Merge topic 'Point-Density-Gradient'

5bbd84e Added capability to compute gradient

Acked-by: Kitware Robot <[email protected]>
Reviewed-by: Andrew Maclean <[email protected]>
Merge-request: !2096

Filters/Points/Testing/Data/Baseline/TestPointDensityGradient.png.md5

@@ -0,0 +1 @@
+c6d0ac5938cc2feb88eed6d499334316

Filters/Points/Testing/Python/CMakeLists.txt

@@ -22,6 +22,7 @@ if ("${VTK_RENDERING_BACKEND}" STREQUAL "OpenGL2")
     TestPCANormalEstimation.py
     TestPCANormalEstimation2.py
     TestPointDensityFilter.py
+    TestPointDensityGradient.py
     TestRadiusOutlierRemoval.py
     TestSignedDistanceFilter.py
     TestStatisticalOutlierRemoval.py

Filters/Points/Testing/Python/TestPointDensityGradient.py

@@ -0,0 +1,159 @@
+#!/usr/bin/env python
+import vtk
+from vtk.test import Testing
+from vtk.util.misc import vtkGetDataRoot
+VTK_DATA_ROOT = vtkGetDataRoot()
+
+# The resolution of the density function volume
+res = 100
+
+# Parameters for debugging
+NPts = 1000000
+math = vtk.vtkMath()
+math.RandomSeed(31415)
+
+# create pipeline
+#
+points = vtk.vtkBoundedPointSource()
+points.SetNumberOfPoints(NPts)
+points.ProduceRandomScalarsOn()
+points.ProduceCellOutputOff()
+points.Update()
+
+# Create a sphere implicit function
+sphere = vtk.vtkSphere()
+sphere.SetCenter(0.0,0.1,0.2)
+sphere.SetRadius(0.75)
+
+# Extract points within sphere
+extract = vtk.vtkFitImplicitFunction()
+extract.SetInputConnection(points.GetOutputPort())
+extract.SetImplicitFunction(sphere)
+extract.SetThreshold(0.005)
+extract.GenerateVerticesOn()
+
+# Clip out some of the points with a plane; requires vertices
+plane = vtk.vtkPlane()
+plane.SetOrigin(sphere.GetCenter())
+plane.SetNormal(1,1,1)
+
+clipper = vtk.vtkClipPolyData()
+clipper.SetInputConnection(extract.GetOutputPort())
+clipper.SetClipFunction(plane);
+
+# Generate density field from points
+# Use fixed radius
+dens0 = vtk.vtkPointDensityFilter()
+dens0.SetInputConnection(clipper.GetOutputPort())
+dens0.SetSampleDimensions(res,res,res)
+dens0.SetDensityEstimateToFixedRadius()
+dens0.SetRadius(0.05)
+#dens0.SetDensityEstimateToRelativeRadius()
+dens0.SetRelativeRadius(2.5)
+#dens0.SetDensityFormToVolumeNormalized()
+dens0.SetDensityFormToNumberOfPoints()
+dens0.ComputeGradientOn()
+dens0.Update()
+vrange = dens0.GetOutput().GetPointData().GetArray("Gradient Magnitude").GetRange()
+
+# Show the gradient magnitude
+assign0 = vtk.vtkAssignAttribute()
+assign0.SetInputConnection(dens0.GetOutputPort())
+assign0.Assign("Gradient Magnitude", "SCALARS", "POINT_DATA")
+
+map0 = vtk.vtkImageSliceMapper()
+map0.BorderOn()
+map0.SliceAtFocalPointOn()
+map0.SliceFacesCameraOn()
+map0.SetInputConnection(assign0.GetOutputPort())
+
+slice0 = vtk.vtkImageSlice()
+slice0.SetMapper(map0)
+slice0.GetProperty().SetColorWindow(vrange[1]-vrange[0])
+slice0.GetProperty().SetColorLevel(0.5*(vrange[0]+vrange[1]))
+
+# Show the region labels
+assign1 = vtk.vtkAssignAttribute()
+assign1.SetInputConnection(dens0.GetOutputPort())
+assign1.Assign("Classification", "SCALARS", "POINT_DATA")
+
+map1 = vtk.vtkImageSliceMapper()
+map1.BorderOn()
+map1.SliceAtFocalPointOn()
+map1.SliceFacesCameraOn()
+map1.SetInputConnection(assign1.GetOutputPort())
+map1.Update()
+
+slice1 = vtk.vtkImageSlice()
+slice1.SetMapper(map1)
+slice1.GetProperty().SetColorWindow(1)
+slice1.GetProperty().SetColorLevel(0.5)
+
+# Show the vectors (gradient)
+assign2 = vtk.vtkAssignAttribute()
+assign2.SetInputConnection(dens0.GetOutputPort())
+assign2.Assign("Gradient", "VECTORS", "POINT_DATA")
+
+plane = vtk.vtkPlane()
+plane.SetNormal(0,0,1)
+plane.SetOrigin(0.0701652, 0.172689, 0.27271)
+
+cut = vtk.vtkFlyingEdgesPlaneCutter()
+cut.SetInputConnection(assign2.GetOutputPort())
+cut.SetPlane(plane)
+cut.InterpolateAttributesOn()
+
+v = vtk.vtkHedgeHog()
+v.SetInputConnection(cut.GetOutputPort())
+v.SetScaleFactor(0.0001)
+
+vMapper = vtk.vtkPolyDataMapper()
+vMapper.SetInputConnection(v.GetOutputPort())
+vMapper.SetScalarRange(vrange)
+
+vectors = vtk.vtkActor()
+vectors.SetMapper(vMapper)
+
+# Create the RenderWindow, Renderer and both Actors
+#
+ren0 = vtk.vtkRenderer()
+ren0.SetViewport(0,0,0.333,1)
+ren1 = vtk.vtkRenderer()
+ren1.SetViewport(0.333,0,0.667,1)
+ren2 = vtk.vtkRenderer()
+ren2.SetViewport(0.667,0,1,1)
+
+renWin = vtk.vtkRenderWindow()
+renWin.AddRenderer(ren0)
+renWin.AddRenderer(ren1)
+renWin.AddRenderer(ren2)
+iren = vtk.vtkRenderWindowInteractor()
+iren.SetRenderWindow(renWin)
+
+# Add the actors to the renderer, set the background and size
+#
+ren0.AddActor(slice0)
+ren0.SetBackground(0,0,0)
+ren1.AddActor(slice1)
+ren1.SetBackground(0,0,0)
+ren2.AddActor(vectors)
+ren2.SetBackground(0,0,0)
+
+renWin.SetSize(450,150)
+
+cam = ren0.GetActiveCamera()
+cam.ParallelProjectionOn()
+cam.SetFocalPoint(0,0,0)
+cam.SetPosition(0,0,1)
+ren0.ResetCamera()
+
+ren1.SetActiveCamera(cam)
+ren2.SetActiveCamera(cam)
+
+iren.Initialize()
+
+# render the image
+#
+renWin.Render()
+
+#iren.Start()

Filters/Points/vtkPointDensityFilter.cxx

@@ -15,6 +15,7 @@
 #include "vtkPointDensityFilter.h"
 
 #include "vtkObjectFactory.h"
+#include "vtkUnsignedCharArray.h"
 #include "vtkFloatArray.h"
 #include "vtkDoubleArray.h"
 #include "vtkPointSet.h"
@@ -201,6 +202,113 @@ struct ComputeWeightedDensity : public ComputePointDensity
   }
 }; //ComputeWeightedDensity
 
+//----------------------------------------------------------------------------
+// Optional kernel to compute gradient of density function. Also the gradient
+// magnitude and function classification is computed.
+struct ComputeGradients
+{
+  int Dims[3];
+  double Origin[3];
+  double Spacing[3];
+  float *Density;
+  float *Gradients;
+  float *GradientMag;
+  unsigned char *FuncClassification;
+
+  ComputeGradients(int dims[3], double origin[3], double spacing[3], float *dens,
+                   float *grad, float *mag, unsigned char *fclass) :
+    Density(dens), Gradients(grad), GradientMag(mag), FuncClassification(fclass)
+
+  {
+    for (int i=0; i < 3; ++i)
+    {
+      this->Dims[i] = dims[i];
+      this->Origin[i] = origin[i];
+      this->Spacing[i] = spacing[i];
+    }
+  }
+
+  void operator() (vtkIdType slice, vtkIdType sliceEnd)
+  {
+    double *spacing=this->Spacing;
+    int *dims=this->Dims;
+    vtkIdType sliceSize=dims[0]*dims[1];
+    float *d = this->Density + slice*sliceSize;
+    float *grad = this->Gradients + 3*slice*sliceSize;
+    float *mag = this->GradientMag + slice*sliceSize;
+    unsigned char *fclass = this->FuncClassification + slice*sliceSize;
+    float f, dp, dm;
+    bool nonZeroComp;
+    int idx[3], incs[3];
+    incs[0] = 1;
+    incs[1] = dims[0];
+    incs[2] = sliceSize;
+
+    for ( ; slice < sliceEnd; ++slice )
+    {
+      idx[2] = slice;
+      for ( int j=0;  j < dims[1]; ++j)
+      {
+        idx[1] = j;
+        for ( int i=0; i < dims[0]; ++i)
+        {
+          idx[0] = i;
+          nonZeroComp = (d[0] != 0.0 || d[1] != 0.0 || d[2] != 0.0 ? true : false);
+          for (int ii=0; ii < 3; ++ii)
+          {
+            if ( idx[ii] == 0 )
+            {
+              dm = *d;
+              dp = *(d + incs[ii]);
+              f = 1.0;
+            }
+            else if ( idx[ii] == dims[ii]-1 )
+            {
+              dm = *(d - incs[ii]);
+              dp = *d;
+              f = 1.0;
+            }
+            else
+            {
+              dm = *(d - incs[ii]);
+              dp = *(d + incs[ii]);
+              f = 0.5;
+            }
+            grad[ii] = f * (dp-dm) / spacing[ii];
+            nonZeroComp = ( dp != 0.0 || dm != 0.0 ? true : nonZeroComp );
+          }
+
+          // magnitude
+          if ( nonZeroComp )
+          {
+            *mag++ = vtkMath::Norm(grad);
+            *fclass++ = vtkPointDensityFilter::NON_ZERO;
+          }
+          else
+          {
+            *mag++ = 0.0;
+            *fclass++ = vtkPointDensityFilter::ZERO;
+          }
+          d++;
+          grad += 3;
+
+        }//over i
+      }//over j
+    }//over slices
+  }
+
+  void Reduce()
+  {
+  }
+
+  static void Execute(int dims[3], double origin[3], double spacing[3],
+                      float *density, float *grad, float *mag,
+                      unsigned char *fclass)
+  {
+    ComputeGradients compGrad(dims, origin, spacing, density, grad, mag, fclass);
+    vtkSMPTools::For(0, dims[2], compGrad);
+  }
+}; //ComputeGradients
 
 } //anonymous namespace
 
@@ -233,6 +341,8 @@ vtkPointDensityFilter::vtkPointDensityFilter()
 
   this->ScalarWeighting = false;
 
+  this->ComputeGradient = false;
+
   this->Locator = vtkStaticPointLocator::New();
 
 }
@@ -500,6 +610,43 @@ int vtkPointDensityFilter::RequestData(
     }
   }
 
+  // If the gradient is requested, compute the vector gradient and magnitude.
+  // Also compute the classification of the gradient value.
+  if ( this->ComputeGradient )
+  {
+    //Allocate space
+    vtkIdType num = density->GetNumberOfTuples();
+
+    vtkFloatArray *gradients = vtkFloatArray::New();
+    gradients->SetNumberOfComponents(3);
+    gradients->SetNumberOfTuples(num);
+    gradients->SetName("Gradient");
+    output->GetPointData()->AddArray(gradients);
+    float *grad = static_cast<float*>(gradients->GetVoidPointer(0));
+    gradients->Delete();
+
+    vtkFloatArray *magnitude = vtkFloatArray::New();
+    magnitude->SetNumberOfComponents(1);
+    magnitude->SetNumberOfTuples(num);
+    magnitude->SetName("Gradient Magnitude");
+    output->GetPointData()->AddArray(magnitude);
+    float *mag = static_cast<float*>(magnitude->GetVoidPointer(0));
+    magnitude->Delete();
+
+    vtkUnsignedCharArray *fclassification =
+      vtkUnsignedCharArray::New();
+    fclassification->SetNumberOfComponents(1);
+    fclassification->SetNumberOfTuples(num);
+    fclassification->SetName("Classification");
+    output->GetPointData()->AddArray(fclassification);
+    unsigned char *fclass =
+      static_cast<unsigned char*>(fclassification->GetVoidPointer(0));
+    fclassification->Delete();
+
+    //Thread the computation over slices
+    ComputeGradients::Execute(dims,origin,spacing,d,grad,mag,fclass);
+  }
+
   return 1;
 }
 
@@ -561,5 +708,8 @@ void vtkPointDensityFilter::PrintSelf(ostream& os, vtkIndent indent)
   os << indent << "Scalar Weighting: "
      << (this->ScalarWeighting ? "On\n" : "Off\n");
 
+  os << indent << "Compute Gradient: "
+     << (this->ComputeGradient ? "On\n" : "Off\n");
+
   os << indent << "Locator: " << this->Locator << "\n";
 }