```Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Triangulation using sphere intersects
Date: Fri, 21 Nov 2008 23:04:02 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 39
Message-ID: <gg7ep2\$dg\$1@fred.mathworks.com>
References: <gg5tb3\$ms3\$1@fred.mathworks.com>
NNTP-Posting-Host: webapp-03-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1227308642 432 172.30.248.38 (21 Nov 2008 23:04:02 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Fri, 21 Nov 2008 23:04:02 +0000 (UTC)
Xref: news.mathworks.com comp.soft-sys.matlab:502434

"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

```