Update the nonlinear subdivision algo to support degenerate cells

Common/DataModel/Testing/Cxx/CMakeLists.txt

@@ -127,6 +127,7 @@ vtk_add_test_cxx(vtkCommonDataModelCxxTests data_tests
   TestPolyhedronDecompose.cxx,NO_DATA
   TestSmoothErrorMetric.cxx
   TestQuadraticPolygonFilters.cxx
+  TestNonlinearSubdivisionOfDegenerateCells.cxx,NO_DATA
   )
 # add to the list but don't define a test
 list(APPEND data_tests

Common/DataModel/Testing/Cxx/TestNonlinearSubdivisionOfDegenerateCells.cxx

+// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
+// SPDX-License-Identifier: BSD-3-Clause
+#include "vtkActor.h"
+#include "vtkDataSetMapper.h"
+#include "vtkDataSetSurfaceFilter.h"
+#include "vtkPoints.h"
+#include "vtkProperty.h"
+#include "vtkRegressionTestImage.h"
+#include "vtkRenderWindow.h"
+#include "vtkRenderWindowInteractor.h"
+#include "vtkRenderer.h"
+#include "vtkUnstructuredGrid.h"
+#include "vtkUnstructuredGridGeometryFilter.h"
+#include <vtkCamera.h>
+
+int TestNonlinearSubdivisionOfDegenerateCells(int argc, char* argv[])
+{
+  vtkNew<vtkUnstructuredGrid> grid;
+  vtkNew<vtkPoints> points;
+  points->SetNumberOfPoints(5);
+  points->SetPoint(0, 0, 0, 0);
+  points->SetPoint(1, 1, 0, 0);
+  points->SetPoint(2, 0, 1, 0);
+  points->SetPoint(3, 1, 1, 0);
+  points->SetPoint(4, 0.5, 0.5, 0.5);
+  grid->SetPoints(points);
+  vtkIdType connectivity1[8] = { 4, 1, 3, 4, 0, 0, 0, 0 };
+  vtkIdType connectivity2[8] = { 0, 0, 0, 0, 3, 4, 4, 2 };
+  grid->InsertNextCell(VTK_LAGRANGE_HEXAHEDRON, 8, connectivity1);
+  grid->InsertNextCell(VTK_LAGRANGE_HEXAHEDRON, 8, connectivity2);
+
+  // extract surface
+  vtkNew<vtkDataSetSurfaceFilter> surfaceFilter;
+  surfaceFilter->SetInputData(grid);
+  surfaceFilter->SetNonlinearSubdivisionLevel(3);
+  surfaceFilter->PassThroughCellIdsOff();
+  surfaceFilter->PassThroughPointIdsOff();
+  surfaceFilter->FastModeOn();
+  surfaceFilter->Update();
+
+  // Create a mapper and actor
+  vtkNew<vtkDataSetMapper> mapper;
+  mapper->SetInputConnection(surfaceFilter->GetOutputPort());
+
+  vtkNew<vtkActor> actor;
+  actor->SetMapper(mapper);
+
+  vtkSmartPointer<vtkProperty> prop = vtkSmartPointer<vtkProperty>::New();
+  prop->LightingOff();
+  prop->SetRepresentationToSurface();
+  prop->EdgeVisibilityOff();
+  prop->SetOpacity(0.5);
+
+  // set property
+  actor->SetProperty(prop);
+
+  // Visualize
+  vtkNew<vtkRenderer> renderer;
+  vtkNew<vtkRenderWindow> renderWindow;
+  renderWindow->AddRenderer(renderer);
+  vtkNew<vtkRenderWindowInteractor> renderWindowInteractor;
+  renderWindowInteractor->SetRenderWindow(renderWindow);
+
+  renderer->AddActor(actor);
+  renderer->ResetCamera();
+  renderer->GetActiveCamera()->Azimuth(30);
+  renderer->GetActiveCamera()->Elevation(10);
+  renderWindow->SetSize(600, 600);
+  renderWindow->Render();
+
+  int retVal = vtkRegressionTestImage(renderWindow);
+  if (retVal == vtkRegressionTester::DO_INTERACTOR)
+  {
+    renderWindowInteractor->Start();
+    retVal = vtkRegressionTester::PASSED;
+  }
+  return (retVal == vtkRegressionTester::PASSED) ? EXIT_SUCCESS : EXIT_FAILURE;
+}

Common/DataModel/Testing/Data/Baseline/TestNonlinearSubdivisionOfDegenerateCells.png.sha512

+1f880bee9007a2ae1fee5bf0fd300c3754692b08f4b299d839afd3b22f1187379368aa07f71c26e2267a0a71075beb712158de6b07da5e5b0fc62f031f910f26

Documentation/release/dev/subdivision-now-supports-degenerate-cells.md

+## Revisit subdivision of nonlinear cells in vtkDataSetSurfaceFilter
+
+In a recent change, the cell->triangulate function has been updated to return local point ids instead of global ids.
+Today, we take advantage of this to use local ids in the nonlinear subdivision algorithms as well.
+
+The main motivation of this change is that working with local ids helps supporting nonlinear subdivisions of degenerate cells.
+Indeed, the subdivision algorithm is now applied in the parametric space ( that is not degenerated), instead of the physical space.
+Once all the subdivisions are done, we move back to the physical space.

Filters/Geometry/vtkDataSetSurfaceFilter.cxx

@@ -292,6 +292,10 @@ public:
   }
   vtkIdType FindEdge(vtkIdType endpoint1, vtkIdType endpoint2)
   {
+    if (endpoint1 == endpoint2)
+    {
+      return endpoint1;
+    }
     if (endpoint1 > endpoint2)
       std::swap(endpoint1, endpoint2);
     MapType::iterator iter = Map.find(std::make_pair(endpoint1, endpoint2));
@@ -1455,7 +1459,7 @@ int vtkDataSetSurfaceFilter::UnstructuredGridExecuteInternal(
 
   // These are for subdividing quadratic cells
   std::vector<double> parametricCoords;
-  std::unique_ptr<vtkEdgeInterpolationMap> localEdgeMap;
+  std::unique_ptr<vtkEdgeInterpolationMap> localEdgeMap(new vtkEdgeInterpolationMap());
   vtkIdList* outPts;
   vtkIdList* pts2;
 
@@ -2028,22 +2032,39 @@ int vtkDataSetSurfaceFilter::UnstructuredGridExecuteInternal(
         }
       }
 
+      bool isDegenerateCell = false;
+      auto isDegeneratedSubTriangle = [&](vtkIdType ii) {
+        return outPts->GetId(pts->GetId(ii)) == outPts->GetId(pts->GetId(ii + 1)) ||
+          outPts->GetId(pts->GetId(ii)) == outPts->GetId(pts->GetId(ii + 2)) ||
+          outPts->GetId(pts->GetId(ii + 1)) == outPts->GetId(pts->GetId(ii + 2));
+      };
+
       // Do any further subdivision if necessary.
       if (this->NonlinearSubdivisionLevel > 1 && pc)
       {
+        for (i = 0; i < pts->GetNumberOfIds(); i += 3)
+        {
+          if (isDegeneratedSubTriangle(i))
+          {
+            isDegenerateCell = true;
+            break;
+          }
+        }
+
         vtkIdType maxNumberOfIds =
           std::pow(4, this->NonlinearSubdivisionLevel - 1) * pts->GetNumberOfIds();
         pts2->Allocate(maxNumberOfIds);
         // We are going to need parametric coordinates to further subdivide.
         parametricCoords.resize(maxNumberOfIds * 3);
-        for (i = 0; i < numFacePts * 3; i++)
-        {
-          parametricCoords[i] = pc[i];
-        }
+        std::copy(&pc[0], &pc[0] + numFacePts * 3, parametricCoords.begin());
 
         // localEdgeMap is simular to this->EdgeMap, but only stores local ids
         localEdgeMap->clear();
 
+        auto isEqualTo1Or0 = [](double a, double e = 1e-10) {
+          return (std::abs(a) <= e) || (std::abs(a - 1) <= e);
+        };
+
         vtkIdType localIdCpt = numFacePts;
         vtkIdType pt1, pt2, id;
         vtkIdType inPts[6];
@@ -2052,6 +2073,37 @@ int vtkDataSetSurfaceFilter::UnstructuredGridExecuteInternal(
         for (j = 1; j < this->NonlinearSubdivisionLevel; j++)
         {
           pts2->Reset();
+          if (isDegenerateCell)
+          {
+            // For degenerate cells, we can have multiple parametric points linked to the same
+            // output point. But we need to select a single one. The rule is to give priority to
+            // the points that are on the contour of the parametric space. This is necessary for
+            // connecting adjacent cells. The way we give this priority is by calling
+            // this->EdgeMap->FindEgde/AddEdge for those points first. So a first iteration over pts
+            // is performed to add those points. During the second iteration (the one not specific
+            // to degenerate cells), when trying to add a duplicate point, the edge map will return
+            // the output id of the already existing point.
+            double coords[3];
+            for (i = 0; i < pts->GetNumberOfIds(); i += 3)
+            {
+              for (k = 0; k < 3; k++)
+              {
+                pt1 = pts->GetId(i + k);
+                pt2 = pts->GetId(i + ((k < 2) ? (k + 1) : 0));
+                {
+                  coords[0] = 0.5 * (parametricCoords[pt1 * 3] + parametricCoords[pt2 * 3]);
+                  coords[1] = 0.5 * (parametricCoords[pt1 * 3 + 1] + parametricCoords[pt2 * 3 + 1]);
+                  coords[2] = 0.5 * (parametricCoords[pt1 * 3 + 2] + parametricCoords[pt2 * 3 + 2]);
+                  if (isEqualTo1Or0(coords[0]) || isEqualTo1Or0(coords[1]))
+                  {
+                    this->GetInterpolatedPointId(outPts->GetId(pt1), outPts->GetId(pt2), input,
+                      cell, coords, weights.data(), newPts, outputPD);
+                  }
+                }
+              }
+            }
+          }
+
           // Each triangle will be split into 4 triangles.
           for (i = 0; i < pts->GetNumberOfIds(); i += 3)
           {
@@ -2104,6 +2156,10 @@ int vtkDataSetSurfaceFilter::UnstructuredGridExecuteInternal(
       }
       for (i = 0; i < pts->GetNumberOfIds(); i += 3)
       {
+        if (isDegenerateCell && isDegeneratedSubTriangle(i))
+        {
+          continue; // Do not record the degenerate triangle
+        }
         newPolys->InsertNextCell(3);
         newPolys->InsertCellPoint(outPts->GetId(pts->GetId(i)));
         newPolys->InsertCellPoint(outPts->GetId(pts->GetId(i + 1)));