diff --git a/main.cpp b/main.cpp
index 06a96e4010dca779e8f042e87e64ba2b313448b4..58c0bdac670dc83400c64cb22453af54a8a16e8a 100644
--- a/main.cpp
+++ b/main.cpp
@@ -56,6 +56,11 @@
#include <vtkSortDataArray.h>
#include "operation_mesh.h"
+#include <vtkGraph.h>
+#include <vtkMutableUndirectedGraph.h>
+#include <Triangulation.h>
+#include <algorithm>
+
int main(int argc, char *argv[])
{
vtkObject::GlobalWarningDisplayOff();
@@ -141,53 +146,157 @@ int main(int argc, char *argv[])
source->GetPointData()->AddArray(shapeIndex);
source->GetPointData()->SetActiveScalars("Shape_Index");
+ /** FILTER **/
+ /** filter based on scalar value of shape index **/
/*
- double pourcent = 0.8;
-
- // Create a vtkThresholdPoints filter for each threshold value
- vtkNew<vtkThresholdPoints> threshold1;
- threshold1->SetInputData(source);
- threshold1->ThresholdBetween(0.0, 0.0 + pourcent*0.25); // Scalar values near 0
- threshold1->Update();
-
- vtkNew<vtkThresholdPoints> threshold2;
- threshold2->SetInputData(source);
- threshold2->ThresholdBetween(0.5 - pourcent*0.25, 0.5 + pourcent*0.25); // Scalar values near 1
- threshold2->Update();
-
- vtkNew<vtkThresholdPoints> threshold3;
- threshold3->SetInputData(source);
- threshold3->ThresholdBetween(1.0 - pourcent*0.25, 1.0); // Scalar values near 1
- threshold3->Update();
-
- // Create a vtkFTRGraph for each thresholded dataset
- vtkNew<ttkFTRGraph> graph1;
- graph1->SetInputData(threshold1->GetOutput());
- graph1->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "Shape_Index");
- graph1->Update();
-
- vtkNew<ttkFTRGraph> graph2;
- graph2->SetInputData(threshold2->GetOutput());
- graph2->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "Shape_Index");
- graph2->Update();
-
- vtkNew<ttkFTRGraph> graph3;
- graph3->SetInputData(threshold3->GetOutput());
- graph3->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "Shape_Index");
- graph3->Update();
-
- // Merge the thresholded datasets together into a single polydata
- vtkNew<vtkAppendPolyData> appendFilter;
- appendFilter->AddInputData(threshold1->GetOutput());
- appendFilter->AddInputData(threshold2->GetOutput());
- appendFilter->AddInputData(threshold3->GetOutput());
- appendFilter->Update();
+ double pourcent = 0.8;
+
+ // Create a vtkThresholdPoints filter for each threshold value
+ vtkNew<vtkThresholdPoints> threshold1;
+ threshold1->SetInputData(source);
+ threshold1->ThresholdBetween(0.0, 0.0 + pourcent * 0.25); // Scalar values near 0
+ threshold1->Update();
+
+ vtkNew<vtkThresholdPoints> threshold2;
+ threshold2->SetInputData(source);
+ threshold2->ThresholdBetween(0.5 - pourcent * 0.25, 0.5 + pourcent * 0.25); // Scalar values near 1
+ threshold2->Update();
+
+ vtkNew<vtkThresholdPoints> threshold3;
+ threshold3->SetInputData(source);
+ threshold3->ThresholdBetween(1.0 - pourcent * 0.25, 1.0); // Scalar values near 1
+ threshold3->Update();
+
+ // Create a vtkFTRGraph for each thresholded dataset
+ vtkNew<ttkFTRGraph> graph1;
+ graph1->SetInputData(threshold1->GetOutput());
+ graph1->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "Shape_Index");
+ graph1->Update();
+
+ vtkNew<ttkFTRGraph> graph2;
+ graph2->SetInputData(threshold2->GetOutput());
+ graph2->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "Shape_Index");
+ graph2->Update();
+
+ vtkNew<ttkFTRGraph> graph3;
+ graph3->SetInputData(threshold3->GetOutput());
+ graph3->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "Shape_Index");
+ graph3->Update();
*/
+ /** Filter based on area **/
+ double threshold = 0.5;
+ vtkCellArray *triangles = source->GetPolys(); // Get triangles of source
+ vtkPoints *points = source->GetPoints(); // Get vertices of source
+
+ std::vector<double> triangleAreas(source->GetNumberOfPoints(), 0.0); // This is where we store triangle's areas
+
+ vtkSmartPointer<vtkIdList> cellPoints = vtkSmartPointer<vtkIdList>::New();
+ triangles->InitTraversal();
+ while (triangles->GetNextCell(cellPoints))
+ {
+ if (cellPoints->GetNumberOfIds() == 3) // Process only triangles
+ {
+ std::array<double, 3> p0, p1, p2;
+
+ // Get coordinates of triangle's points
+ points->GetPoint(cellPoints->GetId(0), p0.data());
+ points->GetPoint(cellPoints->GetId(1), p1.data());
+ points->GetPoint(cellPoints->GetId(2), p2.data());
+
+ double triangleArea = operation::computeTriangleArea(p0, p1, p2);
+
+ // Add triangle's area to corresponding vertices
+ triangleAreas[cellPoints->GetId(0)] += triangleArea;
+ triangleAreas[cellPoints->GetId(1)] += triangleArea;
+ triangleAreas[cellPoints->GetId(2)] += triangleArea;
+ }
+ }
+
+ // Create a vector of indices to keep track of the original vertex indices
+ std::vector<size_t> vertexIndices(source->GetNumberOfPoints());
+ std::iota(vertexIndices.begin(), vertexIndices.end(), 0);
+
+ // Sort the vertex indices based on triangleAreas in descending order
+ std::sort(vertexIndices.begin(), vertexIndices.end(), [&](size_t a, size_t b)
+ { return triangleAreas[a] > triangleAreas[b]; });
+
+ vtkSmartPointer<vtkPolyData> mergedPolyData = vtkSmartPointer<vtkPolyData>::New();
+ mergedPolyData->DeepCopy(source);
+ // Iterate over sorted vertices
+ for (size_t i = 0; i < vertexIndices.size(); i++)
+ {
+ vtkIdType pointId = vertexIndices[i];
+
+ // Check if the vertex is extremal
+ int degree = operation::getDegree(source, pointId);
+ std::cout << "degree: " << degree << std::endl;
+ if (degree == 1)
+ {
+ // Get the weight or area of the current vertex
+ double vertexArea = triangleAreas[pointId];
+
+ std::cout << vertexArea << std::endl;
+
+ // Check if the vertex's area is below the threshold
+ if (vertexArea < threshold)
+ {
+ std::cout << "hi";
+ continue; // ignore the vertex
+ }
+
+ else
+ {
+ // Get the ID of the neighboring vertex
+ vtkIdType neighborId = operation::getNeighborId(mergedPolyData, pointId);
+
+ // Compute the centroid between the current vertex and its neighbor
+ std::array<double, 3> centroid = operation::computeCentroid(mergedPolyData, pointId, neighborId);
+
+ // Replace the two vertices with a new vertex located at the centroid
+ vtkSmartPointer<vtkPoints> points = mergedPolyData->GetPoints();
+ points->SetPoint(pointId, centroid.data()); // Update the position of the current vertex
+ points->SetPoint(neighborId, centroid.data()); // Update the position of the neighbor vertex
+
+ // Update the neighboring triangles to point to the new vertex
+ vtkSmartPointer<vtkIdList> cellIds = vtkSmartPointer<vtkIdList>::New();
+ mergedPolyData->GetPointCells(pointId, cellIds);
+
+ for (vtkIdType i = 0; i < cellIds->GetNumberOfIds(); i++)
+ {
+ vtkIdType cellId = cellIds->GetId(i);
+ vtkSmartPointer<vtkCell> cell = mergedPolyData->GetCell(cellId);
+
+ for (vtkIdType j = 0; j < cell->GetNumberOfPoints(); j++)
+ {
+ vtkIdType cellPointId = cell->GetPointId(j);
+
+ if (cellPointId == pointId || cellPointId == neighborId)
+ {
+ // Update the point IDs of the vertices in the cell to point to the new vertex
+ cell->PointIds->SetId(j, pointId);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // Now vertexIndices contains the sorted vertex indices in descending order of triangleAreas
+ // You can use vertexIndices to access the vertices in the desired order
+
+ // Merge the thresholded datasets together into a single polydata
+ vtkNew<vtkAppendPolyData> appendFilter;
+ appendFilter->AddInputData(mergedPolyData);
+ // appendFilter->AddInputData(threshold1->GetOutput());
+ // appendFilter->AddInputData(threshold2->GetOutput());
+ // appendFilter->AddInputData(threshold3->GetOutput());
+ appendFilter->Update();
+
// Compute the Reeb Graph of the input dataset
vtkNew<ttkFTRGraph> reebGraph;
- // reebGraph->SetInputData(appendFilter->GetOutput()); // ttkFTRGraph for the merged polydata
- reebGraph->SetInputData(source);
+ reebGraph->SetInputData(appendFilter->GetOutput()); // ttkFTRGraph for the merged polydata
+ // reebGraph->SetInputData(source);
reebGraph->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "Shape_Index");
reebGraph->Update();
@@ -255,4 +364,4 @@ int main(int argc, char *argv[])
segWriter->Write();*/
return EXIT_SUCCESS;
-}
+}
\ No newline at end of file
diff --git a/operation_mesh.cpp b/operation_mesh.cpp
index 242b9cb839aa09fbe1756098614b948541408ae5..23480b340049097704249d65bbfe00460b0ee54c 100644
--- a/operation_mesh.cpp
+++ b/operation_mesh.cpp
@@ -1,5 +1,19 @@
#include "operation_mesh.h"
+#include <vtkSmartPointer.h>
+#include <vtkPolyData.h>
+#include <vtkFeatureEdges.h>
+#include <vtkFillHolesFilter.h>
+#include <vtkCleanPolyData.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkActor.h>
+#include <vtkRenderer.h>
+#include <vtkRenderWindow.h>
+#include <vtkRenderWindowInteractor.h>
+
+#include <array>
+#include <cmath>
+
vtkSmartPointer<vtkPolyData> operation::buildValidTriangulation(vtkSmartPointer<vtkPolyData> polyData)
{
vtkSmartPointer<vtkPolyData> triPolyData = vtkSmartPointer<vtkPolyData>::New();
@@ -33,17 +47,6 @@ vtkSmartPointer<vtkPolyData> operation::buildValidTriangulation(vtkSmartPointer<
return triangleFilter->GetOutput();
}
-#include <vtkSmartPointer.h>
-#include <vtkPolyData.h>
-#include <vtkFeatureEdges.h>
-#include <vtkFillHolesFilter.h>
-#include <vtkCleanPolyData.h>
-#include <vtkPolyDataMapper.h>
-#include <vtkActor.h>
-#include <vtkRenderer.h>
-#include <vtkRenderWindow.h>
-#include <vtkRenderWindowInteractor.h>
-
vtkSmartPointer<vtkPolyData> operation::closeMesh(vtkSmartPointer<vtkPolyData> polyData)
{
// Check if the input mesh is open
@@ -58,7 +61,8 @@ vtkSmartPointer<vtkPolyData> operation::closeMesh(vtkSmartPointer<vtkPolyData> p
vtkSmartPointer<vtkPolyData> output = featureEdges->GetOutput();
int numOpenEdges = output->GetNumberOfCells();
- if (numOpenEdges != 0) {
+ if (numOpenEdges != 0)
+ {
std::cout << "There are " << numOpenEdges << " open edges.\n";
// Display the output mesh
@@ -75,7 +79,7 @@ vtkSmartPointer<vtkPolyData> operation::closeMesh(vtkSmartPointer<vtkPolyData> p
renderWindow->Render();
renderWindowInteractor->Start();
- //TODO: not working for the moment...
+ // TODO: not working for the moment...
/*
// Get point coordinates
vtkPoints* points = output->GetPoints();
@@ -131,8 +135,117 @@ vtkSmartPointer<vtkPolyData> operation::closeMesh(vtkSmartPointer<vtkPolyData> p
}
}
return output;*/
- } else {
+ }
+ else
+ {
std::cout << "There are no open edges.\n";
- }
+ }
return polyData;
}
+
+int operation::getDegree(vtkSmartPointer<vtkPolyData> source, vtkIdType pointId)
+{
+ vtkSmartPointer<vtkIdList> cellIds = vtkSmartPointer<vtkIdList>::New();
+ source->GetPointCells(pointId, cellIds);
+
+ int numNeighbors = 0;
+ for (vtkIdType i = 0; i < cellIds->GetNumberOfIds(); i++)
+ {
+ vtkSmartPointer<vtkCell> cell = source->GetCell(cellIds->GetId(i));
+ for (vtkIdType j = 0; j < cell->GetNumberOfPoints(); j++)
+ {
+ vtkIdType neighborId = cell->GetPointId(j);
+ if (neighborId != pointId)
+ {
+ numNeighbors++;
+ }
+ }
+ }
+ return numNeighbors;
+}
+
+void operation::getPointCells(vtkPolyData *polyData, vtkIdType pointId, vtkIdList *cellIds)
+{
+ // Clear the vtkIdList
+ cellIds->Reset();
+
+ // Iterate over all cells in the polyData
+ vtkSmartPointer<vtkIdList> pointIds = vtkSmartPointer<vtkIdList>::New();
+ vtkIdType ncells = polyData->GetNumberOfCells();
+ for (vtkIdType i = 0; i < ncells; i++)
+ {
+ polyData->GetCellPoints(i, pointIds);
+
+ // Check if the pointId is one of the cell's points
+ vtkIdType npoints = pointIds->GetNumberOfIds();
+ for (vtkIdType j = 0; j < npoints; j++)
+ {
+ if (pointIds->GetId(j) == pointId)
+ {
+ // Add the cell to the list of cellIds
+ cellIds->InsertUniqueId(i);
+ break;
+ }
+ }
+ }
+}
+
+double operation::computeTriangleArea(const std::array<double, 3> &p0, const std::array<double, 3> &p1, const std::array<double, 3> &p2)
+{
+ // Calculez les longueurs des côtés du triangle
+ double a = std::sqrt(std::pow(p1[0] - p0[0], 2) + std::pow(p1[1] - p0[1], 2) + std::pow(p1[2] - p0[2], 2));
+ double b = std::sqrt(std::pow(p2[0] - p1[0], 2) + std::pow(p2[1] - p1[1], 2) + std::pow(p2[2] - p1[2], 2));
+ double c = std::sqrt(std::pow(p0[0] - p2[0], 2) + std::pow(p0[1] - p2[1], 2) + std::pow(p0[2] - p2[2], 2));
+
+ // Calculez l'aire du triangle en utilisant la formule de la demi-somme des longueurs des côtés
+ double s = (a + b + c) / 2.0;
+ double area = std::sqrt(s * (s - a) * (s - b) * (s - c));
+
+ return area;
+}
+
+vtkIdType operation::getNeighborId(vtkSmartPointer<vtkPolyData> source, vtkIdType pointId)
+{
+ vtkSmartPointer<vtkIdList> cellIds = vtkSmartPointer<vtkIdList>::New();
+ source->GetPointCells(pointId, cellIds);
+
+ // Browse cell's neighbors
+ for (vtkIdType i = 0; i < cellIds->GetNumberOfIds(); i++)
+ {
+ vtkIdType cellId = cellIds->GetId(i);
+ vtkCell *cell = source->GetCell(cellId);
+
+ // Browse points of each cell
+ for (vtkIdType j = 0; j < cell->GetNumberOfPoints(); j++)
+ {
+ vtkIdType neighborId = cell->GetPointId(j);
+ // Check if the neighbor is not the vertex itself
+ if (neighborId != pointId)
+ {
+ return neighborId;
+ }
+ }
+ }
+ // -1 if can't find any unique neighbor
+ return -1;
+}
+
+std::array<double, 3> operation::computeCentroid(vtkSmartPointer<vtkPolyData> source, vtkIdType pointId1, vtkIdType pointId2)
+{
+ std::array<double, 3> centroid = {0.0, 0.0, 0.0};
+
+ vtkPoints *points = source->GetPoints();
+
+ // get coords of the two points
+ double p1[3], p2[3];
+ points->GetPoint(pointId1, p1);
+ points->GetPoint(pointId2, p2);
+
+ // Calculate the coordinates of the centroid by averaging the coordinates of the two points
+ for (int i = 0; i < 3; i++)
+ {
+ centroid[i] = (p1[i] + p2[i]) / 2.0;
+ }
+
+ return centroid;
+}
diff --git a/operation_mesh.h b/operation_mesh.h
index 0dbec1b848850d114b84d20867fd7d97725d4224..883a4705e5cc3a844b721fedb6f4b58efa89338a 100644
--- a/operation_mesh.h
+++ b/operation_mesh.h
@@ -3,7 +3,7 @@
*
* @author Astrid BEYER
* @date 2023
- * @version 1.0
+ * @version 1.2
*
*/
@@ -12,7 +12,6 @@
#include <vtkTriangleFilter.h>
-
namespace operation
{
/**
@@ -25,9 +24,32 @@ namespace operation
* @brief closeMesh : close the surface of the mesh, if necessary
* @param polyData
* @returns vtkSmartPointer<vtkPolyData>
- */
+ */
vtkSmartPointer<vtkPolyData> closeMesh(vtkSmartPointer<vtkPolyData> polyData);
-
+
+ /**
+ * @brief getDegree : returns number of neighbors a vertex belongs to
+ * @param source, point
+ * @returns int
+ */
+ int getDegree(vtkSmartPointer<vtkPolyData> source, vtkIdType point);
+
+ /**
+ * @brief getPointCells : get the identifiers of the vertices of a given cell
+ * @param polyData, pointId, cellIds
+ * @returns void
+ */
+ void getPointCells(vtkPolyData *polyData, vtkIdType pointId, vtkIdList *cellIds);
+
+ /**
+ * @brief computeTriangleArea : using formula for the half-sum of side lengths to calculate the area of a triangle
+ * @param p0, p1, p2
+ * @returns double
+ */
+ double computeTriangleArea(const std::array<double, 3> &p0, const std::array<double, 3> &p1, const std::array<double, 3> &p2);
+
+ vtkIdType getNeighborId(vtkSmartPointer<vtkPolyData> source, vtkIdType pointId);
+ std::array<double, 3> computeCentroid(vtkSmartPointer<vtkPolyData> source, vtkIdType pointId1, vtkIdType pointId2);
}
#endif