BUG: Handle CellData in Decimation Filter

Currently, cell data is not properly handled in the
SquaredEdgeLengthDecimation filter.  This patch adds entries to
the cell data array as faces are added in the Zip filter, and then
deletes unreferenced entries in a new member function added to
itk::Mesh.

Modules/Core/Mesh/include/itkMesh.h

@@ -359,6 +359,12 @@ class ITK_TEMPLATE_EXPORT Mesh : public PointSet<TPixelType, VDimension, TMeshTr
   const CellDataContainer *
   GetCellData() const;
 
+  /** Delete entries in m_CellDataContainer which do not have a corresponding
+   * entry in m_CellsContainer.
+   */
+  void
+  DeleteUnusedCellData();
+
 #if !defined(ITK_WRAPPING_PARSER)
   /**
    * Set/get the BoundaryAssignmentsContainer for a given dimension.

Modules/Core/Mesh/include/itkMesh.hxx

@@ -1129,6 +1129,30 @@ Mesh<TPixelType, VDimension, TMeshTraits>::Graft(const DataObject * data)
   // test on the container will prevent premature deletion of cells.
   this->m_CellsAllocationMethod = mesh->m_CellsAllocationMethod;
 }
+
+template <typename TPixelType, unsigned int VDimension, typename TMeshTraits>
+void
+Mesh<TPixelType, VDimension, TMeshTraits>::DeleteUnusedCellData()
+{
+
+  if (nullptr == this->GetCellData())
+  {
+    return;
+  }
+
+  std::vector<typename CellDataContainer::ElementIdentifier> cell_data_to_delete;
+  for (auto it = this->GetCellData()->Begin(); it != this->GetCellData()->End(); ++it)
+  {
+    if (!this->GetCells()->IndexExists(it.Index()))
+    {
+      cell_data_to_delete.push_back(it.Index());
+    }
+  }
+  for (const auto c : cell_data_to_delete)
+  {
+    this->GetCellData()->DeleteIndex(c);
+  }
+}
 } // end namespace itk
 
 #endif

Modules/Core/Mesh/test/CMakeLists.txt

@@ -34,6 +34,7 @@ itkCellInterfaceTest.cxx
 itkNewTest.cxx
 itkQuadrilateralCellTest.cxx
 itkTriangleCellTest.cxx
+itkMeshCellDataTest.cxx
 )
 
 set(ITKMesh-Test_LIBRARIES ${ITKMesh-Test_LIBRARIES})
@@ -182,6 +183,7 @@ itk_add_test(NAME itkNewTest COMMAND ITKMeshTestDriver itkNewTest)
 itk_add_test(NAME itkCellInterfaceTest COMMAND ITKMeshTestDriver itkCellInterfaceTest)
 itk_add_test(NAME itkTriangleCellTest COMMAND ITKMeshTestDriver itkTriangleCellTest)
 itk_add_test(NAME itkQuadrilateralCellTest COMMAND ITKMeshTestDriver itkQuadrilateralCellTest)
+itk_add_test(NAME itkMeshCellDataTest COMMAND ITKMeshTestDriver itkMeshCellDataTest)
 
 set_tests_properties(itkVTKPolyDataReaderTest2
    itkVTKPolyDataReaderBadTest0

Modules/Core/Mesh/test/itkMeshCellDataTest.cxx

@@ -0,0 +1,89 @@
+/*=========================================================================
+ *
+ *  Copyright Insight Software Consortium
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+
+#include <itkDefaultDynamicMeshTraits.h>
+#include <itkMesh.h>
+#include <itkTriangleCell.h>
+
+const unsigned int Dimension = 2;
+using TPixel = float;
+using TMeshTraits = itk::DefaultDynamicMeshTraits<TPixel>;
+using TMesh = itk::Mesh<TPixel, Dimension, TMeshTraits>;
+using TPoint = TMesh::PointType;
+using TCell = TMesh::CellType;
+using TTriangle = itk::TriangleCell<TCell>;
+
+int
+itkMeshCellDataTest(int, char *[])
+{
+
+  // 1------3------5
+  //   \    | \    |
+  //      \ |    \ |
+  // 0      2------4
+
+  const auto mesh = TMesh::New();
+  mesh->SetPoint(0, TPoint{ { 0.0, 0.0 } });
+  mesh->SetPoint(1, TPoint{ { 0.0, 1.0 } });
+  mesh->SetPoint(2, TPoint{ { 1.0, 0.0 } });
+  mesh->SetPoint(3, TPoint{ { 1.0, 1.0 } });
+  mesh->SetPoint(4, TPoint{ { 4.0, 0.0 } });
+  mesh->SetPoint(5, TPoint{ { 5.0, 1.0 } });
+
+  {
+    TCell::CellAutoPointer cellpointer;
+    cellpointer.TakeOwnership(new TTriangle);
+    cellpointer->SetPointId(0, 1);
+    cellpointer->SetPointId(1, 2);
+    cellpointer->SetPointId(2, 3);
+    mesh->SetCell(0, cellpointer);
+  }
+  {
+    TCell::CellAutoPointer cellpointer;
+    cellpointer.TakeOwnership(new TTriangle);
+    cellpointer->SetPointId(0, 3);
+    cellpointer->SetPointId(1, 2);
+    cellpointer->SetPointId(2, 4);
+    mesh->SetCell(1, cellpointer);
+  }
+  {
+    TCell::CellAutoPointer cellpointer;
+    cellpointer.TakeOwnership(new TTriangle);
+    cellpointer->SetPointId(0, 3);
+    cellpointer->SetPointId(1, 4);
+    cellpointer->SetPointId(2, 5);
+    mesh->SetCell(2, cellpointer);
+  }
+
+  itkAssertOrThrowMacro(6 == mesh->GetNumberOfPoints(), "Wrong number of points.");
+  itkAssertOrThrowMacro(3 == mesh->GetNumberOfCells(), "Wrong number of cells.");
+
+  // Add one extra entry in the cell data array
+  mesh->SetCellData(0, 0);
+  mesh->SetCellData(1, 1);
+  mesh->SetCellData(2, 2);
+  mesh->SetCellData(3, 3); // Extra entry
+
+  itkAssertOrThrowMacro(4 == mesh->GetCellData()->Size(), "Wrong number of entries in the cell data array.");
+
+  mesh->DeleteUnusedCellData();
+
+  itkAssertOrThrowMacro(3 == mesh->GetCellData()->Size(), "Wrong number of entries in the cell data array.");
+
+  return EXIT_SUCCESS;
+}

Modules/Core/QuadEdgeMesh/include/itkQuadEdgeMeshZipMeshFunction.hxx

@@ -81,6 +81,7 @@ QuadEdgeMeshZipMeshFunction<TMesh, TQEType>::Evaluate(QEType * e)
   OutputType VLeft = e->GetDestination();
   OutputType VRite = b->GetOrigin();
   bool       wasFacePresent = e->IsRightSet();
+  const auto rightFaceID = e->GetRight();
   OutputType resultingPointId = QEType::m_NoPoint;
 
   // We should be cautious and consider the case when the very
@@ -203,6 +204,13 @@ QuadEdgeMeshZipMeshFunction<TMesh, TQEType>::Evaluate(QEType * e)
   if (wasFacePresent)
   {
     this->m_Mesh->AddFace(b);
+    if (nullptr != this->m_Mesh->GetCellData())
+    {
+      if (this->m_Mesh->GetCellData()->IndexExists(rightFaceID))
+      {
+        this->m_Mesh->SetCellData(b->GetLeft(), this->m_Mesh->GetCellData()->ElementAt(rightFaceID));
+      }
+    }
   }
 
   this->m_Mesh->Modified();

Modules/Filtering/QuadEdgeMeshFiltering/include/itkDecimationQuadEdgeMeshFilter.h

@@ -88,6 +88,7 @@ class DecimationQuadEdgeMeshFilter : public QuadEdgeMeshToQuadEdgeMeshFilter<TIn
     } while (!IsCriterionSatisfied());
 
     this->GetOutput()->SqueezePointsIds();
+    this->GetOutput()->DeleteUnusedCellData();
   }
 
   virtual void

Modules/Filtering/QuadEdgeMeshFiltering/test/itkSquaredEdgeLengthDecimationQuadEdgeMeshFilterTest.cxx

@@ -60,6 +60,13 @@ itkSquaredEdgeLengthDecimationQuadEdgeMeshFilterTest(int argc, char * argv[])
 
   MeshType::Pointer mesh = reader->GetOutput();
 
+  for (auto it = mesh->GetCells()->Begin(); it != mesh->GetCells()->End(); ++it)
+  {
+    mesh->SetCellData(it.Index(), 25);
+  }
+  itkAssertOrThrowMacro(mesh->GetNumberOfCells() == mesh->GetCellData()->Size(),
+                        "Incorrect number of elements in the cell data array.");
+
   using CriterionType = itk::NumberOfFacesCriterion<MeshType>;
 
   using DecimationType = itk::SquaredEdgeLengthDecimationQuadEdgeMeshFilter<MeshType, MeshType, CriterionType>;
@@ -77,6 +84,9 @@ itkSquaredEdgeLengthDecimationQuadEdgeMeshFilterTest(int argc, char * argv[])
   decimate->SetCriterion(criterion);
   decimate->Update();
 
+  itkAssertOrThrowMacro(decimate->GetOutput()->GetNumberOfCells() == decimate->GetOutput()->GetCellData()->Size(),
+                        "Incorrect number of elements in the cell data array.");
+
   // ** WRITE OUTPUT **
   WriterType::Pointer writer = WriterType::New();
   writer->SetInput(decimate->GetOutput());