Posted by
ved sharma on
Dec 19, 2010; 8:22pm
URL: http://imagej.273.s1.nabble.com/Re-bending-and-radius-of-curvature-tp3686117.html
Dear Fabrice,
I don't know if you still need it but Michael Doube has written a very nice plugin called FitCircle, which I think will do what you are looking for. Check out his page:
http://doube.org/plugins. It fits circle through the data points using a variety of user selected methods. The original codes for fitting the circle were from Nikolai Chernov's MATLAB scripts (
http://www.math.uab.edu/~chernov/cl/MATLABcircle.html). I took one of these circle fitting methods and rewrote the code in a macro. Check out the attached file.
Ved
/*
This macro fits a circle to the user-selected points in an image.
Reference: Pratt V., Direct least-squares fitting of algebraic surfaces", Computer Graphics, Vol. 21, pages 145-152 (1987).
Original code: Nikolai Chernov's MATLAB script for Newton-based Pratt fit.
(
http://www.math.uab.edu/~chernov/cl/MATLABcircle.html)
Ved Sharma, version: December 18, 2010
*/
requires("1.43k");
if (nImages == 0)
exit("No image is open.");
if (isOpen("ROI Manager")) {
showMessageWithCancel("Circle Fit Macro.", "This macro will reset your ROI Manager.\nDo you want to proceed?");
roiManager("reset");
}
setTool("multipoint");
waitForUser("Circle Fit Macro.", "Multipoint tool selected.\nSelect points in your image, then click OK.");
roiManager("Add"); roiManager("Select", 0); roiManager("Rename", "points");
getSelectionCoordinates(x, y);
n=x.length;
if(n<3)
exit("At least 3 points are required to fit a circle.");
sumx = 0; sumy = 0;
for(i=0;i<n;i++) {
sumx = sumx + x[i];
sumy = sumy + y[i];
}
meanx = sumx/n;
meany = sumy/n;
X = newArray(n); Y = newArray(n);
Mxx=0; Myy=0; Mxy=0; Mxz=0; Myz=0; Mzz=0;
for(i=0;i<n;i++) {
X[i] = x[i] - meanx;
Y[i] = y[i] - meany;
Zi = X[i]*X[i] + Y[i]*Y[i];
Mxy = Mxy + X[i]*Y[i];
Mxx = Mxx + X[i]*X[i];
Myy = Myy + Y[i]*Y[i];
Mxz = Mxz + X[i]*Zi;
Myz = Myz + Y[i]*Zi;
Mzz = Mzz + Zi*Zi;
}
Mxx = Mxx/n;
Myy = Myy/n;
Mxy = Mxy/n;
Mxz = Mxz/n;
Myz = Myz/n;
Mzz = Mzz/n;
Mz = Mxx + Myy;
Cov_xy = Mxx*Myy - Mxy*Mxy;
Mxz2 = Mxz*Mxz;
Myz2 = Myz*Myz;
A2 = 4*Cov_xy - 3*Mz*Mz - Mzz;
A1 = Mzz*Mz + 4*Cov_xy*Mz - Mxz2 - Myz2 - Mz*Mz*Mz;
A0 = Mxz2*Myy + Myz2*Mxx - Mzz*Cov_xy - 2*Mxz*Myz*Mxy + Mz*Mz*Cov_xy;
A22 = A2 + A2;
epsilon=1e-12;
ynew=1e+20;
IterMax=20;
xnew = 0;
// Newton's method starting at x=0
for(iter=1;i<=IterMax;i++) {
yold = ynew;
ynew = A0 + xnew*(A1 + xnew*(A2 + 4.*xnew*xnew));
if (abs(ynew)>abs(yold)) {
print("Newton-Pratt goes wrong direction: |ynew| > |yold|");
xnew = 0;
i = IterMax+1;
}
else {
Dy = A1 + xnew*(A22 + 16*xnew*xnew);
xold = xnew;
xnew = xold - ynew/Dy;
if (abs((xnew-xold)/xnew) < epsilon)
i = IterMax+1;
else {
if (iter >= IterMax) {
print("Newton-Pratt will not converge");
xnew = 0;
}
if (xnew<0) {
// print("Newton-Pratt negative root: x = "+xnew);
xnew = 0;
}
}
}
}
DET = xnew*xnew - xnew*Mz + Cov_xy;
CenterX = (Mxz*(Myy-xnew)-Myz*Mxy)/(2*DET);
CenterY = (Myz*(Mxx-xnew)-Mxz*Mxy)/(2*DET);
radius = sqrt(CenterX*CenterX + CenterY*CenterY + Mz + 2*xnew);
if(isNaN(radius))
exit("Points selected are collinear.");
CenterX = CenterX + meanx;
CenterY = CenterY + meany;
print("Radius = "+d2s(radius,2));
print("Center coordinates: ("+d2s(CenterX,2)+", "+d2s(CenterY,2)+")");
print("(all units in pixel)");
makeOval(CenterX-radius, CenterY-radius, 2*radius, 2*radius);
roiManager("Add"); roiManager("Select", 1); roiManager("Rename", "circle fit");
setTool("point");
makePoint(CenterX, CenterY);
roiManager("Add"); roiManager("Select", 2); roiManager("Rename", "center");
roiManager("Show All without labels");