Fit a circle to 2D coords (least squares method)

Function FitCircleTo2DCoords(w)
	Wave w
	// requires 2D numeric wave with two columnns corresponding to xy coords
	if(Dimsize(w,1) != 2)
		return -1
	endif
	// make two 1D waves for x and y coords
	SplitWave/O/NAME="xW;yW;" w
	Wave xW,yW
	// solve in terms of u and v coordinates
	Duplicate/O xw, uW
	Duplicate/O yw, vW
	uW[] = xw[p] - mean(xw)
	vW[] = yw[p] - mean(yw)
	// sigma calcs - store as variables for clean code
	Variable Su = sum(uW)
	Variable Sv = sum(vW)
	MatrixOp/O/FREE uu = uW * uW
	Variable Suu = sum(uu)
	MatrixOp/O/FREE vv = vW * vW
	Variable Svv = sum(vv)
	MatrixOp/O/FREE uv = uW * vW
	Variable Suv = sum(uv)
	MatrixOp/O/FREE uuu = uW * uW * uW
	Variable Suuu = sum(uuu)
	MatrixOp/O/FREE uvv = uW * vW * vW
	Variable Suvv = sum(uvv)
	MatrixOp/O/FREE vvv = vW * vW * vW
	Variable Svvv = sum(vvv)
	MatrixOp/O/FREE vuu = vW * uW * uW
	Variable Svuu = sum(vuu)
	// linear system
	Make/O/N=(2,2) matA = {{Suu,Suv},{Suv,Svv}}
	Make/O/N=(2,1) matB = {0.5 * (Suuu + Suvv),0.5 * (Svvv + Svuu)}
	// Solve it. matB is overwritten with the result uc,vc
	MatrixLinearSolve/O matA matB
	// transform back to x y coordinate system to get xc,yc
	Duplicate/O matB, matC
	// matC is saved and contains the centre of the circle
	matC[0] = matB[0] + mean(xW)
	matC[1] = matB[1] + mean(yW)
	// alpha is the radius squared
	Variable alpha = matB[0]^2 + matB[1]^2 + ( (Suu + Svv) / numpnts(xW) )
	// clean up 1D waves
	KillWaves/Z xW,yW
	
	// return radius
	return sqrt(alpha)
End

--

This code returns the radius of a circle that fits a set of 2D coordinates. Simple least-squares method. Coordinates may be from an arc or other sparse set of points. Centre of circle is stored as a wave matC.