求两个多边形数据 vtkPolyData 的相交线

vtkIntersectionPolyDataFilter

【简介】

该滤波器计算两个多边形数据 vtkPolyData 的交集。

其第一个输出是 交集线集( a set of lines )；

第二个和第三个输出分别为第一和第二个输入 vtkPolyData.

如果需要，后两个输出可以被相交线分割。

【示例】

#include <vtkActor.h>#include <vtkIntersectionPolyDataFilter.h>#include <vtkPolyDataMapper.h>#include <vtkProperty.h>#include <vtkRenderer.h>#include <vtkRenderWindow.h>#include <vtkRenderWindowInteractor.h>#include <vtkSmartPointer.h>#include <vtkSphereSource.h>int main(int, char *[]){vtkSmartPointer<vtkSphereSource> sphereSource1 = vtkSmartPointer<vtkSphereSource>::New();sphereSource1->SetCenter(0.0, 0.0, 0.0);sphereSource1->SetRadius(2.0f);sphereSource1->Update();vtkSmartPointer<vtkPolyDataMapper> sphere1Mapper = vtkSmartPointer<vtkPolyDataMapper>::New();sphere1Mapper->SetInputConnection( sphereSource1->GetOutputPort() );sphere1Mapper->ScalarVisibilityOff();vtkSmartPointer<vtkActor> sphere1Actor = vtkSmartPointer<vtkActor>::New();sphere1Actor->SetMapper( sphere1Mapper );//sphere1Actor->GetProperty()->SetOpacity(.7);sphere1Actor->GetProperty()->SetColor(1,0,0);sphere1Actor->GetProperty()->SetRepresentationToWireframe();vtkSmartPointer<vtkSphereSource> sphereSource2 = vtkSmartPointer<vtkSphereSource>::New();sphereSource2->SetCenter(1.0, 0.0, 0.0);sphereSource2->SetRadius(2.0f);vtkSmartPointer<vtkPolyDataMapper> sphere2Mapper = vtkSmartPointer<vtkPolyDataMapper>::New();sphere2Mapper->SetInputConnection( sphereSource2->GetOutputPort() );sphere2Mapper->ScalarVisibilityOff();vtkSmartPointer<vtkActor> sphere2Actor = vtkSmartPointer<vtkActor>::New();sphere2Actor->SetMapper( sphere2Mapper );//sphere2Actor->GetProperty()->SetOpacity(.7);sphere2Actor->GetProperty()->SetColor(0,1,0);sphere2Actor->GetProperty()->SetRepresentationToWireframe();vtkSmartPointer<vtkIntersectionPolyDataFilter> intersectionPolyDataFilter =vtkSmartPointer<vtkIntersectionPolyDataFilter>::New();intersectionPolyDataFilter->SetInputConnection( 0, sphereSource1->GetOutputPort() );intersectionPolyDataFilter->SetInputConnection( 1, sphereSource2->GetOutputPort() );intersectionPolyDataFilter->Update();vtkSmartPointer<vtkPolyDataMapper> intersectionMapper = vtkSmartPointer<vtkPolyDataMapper>::New();intersectionMapper->SetInputConnection( intersectionPolyDataFilter->GetOutputPort() );intersectionMapper->ScalarVisibilityOff();vtkSmartPointer<vtkActor> intersectionActor = vtkSmartPointer<vtkActor>::New();intersectionActor->SetMapper( intersectionMapper );vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();renderer->AddViewProp(sphere1Actor);renderer->AddViewProp(sphere2Actor);renderer->AddViewProp(intersectionActor);intersectionActor->GetProperty()->SetColor(0,0,1);vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();renderWindow->AddRenderer( renderer );vtkSmartPointer<vtkRenderWindowInteractor> renWinInteractor =vtkSmartPointer<vtkRenderWindowInteractor>::New();renWinInteractor->SetRenderWindow( renderWindow );renderWindow->Render();renWinInteractor->Start();return EXIT_SUCCESS;}

可以看出使用是非常简单的，

【类实现部分】

vtkSmartPointer<vtkIntersectionPolyDataFilter> intersectionPolyDataFilter =vtkSmartPointer<vtkIntersectionPolyDataFilter>::New();intersectionPolyDataFilter->SetInputConnection( 0, sphereSource1->GetOutputPort() );intersectionPolyDataFilter->SetInputConnection( 1, sphereSource2->GetOutputPort() );intersectionPolyDataFilter->Update();

输入两个多边形数据信息，调用Update() 函数就可以得到像个多边形数据相交的数据信息：

intersectionMapper->SetInputConnection( <span style="background-color: rgb(255, 153, 255);">intersectionPolyDataFilter->GetOutputPort()</span> );

【类实现过程的分析】

其实内部还有很多东西是可以学习的，接下来就分析一下 这个类的源代码：

可以看出该类的基类是 vtkPolyDataAlgorithm ，所以这个类中最为重要的函数就是 

intRequestData (vtkInformation *,vtkInformationVector **, vtkInformationVector *)通过前面文章的分析可以知道，该类调用 Update()  函数的过程中系统会自动调用这个 RequestData() 这个函数。

该类的实现部分：

int vtkIntersectionPolyDataFilter::RequestData(vtkInformation* vtkNotUsed(request),vtkInformationVector** inputVector,vtkInformationVector*  outputVector){vtkInformation* inInfo0 = inputVector[0]->GetInformationObject(0);vtkInformation* inInfo1 = inputVector[1]->GetInformationObject(0);vtkInformation* outIntersectionInfo = outputVector->GetInformationObject(0);vtkInformation* outPolyDataInfo0 = outputVector->GetInformationObject(1);vtkInformation* outPolyDataInfo1 = outputVector->GetInformationObject(2);vtkPolyData *input0 = vtkPolyData::SafeDownCast(inInfo0->Get(vtkDataObject::DATA_OBJECT()));vtkPolyData *input1 = vtkPolyData::SafeDownCast(inInfo1->Get(vtkDataObject::DATA_OBJECT()));vtkPolyData *outputIntersection = vtkPolyData::SafeDownCast(outIntersectionInfo->Get(vtkDataObject::DATA_OBJECT()));vtkSmartPointer< vtkPoints > outputIntersectionPoints =vtkSmartPointer< vtkPoints >::New();outputIntersection->SetPoints(outputIntersectionPoints);vtkPolyData *outputPolyData0 = vtkPolyData::SafeDownCast(outPolyDataInfo0->Get(vtkDataObject::DATA_OBJECT()));vtkPolyData *outputPolyData1 = vtkPolyData::SafeDownCast(outPolyDataInfo1->Get(vtkDataObject::DATA_OBJECT()));// Set up new poly data for the inputs to build cells and links.vtkSmartPointer< vtkPolyData > mesh0 = vtkSmartPointer< vtkPolyData >::New();mesh0->DeepCopy(input0);vtkSmartPointer< vtkPolyData > mesh1 = vtkSmartPointer< vtkPolyData >::New();mesh1->DeepCopy(input1);// Find the triangle-triangle intersections between mesh0 and mesh1vtkSmartPointer< vtkOBBTree > obbTree0 = vtkSmartPointer< vtkOBBTree >::New();obbTree0->SetDataSet(mesh0);obbTree0->SetNumberOfCellsPerNode(10);obbTree0->SetMaxLevel(1000000);obbTree0->SetTolerance(1e-6);obbTree0->AutomaticOn();obbTree0->BuildLocator();vtkSmartPointer< vtkOBBTree > obbTree1 = vtkSmartPointer< vtkOBBTree >::New();obbTree1->SetDataSet(mesh1);obbTree1->SetNumberOfCellsPerNode(10);obbTree1->SetMaxLevel(1000000);obbTree1->SetTolerance(1e-6);obbTree1->AutomaticOn();obbTree1->BuildLocator();// Set up the structure for determining exact triangle-triangle intersections.vtkIntersectionPolyDataFilter::Impl *impl = new vtkIntersectionPolyDataFilter::Impl();impl->Mesh[0]  = mesh0;impl->Mesh[1]  = mesh1;impl->OBBTree1 = obbTree1;vtkSmartPointer< vtkCellArray > lines = vtkSmartPointer< vtkCellArray >::New();outputIntersection->SetLines(lines);impl->IntersectionLines = lines;// Add cell data arrays that map the intersection line to the cells// it splits.impl->CellIds[0] = vtkIdTypeArray::New();impl->CellIds[0]->SetName("Input0CellID");outputIntersection->GetCellData()->AddArray(impl->CellIds[0]);impl->CellIds[0]->Delete();impl->CellIds[1] = vtkIdTypeArray::New();impl->CellIds[1]->SetName("Input1CellID");outputIntersection->GetCellData()->AddArray(impl->CellIds[1]);impl->CellIds[1]->Delete();impl->PointCellIds[0] = vtkIdTypeArray::New();impl->PointCellIds[0]->SetName("PointCellsIDs");impl->PointCellIds[1] = vtkIdTypeArray::New();impl->PointCellIds[1]->SetName("PointCellsIDs");double bounds0[6], bounds1[6];mesh0->GetBounds(bounds0);mesh1->GetBounds(bounds1);for (int i = 0; i < 3; i++){int minIdx = 2*i;int maxIdx = 2*i+1;if (bounds1[minIdx] < bounds0[minIdx]){bounds0[minIdx] = bounds1[minIdx];}if (bounds1[maxIdx] > bounds0[maxIdx]){bounds0[maxIdx] = bounds1[maxIdx];}}vtkSmartPointer< vtkPointLocator > pointMerger = vtkSmartPointer< vtkPointLocator >::New();pointMerger->SetTolerance(1e-6);pointMerger->InitPointInsertion(outputIntersection->GetPoints(), bounds0);impl->PointMerger = pointMerger;// This performs the triangle intersection searchobbTree0->IntersectWithOBBTree(obbTree1, 0, vtkIntersectionPolyDataFilter::Impl::FindTriangleIntersections,impl);// Split the first output if so desiredif ( this->SplitFirstOutput ){// 默认执行mesh0->BuildLinks();impl->SplitMesh(0, outputPolyData0, outputIntersection);}else{outputPolyData0->ShallowCopy( mesh0 );}// Split the second output if desiredif ( this->SplitSecondOutput ){mesh1->BuildLinks();impl->SplitMesh(1, outputPolyData1, outputIntersection);}else{outputPolyData1->ShallowCopy( mesh1 );}impl->PointCellIds[0]->Delete();impl->PointCellIds[1]->Delete();delete impl;return 1;}

分析：

• 对输入的两个面数据 mesh0 + mesh1 分别建立一个 vtkOBBTree
• 两个OBB树调用 IntersectWithOBBTree() ， 遍历两个树的每个节点，并通过调用  FindTriangleIntersections() 函数找所有的相交线,作为第一个输出
• 每个面数据调用 SplitMesh() 函数 将面重新三角化,作为第二三个输出，通过调用 intersectionPolyDataFilter->SetSplitFirstOutput(0); 来决定是否重新三角化，如果用不到重新三角化的数据，可以将 off 掉这个功能，加速。

对于输入面的重新三角化的结果：

图中蓝色为相交面，绿色为源数据信息，红色为重新三角化后的输出数据。

【内部分析】

第一个有趣的地方：

typedef struct _CellEdgeLine {vtkIdType CellId;vtkIdType EdgeId;vtkIdType LineId;} CellEdgeLineType;typedef std::multimap< vtkIdType, CellEdgeLineType > PointEdgeMapType;typedef PointEdgeMapType::iterator                   PointEdgeMapIteratorType;

根据之后代码的分析，我才发现，这个地方主要是为了实现 multimap 的遍历：

        PointEdgeMapType    *PointEdgeMap[2];PointEdgeMapIteratorType iterLower =this->PointEdgeMap[index]->lower_bound( ptId );PointEdgeMapIteratorType iterUpper =this->PointEdgeMap[index]->upper_bound( ptId );while (iterLower != iterUpper){if ( iterLower->second.CellId == cellId ){return;}++iterLower;}

（这样子看，觉得这个挺好的）

第二个有趣的地方：

类的定义出就写着 Private implementation to hide STL.

class vtkIntersectionPolyDataFilter::Impl{public:Impl();virtual ~Impl();static int FindTriangleIntersections(vtkOBBNode *node0, vtkOBBNode *node1,vtkMatrix4x4 *transform, void *arg);int SplitMesh(int inputIndex, vtkPolyData *output, vtkPolyData *intersectionLines);protected:vtkCellArray* SplitCell(vtkPolyData *input, vtkIdType cellId,vtkIdType *cellPts, IntersectionMapType *map,vtkPolyData *interLines);void AddToPointEdgeMap(int index, vtkIdType ptId, double x[3],vtkPolyData *mesh, vtkIdType cellId,vtkIdType edgeId, vtkIdType lineId,vtkIdType triPts[3]);void SplitIntersectionLines(int inputIndex, vtkPolyData *sourceMesh,vtkPolyData *splitLines);public:vtkPolyData         *Mesh[2];vtkOBBTree          *OBBTree1;// Stores the intersection lines.vtkCellArray        *IntersectionLines;// Cell data that indicates in which cell each intersection// lies. One array for each output surface.vtkIdTypeArray      *CellIds[2];// Map from points to the cells that contain them. Used for point// data interpolation. For points on the edge between two cells, it// does not matter which cell is recorded because the interpolation// will be the same.  One array for each output surface.vtkIdTypeArray      *PointCellIds[2];// Merging filter used to convert intersection lines from "line// soup" to connected polylines.vtkPointLocator     *PointMerger;// Map from cell ID to intersection line.IntersectionMapType *IntersectionMap[2];// Map from point to an edge on which it resides, the ID of the// cell, and the ID of the line.PointEdgeMapType    *PointEdgeMap[2];protected:Impl(const Impl&); // purposely not implementedvoid operator=(const Impl&); // purposely not implemented};

第一个公有函数:

int vtkIntersectionPolyDataFilter::Impl::FindTriangleIntersections(vtkOBBNode *node0, vtkOBBNode *node1,vtkMatrix4x4 *transform, void *arg){vtkIntersectionPolyDataFilter::Impl *info =reinterpret_cast<vtkIntersectionPolyDataFilter::Impl*>(arg);vtkPolyData     *mesh0                = info->Mesh[0];vtkPolyData     *mesh1                = info->Mesh[1];vtkOBBTree      *obbTree1             = info->OBBTree1;vtkCellArray    *intersectionLines    = info->IntersectionLines;<span style="white-space:pre"></span>// 相交线的几何信息vtkIdTypeArray  *intersectionCellIds0 = info->CellIds[0];<span style="white-space:pre"></span>// mesh0中每个相交三角形cell 的 cellID vtkIdTypeArray  *intersectionCellIds1 = info->CellIds[1];<span style="white-space:pre"></span>// mesh1中每个相交三角形cell 的 cellIDvtkPointLocator *pointMerger          = info->PointMerger;<span style="white-space:pre"></span>// 相交线的点的信息int numCells0 = node0->Cells->GetNumberOfIds();int retval = 0;// 遍历节点0中所有的CellIDfor (vtkIdType id0 = 0; id0 < numCells0; id0++){vtkIdType cellId0 = node0->Cells->GetId(id0);int type0 = mesh0->GetCellType(cellId0);if (type0 == VTK_TRIANGLE){vtkIdType npts0, *triPtIds0;mesh0->GetCellPoints(cellId0, npts0, triPtIds0);double triPts0[3][3];// 当前Cell 中所有点（也就是3个点）的坐标for (vtkIdType id = 0; id < npts0; id++){mesh0->GetPoint(triPtIds0[id], triPts0[id]);}// 判断当前得到的三角形是不是与 Tree1 的 node1 上的三角形相交// Returns true if triangle (optionally transformed) intersects node.if (obbTree1->TriangleIntersectsNode(node1, triPts0[0], triPts0[1], triPts0[2], transform)){int numCells1 = node1->Cells->GetNumberOfIds();for (vtkIdType id1 = 0; id1 < numCells1; id1++){vtkIdType cellId1 = node1->Cells->GetId(id1);int type1 = mesh1->GetCellType(cellId1);if (type1 == VTK_TRIANGLE){// See if the two cells actually intersect. If they do,// add an entry into the intersection maps and add an// intersection line.vtkIdType npts1, *triPtIds1;mesh1->GetCellPoints(cellId1, npts1, triPtIds1);double triPts1[3][3];     // tree1 中 node1 所包含的 Cells 中所有点（也就是3个点）的坐标for (vtkIdType id = 0; id < npts1; id++){mesh1->GetPoint(triPtIds1[id], triPts1[id]);}// 判断是否两个三角形相交int coplanar = 0;double outpt0[3], outpt1[3];int intersects =vtkIntersectionPolyDataFilter::TriangleTriangleIntersection(triPts0[0], triPts0[1], triPts0[2],triPts1[0], triPts1[1], triPts1[2],coplanar, outpt0, outpt1);// 三角形共面 : coplanar = 1if ( coplanar ){// Coplanar triangle intersection is not handled.// This intersection will not be included in the output.intersects = 0;continue;}// 两个三角形相交线，并且线段有长度不相同if ( intersects &&( outpt0[0] != outpt1[0] ||outpt0[1] != outpt1[1] ||outpt0[2] != outpt1[2] ) ){// 存储相交相交线段的信息vtkIdType lineId = intersectionLines->GetNumberOfCells();intersectionLines->InsertNextCell(2);vtkIdType ptId0, ptId1;pointMerger->InsertUniquePoint(outpt0, ptId0);pointMerger->InsertUniquePoint(outpt1, ptId1);intersectionLines->InsertCellPoint(ptId0);intersectionLines->InsertCellPoint(ptId1);intersectionCellIds0->InsertNextValue(cellId0);// 当前 Tree0.node0 对应的 cellIDintersectionCellIds1->InsertNextValue(cellId1);     // 当前 Tree1.node1 对应的 cellID// Map from points to the cells that contain them, 每个点都是临近两个cellinfo->PointCellIds[0]->InsertValue( ptId0, cellId0 );// Tree0.node0.cell0所对应的点info->PointCellIds[0]->InsertValue( ptId1, cellId0 );info->PointCellIds[1]->InsertValue( ptId0, cellId1 );       // Tree1.node1.cell1所对应的点info->PointCellIds[1]->InsertValue( ptId1, cellId1 );// IntersectionMap[0]: Tree0.node0.cell0 所对应的相交线的IDinfo->IntersectionMap[0]->insert(std::make_pair(cellId0, lineId));info->IntersectionMap[1]->insert(std::make_pair(cellId1, lineId));// Check which edges of cellId0 and cellId1 outpt0 and// outpt1 are on, if any.for (vtkIdType edgeId = 0; edgeId < 3; edgeId++){info->AddToPointEdgeMap(0, ptId0, outpt0, mesh0, cellId0,edgeId, lineId, triPtIds0);info->AddToPointEdgeMap(0, ptId1, outpt1, mesh0, cellId0,edgeId, lineId, triPtIds0);info->AddToPointEdgeMap(1, ptId0, outpt0, mesh1, cellId1,edgeId, lineId, triPtIds1);info->AddToPointEdgeMap(1, ptId1, outpt1, mesh1, cellId1,edgeId, lineId, triPtIds1);}retval++;}}}}}}return retval;}

遍历两个树中某两个节点中的每个 cell ，调用函数 TriangleTriangleIntersection 求出两个三角形相交的线，相交线段的两个端点 outpt0, outpt1

我自己觉得比较麻烦的地方就是为了 vtkIntersectionPolyDataFilter::Impl *info 这个变量各种成员属性的赋值，可能是最近变笨了，所以对于有很多属性的类有点迷茫，会忘记、混淆。。。

对于vtkIdTypeArray的操作也不熟

// Insert data at a specified position in the array. void vtkIdTypeArray::InsertValue (vtkIdType id, vtkIdType f ) // Insert data at the end of the array. Return its location in the array. vtkIdType vtkIdTypeArray::InsertNextValue (vtkIdType f)