Thread Subject:

Triangulation using sphere intersects

Im trying to triangulate locations in xyz space.

I have locations for the centre of spheres and what can be thought of as radii, so what is the best way with Matlab to compute intersects to
1) get and plot the circle where 2 spheres intersect or
2) get and plot the 2 points where 3 spheres may intersect.

Many Thanks
Greg

"gregoire " <s0800314@nus.edu.sg> wrote in message <gg5tb3$ms3$1@fred.mathworks.com>...
> Im trying to triangulate locations in xyz space.
>
> I have locations for the centre of spheres and what can be thought of as radii, so what is the best way with Matlab to compute intersects to
> 1) get and plot the circle where 2 spheres intersect or
> 2) get and plot the 2 points where 3 spheres may intersect.
>
> Many Thanks
> Greg

  I'll just outline how you could proceed in your problem, Greg. For problem 1) requiring plotting the circle of intersection of two spheres, the two equations of the spheres can be expressed:

 (x-x1)^2+(y-y1)^2+(z-z1)^2 = r1^2
 (x-x2)^2+(y-y2)^2+(z-z2)^2 = r2^2

where (x1,y1,z1) and (x2,y2,z2) are the two centers and r1 and r2 the two radii. If you subtract one equation from the other you will have the linear equation satisfied by the plane which contains the desired circle.

  Then the line drawn through the two centers can be expressed parametrically by

 x = t*x1 + (1-t)*x2
 y = t*y1 + (1-t)*y2
 z = t*z1 + (1-t)*z2

These three equations together with the above plane equation are all linear with four unknowns and four equations and can be solved for the point P = (x,y,z) on the connecting line at the intersection with the plane. This is the center of your desired circle.

 Since the vector [x1-x2;y1-y2;z1-z2] is parallel to the above line, you can use the Matlab 'null' function to find two mutually orthogonal unit vectors, u and v, which are both orthogonal to it. By the Pythagoras theorem the square of the circle's radius, R, is the square of one of the two spheres' radii minus the square of the distance from its center to the point P. Hence you can express the circle parametrically as

 Q = P + R*u*cos(s) + R*v*sin(s)

where s is an (angular) parameter and Q an arbitrary point on the circle. You could use 'plot3' to plot it, letting s range from 0 to 2*pi.

  For problem 2) subtract one of the spheres' equations from the other two obtaining two linear equations each of which represents a plane. Then solve for the plane which contains all three centers. The intersection of all these three planes gives a unique point, L. The desired points of intersection of the three spheres must lie along a line through L orthogonal to the plane containing the three centers. Its distance along this line in either direction can again be determined by Pythagoras' theorem using one of the spheres' radii and the distance from L to that sphere's center.

  Of course in both cases if there is no intersection you will be able to tell by the failure of the Pythagoras calculation.

  The solution to 1) and 2) could be reduced to closed formulas in terms of the spheres' centers and radii but I don't have time to work that out for you, so I'll leave that to you.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gg7ep2$dg$1@fred.mathworks.com>...
> .......
> Then solve for the plane which contains all three centers.
> .......

  It occurs to me that I might have taken too large a step when I said to find the equation of the plane through the three sphere centers. To give that step in greater detail, let C1 = [x1;y1;z1], C2 = [x2;y2;z2], and C3 = [x3;y3;z3] be the three centers. Then the cross product, cross(C2-C1,C3-C1), must be orthogonal to the desired plane. Hence any point on the plane P = (x,y,z) must satisfy

 dot(P-C1,cross(C2-C1,C3-C1)) =
 (x-x1)*((y2-y1)*(z3-z1)-(y3-y1)*(z2-z1)) +
 (y-y1)*((z2-z1)*(x3-x1)-(z3-z1)*(x2-x1)) +
 (z-z1)*((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)) = 0

which gives the equation of that plane.

Roger Stafford

Im back on this personal project and let me first say Thank you very much for your advice Roger, it has been very helpful indeed.

It has occured to me that if I had 2 spheres which intersected I'd get a circle, and 3 spheres which intersected I could resolve this down to 2 points, however in coding this problem, I could just find all circle intersects between all pair of intersecting spheres, regardless if there are 2 or 3 intersecting spheres, and then at a later time around, find intersects between these circles (where 3 spheres had intersected) in order to resolve this down to the 2 points? As opposed to treating the case where 2 spheres or 3 sheres intersect, separately. This is correct right?

Working on the above assumption, I have pretty much managed to calculate the centre, and the 2 orthogonal vectors, ie all the information I need to depict for all sphere pair-intersecting circles, in separate rows. So I now have thousands of these circles with the required information in separate rows, and am now having difficulty with a) using plot3 to plot these many many circles together in 3d space and b) attempting to find potential intersections betweens these circles so as to resolve further, some of these circles into just 2 points for increased accuracy, which arise when 3 or more spheres intersect.

Thank you

"gregoire " <dfd@invalid.edu.sg> wrote in message <glmpbp$cte$1@fred.mathworks.com>...
> Im back on this personal project and let me first say Thank you very much for your advice Roger, it has been very helpful indeed.
>
> It has occured to me that if I had 2 spheres which intersected I'd get a circle, and 3 spheres which intersected I could resolve this down to 2 points, however in coding this problem, I could just find all circle intersects between all pair of intersecting spheres, regardless if there are 2 or 3 intersecting spheres, and then at a later time around, find intersects between these circles (where 3 spheres had intersected) in order to resolve this down to the 2 points? As opposed to treating the case where 2 spheres or 3 sheres intersect, separately. This is correct right?
>
> Working on the above assumption, I have pretty much managed to calculate the centre, and the 2 orthogonal vectors, ie all the information I need to depict for all sphere pair-intersecting circles, in separate rows. So I now have thousands of these circles with the required information in separate rows, and am now having difficulty with a) using plot3 to plot these many many circles together in 3d space and b) attempting to find potential intersections betweens these circles so as to resolve further, some of these circles into just 2 points for increased accuracy, which arise when 3 or more spheres intersect.
>
> Thank you

is this a good derivation? seems simple enough...
http://en.wikipedia.org/wiki/Trilateration

google provides many more methods, like:
http://mathforum.org/library/drmath/view/63138.html
http://airglow.csl.uiuc.edu/Teaching/ECE498AL/ECE498AL_HW02_SOLN.pdf
and a bunch more, including some interesting discussions of errors in gps and related locating systems.

Thanks David for your response, some interesting reads there, especially as I may have to look into errors and noise variation in my data at a later stage.

I dont wish to find a new method (though the Newton-Rapshon one seems interesting) for now Im just trying to work out how to

a) use plot3 to plot many circles together in 3d space, which are of the form:-

 Q = P + R*u*cos(s) + R*v*sin(s)
where s is an (angular) parameter and Q an arbitrary point on the circle. u,v are necessary orthogonal vectors,

as that is the way I have already coded and successfully calculated the necessary parts, from my data

and b) attempting to find potential intersections betweens these circles

Thanks again

"gregoire " <dfd@invalid.edu.sg> wrote in message <glmsf6$7b6$1@fred.mathworks.com>...
> Thanks David for your response, some interesting reads there, especially as I may have to look into errors and noise variation in my data at a later stage.
>
> I dont wish to find a new method (though the Newton-Rapshon one seems interesting) for now Im just trying to work out how to
>
> a) use plot3 to plot many circles together in 3d space, which are of the form:-
>
> Q = P + R*u*cos(s) + R*v*sin(s)
> where s is an (angular) parameter and Q an arbitrary point on the circle. u,v are necessary orthogonal vectors,
>
> as that is the way I have already coded and successfully calculated the necessary parts, from my data
>
> and b) attempting to find potential intersections betweens these circles
>
> Thanks again

  I don't think I agree with you, Gregoire, that finding the intersections of the parametric circles is the best way to find the intersection points of triplets of spheres. The way I outlined to you involving three intersecting planes and then a certain distance in either direction along a line orthogonal to the plane of the three centers is a far better method which does not involve iterative methods like Newton-Raphson. You get the two points all in one clean shot without convergence worries or iterations. The fact that the circles will run through these points does not help you in my opinion. I think you will find that things are easier that way.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <glohe8$ai1$1@fred.mathworks.com>...
> "gregoire " <dfd@invalid.edu.sg> wrote in message <glmsf6$7b6$1@fred.mathworks.com>...
> > Thanks David for your response, some interesting reads there, especially as I may have to look into errors and noise variation in my data at a later stage.
> >
> > I dont wish to find a new method (though the Newton-Rapshon one seems interesting) for now Im just trying to work out how to
> >
> > a) use plot3 to plot many circles together in 3d space, which are of the form:-
> >
> > Q = P + R*u*cos(s) + R*v*sin(s)
> > where s is an (angular) parameter and Q an arbitrary point on the circle. u,v are necessary orthogonal vectors,
> >
> > as that is the way I have already coded and successfully calculated the necessary parts, from my data
> >
> > and b) attempting to find potential intersections betweens these circles
> >
> > Thanks again
>
> I don't think I agree with you, Gregoire, that finding the intersections of the parametric circles is the best way to find the intersection points of triplets of spheres. The way I outlined to you involving three intersecting planes and then a certain distance in either direction along a line orthogonal to the plane of the three centers is a far better method which does not involve iterative methods like Newton-Raphson. You get the two points all in one clean shot without convergence worries or iterations. The fact that the circles will run through these points does not help you in my opinion. I think you will find that things are easier that way.
>
> Roger Stafford

i think a more general solution is the newton-raphson approach, especially if he expects to be considering uncertainty in the future. if the spheres are not known exactly, or there are introduced errors where maybe two of the spheres only get very close but don't actually intersect, then an iterative approach may be the only option to find the best possible answer.

"gregoire " <dfd@invalid.edu.sg> wrote in message <glmpbp$cte$1@fred.mathworks.com>...
> Im back on this personal project and let me first say Thank you very much for your advice Roger, it has been very helpful indeed.
>
> It has occured to me that if I had 2 spheres which intersected I'd get a circle, and 3 spheres which intersected I could resolve this down to 2 points, however in coding this problem, I could just find all circle intersects between all pair of intersecting spheres, regardless if there are 2 or 3 intersecting spheres, and then at a later time around, find intersects between these circles (where 3 spheres had intersected) in order to resolve this down to the 2 points? As opposed to treating the case where 2 spheres or 3 sheres intersect, separately. This is correct right?
>
> Working on the above assumption, I have pretty much managed to calculate the centre, and the 2 orthogonal vectors, ie all the information I need to depict for all sphere pair-intersecting circles, in separate rows. So I now have thousands of these circles with the required information in separate rows, and am now having difficulty with a) using plot3 to plot these many many circles together in 3d space and b) attempting to find potential intersections betweens these circles so as to resolve further, some of these circles into just 2 points for increased accuracy, which arise when 3 or more spheres intersect.
>
> Thank you

  Gregoire, I finally decided to stop being lazy and have worked out the matlab code for your problem 2. Given three spheres, you want to find their two points of intersection (if they actually intersect.) As you can see, it gives an answer without any iteration. My claim is that the only time these calculations become subject to inaccuracy is when the problem is inherently ill-conditioned, as for example when the centers are nearly colinear, or when the two points of intersection are very close together.

  Let p1, p2, and p3 be three-element vectors giving the x,y,z coordinates of the spheres' three centers, and r1, r2, and r3 their respective radii. Then do this:

 p21 = p2-p1;
 p31 = p3-p1;
 c = cross(p21,p31);
 c2 = sum(c.^2);
 u1 = cross(((sum(p21.^2)+r1^2-r2^2)*p31 - ...
             (sum(p31.^2)+r1^2-r3^2)*p21)/2,c)/c2;
 v = sqrt(r1^2-sum(u1.^2))*c/sqrt(c2);
 i1 = p1+u1+v;
 i2 = p1+u1-v;

The i1 and i2 will be vectors for the two points of intersection. If there is any possibility of having radii such that no solution is possible (which can happen very easily) you should test that the first square root used in the calculation of v is real. An impossible problem will make it imaginary.

  By way of explanation, the vector c is orthogonal to the plane of the three centers. The vector u1 is the location of the three-plane intersection I mentioned earlier relative to the point p1, and can be obtained by solving the three linear equations in three unknowns I mentioned previously. Vector v stretches from this plane-intersection point along the c direction to one of the three-sphere intersection points, while -v stretches to the other one.

  If the centers are nearly colinear, the cross product c is nearly zero and division by it becomes increasingly inaccurate as colinearity is approached. That is inherent in the situation and no fault of the method. When r1^2 and sum(u1.^2) are nearly equal the square root process becomes inaccurate, and again that is inherent in the geometry of the problem as the two intersection points approach each other.

  David, I disagree with your statement that there is a gain in accuracy from an iterative approach. The very factors that lead to inaccuracy of the above code are the same ones that will cause inaccuracy in any iterative method. Small errors in center locations or radii will lead to disproportionately larger errors in intersection positions in these cases for either type of method.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gg7ibf$iq3$1@fred.mathworks.com>...
> "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gg7ep2$dg$1@fred.mathworks.com>...
> > .......
> > Then solve for the plane which contains all three centers.
> > .......
>
> It occurs to me that I might have taken too large a step when I said to find the equation of the plane through the three sphere centers. To give that step in greater detail, let C1 = [x1;y1;z1], C2 = [x2;y2;z2], and C3 = [x3;y3;z3] be the three centers. Then the cross product, cross(C2-C1,C3-C1), must be orthogonal to the desired plane. Hence any point on the plane P = (x,y,z) must satisfy
>
> dot(P-C1,cross(C2-C1,C3-C1)) =
> (x-x1)*((y2-y1)*(z3-z1)-(y3-y1)*(z2-z1)) +
> (y-y1)*((z2-z1)*(x3-x1)-(z3-z1)*(x2-x1)) +
> (z-z1)*((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)) = 0
>
> which gives the equation of that plane.
>
> Roger Stafford
Greetings Roger
I have followed the steps you specified in this forum to compute the intersection of three spheres. However, I got stuck with some instructions.
In order to finally produce the intercept coordinate points, a line equation must be implemented. This line must contain L (point of interception of the three specified planes) and should be orthogonal to the plane containing the 3 sphere centers. The definition of a line equation is:
<x,y,z> = <x0, y0,z0> + t*<Vx, Vy, Vz>
where <x0, y0, z0> refers to a point in the line, and <Vx, Vy,Vz> refers to a vector parallel to it. The way I figure, if I have this equation defined, I could use it along with L, the computed distance (d) from L to the point of sphere intersection (obtained using Pythagoras theorem), and the distance formula to determine the point of interest:
d^2 = [(x-Lx)^2 + (y-Ly)^2 + (z-Lz)^2].
if (x,y,z) lies along the same line as L, then:
d^2 = [(Lx + tVx-Lx)^2 + (Ly + tVy -Ly)^2 + (Lz - tVz-Lz)^2]
All I would have to do then is to solve for t. Then, the point of interception would be obtained using the Line equation.
It seems dandy and all, but I have a little problem with this. I just do not know how I could obtain the vector needed to define the Line Equation. So, this leads me to the following questions:
Is the approach I have mentioned over here appropriate? If so, what am I doing wrong? What am I missing? Is there a way I could obtain such a vector?
Is there something I could do instead of using the previously mentioned?
I would appreciate if you could help me out with this.

regards, Ricardo
Thanks
<x0, y0, z0> would be L. But I am

"Roger Stafford" wrote in message <gg7ep2$dg$1@fred.mathworks.com>...
> "gregoire " <s0800314@nus.edu.sg> wrote in message <gg5tb3$ms3$1@fred.mathworks.com>...
> > Im trying to triangulate locations in xyz space.
> >
> > I have locations for the centre of spheres and what can be thought of as radii, so what is the best way with Matlab to compute intersects to
> > 1) get and plot the circle where 2 spheres intersect or
> > 2) get and plot the 2 points where 3 spheres may intersect.
> >
> > Many Thanks
> > Greg
>
> I'll just outline how you could proceed in your problem, Greg. For problem 1) requiring plotting the circle of intersection of two spheres, the two equations of the spheres can be expressed:
>
> (x-x1)^2+(y-y1)^2+(z-z1)^2 = r1^2
> (x-x2)^2+(y-y2)^2+(z-z2)^2 = r2^2
>
> where (x1,y1,z1) and (x2,y2,z2) are the two centers and r1 and r2 the two radii. If you subtract one equation from the other you will have the linear equation satisfied by the plane which contains the desired circle.
>
> Then the line drawn through the two centers can be expressed parametrically by
>
> x = t*x1 + (1-t)*x2
> y = t*y1 + (1-t)*y2
> z = t*z1 + (1-t)*z2
>
> These three equations together with the above plane equation are all linear with four unknowns and four equations and can be solved for the point P = (x,y,z) on the connecting line at the intersection with the plane. This is the center of your desired circle.
>
> Since the vector [x1-x2;y1-y2;z1-z2] is parallel to the above line, you can use the Matlab 'null' function to find two mutually orthogonal unit vectors, u and v, which are both orthogonal to it. By the Pythagoras theorem the square of the circle's radius, R, is the square of one of the two spheres' radii minus the square of the distance from its center to the point P. Hence you can express the circle parametrically as
>
> Q = P + R*u*cos(s) + R*v*sin(s)
>
> where s is an (angular) parameter and Q an arbitrary point on the circle. You could use 'plot3' to plot it, letting s range from 0 to 2*pi.
>
> For problem 2) subtract one of the spheres' equations from the other two obtaining two linear equations each of which represents a plane. Then solve for the plane which contains all three centers. The intersection of all these three planes gives a unique point, L. The desired points of intersection of the three spheres must lie along a line through L orthogonal to the plane containing the three centers. Its distance along this line in either direction can again be determined by Pythagoras' theorem using one of the spheres' radii and the distance from L to that sphere's center.
>
> Of course in both cases if there is no intersection you will be able to tell by the failure of the Pythagoras calculation.
>
> The solution to 1) and 2) could be reduced to closed formulas in terms of the spheres' centers and radii but I don't have time to work that out for you, so I'll leave that to you.
>
> Roger Stafford

how will u find the null function of orthogonal vector u and v...
please help me
regards satish

"gregoire" wrote in message <gg5tb3$ms3$1@fred.mathworks.com>...
> Im trying to triangulate locations in xyz space.
>
> I have locations for the centre of spheres and what can be thought of as radii, so what is the best way with Matlab to compute intersects to
> 1) get and plot the circle where 2 spheres intersect or
> 2) get and plot the 2 points where 3 spheres may intersect.
>
> Many Thanks
> Greg

dear gregoire

 I am satish sahu, a student in my thesis i have certain problem in which i have to find and get the plot of circle where two sphere are intersecting.
have you developed any MATLAB code regarding this matter or could you tell me where can i find the code or u have any algorithm.

Please reply soon it would be great help.

Thanks.

"satish sahu" <satishs@iitk.ac.in> wrote in message <jfu4k6$lav$1@newscl01ah.mathworks.com>...
> I am satish sahu, a student in my thesis i have certain problem in which i have to find and get the plot of circle where two sphere are intersecting.
> have you developed any MATLAB code regarding this matter or could you tell me where can i find the code or u have any algorithm.
- - - - - - - - -
  Hello Satish Sahu. I last communicated on this thread three years ago, so it has taken me a while to refresh my memory on this subject. You wish to plot the circle of intersection between two given spheres. The following will accomplish that in matlab code, and derives roughly from the "problem 1)" part of the outline I described in this thread 21 Nov. 2008.

  Let C1 = [x1,y1,z1] and C2 = [x2,y2,z2] be two row vectors defining the spheres' centers and r1 and r2 their respective radii. Then do this:

 C21 = C2-C1;
 d2 = dot(C21,C21);
 C0 = (C1+C2)/2+(r1^2-r2^2)/(2*d2)*C21;
 R = sqrt(((r1+r2)^2-d2)*(d2-(r2-r1)^2)/(4*d2));
 N = null(C21).';
 T = linspace(0,2*pi).';
 V = bsxfun(@plus,C0,R*[cos(T),sin(T)]*N);

  The point C0 is the center of the intersection circle and R is its radius. C0 lies along the line connecting the centers, C1 and C2. The two rows of N, produced with matlab's 'null' function, are mutually orthogonal unit vectors which are also orthogonal to the vector C21 = C2-C1. The column vectors V(:,1), V(:,2), and V(:,3) are the X, Y, and Z vectors you need to do a plot of the circle using 'plot3'. They trace out the path of the circle around a full angle of 2*pi about C0.

  You can verify the accuracy of this code by checking that the distance from points on the circle are all a distance r1 from C1 and r2 from C2. You can also check that vectors pointing from C0 to points on the circle are orthogonal to C2-C1 and are all of length R, which shows that C0 is indeed at the circle's center.

  In case your spheres don't intersect, that will result in taking the square root of a negative number in computing R, so you might want to put a test for that in the above code.

  Note that on plots of this circle with 'plot3' the circle may not appear orthogonal to the line between C1 and C2 unless you have set the plot scaling appropriately.

Roger Stafford

HI Roger Stafford
many many thanks....

"Roger Stafford" wrote in message <glqslm$fjp$1@fred.mathworks.com>...
>
> Let p1, p2, and p3 be three-element vectors giving the x,y,z coordinates of the spheres' three centers, and r1, r2, and r3 their respective radii. Then do this:
>
> p21 = p2-p1;
> p31 = p3-p1;
> c = cross(p21,p31);
> c2 = sum(c.^2);
> u1 = cross(((sum(p21.^2)+r1^2-r2^2)*p31 - ...
> (sum(p31.^2)+r1^2-r3^2)*p21)/2,c)/c2;
> v = sqrt(r1^2-sum(u1.^2))*c/sqrt(c2);
> i1 = p1+u1+v;
> i2 = p1+u1-v;
>

First, thank you Roger!
Your codes work fine. Yet, I'm not quite understand how you manage the vector u1. It seems to contain and simplify some formulas or something? Would you please explain in more details about..
> u1 = cross(((sum(p21.^2)+r1^2-r2^2)*p31 - ...
> (sum(p31.^2)+r1^2-r3^2)*p21)/2,c)/c2;

Many Thanks

"Kus" wrote in message <k0fjsh$l9p$1@newscl01ah.mathworks.com>...
> "Roger Stafford" wrote in message <glqslm$fjp$1@fred.mathworks.com>...
> >
> > Let p1, p2, and p3 be three-element vectors giving the x,y,z coordinates of the spheres' three centers, and r1, r2, and r3 their respective radii. Then do this:
> >
> > p21 = p2-p1;
> > p31 = p3-p1;
> > c = cross(p21,p31);
> > c2 = sum(c.^2);
> > u1 = cross(((sum(p21.^2)+r1^2-r2^2)*p31 - ...
> > (sum(p31.^2)+r1^2-r3^2)*p21)/2,c)/c2;
> > v = sqrt(r1^2-sum(u1.^2))*c/sqrt(c2);
> > i1 = p1+u1+v;
> > i2 = p1+u1-v;
>
> First, thank you Roger!
> Your codes work fine. Yet, I'm not quite understand how you manage the vector u1. It seems to contain and simplify some formulas or something? Would you please explain in more details about..
> > u1 = cross(((sum(p21.^2)+r1^2-r2^2)*p31 - ...
> > (sum(p31.^2)+r1^2-r3^2)*p21)/2,c)/c2;
>
> Many Thanks
- - - - - - - -
  Hello Kus. As you may recall from earlier in this thread, on Jan. 29, 2009 I stated that "The vector u1 is the location of the three-plane intersection I mentioned earlier relative to the point p1, and can be obtained by solving the three linear equations in three unknowns I mentioned previously." Earlier on Nov. 21, 2008 I stated, "subtract one of the spheres' equations from the other two obtaining two linear equations each of which represents a plane." The subtractions p21 = p2-p1 and p31 = p3-p1 represent the centers p2 and p3 in terms of p21 and p31 taking p1 as translated origin.

  From all this you should be able to deduce these three linear equations which represent the intersection of three planes:

 dot(u1,cross(p21,p31)) = 0
 dot(u1,p21) = (sum(p21.^2)+r1^2-r2^2)/2
 dot(u1,p31) = (sum(p31.^2)+r1^2-r3^2)/2

where the unknown components of u1 are cartesian coordinates that have been translated to center p1 as origin.

  It can be shown that the expression which I gave

 u1 = cross(((sum(p21.^2)+r1^2-r2^2)*p31 - ...
             (sum(p31.^2)+r1^2-r3^2)*p21)/2,c)/c2;

gives the unique solution to these equations. For the purpose of demonstrating this, it is useful to recall the three-dimensional vector analysis identities:

 cross(A,cross(B,C)) = dot(A,C)*B - dot(A,B)*C

and

 dot(A,cross(B,C)) = dot(C,cross(A,B)) = dot(B,cross(C,A))

which can be used to show that the above u1 is identically equal to

 ( (a*dot(p31,p31)-b*dot(p21,p31))*p21 + ...
   (b*dot(p21,p21)-a*dot(p21,p31))*p31 ) / ...
 (dot(p21,p21)*dot(p31,p31) - dot(p21,p31)^2)

where

 a = (sum(p21.^2)+r1^2-r2^2)/2
 b = (sum(p31.^2)+r1^2-r3^2)/2

From this second form of u1 it can be readily shown that the second two equations hold. It is obvious from this second form of u1 that it is a linear combination of p21 and p31 and therefore must satisfy the first equation.

  I hope this explanation has helped. Admittedly one could use matlab's matrix division operator to solve the three equations, but I chose to use this explicit form as a solution to Gregoire's original second problem.

Roger Stafford

Roger,
Your explanation has helped me a lot. Thank you so much!

"Kus" wrote in message <k0s7dj$6b0$1@newscl01ah.mathworks.com>...
> Roger,
> Your explanation has helped me a lot. Thank you so much!
- - - - - - -
  I'm afraid I went about demonstrating those three equations the hard way, Kus. It's not necessary to convert u1 to that second form.

  The three equations in question are:

 dot(u1,cross(p21,p31)) = 0
 dot(u1,p21) = a
 dot(u1,p31) = b

where

 a = (sum(p21.^2)+r1^2-r2^2)/2 and
 b = (sum(p31.^2)+r1^2-r3^2)/2

and where u1 can then be expressed in terms of a and b as

 u1 = cross(a*p31-b*p21,cross(p21,p31)) / ...
      dot(cross(p21,p31),cross(p21,p31))

  Applying the "permutation" identity of vector analysis for scalar triple products,

 dot(cross(A,B),C) = dot(cross(B,C),A) = dot(cross(C,A),B) ,

to the numerator of the first equation, the simpler procedure would be

 dot(cross(a*p31-b*p21,cross(p21,p31)),cross(p21,p31)) =
 dot(cross(cross(p21,p31),cross(p21,p31)),a*p31-b*p21) =
 dot(0,a*p31-b*p21) = 0

which implies the first equation.

  The numerator of dot(u1,p21) in the second equation can be expressed

 dot(cross(a*p31-b*p21,cross(p21,p31)),p21) =
 dot(cross(p21,a*p31-b*p21),cross(p21,p31)) =
 dot(a*cross(p21,p31)-b*cross(p21,p21),cross(p21,p31) ) =
 dot(a*cross(p21,p31)-0,cross(p21,p31)) =
 a*dot(cross(p21,p31),cross(p21,p31))

When this is divided by the denominator, the two dot products cancel leaving just 'a' and we get the second equation

 dot(u1,p21) = a.

  A similar reasoning will demonstrate the third equation, dot(u1,p31) = b.

Roger Stafford