Skip to content
Snippets Groups Projects
Commit 27562308 authored by Astrid Beyer's avatar Astrid Beyer
Browse files

wip: filter

parent 04ea05a9
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
main.cpp
+ 152
− 43
View file @ 27562308
......@@ -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
operation_mesh.cpp
+ 128
− 15
View file @ 27562308
#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;
}
operation_mesh.h
+ 26
− 4
View file @ 27562308
......@@ -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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment