# [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[])
{

vtkPoints* points = polydata->GetPoints();

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;
points->GetPoint(i,p);

x[i] = pow(p,2); //x^2
x[i] = pow(p,2); //y^2
x[i] = p * p; //x*y
x[i] = p; // x
x[i] = p; //y
x[i] = 1; // constant
y[i] = p; //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] << " ";
}
std::cout << std::endl;

// create the quadric function definition
quadric->SetCoefficients(m, m, 0, m, 0, 0, m, m, 0, m);

vtkSmartPointer<vtkSampleFunction> sample = vtkSmartPointer<vtkSampleFunction>::New();
sample->SetSampleDimensions(50,50,50);
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->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
```