This post introduces a way to find the biggest plane on the 3D model.
We expand the initial cell by comparing its normal with neighbor’s normal. It indicates they are on a plane if the two normals of cells are approximately equal.
I use vtkOBBTree to find the min axes of the arch to ignore other cells that normals are not equal to the min axes. You can read the article Explore The Orientation Of 3D Object about the min axes.
Finally I mark the cells on the plane red color by scalar to show the algorithm’s result.
#include <vtkVersion.h>
#include <vtkSmartPointer.h>
#include <vtkSurfaceReconstructionFilter.h>
#include <vtkProgrammableSource.h>
#include <vtkContourFilter.h>
#include <vtkReverseSense.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkPolyData.h>
#include <vtkCamera.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkSphereSource.h>
#include <vtkXMLPolyDataReader.h>
#include <vtkUnstructuredGrid.h>
#include <vtkDataSetMapper.h>
#include <vtkPolyDataNormals.h>
#include <vtkPointData.h>
#include <vtkCellData.h>
#include <vtkConeSource.h>
#include <vtkTransform.h>
#include <vtkTransformPolyDataFilter.h>
#include <vtkGlyph3D.h>
#include <vtkSTLReader.h>
#include <vtkLineSource.h>
#include <vtkOBBTree.h>
#include <QList>
#include <vtkTriangleFilter.h>
#include <vtkAppendPolyData.h>
#include <vtkDecimatePro.h>
#include "../share/tool.h"
#define vSP vtkSmartPointer
#define vSPNew(Var, Type) vSP<Type> Var = vSP<Type>::New();
struct CellInfo
{
PointStruct normal;
PointStruct pts[3];
CellInfo operator=(CellInfo anotherCell)
{
normal = anotherCell.normal;
pts[0] = anotherCell.pts[0];
pts[1] = anotherCell.pts[1];
pts[2] = anotherCell.pts[2];
return *this;
}
// heron math way
double GetCellArea()
{
double ds1 = (pts[0] - pts[1]).Length();
double ds2 = (pts[0] - pts[2]).Length();
double ds3 = (pts[1] - pts[2]).Length();
double p = (ds1 + ds2 + ds3) / 2;
double S = sqrt(p * (p - ds1) * (p - ds2) * (p - ds3));
return S;
}
};
void HandleNeighbourCells(vSP<vtkPolyData> mesh,
int seedCellId,
double angleDiffLimit,
std::vector<int>& traveledCellIds,
vSP<vtkDataArray> cellNormals,
bool* cellsVisited)
{
if( cellsVisited[seedCellId] )
{
return ;
}
const vtkIdType numberOfCells = mesh->GetNumberOfCells();
int topStack = 0;
int* TravelStack = new int[numberOfCells];
traveledCellIds.clear();
TravelStack[topStack++] = seedCellId;
vSPNew(cellPointIds, vtkIdList);
vSPNew(neighborCellIds, vtkIdList);
double neighborNormals[2][3];
cellNormals->GetTuple(seedCellId, neighborNormals[0]);
while (topStack > 0)
{
topStack--;
int curSeedId = TravelStack[topStack];
traveledCellIds.push_back(curSeedId);
cellsVisited[curSeedId] = true;
mesh->GetCellPoints(curSeedId, cellPointIds);
int ptsCount = cellPointIds->GetNumberOfIds();
if( 3 != ptsCount )
{
cout << "ptsCount: " << ptsCount << endl;
}
for (int i = 0; i < ptsCount; i++)
{
mesh->GetCellEdgeNeighbors(curSeedId, cellPointIds->GetId(i), cellPointIds->GetId((i + 1) % ptsCount), neighborCellIds);
if (neighborCellIds->GetNumberOfIds() > 0)
{
int neighborCellId = neighborCellIds->GetId(0);
if (cellsVisited[neighborCellId] == false)
{
//cellsVisited[neighborCellId] = true;
cellNormals->GetTuple(neighborCellId, neighborNormals[1]);
double angOfTwoCell = vtkMath::DegreesFromRadians(vtkMath::AngleBetweenVectors(neighborNormals[0], neighborNormals[1]));
if (angOfTwoCell < angleDiffLimit)
{
TravelStack[topStack++] = neighborCellId;
}
}
}
}
}
delete []TravelStack;
}
int main(int argc, char *argv[])
{
vSPNew( STLReader, vtkSTLReader );
STLReader->SetFileName( "~/arch_u.stl" );
STLReader->Update();
vtkPolyData *polyData = STLReader->GetOutput();
vSPNew( decimate, vtkDecimatePro );
decimate->SetTargetReduction( 0.3 );
decimate->SetInputData( polyData );
decimate->Update();
vSPNew( polydata, vtkPolyData );
polydata->DeepCopy( decimate->GetOutput() );
vtkSmartPointer<vtkTriangleFilter> triangle = vtkSmartPointer<vtkTriangleFilter>::New();
triangle->SetInputData( polyData );
triangle->PassLinesOff();
triangle->PassVertsOff();
triangle->Update();
polyData = triangle->GetOutput();
polyData->BuildCells();
polyData->BuildLinks();
PointStruct center( polyData->GetCenter() );
vtkSmartPointer<vtkPolyDataMapper> mapper =
vtkSmartPointer<vtkPolyDataMapper>::New();
mapper->SetInputData( polyData );
mapper->Update();
vtkSmartPointer<vtkActor> surfaceActor =
vtkSmartPointer<vtkActor>::New();
surfaceActor->SetMapper( mapper );
vtkSmartPointer<vtkPolyDataNormals> pdNormals =
vtkSmartPointer<vtkPolyDataNormals>::New();
pdNormals->SetInputData( polyData );
pdNormals->ComputeCellNormalsOn();
pdNormals->ComputePointNormalsOff();
pdNormals->Update();
vtkDataArray *cellNormals = pdNormals->GetOutput()->GetCellData()->GetNormals();
cout << "cellNormals GetNumberOfTuples: " << cellNormals->GetNumberOfTuples() << endl;
cout << polyData->GetNumberOfPoints() << " : " << polyData->GetNumberOfCells() << endl;
const int numberOfCells = cellNormals->GetNumberOfTuples();
PointStruct corner, max, mid, min, size;
PointStruct realSize;
vSPNew( obbTree, vtkOBBTree );
obbTree->ComputeOBB( polyData, corner.point, max.point, mid.point, min.point, size.point );
realSize[0] = max.Length();
realSize[1] = mid.Length();
realSize[2] = min.Length();
max.Unit();
mid.Unit();
min.Unit();
// =========== colorful =============
vtkSmartPointer<vtkUnsignedCharArray> scalar = vtkSmartPointer<vtkUnsignedCharArray>::New();
scalar = vSP<vtkUnsignedCharArray>::New();
scalar->SetNumberOfComponents(3);
scalar->SetName("Colors");
unsigned char white[3] = { 255, 255, 255 };
for (int i = 0; i < numberOfCells; i++)
{
scalar->InsertNextTypedTuple(white);
}
surfaceActor->GetMapper()->GetInput()->GetCellData()->SetScalars(scalar);
// =========== fnished: colorful =============
// =========== start to work ===========
bool* cellsVisited = new bool[numberOfCells];
for (int i = 0; i < numberOfCells; i++)
{
cellsVisited[i] = false;
}
//int maxCellCount = 0;
double maxArea = 0;
std::vector<int> traveledCellIdsForBigPlane;
for( int i = 0; i < numberOfCells; ++i )
{
vtkCell *cell = polyData->GetCell( i );
PointStruct pt( cell->GetPoints()->GetPoint( 0 ) );
PointStruct vec = pt - center;
double dis = fabs( vec.Dot( min ) );
if( dis / realSize[2] < 0.3 )
{
continue;
}
std::vector<int> traveledCellIds;
HandleNeighbourCells( polyData, i, 1, traveledCellIds, cellNormals, cellsVisited );
double area = 0;
for( int j = 0; j < traveledCellIds.size(); ++j )
{
int cellId = traveledCellIds[j];
vtkCell *cell = polyData->GetCell( cellId );
vtkPoints *points = cell->GetPoints();
CellInfo cellInfo;
cellInfo.normal = PointStruct(cellNormals->GetTuple(cellId));
for( int k = 0; k < 3; ++k )
{
cellInfo.pts[k] = PointStruct(points->GetPoint(k));
}
area += cellInfo.GetCellArea();
}
if( maxArea < area )
{
maxArea = area;
traveledCellIdsForBigPlane = traveledCellIds;
}
// if( traveledCellIds.size() > maxCellCount )
// {
// maxCellCount = traveledCellIds.size();
// traveledCellIdsForBigPlane = traveledCellIds;
// }
}
delete []cellsVisited;
// =========== finished: find biggest area ===========
for( int i = 0; i < traveledCellIdsForBigPlane.size(); ++i )
{
scalar->SetTuple3( traveledCellIdsForBigPlane[i], 255, 0, 0 );
}
cout << "maxArea: " << maxArea << endl;
cout << "traveledCellIdsForBigPlane cells count: " << traveledCellIdsForBigPlane.size() << endl;
// Create the RenderWindow, Renderer and both Actors
vtkSmartPointer<vtkRenderer> ren =
vtkSmartPointer<vtkRenderer>::New();
vtkSmartPointer<vtkRenderWindow> renWin =
vtkSmartPointer<vtkRenderWindow>::New();
renWin->AddRenderer(ren);
vtkSmartPointer<vtkRenderWindowInteractor> iren =
vtkSmartPointer<vtkRenderWindowInteractor>::New();
iren->SetRenderWindow(renWin);
// Add the actors to the renderer, set the background and size
ren->AddActor( surfaceActor );
ren->SetBackground(.2, .3, .4);
renWin->Render();
iren->Start();
return EXIT_SUCCESS;
}
