operation_mesh.cpp

#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();

    // if mesh is composed of quads, builds a valid triangulation
    if (polyData->GetMaxCellSize() == 4)
    {
        std::cout << "Mesh is composed of quads\n";
        // we're using a vtkQuadToTriangle filter to convert quads in triangles
        vtkSmartPointer<vtkTriangleFilter> quadToTri = vtkSmartPointer<vtkTriangleFilter>::New();
        quadToTri->SetInputData(polyData);
        quadToTri->Update();

        triPolyData = quadToTri->GetOutput();
    }
    else if (polyData->GetMaxCellSize() == 3) // if mesh is already triangulated, return origin's vtkPolyData
    {
        triPolyData = polyData;
    }
    else
    {
        std::cout << "This mesh is not handled by this program, allows only triangulated and quadrilated meshes.\n";
        std::cout << "Maximum cell size for this mesh: " << polyData->GetMaxCellSize() << ".\n";
        exit(0);
    }
    // adding a vtkTriangleFilter to ensure that the output is a well triangulated mesh
    vtkSmartPointer<vtkTriangleFilter> triangleFilter = vtkSmartPointer<vtkTriangleFilter>::New();
    triangleFilter->SetInputData(triPolyData);
    triangleFilter->Update();

    return triangleFilter->GetOutput();
}

vtkSmartPointer<vtkPolyData> operation::closeMesh(vtkSmartPointer<vtkPolyData> polyData)
{
    // Check if the input mesh is open
    vtkSmartPointer<vtkFeatureEdges> featureEdges = vtkSmartPointer<vtkFeatureEdges>::New();
    featureEdges->SetInputData(polyData);
    // FeatureEdges configuration
    featureEdges->BoundaryEdgesOn();     // Search for boundary edges (only one neighbor)
    featureEdges->NonManifoldEdgesOff(); // Search for non-manifold edges
    featureEdges->FeatureEdgesOff();     // Deactivate search for special edges (projecting, concave, etc)
    featureEdges->Update();

    vtkSmartPointer<vtkPolyData> output = featureEdges->GetOutput();
    int numOpenEdges = output->GetNumberOfCells();

    if (numOpenEdges != 0)
    {
        std::cout << "There are " << numOpenEdges << " open edges.\n";

        // Display the output mesh
        vtkSmartPointer<vtkPolyDataMapper> mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
        mapper->SetInputData(featureEdges->GetOutput());
        vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
        actor->SetMapper(mapper);
        vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
        renderer->AddActor(actor);
        vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
        renderWindow->AddRenderer(renderer);
        vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor = vtkSmartPointer<vtkRenderWindowInteractor>::New();
        renderWindowInteractor->SetRenderWindow(renderWindow);
        renderWindow->Render();
        renderWindowInteractor->Start();

        // TODO: not working for the moment...
        /*
        // Get point coordinates
        vtkPoints* points = output->GetPoints();
        // Loop through open edges and create cells to fill gaps
        for (vtkIdType i = 0; i < output->GetNumberOfCells(); i++) {
            vtkCell* cell = output->GetCell(i);
            if (cell->GetCellType() == VTK_LINE) { // Check if open edge is a line segment
                vtkIdType id0 = cell->GetPointId(0);
                vtkIdType id1 = cell->GetPointId(1);

                vtkIdType newIds[3];
                newIds[0] = id0;
                newIds[1] = id1;

                // Find third point to create triangle
                double point0[3];
                double point1[3];
                double vec[3];
                points->GetPoint(id0, point0);
                points->GetPoint(id1, point1);
                for (vtkIdType j = 0; j < output->GetNumberOfCells(); j++) {
                    if (j == i) {
                        continue;
                    }
                    vtkCell* otherCell = output->GetCell(j);
                    if (otherCell->GetCellType() != VTK_TRIANGLE) {
                        continue;
                    }
                    vtkIdType otherIds[3];
                    otherIds[0] = otherCell->GetPointId(0);
                    otherIds[1] = otherCell->GetPointId(1);
                    otherIds[2] = otherCell->GetPointId(2);
                    if (otherIds[0] == id0 || otherIds[1] == id0 || otherIds[2] == id0) {
                        if (otherIds[0] == id1 || otherIds[1] == id1 || otherIds[2] == id1) {
                            // Found a triangle sharing these two points
                            for (vtkIdType k = 0; k < 3; k++) {
                                if (otherIds[k] != id0 && otherIds[k] != id1) {
                                    points->GetPoint(otherIds[k], vec);
                                    break;
                                }
                            }
                            // Add third point and create new cell
                            newIds[2] = points->InsertNextPoint(vec);
                            vtkSmartPointer<vtkTriangle> newCell = vtkSmartPointer<vtkTriangle>::New();
                            newCell->GetPointIds()->SetId(0, newIds[0]);
                            newCell->GetPointIds()->SetId(1, newIds[1]);
                            newCell->GetPointIds()->SetId(2, newIds[2]);
                            polyData->GetPolys()->InsertNextCell(newCell);
                            break;
                        }
                    }
                }
            }
        }
        return output;*/
    }
    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;
}