36

I have three points on the circumference of a circle:

pt A = (A.x, A.y);
pt B = (B.x, B.y);
pt C = (C.x, C.y);

How do I calculate the center of the circle?

Implementing it in Processing (Java).

I found the answer and implemented a working solution:

 pt circleCenter(pt A, pt B, pt C) {

    float yDelta_a = B.y - A.y;
    float xDelta_a = B.x - A.x;
    float yDelta_b = C.y - B.y;
    float xDelta_b = C.x - B.x;
    pt center = P(0,0);

    float aSlope = yDelta_a/xDelta_a;
    float bSlope = yDelta_b/xDelta_b;  
    center.x = (aSlope*bSlope*(A.y - C.y) + bSlope*(A.x + B.x)
        - aSlope*(B.x+C.x) )/(2* (bSlope-aSlope) );
    center.y = -1*(center.x - (A.x+B.x)/2)/aSlope +  (A.y+B.y)/2;

    return center;
  }
Mark Fisher
  • 965
  • 1
  • 11
  • 30
Russell Strauss
  • 1,160
  • 3
  • 14
  • 30

6 Answers6

21

Here's my Java port, dodging the error condition when the determinant disappears with a very elegant IllegalArgumentException, my approach to coping with the "points are two far apart" or "points lie on a line" conditions. Also, this computes the radius (and copes with exceptional conditions) which your intersecting-slopes approach will not do.

public class CircleThree
{ 
  static final double TOL = 0.0000001;
  
  public static Circle circleFromPoints(final Point p1, final Point p2, final Point p3)
  {
    final double offset = Math.pow(p2.x,2) + Math.pow(p2.y,2);
    final double bc =   ( Math.pow(p1.x,2) + Math.pow(p1.y,2) - offset )/2.0;
    final double cd =   (offset - Math.pow(p3.x, 2) - Math.pow(p3.y, 2))/2.0;
    final double det =  (p1.x - p2.x) * (p2.y - p3.y) - (p2.x - p3.x)* (p1.y - p2.y); 
    
    if (Math.abs(det) < TOL) { throw new IllegalArgumentException("Yeah, lazy."); }

    final double idet = 1/det;
     
    final double centerx =  (bc * (p2.y - p3.y) - cd * (p1.y - p2.y)) * idet;
    final double centery =  (cd * (p1.x - p2.x) - bc * (p2.x - p3.x)) * idet;
    final double radius = 
       Math.sqrt( Math.pow(p2.x - centerx,2) + Math.pow(p2.y-centery,2));
    
    return new Circle(new Point(centerx,centery),radius);
  }
  
  static class Circle
  {
    final Point center;
    final double radius;
    public Circle(Point center, double radius)
    {
      this.center = center; this.radius = radius;
    }
    @Override 
    public String toString()
    {
      return new StringBuilder().append("Center= ").append(center).append(", r=").append(radius).toString();
    }
  }
  
  static class Point
  {
    final double x,y;

    public Point(double x, double y)
    {
      this.x = x; this.y = y;
    }
    @Override
    public String toString()
    {
      return "("+x+","+y+")";
    }
    
  }

  public static void main(String[] args)
  {
    Point p1 = new Point(0.0,1.0);
    Point p2 = new Point(1.0,0.0);
    Point p3 = new Point(2.0,1.0);
    Circle c = circleFromPoints(p1, p2, p3);
    System.out.println(c);
  }
  
}

See algorithm from The Math Forum:

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));
}
Laurel
  • 5,965
  • 14
  • 31
  • 57
andersoj
  • 22,406
  • 7
  • 62
  • 73
17

It can be a rather in depth calculation. There is a simple step-by-step here: http://paulbourke.net/geometry/circlesphere/. Once you have the equation of the circle, you can simply put it in a form involving H and K. The point (h,k) will be the center.

(scroll down a little ways at the link to get to the equations)

Azmisov
  • 6,493
  • 7
  • 53
  • 70
  • This is the page that led me to the answer. I implemented it myself as: – Russell Strauss Nov 05 '10 at 04:09
  • 1
    pt circleCenter(pt A, pt B, pt C) { float yDelta_a = B.y - A.y; float xDelta_a = B.x - A.x; float yDelta_b = C.y - B.y; float xDelta_b = C.x - B.x; pt center = P(0,0); float aSlope = yDelta_a/xDelta_a; float bSlope = yDelta_b/xDelta_b; center.x = (aSlope*bSlope*(A.y - C.y) + bSlope*(A.x + B.x) - aSlope*(B.x+C.x) )/(2* (bSlope-aSlope) ); center.y = -1*(center.x - (A.x+B.x)/2)/aSlope + (A.y+B.y)/2; return center; } – Russell Strauss Nov 05 '10 at 04:10
7

I was looking for a similar algorithm when I hovered over this question. Took your code but found that this will not work in cases when where either of the slope is 0 or infinity (can be true when xDelta_a or xDelta_b is 0).

I corrected the algorithm and here is my code. Note: I used objective-c programming language and am just changing the code for point value initialization, so if that is wrong, I am sure programmer working in java can correct it. The logic, however, is the same for all (God bless algorithms!! :))

Works perfectly fine as far as my own functional testing is concerned. Please let me know if logic is wrong at any point.

pt circleCenter(pt A, pt B, pt C) {

float yDelta_a = B.y - A.y;
float xDelta_a = B.x - A.x;
float yDelta_b = C.y - B.y;
float xDelta_b = C.x - B.x;
pt center = P(0,0);

float aSlope = yDelta_a/xDelta_a;
float bSlope = yDelta_b/xDelta_b;

pt AB_Mid = P((A.x+B.x)/2, (A.y+B.y)/2);
pt BC_Mid = P((B.x+C.x)/2, (B.y+C.y)/2);

if(yDelta_a == 0)         //aSlope == 0
{
    center.x = AB_Mid.x;
    if (xDelta_b == 0)         //bSlope == INFINITY
    {
        center.y = BC_Mid.y;
    }
    else
    {
        center.y = BC_Mid.y + (BC_Mid.x-center.x)/bSlope;
    }
}
else if (yDelta_b == 0)               //bSlope == 0
{
    center.x = BC_Mid.x;
    if (xDelta_a == 0)             //aSlope == INFINITY
    {
        center.y = AB_Mid.y;
    }
    else
    {
        center.y = AB_Mid.y + (AB_Mid.x-center.x)/aSlope;
    }
}
else if (xDelta_a == 0)        //aSlope == INFINITY
{
    center.y = AB_Mid.y;
    center.x = bSlope*(BC_Mid.y-center.y) + BC_Mid.x;
}
else if (xDelta_b == 0)        //bSlope == INFINITY
{
    center.y = BC_Mid.y;
    center.x = aSlope*(AB_Mid.y-center.y) + AB_Mid.x;
}
else
{
    center.x = (aSlope*bSlope*(AB_Mid.y-BC_Mid.y) - aSlope*BC_Mid.x + bSlope*AB_Mid.x)/(bSlope-aSlope);
    center.y = AB_Mid.y - (center.x - AB_Mid.x)/aSlope;
}

return center;
}
IamRKhanna
  • 228
  • 2
  • 10
3

I am sorry my answer is late. Any solution using "slope" will fail when two of the points form a vertical line, because the slope will be infinite.

Here is a simple robust solution for 2019 that always works correctly:

public static boolean circleCenter(double[] p1, double[] p2, double[] p3, double[] center) {
    double ax = (p1[0] + p2[0]) / 2;
    double ay = (p1[1] + p2[1]) / 2;
    double ux = (p1[1] - p2[1]);
    double uy = (p2[0] - p1[0]);
    double bx = (p2[0] + p3[0]) / 2;
    double by = (p2[1] + p3[1]) / 2;
    double vx = (p2[1] - p3[1]);
    double vy = (p3[0] - p2[0]);
    double dx = ax - bx;
    double dy = ay - by;
    double vu = vx * uy - vy * ux;
    if (vu == 0)
        return false; // Points are collinear, so no unique solution
    double g = (dx * uy - dy * ux) / vu;
    center[0] = bx + g * vx;
    center[1] = by + g * vy;
    return true;
}

The above code will return "false" if and only if the 3 points are collinear.

Adam Gawne-Cain
  • 1,347
  • 14
  • 14
  • Why is `ux` set using "y" (`p1[1]`, `p2[1]`) and `uy` set using "x" (`p1[0]`, `p2[0]`)? Same question for `vx` and `vy`. – trans Dec 13 '19 at 17:57
  • The center of circle is where orthogonal bisectors of lines (p1,p2) and (p2,p3) intersect. The orthogonal bisector of line (p1,p2) passes though (ax, ay) and has direction (ux,uy). The orthogonal bisector of line (p2,p3) passes though (bx, by) and has direction (vx,vy). To answer question, swapping X and Y and changing sign of one gives an orthogonal vector. – Adam Gawne-Cain Dec 14 '19 at 21:47
  • I see. Thanks! By the way, your solution came in very handy and saved me a bunch of time. Only thing that could make this solution any better would be a diagram. – trans Dec 16 '19 at 20:53
2
public Vector2 CarculateCircleCenter(Vector2 p1, Vector2 p2, Vector2 p3)
{
    if (
        p2.x - p1.x == 0 ||
        p3.x - p2.x == 0 ||
        p2.y - p1.y == 0 ||
        p3.y - p2.y == 0
    ) return null;

    Vector2 center = new Vector2();
    float ma = (p2.y - p1.y) / (p2.x - p1.x);
    float mb = (p3.y - p2.y) / (p3.x - p2.x);
    center.x = (ma * mb * (p1.y - p3.y) + mb * (p1.x - p2.x) - ma * (p2.x + p3.x)) / (2 * (mb - ma));
    center.y = (-1 / ma) * (center.x - (p1.x + p2.x) * 0.5) + (p1.y + p2.y) * 0.5;
    return center;
}
Ernesto Alfonso
  • 650
  • 10
  • 30
0
def circle_mid_point(list_b):
list_k = []
for i in range(3):
    for j in range(1,4):
        if i+j <=3:
            midpoint_x1 =  (list_b[i][0] + list_b[i+j][0])/2
            midpoint_y1 = (list_b[i][1] + list_b[i + j][1]) / 2
            list_k.append([midpoint_x1,midpoint_y1]) # list of all the midpoints of the lines

for k in range(len(list_k)):
    for j in range(1,len(list_k)-1):
        if list_k[k] == list_k[k+j]: #at centre only two midpoints will have the same value
            return list_k[k]

k = circle_mid_point([[0,1],[1,0],[-1,0],[0,-1]])
print(k)
deadshot
  • 8,881
  • 4
  • 20
  • 39
  • use midpoints of all the points given on the circumference. only single pair will match and that is the midpoint of the circle. – Nitin Jindal May 26 '20 at 05:54