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

ThreadSafe Function CalcMaxDistance(wx,wy,np,pp)
	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