[vtkusers] vtkDelaunay3D: wrong convex hull

[David Doria]

2011/7/26 Miguel Sotaquirá <msotaquira at gmail.com>

> Hi David!
>
> Here's the hardcoded version. Just to clarify: I just want to know if the
> current 3-D vertex ("c" in the code) is contained in the convex hull of its
> 3-D neighbors ("pts" in the code). I've used vtkDelaunay3D and vtkDelaunay2D
> currently without any success.
> Miguel
>
> #include <vtkSmartPointer.h>
> #include <vtkPolyData.h>
> #include <vtkPoints.h>
> #include <vtkUnstructuredGrid.h>
> #include <vtkDelaunay3D.h>
> #include <vtkPlane.h>
> #include <vtkDelaunay2D.h>
> #include <vtkXMLPolyDataWriter.h>
> #include <vtkXMLUnstructuredGridWriter.h>
>
> int main (int argc, char* argv[])
> {
> // Current vertex coordinates (this data is not used. serves only as a
> reference)
> double c[3] = {-0.59842, 2.12762, 5.17183};
>  // Current vertex normal
> double N[3] = {-0.604189, 0.796814, -0.00654253};
>  // Original 3-D neighbors coordinates
> double p0[3] = {-0.53592, 2.21902, 5.10933};
>  double p1[3] = {-0.59842, 2.17465, 5.10933};
> double p2[3] = {-0.53592, 2.17465, 5.17183};
>  double p3[3] = {-0.66092, 2.07811, 5.17183};
> double p4[3] = {-0.59842, 2.17465, 5.23433};
>  double p5[3] = {-0.66092, 2.12762, 5.23433};
>  // Insert 3-D neighbors in a polydata
>  vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
> pts->InsertNextPoint( p0 );
>  pts->InsertNextPoint( p1 );
> pts->InsertNextPoint( p2 );
> pts->InsertNextPoint( p3 );
>  pts->InsertNextPoint( p4 );
> pts->InsertNextPoint( p5 );
>  /*
> // 3-D convex hull: NOT WORKING
> // point cloud of neighbors
> vtkSmartPointer<vtkUnstructuredGrid> points =
> vtkSmartPointer<vtkUnstructuredGrid>::New();
>  points->SetPoints( pts );
> points->Update();
>  // convex hull
> vtkSmartPointer<vtkDelaunay3D> delaunay3D =
> vtkSmartPointer<vtkDelaunay3D>::New();
>  delaunay3D->SetInput( points );
> delaunay3D->Update();
> */
>  //
> // A 2D triangulation (using vtkDelaunay2D) will be obtained by projecting
>  // 3-D points onto a plane
> //
>  // create a plane, seting its origin in (0,0,0) and its normal N
> vtkSmartPointer<vtkPlane> plane = vtkSmartPointer<vtkPlane>::New();
>  plane->SetOrigin(0.0, 0.0, 0.0);
> plane->SetNormal(N[0], N[1], N[2]);
>  // project each point in the point cloud onto the plane
> vtkSmartPointer<vtkPoints> ptsProj = vtkSmartPointer<vtkPoints>::New();
>  double origin[3] = {0.0, 0.0, 0.0};
> double normal[3] = {N[0], N[1], N[2]};
>  double projected[3];
> double p[3];
> for ( int i = 0; i < pts->GetNumberOfPoints(); i++ ) {
>  pts->GetPoint(i,p);
> plane->ProjectPoint( p, origin, normal, projected );
>  ptsProj->InsertNextPoint(projected);
> }
>  // apply delaunay2D on the projected points
> vtkSmartPointer<vtkPolyData> ptsForTriang =
> vtkSmartPointer<vtkPolyData>::New();
>  ptsForTriang->SetPoints(ptsProj);
>  vtkSmartPointer<vtkDelaunay2D> delaunay2D =
> vtkSmartPointer<vtkDelaunay2D>::New();
>  delaunay2D->SetInput(ptsForTriang);
> delaunay2D->Update();
>
>  // export points (original and projected) and triangulation for viewing
> in Paraview
> vtkSmartPointer<vtkUnstructuredGrid> points2D =
> vtkSmartPointer<vtkUnstructuredGrid>::New();
>  vtkSmartPointer<vtkUnstructuredGrid> points3D =
> vtkSmartPointer<vtkUnstructuredGrid>::New();
> points2D->SetPoints(ptsProj);
>  points3D->SetPoints(pts);
>  vtkSmartPointer<vtkXMLUnstructuredGridWriter> writerpts2D =
> vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
>  writerpts2D->SetFileName("points2D.vtu");
> writerpts2D->SetInput(points2D);
>  writerpts2D->Write();
>  vtkSmartPointer<vtkXMLUnstructuredGridWriter> writerpts3D =
> vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
>  writerpts3D->SetFileName("points3D.vtu");
> writerpts3D->SetInput(points3D);
>  writerpts3D->Write();
>  vtkSmartPointer<vtkXMLPolyDataWriter> writer =
> vtkSmartPointer<vtkXMLPolyDataWriter>::New();
>  writer->SetFileName("delaunay2D.vtp");
> writer->SetInput(delaunay2D->GetOutput());
>  writer->Update();
>  return 0;
> }
>

You are projecting the points onto an arbitrary plane, but Delaunay2D
expected them to be on the XY plane:

"The input to this filter is a list of points specified in 3D, even though
the triangulation is 2D. Thus the triangulation is constructed in the x-y
plane, and the z coordinate is ignored (although carried through to the
output)"

I've done this by finding a local coordinate system on the plane and mapping
the points into that coordinate system:
https://github.com/daviddoria/vtkPlanarPatch/blob/master/vtkPlanarPatch.cxx

You can also do this by passing a transform to Delaunay2D with
SetTransform().

The idea is to apply the topology that results from the triangulation to
your original points.

David
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20110727/775eaa04/attachment.htm>