Threaded and templated elevation filters

Used vtkSMPTools and templates to speed up the elevation filters.
A fast path was created for processing vtkPointSet input data.

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

+3e99ef9b2f4b4438084874531b97b5f5

Filters/Core/Testing/Python/CMakeLists.txt

@@ -15,6 +15,7 @@ vtk_add_test_python(
   Delaunay3DAlphaTest.py
   QuadricDecimation.py
   StreamPolyData.py
+  TestElevationFilter.py
   TestGridSynchronizedTemplates3D.py
   TestMarchingSquares.py
   TestRectilinearSynchronizedTemplates.py

Filters/Core/Testing/Python/TestElevationFilter.py

+#!/usr/bin/env python
+
+import math
+import vtk
+import vtk.test.Testing
+from vtk.util.misc import vtkGetDataRoot
+VTK_DATA_ROOT = vtkGetDataRoot()
+
+res = 200
+
+# Create the RenderWindow, Renderer and both Actors
+#
+ren1 = vtk.vtkRenderer()
+ren2 = vtk.vtkRenderer()
+renWin = vtk.vtkRenderWindow()
+renWin.AddRenderer(ren1)
+renWin.AddRenderer(ren2)
+
+# Create the data -- a plane
+#
+plane = vtk.vtkPlaneSource()
+plane.SetXResolution(res)
+plane.SetYResolution(res)
+plane.SetOrigin(-10,-10,0)
+plane.SetPoint1(10,-10,0)
+plane.SetPoint2(-10,10,0)
+
+tf = vtk.vtkTriangleFilter()
+tf.SetInputConnection(plane.GetOutputPort())
+tf.Update()
+
+pd = vtk.vtkPolyData()
+pd.CopyStructure(tf.GetOutput())
+numPts = pd.GetNumberOfPoints()
+oldPts = tf.GetOutput().GetPoints()
+newPts = vtk.vtkPoints()
+newPts.SetNumberOfPoints(numPts)
+for i in range(0, numPts):
+    pt = oldPts.GetPoint(i)
+    r = math.sqrt(pt[0]*pt[0] + pt[1]*pt[1])
+    z = 1.5*math.cos(2*r)
+    newPts.SetPoint(i, pt[0], pt[1], z)
+pd.SetPoints(newPts)
+
+ele = vtk.vtkSimpleElevationFilter()
+ele.SetInputData(pd)
+
+mapper = vtk.vtkPolyDataMapper()
+mapper.SetInputConnection(ele.GetOutputPort())
+
+actor = vtk.vtkActor()
+actor.SetMapper(mapper)
+
+ele2 = vtk.vtkElevationFilter()
+ele2.SetInputData(pd)
+ele2.SetLowPoint(0,0,-1.5)
+ele2.SetHighPoint(0,0,1.5)
+ele2.SetScalarRange(-1.5,1.5)
+
+mapper2 = vtk.vtkPolyDataMapper()
+mapper2.SetInputConnection(ele2.GetOutputPort())
+
+actor2 = vtk.vtkActor()
+actor2.SetMapper(mapper2)
+
+# Add the actors to the renderer, set the background and size
+#
+ren1.SetViewport(0, 0, .5, 1)
+ren2.SetViewport(.5, 0, 1, 1)
+
+ren1.AddActor(actor)
+ren2.AddActor(actor2)
+
+camera = vtk.vtkCamera()
+camera.SetPosition(1,1,1)
+ren1.SetActiveCamera(camera)
+ren2.SetActiveCamera(camera)
+
+ren1.SetBackground(0, 0, 0)
+ren2.SetBackground(0, 0, 0)
+
+renWin.SetSize(500, 250)
+
+# render and interact with data
+
+iRen = vtk.vtkRenderWindowInteractor()
+iRen.SetRenderWindow(renWin);
+ren1.ResetCamera()
+renWin.Render()
+
+iRen.Initialize()
+#iRen.Start()

Filters/Core/vtkElevationFilter.cxx

@@ -16,6 +16,7 @@
 
 #include "vtkCellData.h"
 #include "vtkDataSet.h"
+#include "vtkPointSet.h"
 #include "vtkFloatArray.h"
 #include "vtkInformation.h"
 #include "vtkInformationVector.h"
@@ -23,10 +24,103 @@
 #include "vtkObjectFactory.h"
 #include "vtkPointData.h"
 #include "vtkSmartPointer.h"
+#include "vtkSMPTools.h"
 
 vtkStandardNewMacro(vtkElevationFilter);
 
+// The heart of the algorithm plus interface to the SMP tools. Double templated
+// over point and scalar types.
+template <class TP>
+class vtkElevationAlgorithm
+{
+public:
+  vtkIdType NumPts;
+  double LowPoint[3];
+  double HighPoint[3];
+  double ScalarRange[2];
+  const TP *Points;
+  float *Scalars;
+  const double *V;
+  double L2;
+
+  // Contructor
+  vtkElevationAlgorithm();
+
+  // Interface between VTK and templated functions
+  static void Elevate(vtkElevationFilter *self, vtkIdType numPts,
+                      double v[3], double l2, TP *points, float *scalars);
+
+  // Interface implicit function computation to SMP tools.
+  template <class TP> class ElevationOp
+    {
+    public:
+      ElevationOp(vtkElevationAlgorithm<TP> *algo)
+        { this->Algo = algo;}
+      vtkElevationAlgorithm *Algo;
+      void  operator() (vtkIdType k, vtkIdType end)
+        {
+        double ns, vec[3];
+        const double *range = this->Algo->ScalarRange;
+        const double diffScalar = range[1] - range[0];
+        const double *v = this->Algo->V;
+        const double l2 = this->Algo->L2;
+        const double *lp = this->Algo->LowPoint;
+        const TP *p = this->Algo->Points + 3*k;
+        float *s = this->Algo->Scalars + k;
+        for ( ; k < end; ++k)
+          {
+          vec[0] = p[0] - lp[0];
+          vec[1] = p[1] - lp[1];
+          vec[2] = p[2] - lp[2];
+          ns = (vec[0]*v[0] + vec[1]*v[1] + vec[2]*v[2]) / l2;
+          ns = (ns < 0.0 ? 0.0 : ns > 1.0 ? 1.0 : ns);
+
+          // Store the resulting scalar value.
+          *s = range[0] + ns*diffScalar;
+
+          p+=3;
+          ++s;
+          }
+        }
+    };
+};
+
+//----------------------------------------------------------------------------
+// Initialized mainly to eliminate compiler warnings.
+template <class TP> vtkElevationAlgorithm<TP>::
+vtkElevationAlgorithm():Points(NULL),Scalars(NULL)
+{
+  this->LowPoint[0] = this->LowPoint[1] = this->LowPoint[2] = 0.0;
+  this->HighPoint[0] = this->HighPoint[1] = 0.0;
+  this->HighPoint[2] = 1.0;
+  this->ScalarRange[0] = 0.0;
+  this->ScalarRange[1] = 1.0;
+}
+
 //----------------------------------------------------------------------------
+// Templated class is glue between VTK and templated algorithms.
+template <class TP> void vtkElevationAlgorithm<TP>::
+Elevate(vtkElevationFilter *self, vtkIdType numPts,
+        double *v, double l2, TP *points, float *scalars)
+{
+  // Populate data into local storage
+  vtkElevationAlgorithm<TP> algo;
+  algo.NumPts = numPts;
+  self->GetLowPoint(algo.LowPoint);
+  self->GetHighPoint(algo.HighPoint);
+  self->GetScalarRange(algo.ScalarRange);
+  algo.Points = points;
+  algo.Scalars = scalars;
+  algo.V = v;
+  algo.L2 = l2;
+
+  // Okay now generate samples using SMP tools
+  ElevationOp<TP> values(&algo);
+  vtkSMPTools::For(0,algo.NumPts, values);
+}
+
+//----------------------------------------------------------------------------
+// Begin the class proper
 vtkElevationFilter::vtkElevationFilter()
 {
   this->LowPoint[0] = 0.0;
@@ -100,34 +194,57 @@ int vtkElevationFilter::RequestData(vtkInformation*,
     length2 = 1.0;
     }
 
-  // Support progress and abort.
-  vtkIdType tenth = (numPts >= 10? numPts/10 : 1);
-  double numPtsInv = 1.0/numPts;
-  int abort = 0;
-
-  // Compute parametric coordinate and map into scalar range.
-  double diffScalar = this->ScalarRange[1] - this->ScalarRange[0];
   vtkDebugMacro("Generating elevation scalars!");
-  for(vtkIdType i=0; i < numPts && !abort; ++i)
+
+  // Create a fast path for point set input
+  //
+  vtkPointSet *ps = vtkPointSet::SafeDownCast(input);
+//  if ( ps )
+  if ( 0 )
     {
-    // Periodically update progress and check for an abort request.
-    if(i % tenth == 0)
+    float *scalars =
+      static_cast<float*>(newScalars->GetVoidPointer(0));
+    vtkPoints *points = ps->GetPoints();
+    void *pts = points->GetData()->GetVoidPointer(0);
+    switch ( points->GetDataType() )
       {
-      this->UpdateProgress((i+1)*numPtsInv);
-      abort = this->GetAbortExecute();
+      vtkTemplateMacro(
+        vtkElevationAlgorithm<VTK_TT>::Elevate(this,numPts,diffVector,length2,
+                                               (VTK_TT *)pts,scalars));
       }
+    }//fast path
+
+  else
+    {
+    // Too bad, got to take the scenic route.
+    // Support progress and abort.
+    vtkIdType tenth = (numPts >= 10? numPts/10 : 1);
+    double numPtsInv = 1.0/numPts;
+    int abort = 0;
 
-    // Project this input point into the 1D system.
-    double x[3];
-    input->GetPoint(i, x);
-    double v[3] = { x[0] - this->LowPoint[0],
-                    x[1] - this->LowPoint[1],
-                    x[2] - this->LowPoint[2] };
-    double s = vtkMath::Dot(v, diffVector) / length2;
-    s = (s < 0.0 ? 0.0 : s > 1.0 ? 1.0 : s);
-
-    // Store the resulting scalar value.
-    newScalars->SetValue(i, this->ScalarRange[0] + s*diffScalar);
+    // Compute parametric coordinate and map into scalar range.
+    double diffScalar = this->ScalarRange[1] - this->ScalarRange[0];
+    for(vtkIdType i=0; i < numPts && !abort; ++i)
+      {
+      // Periodically update progress and check for an abort request.
+      if(i % tenth == 0)
+        {
+        this->UpdateProgress((i+1)*numPtsInv);
+        abort = this->GetAbortExecute();
+        }
+
+      // Project this input point into the 1D system.
+      double x[3];
+      input->GetPoint(i, x);
+      double v[3] = { x[0] - this->LowPoint[0],
+                      x[1] - this->LowPoint[1],
+                      x[2] - this->LowPoint[2] };
+      double s = vtkMath::Dot(v, diffVector) / length2;
+      s = (s < 0.0 ? 0.0 : s > 1.0 ? 1.0 : s);
+
+      // Store the resulting scalar value.
+      newScalars->SetValue(i, this->ScalarRange[0] + s*diffScalar);
+      }
     }
 
   // Copy all the input geometry and data to the output.

Filters/Core/vtkElevationFilter.h

@@ -20,6 +20,17 @@
 // a line. The line can be oriented arbitrarily. A typical example is
 // to generate scalars based on elevation or height above a plane.
 
+// .SECTION Caveats
+// vtkSimpleElevationFilter may be easier to use in many cases; e.g.,
+// compute vertical elevation above zero z-point.
+//
+// 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
+// vtkSimpleElevationFilter
+
 #ifndef vtkElevationFilter_h
 #define vtkElevationFilter_h
 
@@ -61,6 +72,7 @@ protected:
   double LowPoint[3];
   double HighPoint[3];
   double ScalarRange[2];
+
 private:
   vtkElevationFilter(const vtkElevationFilter&);  // Not implemented.
   void operator=(const vtkElevationFilter&);  // Not implemented.

Filters/Core/vtkSimpleElevationFilter.cxx

@@ -15,18 +15,87 @@
 #include "vtkSimpleElevationFilter.h"
 
 #include "vtkCellData.h"
-#include "vtkDataSet.h"
+#include "vtkPointSet.h"
 #include "vtkFloatArray.h"
 #include "vtkInformation.h"
 #include "vtkInformationVector.h"
 #include "vtkMath.h"
 #include "vtkObjectFactory.h"
 #include "vtkPointData.h"
+#include "vtkSMPTools.h"
 
 vtkStandardNewMacro(vtkSimpleElevationFilter);
 
-// Construct object with LowPoint=(0,0,0) and HighPoint=(0,0,1). Scalar
-// range is (0,1).
+// The heart of the algorithm plus interface to the SMP tools. Double templated
+// over point and scalar types.
+template <class TP>
+class vtkSimpleElevationAlgorithm
+{
+public:
+  vtkIdType NumPts;
+  double Vector[3];
+  const TP *Points;
+  float *Scalars;
+
+  // Contructor
+  vtkSimpleElevationAlgorithm();
+
+  // Interface between VTK and templated functions
+  static void Elevate(vtkSimpleElevationFilter *self, vtkIdType numPts,
+                      const TP *points, float *scalars);
+
+  // Interface implicit function computation to SMP tools.
+  template <class TP> class ElevationOp
+    {
+    public:
+      ElevationOp(vtkSimpleElevationAlgorithm<TP> *algo)
+        { this->Algo = algo;}
+      vtkSimpleElevationAlgorithm *Algo;
+      void  operator() (vtkIdType k, vtkIdType end)
+        {
+        const double *v = this->Algo->Vector;
+        const TP *p = this->Algo->Points + 3*k;
+        float *s = this->Algo->Scalars + k;
+        for ( ; k < end; ++k)
+          {
+          *s = v[0]*p[0] + v[1]*p[1] + v[2]*p[2];
+          p+=3;
+          ++s;
+          }
+        }
+    };
+};
+
+//----------------------------------------------------------------------------
+// Initialized mainly to eliminate compiler warnings.
+template <class TP> vtkSimpleElevationAlgorithm<TP>::
+vtkSimpleElevationAlgorithm():Points(NULL),Scalars(NULL)
+{
+  this->Vector[0] = this->Vector[1] = this->Vector[2] = 0.0;
+}
+
+//----------------------------------------------------------------------------
+// Templated class is glue between VTK and templated algorithms.
+template <class TP> void vtkSimpleElevationAlgorithm<TP>::
+Elevate(vtkSimpleElevationFilter *self, vtkIdType numPts,
+        const TP *points, float *scalars)
+{
+  // Populate data into local storage
+  vtkSimpleElevationAlgorithm<TP> algo;
+  algo.NumPts = numPts;
+  self->GetVector(algo.Vector);
+  algo.Points = points;
+  algo.Scalars = scalars;
+
+  // Okay now generate samples using SMP tools
+  ElevationOp<TP> values(&algo);
+  vtkSMPTools::For(0,algo.NumPts, values);
+}
+
+
+// Okay begin the class proper
+//----------------------------------------------------------------------------
+// Construct object with Vector=(0,0,1).
 vtkSimpleElevationFilter::vtkSimpleElevationFilter()
 {
   this->Vector[0] = 0.0;
@@ -34,7 +103,8 @@ vtkSimpleElevationFilter::vtkSimpleElevationFilter()
   this->Vector[2] = 1.0;
 }
 
-// Convert position along ray into scalar value.  Example use includes
+//----------------------------------------------------------------------------
+// Convert position along the ray into scalar value.  Example use includes
 // coloring terrain by elevation.
 //
 int vtkSimpleElevationFilter::RequestData(
@@ -60,7 +130,7 @@ int vtkSimpleElevationFilter::RequestData(
   //
   vtkDebugMacro(<<"Generating elevation scalars!");
 
- // First, copy the input to the output as a starting point
+  // First, copy the input to the output as a starting point
   output->CopyStructure( input );
 
   if ( ((numPts=input->GetNumberOfPoints()) < 1) )
@@ -74,7 +144,7 @@ int vtkSimpleElevationFilter::RequestData(
   newScalars = vtkFloatArray::New();
   newScalars->SetNumberOfTuples(numPts);
 
-  // Set up 1D parametric system
+  // Ensure that there is a valid vector
   //
   if ( vtkMath::Dot(this->Vector,this->Vector) == 0.0)
     {
@@ -82,21 +152,42 @@ int vtkSimpleElevationFilter::RequestData(
     this->Vector[0] = this->Vector[1] = 0.0; this->Vector[2] = 1.0;
     }
 
-  // Compute dot product
+  // Create a fast path for point set input
   //
-  int abort=0;
-  vtkIdType progressInterval=numPts/20 + 1;
-  for (i=0; i<numPts && !abort; i++)
+  vtkPointSet *ps = vtkPointSet::SafeDownCast(input);
+//  if ( ps )
+  if ( 0 )
     {
-    if ( ! (i % progressInterval) )
+    float *scalars =
+      static_cast<float*>(newScalars->GetVoidPointer(0));
+    vtkPoints *points = ps->GetPoints();
+    void *pts = points->GetData()->GetVoidPointer(0);
+    switch ( points->GetDataType() )
       {
-      this->UpdateProgress ((double)i/numPts);
-      abort = this->GetAbortExecute();
+      vtkTemplateMacro(vtkSimpleElevationAlgorithm<VTK_TT>::
+                       Elevate(this,numPts,(VTK_TT *)pts,scalars));
       }
+    }//fast path
 
-    input->GetPoint(i,x);
-    s = vtkMath::Dot(this->Vector,x);
-    newScalars->SetComponent(i,0,s);
+  else
+    {
+    // Too bad, got to take the scenic route.
+    // Compute dot product.
+    //
+    int abort=0;
+    vtkIdType progressInterval=numPts/20 + 1;
+    for (i=0; i<numPts && !abort; i++)
+      {
+      if ( ! (i % progressInterval) )
+        {
+        this->UpdateProgress ((double)i/numPts);
+        abort = this->GetAbortExecute();
+        }
+
+      input->GetPoint(i,x);
+      s = vtkMath::Dot(this->Vector,x);
+      newScalars->SetComponent(i,0,s);
+      }
     }
 
   // Update self
@@ -114,6 +205,7 @@ int vtkSimpleElevationFilter::RequestData(
   return 1;
 }
 
+//----------------------------------------------------------------------------
 void vtkSimpleElevationFilter::PrintSelf(ostream& os, vtkIndent indent)
 {
   this->Superclass::PrintSelf(os,indent);

Filters/Core/vtkSimpleElevationFilter.h

@@ -18,6 +18,18 @@
 // dataset.  The scalar values are generated by dotting a user-specified
 // vector against a vector defined from the input dataset points to the
 // origin.
+
+// .SECTION Caveats
+// 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.  Note
+// however that the associated threading/templating implements a fast path,
+// which will only kick in when the input has explicit point representation
+// (e.g. a vtkPointSet).
+//
+// See also vtkElevationFilter provides more control over the operation,
+// including clamping the output scalars within a range.
+
 // .SECTION See Also
 // vtkElevationFilter
 
@@ -34,11 +46,11 @@ public:
   void PrintSelf(ostream& os, vtkIndent indent);
 
   // Description:
-  // Construct object with Vector=(0,0,1);
+  // Construct object with Vector=(0,0,1).
   static vtkSimpleElevationFilter *New();
 
   // Description:
-  // Define one end of the line (small scalar values).
+  // Define the vector with which to dot against.
   vtkSetVector3Macro(Vector,double);
   vtkGetVectorMacro(Vector,double,3);
 
@@ -54,5 +66,3 @@ private:
 };
 
 #endif
-
-