Find 3D Model’s Boundary In VTK

We define all edges in the boundary have no neighbor cells. It indicates that there is no boundary on the closed surface.

In the following example, I read an STL file and get an arch.
Step 1, clip the original arch with a plane that has normal (0, 1, 0) to remove the upper part.
Step 2, find all lines which are used by only one cell.
Step 3, connect these lines to form long lists.
Step 4, find the longest list and show it.
Here are the code snippets.

main.cpp

#include <iostream>
#include <vtkSmartPointer.h>
#include <vtkSphereSource.h>
#include <vtkActor.h>
#include <vtkConeSource.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkPolyDataMapper.h>
#include <vtkSphereSource.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkProperty.h>
#include <vtkLineSource.h>
#include <vtkSTLReader.h>
#include <vtkPlaneSource.h>
#include <vtkCubeSource.h>
#include <vtkIdList.h>
#include <vtkCharArray.h>
#include <vtkCamera.h>
#include <vtkSTLWriter.h>
#include <vtkClipPolyData.h>
#include <vtkColorTransferFunction.h>
#include <vtkSTLReader.h>
#include <vtkPointData.h>
#include <vtkCutter.h>
#include <set>
#include <vtkPlane.h>

#include "ConnectedEdgeFilter.hpp"
#include "../share/tool.h"

using namespace std;

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

int main()
{
    setbuf( stdout, nullptr );
    vtkSPtrNew( reader, vtkSTLReader );
    reader->SetFileName( "/Users/weiyang/Desktop/u.stl" );
    reader->Update();

    vtkPolyData *polyData = reader->GetOutput();

    // Clip the top to make arch unclosed.
    vtkSPtrNew( plane, vtkPlane );
    polyData->ComputeBounds();
    double bounds[6];
    polyData->GetBounds( bounds );
    PointStruct origin( polyData->GetCenter() );
    origin[1] += (bounds[3] - bounds[2] ) / 2 - 2;
    plane->SetOrigin( origin.point );
    plane->SetNormal( 0, 1, 0 );
    vtkSPtrNew( clipper, vtkClipPolyData );
    clipper->SetInputData( polyData );
    clipper->SetClipFunction( plane );
    clipper->InsideOutOn();
    clipper->Update();

    polyData = clipper->GetOutput();

    vtkSPtrNew( mapper, vtkPolyDataMapper );
    mapper->SetInputData( polyData );

    vtkSPtrNew( actor, vtkActor );
    actor->SetMapper( mapper );

    // find all lines which has not neighbors
    vtkSPtrNew( boundaryArray, vtkCellArray );
    FindIndependentEdges( boundaryArray, polyData );

    // Define viewport ranges
    // (xmin, ymin, xmax, ymax) left-top and right-bottom point
    double leftViewport[4] = {0.0, 0.0, 0.5, 1.0};
    double rightViewport[4] = {0.5, 0.0, 1.0, 1.0};

    // Setup renderers
    vtkSmartPointer<vtkRenderer> leftRenderer =
    vtkSmartPointer<vtkRenderer>::New();
    leftRenderer->SetViewport( leftViewport );
    leftRenderer->SetBackground( .6, .5, .4 );

    vtkSmartPointer<vtkRenderer> rightRenderer =
    vtkSmartPointer<vtkRenderer>::New();
    rightRenderer->SetViewport( rightViewport );
    rightRenderer->SetBackground( .4, .5, .6);

    leftRenderer->AddActor( actor );
    leftRenderer->SetBackground( 0, 0, 0 );

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

    vtkSPtr<vtkIdList> longestList = connectFilter->GetLongestList();
    std::cout << "longestList->GetNumberOfIds: " << longestList->GetNumberOfIds() << std::endl;

    vtkSPtrNew( edge, vtkPolyData );
    vtkSPtrNew( points, vtkPoints );
    vtkSPtrNew( lines, vtkCellArray );
    const int ptsCount = longestList->GetNumberOfIds();
    vtkIdType pts1[ptsCount];
    bool valid = true;
    for( int i = 0; i < longestList->GetNumberOfIds(); ++i )
    {
        vtkIdType id = longestList->GetId( i );
        if( id > polyData->GetNumberOfPoints() )
        {
            std::cout << "i: "<< i << ", bad id: " << id << std::endl;
            valid = false;
        }
        PointStruct pt( polyData->GetPoint( id ) );
        points->InsertNextPoint( pt.point );
        pts1[i] = i;
    }
    if( valid )
    {
        edge->SetPoints( points );
        lines->InsertNextCell( ptsCount, pts1 );
        edge->SetPolys( lines );
        edge->BuildCells();
        edge->BuildLinks();

        vtkSPtrNew( edgeMapper, vtkPolyDataMapper );
        edgeMapper->SetInputData( edge );

        vtkSPtrNew( edgeActor, vtkActor );
        edgeActor->SetMapper( edgeMapper );
        edgeActor->GetProperty()->SetColor( 1, 0, 0 );
        edgeActor->GetProperty()->SetLineWidth( 10 );

        rightRenderer->AddActor( edgeActor );
    }

    std::cout << polyData->GetNumberOfPoints() << " : " << longestList->GetNumberOfIds() << std::endl;

    vtkSPtrNew( renderWindow, vtkRenderWindow );
    renderWindow->AddRenderer( leftRenderer );
    renderWindow->AddRenderer( rightRenderer );

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

    vtkSPtrNew( camera, vtkCamera );
    leftRenderer->SetActiveCamera( camera );
    rightRenderer->SetActiveCamera( camera );
    leftRenderer->ResetCamera();
    rightRenderer->ResetCamera();

    renderWindow->SetSize(600, 300);
    renderWindow->Render();
    renderWindowInteractor->Start();

    delete connectFilter;
    return 0;
}

ConnectedEdgeFilter.hpp

#pragma once

#include <vtkCellArray.h>
#include <vtkIdList.h>

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

class ConnectedEdgeFilter
{
public:
    ConnectedEdgeFilter(){}
    ~ConnectedEdgeFilter()
    {
        m_ConnectedLists.clear();
    }
    void Initialise( vtkSPtr<vtkCellArray> edgeSegments )
    {
        m_EdgeSegments = edgeSegments;
        m_EdgeSegments->InitTraversal();
        m_ConnectedLists.clear();
        vtkSPtrNew( cell, vtkIdList );
        while( m_EdgeSegments->GetNextCell( cell ) )
        {
            m_EdgeStartPointIds.push_back( cell->GetId( 0 ) );
            m_EdgeEndPointIds.push_back( cell->GetId( 1 ) );
            m_EgdeHandleFlag.push_back( false );
        }
    }
    void BackSearch(int startIndex, std::vector<int>& connectedList)
    {
        if ( nullptr == m_EdgeSegments || m_EdgeSegments->GetNumberOfCells() < 1 )
        {
            printf( "%s, %d: bad data", __FILE__, __LINE__ );
            return;
        }
        int segmentStartPointId = m_EdgeStartPointIds[startIndex];
        int segmentEndPointId = m_EdgeEndPointIds[startIndex];
        connectedList.push_back( segmentStartPointId );
        connectedList.push_back( segmentEndPointId );
        m_EgdeHandleFlag[ startIndex ] = true;

        int nextId = segmentEndPointId;
        bool foundNext = true;
        while (foundNext)
        {
            foundNext = false;
            for (int i = 0; i< m_EdgeSegments->GetNumberOfCells(); i++)
            {
                if (m_EgdeHandleFlag[i])
                {
                    continue;
                }
                if (nextId == m_EdgeEndPointIds[i])
                {
                    nextId = m_EdgeStartPointIds[i];
                    connectedList.push_back( nextId );
                    m_EgdeHandleFlag[i] = true;
                    foundNext = true;
                    break;

                }
                else if (nextId == m_EdgeStartPointIds[i])
                {
                    nextId = m_EdgeEndPointIds[i];
                    connectedList.push_back( nextId );
                    m_EgdeHandleFlag[i] = true;
                    foundNext = true;
                    break;
                }
            }
        }
    }
    void FrontSearch(int startIndex, std::vector<int>& connectedList)
    {
        if ( nullptr == m_EdgeSegments || m_EdgeSegments->GetNumberOfCells() < 1 )
        {
            printf( "%s, %d: bad data", __FILE__, __LINE__ );
            return;
        }

        bool findPre = true;
        int preId = m_EdgeStartPointIds[startIndex];
        while ( findPre )
        {
            findPre = false;
            for (int i = 0; i < m_EdgeSegments->GetNumberOfCells(); i++ )
            {
                if (m_EgdeHandleFlag[i])
                {
                    continue;
                }
                if (preId == m_EdgeStartPointIds[i])
                {
                    preId = m_EdgeEndPointIds[i];
                    connectedList.push_back(preId);
                    m_EgdeHandleFlag[i] = true;
                    findPre = true;
                    break;

                }
                else if (preId == m_EdgeEndPointIds[i])
                {
                    preId = m_EdgeStartPointIds[i];
                    connectedList.push_back(preId);
                    m_EgdeHandleFlag[i] = true;
                    findPre = true;
                    break;
                }
            }
        }
    }
    void HandleEdges()
    {
        if ( nullptr == m_EdgeSegments || m_EdgeSegments->GetNumberOfCells() < 1 )
        {
            printf( "%s, %d: bad data", __FILE__, __LINE__ );
            return;
        }

        int edgeIndex = 0;
        while ( edgeIndex !=-1 )
        {
            std::vector<int> connectedEdge0;
            BackSearch( edgeIndex, connectedEdge0 );

            std::vector<int> connectedEdge1;
            FrontSearch( edgeIndex, connectedEdge1 );

            // connect both list
            vtkSPtrNew( connectedList, vtkIdList );
            for( int i = connectedEdge0.size()-1; i >= 0; --i )
            {
                connectedList->InsertNextId( connectedEdge0[i] );
            }
            for( int i = 0; i < connectedEdge1.size(); ++i )
            {
                connectedList->InsertNextId( connectedEdge1[i] );
            }

            m_ConnectedLists.push_back( connectedList );
            edgeIndex = -1;
            int numberOfSegments = m_EgdeHandleFlag.size();
            for (int i = 0; i<numberOfSegments; i++)
            {
                if (m_EgdeHandleFlag[i] == false)
                {
                    edgeIndex = i;
                    break;
                }
            }
        }
    }
    vtkSPtr<vtkIdList> GetLongestList()
    {
        vtkSPtrNew( longestList, vtkIdList );
        int max = -1;
        for (int i = 0; i < m_ConnectedLists.size(); i++)
        {
            vtkIdList *list = m_ConnectedLists[i];
            vtkIdType listSize = list->GetNumberOfIds();
            if ( listSize > max )
            {
                max = listSize;
                longestList = list;
            }
        }
        return longestList;
    }
protected:
    vtkSPtr<vtkCellArray>               m_EdgeSegments;

    std::vector<bool>           m_EgdeHandleFlag;
    std::vector<int>            m_EdgeStartPointIds;
    std::vector<int>            m_EdgeEndPointIds;

    std::vector<vtkSPtr<vtkIdList>>     m_ConnectedLists;

};

We can compare the original arch and boundary in both renderers.

Move the cursor to the right renderer and tap the key w to switch it to modify the representation of all actors so that they are wireframe.

We can show all lists in a simple but rude way like the following code.

    //check all segments that contain longest list

    for( int i = 0; i < boundaryArray->GetNumberOfCells(); ++i )
    {
        vtkSPtrNew( pts, vtkIdList );
        boundaryArray->GetCell( i, pts );
        std::cout << "GetNumberOfIds: " << pts->GetNumberOfIds() << std::endl;
        vtkSPtrNew( edge, vtkPolyData );
        vtkSPtrNew( points, vtkPoints );
        vtkSPtrNew( lines, vtkCellArray );
        const int ptsCount = pts->GetNumberOfIds();
        vtkIdType pts1[ptsCount];
        bool valid = true;
        for( int j = 0; j < pts->GetNumberOfIds(); ++j )
        {
            vtkIdType id = pts->GetId( j );
            if( id > polyData->GetNumberOfPoints() )
            {
                std::cout << "i: "<< j << ", bad id: " << id << std::endl;
                valid = false;
                break;
            }
            points->InsertNextPoint( polyData->GetPoint( id ) );
            pts1[j] = j;
        }
        if( !valid )
        {
            continue;
        }
        edge->SetPoints( points );
        lines->InsertNextCell( ptsCount, pts1 );
        edge->SetPolys( lines );

        vtkSPtrNew( edgeMapper, vtkPolyDataMapper );
        edgeMapper->SetInputData( edge );

        vtkSPtrNew( edgeActor, vtkActor );
        edgeActor->SetMapper( edgeMapper );
        edgeActor->GetProperty()->SetColor( 1, 0, 0 );

        rightRenderer->AddActor( edgeActor );
    }

Output:

longestList->GetNumberOfIds: 545
GetNumberOfIds: 2
GetNumberOfIds: 9253
GetNumberOfIds: 9254
i: 3517, bad id: 4371982720
GetNumberOfIds: 2
GetNumberOfIds: 9254
i: 3515, bad id: 4371982720
GetNumberOfIds: 9255
i: 3514, bad id: 4371982720
GetNumberOfIds: 2
GetNumberOfIds: 9256
i: 3512, bad id: 4371982720
GetNumberOfIds: 9257
i: 3511, bad id: 4371982720
GetNumberOfIds: 2
GetNumberOfIds: 9257
i: 3509, bad id: 4371982720
GetNumberOfIds: 9258
i: 3508, bad id: 4371982720
...