From: rusin@washington.math.niu.edu (Dave Rusin)
Newsgroups: comp.soft-sys.math.mathematica
Subject: Re: Triangulation Problem
Date: 20 Jan 1996 08:28:12 GMT
In article <4diafm$84a@dragonfly.wri.com>, david hoare wrote:
>
>I am looking for the trigonometric / algebreic formulae one would use
>to triangulate an unknown position, knowing 3 fixed points and the
>respective distances to the unknown point. (make sence?)
>Oh - in 3 dimensional space.
So you want to know what point is of distance D1 from point P1, and likewise
of distances D2 and D3 from P2 and P3 respectively.
Think geometrically!
The set of all points meeting the first condition is a sphere of radius D1
centered at P1. The set of points meeting the first two conditions is
then the intersection of two such spheres, which is clearly a
circle, perpendicular to the line segment joining P1 and P2, and of some
radius r, where I can compute r knowing the distance D12 between
P1 and P2. (You'll have to draw your own picture since I failed my
kidergarten ASCII-art class): if x is the distance from P1 to the
center of this circle, and thus D12-x the distance from there to P2,
we need r^2+x^2=D1^2 and r^2+(D12-x)^2=D2^2, forcing
x = [D12^2+D1^2-D2^2]/(2 D12)
and thus r = sqrt(D1^2-x^2) is known.
At this point I think it's easiest to parameterize the circle: Let
v0 = (P2 - P1)/ D12 be the unit vector from P1 to P2. If v is any
vector at all not parallel to v0, then the cross-products
v1 = v0 x v and v2 = v1 x v0 are perpendicular to each other and
to v0. Scale v1 and v2 to make them also of unit length.
Then the circle described in the previous paragraph is the set of points
P = r cos(t) v1 + r sin(t) v2 + (P1 + x v0)
Now apply your third condition: you want dist(P, P3) = D3. If the vector
from P1 to P3 is decomposed in the current coordinate system as
a v1+ b v2 + c v0, then that means you need
D3^2 = (r cos(t)-a)^2 + (r sin(t)-b)^2 + (x-c)^2.
This is a relatively simple trigonometric equation to solve; once the
smoke clears it looks like
A cos(t) + B sin(t) = C.
I usually solve such equations by dividing by sqrt(A^2+B^2). Choosing an
angle u so that sin(u) = A/sqrt(A^2+B^2) and cos(u)= B/sqrt(A^2+B^2),
we have just the equation
sin(t+u) = C/sqrt(A^2+B^2).
Take arcsines and subtract u to discover t; plug back into the
parameterization of the circle to find your point.
Alternatively you can replace the trig of the previous paragraph with
some algebra by parameterizing the unit circle as the set of points
of the form ( (1-t^2)/(1+t^2), (2t)/(1+t^2) ), instead of the usual
(cos(t), sin(t)). I don't think this is particularly easier, unless you
have an aversion to trig functions and their inverses.
>And actually, the reverse operations would be handy too...ie: given a
>known central point, AND 3 other fixed points, how to find the lengths
>of the lines required to meet at that first point, in 3D.
The lengths of these lines are simply given by the distance formula
(Pythagorean theorem).
dave
==============================================================================
From: Preston Nichols
Newsgroups: comp.soft-sys.math.mathematica
Subject: Re: [mg2968] Triangulation Problem
Date: 30 Jan 1996 06:37:45 GMT
David Hoare wrote:
<<<<
I am looking for the trigonometric / algebreic formulae one would
use to triangulate an unknown position, knowing 3 fixed points and
the respective distances to the unknown point. (make sense?) Oh - in
3 dimensional space.
[....]
It's for a program I'm writing, and it's been a while since I've done
this sort of thing - - my math has left me!
Any Help would be GREATLY appreciated!!
>>>>
I have seen essentially two responses to this, and I have a few
comments on each.
1. Anthony D. Gleckler's approach [mg3007], slightly paraphrased
(and the approach I used):
If a and b are points (given as lists of coordinates), then the
distance between them is, in Mathematica,
Sqrt[(a-b).(a-b)]
Let the three given points be
a = {a1,a2,a3};
b = {b1,b2,b3};
c = {c1,c2,c3};
let the associated distances be A,B,C, and let the unknown point be
X = {x,y,z};
Then the desired point is the simultaneous solution of these three
equations:
spherea := ( (X-a).(X-a) == A^2 );
sphereb := ( (X-b).(X-b) == B^2 );
spherec := ( (X-c).(X-c) == C^2 );
each of which just states the desired equality of distances, but
with both sides squared. This is three quadratic equations in three
unknowns;
sols = Solve[{eqa,eqb,eqc}, X]
will return the solution as replacement rules, and
soughtpoints = X/.sols
will give the unknown point as a list of coordinates -- but not
really! First of all, the Mathematica commands given above work
fine for numerical values for a, b, c, A, B, and C. If you try to
get a general symbolic solution, however, the expressions get
unmanageably big.
But there's another issue where writing the necessary code is not
hard, but deciding just what you want the code to do requires some
reflection. Unless the measurements of the given data (a, b, c, A,
B, and C) are absolutely perfect, the system of equations will have
TWO DISTINCT SOLUTIONS, possibly complex (having equal and opposite
imaginary parts), and the soughtpoints as computed above will
actually be a list of TWO DISTINCT POINTS, possibly with complex
coordinates. Each of the equations shperea, sphereb, and spherec
define a sphere; two of the spheres generically intersect in a
circle, which will generically intersect the third sphere in two
points. These two points (if real) determine a line perpendicular
to the plane containing the given points a,b,c, and the two points
are equidistant from that plane.
If the measurements of a, b, c, A, B, and C are only slightly off,
these two points will be very close together, and simply averaging
them will give a reasonable, usable result:
thepoint = Apply[Plus, soughtpoints]/2
(Conveniently, this will give a real result even if complex
solutions are involved!)
You might, however, wish to use the distance between the two
solutions as an indicator of the quality of your data (i.e. the
given points and distances). At the simplest level, a relatively
"large" distance between the two solution points indicates that your
values for A, B, and C *can't* be accepted as distances of a, b,
and c from a common fourth point; you could use that "largeness"
simply as a tipoff that something's wrong with the data. With some
more effort, you could also convert the size of that distance into a
numerical measurement of how well all your pieces are fitting
together.
---------------------------------------------------
2. Dave Rusin's approach [mg3015]
I really enjoyed Dave's method, because it so closely follows what
we might do if we could directly manipulate abstract geometrical
objects in 3-space with our hands. He suggests locating and
explicitly parametrizing the circle where the first two spheres
intersect. (I also think parametrizing curves is underappreciated
in general.) This can be done in Mathematica by using some tricks
with NullSpace and LinearAlgebra`Orthogonalization`GramSchmidt.
Dave's method will lead to the same two solutions for less than
perfect data, but this time they appear in the (ususally) multiple
solutions of the trigonometric equation
> D3^2 = (r cos(t)-a)^2 + (r sin(t)-b)^2 + (x-c)^2
the solving of which involves taking multi-valued arcsines. (In 1,
we might have been alerted to multiple solutions by the fact that
we were solving quadratic equations.)
---------------------------------------------------
Preston Nichols
Visiting Assistant Professor
Department of Mathematics
Carnegie Mellon University
(seeking employment for Fall 1996)