Re: bending and radius of curvature

[ved sharma]

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");

[Michael Doube]

Hi Ved, Fabrice,

The FitCircle class is now being maintained at:
https://github.com/mdoube/BoneJ/blob/master/src/org/doube/geometry/FitCircle.java

And there are unit tests at:
https://github.com/mdoube/BoneJ/blob/master/test/org/doube/geometry/FitCircleTest.java

Also, it would be better to use the public methods in the class, rather
than reimplementing them in macro language, although the call() macro
function doesn't seem to handle non-string arguments (we need to pass a
double[] argument to the circle fitting methods).  Perhaps another
scripting language (javascript?) or Java plugin would be better.

Michael

> 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

[Rasband, Wayne (NIH/NIMH) [E]]

On Dec 19, 2010, at 3:22 PM, Ved Sharma wrote:

> 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
> <CircleFit_macro.txt>

Michael Doube's FitCircle.prattNewton() circle fitting method is built into the ImageJ 1.44m daily build as the Edit>Selection>Fit Circle command. Use the Measure command, with "Centroid" and "Bounding rectangle" enabled in Set Measurements, to get the center and diameter of the fitted circle.

-wayne

[vischer]

Hi Wayne and others,

For closed shapes, I rather would associate "Fit Circle" to be similar to "Fit Ellipse", i.e. same area and centroid. This is currently not the case: a longish, symmetrical particle outline could result in a very large radius (see macro below). Also if I remove the "Fit Spline" in the macro, I don't understand the meaning of the resulting circle.

What about applying the "Fit Circle" only to open shapes (MultiPoint, Segmented Line), while closed shapes just would return a circle with same area and centroid?

Norbert Vischer


newImage("A", "8-bit White", 550, 550, 1);
makeOval(80, 200, 200, 50);
setKeyDown("shift");
makeOval(90, 220, 10, 10);//so it behaves like a traced roi
run("Fit Spline");//optional
run("Invert")
wait(1000);
run("Fit Circle");

[Gabriel Landini]

On Tuesday 21 December 2010 11:49:10 you wrote:
> For closed shapes, I rather would associate "Fit Circle" to be similar to
> "Fit Ellipse", i.e. same area and centroid. This is currently not the case:
> a longish, symmetrical particle outline could result in a very large radius
> (see macro below). Also if I remove the "Fit Spline" in the macro, I don't
> understand the meaning of the resulting circle.

I guess that it is assuming a cloud of points (in this case the points forming
the ellipse) as defining an arc of the circle (instead of the whole circle).
Not sure how the algorithm works but maybe it looks at minimising residuals.
If so, the arc might have smaller residuals than a circle which intuitively
looks similarly sized to the cloud?


Cheers

Gabriel

[Rasband, Wayne (NIH/NIMH) [E]]

On Dec 21, 2010, at 5:49 AM, Norbert Vischer wrote:

> Hi Wayne and others,
>
> For closed shapes, I rather would associate "Fit Circle" to be similar to "Fit Ellipse", i.e. same area and centroid. This is currently not the case: a longish, symmetrical particle outline could result in a very large radius (see macro below). Also if I remove the "Fit Spline" in the macro, I don't understand the meaning of the resulting circle.
>
> What about applying the "Fit Circle" only to open shapes (MultiPoint, Segmented Line), while closed shapes just would return a circle with same area and centroid?
>
> Norbert Vischer
>
>
> newImage("A", "8-bit White", 550, 550, 1);
> makeOval(80, 200, 200, 50);
> setKeyDown("shift");
> makeOval(90, 220, 10, 10);//so it behaves like a traced roi
> run("Fit Spline");//optional
> run("Invert")
> wait(1000);
> run("Fit Circle");

In the 1.44m5 daily build, Edit>Selection>Fit Circle only fits a circle to open shapes (lines and points). For closed shapes, it creates a circle with the same area and centroid.

-wayne for