查找兩個球體的交點

Question

考慮以下問題 - 我給出了長度爲L0和L1的2個鏈接。 P0是第一個鏈接開始的點，P1是我希望第二個鏈接的結束位於三維空間的點。我應該編寫一個函數，將這些三維點（P0和P1）作爲輸入，並且應該找到將第二個鏈接的終點放在P1處的所有鏈接配置。

我對如何去做的理解是 - 每個鏈接L0和L1將圍繞它自己創建一個球體S0和S1。我應該找出這兩個球體（這將是一個圓圈）的交點，並打印出該圓周上的所有點。

我看到gmatt在Finding intersection points between 3 spheres上的第一個回覆，但由於圖像沒有顯示，所以無法正確理解。我也看到了一個公式，找出在http://mathworld.wolfram.com/Sphere-SphereIntersection.html

十字路口我可以找到數學世界給出的方法的交點半徑。我也可以找到該圓的中心，然後使用圓的參數方程來查找點。我唯一的疑問就是這種方法適用於上面提到的P0和P1點嗎？

請評論並讓我知道您的想法。

Answer 1

兩個球體的方程可以寫成：

<X-P0,X-P0> - L0^2 = 0 (Eq0) 
<X-P1,X-P1> - L1^2 = 0 (Eq1) 

凡<U,V>表示點積。交點圓的中心（如果定義的話）是線P0，P1和由Eq0-Eq1（支持該圓）定義的平面之間的交點。這架飛機被稱爲兩架飛機的激進飛機。這個平面的方程爲（E）=（EQ0） - （EQ1）：

<P0,P0> - <P1,P1> + 2*<X,P1-P0> - L0^2 + L1^2 = 0 (E) 

表示由X線P0，P1的點（A）= A * P0 +（1-α）* P1，並注入（E）以獲得a中的線性方程。解是a0，圓的中心是C = X（a0）。請注意，C可能位於段P0，P1的外部（當一個球體的中心位於另一個球體的內部時）。我們得到：

2*a0 = 1 - (L0^2-L1^2)/dist(P0,P1)^2 

半徑圓的r然後得到解決：

dist(C,P0)^2+r^2=L0^2, or equivalently 
dist(C,P1)^2+r^2=L1^2 

它可能沒有解決辦法，如果球沒有交集。

Answer 2

我附上解決方案的代碼。 P0被認爲是肩膀點，P1被認爲是位於太空中的點（根據我的上臂和前臂的長度，我應該抓住P1）。請評論它，讓我知道你的想法。

#include <stdio.h> 
#include <math.h> 
#include <stdlib.h> 

struct point { 
    double x, y, z; 
}; 

/* 
* Used to represent a vector in 
* Xi + Yj + Zk format 
*/ 
struct vector { 
    double i, j, k; 
}; 

/* 
* Used to represent a plane in 
* Ax + By + Cz + D = 0 format 
*/ 
struct plane { 
    double A, B, C, D; 
}; 

/* 
* Represents the final assembly of the configuration. When two spheres 
* intersect they form a circle whose center is stored at "center" and radius is 
* stored at "radius". The circle also has a support which is defined by a plane 
* called as "radical plane" and its equation is stored at "p". 
*/ 
struct configuration { 
    struct point center; 
    double radius; 
    struct plane p; 
}; 

/* 
* Conversion functions between vector and point 
*/ 
struct vector get_vector_from_point(struct point p) { 
    static struct vector v; 
    v.i = p.x; 
    v.j = p.y; 
    v.k = p.z; 
    return v; 
} 

struct point get_point_from_vector(struct vector v) { 
    static struct point p; 
    p.x = v.i; 
    p.y = v.j; 
    p.z = v.k; 
    return p; 
} 

int check_if_same_points(struct point p1, struct point p2) { 
    return ((p1.x == p2.x) && (p1.y == p2.y) && (p1.z == p2.z)); 
} 
/* 
* Distance formula 
*/ 
double distance(struct point p0, struct point p1) { 
    return sqrt(pow((fabs(p1.z - p0.z)), 2) + pow((fabs(p1.y - p0.y)), 2) + pow((fabs(p1.x - p0.x)), 2)); 
} 

/* 
* Customized utility functions used for vector mathematics 
*/ 
double dot_product(struct vector p0, struct vector p1) { 
    return (p0.i * p1.i + p0.j * p1.j + p0.k * p1.k); 
} 
struct vector scale_vector(double scale, struct vector v) { 
    static struct vector scaled_vec; 
    scaled_vec.i = scale * v.i; 
    scaled_vec.j = scale * v.j; 
    scaled_vec.k = scale * v.k; 
    return scaled_vec; 
} 

struct vector add_vectors(struct vector v1, struct vector v2) { 
    static struct vector v; 
    v.i = v1.i + v2.i; 
    v.j = v1.j + v2.j; 
    v.k = v1.k + v2.k; 
    return v; 
} 

struct vector subtract_vectors(struct vector v1, struct vector v2) { 
    static struct vector v; 
    v.i = v1.i - v2.i; 
    v.j = v1.j - v2.j; 
    v.k = v1.k - v2.k; 
    return v; 
} 

/* 
* Takes the given assembly of points and links. Returns object of configuration 
* structure with necessary information. The center and radius from the returned 
* structure can be used find possible locations of elbow. 
* Client can use following parametric equation of circle in 3-D 
* X(t) = C + r (cos(t) * U + sin(t) * V) 
* where 0 <= t < 2 * pie, C is the center, r is the radius, U and V are unit 
* normals to the plane such that if N is a unit length plane normal then 
* {U,V,N} are mutually orthogonal. 
*/ 
struct configuration return_config(struct point p0, double l0, struct point p1, double l1) { 

    struct vector p0_v = get_vector_from_point(p0); 
    struct vector p1_v = get_vector_from_point(p1); 

    double dot_prd_p0 = dot_product(p0_v, p0_v); 
    double dot_prd_p1 = dot_product(p1_v, p1_v); 

    struct vector sub_vec = subtract_vectors(p1_v, p0_v); 
    double D = ((l0 * l0) - (l1 * l1) + dot_prd_p1 - dot_prd_p0)/2.0f; 

    static struct plane p; 
    p.A = sub_vec.i; p.B = sub_vec.j; p.C = sub_vec.k; p.D = D; 

    static struct configuration c; 

    /* 
     * Special case when object point and shoulder point are same. 
     */ 
    if(check_if_same_points(p0, p1)) { 
      printf("object and shoulder are at same location \n"); 
      c.center.x = p0.x; c.center.y = p0.y; c.center.z = p0.z; 
      c.radius = l0; 
      c.p.A = c.p.B = c.p.C = c.p.D = 0.0f; 
      return c; 
    } 

    double a0 = (1.0f - (((l0 * l0) - (l1 * l1))/(distance(p0, p1) * distance(p0, p1))))/2.0f; 


    struct vector lhs = scale_vector(a0,p0_v); 
    struct vector rhs = scale_vector(1.0f - a0, p1_v); 
    struct vector ans = add_vectors(lhs, rhs); 

    struct point center = get_point_from_vector(ans); 
    double radius = sqrt((l0 * l0) - (distance(center, p0) * distance(center, p0))); 

    c.center.x = center.x; c.center.y = center.y; c.center.z = center.z; 
    c.radius = radius; 
    c.p.A = p.A; c.p.B = p.B; c.p.C = p.C; c.p.D = D; 
    return c; 
} 

/* 
* The logic is as follows - Point P0 generates a sphere of radius L0 around it, 
* P1 generates another sphere of radius L1 around it. The intersection(if any) 
* will be a circle with a plane. If I can return the center and radius of that 
* circle and equation of the plane, then the client can find out any possible 
* location of the elbow by varying the value of theta in the parametric 
* equation of the circle. Thus if the spheres meet to generate a circle then 
* there will be infinite elbow positions. If it does not generate a circle and 
* meet "externally" then there will be only single elbow position. Otherwise 
* there will be no solutions at all. 
*/ 
int main() { 

    struct point p0, p1; 
    p0.x = 0, p0.y = 0, p0.z = 0; 
    p1.x = 50, p1.y = 50, p1.z = 0; 
    double l0 = 50, l1 = 50; 

    printf("Shoulder coordinates : (%lf, %lf, %lf) \n", p0.x, p0.y, p0.z); 
    printf("Object coordinates: (%lf, %lf, %lf) \n", p1.x, p1.y, p1.z); 
    printf("link0 = %lf, link1 = %lf \n", l0, l1); 

    if(distance(p0, p1) > (l0 + l1)) { 
      printf("The given combination of the points and links cannot make a valid configuration"); 
      return -1; 
    } 
    struct configuration c = return_config(p0, l0, p1, l1); 
    printf("Center = (%lf, %lf, %lf), radius = %lf \n", c.center.x, c.center.y, c.center.z, c.radius); 
    printf("Equation of the radical plane = %lfA + (%lf)B + (%lf)C + %lf = 0 \n", c.p.A, c.p.B, c.p.C, c.p.D); 
    return 0; 
}