[vtk-developers] Problem (and proposed fix) in vtkDelaunay3D

[Wilson, Andrew T]

I believe I've found a problem in vtkDelaunay3D.  It will sometimes (often)
create tetrahedra whose points are arranged incorrectly.  This causes
problems when I'm trying to compute face normals for the cells in an
unstructured grid: if the points are in the wrong order, I can get an
inward-facing normal when I expected an outward-facing one.

In the file format specification
(http://public.kitware.com/VTK/pdf/file-formats.pdf), page 9, it seems to me
that VTK_TETRA has the following requirement.  When looking downward from
point 3, points 0, 1, and 2 must appear in counterclockwise order.  This is
not enforced in the meshes that come out of vtkDelaunay3D.

At the end of this message I've pasted in a diff against the current
revision of vtkDelaunay3D.cxx that I believe fixes the problem.  My solution
was to create the initial bounding tetrahedra with the correct point order,
then ensure that whenever a new tetrahedron is added to the mesh it gets the
order right as well.  The test results for vtkDelaunay3D
(/Graphics/Testing/Tcl/{Delaunay3D.tcl,cylMap.tcl}) look the same both with
and without this patch in place, although if there are more stringent tests
I'm not aware of them.  I've also run my own tests to examine the point
ordering in the output.

Should this patch go into the repository?

-- Andy


Index: vtkDelaunay3D.cxx
===================================================================
RCS file: /cvsroot/VTK/VTK/Graphics/vtkDelaunay3D.cxx,v
retrieving revision 1.64
diff -r1.64 vtkDelaunay3D.cxx
263,265c263,277
<       p1 = tetraPts[j];
<       p2 = tetraPts[(j+1)%4];
<       p3 = tetraPts[(j+2)%4];
---
>       // Make sure to arrange these points so that they're in
>       // counterclockwise order when viewed from the center of the
>       // cell
>       switch (j)
> 	{
> 	case 0: // face 0: points 0, 1, 2
> 	  p1 = tetraPts[0]; p2 = tetraPts[1]; p3 = tetraPts[2]; break;
> 	case 1: // face 1: points 1, 2, 3 (must flip order!)
> 	  p1 = tetraPts[1]; p2 = tetraPts[3]; p3 = tetraPts[2]; break;
> 	case 2: // face 2: points 2, 3, 0
> 	  p1 = tetraPts[2]; p2 = tetraPts[3]; p3 = tetraPts[0]; break;
> 	case 3: // face 3: points 3, 0, 1 (must flip order!)
> 	  p1 = tetraPts[3]; p2 = tetraPts[1]; p3 = tetraPts[0]; break;
> 	}
> 
488a501
> 
777,778c790,791
<   pts[0] = numPtsToInsert + 4; pts[1] = numPtsToInsert + 5; 
<   pts[2] = numPtsToInsert ; pts[3] = numPtsToInsert + 2;
---
>   pts[0] = numPtsToInsert + 4; pts[1] = numPtsToInsert + 5;
>   pts[2] = numPtsToInsert;     pts[3] = numPtsToInsert + 2;
843,847c856,866
<       //define tetrahedron
<       nodes[0] = ptId;
<       nodes[1] = this->Faces->GetId(3*tetraNum);
<       nodes[2] = this->Faces->GetId(3*tetraNum+1);
<       nodes[3] = this->Faces->GetId(3*tetraNum+2);
---
>       // Define tetrahedron.  The order of the points matters: points
>       // 0, 1, and 2 must appear in counterclockwise order when seen
>       // from point 3.  When we get here, point ptId is inside the
>       // tetrahedron whose faces we're considering and we've
>       // guaranteed that the 3 points in this face are
>       // counterclockwise wrt the new point.  That lets us create a
>       // new tetrahedron with the right ordering.
>       nodes[0] = this->Faces->GetId(3*tetraNum);
>       nodes[1] = this->Faces->GetId(3*tetraNum+1);
>       nodes[2] = this->Faces->GetId(3*tetraNum+2);
>       nodes[3] = ptId;

(end of diff)