# 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 = matB + mean(xW)
matC = matB + mean(yW)
// alpha is the radius squared
Variable alpha = matB^2 + matB^2 + ( (Suu + Svv) / numpnts(xW) )
// clean up 1D waves
KillWaves/Z xW,yW

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. Forum Support Gallery