Open CV trivial circle detection -- how to get least squares instead of a contour?

Here is another way to fit a circle by getting the equivalent circle center and radius from the binary image using connected components and drawing a circle from that using Python/OpenCV/Skimage.

Input:

import cv2
import numpy as np
from skimage import measure

# load image and set the bounds
img = cv2.imread("dark_circle.png")

# convert to grayscale
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# blur
blur = cv2.GaussianBlur(gray, (3,3), 0)

# threshold
thresh = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)[1]

# apply morphology open with a circular shaped kernel
kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (5,5))
binary = cv2.morphologyEx(thresh, cv2.MORPH_OPEN, kernel, iterations=2)

# find contour and draw on input (for comparison with circle)
cnts = cv2.findContours(binary, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
cnts = cnts[0] if len(cnts) == 2 else cnts[1]
c = cnts[0]
result = img.copy()
cv2.drawContours(result, [c], -1, (0, 255, 0), 1)

# find radius and center of equivalent circle from binary image and draw circle
# see https://scikit-image.org/docs/dev/api/skimage.measure.html#skimage.measure.regionprops
# Note: this should be the same as getting the centroid and area=cv2.CC_STAT_AREA from cv2.connectedComponentsWithStats and computing radius = 0.5*sqrt(4*area/pi) or approximately from the area of the contour and computed centroid via image moments.
regions = measure.regionprops(binary)
circle = regions[0]
yc, xc = circle.centroid
radius = circle.equivalent_diameter / 2.0
print("radius =",radius, "  center =",xc,",",yc)
xx = int(round(xc))
yy = int(round(yc))
rr = int(round(radius))
cv2.circle(result, (xx,yy), rr, (0, 0, 255), 1)

# write result to disk
cv2.imwrite("dark_circle_fit.png", result)

# display it
cv2.imshow("image", img)
cv2.imshow("thresh", thresh)
cv2.imshow("binary", binary)
cv2.imshow("result", result)
cv2.waitKey(0)

Result showing contour (green) compared to circle fit (red):

Circle Radius and Center:

radius = 117.6142467296168   center = 220.2169911178609 , 150.26823599797507

A least squares fit method (between the contour points and a circle) can be obtained using Scipy. For example, see:

https://gist.github.com/lorenzoriano/6799568

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html

I would suggest computing a mask as in nathancy's answer, but then simply counting the number of pixels in the mask opening that he computed (which is an unbiased estimate of the area of the hole), and then translating the area to a radius using radius = sqrt(area/pi). This will give you the radius of the circle with the same area as the hole, and corresponds to one method to obtain a best fit circle.

A different way of obtaining a best fit circle is to take the contour of the hole (as returned in cnts by cv.findContours in nethancy's answer), finding its centroid, and then computing the mean distance of each vertex to the centroid. This would correspond approximately* to a least squares fit of a circle to the hole perimeter.

^{* I say approximately because the vertices of the contour are an approximation to the contour, and the distances between these vertices is likely not uniform. The error should be really small though.}

---

Here's code example using DIPlib (disclosure: I'm an author) (note: the import PyDIP statement below requires you install DIPlib, and you cannot install it with pip, there is a binary release for Windows on the GitHub page, or otherwise you need to build it from sources).

import PyDIP as dip
import imageio
import math

img = imageio.imread('https://i.stack.imgur.com/szvc2.jpg')
img = dip.Image(img[:,2600:-1])
img.SetPixelSize(0.01, 'mm')      # Use your actual values!
bin = ~dip.OtsuThreshold(dip.Gauss(img, [3]))
bin = dip.Opening(bin, 25)
#dip.Overlay(img, bin - dip.BinaryErosion(bin, 1, 3)).Show()

msr = dip.MeasurementTool.Measure(dip.Label(bin), features=['Size', 'Radius'])
#print(msr)

print('Method 1:', math.sqrt(msr[1]['Size'][0] / 3.14), 'mm')
print('Method 2:', msr[1]['Radius'][1], 'mm')

The MeasurementTool.Measure function computes 'Size', which is the area; and 'Radius', which returns the max, mean, min and standard deviation of the distances between each boundary pixel and the centroid. From 'Radius', we take the 2nd value, the mean radius.

This outputs:

Method 1: 7.227900647539411 mm
Method 2: 7.225178113501325 mm

But do note that I assigned a random pixel size (0.01mm per pixel), you'll need to fill in the right pixels-to-mm conversion value.

Note how the two estimates are very close. Both methods are good, unbiased estimates. The first method is computationally cheaper.