[vtkusers] Best fitting plane....

[Malcolm Drummond]

Hi Andres

Here's the execute method of a filter I wrote to return origin, dip and
dip-direction - like declination and azimuth but below horizon. It shows how
to get the result using Jacobi method. The normal you require is in
evecs[*][2].

The filter would be a lot more useful if it also returned the plane-normal
(I just derive dip and dip-dir from it), I'll add it as soon as I'm done
with financial year-end stuff - probably end of the week. Sorry I'm a bit
snowed under right now and can't check your code. I can send you all the
source for this filter if that's the way you'd like to go - you'd just have
to add one variable and a Get macro.

Regards
Malcolm

//----------------------------------------------------------

void vtkPlaneStats::Execute()
{
 vtkPolyData *pd = this->GetInput();
 vtkPoints *pts = pd->GetPoints();
 int npts = pts->GetNumberOfPoints();
 // todo - error if no points or nPoints < 3
 float sumx=0.0, sumy=0.0, sumz=0.0;
 for (int i=0; i<npts; i++)
 {
  float *pt = pts->GetPoint(i);
  sumx += pt[0];
  sumy += pt[1];
  sumz += pt[2];
 }
 float ox = sumx/npts; float oy = sumy/npts; float oz = sumz/npts;
 float sxx=0.0, syy=0.0, szz=0.0, sxy=0.0, sxz=0.0, syz=0.0;
 for (i=0; i<npts; i++)
 {
  float* pt = pts->GetPoint(i);
  float dx = ox-pt[0], dy = oy-pt[1], dz = oz-pt[2];
  sxx+=dx*dx; syy+=dy*dy; szz+=dz*dz;
  sxy+=dx*dy; sxz+=dx*dz; syz+=dy*dz;
 }
 float *mat[3],evals[3],*evecs[3];
 float mat0[3],mat1[3],mat2[3];
 float evec0[3],evec1[3],evec2[3];
 //setup
 mat[0] = mat0; mat[1] = mat1; mat[2] = mat2;
 evecs[0] = evec0; evecs[1] = evec1; evecs[2] = evec2;
 mat[0][0] = sxx; mat[0][1] = sxy; mat[0][2] = sxz;
 mat[1][0] = sxy; mat[1][1] = syy; mat[1][2] = syz;
 mat[2][0] = sxz; mat[2][1] = syz; mat[2][2] = szz;
 vtkMath::Jacobi(mat,evals,evecs);
 // vtkMath::Jacobi sorts results in dec order, vecs normalized
 // min evec is normal to plane - derive dip, dipdir
 // NB: evecs in columns
 float evec[3];
 evec[0] = evecs[0][2]; evec[1] = evecs[1][2]; evec[2] = evecs[2][2];
 // want plane normal pointing UP
 if (evec[2] < 0.0)
 {
  evec[0] *= -1;
  evec[1] *= -1;
  evec[2] *= -1;
 }
 float dip = acos(evec[2])*180.0/vtkMath::Pi();
 float dipdir = atan2(evec[1],evec[0])*180.0/vtkMath::Pi();
 this->Dip = dip;
 this->DipDir = fmod(dipdir+360.0,360.0);
 this->Origin[0] = ox;
 this->Origin[1] = oy;
 this->Origin[2] = oz;
}

----- Original Message -----
From: "Andres Barrera" <andresba at hotmail.com>
To: <malcolm at geovision.co.za>
Sent: Monday, March 29, 2004 11:28 PM
Subject: Re: [vtkusers] Best fitting plane....


> Hello Malcom,
>
>     You told me to let you know if I had any problem..... well.. here I
am!
> ;)
>
>    I tryed to imlpement what you said, but I might be doing something
> wrong... I attached the code ("Tests.cpp"), and a picture showing what I
got
> ("ViewSample.jpg"). May be it you can help me finding out the problem, and
> send me the demo code you offered me before for checking myself.
>
>     Once again, thant you very much for your time and your help.
>
>
>          Andres
>
>
>
>
> >
> >Hi Andres
> >
> >Leo van Ruijven posted a response to a similar question some time back
> >describing an eigenvector method.
> >
> >For every x,y,z calculate products xx,xy,xz,yy,yz,zz. Put sums of
products
> >in 3*3 matrix (ie. sxx,sxy,sxz;sxy,syy,syz;sxz,syz,szz) which now
> >represents
> >an ellipsoid bounding your point-set. Use vtkMath to calculate the
> >eigenvectors and eigenvalues. The two largest values are in the plane,
the
> >smallest is normal to the plane.
> >
> >Let me know if you have any problems, I used this method recently and
could
> >dig up some demo code if you need it.
> >
> >HTH
> >Malcolm
> >
> >----- Original Message -----
> >From: "Andres Barrera" <andresba at hotmail.com>
> >To: <vtkusers at vtk.org>
> >Sent: Wednesday, March 24, 2004 7:06 PM
> >Subject: [vtkusers] Best fitting plane....
> >
> >
> > > Dear list,
> > >
> > > On Fri, 22 Mar 2002 16:39:55 +0100, Anner Adrian wrote:
> > >
> > > >>Hello,
> > > >>
> > > >>is in vtk a function to find the best fit for a plane with n-points?
> > > >>I have n-points and now I want place a vtkPlaneSource so that the
> > > >>distance from the plane to the points is minimal.
> > >
> > > I need to solve a similar problem now, but I could find any answer to
> >that
> > > in the User's List archives.
> > >
> > > Is there any defined way of calculating the best fitting plane through
a
> > > series of points in vtk?
> > >
> > >    Thank you in advance
> > >
> > >          Andres
>
> _________________________________________________________________
> MSN Toolbar provides one-click access to Hotmail from any Web page - FREE
> download! http://toolbar.msn.com/go/onm00200413ave/direct/01/
>