求两个多边形数据 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 ，所以这个类中最为重要的函数就是 

int	RequestData (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 mesh1
vtkSmartPointer< 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 search
obbTree0->IntersectWithOBBTree
(obbTree1, 0, vtkIntersectionPolyDataFilter::Impl::FindTriangleIntersections,
impl);

// Split the first output if so desired
if ( this->SplitFirstOutput )
{	// 默认执行
mesh0->BuildLinks();
impl->SplitMesh(0, outputPolyData0, outputIntersection);
}
else
{
outputPolyData0->ShallowCopy( mesh0 );
}

// Split the second output if desired
if ( 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 implemented
void 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 的 cellID
vtkPointLocator *pointMerger          = info->PointMerger;<span style="white-space:pre">		</span>// 相交线的点的信息

int numCells0 = node0->Cells->GetNumberOfIds();
int retval = 0;
// 遍历节点0中所有的CellID
for (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 = 1
if ( 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 对应的 cellID
intersectionCellIds1->InsertNextValue(cellId1);     // 当前 Tree1.node1 对应的 cellID
// Map from points to the cells that contain them, 每个点都是临近两个cell
info->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 所对应的相交线的ID
info->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)