Re: Sphere through 4 points
- To: mathgroup at smc.vnet.net
- Subject: [mg85542] Re: Sphere through 4 points
- From: dh <dh at metrohm.ch>
- Date: Wed, 13 Feb 2008 04:27:40 -0500 (EST)
- References: <fojrcq$gme$1@smc.vnet.net>
Hi Steve,
by reduction of parameters and precomputing expressions you can save
time. The following takes a list of quadruples of 3-vectors and returns
a list of centers and radius:
getCircles[v : {{_, _, _, _} ..}] :=
Module[{x1, x2, x3, y1, y2, y3, z1, z2, z3, v1, v2, v3},
subs := {x1 -> v1[[1]], x2 -> v2[[1]], x3 -> v3[[1]], y1 -> v1[[2]],
y2 -> v2[[2]], y3 -> v3[[2]], z1 -> v1[[3]], z2 -> v2[[3]],
z3 -> v3[[3]]};
( {v1, v2,
v3} = {#[[1]] - #[[4]], #[[2]] - #[[4]], #[[3]] - #[[4]]} ;
x = (x3^2 (y2 z1 - y1 z2) + x2^2 (-y3 z1 + y1 z3) +
y2^2 (-y3 z1 + y1 z3) +
z2 (y3 (x1^2 + y1^2 - y1 y3 + z1 (z1 - z2)) + y1 z2 z3 -
y1 z3^2) +
y2 (y3^2 z1 - (x1^2 + y1^2 + z1^2) z3 +
z1 z3^2))/(2 (x3 y2 z1 - x2 y3 z1 - x3 y1 z2 + x1 y3 z2 +
x2 y1 z3 - x1 y2 z3)) /. subs;
y = (x1 x3^2 z2 +
x3 (y2^2 z1 - (x1^2 + y1^2 + z1^2) z2 + z1 z2^2) +
x2^2 (x3 z1 - x1 z3) +
x2 (-(x3^2 + y3^2) z1 + (x1^2 + y1^2 + z1^2) z3 - z1 z3^2) +
x1 (y3^2 z2 - (y2^2 + z2^2) z3 + z2 z3^2))/(2 (x3 y2 z1 -
x2 y3 z1 - x3 y1 z2 + x1 y3 z2 + x2 y1 z3 - x1 y2 z3)) /.
subs;
z = (x1^2 x3 y2 + x2^2 (-x3 y1 + x1 y3) +
x3 (y1^2 y2 + y2 z1^2 - y1 (y2^2 + z2^2)) +
x2 (x3^2 y1 - y3 (x1^2 + y1^2 - y1 y3 + z1^2) + y1 z3^2) -
x1 (x3^2 y2 - y3 (y2^2 - y2 y3 + z2^2) +
y2 z3^2))/(2 (x3 y2 z1 - x2 y3 z1 - x3 y1 z2 + x1 y3 z2 +
x2 y1 z3 - x1 y2 z3)) /. subs;
r = Sqrt[x^2 + y^2 + z^2];
{{x, y, z} + #[[4]], r}
) & /@ v
]
here we compute 1000 circles:
dat=RandomReal[{-1,1},{1000,4,3}];
Timing[getCircles @ dat; ]
{1.188,Null}
hope this helps, Daniel
Steve Gray wrote:
> Given 4 points in 3-space not lying in a plane, I want the
> center and radius of the incident sphere. I know about the elegant
> determinant solution (4 or 5 determinants, each 4x4, plus a few
> arithmetic operations) and I'm using that, but maybe someone has coded
> a faster version.
> Thanks for any information.
>
> Steve Gray
>