Create a sphere and clip it as the following image.
We will pull down the boundary and close the opening to get a new model.
How to pull down the boundary of the hemispheres?
Step 1. Find all independent edges which are used by only one cell.
Step 2. Connect all independent edges and forms a list.
Step 3. Project all points on the list to the bottom plane.
Step 4. Link the points on the list and projected points to generate the surface.
void BuildSurface
(
vtkSPtr<vtkPolyData> polyData,
vtkSPtr<vtkIdList> holeList,
vtkSPtr<vtkPoints> projectedPts
)
{
if( polyData->NeedToBuildCells() )
{
polyData->BuildCells();
polyData->BuildLinks();
}
vtkIdType proPtsCount = projectedPts->GetNumberOfPoints();
vtkIdType dataPtsCount = polyData->GetNumberOfPoints();
bool isClosed = false;
if( holeList->GetId(0) == holeList->GetId( holeList->GetNumberOfIds() - 1 ) )
{
isClosed = true;
proPtsCount = proPtsCount - 1;
}
// add points
for( int i = 0; i < proPtsCount; ++i )
{
double ptPos[3];
projectedPts->GetPoint( i, ptPos );
polyData->InsertNextLinkedPoint( ptPos, 1 );
}
// add triangle as cell
for( int i = 0; i < proPtsCount - 1; ++i )
{
vtkIdType preId1 = holeList->GetId( i );
vtkIdType preId2 = holeList->GetId( i+1 );
vtkIdType nextId1 = dataPtsCount + i;
vtkIdType nextId2 = nextId1 + 1;
vtkIdType tri1[3] = { preId1, preId2, nextId1 };
vtkIdType tri2[3] = { nextId2, nextId1, preId2 };
polyData->InsertNextLinkedCell( VTK_TRIANGLE, 3, tri1 );
polyData->InsertNextLinkedCell( VTK_TRIANGLE, 3, tri2 );
}
if( isClosed )
{
vtkIdType preId1 = holeList->GetId( proPtsCount - 1 );
vtkIdType preId2 = holeList->GetId( proPtsCount );
vtkIdType nextId1 = dataPtsCount + proPtsCount - 1;
vtkIdType nextId2 = dataPtsCount;
vtkIdType tri1[3] = { preId1, preId2, nextId1 };
vtkIdType tri2[3] = { nextId2, nextId1, preId2 };
polyData->InsertNextLinkedCell( VTK_TRIANGLE, 3, tri1 );
polyData->InsertNextLinkedCell( VTK_TRIANGLE, 3, tri2 );
}
polyData->Modified();
}
int main()
{
// ...
vtkSPtrNew( bottomPlane, vtkPlane );
bottomPlane->SetOrigin( 0, -2, 0 );
bottomPlane->SetNormal( 0, 1, 0 );
vtkSPtrNew( projectPoints, vtkPoints );
for( int i = 0; i < longestList->GetNumberOfIds(); ++i )
{
vtkIdType id = longestList->GetId( i );
double *pt = originalPd->GetPoint( id );
double projectedPt[3];
bottomPlane->ProjectPoint( pt, projectedPt );
projectPoints->InsertNextPoint( projectedPt );
}
//...
Finally, we have to fill the hole. VTK provides the algorithm vtkFillHolesFilter to help us to fulfill the task.
All implementation details can be found in the following code.
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 <vtkRenderWindowInteractor.h>
#include <vtkClipPolyData.h>
#include <vtkPlane.h>
#include <vtkPoints.h>
#include <vtkSTLWriter.h>
#include <vtkFillHolesFilter.h>
#include <vtkProperty.h>
#include "ConnectedEdgeFilter.hpp"
#define vtkSPtr vtkSmartPointer
#define vtkSPtrNew(Var, Type) vtkSPtr<Type> Var = vtkSPtr<Type>::New();
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;
}
}
}
}
void BuildSurface
(
vtkSPtr<vtkPolyData> polyData,
vtkSPtr<vtkIdList> holeList,
vtkSPtr<vtkPoints> projectedPts
)
{
if( polyData->NeedToBuildCells() )
{
polyData->BuildCells();
polyData->BuildLinks();
}
vtkIdType proPtsCount = projectedPts->GetNumberOfPoints();
vtkIdType dataPtsCount = polyData->GetNumberOfPoints();
bool isClosed = false;
if( holeList->GetId(0) == holeList->GetId( holeList->GetNumberOfIds() - 1 ) )
{
isClosed = true;
proPtsCount = proPtsCount - 1;
}
// add points
for( int i = 0; i < proPtsCount; ++i )
{
double ptPos[3];
projectedPts->GetPoint( i, ptPos );
polyData->InsertNextLinkedPoint( ptPos, 1 );
}
// add triangle as cell
for( int i = 0; i < proPtsCount - 1; ++i )
{
vtkIdType preId1 = holeList->GetId( i );
vtkIdType preId2 = holeList->GetId( i+1 );
vtkIdType nextId1 = dataPtsCount + i;
vtkIdType nextId2 = nextId1 + 1;
vtkIdType tri1[3] = { preId1, preId2, nextId1 };
vtkIdType tri2[3] = { nextId2, nextId1, preId2 };
polyData->InsertNextLinkedCell( VTK_TRIANGLE, 3, tri1 );
polyData->InsertNextLinkedCell( VTK_TRIANGLE, 3, tri2 );
}
if( isClosed )
{
vtkIdType preId1 = holeList->GetId( proPtsCount - 1 );
vtkIdType preId2 = holeList->GetId( proPtsCount );
vtkIdType nextId1 = dataPtsCount + proPtsCount - 1;
vtkIdType nextId2 = dataPtsCount;
vtkIdType tri1[3] = { preId1, preId2, nextId1 };
vtkIdType tri2[3] = { nextId2, nextId1, preId2 };
polyData->InsertNextLinkedCell( VTK_TRIANGLE, 3, tri1 );
polyData->InsertNextLinkedCell( VTK_TRIANGLE, 3, tri2 );
}
polyData->Modified();
}
int main()
{
vtkSPtrNew( sphere, vtkSphereSource );
sphere->SetRadius( 1 );
sphere->Update();
vtkSPtrNew( cutPlane, vtkPlane );
cutPlane->SetOrigin( 0, 0, 0 );
cutPlane->SetNormal( 0, 1, 0 );
vtkSPtrNew( clipFilter, vtkClipPolyData );
clipFilter->SetInputData( sphere->GetOutput() );
clipFilter->SetClipFunction( cutPlane );
clipFilter->Update();
vtkPolyData *originalPd = clipFilter->GetOutput();
// ============ find all lines which has not neighbors ==============
vtkSPtrNew( boundaryArray, vtkCellArray );
FindIndependentEdges( boundaryArray, originalPd );
ConnectedEdgeFilter *connectFilter = new ConnectedEdgeFilter();
connectFilter->Initialise( boundaryArray );
connectFilter->HandleEdges();
vtkSPtr<vtkIdList> longestList = connectFilter->GetLongestList();
// =============== pull down the boundary ===================
vtkSPtrNew( bottomPlane, vtkPlane );
bottomPlane->SetOrigin( 0, -2, 0 );
bottomPlane->SetNormal( 0, 1, 0 );
vtkSPtrNew( projectPoints, vtkPoints );
for( int i = 0; i < longestList->GetNumberOfIds(); ++i )
{
vtkIdType id = longestList->GetId( i );
double *pt = originalPd->GetPoint( id );
double projectedPt[3];
bottomPlane->ProjectPoint( pt, projectedPt );
projectPoints->InsertNextPoint( projectedPt );
}
BuildSurface( originalPd, longestList, projectPoints );
// =============== fill the hole ===================
vtkSPtrNew( fillHoleFilter, vtkFillHolesFilter );
fillHoleFilter->SetInputData( originalPd );
fillHoleFilter->SetHoleSize( 10000 );
fillHoleFilter->Update();
vtkSPtrNew( mapper, vtkPolyDataMapper );
mapper->SetInputData( fillHoleFilter->GetOutput() );
vtkSPtrNew( actor, vtkActor );
actor->SetMapper( mapper );
vtkSPtrNew( renderer, vtkRenderer );
renderer->AddActor(actor);
renderer->SetBackground( 0, 0, 0 );
vtkSPtrNew( renderWindow, vtkRenderWindow );
renderWindow->AddRenderer( renderer );
vtkSPtrNew( renderWindowInteractor, vtkRenderWindowInteractor );
renderWindowInteractor->SetRenderWindow( renderWindow );
renderer->ResetCamera();
renderWindow->Render();
renderWindowInteractor->Start();
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;
};
