Used 200grit to polish the code (bugs fixes)

A series of bug fixes including: 1) gradient calculation error;
2) addressing boundary edges interpolation in certain rare cases;
3) setting scalar values (if requested) in one fill_n operation;
4) correct treatment of normals.

Also the documentation was improved, and a new .cxx test case was
added.

Filters/Core/Testing/Cxx/CMakeLists.txt

@@ -16,6 +16,7 @@ vtk_add_test_cxx(${vtk-module}CxxTests tests
   TestDelaunay3D.cxx,NO_VALID
   TestExecutionTimer.cxx,NO_VALID
   TestFeatureEdges.cxx,NO_VALID
+  TestFlyingEdges.cxx
   TestGlyph3D.cxx
   TestHedgeHog.cxx,NO_VALID
   TestImplicitPolyDataDistance.cxx

Filters/Core/Testing/Cxx/TestFlyingEdges.cxx

+/*=========================================================================
+
+  Program:   Visualization Toolkit
+  Module:    TestFlyingEdges.cxx
+
+  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
+  All rights reserved.
+  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
+
+  This software is distributed WITHOUT ANY WARRANTY; without even
+  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+  PURPOSE.  See the above copyright notice for more information.
+
+=========================================================================*/
+
+// Description
+// This test creates a wavelet dataset and creates isosurfaces using
+// vtkFlyingEdges3D
+
+#include "vtkActor.h"
+#include "vtkFlyingEdges3D.h"
+#include "vtkNew.h"
+#include "vtkPolyDataMapper.h"
+#include "vtkRegressionTestImage.h"
+#include "vtkRenderer.h"
+#include "vtkRenderWindow.h"
+#include "vtkRenderWindowInteractor.h"
+#include "vtkRTAnalyticSource.h"
+#include "vtkTesting.h"
+
+int TestFlyingEdges(int argc, char *argv[])
+{
+  cout << "CTEST_FULL_OUTPUT (Avoid ctest truncation of output)" << endl;
+
+  // Create the sample dataset
+  vtkNew<vtkRTAnalyticSource> wavelet;
+  wavelet->SetWholeExtent(-63, 64,
+                          -63, 64,
+                          -63, 64);
+  wavelet->SetCenter(0.0, 0.0, 0.0);
+  wavelet->Update();
+
+  vtkNew<vtkFlyingEdges3D> flyingEdges;
+  flyingEdges->SetInputConnection(wavelet->GetOutputPort());
+  flyingEdges->GenerateValues(54, 128.0, 225.0);
+  flyingEdges->ComputeNormalsOn();
+  flyingEdges->ComputeGradientsOn();
+  flyingEdges->ComputeScalarsOn();
+  flyingEdges->SetArrayComponent(0);
+
+  vtkNew<vtkPolyDataMapper> mapper;
+  mapper->SetInputConnection(flyingEdges->GetOutputPort());
+  vtkNew<vtkActor> actor;
+  actor->SetMapper(mapper.GetPointer());
+  vtkNew<vtkRenderer> ren;
+  ren->AddActor(actor.GetPointer());
+
+  vtkNew<vtkRenderWindow> renWin;
+  renWin->SetSize(399, 401);
+  renWin->SetMultiSamples(0);
+  renWin->AddRenderer(ren.GetPointer());
+  vtkNew<vtkRenderWindowInteractor> iren;
+  iren->SetRenderWindow(renWin.GetPointer());
+
+  ren->ResetCamera();
+  renWin->Render();
+  iren->Initialize();
+
+  int retVal = vtkRegressionTestImage(renWin.GetPointer());
+  if (retVal == vtkRegressionTester::DO_INTERACTOR)
+    {
+    iren->Start();
+    }
+  return !(retVal == vtkTesting::DO_INTERACTOR);
+}

Filters/Core/Testing/Data/Baseline/TestFlyingEdges.png.md5

+9f573000b865ab2da5383a392c32118b

Filters/Core/Testing/Data/Baseline/TestFlyingEdges3D.png.md5

-4f0d7b8f24a4f3f3b7a7f2485409fd93
+4650004e30a3cd6767079b7b64829659

Filters/Core/Testing/Python/TestFlyingEdges2D.py

-
+#!/usr/bin/env python
 import vtk
 from vtk.test import Testing
 from vtk.util.misc import vtkGetDataRoot
@@ -47,3 +47,4 @@ iren.Initialize()
 
 renWin.Render()
 # --- end of script --
+iren.Start()

Filters/Core/Testing/Python/TestFlyingEdges3D.py

-
+#!/usr/bin/env python
 import vtk
 from vtk.test import Testing
 from vtk.util.misc import vtkGetDataRoot
@@ -18,15 +18,17 @@ sphere = vtk.vtkSphere()
 sphere.SetCenter( 0.0,0.0,0.0)
 sphere.SetRadius(0.25)
 
-# iso-surface to create geometry
+# iso-surface to create geometry. This uses a very small volume to stress
+# boundary conditions in vtkFlyingEdges; i.e., we want the sphere isosurface
+# to intersect the boundary.
 sample = vtk.vtkSampleFunction()
 sample.SetImplicitFunction(sphere)
 sample.SetModelBounds(-0.5,0.5, -0.5,0.5, -0.5,0.5)
-sample.SetSampleDimensions(25,35,40)
+sample.SetSampleDimensions(3,3,3)
 
 iso = vtk.vtkFlyingEdges3D()
 iso.SetInputConnection(sample.GetOutputPort())
-iso.GenerateValues(3,-.11,.11)
+iso.SetValue(0,0.25)
 iso.ComputeNormalsOn()
 iso.ComputeGradientsOn()
 
@@ -36,7 +38,7 @@ isoMapper.ScalarVisibilityOff()
 isoActor = vtk.vtkActor()
 isoActor.SetMapper(isoMapper)
 isoActor.GetProperty().SetColor(1,1,1)
-isoActor.GetProperty().SetOpacity(0.5)
+isoActor.GetProperty().SetOpacity(1)
 
 outline = vtk.vtkOutlineFilter()
 outline.SetInputConnection(sample.GetOutputPort())

Filters/Core/vtkFlyingEdges2D.cxx

@@ -106,7 +106,7 @@ public:
   // The three passes of the algorithm.
   void ProcessXEdge(double value, T* inPtr, vtkIdType row); //PASS 1
   void ProcessYEdges(vtkIdType row); //PASS 2
-  void GenerateOutput(double value, T* inPtr, vtkIdType row); //PASS 3
+  void GenerateOutput(double value, T* inPtr, vtkIdType row); //PASS 4
 
   // Place holder for now in case fancy bit fiddling is needed later.
   void SetXEdge(unsigned char *ePtr, unsigned char edgeCase)
@@ -230,10 +230,10 @@ public:
           }//for all rows in this batch
         }
     };
-  template <class TT> class Pass3
+  template <class TT> class Pass4
     {
     public:
-      Pass3(vtkFlyingEdges2DAlgorithm<TT> *algo, double value)
+      Pass4(vtkFlyingEdges2DAlgorithm<TT> *algo, double value)
         { this->Algo = algo; this->Value = value;}
       vtkFlyingEdges2DAlgorithm<TT> *Algo;
       double Value;
@@ -525,8 +525,8 @@ ProcessYEdges(vtkIdType row)
 }
 
 //----------------------------------------------------------------------------
-// PASS 3: Process the x-row cells to generate output primitives, including
-// point coordinates and line segments. This is the third pass of the
+// PASS 4: Process the x-row cells to generate output primitives, including
+// point coordinates and line segments. This is the fourth pass of the
 // algorithm.
 template <class T> void vtkFlyingEdges2DAlgorithm<T>::
 GenerateOutput(double value, T* rowPtr, vtkIdType row)
@@ -549,14 +549,6 @@ GenerateOutput(double value, T* rowPtr, vtkIdType row)
   ePtr0 = this->XCases + row*(this->Dims[0]-1) + xL;
   ePtr1 = ePtr0 + this->Dims[0]-1;
 
-  // Update scalars along this x-row if necessary
-  vtkIdType numNewPts = eMD1[0] - eMD0[0];
-  if ( this->NewScalars && numNewPts > 0 )
-    {
-    T TValue = static_cast<T>(value);
-    std::fill_n(this->NewScalars+eMD0[0], numNewPts, TValue);
-    }
-
   // Traverse all pixels in this row, those containing the contour are
   // further identified for processing, meaning generating points and
   // triangles. Begin by setting up point ids on pixel edges.
@@ -751,11 +743,13 @@ ContourImage(vtkFlyingEdges2D *self, T *scalars, vtkPoints *newPts,
       {
       newScalars->WriteVoidPointer(0,numOutXPts+numOutYPts);
       algo.NewScalars = static_cast<T*>(newScalars->GetVoidPointer(0));
+      T TValue = static_cast<T>(value);
+      std::fill_n(algo.NewScalars, numOutXPts+numOutYPts, TValue);
       }
 
-    // Now process each x-row and produce the output primitives.
-    Pass3<T> pass3(&algo,value);
-    vtkSMPTools::For(0,algo.Dims[1]-1, pass3);
+    // PASS 4: Now process each x-row and produce the output primitives.
+    Pass4<T> pass4(&algo,value);
+    vtkSMPTools::For(0,algo.Dims[1]-1, pass4);
 
     // Handle multiple contours
     startXPts = numOutXPts;

Filters/Core/vtkFlyingEdges2D.h

@@ -25,19 +25,28 @@
 // which partial computational results can be used to eliminate future
 // computations.
 //
-// This is a three-pass algorithm. The first pass processes all x-edges and
+// This is a four-pass algorithm. The first pass processes all x-edges and
 // builds x-edge case values (which, when the two x-edges defining a pixel
 // are combined, are equivalent to vertex-based case table except edge-based
 // approaches are separable to parallel computing). Next x-pixel rows are
 // processed to gather information from y-edges (basically to count the
-// number of edge intersections and lines generated). Finally in the
-// third pass output primitives are generated into pre-allocated arrays. This
-// implementation uses pixel cell axes (a x-y dyad located at the pixel
-// origin) to ensure that each edge is intersected at most one time.
+// number of edge intersections and lines generated). In the third pass a
+// prefix sum is used to count and allocate memory for the output
+// primitives. Finally in the fourth pass output primitives are generated into
+// pre-allocated arrays. This implementation uses pixel cell axes (a x-y dyad
+// located at the pixel origin) to ensure that each edge is intersected at
+// most one time.
+//
+// See the paper "Flying Edges: A High-Performance Scalable Isocontouring
+// Algorithm" by Schroeder, Maynard, Geveci. Proc. of LDAV 2015. Chicago, IL.
 
 // .SECTION Caveats
 // This filter is specialized to 2D images. This implementation can produce
 // degenerate line segments (i.e., zero-length line segments).
+//
+// This class has been threaded with vtkSMPTools. Using TBB or other
+// non-sequential type (set in the CMake variable
+// VTK_SMP_IMPLEMENTATION_TYPE) may improve performance significantly.
 
 // .SECTION See Also
 // vtkContourFilter vtkFlyingEdges3D vtkSynchronizedTemplates2D

Filters/Core/vtkFlyingEdges3D.cxx

@@ -212,6 +212,7 @@ public:
       x[0] = x0[0] + t*(x1[0]-x0[0]);
       x[1] = x0[1] + t*(x1[1]-x0[1]);
       x[2] = x0[2] + t*(x1[2]-x0[2]);
+
       if ( this->NeedGradients )
         {
         float gTmp[3], g1[3];
@@ -337,10 +338,10 @@ public:
           }//for all slices in this batch
         }
     };
-  template <class TT> class Pass3
+  template <class TT> class Pass4
     {
     public:
-      Pass3(vtkFlyingEdges3DAlgorithm<TT> *algo, double value)
+      Pass4(vtkFlyingEdges3DAlgorithm<TT> *algo, double value)
         {this->Algo = algo; this->Value = value;}
       vtkFlyingEdges3DAlgorithm<TT> *Algo;
       double Value;
@@ -460,14 +461,10 @@ vtkFlyingEdges3DAlgorithm():XCases(NULL),EdgeMetaData(NULL),NewScalars(NULL),
             *edgeCase++ = numTris;
             for ( edge = triCase->edges; edge[0] > -1; edge += 3, edgeCase+=3 )
               {
-              // Build new case table. You're probably wondering why the
-              // crazy (0,2,1) edge order below. Simple: as originally
-              // presented the MC algorithm used a left-handed coordinate
-              // system, so we have to reverse the ordering of the triangle
-              // to make it consistent with any generated normals.
+              // Build new case table.
               edgeCase[0] = this->EdgeMap[edge[0]];
-              edgeCase[2] = this->EdgeMap[edge[1]];
-              edgeCase[1] = this->EdgeMap[edge[2]];
+              edgeCase[1] = this->EdgeMap[edge[1]];
+              edgeCase[2] = this->EdgeMap[edge[2]];
               }
             }
           }//x-edges
@@ -649,6 +646,7 @@ InterpolateEdge(double value, vtkIdType ijk[3],
   xPtr[0] = x0[0] + t*(x1[0]-x0[0]);
   xPtr[1] = x0[1] + t*(x1[1]-x0[1]);
   xPtr[2] = x0[2] + t*(x1[2]-x0[2]);
+
   if ( this->NeedGradients )
     {
     float gTmp[3], g0[3], g1[3];
@@ -700,6 +698,7 @@ GeneratePoints(double value, unsigned char loc, vtkIdType ijk[3],
                           g0);
     }
 
+  // Interpolate the cell axes edges
   for(int i=0; i < 3; ++i)
     {
     if(edgeUses[i*4])
@@ -708,7 +707,7 @@ GeneratePoints(double value, unsigned char loc, vtkIdType ijk[3],
       //edgesUses[4] == y axes edge
       //edgesUses[8] == z axes edge
       float x1[3] = {x[0], x[1], x[2] }; x1[i] += this->Spacing[i];
-      vtkIdType ijk1[3] = { ijk[0], ijk[1], ijk[2] }; ++ijk[i];
+      vtkIdType ijk1[3] = { ijk[0], ijk[1], ijk[2] }; ++ijk1[i];
 
       T const * const sPtr2 = (sPtr+incs[i]);
       double t = (value - *sPtr) / (*sPtr2 - *sPtr);
@@ -716,40 +715,43 @@ GeneratePoints(double value, unsigned char loc, vtkIdType ijk[3],
       }
     }
 
-  // Otherwise do more general gyrations. These are boundary situations where
-  // the voxel axes is not fully formed. These situations occur on the
-  // +x,+y,+z volume boundaries. (The other cases are handled by the default:
-  // case and are expected.)
-  switch (loc) //location is one of 27 regions in the volume
+  // On the boundary cells special work has to be done to cover the partial
+  // cell axes. These are boundary situations where the voxel axes is not
+  // fully formed. These situations occur on the +x,+y,+z volume
+  // boundaries. (The other cases fall through the default: case which is
+  // expected.)
+  //
+  // Note that loc is one of 27 regions in the volume, with (0,1,2)
+  // indicating (interior, min, max) along coordinate axes.
+  switch (loc)
     {
-    case 2: case 6: case 18:
-    case 22: case 26: //+x & +x -y & +x -z & +x -y -z +x +y -z
+    case 2: case 6: case 18: case 22: //+x
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 5, edgeUses, eIds);
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 9, edgeUses, eIds);
       break;
-    case 8: case 24: case 25: //+y & +y -z & +y -x -z
+    case 8: case 9: case 24: case 25: //+y
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 1, edgeUses, eIds);
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 10, edgeUses, eIds);
       break;
-    case 10://+x +y
+    case 32: case 33: case 36: case 37: //+z
+      this->InterpolateEdge(value, ijk, sPtr, incs, x, 2, edgeUses, eIds);
+      this->InterpolateEdge(value, ijk, sPtr, incs, x, 6, edgeUses, eIds);
+      break;
+    case 10: case 26: //+x +y
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 1, edgeUses, eIds);
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 5, edgeUses, eIds);
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 9, edgeUses, eIds);
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 10, edgeUses, eIds);
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 11, edgeUses, eIds);
       break;
-    case 32: case 33: case 36: case 37: //+z & -x +z & -y +z & -x -y +z
-      this->InterpolateEdge(value, ijk, sPtr, incs, x, 2, edgeUses, eIds);
-      this->InterpolateEdge(value, ijk, sPtr, incs, x, 6, edgeUses, eIds);
-      break;
-    case 34: case 38: //+x +z & +x -y +z
+    case 34: case 38: //+x +z
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 2, edgeUses, eIds);
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 5, edgeUses, eIds);
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 9, edgeUses, eIds);
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 6, edgeUses, eIds);
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 7, edgeUses, eIds);
       break;
-    case 9: case 40: case 41: //-x +y & +y +z & -x + y + z
+    case 40: case 41: //+y +z
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 1, edgeUses, eIds);
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 2, edgeUses, eIds);
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 3, edgeUses, eIds);
@@ -767,7 +769,7 @@ GeneratePoints(double value, unsigned char loc, vtkIdType ijk[3],
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 6, edgeUses, eIds);
       this->InterpolateEdge(value, ijk, sPtr, incs, x, 7, edgeUses, eIds);
       break;
-    default: //interior, or -x,-y,-z boundary
+    default: //interior, or -x,-y,-z boundaries
       return;
     }
 }
@@ -941,8 +943,9 @@ ProcessYZEdges(vtkIdType row, vtkIdType slice)
 }
 
 //----------------------------------------------------------------------------
-// PASS 3: Process the x-row cells to generate output primitives, including
-// point coordinates and triangles. This is the third pass of the algorithm.
+// PASS 4: Process the x-row cells to generate output primitives, including
+// point coordinates and triangles. This is the fourth and final pass of the
+// algorithm.
 template <class T> void vtkFlyingEdges3DAlgorithm<T>::
 GenerateOutput(double value, T* rowPtr, vtkIdType row, vtkIdType slice)
 {
@@ -976,14 +979,6 @@ GenerateOutput(double value, T* rowPtr, vtkIdType row, vtkIdType slice)
   ePtr[2] = ePtr[0] + this->SliceOffset;
   ePtr[3] = ePtr[2] + this->Dims[0]-1;
 
-  // Update scalars along this x-row if necessary
-  vtkIdType numNewPts = eMD[1][0] - eMD[0][0];
-  if ( this->NewScalars && numNewPts > 0 )
-    {
-    T TValue = static_cast<T>(value);
-    std::fill_n(this->NewScalars+eMD[0][0], numNewPts, TValue);
-    }
-
   // Traverse all voxels in this row, those containing the contour are
   // further identified for processing, meaning generating points and
   // triangles. Begin by setting up point ids on voxel edges.
@@ -1121,7 +1116,10 @@ Contour(vtkFlyingEdges3D *self, vtkImageData *input, int extent[6],
     // independent threads can write without collisions. Once allocation is
     // complete, the volume is processed on a voxel row by row basis to
     // produce output points and triangles, and interpolate point attribute
-    // data (as necessary).
+    // data (as necessary). NOTE: This implementation is serial. It is
+    // possible to use a threaded prefix sum to make it even faster. Since
+    // this pass usually takes a small amount of time, we choose simplicity
+    // over performance.
     numOutXPts = startXPts;
     numOutYPts = startYPts;
     numOutZPts = startZPts;
@@ -1158,6 +1156,8 @@ Contour(vtkFlyingEdges3D *self, vtkImageData *input, int extent[6],
       {
       newScalars->WriteVoidPointer(0,(numOutXPts+numOutYPts+numOutZPts));
       algo.NewScalars = static_cast<T*>(newScalars->GetVoidPointer(0));
+      T TValue = static_cast<T>(value);
+      std::fill_n(algo.NewScalars, numOutXPts+numOutYPts+numOutZPts, TValue);
       }
     if (newGradients)
       {
@@ -1171,9 +1171,12 @@ Contour(vtkFlyingEdges3D *self, vtkImageData *input, int extent[6],
       }
     algo.NeedGradients = (algo.NewGradients || algo.NewNormals ? 1 : 0);
 
-    // Process voxel rows and generate output
-    Pass3<T> pass3(&algo,value);
-    vtkSMPTools::For(0,algo.Dims[2]-1, pass3);
+    // PASS 4: Fourth and final pass: Process voxel rows and generate output.
+    // Note that we are simultaneously generating triangles and interpolating
+    // points. These could be split into separate, parallel operations for
+    // maximum performance.
+    Pass4<T> pass4(&algo,value);
+    vtkSMPTools::For(0,algo.Dims[2]-1, pass4);
 
     // Handle multiple contours
     startXPts = numOutXPts;
@@ -1320,7 +1323,7 @@ int vtkFlyingEdges3D::RequestData(
     {
     newNormals = vtkFloatArray::New();
     newNormals->SetNumberOfComponents(3);
-    newNormals->SetName("Gradients");
+    newNormals->SetName("Normals");
     }
   if (this->ComputeGradients)
     {
@@ -1367,7 +1370,8 @@ int vtkFlyingEdges3D::RequestData(
 
   if (newGradients)
     {
-    output->GetPointData()->AddArray(newGradients);
+    int idx = output->GetPointData()->AddArray(newGradients);
+    output->GetPointData()->SetActiveAttribute(idx, vtkDataSetAttributes::VECTORS);
     newGradients->Delete();
     }
 

Filters/Core/vtkFlyingEdges3D.h

@@ -25,22 +25,30 @@
 // which partial computational results can be used to eliminate future
 // computations.
 //
-// This is a three-pass algorithm. The first pass processes all x-edges and
+// This is a four-pass algorithm. The first pass processes all x-edges and
 // builds x-edge case values (which, when the four x-edges defining a voxel
 // are combined, are equivalent to vertex-based case table except edge-based
 // approaches are separable in support of parallel computing). Next x-voxel
 // rows are processed to gather information from yz-edges (basically to count
-// the number of edge intersections and triangles generated). Finally in the
-// third pass output primitives are generated into pre-allocated arrays. This
-// implementation uses voxel cell axes (a x-y-z triad located at the voxel
-// origin) to ensure that each edge is intersected at most one time. Note
-// that this implementation also reuses the VTK Marching Cubes case table,
-// although the vertex-based MC table is transformed into an edge-based
-// table.
+// the number of y-z edge intersections and triangles generated). In the third
+// pass a prefix sum is used to count and allocate memory for the output
+// primitives. Finally in the fourth pass output primitives are generated into
+// pre-allocated arrays. This implementation uses voxel cell axes (a x-y-z
+// triad located at the voxel origin) to ensure that each edge is intersected
+// at most one time. Note that this implementation also reuses the VTK
+// Marching Cubes case table, although the vertex-based MC table is
+// transformed into an edge-based table on object instantiation.
+//
+// See the paper "Flying Edges: A High-Performance Scalable Isocontouring
+// Algorithm" by Schroeder, Maynard, Geveci. Proc. of LDAV 2015. Chicago, IL.
 
 // .SECTION Caveats
 // This filter is specialized to 3D volumes. This implementation can produce
 // degenerate triangles (i.e., zero-area triangles).
+//
+// This class has been threaded with vtkSMPTools. Using TBB or other
+// non-sequential type (set in the CMake variable
+// VTK_SMP_IMPLEMENTATION_TYPE) may improve performance significantly.
 
 // .SECTION See Also
 // vtkContourFilter vtkFlyingEdges2D vtkSynchronizedTemplates3D