"sliding" measure of radius of curvature

[fabrice senger-2]

Hi,

in addition to my previous mail I send a new code wich works in my hands.
With this code you create evenly spaced spot selection along a line you
draw and then calculates the radius of curvature for each consecutive
triplet of points.

For those who are interested I welcome your comments.

This is some code put together from other sources, so, thank you all.
(especially the implementation taken from M. Doube for the calculation of
radius of curvature)

Fabrice.


// Generate an array of selections along a line.


setTool("polyline");
waitForUser("Draw your line selection");


 inc = 30;

 run("Fit Spline", "straighten");
 wait(2000);
 getSelectionCoordinates(x, y);
 for (h=0; h<x.length; h+=inc) {
    makePoint(x[h]-inc, y[h]-inc);
    run("Add to Manager");

 }

 // cycle trough ROIs, select 3 points and fit circle to output radius of
curvature.
myID = getImageID();

setBatchMode(true);
 p = roiManager("count");
 for (j=0; j<p-2; j++) {

a=j;
b=j+1;
c=j+2;
a1=newArray(a,b,c);

roiManager("Select",a1);

     roiManager("Combine"); //unfortunatly, at this point a message
pop's-up...
     roiManager("Add"); // a new 3 point selection is added to the manager,
at the end of the list
     m = roiManager("count");

roiManager("select", m-1);// this is the last element of the roi list, the
generated 3 point selection
getSelectionCoordinates(x, y); //this part of code has been initially
developed by M. DOUBE...
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."); I suppressed this part of
code to not interupt the macro if 3 points are colinear...hope this will
not generate problems...
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)");

roiManager("select", m-1); //here I suppress the last ROI (the 3 point
selection...I dont want to populate the manager too much)
roiManager("Delete");
run("Select None");


 }
 setBatchMode(false);

[ckhripin]

Fabrice,

I tried your macro and had a problem. When the object I am trying to get the curvature of is close to horizontal or vertical, even if its straight, I get some low curvature segments. I modified your code so that inc = 5 and my curves are about 50 pixels in length. Any help would be appreciated.

Constantine