VTK – Increase Density For Model By Add Edge Points

Here is a way to remesh 3D model, for every edge we can add a few points on it and create new small cells. The old big triangle will be removed from the model.
You can read other relative posts about remesh model, VTK – Remesh Model By Flip Operation and VTK – Remesh Model By Add Center Point.

#include <iostream>
#include <vtkSmartPointer.h>
#include <vtkSphereSource.h>
#include <vtkActor.h>
#include <vtkConeSource.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkPolyDataMapper.h>
#include <vtkXMLPolyDataWriter.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkPolyData.h>
#include <vtkSTLWriter.h>
#include <vtkPlane.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;

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;
                }
            }
        }
    }
}

bool Split(vtkSPtr<vtkPolyData> data, int pointId1, int pointId2, double edgeLenMax)
{
    if (data->IsEdge(pointId1, pointId2) == false)
        return false;
    Point edgePointPos, edgePointPos1;
    data->GetPoint(pointId1, edgePointPos.point);
    data->GetPoint(pointId2, edgePointPos1.point);

    Point vec = edgePointPos1-edgePointPos;
    double length = vec.Length();
    if ( length < edgeLenMax )
        return false;

    vtkIdType thirdPointId1 = -1, thirdPointId2 = -1;
    vtkIdType edgeCellId1 = -1, edgeCellId2 = -1;
    FindEdgeCellsAndThirdPts(data, pointId1, pointId2, edgeCellId1, edgeCellId2,
                             thirdPointId1, thirdPointId2);
    if (thirdPointId1 == -1)
        return false;
    if (data->IsEdge(thirdPointId1, thirdPointId2))
        return false;
    if (data->GetCellType(edgeCellId1) != VTK_TRIANGLE)
        return false;
    if (edgeCellId2 != -1 && data->GetCellType(edgeCellId2) != VTK_TRIANGLE)
        return false;

    int numberOfPointToAdd = length / edgeLenMax;
    Point moveVecEverySeg = vec*1.0 / (numberOfPointToAdd + 1);
    std::vector<vtkIdType> pointIdsOnEdge;
    pointIdsOnEdge.push_back(pointId1);
    for (int i = 0; i < numberOfPointToAdd; i++)
    {
        Point addedPointPos = edgePointPos + (i+1)*moveVecEverySeg;
        int addedPointId;
        if (thirdPointId2 != -1)
            addedPointId = data->InsertNextLinkedPoint(addedPointPos.point, 4);
        else
            addedPointId = data->InsertNextLinkedPoint(addedPointPos.point, 2);
        pointIdsOnEdge.push_back(addedPointId);
    }
    pointIdsOnEdge.push_back(pointId2);

    int edgeOrder = 0;          //0: m_PointId1 to m_PointId2;1:m_PointId2 to m_PointId1
    if (edgeCellId1 != -1)
    {
        vtkIdType nPoints;
        const vtkIdType *pts;
        data->GetCellPoints(edgeCellId1, nPoints, pts);
        for (int i = 0; i < nPoints; i++)
        {
            if (pts[i] != pointId1)
                continue;
            if (pts[(i + 2) % 3] == pointId2)
                edgeOrder = 1;
        }
    }


    int edgeOrder1 = 0;         // 0: PointId1 to PointId2
                                // 1: PointId2 to PointId1
    if (edgeCellId2 != -1)
    {
        vtkIdType nPoints;
        const vtkIdType *pts;
        data->GetCellPoints(edgeCellId2, nPoints, pts);
        for (int i = 0; i < nPoints; i++)
        {
            if (pts[i] != pointId1)
                continue;
            if (pts[(i + 2) % 3] == pointId2)
                edgeOrder1 = 1;
        }
    }

    for (size_t i = 1; i < pointIdsOnEdge.size(); i++)
    {
        vtkIdType prePointId = pointIdsOnEdge[i - 1];
        vtkIdType curPointId = pointIdsOnEdge[i];

        //for edgeCellId1
        if (thirdPointId1 != -1 && edgeCellId1 != -1)
        {
            vtkIdType tri[3] = { prePointId, curPointId, thirdPointId1 };
            if (edgeOrder == 1)     //switch order
            {
                tri[0] = curPointId;
                tri[1] = prePointId;
            }

            if (i == 1)
            {
                data->RemoveReferenceToCell(pointId2, edgeCellId1);
                data->ReplaceCell(edgeCellId1, 3, tri);
                data->ResizeCellList(curPointId, 1);
                data->AddReferenceToCell(curPointId, edgeCellId1);
            }
            else
            {
                data->InsertNextLinkedCell(VTK_TRIANGLE, 3, tri);
            }
        }

        //for edgeCellId2
        if (thirdPointId2 != -1 && edgeCellId2 != -1)
        {
            vtkIdType tri[3] = { prePointId, curPointId, thirdPointId2 };
            if (edgeOrder1 == 1)        //switch order
            {
                tri[0] = curPointId;
                tri[1] = prePointId;
            }

            if (i == 1)
            {
                data->RemoveReferenceToCell(pointId2, edgeCellId2);
                data->ReplaceCell(edgeCellId2, 3, tri);
                data->ResizeCellList(curPointId, 1);
                data->AddReferenceToCell(curPointId, edgeCellId2);
            }
            else
            {
                data->InsertNextLinkedCell(VTK_TRIANGLE, 3, tri);
            }
        }
    }
}

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

    vtkSPtrNew(neighborCells, vtkIdList);
    int cellCount = polyData->GetNumberOfCells();
    cout << "polyData->GetNumberOfCells(): " << polyData->GetNumberOfCells() << endl;
    for (vtkIdType cellId = 0; cellId < cellCount; 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;
            }
            cout << "cellId: " << cellId << endl;
            Split( polyData, edgePointId1, edgePointId2, 1 );
            break;
        }
    }

    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;
}