From ada8a3025551c3ea8cc6d7272023b74e96d63f56 Mon Sep 17 00:00:00 2001
From: Astrid <astrid.beyer@etu.univ-amu.fr>
Date: Thu, 25 May 2023 17:50:55 +0200
Subject: [PATCH] wip: filter / error on paraview...
---
main.cpp | 176 +++++++++++++++++++------------------------------------
1 file changed, 59 insertions(+), 117 deletions(-)
diff --git a/main.cpp b/main.cpp
index 58c0bda..3ff1693 100644
--- a/main.cpp
+++ b/main.cpp
@@ -182,121 +182,63 @@ int main(int argc, char *argv[])
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();
*/
- /** Filter based on area **/
- double threshold = 0.5;
- vtkCellArray *triangles = source->GetPolys(); // Get triangles of source
- vtkPoints *points = source->GetPoints(); // Get vertices of source
+ /** Another filter **/
+ double threshold = 0.9;
- std::vector<double> triangleAreas(source->GetNumberOfPoints(), 0.0); // This is where we store triangle's areas
+ vtkSmartPointer<vtkPoints> filteredPoints = vtkSmartPointer<vtkPoints>::New();
+ vtkSmartPointer<vtkDoubleArray> filteredShapeIndex = vtkSmartPointer<vtkDoubleArray>::New();
+ filteredShapeIndex->SetName("Filtered_Shape_Index");
- vtkSmartPointer<vtkIdList> cellPoints = vtkSmartPointer<vtkIdList>::New();
- triangles->InitTraversal();
- while (triangles->GetNextCell(cellPoints))
+ for (vtkIdType i = 0; i < k2->GetNumberOfPoints(); ++i)
{
- 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);
+ double kmax = k1Array->GetVariantValue(i).ToDouble();
+ double kmin = k2Array->GetVariantValue(i).ToDouble();
+ double si = 0.5 - (1 / M_PI) * atan((kmax + kmin) / (kmin - kmax)); // Calcul de l'indice de forme
- // Add triangle's area to corresponding vertices
- triangleAreas[cellPoints->GetId(0)] += triangleArea;
- triangleAreas[cellPoints->GetId(1)] += triangleArea;
- triangleAreas[cellPoints->GetId(2)] += triangleArea;
+ if (si > 0 && si < 1) // Vérifiez si l'indice de forme est différent de 0 et 1
+ {
+ if (si >= threshold) // Vérifiez si l'indice de forme dépasse le seuil
+ {
+ filteredPoints->InsertNextPoint(source->GetPoint(i));
+ filteredShapeIndex->InsertNextValue(si);
+ }
}
}
- // 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]; });
+ // Créez un nouveau maillage à partir des points filtrés
+ vtkSmartPointer<vtkPolyData> filteredMesh = vtkSmartPointer<vtkPolyData>::New();
+ filteredMesh->SetPoints(filteredPoints);
+ filteredMesh->GetPointData()->AddArray(filteredShapeIndex);
- vtkSmartPointer<vtkPolyData> mergedPolyData = vtkSmartPointer<vtkPolyData>::New();
- mergedPolyData->DeepCopy(source);
- // Iterate over sorted vertices
- for (size_t i = 0; i < vertexIndices.size(); i++)
+ for (vtkIdType i = 0; i < filteredShapeIndex->GetNumberOfTuples(); ++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);
- }
- }
- }
- }
- }
+ double value = filteredShapeIndex->GetValue(i);
+ std::cout << "Shape_Index[" << i << "] = " << value << std::endl;
}
- // Now vertexIndices contains the sorted vertex indices in descending order of triangleAreas
- // You can use vertexIndices to access the vertices in the desired order
+ source->GetPointData()->AddArray(filteredShapeIndex);
+ source->GetPointData()->SetActiveScalars("Shape_Index");
- // 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();
+ // vtkNew<ttkFTRGraph> graph;
+ // graph->SetInputData(filteredMesh);
+ // graph->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "Shape_Index");
+ // graph->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->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "Shape_Index");
+ reebGraph->SetInputData(source);
reebGraph->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "Shape_Index");
reebGraph->Update();
@@ -307,24 +249,24 @@ int main(int argc, char *argv[])
ttkIcospheresFromPoints->SetNumberOfSubdivisions(1);
ttkIcospheresFromPoints->Update();
- // Filters dependencies with Tubes : applies smooth to the geometry of the Reeb Graph
- vtkNew<ttkGeometrySmoother> ttkGeometrySmoother{};
- ttkGeometrySmoother->SetInputConnection(reebGraph->GetOutputPort(1));
- ttkGeometrySmoother->Update();
+ // // Filters dependencies with Tubes : applies smooth to the geometry of the Reeb Graph
+ // vtkNew<ttkGeometrySmoother> ttkGeometrySmoother{};
+ // ttkGeometrySmoother->SetInputConnection(reebGraph->GetOutputPort(1));
+ // ttkGeometrySmoother->Update();
- vtkNew<vtkGeometryFilter> geometryFilter{};
- geometryFilter->SetInputConnection(ttkGeometrySmoother->GetOutputPort());
- geometryFilter->Update();
+ // vtkNew<vtkGeometryFilter> geometryFilter{};
+ // geometryFilter->SetInputConnection(ttkGeometrySmoother->GetOutputPort());
+ // geometryFilter->Update();
- vtkNew<vtkPolyData> polyData{};
- polyData->ShallowCopy(geometryFilter->GetOutput());
+ // vtkNew<vtkPolyData> polyData{};
+ // polyData->ShallowCopy(geometryFilter->GetOutput());
- // Create Tubes to represent edges in the Reeb Graph
- vtkNew<vtkTubeFilter> tube{};
- tube->SetInputData(polyData); // smoothed geometry of the Reeb graph
- tube->SetRadius(0.004);
- // tube->SetNumberOfSides(1);
- tube->Update();
+ // // Create Tubes to represent edges in the Reeb Graph
+ // vtkNew<vtkTubeFilter> tube{};
+ // tube->SetInputData(polyData); // smoothed geometry of the Reeb graph
+ // tube->SetRadius(0.004);
+ // // tube->SetNumberOfSides(1);
+ // tube->Update();
// Save the graph nodes with IcospheresFromPoints
vtkNew<vtkXMLPolyDataWriter> writerNodes{};
@@ -332,11 +274,11 @@ int main(int argc, char *argv[])
writerNodes->SetInputData(ttkIcospheresFromPoints->GetOutput());
writerNodes->Write();
- // Save the graph edges with Tubes
- vtkNew<vtkXMLPolyDataWriter> writerArcs{};
- writerArcs->SetFileName("ReebGraphArcs.vtp");
- writerArcs->SetInputData(tube->GetOutput());
- writerArcs->Write();
+ // // Save the graph edges with Tubes
+ // vtkNew<vtkXMLPolyDataWriter> writerArcs{};
+ // writerArcs->SetFileName("ReebGraphArcs.vtp");
+ // writerArcs->SetInputData(tube->GetOutput());
+ // writerArcs->Write();
// Save the graph coloration
vtkNew<vtkXMLPolyDataWriter> segWriter{};
--
GitLab