CDR_Calculation.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "CDR Calculation",
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/kaylaque/glauco-scan/blob/main/CDR_Calculation.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1hnHuYe3zNCX"
      },
      "source": [
        "import cv2\n",
        "import sys\n",
        "import glob\n",
        "import pandas as pd\n",
        "import os\n",
        "from PIL import Image\n",
        "import math\n",
        "from pylab import*\n",
        "from scipy import signal\n",
        "import numpy as np\n",
        "from matplotlib import pyplot as plt\n",
        "from skimage.feature import greycomatrix, greycoprops\n",
        "from skimage import data\n",
        "clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(10,10))\n",
        "def cdr(cup,disc,plot):\n",
        "    \n",
        "    #morphological closing and opening operations\n",
        "    R1 = cv2.morphologyEx(cup, cv2.MORPH_CLOSE, cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(2,2)), iterations = 1)\n",
        "    r1 = cv2.morphologyEx(R1, cv2.MORPH_OPEN, cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(7,7)), iterations = 1)\n",
        "    R2 = cv2.morphologyEx(r1, cv2.MORPH_CLOSE, cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(1,21)), iterations = 1)\n",
        "    r2 = cv2.morphologyEx(R2, cv2.MORPH_OPEN, cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(21,1)), iterations = 1)\n",
        "    R3 = cv2.morphologyEx(r2, cv2.MORPH_CLOSE, cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(33,33)), iterations = 1)\t\n",
        "    r3 = cv2.morphologyEx(R3, cv2.MORPH_OPEN, cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(43,43)), iterations = 1)\n",
        "\n",
        "    img = clahe.apply(r3)\n",
        "    \n",
        "    \n",
        "    ret,thresh = cv2.threshold(cup,127,255,0)\n",
        "    # print(thresh)\n",
        "    # print(len(cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)))\n",
        "    contours,hierarchy = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE) #Getting all possible contours in the segmented image\n",
        "    cup_diameter = 0\n",
        "    largest_area = 0\n",
        "    el_cup = contours[0]\n",
        "    if len(contours) != 0:\n",
        "        for i in range(len(contours)):\n",
        "            if len(contours[i]) >= 5:\n",
        "                area = cv2.contourArea(contours[i]) #Getting the contour with the largest area\n",
        "                if (area>largest_area):\n",
        "                    largest_area=area\n",
        "                    index = i\n",
        "                    el_cup = cv2.fitEllipse(contours[i])\n",
        "                \n",
        "    cv2.ellipse(img,el_cup,(140,60,150),3)  #fitting ellipse with the largest area\n",
        "    x,y,w,h = cv2.boundingRect(contours[index]) #fitting a rectangle on the ellipse to get the length of major axis\n",
        "    cup_diameter = max(w,h) #major axis is the diameter\n",
        "\n",
        "    #morphological closing and opening operations\n",
        "    R1 = cv2.morphologyEx(disc, cv2.MORPH_CLOSE, cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(2,2)), iterations = 1)\n",
        "    r1 = cv2.morphologyEx(R1, cv2.MORPH_OPEN, cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(7,7)), iterations = 1)\n",
        "    R2 = cv2.morphologyEx(r1, cv2.MORPH_CLOSE, cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(1,21)), iterations = 1)\n",
        "    r2 = cv2.morphologyEx(R2, cv2.MORPH_OPEN, cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(21,1)), iterations = 1)\n",
        "    R3 = cv2.morphologyEx(r2, cv2.MORPH_CLOSE, cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(33,33)), iterations = 1)\n",
        "    r3 = cv2.morphologyEx(R3, cv2.MORPH_OPEN, cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(43,43)), iterations = 1)\n",
        "\n",
        "    img2 = clahe.apply(r3)\n",
        "    \n",
        "    ret,thresh = cv2.threshold(disc,127,255,0)\n",
        "    contours,hierarchy = cv2.findContours(thresh, cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_SIMPLE) #Getting all possible contours in the segmented image\n",
        "    disk_diameter = 0\n",
        "    largest_area = 0\n",
        "    el_disc = el_cup\n",
        "    if len(contours) != 0:\n",
        "          for i in range(len(contours)):\n",
        "            if len(contours[i]) >= 5:\n",
        "                area = cv2.contourArea(contours[i]) #Getting the contour with the largest area\n",
        "                if (area>largest_area):\n",
        "                    largest_area=area\n",
        "                    index = i\n",
        "                    el_disc = cv2.fitEllipse(contours[i])\n",
        "                    \n",
        "    cv2.ellipse(img2,el_disc,(140,60,150),3) #fitting ellipse with the largest area\n",
        "    x,y,w,h = cv2.boundingRect(contours[index]) #fitting a rectangle on the ellipse to get the length of major axis\n",
        "    disk_diameter = max(w,h) #major axis is the diameter\n",
        "                \n",
        "    if plot:\n",
        "        plt.imshow(img2, cmap = 'gray', interpolation = 'bicubic')\n",
        "        plt.axis(\"off\")\n",
        "        plt.title(\"Optic Disk\")\n",
        "        plt.show()\n",
        "        plt.imshow(img)\n",
        "        plt.axis(\"off\")\n",
        "        plt.title(\"Optic Cup\")\n",
        "        plt.show()\n",
        "        \n",
        "    if(disk_diameter == 0): return 1 # if disc not segmented properly then cdr might be infinity\n",
        "    cdr = cup_diameter/disk_diameter #ration of major axis of cup and disc\n",
        "    return cdr\n",
        "'''\n",
        "# use this if the input from a directory\n",
        "cup = cv2.imread('./cup.png',0) \n",
        "disc = cv2.imread('./disk.png',0) \n",
        "'''\n",
        "cdr_cal = cdr(cup,disc,True)\n",
        "print(cdr_cal)"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}