Re: 3D data set
- To: mathgroup at smc.vnet.net
- Subject: [mg51077] Re: 3D data set
- From: koopman at sfu.ca (Ray Koopman)
- Date: Mon, 4 Oct 2004 06:17:35 -0400 (EDT)
- References: <cjdq0q$ar3$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Giovanni Bellesia <giovanni.bellesia at ucd.ie> wrote in message
news:<cjdq0q$ar3$1 at smc.vnet.net>...
> [...]
> Does anybody knows a clear and efficient procedure to fit these points
> to a regular, circular helix.
> [...]
To the extent that fitting a cylinder prior to fitting a helix is not
a waste of time, the following may be of some use.
Given n points in two dimensions, {{x1,y1},...,{xn,yn}}, the following
routine fits a circle thru them by minimizing the variance of the
squared distances of the points from the putative center. It returns
{variance, center}. It makes no attempt to get an explicit estimate of
the radius, but the way the center is found implicitly estimates the
radius as the rms distance from the center.
fitcircle[xy_] := Block[{m = -Mean@xy, x,y, u,v},
{x,y} = Transpose@xy + m;
{u,v} = -.5 LinearSolve[{x,y}.xy,{x,y}.(x^2 + y^2)];
{Variance[(x + u)^2 + (y + v)^2], -({u,v} + m)}]
To fit a cylinder to a set of points in three dimensions, postmultiply
the data {{x1,y1,z1},...,{xn,yn,zn}} by a 3 x 3 rotation matrix to
get a new matrix in which one of the columns -- the third, say -- is
supposed to contain the projections of the points onto the axis of the
cylinder. Then the first two columns contain the projections of the
points onto the plane that is orthogonal to the axis, and those
projections should lie on a circle in that plane.
The problem is to find the rotation matrix that minimizes the variance
returned by fitcircle when it is given the first two columns of the
rotated matrix. In principle this is something that FindMinimum should
be able to do, but so far I have been unable to coax a solution out of
it. My fallback has been to use a sequence of ContourPlots, specifying
successively narrower ranges for the angles of rotation.
T13[angle_] := With[{c = N@Cos@angle, s = N@Sin@angle},
{{c,0,-s},{0,1.,0},{s,0,c}}]
T23[angle_] := With[{c = N@Cos@angle, s = N@Sin@angle},
{{1.,0,0},{0,c,-s},{0,s,c}}]
ContourPlot[-Sqrt@First@fitcircle[xyz.(T23[a].T13[b])[[All,{1,2}]]],
{a,-Pi/2,Pi/2},{b,-Pi/2,Pi/2}, PlotPoints->50];
I change sign so as to be looking for white, not black.
I use Sqrt because it gives sharper peaks initially.
I adjust PlotPoints as necessary.