views:

513

answers:

2

Hi,

Consider the following problem - I am given 2 links of length L0 and L1. P0 is the point that the first link starts at and P1 is the point that I want the end of second link to be at in 3-D space. I am supposed to write a function that should take in these 3-D points (P0 and P1) as inputs and should find all configurations of the links that put the second link's end point at P1.

My understanding of how to go about it is - Each link L0 and L1 will create a sphere S0 and S1 around itself. I should find out the intersection of those two spheres (which will be a circle) and print all points that are on the circumference of that circle.

I saw gmatt's first reply on the http://stackoverflow.com/questions/1406375/finding-intersection-points-between-3-spheres but could not understand it properly since the images did not show up. I also saw a formula for finding out the intersection at http://mathworld.wolfram.com/Sphere-SphereIntersection.html

I could find the radius of intersection by the method given on mathworld. Also I can find the center of that circle and then use the parametric equation of circle to find the points. The only doubt that I have is will this method work for the points P0 and P1 mentioned above ?

Please comment and let me know your thoughts.

+1  A: 

The equations of the two spheres may be written:

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

Where <U,V> denotes the dot product. The center of the intersection circle, if defined, is the intersection between line P0,P1 and the plane defined by Eq0-Eq1 (support of the circle). This plane is known as the radical plane of the two spheres. The equation of this plane is (E)=(Eq0)-(Eq1):

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

Represent a point on line P0,P1 by X(a)=a*P0+(1-a)*P1, and inject in (E) to obtain a linear equation in a. Solution is a0, and the center of the circle is C=X(a0). Note that C may be outside segment P0,P1 (when the center of one sphere is inside the other). We get:

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

The radius r of the circle is then obtained solving:

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

It may not have solutions if the spheres have no intersections.

Eric Bainville
Thanks a lot Eric for your reply. I think I got the solution and I am posting it as an answer.
Onkar Deshpande
A: 

I am attaching the code for solution. P0 is considered as a shoulder point and P1 is considered as the point that lies in space (which I am supposed to grab given the length of my upper arm and forearm). Please comment on it and let me know what you think.

#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;
}
Onkar Deshpande