VTK – Fill Hole By Center Point Suture Method

Create a sphere and use a plane to cut it to generate a big hole. This is our test data.

Then we will fill the hole by inserting a center point of the boundary and creating new triangle cells to form a new surface.

We will use ConnectedEdgeFilter.hpp to handles edges on the boundary to get a connected list. You can find the file at: https://www.weiy.city/2020/03/find-3d-models-boundary-in-vtk/.

Finally, fill all holes with a circumference greater than 1.5.

Here are all code:

#include <iostream>
#include <vtkSmartPointer.h>
#include <vtkTransform.h>
#include <vtkActor.h>
#include <vtkConeSource.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkPolyDataMapper.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkAxesActor.h>
#include <vtkLine.h>
#include <vtkPlane.h>
#include <vtkSTLWriter.h>
#include <vtkPointLocator.h>
#include <vtkProperty.h>
#include <vtkSphereSource.h>
#include <vtkCubeSource.h>
#include <vtkAxesActor.h>
#include <vtkCallbackCommand.h>
#include <vtkPlaneSource.h>
#include <vtkClipPolyData.h>

#include "point.hpp"
#include "ConnectedEdgeFilter.hpp"

using namespace std;

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

void FindIndependentEdges( vtkSPtr<vtkCellArray> edges, vtkSPtr<vtkPolyData> originPolyData )
{
    vtkIdType cellId;
    vtkIdType npts;
    vtkIdType p1, p2;
    vtkIdType *pts = nullptr;
    edges->Initialize();
    originPolyData->BuildLinks();
    vtkSPtrNew( neighbors, vtkIdList );
    neighbors->Allocate(VTK_CELL_SIZE);
    vtkCellArray *polys = originPolyData->GetPolys();
    for (cellId=0, polys->InitTraversal();
         polys->GetNextCell(npts,pts);
         cellId++)
    {
        for (int i = 0; i < npts; i++ )
        {
            p1 = pts[ i ];
            p2 = pts[ (i+1) %npts ];

            originPolyData->GetCellEdgeNeighbors( cellId, p1, p2, neighbors );
            vtkIdType numNei = neighbors->GetNumberOfIds();
            if (numNei < 1)
            {
                vtkIdType edgeIds[2] = { p1, p2 };
                edges->InsertNextCell( 2, edgeIds );
            }
            else
            {
                continue;
            }
        }
    }
}

void CenterPointFillHole(vtkSmartPointer<vtkPolyData> polyData, vtkSmartPointer<vtkIdList> list)
{
    Point center( 0, 0, 0 );
    for( int i = 0; i < list->GetNumberOfIds() - 1; ++i )
    {
        int id = list->GetId( i );
        Point pt( polyData->GetPoint( id ) );
        center = center + pt;
    }
    center /= (list->GetNumberOfIds() - 1);
    vtkIdType cellCenterPointId = polyData->InsertNextLinkedPoint(center.point, list->GetNumberOfIds() - 1);
    for( int i = 0; i < list->GetNumberOfIds() - 1; ++i )
    {
        int ptId0 = list->GetId( i );
        int ptId1 = list->GetId( i+1 );
        vtkIdType tri1[3] = {  ptId0, ptId1, cellCenterPointId };
        int cellId1 = polyData->InsertNextLinkedCell(VTK_TRIANGLE, 3, tri1);
    }
    polyData->Modified();
}

void CenterPointFillAllBigHoles(vtkSmartPointer<vtkPolyData> polyData, double holeLengthLimit)
{
    // find all lines which has not neighbors
    vtkSPtrNew( boundaryArray, vtkCellArray );
    FindIndependentEdges( boundaryArray, polyData );

    ConnectedEdgeFilter *connectFilter = new ConnectedEdgeFilter();
    connectFilter->Initialise( boundaryArray );
    connectFilter->HandleEdges();

    for( int i = 0; i < connectFilter->GetListsCount(); ++i )
    {
        vtkSmartPointer<vtkIdList> list = connectFilter->GetList(i);
        double edgeLength = 0.0;
        if (list->GetNumberOfIds() < 3)
            continue;

        // it's a closed list
        std::cout << list->GetId( 0 ) << ", " << list->GetId( list->GetNumberOfIds() - 1 ) << std::endl;

        for (vtkIdType j = 1; j < list->GetNumberOfIds(); j++)
        {
            Point prePoint, point;
            polyData->GetPoint(list->GetId(j), point.point);
            polyData->GetPoint(list->GetId(j - 1), prePoint.point);

            double length = (point - prePoint).Length();
            edgeLength += length;
        }
        if( edgeLength > holeLengthLimit )
        {
            CenterPointFillHole( polyData, list );
        }
    }
}


int main()
{
    setbuf( stdout, nullptr );

    vtkSmartPointer<vtkSphereSource> source =
            vtkSmartPointer<vtkSphereSource>::New();
    source->Update();

    vtkSmartPointer<vtkPlane> cutPlane =
            vtkSmartPointer<vtkPlane>::New();
    cutPlane->SetOrigin( 0, 0, 0 );
    cutPlane->SetNormal( 1, 0, 0 );

    vtkSmartPointer<vtkClipPolyData> clip =
            vtkSmartPointer<vtkClipPolyData>::New();
    clip->SetClipFunction( cutPlane );
    clip->SetInputData( source->GetOutput() );
    clip->Update();

    vtkSmartPointer<vtkPolyData> clippedData = clip->GetOutput();
    clippedData->BuildCells();
    clippedData->BuildLinks();

    CenterPointFillAllBigHoles( clippedData, 1.5 );

    vtkSmartPointer<vtkSTLWriter> writer = vtkSmartPointer<vtkSTLWriter>::New();
    writer->SetInputData( clippedData );
    writer->SetFileName( "out.stl" );
    writer->Write();

    vtkSmartPointer<vtkPolyDataMapper> mapper =
            vtkSmartPointer<vtkPolyDataMapper>::New();
    mapper->SetInputData( clippedData );

    vtkSmartPointer<vtkActor> actor =
            vtkSmartPointer<vtkActor>::New();
    actor->SetMapper( mapper );

    vtkSmartPointer<vtkAxesActor> axesActor =
            vtkSmartPointer<vtkAxesActor>::New();

    vtkSmartPointer<vtkRenderer> renderer =
            vtkSmartPointer<vtkRenderer>::New();
    renderer->AddActor( actor );
    //renderer->AddActor( axesActor );

    renderer->SetBackground( 0, 0, 0 );

    vtkSmartPointer<vtkRenderWindow> renderWindow =
            vtkSmartPointer<vtkRenderWindow>::New();
    renderWindow->AddRenderer( renderer );

    vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
            vtkSmartPointer<vtkRenderWindowInteractor>::New();
    renderWindowInteractor->SetRenderWindow( renderWindow );

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

    renderWindowInteractor->Start();
    return 0;
}

Result:

Shows all normals on Paraview by NormalGlyphs filter: