From ce897a721afd192a247f1fd51f0c6f63253e2c98 Mon Sep 17 00:00:00 2001
From: Thomas Galland <thomas.galland@kitware.com>
Date: Thu, 2 Mar 2023 08:42:26 +0100
Subject: [PATCH] Add vtkPolyhedronUtilities and Decompose method
Add a new method allowing to decompose polyhedron into tets.
This method will generate new points on polyhedron and it's purpose
is to improve the result of subsequent filters (e.g. contours) on
polyhedrons with concave faces.
The vtkPolyhedronUtilities class is added to contain this method, and
potential future ones working on polyhedrons for the same kind of purpose.
This class can be seen as a toolbox for vtkPolyhedron.
---
Common/DataModel/CMakeLists.txt | 1 +
Common/DataModel/Testing/Cxx/CMakeLists.txt | 1 +
.../Testing/Cxx/TestPolyhedronDecompose.cxx | 376 ++++++++++++++++++
.../TestPolyhedronDecompose.png.sha512 | 1 +
Common/DataModel/vtkPolyhedron.cxx | 18 +-
Common/DataModel/vtkPolyhedron.h | 5 +-
Common/DataModel/vtkPolyhedronUtilities.cxx | 356 +++++++++++++++++
Common/DataModel/vtkPolyhedronUtilities.h | 65 +++
.../release/dev/add-polyhedron-utilities.md | 8 +
9 files changed, 820 insertions(+), 11 deletions(-)
create mode 100644 Common/DataModel/Testing/Cxx/TestPolyhedronDecompose.cxx
create mode 100644 Common/DataModel/Testing/Data/Baseline/TestPolyhedronDecompose.png.sha512
create mode 100644 Common/DataModel/vtkPolyhedronUtilities.cxx
create mode 100644 Common/DataModel/vtkPolyhedronUtilities.h
create mode 100644 Documentation/release/dev/add-polyhedron-utilities.md
diff --git a/Common/DataModel/CMakeLists.txt b/Common/DataModel/CMakeLists.txt
index eb84a44d9a4..3b34d51a08c 100644
--- a/Common/DataModel/CMakeLists.txt
+++ b/Common/DataModel/CMakeLists.txt
@@ -191,6 +191,7 @@ set(classes
vtkPolyVertex
vtkPolygon
vtkPolyhedron
+ vtkPolyhedronUtilities
vtkPyramid
vtkQuad
vtkQuadraticEdge
diff --git a/Common/DataModel/Testing/Cxx/CMakeLists.txt b/Common/DataModel/Testing/Cxx/CMakeLists.txt
index eed3ae9d035..4ed4077cde4 100644
--- a/Common/DataModel/Testing/Cxx/CMakeLists.txt
+++ b/Common/DataModel/Testing/Cxx/CMakeLists.txt
@@ -120,6 +120,7 @@ vtk_add_test_cxx(vtkCommonDataModelCxxTests data_tests
TestMeanValueCoordinatesInterpolation1.cxx
TestMeanValueCoordinatesInterpolation2.cxx
TestPolyhedron4.cxx,NO_VALID,NO_OUTPUT
+ TestPolyhedronDecompose.cxx,NO_DATA
TestSmoothErrorMetric.cxx
TestQuadraticPolygonFilters.cxx
)
diff --git a/Common/DataModel/Testing/Cxx/TestPolyhedronDecompose.cxx b/Common/DataModel/Testing/Cxx/TestPolyhedronDecompose.cxx
new file mode 100644
index 00000000000..7af48b2151d
--- /dev/null
+++ b/Common/DataModel/Testing/Cxx/TestPolyhedronDecompose.cxx
@@ -0,0 +1,376 @@
+/*=========================================================================
+
+Program: Visualization Toolkit
+Module: TestPolyhedronDecompose.cxx
+
+Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
+All rights reserved.
+See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
+
+This software is distributed WITHOUT ANY WARRANTY; without even
+the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+PURPOSE. See the above copyright notice for more information.
+
+=========================================================================*/
+
+#include "vtkRegressionTestImage.h"
+#include <vtkActor.h>
+#include <vtkBitArray.h>
+#include <vtkCamera.h>
+#include <vtkCellData.h>
+#include <vtkContourFilter.h>
+#include <vtkDoubleArray.h>
+#include <vtkGeometryFilter.h>
+#include <vtkNew.h>
+#include <vtkPointData.h>
+#include <vtkPoints.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkPolyhedron.h>
+#include <vtkPolyhedronUtilities.h>
+#include <vtkProperty.h>
+#include <vtkRenderWindow.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkRenderer.h>
+#include <vtkSmartPointer.h>
+#include <vtkStringArray.h>
+#include <vtkUnstructuredGrid.h>
+
+#include <cmath>
+#include <limits>
+
+namespace
+{
+vtkSmartPointer<vtkPolyhedron> MakePolyhedron1();
+vtkSmartPointer<vtkPolyhedron> MakePolyhedron2();
+
+template <typename T>
+bool testValue(const std::string& value, T gotVal, T expectedVal);
+}
+
+//------------------------------------------------------------------------------
+int TestPolyhedronDecompose(int argc, char* argv[])
+{
+ ////////// Setup data objects //////////
+
+ // Create neighboring polyhedrons
+ auto polyhedron1 = ::MakePolyhedron1();
+ auto polyhedron2 = ::MakePolyhedron2();
+
+ // Add some cell data
+ vtkNew<vtkDoubleArray> cellArray;
+ cellArray->SetNumberOfValues(1);
+ cellArray->SetName("Cell array");
+ cellArray->SetValue(0, 1.5);
+
+ vtkNew<vtkCellData> cellData;
+ cellData->AddArray(cellArray);
+
+ // Add some point data
+ constexpr std::array<double, 8> doubleValues = { 2, 5, 2, 2, 2, 3, 2, 3 };
+ const std::array<std::string, 8> stringValues = { "A", "A", "A", "A", "A", "A", "A", "A" };
+ constexpr std::array<int, 8> bitValues = { 0, 0, 0, 0, 0, 0, 0, 0 };
+
+ vtkNew<vtkDoubleArray> pointArrayDouble; // Will be dispatched
+ pointArrayDouble->SetNumberOfValues(8);
+ pointArrayDouble->SetName("Doubles");
+
+ vtkNew<vtkStringArray> pointArrayString; // Will not
+ pointArrayString->SetNumberOfValues(8);
+ pointArrayString->SetName("Strings");
+
+ vtkNew<vtkBitArray> pointArrayBits; // Will not
+ pointArrayBits->SetNumberOfValues(8);
+ pointArrayBits->SetName("Bits");
+
+ for (int i = 0; i < 8; i++)
+ {
+ pointArrayDouble->SetValue(i, doubleValues[i]);
+ pointArrayString->SetValue(i, stringValues[i]);
+ pointArrayBits->SetValue(i, bitValues[i]);
+ }
+
+ vtkNew<vtkPointData> pointData;
+ pointData->AddArray(pointArrayDouble);
+ pointData->AddArray(pointArrayString);
+ pointData->AddArray(pointArrayBits);
+
+ // Decompose polyhedra
+ auto decomposedUG1 = vtkPolyhedronUtilities::Decompose(polyhedron1, pointData, 0, cellData);
+ auto decomposedUG2 = vtkPolyhedronUtilities::Decompose(polyhedron2, pointData, 0, cellData);
+
+ ////////// Test geometry //////////
+
+ // New number of pts = original pts + face barycenters + cell barycenter
+ auto numberOfPts = decomposedUG1->GetNumberOfPoints();
+ if (!testValue("number of points", numberOfPts, vtkIdType(8 + 6 + 1)))
+ {
+ return EXIT_FAILURE;
+ }
+
+ // New number of cells = original nb of faces * 4
+ auto numberOfCells = decomposedUG1->GetNumberOfCells();
+ if (!testValue("number of cells", numberOfCells, vtkIdType(6 * 4)))
+ {
+ return EXIT_FAILURE;
+ }
+
+ ////////// Test data //////////
+
+ auto pointDataDec = decomposedUG1->GetPointData();
+ auto cellDataDec = decomposedUG1->GetCellData();
+
+ // Test barycenters point data
+ // Face barycenter: mean value of face point data
+ // Cell barycenter (last one): mean value of face barycenters point data
+ vtkDoubleArray* doubleArray =
+ vtkDoubleArray::SafeDownCast(pointDataDec->GetAbstractArray("Doubles"));
+ if (!doubleArray)
+ {
+ std::cerr << "Unable to retrieve \"Doubles\" point data." << std::endl;
+ }
+
+ if (!testValue("point data (\"Doubles\") nb of tuples", doubleArray->GetNumberOfTuples(),
+ decomposedUG1->GetNumberOfPoints()))
+ {
+ return EXIT_FAILURE;
+ }
+
+ constexpr double expectedValues[7] = { 2.75, 3, 2, 3.25, 2.25, 2.5, 2.625 };
+ for (int i = 0; i < 7; i++)
+ {
+ if (!testValue("point data (\"Doubles\") at index " + i, doubleArray->GetValue(i + 8),
+ expectedValues[i]))
+ {
+ return EXIT_FAILURE;
+ }
+ }
+
+ // vtkStringArray is not dispatched, check that the fallback initialized
+ // the values as an empty string
+ vtkStringArray* stringArray =
+ vtkStringArray::SafeDownCast(pointDataDec->GetAbstractArray("Strings"));
+ if (!stringArray)
+ {
+ std::cerr << "Unable to retrieve \"Strings\" point data." << std::endl;
+ }
+
+ if (!testValue("point data (\"Strings\") nb of tuples", stringArray->GetNumberOfTuples(),
+ decomposedUG1->GetNumberOfPoints()))
+ {
+ return EXIT_FAILURE;
+ }
+
+ for (int i = 0; i < 7; i++)
+ {
+ if (!testValue(
+ "point data (\"Strings\") at index " + i, stringArray->GetValue(i + 8), vtkStdString()))
+ {
+ return EXIT_FAILURE;
+ }
+ }
+
+ // vtkBitArray is not dispatched, check that the fallback initialized
+ // the values with 0
+ vtkBitArray* bitArray = vtkBitArray::SafeDownCast(pointDataDec->GetAbstractArray("Bits"));
+ if (!bitArray)
+ {
+ std::cerr << "Unable to retrieve \"Bits\" point data." << std::endl;
+ }
+
+ if (!testValue("point data (\"Bits\") nb of tuples", bitArray->GetNumberOfTuples(),
+ decomposedUG1->GetNumberOfPoints()))
+ {
+ return EXIT_FAILURE;
+ }
+
+ for (int i = 0; i < 7; i++)
+ {
+ if (!testValue("point data (\"Bits\") at index " + i, bitArray->GetValue(i + 8), 0))
+ {
+ return EXIT_FAILURE;
+ }
+ }
+
+ // Cell data should be copied to all cells
+ vtkDoubleArray* doubleArrayCells =
+ vtkDoubleArray::SafeDownCast(cellDataDec->GetAbstractArray("Cell array"));
+ if (!doubleArrayCells)
+ {
+ std::cerr << "Unable to retrieve \"Cell array\" cell data." << std::endl;
+ }
+
+ for (vtkIdType cellId = 0; cellId < decomposedUG1->GetNumberOfCells(); cellId++)
+ {
+ if (!testValue("Cell array", doubleArrayCells->GetValue(cellId), 1.5))
+ {
+ return EXIT_FAILURE;
+ }
+ }
+
+ ////////// Test contour //////////
+
+ // Extract contours from decomposed ugs
+ vtkNew<vtkContourFilter> contourFilter;
+ contourFilter->SetInputData(decomposedUG1);
+ contourFilter->SetInputArrayToProcess(
+ 0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "Doubles");
+ contourFilter->SetNumberOfContours(1);
+ contourFilter->SetValue(0, 3.5);
+ contourFilter->Update();
+
+ vtkNew<vtkPolyData> contour1;
+ contour1->DeepCopy(vtkPolyData::SafeDownCast(contourFilter->GetOutputDataObject(0)));
+
+ contourFilter->SetInputData(decomposedUG2);
+ contourFilter->Update();
+
+ vtkNew<vtkPolyData> contour2;
+ contour2->DeepCopy(vtkPolyData::SafeDownCast(contourFilter->GetOutputDataObject(0)));
+
+ // Extract surface from decomposed ugs for rendering
+ vtkNew<vtkGeometryFilter> filter;
+ filter->SetInputDataObject(decomposedUG1);
+ filter->Update();
+
+ vtkNew<vtkPolyData> ugSurface1;
+ ugSurface1->DeepCopy(vtkPolyData::SafeDownCast(filter->GetOutputDataObject(0)));
+
+ filter->SetInputDataObject(decomposedUG2);
+ filter->Update();
+
+ vtkNew<vtkPolyData> ugSurface2;
+ ugSurface2->DeepCopy(vtkPolyData::SafeDownCast(filter->GetOutputDataObject(0)));
+ ugSurface2->DeepCopy(vtkPolyData::SafeDownCast(filter->GetOutputDataObject(0)));
+
+ // Mappers
+ vtkNew<vtkPolyDataMapper> ugMapper1;
+ ugMapper1->SetInputData(ugSurface1);
+ vtkNew<vtkPolyDataMapper> ugMapper2;
+ ugMapper2->SetInputData(ugSurface2);
+ vtkNew<vtkPolyDataMapper> contourMapper1;
+ contourMapper1->SetInputData(contour1);
+ vtkNew<vtkPolyDataMapper> contourMapper2;
+ contourMapper2->SetInputData(contour2);
+
+ // Actors
+ vtkNew<vtkActor> ugActor1;
+ ugActor1->SetMapper(ugMapper1);
+ ugActor1->GetProperty()->SetOpacity(0.1);
+ vtkNew<vtkActor> ugActor2;
+ ugActor2->SetMapper(ugMapper2);
+ ugActor2->GetProperty()->SetOpacity(0.1);
+ vtkNew<vtkActor> contourActor1;
+ contourActor1->SetMapper(contourMapper1);
+ vtkNew<vtkActor> contourActor2;
+ contourActor2->SetMapper(contourMapper2);
+
+ // Renderer
+ vtkNew<vtkRenderer> renderer;
+ renderer->AddActor(ugActor1);
+ renderer->AddActor(ugActor2);
+ renderer->AddActor(contourActor1);
+ renderer->AddActor(contourActor2);
+
+ // Camera
+ renderer->GetActiveCamera()->Azimuth(135);
+ renderer->ResetCamera();
+
+ // Render window
+ vtkNew<vtkRenderWindow> renderWindow;
+ renderWindow->AddRenderer(renderer);
+ renderWindow->SetSize(300, 300);
+
+ // Interactor
+ vtkNew<vtkRenderWindowInteractor> renderWindowInteractor;
+ renderWindowInteractor->SetRenderWindow(renderWindow);
+
+ // Regression image testing
+ renderWindow->Render();
+ int retVal = vtkRegressionTestImage(renderWindow);
+
+ if (retVal == vtkRegressionTester::DO_INTERACTOR)
+ {
+ renderWindowInteractor->Start();
+ retVal = vtkRegressionTester::PASSED;
+ }
+
+ return (retVal == vtkRegressionTester::PASSED) ? EXIT_SUCCESS : EXIT_FAILURE;
+}
+
+namespace
+{
+//------------------------------------------------------------------------------
+vtkSmartPointer<vtkPolyhedron> MakePolyhedron1()
+{
+ vtkSmartPointer<vtkPolyhedron> polyhedron = vtkSmartPointer<vtkPolyhedron>::New();
+
+ // Point Ids
+ for (int i = 0; i < 8; ++i)
+ {
+ polyhedron->GetPointIds()->InsertNextId(i);
+ }
+
+ // Points
+ polyhedron->GetPoints()->InsertNextPoint(2.5, -7.5, 2.5);
+ polyhedron->GetPoints()->InsertNextPoint(5.31, -5.31, 4.68);
+ polyhedron->GetPoints()->InsertNextPoint(2.5, -2.5, 2.5);
+ polyhedron->GetPoints()->InsertNextPoint(7.5, -2.5, 2.5);
+ polyhedron->GetPoints()->InsertNextPoint(2.5, -7.5, 7.5);
+ polyhedron->GetPoints()->InsertNextPoint(6.25, -6.25, 6.25);
+ polyhedron->GetPoints()->InsertNextPoint(2.5, -2.5, 7.5);
+ polyhedron->GetPoints()->InsertNextPoint(6.25, -3.75, 6.25);
+
+ // Faces
+ vtkIdType faces[31] = { 6, 4, 0, 1, 3, 2, 4, 0, 4, 5, 1, 4, 0, 2, 6, 4, 4, 1, 5, 7, 3, 4, 3, 7, 6,
+ 2, 4, 4, 6, 7, 5 };
+
+ polyhedron->SetFaces(faces);
+ polyhedron->Initialize();
+
+ return polyhedron;
+}
+
+//------------------------------------------------------------------------------
+vtkSmartPointer<vtkPolyhedron> MakePolyhedron2()
+{
+ vtkSmartPointer<vtkPolyhedron> polyhedron = vtkSmartPointer<vtkPolyhedron>::New();
+
+ // Point Ids
+ for (int i = 0; i < 8; ++i)
+ {
+ polyhedron->GetPointIds()->InsertNextId(i);
+ }
+
+ // Points
+ polyhedron->GetPoints()->InsertNextPoint(2.5, -7.5, 2.5);
+ polyhedron->GetPoints()->InsertNextPoint(5.31, -5.31, 4.68);
+ polyhedron->GetPoints()->InsertNextPoint(2.5, -12.5, 2.5);
+ polyhedron->GetPoints()->InsertNextPoint(7.5, -12.5, 2.5);
+ polyhedron->GetPoints()->InsertNextPoint(2.5, -7.5, 7.5);
+ polyhedron->GetPoints()->InsertNextPoint(6.25, -6.25, 6.25);
+ polyhedron->GetPoints()->InsertNextPoint(2.5, -12.5, 7.5);
+ polyhedron->GetPoints()->InsertNextPoint(6.25, -13.75, 6.25);
+
+ // Faces
+ vtkIdType faces[31] = { 6, 4, 2, 3, 1, 0, 4, 1, 5, 4, 0, 4, 4, 6, 2, 0, 4, 3, 7, 5, 1, 4, 2, 6, 7,
+ 3, 4, 5, 7, 6, 4 };
+
+ polyhedron->SetFaces(faces);
+ polyhedron->Initialize();
+
+ return polyhedron;
+}
+
+//------------------------------------------------------------------------------
+template <typename T>
+bool testValue(const std::string& value, T gotVal, T expectedVal)
+{
+ if (gotVal != expectedVal)
+ {
+ std::cerr << "Wrong " << value << ". Got " << gotVal << ", expected is : " << expectedVal
+ << std::endl;
+ return false;
+ }
+ return true;
+}
+}
diff --git a/Common/DataModel/Testing/Data/Baseline/TestPolyhedronDecompose.png.sha512 b/Common/DataModel/Testing/Data/Baseline/TestPolyhedronDecompose.png.sha512
new file mode 100644
index 00000000000..f5810df3b8b
--- /dev/null
+++ b/Common/DataModel/Testing/Data/Baseline/TestPolyhedronDecompose.png.sha512
@@ -0,0 +1 @@
+8d5245890941c61ab90317ad05068fcd6e26e18a14a015f2b7d71237d1995802e96fd9cd8a05fa44db0700b71e001178aa3f1b3f42c2b3c5053bacc4bbf2f09f
diff --git a/Common/DataModel/vtkPolyhedron.cxx b/Common/DataModel/vtkPolyhedron.cxx
index d5f4d0ca2ff..dc896018aa3 100644
--- a/Common/DataModel/vtkPolyhedron.cxx
+++ b/Common/DataModel/vtkPolyhedron.cxx
@@ -46,9 +46,6 @@ vtkStandardNewMacro(vtkPolyhedron);
// Special typedef
typedef std::vector<vtkIdType> vtkIdVectorType;
-class vtkPointIdMap : public std::map<vtkIdType, vtkIdType>
-{
-};
// an edge consists of two id's and their order
// is *not* important. To that end special hash and
@@ -1349,8 +1346,9 @@ vtkIdType vtkPolyhedron::GetPointToIncidentFaces(vtkIdType pointId, const vtkIdT
return this->ValenceAtPoint[pointId];
}
-bool IntersectWithContour(vtkCell* cell, vtkDataArray* pointScalars, vtkPointIdMap* pointIdMap,
- double value, std::function<bool(double, double)>& compare, bool& allTrue)
+bool IntersectWithContour(vtkCell* cell, vtkDataArray* pointScalars,
+ vtkPolyhedron::vtkPointIdMap* pointIdMap, double value,
+ std::function<bool(double, double)>& compare, bool& allTrue)
{
allTrue = true;
bool allFalse = true;
@@ -1571,8 +1569,8 @@ int TriangulatePolygonAt(vtkCell* polygon, int offset, vtkIdList* triIds)
return nPoints - 2;
}
-void CalculateAngles(const vtkIdType* tri, vtkPoints* phPoints, const vtkPointIdMap* pointIdMap,
- double& minAngle, double& maxAngle)
+void CalculateAngles(const vtkIdType* tri, vtkPoints* phPoints,
+ const vtkPolyhedron::vtkPointIdMap* pointIdMap, double& minAngle, double& maxAngle)
{
vtkIdType idx0 = tri[0];
vtkIdType idx1 = tri[1];
@@ -1620,7 +1618,7 @@ void CalculateAngles(const vtkIdType* tri, vtkPoints* phPoints, const vtkPointId
}
void TriangulatePolygon(vtkCell* polygon, FaceVector& faces, vtkIdList* triIds, vtkPoints* phPoints,
- vtkPointIdMap* pointIdMap)
+ vtkPolyhedron::vtkPointIdMap* pointIdMap)
{
// attempt a fan triangulation for each point on the polygon and choose the
// fan triangulation with the lowest range in internal angles differing from 60 degrees
@@ -1666,7 +1664,7 @@ void TriangulatePolygon(vtkCell* polygon, FaceVector& faces, vtkIdList* triIds,
}
void TriangulateFace(vtkCell* face, FaceVector& faces, vtkIdList* triIds, vtkPoints* phPoints,
- vtkPointIdMap* pointIdMap)
+ vtkPolyhedron::vtkPointIdMap* pointIdMap)
{
switch (face->GetCellType())
{
@@ -1710,7 +1708,7 @@ bool CheckNonManifoldTriangulation(EdgeFaceSetMap& edgeFaceMap)
}
bool GetContourPoints(double value, vtkPolyhedron* cell,
- vtkPointIdMap* pointIdMap, // from global id to local cell id
+ vtkPolyhedron::vtkPointIdMap* pointIdMap, // from global id to local cell id
FaceEdgesVector& faceEdgesVector, EdgeFaceSetMap& edgeFaceMap, EdgeSet& originalEdges,
std::vector<std::vector<vtkIdType>>& oririginalFaceTriFaceMap,
PointIndexEdgeMultiMap& contourPointEdgeMultiMap, EdgePointIndexMap& edgeContourPointMap,
diff --git a/Common/DataModel/vtkPolyhedron.h b/Common/DataModel/vtkPolyhedron.h
index 3ea43540f64..dfc6aecb07f 100644
--- a/Common/DataModel/vtkPolyhedron.h
+++ b/Common/DataModel/vtkPolyhedron.h
@@ -46,7 +46,6 @@ class vtkQuad;
class vtkTetra;
class vtkPolygon;
class vtkLine;
-class vtkPointIdMap;
class vtkIdToIdVectorMapType;
class vtkIdToIdMapType;
class vtkEdgeTable;
@@ -58,6 +57,8 @@ class vtkPointLocator;
class VTKCOMMONDATAMODEL_EXPORT vtkPolyhedron : public vtkCell3D
{
public:
+ typedef std::map<vtkIdType, vtkIdType> vtkPointIdMap;
+
///@{
/**
* Standard new methods.
@@ -345,6 +346,8 @@ protected:
private:
vtkPolyhedron(const vtkPolyhedron&) = delete;
void operator=(const vtkPolyhedron&) = delete;
+
+ friend class vtkPolyhedronUtilities;
};
//----------------------------------------------------------------------------
diff --git a/Common/DataModel/vtkPolyhedronUtilities.cxx b/Common/DataModel/vtkPolyhedronUtilities.cxx
new file mode 100644
index 00000000000..50e397b1fe9
--- /dev/null
+++ b/Common/DataModel/vtkPolyhedronUtilities.cxx
@@ -0,0 +1,356 @@
+/*=========================================================================
+
+ Program: Visualization Toolkit
+ Module: vtkPolyhedronUtilities.cxx
+
+ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
+ All rights reserved.
+ See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notice for more information.
+
+ =========================================================================*/
+
+#include "vtkPolyhedronUtilities.h"
+#include "vtkArrayDispatch.h"
+#include "vtkCellData.h"
+#include "vtkDataArray.h"
+#include "vtkPointData.h"
+#include "vtkPolyhedron.h"
+#include "vtkStringArray.h"
+#include "vtkUnstructuredGrid.h"
+
+namespace
+{
+//------------------------------------------------------------------------------
+struct InitWorker
+{
+ // Insert new tuple to the array and initialize all its components to 0
+ template <typename ArrayType>
+ void operator()(ArrayType* outArray)
+ {
+ using T = vtk::GetAPIType<ArrayType>;
+
+ std::vector<T> nextTuple(outArray->GetNumberOfComponents(), 0);
+ outArray->InsertNextTypedTuple(nextTuple.data());
+ }
+};
+
+//------------------------------------------------------------------------------
+// Fallback for not dispatched data arrays
+template <>
+void InitWorker::operator()(vtkDataArray* outArray)
+{
+ std::vector<double> nextTuple(outArray->GetNumberOfComponents(), 0);
+ outArray->InsertNextTuple(nextTuple.data());
+}
+
+//------------------------------------------------------------------------------
+// Fallback for string array
+template <>
+void InitWorker::operator()(vtkStringArray* outArray)
+{
+ vtkIdType nbOfComponents = outArray->GetNumberOfComponents();
+ vtkIdType nbOfValuesOld = outArray->GetNumberOfValues();
+ outArray->SetNumberOfValues(nbOfValuesOld + nbOfComponents);
+
+ for (int i = 0; i < nbOfComponents; i++)
+ {
+ outArray->SetValue(nbOfValuesOld + i, "");
+ }
+}
+
+//------------------------------------------------------------------------------
+// Fallback for all other arrays
+template <>
+void InitWorker::operator()(vtkAbstractArray* outArray)
+{
+ // Insert one "uninitialized" tuple
+ vtkIdType nbOfComponents = outArray->GetNumberOfComponents();
+ vtkIdType nbOfValuesOld = outArray->GetNumberOfValues();
+ outArray->SetNumberOfValues(nbOfValuesOld + nbOfComponents);
+}
+
+//------------------------------------------------------------------------------
+struct AccuWorker
+{
+ // Add the components of the given tuple of inArray to the components
+ // of the last tuple of outArray
+ template <typename ArrayType1, typename ArrayType2>
+ void operator()(ArrayType1* inArray, ArrayType2* outArray, vtkIdType inPtId)
+ {
+ using T = vtk::GetAPIType<ArrayType1>;
+
+ int nbOfComp = outArray->GetNumberOfComponents();
+ vtkIdType nbOfValues = outArray->GetNumberOfValues();
+
+ const auto inRange =
+ vtk::DataArrayValueRange(inArray, inPtId * nbOfComp, (inPtId + 1) * nbOfComp);
+ auto outRange = vtk::DataArrayValueRange(outArray, nbOfValues - nbOfComp, nbOfValues);
+
+ std::transform(inRange.cbegin(), inRange.cend(), outRange.cbegin(), outRange.begin(),
+ [=](T inValue, T outValue) { return outValue + inValue; });
+ }
+};
+
+//------------------------------------------------------------------------------
+struct DivWorker
+{
+ // Divide all the components of the last tuple of outArray with div
+ template <typename ArrayType>
+ void operator()(ArrayType* outArray, vtkIdType div)
+ {
+ using T = vtk::GetAPIType<ArrayType>;
+
+ int nbOfComp = outArray->GetNumberOfComponents();
+ vtkIdType nbOfValues = outArray->GetNumberOfValues();
+
+ auto range = vtk::DataArrayValueRange(outArray, nbOfValues - nbOfComp, nbOfValues);
+
+ std::transform(
+ range.cbegin(), range.cend(), range.begin(), [=](T inValue) { return inValue / div; });
+ }
+};
+}
+
+VTK_ABI_NAMESPACE_BEGIN
+
+//------------------------------------------------------------------------------
+vtkSmartPointer<vtkUnstructuredGrid> vtkPolyhedronUtilities::Decompose(
+ vtkPolyhedron* polyhedron, vtkPointData* inPd, vtkIdType cellId, vtkCellData* inCd)
+{
+ if (!polyhedron->GetPoints() || !polyhedron->GetNumberOfPoints())
+ {
+ return 0;
+ }
+
+ vtkNew<vtkUnstructuredGrid> outputGrid;
+ ::InitWorker initWorker;
+ ::AccuWorker accuWorker;
+ ::DivWorker divWorker;
+ typedef vtkArrayDispatch::DispatchByValueType<vtkArrayDispatch::AllTypes> Dispatcher;
+ typedef vtkArrayDispatch::Dispatch2BySameValueType<vtkArrayDispatch::AllTypes> Dispatcher2;
+
+ ////////// Compute barycenters and barycenters data //////////
+ // Here we iterate over each face and generate a new point (barycenter of the face).
+ // We also add new point data for the barycenter, that is the mean value of the face points data.
+ // XXX Consider rework this code in order to include the face and face points interations
+ // inside the workers in order to reduce the number of dispatchs (that are costly)
+
+ // Copy the point data to the output
+ vtkPointData* outPd = outputGrid->GetPointData();
+ outPd->DeepCopy(inPd);
+
+ // Global faces are faces with global point indexes
+ vtkIdType* globalFaces = polyhedron->GetFaces();
+ vtkIdType facesNb = globalFaces[0];
+ vtkIdType* globalFace = globalFaces + 1;
+ vtkIdType numberOfNewCells = 0; // Account for the number of cells of the output UG
+ vtkPolyhedron::vtkPointIdMap* pointIdMap = polyhedron->PointIdMap;
+
+ vtkNew<vtkPoints> barycenters;
+ // Iterate on each face to compute face barycenters and barycenters data (point data)
+ for (vtkIdType faceCount = 0; faceCount < facesNb; faceCount++)
+ {
+ vtkIdType nbFacePts = globalFace[0];
+
+ // Add a new value for each output array, init to 0.0
+ for (vtkIdType arrayId = 0; arrayId < outPd->GetNumberOfArrays(); arrayId++)
+ {
+ vtkAbstractArray* abstArray = outPd->GetAbstractArray(arrayId);
+ if (auto dataArray = vtkDataArray::SafeDownCast(abstArray))
+ {
+ if (!Dispatcher::Execute(dataArray, initWorker))
+ {
+ initWorker(dataArray); // Fallback for vtkDataArray subtypes
+ }
+ }
+ else if (auto stringArray = vtkStringArray::SafeDownCast(abstArray))
+ {
+ initWorker(stringArray); // Fallback for vtkStringArray
+ }
+ else
+ {
+ vtkWarningWithObjectMacro(nullptr,
+ "" << abstArray->GetName()
+ << ": array type is not supported. Values on new points will be undefined.");
+ initWorker(abstArray); // Fallback for all other types
+ }
+ }
+
+ std::array<double, 3> barycenter = { 0.0, 0.0, 0.0 };
+
+ for (vtkIdType i = 1; i <= nbFacePts; i++)
+ {
+ // Accumulate face points coordinates
+ auto globalPtId = globalFace[i];
+ auto localPtId = (*pointIdMap)[globalPtId];
+
+ std::array<double, 3> pt = { 0.0, 0.0, 0.0 };
+ polyhedron->GetPoints()->GetPoint(localPtId, pt.data());
+ barycenter[0] += pt[0];
+ barycenter[1] += pt[1];
+ barycenter[2] += pt[2];
+
+ // Accumulate barycenter new point data
+ for (vtkIdType arrayId = 0; arrayId < outPd->GetNumberOfArrays(); arrayId++)
+ {
+ vtkDataArray* inArray = vtkDataArray::SafeDownCast(inPd->GetAbstractArray(arrayId));
+ vtkDataArray* outArray = vtkDataArray::SafeDownCast(outPd->GetAbstractArray(arrayId));
+ if (inArray && outArray)
+ {
+ if (!Dispatcher2::Execute(inArray, outArray, accuWorker, globalPtId))
+ {
+ // We only need to fallback for vtkDataArray subtypes, other types
+ // are just initialised and no mean value is computed.
+ accuWorker(inArray, outArray, globalPtId);
+ }
+ }
+ }
+ }
+
+ // Compute the barycenter of the face
+ for (int i = 0; i < 3; i++)
+ {
+ barycenter[i] /= nbFacePts;
+ }
+
+ barycenters->InsertNextPoint(barycenter.data());
+
+ // Compute barycenter point data
+ for (vtkIdType arrayId = 0; arrayId < outPd->GetNumberOfArrays(); arrayId++)
+ {
+ if (auto array = vtkDataArray::SafeDownCast(outPd->GetAbstractArray(arrayId)))
+ {
+ if (!Dispatcher::Execute(array, divWorker, nbFacePts))
+ {
+ // Fallback for vtkDataArray subtypes only
+ divWorker(array, nbFacePts);
+ }
+ }
+ }
+
+ numberOfNewCells += nbFacePts;
+ globalFace += nbFacePts + 1; // Go to next face
+ }
+
+ // Compute polyhedron barycenter from faces barycenters
+ std::array<double, 3> polyBarycenter = { 0.0, 0.0, 0.0 };
+ for (vtkIdType ptId = 0; ptId < barycenters->GetNumberOfPoints(); ptId++)
+ {
+ std::array<double, 3> pt = { 0.0, 0.0, 0.0 };
+ barycenters->GetPoint(ptId, pt.data());
+ polyBarycenter[0] += pt[0];
+ polyBarycenter[1] += pt[1];
+ polyBarycenter[2] += pt[2];
+ }
+ for (int i = 0; i < 3; i++)
+ {
+ polyBarycenter[i] /= barycenters->GetNumberOfPoints();
+ }
+
+ // Compute polyhedron barycenter point data
+ for (vtkIdType arrayId = 0; arrayId < outPd->GetNumberOfArrays(); arrayId++)
+ {
+ vtkAbstractArray* abstArray = outPd->GetAbstractArray(arrayId);
+ vtkDataArray* array = vtkDataArray::SafeDownCast(abstArray);
+ if (array)
+ {
+ if (!Dispatcher::Execute(array, initWorker))
+ {
+ initWorker(array); // Fallback for vtkDataArray subtypes
+ }
+ }
+ else if (auto stringArray = vtkStringArray::SafeDownCast(abstArray))
+ {
+ initWorker(stringArray); // Fallback for vtkStringArray
+ }
+ else
+ {
+ vtkWarningWithObjectMacro(nullptr,
+ "" << abstArray->GetName()
+ << ": array type is not supported. Values on new points will be undefined.");
+ initWorker(abstArray); // Fallback for all other types
+ }
+
+ for (vtkIdType pointId = polyhedron->GetNumberOfPoints();
+ pointId < polyhedron->GetNumberOfPoints() + barycenters->GetNumberOfPoints(); pointId++)
+ {
+ if (array)
+ {
+ if (!Dispatcher2::Execute(array, array, accuWorker, pointId))
+ {
+ // Fallback for vtkDataArray subtypes only
+ accuWorker(array, array, pointId);
+ }
+ }
+ }
+ if (array)
+ {
+ if (!Dispatcher::Execute(array, divWorker, barycenters->GetNumberOfPoints()))
+ {
+ // Fallback for vtkDataArray subtypes only
+ divWorker(array, barycenters->GetNumberOfPoints());
+ }
+ }
+ }
+
+ /////////////////// Construct output UG ///////////////////
+ // Here we construct the output UG (geometry and topology).
+ // For each concave face, we generate a tetrahedron for each face edge (2 edge points + face
+ // barycenter + cell barycenter) We also fill the cell data here.
+
+ // Copy the original points to the output UG
+ vtkNew<vtkPoints> outputPoints;
+ outputGrid->SetPoints(outputPoints);
+ outputPoints->DeepCopy(polyhedron->GetPoints());
+
+ // Add the new points (barycenters) to the output UG
+ outputPoints->Resize(polyhedron->GetNumberOfPoints() + barycenters->GetNumberOfPoints() + 1);
+ for (vtkIdType newPtId = 0; newPtId < barycenters->GetNumberOfPoints(); newPtId++)
+ {
+ outputPoints->InsertNextPoint(barycenters->GetPoint(newPtId));
+ }
+ outputPoints->InsertNextPoint(polyBarycenter.data());
+
+ // Prepare output cell data
+ vtkCellData* outCd = outputGrid->GetCellData();
+ outCd->CopyAllocate(inCd, numberOfNewCells);
+
+ vtkIdType barycenterId = polyhedron->GetNumberOfPoints();
+ vtkIdType polyBarycenterId = polyhedron->GetNumberOfPoints() + barycenters->GetNumberOfPoints();
+
+ // Insert to UG a new tetra. Tetra points:
+ // ptId1, ptId2 (forming one face edge), face barycenter, polyhedron barycenter
+ auto insertTetra = [&](vtkCell* face, vtkIdType ptId1, vtkIdType ptId2) {
+ vtkIdType ptIds[4] = { 0 };
+ ptIds[0] = face->GetPointId(ptId1);
+ ptIds[1] = barycenterId;
+ ptIds[2] = face->GetPointId(ptId2);
+ ptIds[3] = polyBarycenterId;
+ vtkIdType newCellId = outputGrid->InsertNextCell(VTK_TETRA, 4, ptIds);
+ outCd->CopyData(inCd, cellId, newCellId);
+ };
+
+ // Add cells to output UG. Each new cell will contain the same data that the current polyhedron.
+ // This can be potentially improved by finding a way to insert index at the same time we
+ // create the barycenters, in order to avoid re-iterating over all the faces again
+ for (vtkIdType faceId = 0; faceId < polyhedron->GetNumberOfFaces(); ++faceId)
+ {
+ vtkCell* face = polyhedron->GetFace(faceId);
+
+ for (vtkIdType ptId = 0; ptId < face->GetNumberOfPoints() - 1; ptId++)
+ {
+ insertTetra(face, ptId, ptId + 1);
+ }
+
+ insertTetra(face, face->GetNumberOfPoints() - 1, 0);
+ barycenterId++;
+ }
+
+ return outputGrid;
+}
+
+VTK_ABI_NAMESPACE_END
diff --git a/Common/DataModel/vtkPolyhedronUtilities.h b/Common/DataModel/vtkPolyhedronUtilities.h
new file mode 100644
index 00000000000..71ed667d071
--- /dev/null
+++ b/Common/DataModel/vtkPolyhedronUtilities.h
@@ -0,0 +1,65 @@
+/*=========================================================================
+
+ Program: Visualization Toolkit
+ Module: vtkPolyhedronUtilities.h
+
+ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
+ All rights reserved.
+ See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notice for more information.
+
+=========================================================================*/
+/**
+ * @class vtkPolyhedronUtilities
+ * @brief vtkPolyhedron utilities
+ *
+ * This class contains specific methods used to process vtkPolyhedron.
+ * These methods are intended to improve filters behavior over bad-shaped or degenerated
+ * polyhedrons (for example, non-planar ones).
+ *
+ * @sa
+ * vtkPolyhedron
+ */
+
+#ifndef vtkPolyhedronUtilities_h
+#define vtkPolyhedronUtilities_h
+
+#include "vtkCommonDataModelModule.h" // For export macro
+#include "vtkObject.h"
+#include "vtkSetGet.h" // For vtkTypeMacro
+#include "vtkSmartPointer.h" // For vtkSmartPointer
+#include "vtkType.h" // For vtkIdType
+
+VTK_ABI_NAMESPACE_BEGIN
+class vtkCellData;
+class vtkPointData;
+class vtkPolyhedron;
+class vtkUnstructuredGrid;
+
+class VTKCOMMONDATAMODEL_EXPORT vtkPolyhedronUtilities
+{
+public:
+ /**
+ * Decompose the input polyhedron into tetrahedrons.
+ * This method will generate new points on each faces (faces barycenters) and another that is the
+ * barycenter of the cell. These new points are used to create the tetrahedrons and will lead to
+ * better result when applying filters (for example contours) on the output if the input
+ * polyhedron contains concave faces. The user can give point data and cell data to be passed
+ * through the decomposition. The point data on the new points (barycenters) correspond to the
+ * mean value of the respective data on the faces points. The point data on the barycenter of the
+ * cell correspond to the mean value of the respective data on all points. The cell data at given
+ * cellId will simply be copied to each output tetrahedron.
+ */
+ static vtkSmartPointer<vtkUnstructuredGrid> Decompose(
+ vtkPolyhedron* polyhedron, vtkPointData* inPd, vtkIdType cellId, vtkCellData* inCd);
+
+private:
+ vtkPolyhedronUtilities() = default;
+ ~vtkPolyhedronUtilities() = default;
+};
+
+VTK_ABI_NAMESPACE_END
+#endif
diff --git a/Documentation/release/dev/add-polyhedron-utilities.md b/Documentation/release/dev/add-polyhedron-utilities.md
new file mode 100644
index 00000000000..ce42a6e3730
--- /dev/null
+++ b/Documentation/release/dev/add-polyhedron-utilities.md
@@ -0,0 +1,8 @@
+## Add vtkPolyhedronUtilities and Decompose method
+
+Add a new method allowing to decompose a polyhedron into tetrahedrons.
+This method will generate an unstructured grid containing more points that the input polyhedron.
+Its purpose is to improve the result of subsequent filters (e.g. contours) on polyhedrons with concave faces.
+
+The vtkPolyhedronUtilities class is added to contain this method,
+and potentially future ones used for the same kind of purpose (adding a better support of bad-shaped polyhedrons).
--
GitLab