[vtkusers] ThinPlateSpline regression

[Brehler, Michael]

Hi David,

thank you for your reply. I found this (broken) example http://www.vtk.org/Wiki/VTK/Examples/Broken/QuadricSurfaceFitting which some time ago probably did what I want to do. I adjusted the main method to work with the VTK6 input and output and added the correct matrices to be freed:

 int main (int argc, char *argv[])
{
  vtkSmartPointer<vtkXMLPolyDataReader> reader = vtkSmartPointer<vtkXMLPolyDataReader>::New();
  reader->SetFileName(argv[1]);
  reader->Update();
 
  vtkPolyData* polydata = reader->GetOutput();
  vtkPoints* points = polydata->GetPoints();
 
  //fit a quadric surface
  int numberOfSamples = points->GetNumberOfPoints();
  int numberOfVariables = 6;

  double **x = create_matrix<double> (numberOfSamples, numberOfVariables);
  double **y = create_matrix<double> ( numberOfSamples, 1 );
 
  for(unsigned int i = 0; i < numberOfSamples; i++)
    {
    double p[3];
    points->GetPoint(i,p);
 
    x[i][0] = pow(p[0],2); //x^2
    x[i][1] = pow(p[1],2); //y^2
    x[i][2] = p[0] * p[1]; //x*y
    x[i][3] = p[0]; // x
    x[i][4] = p[1]; //y
    x[i][5] = 1; // constant
    y[i][0] = p[2]; //z
  }
 
  double **m = create_matrix<double> ( numberOfVariables, 1 );
  vtkMath::SolveLeastSquares(numberOfSamples, x, numberOfVariables, y, 1, m);
 
  std::cout << "Solution is: " << std::endl;
  for(unsigned int i = 0; i < numberOfVariables; i++)
    {
    std::cout << m[i][0] << " ";
    }
  std::cout << std::endl;
 
  // create the quadric function definition
  vtkSmartPointer<vtkQuadric> quadric = vtkSmartPointer<vtkQuadric>::New();
  quadric->SetCoefficients(m[0][0], m[1][0], 0, m[2][0], 0, 0, m[3][0], m[4][0], 0, m[5][0]);
 
  // sample the quadric function
  vtkSmartPointer<vtkSampleFunction> sample = vtkSmartPointer<vtkSampleFunction>::New();
  sample->SetSampleDimensions(50,50,50);
  sample->SetImplicitFunction(quadric);
  double xmin = -5, xmax=5, ymin=-5, ymax=5, zmin=-5, zmax=5;
  sample->SetModelBounds(xmin, xmax, ymin, ymax, zmin, zmax);

  //create the 0 isosurface
  vtkSmartPointer<vtkContourFilter> contours = vtkSmartPointer<vtkContourFilter>::New();
  contours->SetInputConnection(sample->GetOutputPort());
  contours->GenerateValues(1, 1, 1);
  contours->Update();
  vtkSmartPointer<vtkPolyData> polyD = contours->GetOutput();

  vtkSmartPointer<vtkXMLPolyDataWriter> surfaceWriter = vtkSmartPointer<vtkXMLPolyDataWriter>::New();
  surfaceWriter->SetFileName("Quadratic.vtp");
  surfaceWriter->SetInputData(polyD);
  surfaceWriter->Write();  
 
  free_matrix(x);
  free_matrix(y);
  free_matrix(m);

  return 0;
}

But using the input points from the example page the output is a elliptic cylinder centered around the input points, which is not the desired result. Could somebody tell me what went wrong?

Best,
Michael


-----Ursprüngliche Nachricht-----
Von: David Gobbi [mailto:david.gobbi at gmail.com] 
Gesendet: Freitag, 29. August 2014 16:58
An: Brehler, Michael
Cc: vtkusers at vtk.org
Betreff: Re: [vtkusers] ThinPlateSpline regression

Hi Michael,

No, I don't think that ThinPlateSplineTransform provides a way to do this.
I see a ComputeBestFittingPlane() method in vtkDelaunay2D, but it is a protected method and, looking at the method that it uses, I don't think it computes the "best" plane, anyway.

Hopefully other people who have used VTK to compute the best-fit plane can give you a better answer, but my suggestion is to search the internet for a method that computes the best plane via an eigenvalue decomposition (my suspicion is that for the best-fit plane, the eigenvector for the smallest eigenvalue is the normal of the plane, and the "origin" of the plane is the centre-of-mass).

The vtkMath::Jacobi() function can compute the eigenvalues.

 - David


On Fri, Aug 29, 2014 at 3:29 AM, Brehler, Michael <m.brehler at dkfz-heidelberg.de> wrote:
> Hi,
>
>
>
> I am using VTK to find a best fitting hyperplane (or regression plane) 
> for a given 3D point cloud (polydata surface) and was wondering if the 
> ThinPlateSplineTransform could be of any help. Can I use an equally 
> distributed 'grid' of points to calculate the matching hyperplane with 
> the ThinPlateSplineTransform or is there a better way?
>
>
>
> Regards,
>
> Michael