# Find Smallest Circle Fit to Image Blob

This will seek the smallest circle to fit to a given shape. The inputs wx, wy are the boundaries to the shape.

Further optimizations are likely possible.

Reference Discussions

https://www.wavemetrics.com/forum/general/iso-algorithm-apply-smallest-circle-problem

https://www.wavemetrics.com/forum/general/circle-cut-edges

// calculate the smallest circle to fit boundary defined by wy, wx
Function calc_SmallestCircle(wy,wx)
wave wy, wx

variable np, npl, npr
variable px1, py1, px2, py2, px3, py3
variable vm1, vm2

// restrict to under 1000 points as needed to avoid long times
np = numpnts(wx)
npl = log(np)
if (npl > 3)
npl = trunc(npl)
npr = trunc(np/10^npl)
npl = (npr*10^(npl - 2))/5
duplicate/O wx wdx
duplicate/O wy wdy
Resample/DOWN=(npl) wdx
Resample/DOWN=(npl) wdy
np = numpnts(wdx)
else
duplicate/O wx, wdx
duplicate/O wy, wdy
np = numpnts(wx)
endif

// calculate the distances
make/N=(np)/FREE/D dmaxarray
dmaxarray = 0

MultiThread dmaxarray = CalcMaxDistance(wdx, wdy, np, p)

// get three points
make/N=3/D/FREE wpx, wpy
npr = 0.03*np
WaveStats/Q dmaxarray
vm1 = v_maxloc
wpx[0] = wdx[vm1]
wpy[0] = wdy[vm1]
npl = (vm1 - npr) < 0 ? 0 : vm1 - npr
dmaxarray[npl,vm1 + npr] = 0
WaveStats/Q/R=(inf,vm1) dmaxarray
vm2 = v_maxloc
wpx[1] = wdx[vm2]
wpy[1] = wdy[vm2]
npl = (vm2 + npr) > (np - 1) ? np - 1 : vm2 + npr
dmaxarray[vm2 - npr, npl] = 0
if ((vm2 - vm1) > (np - vm2))
WaveStats/Q/R=(vm1,vm2) dmaxarray
else
WaveStats/Q/R=(vm2,np-1) dmaxarray
endif
wpx[2] = wdx[v_maxloc]
wpy[2] = wdy[v_maxloc]
dmaxarray[v_maxloc] = 0

GetCircumCenter(wpx,wpy)
return 0
end

wave wx, wy
variable np, pp

variable rmax

// create 2-D wave
make/N=(np)/FREE dd = 0
// calculate distances and return maximum
dd = sqrt((wx[p] - wx[pp])^2 + (wy[p] - wy[pp])^2)
rmax = wavemax(dd)

return rmax
end

Function GetCircumCenter(wpx,wpy)
wave wpx, wpy

variable dlen, cx, cy, rc

dlen = 2*(wpy[0]*wpx[2] + wpy[1]*wpx[0] - wpy[1]*wpx[2] - wpy[0]*wpx[1] - wpy[2]*wpx[0] + wpy[2]*wpx[1])

cx = wpy[1]*wpx[0]^2 - wpy[2]*wpx[0]^2 - (wpy[1]^2)*wpy[0] + (wpy[2]^2)*wpy[0]
cx += (wpx[1]^2)*wpy[2] + (wpy[0]^2)*wpy[1] + wpx[2]^2*wpy[0] - (wpy[2]^2)*wpy[1]
cx += - (wpx[2]^2)*wpy[1] - (wpx[1]^2)*wpy[0] + (wpy[1]^2)*wpy[2] - (wpy[0]^2)*wpy[2]
cx /= dlen

cy = (wpx[0]^2)*wpx[2] + (wpy[0]^2)*wpx[2] + (wpx[1]^2)*wpx[0] - (wpx[1]^2)*wpx[2]
cy += + (wpy[1]^2)*wpx[0] - (wpy[1]^2)*wpx[2] - (wpx[0]^2)*wpx[1] - (wpy[0]^2)*wpx[1]
cy += - (wpx[2]^2)*wpx[0] + (wpx[2]^2)*wpx[1] - (wpy[2]^2)*wpx[0] + (wpy[2]^2)*wpx[1]
cy /= dlen

rc = sqrt((wpx[0] - cx)^2 + (wpy[0] - cy)^2)

print "Values (cx, cy, rc) are: ", cx, cy, rc

return 0
end

Forum

Support

Gallery