MathGroup Archive 2004

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

Search the Archive

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.


  • Prev by Date: Re: Problem with Maximize and conditions.
  • Next by Date: Re: Derivatives of user-defined control-flow functions
  • Previous by thread: Re: 3D data set
  • Next by thread: Re: Partial derivation with regards to variables in a list