algorithm - Puntos de intersección círculo-círculo
math geometry (5)
Intersección de dos círculos
Escrito por Paul Bourke
La siguiente nota describe cómo encontrar los puntos de intersección entre dos círculos en un plano, se usa la siguiente notación. El objetivo es encontrar los dos puntos P 3 = (x 3 , y 3 ) si existen.
Primero calcule la distancia d entre el centro de los círculos. d = || P 1 - P 0 ||.
- Si d> r 0 + r 1, entonces no hay soluciones, los círculos están separados.
- Si d <| r 0 - r 1 | entonces no hay soluciones porque un círculo está contenido dentro del otro.
- Si d = 0 y r 0 = r 1, entonces los círculos son coincidentes y hay un número infinito de soluciones.
Teniendo en cuenta los dos triángulos P 0 P 2 P 3 y P 1 P 2 P 3 podemos escribir
a 2 + h 2 = r 0 2 y b 2 + h 2 = r 1 2
Usando d = a + b podemos resolver por a,
a = (r 0 2 - r 1 2 + d 2 ) / (2 d)
Se puede demostrar fácilmente que esto se reduce a r 0 cuando los dos círculos se tocan en un punto, es decir: d = r 0 + r 1 Resuelve para h sustituyendo a en la primera ecuación, h 2 = r 0 2 - a 2Asi que
P 2 = P 0 + a (P 1 - P 0 ) / d
Y finalmente, P 3 = (x 3 , y 3 ) en términos de P 0 = (x 0 , y 0 ), P 1 = (x 1 , y 1 ) y P 2 = (x 2 , y 2 ), es
x 3 = x 2 + - h (y 1 - y 0 ) / d
y 3 = y 2 - + h (x 1 - x 0 ) / d
¿Cómo puedo calcular los puntos de intersección de dos círculos? Esperaría que hubiera dos, uno o ningún punto de intersección en todos los casos.
Tengo las coordenadas xey del punto central y el radio de cada círculo.
Se preferiría una respuesta en python, pero cualquier algoritmo de trabajo sería aceptable.
¿Por qué no usar 7 líneas de su lenguaje de procedimiento favorito (o calculadora programable) como se muestra a continuación?
Suponiendo que se le otorgan coordenadas P0 (x0, y0), colas P1 (x1, y1), r0 y r1 y quiere encontrar colas 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
Aquí está mi implementación de C ++ basada en el artículo de Paul Bourke . Solo funciona si hay dos intersecciones, de lo contrario probablemente devuelva 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));
}
};
Aquí hay una implementación en JavaScript usando vectores. El código está bien documentado, debería poder seguirlo. Aquí está la fuente original
Vea la demostración en vivo 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];
}
Prueba esto;
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)