Detecting Elliptical Shape in Noisy Image

I’m trying to detect the region here (circled in red) more effectively. As it currently stands, I have a few steps to get the area:

  1. Brighten the input image

2

to increase contrast and likelihood to pick up edges of the image file to get this image

3

  1. Crop and threshold the region of interest, and add Gaussian blur to get this image:

4

  1. Use OpenCV to detect Hough Circles on the thresholded image
  2. Select the top 10 largest circles found, then choose the one closest to the grid intersection (in the code as fv_cx and fv_cy) vertically and closest to the edge of the image horizontally.

While this works well in general

1

often times it misses the right circle

5

or it encircles an area too small

6.

Using these input images, is there a better way to work on this problem?

This is my code so far:

from __future__ import print_function
import pandas as pd
from pandas.api.types import is_numeric_dtype
import os
from PIL import Image, ImageDraw, ImageFont
import math
import cv2
import matplotlib.pyplot as plt
import time
import re
import csv
from skimage import data, color, io, img_as_ubyte
from skimage.transform import hough_circle, hough_circle_peaks, hough_ellipse
from skimage.feature import canny
from skimage.draw import circle_perimeter, ellipse_perimeter
from skimage.util import img_as_ubyte
from builtins import input
import numpy as np
   
def append_list_as_row(file_name, list_of_elem):
    with open(file_name, 'a+', newline='', encoding='utf-8') as write_obj:
        csv_writer = csv.writer(write_obj, dialect='excel')
        csv_writer.writerow(list_of_elem)

def round_up_to_odd(f):
    return int(np.ceil(f) // 2 * 2 + 1)
# Folder path here 
folder = r""
csv_file = folder + os.sep + "Measurements.csv"
csv_file2 = folder + os.sep + "Measurements2.csv"
df2 = pd.DataFrame(columns = ["filepath","od_cx","od_cy", "fv_x", "fv_y"])

for subdir, dirs, files in os.walk(folder):
    for file in files:
        #print os.path.join(subdir, file)
        filepath = subdir + os.sep + file
        if filepath.endswith(".jpeg") or filepath.endswith(".tiff") and not filepath.endswith("_OD.tiff") and not filepath.endswith("_bright.tiff") and not filepath.endswith("_FV.tiff") and not filepath.endswith("_mask.tiff"):
            og_cv = cv2.imread(filepath, cv2.IMREAD_COLOR)
            if "left" in str(filepath):
                od = "left"
            elif "right" in str(filepath):
                od = "right"
            
            OD_path = subdir + os.sep + "OD"
            if not os.path.exists(str(OD_path)):
                os.mkdir(str(OD_path))

            OD = OD_path + os.sep + str(os.path.splitext(file)[0]) + "_OD.tiff"
            fovea_path = OD_path + os.sep + str(os.path.splitext(file)[0]) + "_FV.tiff"

            temp_path = subdir + os.sep + "Temp"
            if not os.path.exists(str(temp_path)):
                os.mkdir(str(temp_path))

            bright = temp_path + os.sep + str(os.path.splitext(file)[0]) + "_bright.tiff"
            thresholded_od = temp_path + os.sep + str(os.path.splitext(file)[0]) + "_thresholded_OD.tiff"
            thresholded_fv = temp_path + os.sep + str(os.path.splitext(file)[0]) + "_thresholded_FV.tiff"
            mask_file = temp_path + os.sep + str(os.path.splitext(file)[0]) + "_mask.tiff"

            ## Fovea
            image = cv2.imread(filepath)
            image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)
            h = image.shape[0]
            w = image.shape[1]

            # loop over the image
            fv_cx = []
            fv_cy = []
            for y in range(0, h):
                for x in range(0, w):
                # threshold the pixel
                    if np.all(image[y, x] == (255, 0, 255)) and np.all(image[y, x+3] == (0, 255, 255)) and np.all(image[y, x-3] == (0, 255, 255)):
                        print("Found fovea")
                        fv_cx.append(x)
                        fv_cy.append(y)
                        # Draw them
                        # fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(10, 4))
                        # image = color.gray2rgb(image)
                        image_draw = image
                        if image.shape[2] == 3:
                            image[fv_cy, fv_cx] = (220, 20, 20)
                        if image.shape[2] == 4:
                            image[fv_cy, fv_cx] = (220, 20, 20, 0)
                        plt.imsave(OD, image)
                        print(fv_cx, fv_cy)
                    # else:
                    #     fv_cx = "No fv_cx"
                    #     fv_cy = "No fv_cy"

            ## Find image dimensions
            source_img = Image.open(filepath)
            width, height = source_img.size
            x_max = int(width)
            y_max = int(height)
            print(x_max)
            print(y_max)
                        
            #Load image
            im  = cv2.imread(filepath, cv2.IMREAD_COLOR)
            background = Image.open(filepath).convert('RGB')
            width, height = background.size
            x_max = int(width)
            y_max = int(height)

            # Brightness adjustment - https://docs.opencv.org/3.4/d3/dc1/tutorial_basic_linear_transform.html
            new_image = np.zeros(im.shape, im.dtype)
            alpha = 1.0 # contrast control
            beta = 0 # brightness control
            new_image = cv2.convertScaleAbs(im, alpha=alpha, beta=100)
            # cv2.imshow('New Image', new_image)
            # cv2.waitKey(0)
            cv2.imwrite(bright, new_image)
            new_image = cv2.imread(bright, cv2.IMREAD_COLOR)

            ## OD
            #Convert to HLS, so we can remove the saturated fovea
            HLS = cv2.cvtColor(new_image,cv2.COLOR_BGR2HLS)
            Schannel = HLS[:,:,2]
            mask = cv2.inRange(Schannel, 0, 0)
            # res = cv2.bitwise_and(new_image,new_image, mask= mask)
            new_image = cv2.cvtColor(new_image,cv2.COLOR_BGR2GRAY)
            
            thresh_x = round_up_to_odd((21/1033) * width)
            thresh_x = 21

            #### Thresholding Example Options
            # img = cv2.bitwise_and(new_image,new_image, mask= mask)
            img = cv2.medianBlur(new_image,5)
            
            pil_im  = Image.fromarray(mask)
            # mask_width, mask_height = pil_im.size
            mask_width, mask_height = (165 * (width/290)), (165 * (width/290))
            print(width, height)
            print(mask_width, mask_height)
            margin = 10

            if "_L" in filepath or "OS" in filepath:
                x_center = width/2
                crop_x_start = 0
                crop_x_stop = int(x_center-(mask_width/2)) + margin
                crop_img = img[0:height, crop_x_start:crop_x_stop]
                # cv2.imshow("cropped", crop_img)
                cv2.waitKey()
            if "_R" in filepath or "OD" in filepath:
                x_center = width/2
                crop_x_start = int((x_center+(mask_width/2))) - margin
                crop_x_stop = width
                crop_img = img[0:height, crop_x_start:crop_x_stop]
                # cv2.imshow("cropped", crop_img)
                cv2.waitKey()

            th2 = cv2.adaptiveThreshold(crop_img,255,cv2.ADAPTIVE_THRESH_MEAN_C,\
                        cv2.THRESH_BINARY,thresh_x,2)
            th2 = cv2.GaussianBlur(th2,(21,21),0)
            # cv2.imshow("cropped", th2)
            cv2.waitKey()
            cv2.imwrite(thresholded_od, th2)

            ## Hough Circle
            # Load picture and detect edges
            image = img_as_ubyte(th2)
            edges = canny(image, sigma=3, low_threshold=10, high_threshold=50)

            # Detect two radii
            x=50
            y=500
            z=2
            start = math.ceil((x/1033) * width)
            stop = math.ceil((y/1033) * width)
            step = math.ceil((z/1033) * width)
            
            hough_radii = np.arange(start, stop, step)
            hough_res = hough_circle(edges, hough_radii)

            if fv_cy != []:
                # Select the most prominent 3 circles
                accums, cx, cy, radii = hough_circle_peaks(hough_res, hough_radii,
                                                        total_num_peaks=10)

                df = pd.DataFrame(columns = ["index", "distance", "area", "cX", "cY"])
                idx = (0)
                # Draw them
                fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(10, 4))
                # image = color.gray2rgb(image)
                image = io.imread(filepath)
                cx = (cx + crop_x_start)
                for center_y, center_x, radius in zip(cy, cx, radii):
                    # d = math.sqrt(((center_x-fv_cx)**2) + ((center_y-fv_cy)**2))
                    p = abs(center_y - fv_cy)
                    q = -1 * abs(center_x - fv_cx)
                    d = p + q
                    print(d)
                    area = math.pi * (radius**2)
                    df.loc[idx, 'index'] = idx
                    df.loc[idx, 'area'] = int(area)
                    df.loc[idx, 'distance'] = int(d)
                    df.loc[idx, 'cX'] = int(center_x)
                    df.loc[idx, 'cY'] = int(center_y)
                    df.loc[idx, 'radius'] = int(radius)
                    idx += 1
                df['distance'] = pd.to_numeric(df['distance'])
                df['radius'] = pd.to_numeric(df['radius'])
                print("DF?")
                print(df)
                if len(df["distance"]) > 0:
                    print("pass")
                    df_radius = df.nsmallest(3, 'distance')
                    print(df_radius)
                    if (df_radius['radius'].max()-df_radius['radius'].min()) < 3:
                        idx = df_radius['radius'].idxmax()
                    else:
                        idx = df['distance'].idxmin()
                    center_y = int(df.loc[idx, 'cY'])
                    center_x = int(df.loc[idx, 'cX'])
                    radius = int(df.loc[idx, 'radius'])
                    print(center_y, center_x, radius)
                    # Draw them
                    # fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(10, 4))
                    # image = color.gray2rgb(image)
                    image = io.imread(filepath)
                    image = image_draw
                    circy, circx = circle_perimeter(center_y, center_x, radius,
                                                    shape=image.shape)
                    print(image.shape)
                    print(image.shape[2])
                    if image.shape[2] == 3:
                        image_draw[circy, circx] = (220, 20, 20)
                        circy, circx = circle_perimeter(center_y, center_x, 0,
                                                        shape=image.shape)
                        image_draw[circy, circx] = (220, 20, 20)
                    if image.shape[2] == 4:
                        image_draw[circy, circx] = (220, 20, 20, 0)
                        circy, circx = circle_perimeter(center_y, center_x, 0,
                                                        shape=image.shape)
                        image_draw[circy, circx] = (220, 20, 20, 0)
                    # final = ax.imshow(image, cmap=plt.cm.gray)
                    # fig = plt.show()
                    ## Need to fix saving
                    plt.imsave(OD, image_draw)
            else:
                hough_radii = np.arange(start, stop, step)
                hough_res = hough_circle(edges, hough_radii)

                # Select the most prominent 3 circles
                accums, cx, cy, radii = hough_circle_peaks(hough_res, hough_radii,
                                                        total_num_peaks=1)
                if cx != None:
                    print("Found OD")
                    od_cx = (re.search(r"\[([A-Za-z0-9_]+)\]",str(cx))).group(1)
                    od_cy = (re.search(r"\[([A-Za-z0-9_]+)\]",str(cy))).group(1)
                else:
                    od_cx = "Not found"
                    od_cy = "Not found"

                # Draw them
                #fig, ax = plt.subplots(ncols=1, nrows=1, #figsize=(10, 4))
                # image = color.gray2rgb(image)
                image = io.imread(filepath)
                image = image_draw
                for center_y, center_x, radius in zip(cy, cx, radii):
                    circy, circx = circle_perimeter(center_y, center_x, radius,
                                                    shape=image.shape)
                    print(image.shape)
                    print(image.shape[2])
                    if image.shape[2] == 3:
                        image_draw[circy, circx] = (220, 20, 20)
                        circy, circx = circle_perimeter(center_y, center_x, 0,
                                                        shape=image.shape)
                        image_draw[circy, circx] = (220, 20, 20)
                    if image.shape[2] == 4:
                        image_draw[circy, circx] = (220, 20, 20, 0)
                        circy, circx = circle_perimeter(center_y, center_x, 0,
                                                        shape=image.shape)
                        image_draw[circy, circx] = (220, 20, 20, 0)
                # final = ax.imshow(image, cmap=plt.cm.gray)
                # fig = plt.show()
                ## Need to fix saving
                plt.imsave(OD, image_draw)
            append_list_as_row(csv_file,[filepath,center_x,center_y, fv_cx, fv_cy])
            plt.close('all')
            df2 = df2.append({"filepath":filepath,"od_cx":center_x, "od_cy":center_y, "fv_x":fv_cx, "fv_y":fv_cy}, ignore_index=True)
            print(df2)
df2.to_csv(csv_file2)

Hi, sorry for the late reply !

Just to be sure, you have a grayscale image that you load as RGB and draw this central square right?
Do you need to draw this square for the processing ?
You would potentially get a better signal loading the image as grayscale rather than RGB (at least if the image is 16-bit).
image = cv2.imread(filepath, cv2.IMREAD_UNCHANGED)

I see that you detect the edges with Canny before the Hough circle.
I am not sure if you do it on the gray image (which is what I would do) or on the thresholded version. Maybe you can try to adjust the values below to avoid too small circles. Also checking the result of the edge map is important, there are a few parameters for the canny edge detection, so you have to make sure you can distinguish the circle in this edge map.

A different solution could be to try some sort of spot detection, with a DoG filter (Difference of Gaussians). With a large Gaussian kernel (for one of the 2 gaussians) that would smooth out everything and emphasize round shapes. From there with a maxima detector you could get the center of the circle.
However to recover the outline maybe a seeded watershed from the center or an active contour could grow a mask from the detected center to the actual circle perimenter.

Good luck!