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