Re: [mg 1129] Find the center!
- To: mathgroup at christensen.cybernetics.net
- Subject: [mg1170] Re: [mg1120] [mg 1129] Find the center!
- From: Allan Hayes <hay at haystack.demon.co.uk>
- Date: Sat, 20 May 1995 03:33:27 -0400
Gottfried Mayer-Kress <gmk at pegasos.ccsr.uiuc.edu> wrote in
[mg1120] Find the center!
> There are a bunch of points in 3-space and they are supposed to be
> scattered around the surface of a sphere. You want to find the
> center of the sphere and the radius.
Jorma Virtamo ([mg 1129]) gave a function for doing this.
In[1]:=
(*Virtamo*)
center[pts_] := Module[ {R,r=Plus@@pts/Length[pts]},
{ R = Inverse[ Plus@@(Outer[Times,#-r,#-r]& /@ pts) ].
Plus@@(#.#(#-r)& /@ pts)/2,
Sqrt[Plus@@((#-R).(#-R)& /@ pts)/Length[pts]] } ]
We can use Fit and the standard equation for a circle to get a
different function.
In[2]:=
CenterRad[pts_]:=
Block[{data, f,g,h,c},
data = {Sequence@@#,#.#}&/@pts;
{c,f,g,h} = List@@Fit[data, {1,x,y,z}, {x,y,z}]/{1,2x,2y,2z};
{{f,g,h}, Sqrt[f^2+g^2+h^2 +c]}
]
This is faster than Virtamo's function on the test that he uses
TEST
First generate some points, roughly on a sphers center {1,2,3}, radius 1.
In[3]:=
rndnrm[x_] := (.98+.04Random[])x/Sqrt[x.x]
In[4]:=
points = ({1,2,3}+#)& /@ rndnrm /@ Table[Random[]-.5,{10},{3}];
Now try to find the center and radius
In[5]:=
CenterRad[points]
Do[CenterRad[points],{10}]//Timing//First
Out[5]=
{{0.995855, 1.98329, 3.00476}, 1.00521}
Out[6]=
0.55 Second
In[7]:=
center[points]
Do[center[points],{10}]//Timing//First
Out[7]=
{{0.995855, 1.98329, 3.00476}, 1.00521}
Out[8]=
1.46667 Second
Allan Hayes
hay at haystack.demon.co.uk