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