VTK – Remesh Model By Flip Operation

A edge pt1-pt2 is used by two cells as in the following image, we need to flip the cells if the point opt2 is in the circumscribed sphere of points set (pt1, pt2, opt1).

We have to find all edges and judge whether the point is in the circumscribed sphere of another three points, the target edge is used by two cells, it has only one neighbor cell in the other words. Finally, flip the cells as the above image.

I use multimap to store all edges, the data structure maybe not the best one, container vector with the custom C++ structure can also work and be faster.

The original model looks like this scene.

#include <iostream>
#include <vtkSmartPointer.h>
#include <vtkSphereSource.h>
#include <vtkActor.h>
#include <vtkConeSource.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkPolyDataMapper.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkPolyData.h>
#include <vtkSTLWriter.h>
#include <vtkImplicitPolyDataDistance.h>
#include <vtkProperty.h>
#include <vtkSelectEnclosedPoints.h>
#include <vtkTriangleFilter.h>
#include <vtkSTLReader.h>
#include <vtkPointData.h>
#include <vtkColorTransferFunction.h>
#include <vtkPolyDataConnectivityFilter.h>
#include <vtkLinearSubdivisionFilter.h>
#include <vtkButterflySubdivisionFilter.h>
#include <vtkNamedColors.h>
#include <vtkClipPolyData.h>

#include <map>
#include "point.hpp"

#define vtkSPtr vtkSmartPointer
#define vtkSPtrNew(Var, Type) vtkSPtr<Type> Var = vtkSPtr<Type>::New();

using namespace std;

class EdgeList
{
public:
    void Init(){ m_List.clear();  };
    void AddEdge(vtkIdType ptId1, vtkIdType ptId2){
        if( ptId1 > ptId2 ){
            std::swap( ptId1, ptId2 );
        }
        std::multimap<int, int>::iterator itr = m_List.begin();
        for( const auto&[key, value]: m_List )
        {
            if( key == ptId1 && value == ptId2 )
            {
                return ;
            }
            itr++;
        }

        m_List.insert( itr, std::make_pair( ptId1, ptId2 ) );
    }
    void RemoveEdge(vtkIdType ptId1, vtkIdType ptId2){
        if( ptId1 > ptId2 ){
            std::swap( ptId1, ptId2 );
        }
        std::multimap<int, int>::iterator itr = m_List.begin();
        for( const auto&[key, value]: m_List )
        {
            if( key == ptId1 && value == ptId2 )
            {
                break;
            }
            itr++;
        }
        if( itr != m_List.end() )
        {
            m_List.erase( itr );
        }
    }
    std::multimap<int, int> Getm_List(){ return m_List; };
protected:
    std::multimap<int, int> m_List;
};

bool IsInsideCircumcircle
    (
    double a[3], 
    double b[3], 
    double c[3], 
    double point[3]
    )
{
    double center[3];
    vtkMath::Solve3PointCircle(a, b, c, center);
    double xCenterVec[3], aCenterVec[3];
    vtkMath::Subtract(point, center, xCenterVec);
    vtkMath::Subtract(a, center, aCenterVec);
    double lengthXCenter = vtkMath::Norm(xCenterVec);
    double lengthACenter = vtkMath::Norm(aCenterVec);
    return (lengthXCenter < lengthACenter);
}

void FindEdgeCellsAndThirdPts
    (
    vtkSPtr<vtkPolyData> data,
    int pointId1,
    int pointId2,
    vtkIdType &edgeCellId1,
    vtkIdType &edgeCellId2,
    vtkIdType &thirdPointId1,
    vtkIdType &thirdPointId2
    )
{
    if (data->IsEdge(pointId1, pointId2) == false)
        return;

    edgeCellId1 = -1;
    edgeCellId2 = -1;
    thirdPointId1 = -1;
    thirdPointId2 = -1;

    vtkSPtrNew( cellPtIds1, vtkIdList );
    vtkSPtrNew( cellPtIds2, vtkIdList );

    data->GetPointCells(pointId1, cellPtIds1);
    data->GetPointCells(pointId2, cellPtIds2);

    for (int i = 0; i < cellPtIds1->GetNumberOfIds(); i++)
    {
        int cellId = cellPtIds1->GetId(i);
        if (data->GetCellType(cellId) != VTK_TRIANGLE)
            continue;
        for (int j = 0; j < cellPtIds2->GetNumberOfIds(); j++)
        {
            int cellId1 = cellPtIds2->GetId(j);
            if (data->GetCellType(cellId1) != VTK_TRIANGLE)
                continue;
            if (cellId == cellId1)
            {
                vtkIdType nPoints;
                const vtkIdType *pts;
                data->GetCellPoints(cellId, nPoints, pts);
                int pointId = -1;
                for (int k = 0; k < nPoints; k++)
                {
                    if (pts[k] != pointId1&&pts[k] != pointId2)
                    {
                        pointId = pts[k];
                        break;
                    }
                }
                if (thirdPointId1 == -1)
                {
                    thirdPointId1 = pointId;
                    edgeCellId1 = cellId;
                    break;
                }
                else if (thirdPointId2 == -1)
                {
                    thirdPointId2 = pointId;
                    edgeCellId2 = cellId;
                    break;
                }
            }
        }
    }
}

void Flip
    (
    vtkSPtr<vtkPolyData> data, 
    int pointId1, 
    int pointId2, 
    vtkIdType edgeCellId1, 
    vtkIdType edgeCellId2, 
    vtkIdType thirdPointId, 
    vtkIdType thirdPointId1
    )
{
    vtkSPtrNew( cellPtIds1, vtkIdList );
    vtkSPtrNew( cellPtIds2, vtkIdList );

    data->GetCellPoints(edgeCellId1, cellPtIds1);
    data->GetCellPoints(edgeCellId2, cellPtIds2);

    vtkIdType tri[3];
    for (int i = 0; i < 3; i++)
    {
        if (cellPtIds1->GetId(i) == pointId2)
            tri[i] = thirdPointId1;
        else
            tri[i] = cellPtIds1->GetId(i);
    }
    data->RemoveReferenceToCell(pointId2, edgeCellId1);
    data->ReplaceCell(edgeCellId1, 3, tri);
    data->ResizeCellList(thirdPointId1, 1);
    data->AddReferenceToCell(thirdPointId1, edgeCellId1);

    vtkIdType tri1[3];
    for (int i = 0; i < 3; i++)
    {
        if (cellPtIds2->GetId(i) == pointId1)
            tri1[i] = thirdPointId;
        else
            tri1[i] = cellPtIds2->GetId(i);
    }
    data->RemoveReferenceToCell(pointId1, edgeCellId2);
    data->ReplaceCell(edgeCellId2, 3, tri1);
    data->ResizeCellList(thirdPointId, 1);
    data->AddReferenceToCell(thirdPointId, edgeCellId2);
}


int main()
{
    vtkSPtrNew( points, vtkPoints );
    points->InsertNextPoint( 0, 0, -2 );
    points->InsertNextPoint( 2, 0, 2 );
    points->InsertNextPoint( -2, 0, 2 );
    points->InsertNextPoint( 0, 2, 0 );
    vtkSPtrNew( cells, vtkCellArray );
    vtkIdType cell0[3] = { 0, 1, 2 };
    cells->InsertNextCell( 3, cell0 );
    vtkIdType cell1[3] = { 3, 2, 1 };
    cells->InsertNextCell( 3, cell1 );
    vtkSPtrNew( polyData, vtkPolyData );
    polyData->SetPoints( points );
    polyData->SetPolys( cells );
    polyData->BuildCells();
    polyData->BuildLinks();
    polyData->Modified();

    EdgeList edges;
    // ================== get edge list ===================
    {
        edges.Init();

        vtkSPtrNew(neighborCells, vtkIdList);
        cout << "polyData->GetNumberOfCells(): " << polyData->GetNumberOfCells() << endl;
        for (vtkIdType cellId = 0; cellId < polyData->GetNumberOfCells(); cellId++)
        {
            vtkIdType nPoints = 0;
            const vtkIdType *pts;
            polyData->GetCellPoints(cellId, nPoints, pts);
            if (nPoints != 3)    //only process triangles
                continue;
            for (vtkIdType j = 0; j < 3; j++)
            {
                vtkIdType edgePointId1 = pts[j];
                vtkIdType edgePointId2 = pts[(j + 1) % 3];
                polyData->GetCellEdgeNeighbors(cellId, edgePointId1, edgePointId2, neighborCells);
                vtkIdType numberOfNgbCells = neighborCells->GetNumberOfIds();
                if (numberOfNgbCells != 1)
                {
                    cout << "continue, numberOfNgbCells: " << numberOfNgbCells << endl;
                    continue;
                }
                vtkIdType ngbCellId = neighborCells->GetId(0);
                if (polyData->GetCellType(ngbCellId) != VTK_TRIANGLE)
                {
                    cout << "continue, GetCellType: " << polyData->GetCellType(ngbCellId) << endl;
                    continue;
                }
                edges.AddEdge(edgePointId1, edgePointId2);
            }
        }
    }
    // ================== finish: ===================

    // ================== find the edge and flip cells ================
    auto edgesList = edges.Getm_List();
    std::multimap<int, int>::iterator itr = edgesList.begin();
    while( itr != edgesList.end() )
    {
        vtkIdType edgeCellId[2] = {-1, -1}, thirdPts[2] = {-1, -1};
        vtkIdType edgePtIds[2] = { itr->first, itr->second };
        FindEdgeCellsAndThirdPts( polyData, edgePtIds[0], edgePtIds[1],
                                  edgeCellId[0], edgeCellId[1], thirdPts[0], thirdPts[1] );
        if( -1 == thirdPts[0] || -1 == thirdPts[1] ) continue;
        if( polyData->IsEdge( thirdPts[0], thirdPts[1] ) ) continue;

        Point pt[4];
        polyData->GetPoint( edgePtIds[0], pt[0].point );
        polyData->GetPoint( edgePtIds[1], pt[1].point );
        polyData->GetPoint( thirdPts[0], pt[2].point );
        polyData->GetPoint( thirdPts[1], pt[3].point );

        if( IsInsideCircumcircle( pt[0].point, pt[1].point, pt[2].point, pt[3].point ) )
        {
            Flip( polyData, edgePtIds[0], edgePtIds[1],
                    edgeCellId[0], edgeCellId[1], thirdPts[0], thirdPts[1] );
        }
        itr++;
    }
    // ================== finish: ===================

    vtkSPtrNew( writer, vtkSTLWriter );
    writer->SetInputData( polyData );
    writer->SetFileName( "/Users/weiyang/Desktop/flipped.stl" );
    writer->Write();

    vtkSPtrNew( mapper0, vtkPolyDataMapper );
    mapper0->SetInputData( polyData );
    vtkSPtrNew( actor0, vtkActor );
    actor0->SetMapper( mapper0 );

    vtkSPtrNew( renderer, vtkRenderer );
    renderer->AddActor( actor0 );

    vtkSPtrNew( renderWindow, vtkRenderWindow );
    renderWindow->AddRenderer( renderer );

    vtkSPtrNew( renderWindowInteractor, vtkRenderWindowInteractor );
    renderWindowInteractor->SetRenderWindow( renderWindow );

    renderer->ResetCamera();
    renderWindow->Render();

    renderWindowInteractor->Start();
    return 0;
}

We will get the new remeshed model.