macro to draw a circle from three points

[Christopher Coulon-2]

Before I reinvent the wheel, is there anyone who can share the code for drawing a circle from three points?  

Thanks in advance,

Chris
--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

[Jerome Mutterer-3]

Hi Christopher,
The "Fit Circle" command does this. It's code is at:
https://github.com/imagej/imagej1/blob/master/ij/plugin/Selection.java#L103
Sincerely,
Jerome.


On 13 February 2016 at 04:30, Christopher Coulon <
[hidden email]> wrote:
>
> Before I reinvent the wheel, is there anyone who can share the code for
drawing a circle from three points?
>
> Thanks in advance,
>
> Chris
> --
> ImageJ mailing list: http://imagej.nih.gov/ij/list.html




--
Jerome Mutterer
CNRS - Institut de biologie moléculaire des plantes
12, rue du Général Zimmer
67084 Strasbourg Cedex
www.ibmp.cnrs.fr

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

[Christopher Coulon-2]

Thank you Jerome!  Last night, as I was waiting to go to sleep, I got the inspiration that I could do this by looking through the entire image for the pixel location equidistant from all three points, so I wrote the macro below to test the theory, and it worked!  :-)

newImage("Untitled", "8-bit white", 1000, 1000, 1);
run("Colors...", "foreground=black background=white selection=yellow");
run("Set Measurements...", "area centroid limit redirect=None decimal=8");
makeLine(234, 156, 240, 156);
run("Draw", "slice");
makeLine(369, 120, 372, 120);
run("Draw", "slice");
makeLine(466, 180, 470, 182);
run("Draw", "slice");
setThreshold(0, 128);
run("Analyze Particles...", "size=2-Infinity display exclude clear record");
run("Select None");
resetThreshold();

x0 = getResult("X", 0);
y0 = getResult("Y", 0);
x1 = getResult("X", 1);
y1 = getResult("Y", 1);
x2 = getResult("X", 2);
y2 = getResult("Y", 2);
l0 = 0;
l1 = 0;
l2 = 0;
top = getHeight();
left = getWidth();

center = newArray(2);;
radius = 0;
go = true;

for(x=0; x<getWidth(); x++) {
        for(y = 0; y < getHeight(); y++) {
                l1 = round(sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0)));
                l2 = round(sqrt((x-x1)*(x-x1)+(y-y1)*(y-y1)));
                l3 = round(sqrt((x-x2)*(x-x2)+(y-y2)*(y-y2)));
                if(go && l1 == l2 && l1 == l3) {
                        center[0] = x;
                        center[1] = y;
                        radius = l1;
                        if(y > 0) go = false;
                }
        }
}

print("center x, y = " + center[0] + ", " + center[1]);
print("radius = " + radius);

left = center[0] - radius;
top = center[1] - radius;
diam = radius * 2;
drawOval(left, top, diam, diam);
print("left = " + left + "  top = " + top);

> On Feb 13, 2016, at 1:20 AM, Jerome Mutterer <[hidden email]> wrote:
>
> Hi Christopher,
> The "Fit Circle" command does this. It's code is at:
> https://github.com/imagej/imagej1/blob/master/ij/plugin/Selection.java#L103
> Sincerely,
> Jerome.
>
>
> On 13 February 2016 at 04:30, Christopher Coulon <
> [hidden email]> wrote:
>>
>> Before I reinvent the wheel, is there anyone who can share the code for
> drawing a circle from three points?
>>
>> Thanks in advance,
>>
>> Chris
>> --
>> ImageJ mailing list: http://imagej.nih.gov/ij/list.html
>
>
>
>
> --
> Jerome Mutterer
> CNRS - Institut de biologie moléculaire des plantes
> 12, rue du Général Zimmer
> 67084 Strasbourg Cedex
> www.ibmp.cnrs.fr
>
> --
> ImageJ mailing list: http://imagej.nih.gov/ij/list.html

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

[Peter Haub]

Hi Christopher,

3 points definitely define a unique circle.

Here is the simple macro code to calculate the center and radius of a
circle defined by a Multipoint selection of 3 points:


//  ****************   Macro: Calculate Circle from 3 points
****************


if ( selectionType() == 10 ){ // is MultiPoint selection
    getSelectionCoordinates(x, y);

    if (lengthOf(x) != 3){
       print("Must be 3 points !");
       return;
    }

    b1 = x[0]* x[0] + y[0]*y[0];
    b2 = x[1]* x[1] + y[1]*y[1];
    b3 = x[2]* x[2] + y[2]*y[2];

    detA = 4*( x[0]*(y[1]-y[2]) - y[0]*(x[1]-x[2]) + (x[1]*y[2]-
x[2]*y[1]) );

    xc = 2*( b1*(y[1]-y[2]) - y[0]*(b2-b3) + (b2*y[2]-b3*y[1]) )/detA;
    yc = 2*( x[0]*(b2-b3) - b1*(x[1]-x[2]) + (x[1]*b3-x[2]*b2) )/detA;
    K = 4*( x[0]*(y[1]*b3-y[2]*b2) - y[0]*(x[1]*b3-x[2]*b2) +
b1*(x[1]*y[2]-x[2]*y[1]) )/detA;
    r = sqrt(abs(K + xc*xc + yc*yc));

    print("xc : " + xc);
    print("yc : " + yc);
    print("r : " + r);

   drawOval(xc-r, yc-r, 2*r, 2*r);

}

// *************** end of macro  *******************

Instead of point coordinates you can use your x/y measurements from the
result table. Adapt the code to your needs.


Here is the math explanations:

The general equation of a circle
    0 = sqr(x-xc) + sqr(y-yc)- sqr(r)
can be reformulated as
    0 = sqr(x) + sqr(y) + sqr(xc) + sqr(yc) -2x*xc - 2y*yc - sqr(r)
and further
    sqr(x) + sqr(y) = b = 2+x*xc + 2*y*yc + K

where:
r - Radius
x,y - Point coordinates
xc, yc - Center coordinates
sqr(a) = a*a
K = sqr(r) - sqr(xc) - sqr(yc)

For 3 points the systems of 3 linear equations can be solved for the
unknowns xc, yc, A.
b1 = 2*x1*xc + 2*y1*yc + K
b2 = 2*x2*xc + 2*y2*yc + K
b3 = 2*x3*xc + 2*y3*yc + K

This is in matrix notation
B =  Matrix ( (2*x1,   2*y1,   1,), ( 2*x2, 2*y3, 1 ), ( 2*x3, 2*y3, 1)
) * Transpose(xc, yc, K)
or
B = A * s

The solution is
          | b1 2*y1 1 |
          | b2 2*y2 1 |
          | b3 2*y3 1 |
xc =   --------------
         | 2*x1 2*y1 1 |
         | 2*x2 2*y2 1 |
         | 2*x3 2*y3 1 |

          | 2*x1 b1 1 |
          | 2*x2 b2 1 |
          | 2*x3 b3 1 |
yc =   --------------
         | 2*x1 2*y1 1 |
         | 2*x2 2*y2 1 |
         | 2*x3 2*y3 1 |

         | 2*x1 2*y1 b1 |
         | 2*x2 2*y2 b2 |
         | 2*x3 2*y3 b3 |
K =     --------------
         | 2*x1 2*y1 1 |
         | 2*x2 2*y2 1 |
         | 2*x3 2*y3 1 |

| 2*x1 2*y1 1 |
| 2*x2 2*y2 1 | = 4*( x1*(y2-y3) - y1*(x2-x3) + (x2*y3 - x3*y2) ) = detA
| 2*x3 2*y3 1 |

So the final solution is :

xc = 2*( b1*(y2-y3) - y1*(b2-b3) + (b2*y3-b3*y2) ) / detA
yc = 2*( x1*(b2-b3) - b1*(x2-x3) + (x2*b3-x3*b2) ) / detA
K = 4*( x1*(y2*b3-y3*b2) - y1*(x2*b3-x3*b2) + b1*(x2*y3-x3*y2) )/detA
r = sqrt(K + sqr(x) + sqr(y))

... or simply use the
Best regards,
Peter


On 13.02.2016 18:59, Christopher Coulon wrote:

> Thank you Jerome!  Last night, as I was waiting to go to sleep, I got the inspiration that I could do this by looking through the entire image for the pixel location equidistant from all three points, so I wrote the macro below to test the theory, and it worked!  :-)
>
> newImage("Untitled", "8-bit white", 1000, 1000, 1);
> run("Colors...", "foreground=black background=white selection=yellow");
> run("Set Measurements...", "area centroid limit redirect=None decimal=8");
> makeLine(234, 156, 240, 156);
> run("Draw", "slice");
> makeLine(369, 120, 372, 120);
> run("Draw", "slice");
> makeLine(466, 180, 470, 182);
> run("Draw", "slice");
> setThreshold(0, 128);
> run("Analyze Particles...", "size=2-Infinity display exclude clear record");
> run("Select None");
> resetThreshold();
>
> x0 = getResult("X", 0);
> y0 = getResult("Y", 0);
> x1 = getResult("X", 1);
> y1 = getResult("Y", 1);
> x2 = getResult("X", 2);
> y2 = getResult("Y", 2);
> l0 = 0;
> l1 = 0;
> l2 = 0;
> top = getHeight();
> left = getWidth();
>
> center = newArray(2);;
> radius = 0;
> go = true;
>
> for(x=0; x<getWidth(); x++) {
> for(y = 0; y < getHeight(); y++) {
> l1 = round(sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0)));
> l2 = round(sqrt((x-x1)*(x-x1)+(y-y1)*(y-y1)));
> l3 = round(sqrt((x-x2)*(x-x2)+(y-y2)*(y-y2)));
> if(go && l1 == l2 && l1 == l3) {
> center[0] = x;
> center[1] = y;
> radius = l1;
> if(y > 0) go = false;
> }
> }
> }
>
> print("center x, y = " + center[0] + ", " + center[1]);
> print("radius = " + radius);
>
> left = center[0] - radius;
> top = center[1] - radius;
> diam = radius * 2;
> drawOval(left, top, diam, diam);
> print("left = " + left + "  top = " + top);
>
>> On Feb 13, 2016, at 1:20 AM, Jerome Mutterer <[hidden email]> wrote:
>>
>> Hi Christopher,
>> The "Fit Circle" command does this. It's code is at:
>> https://github.com/imagej/imagej1/blob/master/ij/plugin/Selection.java#L103
>> Sincerely,
>> Jerome.
>>
>>
>> On 13 February 2016 at 04:30, Christopher Coulon <
>> [hidden email]> wrote:
>>> Before I reinvent the wheel, is there anyone who can share the code for
>> drawing a circle from three points?
>>> Thanks in advance,
>>>
>>> Chris
>>> --
>>> ImageJ mailing list: http://imagej.nih.gov/ij/list.html
>>
>>
>>
>> --
>> Jerome Mutterer
>> CNRS - Institut de biologie moléculaire des plantes
>> 12, rue du Général Zimmer
>> 67084 Strasbourg Cedex
>> www.ibmp.cnrs.fr
>>
>> --
>> ImageJ mailing list: http://imagej.nih.gov/ij/list.html
>
> --
> ImageJ mailing list: http://imagej.nih.gov/ij/list.html
>

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html