Created
May 24, 2022 09:08
-
-
Save quantixed/3249ead928e66f29f554cf7dda8fb10e to your computer and use it in GitHub Desktop.
Igor function to fit a circle to 2D coords (sparse, irregular, arc). Returns the radius in units of xy
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment