Modify vtkDecimatePolylineFilter to accommodate multiple polylines.

This fix is in reference to bug report 0015055. vtkDecimatePolylineFilter can
now be run on polydata containing multiple distinct polylines, error values
for each point are updated during decimation, and common points within the
same polyline are maintained as common points (rather than as duplicates).
Also, the test for this filter has been updated to reflect the augmented
feature set, and a comparative image was added.

Filters/Core/Testing/Cxx/CMakeLists.txt

@@ -10,7 +10,7 @@ vtk_add_test_cxx(${vtk-module}CxxTests tests
   TestClipPolyData.cxx,NO_VALID
   TestConnectivityFilter.cxx,NO_VALID
   TestCutter.cxx,NO_VALID
-  TestDecimatePolylineFilter.cxx,NO_VALID
+  TestDecimatePolylineFilter.cxx
   TestDecimatePro.cxx,NO_VALID
   TestDelaunay2D.cxx
   TestDelaunay3D.cxx,NO_VALID

Filters/Core/Testing/Cxx/TestDecimatePolylineFilter.cxx

@@ -29,42 +29,62 @@
 #include <vtkPolyData.h>
 #include <vtkPolyDataMapper.h>
 #include <vtkProperty.h>
+#include "vtkRegressionTestImage.h"
 #include <vtkRenderer.h>
 #include <vtkRenderWindow.h>
 #include <vtkRenderWindowInteractor.h>
 #include <vtkSmartPointer.h>
 
-int TestDecimatePolylineFilter(int vtkNotUsed(argc), char *vtkNotUsed(argv)[])
+int TestDecimatePolylineFilter(int argc, char *argv[])
 {
-  const unsigned int numberOfPoints = 100;
+  const unsigned int numberOfPointsInCircle = 100;
 
   vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
   points->SetDataType(VTK_FLOAT);
 
-  vtkIdType *lineIds = new vtkIdType[numberOfPoints+1];
+  // We will create two polylines: one complete circle, and one circular arc
+  // subtending 3/4 of a circle.
+  vtkIdType *lineIds = new vtkIdType[(numberOfPointsInCircle*7)/4+1];
 
-  for(unsigned int i = 0; i < numberOfPoints; ++i)
+  vtkIdType lineIdCounter = 0;
+
+  // First circle:
+  for(unsigned int i = 0; i < numberOfPointsInCircle; ++i)
     {
     const double angle = 2.0 * vtkMath::Pi() * static_cast<double>(i)
-      / static_cast<double>(numberOfPoints);
+      / static_cast<double>(numberOfPointsInCircle);
     points->InsertPoint(static_cast<vtkIdType>(i), std::cos(angle),
       std::sin(angle), 0.0);
-    lineIds[i] = static_cast<vtkIdType>(i);
+    lineIds[i] = lineIdCounter++;
     }
+  lineIds[numberOfPointsInCircle] = 0;
 
-  lineIds[numberOfPoints] = 0;
+  // Second circular arc:
+  for(unsigned int i = 0; i < (numberOfPointsInCircle*3)/4; ++i)
+    {
+    const double angle = 3.0 / 2.0 * vtkMath::Pi() * static_cast<double>(i)
+      / static_cast<double>((numberOfPointsInCircle*3)/4);
+    points->InsertPoint(static_cast<vtkIdType>(i+numberOfPointsInCircle),
+                        std::cos(angle), std::sin(angle), 1.0);
+    lineIds[numberOfPointsInCircle + 1 + i] = lineIdCounter++;
+    }
 
+  // Construct associated cell array, containing both polylines.
   vtkSmartPointer<vtkCellArray> lines = vtkSmartPointer<vtkCellArray>::New();
-  lines->InsertNextCell(numberOfPoints + 1, lineIds);
+  // 1st:
+  lines->InsertNextCell(numberOfPointsInCircle + 1, lineIds);
+  // 2nd:
+  lines->InsertNextCell((numberOfPointsInCircle*3)/4,
+                        &lineIds[numberOfPointsInCircle+1]);
   delete[] lineIds;
 
-  vtkSmartPointer<vtkPolyData> circle = vtkSmartPointer<vtkPolyData>::New();
-  circle->SetPoints(points);
-  circle->SetLines(lines);
+  vtkSmartPointer<vtkPolyData> circles = vtkSmartPointer<vtkPolyData>::New();
+  circles->SetPoints(points);
+  circles->SetLines(lines);
 
   vtkSmartPointer<vtkPolyDataMapper> circleMapper
     = vtkSmartPointer<vtkPolyDataMapper>::New();
-  circleMapper->SetInputData(circle);
+  circleMapper->SetInputData(circles);
 
   vtkSmartPointer<vtkActor> circleActor = vtkSmartPointer<vtkActor>::New();
   circleActor->SetMapper(circleMapper);
@@ -72,8 +92,8 @@ int TestDecimatePolylineFilter(int vtkNotUsed(argc), char *vtkNotUsed(argv)[])
   vtkSmartPointer<vtkDecimatePolylineFilter> decimatePolylineFilter
     = vtkSmartPointer<vtkDecimatePolylineFilter>::New();
   decimatePolylineFilter->SetOutputPointsPrecision(vtkAlgorithm::DEFAULT_PRECISION);
-  decimatePolylineFilter->SetInputData(circle);
-  decimatePolylineFilter->SetTargetReduction(0.95);
+  decimatePolylineFilter->SetInputData(circles);
+  decimatePolylineFilter->SetTargetReduction(0.9);
   decimatePolylineFilter->Update();
 
   if(decimatePolylineFilter->GetOutput()->GetPoints()->GetDataType() != VTK_FLOAT)
@@ -114,13 +134,19 @@ int TestDecimatePolylineFilter(int vtkNotUsed(argc), char *vtkNotUsed(argv)[])
   vtkSmartPointer<vtkRenderWindow> renderWindow
     = vtkSmartPointer<vtkRenderWindow>::New();
   renderWindow->AddRenderer(renderer);
+  renderWindow->SetSize(300,300);
 
   vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor
     = vtkSmartPointer<vtkRenderWindowInteractor>::New();
   renderWindowInteractor->SetRenderWindow(renderWindow);
 
   renderWindow->Render();
-  renderWindowInteractor->CreateOneShotTimer(1);
 
-  return EXIT_SUCCESS;
+  int retVal = vtkRegressionTestImage(renderWindow.GetPointer());
+  if (retVal == vtkRegressionTester::DO_INTERACTOR)
+    {
+    renderWindowInteractor->Start();
+    }
+
+  return !retVal;
 }

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

+aaa862b5ab519aa449c473f8095ed2b2

Filters/Core/vtkDecimatePolylineFilter.cxx

@@ -27,46 +27,85 @@
 #include "vtkPointData.h"
 #include "vtkCellData.h"
 #include "vtkPriorityQueue.h"
-#include <vector>
-#include <queue>
+
 #include <map>
 
-struct vtkDecimatePolylineFilter::vtkDecimatePolylineVertexErrorSTLMap
+vtkStandardNewMacro(vtkDecimatePolylineFilter);
+
+// Representation of a polyline as a doubly linked list of vertices.
+class vtkDecimatePolylineFilter::Polyline
 {
-  std::map< int, double > VertexErrorMap;
-};
+public:
+  struct Vertex
+  {
+    unsigned index;
+    vtkIdType id;
+    Vertex* prev;
+    Vertex* next;
+    bool removable;
+  };
+
+  Polyline(vtkIdType* vertexOrdering, vtkIdType size)
+  {
+    this->Size = size;
+    Vertices = new Vertex[size];
+    for (unsigned idx=0;idx<size;idx++)
+      {
+      Vertices[idx].index = idx;
+      Vertices[idx].id = vertexOrdering[idx];
+      Vertices[idx].prev = (idx > 0 ? &Vertices[idx-1] : NULL);
+      Vertices[idx].next = (idx < size-1 ? &Vertices[idx+1] : NULL);
+      Vertices[idx].removable = true;
+      }
+    Vertices[0].removable = Vertices[size-1].removable = false;
+  }
 
+  ~Polyline()
+  {
+    if (Vertices)
+      {
+      delete [] Vertices;
+      Vertices = NULL;
+      }
+  }
 
-vtkStandardNewMacro(vtkDecimatePolylineFilter);
+  void Remove(vtkIdType vertexIdx)
+  {
+    this->Size--;
+    (*(this->Vertices[vertexIdx].prev)).next = this->Vertices[vertexIdx].next;
+    (*(this->Vertices[vertexIdx].next)).prev = this->Vertices[vertexIdx].prev;
+  }
+
+  vtkIdType Size;
+  Vertex* Vertices;
+};
 
 //---------------------------------------------------------------------
 // Create object with specified reduction of 90%.
 vtkDecimatePolylineFilter::vtkDecimatePolylineFilter()
 {
   this->TargetReduction = 0.90;
-  this->Closed = true;
   this->PriorityQueue = vtkSmartPointer< vtkPriorityQueue >::New();
-  this->ErrorMap = new vtkDecimatePolylineVertexErrorSTLMap;
   this->OutputPointsPrecision = vtkAlgorithm::DEFAULT_PRECISION;
 }
 
 //---------------------------------------------------------------------
 vtkDecimatePolylineFilter::~vtkDecimatePolylineFilter()
 {
-  delete this->ErrorMap;
 }
 
 //---------------------------------------------------------------------
 double
 vtkDecimatePolylineFilter::ComputeError( vtkPolyData* input,
-  int prev, int id, int next )
+                                         Polyline* polyline,
+                                         vtkIdType idx )
 {
   vtkPoints * inputPoints = input->GetPoints();
 
   double x1[3], x[3], x2[3];
-  inputPoints->GetPoint( prev, x1 );
-  inputPoints->GetPoint( id, x );
-  inputPoints->GetPoint( next, x2 );
+  inputPoints->GetPoint( polyline->Vertices[idx].prev->id, x1 );
+  inputPoints->GetPoint( polyline->Vertices[idx].id, x );
+  inputPoints->GetPoint( polyline->Vertices[idx].next->id, x2 );
 
   if ( vtkMath::Distance2BetweenPoints( x1, x2 ) == 0.0 )
     {
@@ -137,58 +176,90 @@ int vtkDecimatePolylineFilter::RequestData(
   outPD->CopyAllocate(inPD);
   outCD->CopyAllocate(inCD);
 
-  // Loop over all polylines, decimating each independently
-  vtkIdType i, cellId = 0, newId;
-  double error;
+  vtkIdType *linePtr = inputLines->GetPointer();
+  vtkIdType firstVertexIndex = 0;
+  vtkIdType polylineSize = 0;
 
-  for ( i = 0; i < numPts; ++i )
+  // Decimate each polyline (represented as a single cell) in series
+  for ( vtkIdType lineId=0;lineId<numLines;lineId++,
+          firstVertexIndex+=polylineSize)
     {
-    this->ErrorMap->VertexErrorMap[i] = 0.;
-    }
+    polylineSize = linePtr[firstVertexIndex + lineId];
 
-  for ( i = 0; i < numPts; ++i )
-    {
-    error = ComputeError( input, GetPrev(i), i, GetNext(i) );
-    this->ErrorMap->VertexErrorMap[i] = error;
-    this->PriorityQueue->Insert( error, i );
-    }//for all points in polyline
-
-  // Now process structures,
-  // deleting points until the decimation target is met.
-  vtkIdType currentNumPts = this->PriorityQueue->GetNumberOfItems();
-  while ( 1.0 - ( static_cast<double>( currentNumPts ) / static_cast<double>(
-                    numPts ) )
-          < this->TargetReduction &&
-          currentNumPts > 2 )
-    {
-    i = this->PriorityQueue->Pop( );
-    --currentNumPts;
-    UpdateError( input, i );
-    this->ErrorMap->VertexErrorMap.erase( i );
-    }
+    // construct a polyline as a doubly linked list
+    vtkDecimatePolylineFilter::Polyline* polyline =
+      new vtkDecimatePolylineFilter::Polyline(&(linePtr[firstVertexIndex +
+                                                        lineId + 1]),
+                                              polylineSize);
 
-  // What's left over is now spit out as a new polyline
-  newId = newLines->InsertNextCell( currentNumPts + 1);
-  outCD->CopyData( inCD, cellId, newId );
+    double error;
+    for (vtkIdType vertexIdx=0;vertexIdx<polyline->Size;++vertexIdx)
+      {
+      // only vertices that are removable have associated error values
+      if (polyline->Vertices[vertexIdx].removable)
+        {
+        error = this->ComputeError(input,polyline,vertexIdx);
+        this->PriorityQueue->Insert(error,vertexIdx);
+        }
+      }
 
-  std::map< int, double >::iterator it = this->ErrorMap->VertexErrorMap.begin();
+    // Now process structures,
+    // deleting vertices until the decimation target is met.
+    vtkIdType currentNumPts = polylineSize;
+    while ( 1.0 - ( static_cast<double>( currentNumPts ) / static_cast<double>(
+                      polylineSize ) )
+            < this->TargetReduction &&
+            currentNumPts > 2 )
+      {
+      --currentNumPts;
+      vtkIdType poppedIdx = this->PriorityQueue->Pop();
+      polyline->Remove(poppedIdx);
+      vtkIdType prevIdx = polyline->Vertices[poppedIdx].prev->index;
+      vtkIdType nextIdx = polyline->Vertices[poppedIdx].next->index;
+
+      // again, only vertices that are removable have associated error values
+      if (polyline->Vertices[poppedIdx].prev->removable)
+        {
+        error = this->ComputeError(input,polyline,prevIdx);
+        this->PriorityQueue->DeleteId(prevIdx);
+        this->PriorityQueue->Insert(error,prevIdx);
+        }
+
+      if (polyline->Vertices[poppedIdx].next->removable)
+        {
+        error = this->ComputeError(input,polyline,nextIdx);
+        this->PriorityQueue->DeleteId(nextIdx);
+        this->PriorityQueue->Insert(error,nextIdx);
+        }
+      }
 
-  while( it != this->ErrorMap->VertexErrorMap.end() )
-    {
-    newId = newPts->InsertNextPoint( inputPoints->GetPoint( it->first ) );
-    newLines->InsertCellPoint( newId );
-    outPD->CopyData( inPD, it->first, newId );
-    ++it;
-    }
-  if( this->Closed )
-    {
-    newId = newPts->InsertNextPoint( newPts->GetPoint( 0 ) );
-    newLines->InsertCellPoint( 0 );
-    outPD->CopyData( inPD, this->ErrorMap->VertexErrorMap.begin()->first, newId );
-    }
+    // What's left over is now spit out as a new polyline
+    vtkIdType newId = newLines->InsertNextCell(currentNumPts);
+    outCD->CopyData(inCD,firstVertexIndex,newId);
+
+    std::map<vtkIdType,vtkIdType> pointIdMap;
+    std::map<vtkIdType,vtkIdType>::iterator it;
 
-  // Clean up in preparation for the next line
-  this->PriorityQueue->Reset();
+    Polyline::Vertex* vertex = &(polyline->Vertices[0]);
+    while (vertex != NULL)
+      {
+      // points that are repeated within a single polyline are represented by
+      // only one point instance
+      it = pointIdMap.find(vertex->id);
+      if (it == pointIdMap.end())
+        {
+        newId = newPts->InsertNextPoint( inputPoints->GetPoint( vertex->id ) );
+        newLines->InsertCellPoint( newId );
+        outPD->CopyData( inPD, vertex->id, newId );
+        }
+      else
+        newLines->InsertCellPoint( it->second );
+
+      vertex = vertex->next;
+      }
+
+    delete polyline;
+    }
 
   // Create output and clean up
   output->SetPoints( newPts );
@@ -199,75 +270,6 @@ int vtkDecimatePolylineFilter::RequestData(
 
   return 1;
 }
-//---------------------------------------------------------------------
-int vtkDecimatePolylineFilter::GetPrev( int iId )
-{
-  std::map< int, double >::iterator it = this->ErrorMap->VertexErrorMap.find( iId );
-
-  if( it == this->ErrorMap->VertexErrorMap.begin() )
-    {
-    if( this->Closed )
-      {
-      it = this->ErrorMap->VertexErrorMap.end();
-      --it;
-      return it->first;
-      }
-    else
-      {
-      return iId;
-      }
-    }
-  else
-    {
-    --it;
-    return it->first;
-    }
-}
-//---------------------------------------------------------------------
-int vtkDecimatePolylineFilter::GetNext( int iId )
-{
-  std::map< int, double >::iterator
-    it = this->ErrorMap->VertexErrorMap.find( iId );
-  std::map< int, double >::iterator
-    end_it = this->ErrorMap->VertexErrorMap.end();
-  --end_it;
-  if( it == end_it )
-    {
-    if( this->Closed )
-      {
-      return this->ErrorMap->VertexErrorMap.begin()->first;
-      }
-    else
-      {
-      return iId;
-      }
-    }
-  else
-    {
-    ++it;
-    return it->first;
-    }
-}
-//---------------------------------------------------------------------
-void vtkDecimatePolylineFilter::UpdateError( vtkPolyData* input, int iId )
-{
-  int prev = GetPrev( iId );
-  int prev_prev = GetPrev( prev );
-  int next = GetNext( iId );
-  int next_next = GetNext( next );
-
-  double prev_error = ComputeError( input, prev_prev, prev, next );
-  this->ErrorMap->VertexErrorMap[prev] = prev_error;
-  this->PriorityQueue->DeleteId( prev );
-  this->PriorityQueue->Insert( prev_error, prev );
-
-  double next_error = ComputeError( input, prev, next, next_next );
-  this->ErrorMap->VertexErrorMap[next] = next_error;
-  this->PriorityQueue->DeleteId( next );
-  this->PriorityQueue->Insert( next_error, next );
-}
-//---------------------------------------------------------------------
-
 //---------------------------------------------------------------------
 void vtkDecimatePolylineFilter::PrintSelf(ostream& os, vtkIndent indent)
 {

Filters/Core/vtkDecimatePolylineFilter.h

@@ -23,12 +23,10 @@
 //
 // .SECTION Caveats
 // This algorithm is a very simple implementation that overlooks some
-// potential complexities. First, if a vertex is multiply connected,
-// meaning that it is used by multiple polylines, then the extra
-// topological constraints are ignored. Second, the error is not updated
-// as vertices are deleted (similar to iteratively computing a quadric
-// error metric). Thus, once calculated, the error is used to determine
-// which vertices are removed. This can produce less than optimal results.
+// potential complexities. For example, if a vertex is multiply connected,
+// meaning that it is used by multiple distinct polylines, then the extra
+// topological constraints are ignored. This can produce less than optimal
+// results.
 //
 // .SECTION See Also
 // vtkDecimate vtkDecimateProp vtkQuadricClustering vtkQuadricDecimation
@@ -77,17 +75,12 @@ protected:
 
   int RequestData(vtkInformation *, vtkInformationVector **, vtkInformationVector *);
 
-  double ComputeError( vtkPolyData* input, int prev, int id, int next );
-  void UpdateError( vtkPolyData* input, int iId );
-
-  int GetPrev( int iId );
-  int GetNext( int iId );
-
-  struct    vtkDecimatePolylineVertexErrorSTLMap;
-  vtkDecimatePolylineVertexErrorSTLMap*  ErrorMap;
+  class Polyline;
+  double ComputeError( vtkPolyData* input,
+                       Polyline* polyline,
+                       vtkIdType id );
 
   vtkSmartPointer< vtkPriorityQueue >   PriorityQueue;
-  bool                                  Closed;
   double                                TargetReduction;
   int                                   OutputPointsPrecision;