Каков алгоритм нахождения центра окружности по трем точкам?

У меня есть три точки на окружности круга:

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

Как вычислить центр круга?

Реализация в Processing (Java).

Я нашел ответ и реализовал рабочее решение:

 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;
  }

person Russell Strauss    schedule 05.11.2010    source источник
comment
mathforum.org/library/drmath/view/55233.html   -  person John C Earls    schedule 05.11.2010
comment
Большое спасибо за публикацию ответа на ваш вопрос.   -  person SpoiledTechie.com    schedule 15.04.2011
comment
Ваше решение дает странные результаты, если B.x - A.x или C.x - B.x равно нулю, потому что вы затем делите на ноль.   -  person gexicide    schedule 14.06.2015


Ответы (6)


Это может быть довольно глубокий расчет. Здесь есть простая пошаговая инструкция: http://paulbourke.net/geometry/circlesphere/. Когда у вас есть уравнение окружности, вы можете просто представить его в форме, включающей H и K. Точка (h,k) будет центром.

(прокрутите немного вниз по ссылке, чтобы перейти к уравнениям)

person Azmisov    schedule 05.11.2010
comment
Это страница, которая привела меня к ответу. Я реализовал это сам как: - person Russell Strauss; 05.11.2010
comment
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 центр = P (0,0); float aSlope = yDelta_a/xDelta_a; float bSlope = yDelta_b/xDelta_b; center.x = (aSlopebSlope*(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; центр возврата; } - person Russell Strauss; 05.11.2010

Вот мой порт Java, уклоняющийся от ошибки, когда определитель исчезает с очень элегантным IllegalArgumentException, мой подход к работе с условиями «точки находятся на расстоянии двух друг от друга» или «точки лежат на линии». Кроме того, это вычисляет радиус (и справляется с исключительными условиями), чего не будет делать ваш подход с пересекающимися склонами.

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

}

См. алгоритм здесь:

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));
}
person andersoj    schedule 05.11.2010
comment
Мне трудно увидеть, какие 3 вершины являются исходными. v1, v2 и v3? - person Russell Strauss; 05.11.2010
comment
Да, это не лучший код; Я записал это. v1,2,3 — исходные вершины. (bx,by), (cx,cy), (dx,dy) — координаты. - person andersoj; 05.11.2010
comment
@Russell Strauss: я предоставил Java-порт этого кода, что делает поток более понятным. - person andersoj; 05.11.2010

Я искал подобный алгоритм, когда зависал над этим вопросом. Взял ваш код, но обнаружил, что это не будет работать в случаях, когда любой из наклонов равен 0 или бесконечности (может быть правдой, когда xDelta_a или xDelta_b равно 0).

Я поправил алгоритм и вот мой код. Примечание. Я использовал язык программирования Objective-C и просто меняю код для инициализации значения точки, поэтому, если это неправильно, я уверен, что программист, работающий на Java, может это исправить. Логика, однако, у всех одинаковая (да благословит Бог алгоритмы!! :))

Работает отлично, что касается моего собственного функционального тестирования. Пожалуйста, дайте мне знать, если логика неверна в любой момент.

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;
}
person IamRKhanna    schedule 10.10.2013

Мне жаль, что мой ответ запоздал. Любое решение, использующее «наклон», потерпит неудачу, если две точки образуют вертикальную линию, потому что наклон будет бесконечным.

Вот простое надежное решение для 2019 года, которое всегда работает правильно:

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;
}

Приведенный выше код вернет «false» тогда и только тогда, когда 3 точки лежат на одной прямой.

person Adam Gawne-Cain    schedule 08.08.2019
comment
Почему ux задается с помощью y (p1[1], p2[1]), а uy задается с помощью x (p1[0], p2[0])? Тот же вопрос для vx и vy. - person trans; 13.12.2019
comment
Центр окружности находится там, где пересекаются ортогональные биссектрисы прямых (p1,p2) и (p2,p3). Ортогональная биссектриса прямой (p1,p2) проходит через (ax,y) и имеет направление (ux,uy). Ортогональная биссектриса прямой (p2,p3) проходит через (bx,by) и имеет направление (vx,vy). Чтобы ответить на вопрос, поменяв местами X и Y и изменив знак единицы, вы получите ортогональный вектор. - person Adam Gawne-Cain; 15.12.2019
comment
Я понимаю. Спасибо! Кстати, ваше решение оказалось очень кстати и сэкономило мне кучу времени. Единственное, что может сделать это решение лучше, — это диаграмма. - person trans; 16.12.2019

person    schedule
comment
Опечатка: mb * (p1.x - p2.x) должно быть mb * (p1.x + p2.x) - person tim_hutton; 15.01.2016
comment
Это решение не работает с zero_divide, если две точки имеют одинаковую координату X. - person Adam Gawne-Cain; 08.08.2019
comment
Да, вы правы, это частный случай, который нужно учитывать - person Ernesto Alfonso; 08.08.2019

person    schedule
comment
использовать середины всех точек, данных на окружности. совпадет только одна пара, и это будет середина круга. - person Nitin Jindal; 26.05.2020