[vtkusers] planar points to get a plane and its normal vector

[David Doria]

I'm not sure why there isn't a function for this (that I can find at
least!). You can find a function that does what you want here:

> http://www.vtk.org/pipermail/vtkusers/2004-March/072916.html
> or here
>
> http://markmail.org/message/52odxars4ijopqfe#query:fit%20plane%20vtk+page:1+mid:52odxars4ijopqfe+state:results
>
> The idea is that the eigen vector corresponding to the smallest eigen value
> of the scatter matrix of the points is the normal of the plane. It seems
> like there should be a function somewhere in vtk:
> vtkPlane* FitPlane(vtkPoints* Points);
>
> GC, when you get this working how you want it can you post it here and
> maybe someone can add it somewhere in vtk?
>
> Thanks,
>
> David
>

I extracted this function mostly from vtkDelaunay2d. Maybe someone can add
this to vtkPlane ?

vtkPlane* BestFitPlane(vtkPoints *points)
    {
        vtkIdType NumPoints = points->GetNumberOfPoints();

        //find the center of mass of the points
        double Center[3] = {0.0, 0.0, 0.0};

        for(vtkIdType i = 0; i < NumPoints; i++)
        {
            double point[3];
            points->GetPoint(i, point);

            Center[0] += point[0];
            Center[1] += point[1];
            Center[2] += point[2];
        }

        Center[0] = Center[0]/static_cast<double>(NumPoints);
        Center[1] = Center[1]/static_cast<double>(NumPoints);
        Center[2] = Center[2]/static_cast<double>(NumPoints);

        //Compute Sample Covariance Matrix (with points centered at the
origin)

        double *a[3], a0[3], a1[3], a2[3];
        a[0] = a0; a[1] = a1; a[2] = a2;
        for(unsigned int i = 0; i < 3; i++)
        {
            a0[i] = a1[i] = a2[i] = 0.0;
        }

        for(unsigned int pointId = 0; pointId < NumPoints; pointId++ )
        {
            double x[3], xp[3];
            points->GetPoint(pointId, x);
            xp[0] = x[0] - Center[0];
            xp[1] = x[1] - Center[1];
            xp[2] = x[2] - Center[2];
            for (unsigned int i = 0; i < 3; i++)
            {
                a0[i] += xp[0] * xp[i];
                a1[i] += xp[1] * xp[i];
                a2[i] += xp[2] * xp[i];
            }
        }

        for(unsigned int i = 0; i < 3; i++)
        {
            a0[i] /= static_cast<double>(NumPoints);
            a1[i] /= static_cast<double>(NumPoints);
            a2[i] /= static_cast<double>(NumPoints);
        }

        // Extract eigenvectors from covariance matrix
        double *v[3], v0[3], v1[3], v2[3];
        v[0] = v0; v[1] = v1; v[2] = v2;
        double eigval[3];
        vtkMath::Jacobi(a,eigval,v);
        //Jacobi iteration for the solution of eigenvectors/eigenvalues of a
3x3 real symmetric matrix. Square 3x3 matrix a; output eigenvalues in w; and
output eigenvectors in v. Resulting eigenvalues/vectors are sorted in
decreasing order; eigenvectors are normalized.

        vtkPlane* BestPlane = vtkPlane::New();
        //Set the plane normal to the smallest eigen vector
        BestPlane->SetNormal(v2[0], v2[1], v2[2]);

        //Set the plane origin to the center of mass
        BestPlane->SetOrigin(Center[0], Center[1], Center[2]);

        return BestPlane;
    }





and a test call:

    vtkPoints* Points = vtkPoints::New();

    //create 3 known points
    Points->InsertNextPoint(0.0, 0.0, 0.0);
    Points->InsertNextPoint(1.0, 0.0, 0.0);
    Points->InsertNextPoint(0.0, 1.0, 0.0);

    vtkPlane* BestPlane = VTKHelpers::BestFitPlane(Points);

    std::cout << "Best Plane:\n" << *BestPlane << "\n";



Thanks,

David
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20090707/5ad31a79/attachment.htm>