Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Re: Sphere through 4 points

• To: mathgroup at smc.vnet.net
• Subject: [mg85539] Re: [mg85415] Re: [mg85401] Sphere through 4 points
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Wed, 13 Feb 2008 04:26:08 -0500 (EST)
• References: <200802090915.EAA16803@smc.vnet.net> <200802101010.FAA17539@smc.vnet.net>

danl at wolfram.com 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
>>
> 
> 
> Not sure if this is faster, but could do:
> 
> centerAndRadius[data_] := Module[
>   {x, y, z, vars, rsq, coeffs, polys, lpolys, cen},
>   vars = {x, y, z};
>   polys = (vars.vars - rsq) /.
>     Thread[vars -> (vars - Transpose[data])];
>   lpolys = Expand[Rest[polys] - First[polys]];
>   coeffs = CoefficientArrays[lpolys, vars];
>   cen = LinearSolve[Last[coeffs], First[coeffs]];
>   {cen, Sqrt[rsq]} /.
>    First[Solve[(First[polys] /. Thread[vars -> cen]) == 0]]
>   ]
> 
> In[147]:= SeedRandom[1111];
> data = RandomReal[{-5, 5}, {4, 3}];
> centerAndRadius[data]
> 
> Out[149]= {{-2.25666, -13.9661, 7.0018}, 15.1774}
> 
> Can do around 200+ per second, on my laptop.
> 
> In[150]:= Timing[Do[centerAndRadius[data], {100}]]
> Out[150]= {0.421, Null}
> 
> Daniel Lichtblau
> Wolfram Research

It was pointed out to me that this code does not work correctly. Ther 
eason is I forgot to negate the constant terms, in effect moving them to 
the other side to make equations from linear polynomials. Here is the 
correct form.

centerAndRadius[data_] :=
  Module[{x, y, z, vars, rsq, coeffs, polys, lpolys, cen},
   vars = {x, y, z};
   polys = (vars.vars - rsq) /.
     Thread[vars -> (vars - Transpose[data])];
   lpolys = Expand[Rest[polys] - First[polys]];
   coeffs = CoefficientArrays[lpolys, vars];
   cen = LinearSolve[Last[coeffs], -First[coeffs]];
   {cen, Sqrt[rsq]} /.
    First[Solve[(First[polys] /. Thread[vars -> cen]) == 0]]]

Check:

In[85]:= SeedRandom[1111];
data = RandomReal[{-5, 5}, {4, 3}];
{cen, rad} = centerAndRadius[data]

Out[87]= {{2.25666, 13.9661, -7.0018}, 17.9779}

In[89]:= Map[Norm[# - cen] - rad &, data]
Out[89]= {0., 0., 0., 0.}

If raw speed is the issue, the variant below is faster.

centerAndRadius2[data_] :=
  Module[{x, y, z, x0, y0, z0, rsqr, coeffs, rules, polys, lpolys, cen},
   rules = (Thread[{x, y, z} -> #1] &) /@ data;
   polys = x^2 - 2*x*x0 + y^2 - 2*y*y0 + z^2 - 2*z*z0 /. rules;
   lpolys = Rest[polys] - First[polys];
   coeffs = Normal[CoefficientArrays[lpolys, {x0, y0, z0}]];
   cen = LinearSolve[Last[coeffs], -First[coeffs]];
   {cen, Norm[cen - First[data]]}]

Daniel Lichtblau
Wolfram Research

• Follow-Ups:
  • Re: Re: Re: Sphere through 4 points
    • From: Daniel Lichtblau <danl@wolfram.com>

• References:
  • Sphere through 4 points
    • From: Steve Gray <stevebg@roadrunner.com>
  • Re: Sphere through 4 points
    • From: danl@wolfram.com