Completing contour for "broken" image

Question

I have a noisy gray-scale image for which I want to segment/mask the large arc spanning the image from the rest. I intend to mask the arc and all of the pixels above the arc.

To do this, I have thresholded the image to create a binary image, and used cv2.findContours() to trace the outline of the arc.

Original Image:

Image after Otsu threshold:

Threshold + Closing:

Contours of closed image:

As you can see, the closed image does not create a solid arc. Closing further causes the arc to lose it's shape. The green line is a contour of the closed image. The blue line is created with approxpolyDP() but I can't get it to work. Are there better ways to mask the arc in the image perhaps?

Here is my code:

import cv2, matplotlib
import numpy as np
import matplotlib.pyplot as plt

# read an image
img = cv2.imread('oct.png')

# get gray image and apply Gaussian blur
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
blur = cv2.GaussianBlur(gray, (5, 5), 0)

# get binary image
ret, thresh = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)

# close image to "solidify" it
kernel = np.ones((3,3),np.uint8)
closing = cv2.morphologyEx(thresh,cv2.MORPH_CLOSE,kernel, iterations = 3)


# find contours
(_, contours, _) = cv2.findContours(closing, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)

cnt = contours[0]
max_area = cv2.contourArea(cnt)

for cont in contours:
    if cv2.contourArea(cont) > max_area:
        cnt = cont
        max_area = cv2.contourArea(cont)

# define main arc contour approx. and hull
perimeter = cv2.arcLength(cnt, True)
epsilon = 0.1 * cv2.arcLength(cnt, True)
approx = cv2.approxPolyDP(cnt, epsilon, True)

# hull = cv2.convexHull(cnt)

# cv2.isContourConvex(cnt)
imgcopy = np.copy(img)
cv2.drawContours(imgcopy, [cnt], -1, (0, 255, 0), 3)
cv2.drawContours(imgcopy, [approx], -1, (0, 0, 255), 3)

# plot figures
plt.figure(1)
plt.imshow(imgcopy, cmap="gray")
plt.figure(2)
plt.imshow(thresh, cmap="gray")
plt.figure(3)
plt.imshow(closing, cmap="gray")

Answer 1

I would suggest to use a RANSAC method to fit 2 ellipses using edge information of the arc. Edge can be obtain simply by using canny or any other method you see fit. Of course this method can only work if the arc is elliptical. If its a straight line, you can replace the ellipse fitting part with a line fitting part.

Here is the result:

Here is the code:

import numpy as np
import cv2
import random as rp


def ransac_ellipse(iter, srcimg, x, y):

    x_size = np.size(x)
    best_count = x_size
    for i in range(iter):

        base = srcimg.copy()

        # get 5 random points
        r1 = int(rp.random() * x_size)
        r2 = int(rp.random() * x_size)  
        r3 = int(rp.random() * x_size)
        r4 = int(rp.random() * x_size)
        r5 = int(rp.random() * x_size)  

        p1 = (x[r1],y[r1])
        p2 = (x[r2],y[r2])
        p3 = (x[r3],y[r3])
        p4 = (x[r4],y[r4])
        p5 = (x[r5],y[r5])

        p_set = np.array((p1,p2,p3,p4,p5))

        # fit ellipse
        ellipse = cv2.fitEllipse(p_set)

        # remove intersected ellipse
        cv2.ellipse(base,ellipse,(0),1)

        # count remain
        local_count = cv2.countNonZero(base)

        # if count is smaller than best, update
        if local_count < best_count:
            best_count = local_count
            best_ellipse = ellipse


    return best_ellipse

img = cv2.imread('arc.jpg',0)

# Speed up and remove noise
small = cv2.resize(img,(0,0),fx = 0.25,fy = 0.25)

# remove remaining noise
median = cv2.medianBlur(small,21)

# get canny edge
edge = cv2.Canny(median,180,20)


cv2.imshow("Edge",edge)


# obtain the non zero locations
y, x = np.where(edge > 0)
# ransac ellipse to get the outter circle
ellipse1 = ransac_ellipse(10000,edge,x,y)


# remove the outter circle
cv2.ellipse(edge,ellipse1,(0),2)

# ransac ellipse to get the inner circle
y, x = np.where(edge > 0)
ellipse2 = ransac_ellipse(10000,edge,x,y)

disp = cv2.cvtColor(small,cv2.COLOR_GRAY2BGR)
cv2.ellipse(disp,ellipse1,(0,0,255),1)
cv2.ellipse(disp,ellipse2,(0,0,255),1)


cv2.imshow("result",disp)
cv2.waitKey(0)

Answer 2

You are on the right path. Your closing will likely work better if you smooth the image a little bit first. I like to apply the thresholding at the end, after the morphological operations. In this case, it doesn't really matter which the order for closing and thresholding is, but keeping thresholding at the end helps later when refining the pre-processing. Once you threshold you loose a lot of information, you need to make sure you preserve all the information you will need, and thus filtering the image properly before thresholding is important.

Here is a quick attempt, I'm sure it can be refined:

import matplotlib.pyplot as pp
import PyDIP as dip

img = pp.imread('/Users/cris/Downloads/MipBB.jpg')
img = img[:,:,0]
smooth = dip.Gauss(img, [3])        # Gaussian smoothing with sigma=3
smooth = dip.Closing(smooth, 25)    # Note! This uses a disk SE with diameter 25 pixels
out = dip.Threshold(smooth, 'triangle')[0]
pp.imsave('/Users/cris/Downloads/MipBB_out.png', out)

I used the triangle threshold method (also known as the chord method, or skewed bi-modality threshold, see P.L. Rosin, "Unimodal thresholding", Pattern Recognition 34(11):2083-2096, 2001) because it works better in this case.

The code uses PyDIP, but I'm sure you can re-create the same process using OpenCV.