[vtkusers] ThinPlateSpline regression

Brehler, Michael m.brehler at dkfz-heidelberg.de
Tue Sep 2 04:46:22 EDT 2014

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();
  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];
    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();
  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->GenerateValues(1, 1, 1);
  vtkSmartPointer<vtkPolyData> polyD = contours->GetOutput();

  vtkSmartPointer<vtkXMLPolyDataWriter> surfaceWriter = vtkSmartPointer<vtkXMLPolyDataWriter>::New();

  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?


-----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

More information about the vtkusers mailing list