Radius, Center of Circle Given 3 Points

Date: 8/28/96 at 17:41:17
From: Stuart Simmons
Subject: Radius, Center of Circle Given 3 Points

I'm writing a small piece of software for personal use and I have come
across a problem. 

I need an easy formula to calculate the centre point and radius of a 
circle given three points around the circumference of the circle. 
(x1,y1),(x2,y2),(x3,y3)

Regards
Stuart

Date: 8/30/96 at 13:27:38
From: Doctor Tom
Subject: Re: Radius, Center of Circle Given 3 Points

Here's how I'd do it:

The perpendicular bisector of the segment between two points is the
set of all points equidistant from both.  So if you take the
perpendicular bisector of (x1,y1) and (x2,y2) and the perpendicular
bisector of the segment from (x2,y2) to (x3,y3) and find the
intersection of those lines, that point will be the center.  

To find the radius, find the distance from that point to any of the
three points.  To find the equation of the perpendicular bisector
of (x1,y1) to (x2,y2), you know that it passes through the midpoint of 
the segment: ((x1+x2)/2,(y1+y2)/2), and if the slope of the line 
connecting (x1,y1) to (x2,y2) is m, the slope of the perpendicular 
bisector is -1/m.  Work out the equations for the two lines, find 
their intersection, and bingo!  You've got the coordinates of the 
center.

An error should occur if the three points lie on a line, and you'll
need special code to check for the case where one of the slopes is
zero.

Otherwise it's pretty straightforward.

By a miracle, I wrote some code (in C++) to do something like this
once.  I'll include it, and hope you can make some use of it.  Don't
worry about the "w" coordinates; assume they're all 1.0.  I wanted
my code to work in the perspective case and I had to worry about
"points at infinity":

void circle_vvv(circle *c)
{
    c->center.w = 1.0;
    vertex *v1 = (vertex *)c->c.p1;
    vertex *v2 = (vertex *)c->c.p2;
    vertex *v3 = (vertex *)c->c.p3;
    float bx = v1->xw; float by = v1->yw;
    float cx = v2->xw; float cy = v2->yw;
    float dx = v3->xw; float dy = v3->yw;
    float temp = cx*cx+cy*cy;
    float bc = (bx*bx + by*by - temp)/2.0;
    float cd = (temp - dx*dx - dy*dy)/2.0;
    float det = (bx-cx)*(cy-dy)-(cx-dx)*(by-cy);
    if (fabs(det) < 1.0e-6) {
        c->center.xw = c->center.yw = 1.0;
        c->center.w = 0.0;
        c->v1 = *v1;
        c->v2 = *v2;
        c->v3 = *v3;
        return;
        }
    det = 1/det;
    c->center.xw = (bc*(cy-dy)-cd*(by-cy))*det;
    c->center.yw = ((bx-cx)*cd-(cx-dx)*bc)*det;
    cx = c->center.xw; cy = c->center.yw;
    c->radius = sqrt((cx-bx)*(cx-bx)+(cy-by)*(cy-by));
}

-Doctor Tom,  The Math Forum
 Check out our web site!  http://mathforum.org/dr.math/