Точки пересечения окружности - PullRequest
49 голосов
/ 28 июля 2010

Как рассчитать точки пересечения двух окружностей.Я ожидал бы, что во всех случаях будут две, одна или ни одна точка пересечения.

У меня есть координаты x и y центральной точки и радиус для каждого круга.

Ответ на python был бы предпочтительнее, но любой рабочий алгоритм был бы приемлемым.

Ответы [ 5 ]

85 голосов
/ 28 июля 2010

Пересечение двух окружностей

Автор Пол Бурк

В следующем примечании описывается, как найти точку (точки) пересечения двух окружностей наплоскость, используются следующие обозначения.Цель состоит в том, чтобы найти две точки P 3 = (x 3 , y 3 ), если они существуют.

Intersection of 2 circles

Сначала вычислите расстояние d между центром окружностей.d = || P 1 - P 0 ||.

  • Если d> r 0 + r 1 , то решений нет, кружки разделены.

  • Если d <| r <sub>0 - r 1 |тогда нет решений, потому что один круг содержится внутри другого.

  • Если d = 0 и r 0 = r 1 , то круги совпадают и существуютбесконечное число решений.

С учетом двух треугольников P 0 P 2 P 3 и P 1 P 2 P 3 мы можем написать

a 2 + h 2 = r 0 2 и b 2 + h 2 = r 1 2

Использование d =a + b мы можем найти для a,

a = (r 0 2 - r 1 2 + d 2 ) / (2 d)

Легко показать, что это уменьшается до r
0 , когда два круга касаются одной точки, то есть: d =r 0 + r 1 Решите для h, подставив в первое уравнение h 2 = r 0 2 - a 2

Итак

P 2 = P 0 + a (P 1 - P 0 ) / d

И, наконец, P 3 = (x 3 , y 3 ) вусловия P 0 = (x 0 , y 0 ), P 1 = (x 1 ,y 1 ) и P 2 = (x 2 , y 2 ), составляет

x 3 = x 2 + - ч (у 1 - у 0 ) / д

у 3 = у 2 - + ч (х 1 - х 0 ) / д

Источник: http://paulbourke.net/geometry/circlesphere/

15 голосов
/ 03 января 2013

Вот моя реализация C ++, основанная на статье Пола Бурка .Он работает только при наличии двух пересечений, в противном случае он, вероятно, возвращает NaN NAN NAN NAN.

class Point{
    public:
        float x, y;
        Point(float px, float py) {
            x = px;
            y = py;
        }
        Point sub(Point p2) {
            return Point(x - p2.x, y - p2.y);
        }
        Point add(Point p2) {
            return Point(x + p2.x, y + p2.y);
        }
        float distance(Point p2) {
            return sqrt((x - p2.x)*(x - p2.x) + (y - p2.y)*(y - p2.y));
        }
        Point normal() {
            float length = sqrt(x*x + y*y);
            return Point(x/length, y/length);
        }
        Point scale(float s) {
            return Point(x*s, y*s);
        }
};

class Circle {
    public:
        float x, y, r, left;
        Circle(float cx, float cy, float cr) {
            x = cx;
            y = cy;
            r = cr;
            left = x - r;
        }
        pair<Point, Point> intersections(Circle c) {
            Point P0(x, y);
            Point P1(c.x, c.y);
            float d, a, h;
            d = P0.distance(P1);
            a = (r*r - c.r*c.r + d*d)/(2*d);
            h = sqrt(r*r - a*a);
            Point P2 = P1.sub(P0).scale(a/d).add(P0);
            float x3, y3, x4, y4;
            x3 = P2.x + h*(P1.y - P0.y)/d;
            y3 = P2.y - h*(P1.x - P0.x)/d;
            x4 = P2.x - h*(P1.y - P0.y)/d;
            y4 = P2.y + h*(P1.x - P0.x)/d;

            return pair<Point, Point>(Point(x3, y3), Point(x4, y4));
        }

};
8 голосов
/ 06 июля 2017

Вот реализация в Javascript с использованием векторов.Код хорошо документирован, вы должны быть в состоянии следовать ему.Вот исходный код

Смотрите демонстрацию здесь : enter image description here

// Let EPS (epsilon) be a small value
var EPS = 0.0000001;

// Let a point be a pair: (x, y)
function Point(x, y) {
  this.x = x;
  this.y = y;
}

// Define a circle centered at (x,y) with radius r
function Circle(x,y,r) {
  this.x = x;
  this.y = y;
  this.r = r;
}

// Due to double rounding precision the value passed into the Math.acos
// function may be outside its domain of [-1, +1] which would return
// the value NaN which we do not want.
function acossafe(x) {
  if (x >= +1.0) return 0;
  if (x <= -1.0) return Math.PI;
  return Math.acos(x);
}

// Rotates a point about a fixed point at some angle 'a'
function rotatePoint(fp, pt, a) {
  var x = pt.x - fp.x;
  var y = pt.y - fp.y;
  var xRot = x * Math.cos(a) + y * Math.sin(a);
  var yRot = y * Math.cos(a) - x * Math.sin(a);
  return new Point(fp.x+xRot,fp.y+yRot);
}

// Given two circles this method finds the intersection
// point(s) of the two circles (if any exists)
function circleCircleIntersectionPoints(c1, c2) {

  var r, R, d, dx, dy, cx, cy, Cx, Cy;

  if (c1.r < c2.r) {
    r  = c1.r;  R = c2.r;
    cx = c1.x; cy = c1.y;
    Cx = c2.x; Cy = c2.y;
  } else {
    r  = c2.r; R  = c1.r;
    Cx = c1.x; Cy = c1.y;
    cx = c2.x; cy = c2.y;
  }

  // Compute the vector <dx, dy>
  dx = cx - Cx;
  dy = cy - Cy;

  // Find the distance between two points.
  d = Math.sqrt( dx*dx + dy*dy );

  // There are an infinite number of solutions
  // Seems appropriate to also return null
  if (d < EPS && Math.abs(R-r) < EPS) return [];

  // No intersection (circles centered at the 
  // same place with different size)
  else if (d < EPS) return [];

  var x = (dx / d) * R + Cx;
  var y = (dy / d) * R + Cy;
  var P = new Point(x, y);

  // Single intersection (kissing circles)
  if (Math.abs((R+r)-d) < EPS || Math.abs(R-(r+d)) < EPS) return [P];

  // No intersection. Either the small circle contained within 
  // big circle or circles are simply disjoint.
  if ( (d+r) < R || (R+r < d) ) return [];

  var C = new Point(Cx, Cy);
  var angle = acossafe((r*r-d*d-R*R)/(-2.0*d*R));
  var pt1 = rotatePoint(C, P, +angle);
  var pt2 = rotatePoint(C, P, -angle);
  return [pt1, pt2];

}
6 голосов
/ 15 марта 2017

Почему бы просто не использовать 7 строк вашего любимого процедурного языка (или программируемого калькулятора!), Как показано ниже.

При условии, что вам заданы координаты P0 (x0, y0), координаты P1 (x1, y1), r0 и r1, и вы хотите найти координаты P3 (x3, y3):

d=sqr((x1-x0)^2 + (y1-y0)^2)
a=(r0^2-r1^2+d^2)/(2*d)
h=sqr(r0^2-a^2)
x2=x0+a*(x1-x0)/d   
y2=y0+a*(y1-y0)/d   
x3=x2+h*(y1-y0)/d       // also x3=x2-h*(y1-y0)/d
y3=y2-h*(x1-x0)/d       // also y3=y2+h*(x1-x0)/d
0 голосов
/ 30 мая 2018

Попробуйте это;

    def ri(cr1,cr2,cp1,cp2):
        int1=[]
        int2=[]
        ori=0
        if cp1[0]<cp2[0] and cp1[1]!=cp2[1]:
            p1=cp1
            p2=cp2
            r1=cr1
            r2=cr2
            if cp1[1]<cp2[1]:
                ori+=1
            elif cp1[1]>cp2[1]:
                ori+=2        
        elif cp1[0]>cp2[0] and cp1[1]!=cp2[1]:
            p1=cp2
            p2=cp1
            r1=cr2
            r2=cr1
            if p1[1]<p2[1]:
                ori+=1
            elif p1[1]>p2[1]:
                ori+=2
        elif cp1[0]==cp2[0]:
            ori+=4
            if cp1[1]>cp2[1]:
                p1=cp1
                p2=cp2
                r1=cr1
                r2=cr2
            elif cp1[1]<cp2[1]:
                p1=cp2
                p2=cp1
                r1=cr2
                r2=cr1
        elif cp1[1]==cp2[1]:
            ori+=3
            if cp1[0]>cp2[0]:
                p1=cp2
                p2=cp1
                r1=cr2
                r2=cr1
            elif cp1[0]<cp2[0]:
                p1=cp1
                p2=cp2
                r1=cr1
                r2=cr2
        if ori==1:#+
            D=calc_dist(p1,p2)
            tr=r1+r2
            el=tr-D
            a=r1-el
            b=r2-el
            A=a+(el/2)
            B=b+(el/2)
            thta=math.degrees(math.acos(A/r1))
            rs=p2[1]-p1[1]
            rn=p2[0]-p1[0]
            gd=rs/rn
            yint=p1[1]-((gd)*p1[0])
            dty=calc_dist(p1,[0,yint])

            aa=p1[1]-yint
            bb=math.degrees(math.asin(aa/dty))
            d=90-bb
            e=180-d-thta
            g=(dty/math.sin(math.radians(e)))*math.sin(math.radians(thta))
            f=(g/math.sin(math.radians(thta)))*math.sin(math.radians(d))
            oty=yint+g
            h=f+r1
            i=90-e
            j=180-90-i
            l=math.sin(math.radians(i))*h
            k=math.cos(math.radians(i))*h
            iy2=oty-l
            ix2=k
            int2.append(ix2)
            int2.append(iy2)

            m=90+bb
            n=180-m-thta
            p=(dty/math.sin(math.radians(n)))*math.sin(math.radians(m))
            o=(p/math.sin(math.radians(m)))*math.sin(math.radians(thta))
            q=p+r1
            r=90-n
            s=math.sin(math.radians(r))*q
            t=math.cos(math.radians(r))*q
            otty=yint-o
            iy1=otty+s
            ix1=t
            int1.append(ix1)
            int1.append(iy1)
        elif ori==2:#-
            D=calc_dist(p1,p2)
            tr=r1+r2
            el=tr-D
            a=r1-el
            b=r2-el
            A=a+(el/2)
            B=b+(el/2)
            thta=math.degrees(math.acos(A/r1))
            rs=p2[1]-p1[1]
            rn=p2[0]-p1[0]
            gd=rs/rn
            yint=p1[1]-((gd)*p1[0])
            dty=calc_dist(p1,[0,yint])

            aa=yint-p1[1]
            bb=math.degrees(math.asin(aa/dty))
            c=180-90-bb
            d=180-c-thta
            e=180-90-d
            f=math.tan(math.radians(e))*p1[0]
            g=math.sqrt(p1[0]**2+f**2)
            h=g+r1
            i=180-90-e
            j=math.sin(math.radians(e))*h
            jj=math.cos(math.radians(i))*h
            k=math.cos(math.radians(e))*h
            kk=math.sin(math.radians(i))*h
            l=90-bb
            m=90-e
            tt=l+m+thta
            n=(dty/math.sin(math.radians(m)))*math.sin(math.radians(thta))
            nn=(g/math.sin(math.radians(l)))*math.sin(math.radians(thta))
            oty=yint-n
            iy1=oty+j
            ix1=k
            int1.append(ix1)
            int1.append(iy1)

            o=bb+90
            p=180-o-thta
            q=90-p
            r=180-90-q
            s=(dty/math.sin(math.radians(p)))*math.sin(math.radians(o))
            t=(s/math.sin(math.radians(o)))*math.sin(math.radians(thta))
            u=s+r1
            v=math.sin(math.radians(r))*u
            vv=math.cos(math.radians(q))*u
            w=math.cos(math.radians(r))*u
            ww=math.sin(math.radians(q))*u
            ix2=v
            otty=yint+t
            iy2=otty-w
            int2.append(ix2)
            int2.append(iy2)

        elif ori==3:#y
            D=calc_dist(p1,p2)
            tr=r1+r2
            el=tr-D
            a=r1-el
            b=r2-el
            A=a+(el/2)
            B=b+(el/2)
            b=math.sqrt(r1**2-A**2)
            int1.append(p1[0]+A)
            int1.append(p1[1]+b)
            int2.append(p1[0]+A)
            int2.append(p1[1]-b)
        elif ori==4:#x
            D=calc_dist(p1,p2)
            tr=r1+r2
            el=tr-D
            a=r1-el
            b=r2-el
            A=a+(el/2)
            B=b+(el/2)
            b=math.sqrt(r1**2-A**2)
            int1.append(p1[0]+b)
            int1.append(p1[1]-A)
            int2.append(p1[0]-b)
            int2.append(p1[1]-A)
        return [int1,int2]
    def calc_dist(p1,p2):
        return math.sqrt((p2[0] - p1[0]) ** 2 +
                         (p2[1] - p1[1]) ** 2)
...