Re: Re: Sphere through 4 points

```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:
>
>   {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}];
>
> Out[149]= {{-2.25666, -13.9661, 7.0018}, 15.1774}
>
> Can do around 200+ per second, on my laptop.
>
> 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.

Module[{x, y, z, vars, rsq, coeffs, polys, lpolys, cen},
vars = {x, y, z};
polys = (vars.vars - rsq) /.
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}];

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.

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

```

• Prev by Date: Re: Re: How should I start with mathematica?
• Next by Date: RE: Fitting experimental data
• Previous by thread: Re: Sphere through 4 points
• Next by thread: Re: Re: Re: Sphere through 4 points