VTK知识学习（42）-基本的图形操作（三）

1、网格平滑

1）概述

现代扫描技术的发展使得获取点云数据不再困难，通过曲面重建技术可以获取表面网格来表示各种复杂的实体。但是点云数据中往往存在噪声，这样得到的重建网格通常都需要进行平滑处理。
拉普拉斯平滑是一种常用的网格平滑算法。该方法的原理比较简单，如图所示，将每个点用其邻域点的中心来代替。通过不断地迭代，可以得到较为光滑的网格。

VTK中的 vtkSmoothPolyDataFilter 类实现了网格的拉普拉斯平滑算法。

2）代码

private void TestLasianSmooth()
{vtkPolyDataReader reader = vtkPolyDataReader.New();reader.SetFileName("F:\\code\\VTK\\TestActiViz\\data\\fran_cut.vtk");reader.Update();//拉普拉斯平滑vtkSmoothPolyDataFilter smoothFilter = vtkSmoothPolyDataFilter.New();smoothFilter.SetInputConnection(reader.GetOutputPort());smoothFilter.SetNumberOfIterations(200);   //平滑次数 次数越大平滑越厉害smoothFilter.BoundarySmoothingOn();    //控制对边界点平滑开启smoothFilter.Update();vtkWindowedSincPolyDataFilter sincFilter = vtkWindowedSincPolyDataFilter.New();sincFilter.SetInputConnection(reader.GetOutputPort());sincFilter.SetNumberOfIterations(200);sincFilter.Update();ShowImageBy3(reader.GetOutput(), smoothFilter.GetOutput(), sincFilter.GetOutput());
}

3）效果

4）说明

vtkSmoothPolyDataFilter::SetNumberOfterations()控制平滑次数，次数越大平滑越厉害。平滑结果如图所示，左图为原始模型，右图为平滑后的模型。从图中可以看出，在经过200次拉普拉斯平滑后，模型变得非常光滑，但是在平滑的同时也损失了一些细节信息，例如在眼窝处的细节已经平滑掉，在使用该方法时需要注意这点。

该类中还有多个变量来控制平滑过程,利用这些变量在一定程度上可以控制细节的损失。BoundarySmoothing控制是否对边界点平滑。这里需要理解边界点的概念:如果一个网格模型中，一条边只被一个单元包含，那么这条边就是边界边，而边界边上的点则为边界点。如果一个模型中含有边界边，则说明该模型不是封闭的，如本例中的模型。
FeatureEdgeSmoothing控制是否对特征边上的点的平滑。如果一条边被两个邻近的多边形共用，若这两个多边形法向量的夹角(特征角)大于定义的值，则说明该边为一条特征边。因此，FeatureEdgeSmoothing设置开始时，需要调用SetFeatureAngle()函数设置特征角的阈值。

如图所示，特征边e的两个相邻面片的特征角为θ，特征角值越大，说明该边越尖锐。特征边/点往往表示模型的细节,平滑过程中不进行处理，以保持细节。

虽然通过特征边平滑设置可以降低一部分细节损失，但并不能完全避免，且随着拉普拉斯平滑的不断迭代，模型会逐渐向网格的中心收缩。所以，视窗多数据过滤器(详见参考文献[11)是一个更好的选择，该算法使用窗口辛克函数实现网格平滑，能够最小程度地避免收缩，其使用方法与vtkSmooth Polydata Filter相同。

2、封闭性检测

1）概述

由于受原始数据、重建方法的限制，得到的网格模型并不是封闭的。有时为了显示或者处理某些要求，需要网格必须是封闭的。封闭性网格应该比较好理解，比如一个球面网格。上节内容中提到了边界边的概念:如果一条边只被一个多边形包含，那么这条边就是边界边。是否存在边界边是检测一个网格模型是否封闭的重要特征。

vtkFeatureEdges是一个非常重要的类，该类能够提取多边形网格模型中四种类型的边。

①边界边。即只被一个多边形或者一条边包含的边。
②非流形边。被3个或者3个以上的多边形包含的边即为非流形边。
③特征边。上节中也提到，需要设置一个特征角的阈值，当包含同一条边的两个三角形的法向量的夹角大于该阀值时，即为一个特征边。
④流形边。只被两个多边形包含的边即为流形边。

可以使用该类来检测是否存在边界边，并依此来判断网格是否封闭。

2）代码

void GenerateData(vtkSmartPointer<vtkPolyData> input)
{vtkSmartPointer<vtkSphereSource> sphereSource =vtkSmartPointer<vtkSphereSource>::New();sphereSource->Update();vtkSmartPointer<vtkIdTypeArray> ids =vtkSmartPointer<vtkIdTypeArray>::New();ids->SetNumberOfComponents(1);ids->InsertNextValue(2);ids->InsertNextValue(10);vtkSmartPointer<vtkSelectionNode> selectionNode =vtkSmartPointer<vtkSelectionNode>::New();selectionNode->SetFieldType(vtkSelectionNode::CELL);selectionNode->SetContentType(vtkSelectionNode::INDICES);selectionNode->SetSelectionList(ids);selectionNode->GetProperties()->Set(vtkSelectionNode::INVERSE(), 1);vtkSmartPointer<vtkSelection> selection =vtkSmartPointer<vtkSelection>::New();selection->AddNode(selectionNode);vtkSmartPointer<vtkExtractSelection> extractSelection =vtkSmartPointer<vtkExtractSelection>::New();extractSelection->SetInputData(0, sphereSource->GetOutput());extractSelection->SetInputData(1, selection);extractSelection->Update();vtkSmartPointer<vtkDataSetSurfaceFilter> surfaceFilter =vtkSmartPointer<vtkDataSetSurfaceFilter>::New();surfaceFilter->SetInputConnection(extractSelection->GetOutputPort());surfaceFilter->Update();input->ShallowCopy(surfaceFilter->GetOutput());
}void  TestDataClosed()
{vtkSmartPointer<vtkPolyData> input = vtkSmartPointer<vtkPolyData>::New();GenerateData(input);vtkSmartPointer<vtkFeatureEdges> featureEdges =vtkSmartPointer<vtkFeatureEdges>::New();featureEdges->SetInputData(input);featureEdges->BoundaryEdgesOn();featureEdges->FeatureEdgesOff();featureEdges->ManifoldEdgesOff();featureEdges->NonManifoldEdgesOff();featureEdges->Update();int numberOfOpenEdges = featureEdges->GetOutput()->GetNumberOfCells();if (numberOfOpenEdges){std::cout << "该网格模型不是封闭的..." << std::endl;}else{std::cout << "该网格模型是封闭的..." << std::endl;return;}vtkSmartPointer<vtkFillHolesFilter> fillHolesFilter =vtkSmartPointer<vtkFillHolesFilter>::New();fillHolesFilter->SetInputData(input);fillHolesFilter->Update();vtkSmartPointer<vtkPolyDataNormals> normals =vtkSmartPointer<vtkPolyDataNormals>::New();normals->SetInputConnection(fillHolesFilter->GetOutputPort());normals->ConsistencyOn();normals->SplittingOff();normals->Update();//double leftViewport[4] = { 0.0, 0.0, 0.5, 1.0 };double rightViewport[4] = { 0.5, 0.0, 1.0, 1.0 };vtkSmartPointer<vtkPolyDataMapper> originalMapper =vtkSmartPointer<vtkPolyDataMapper>::New();originalMapper->SetInputData(input);vtkSmartPointer<vtkProperty> backfaceProp =vtkSmartPointer<vtkProperty>::New();backfaceProp->SetDiffuseColor(0.89, 0.81, 0.34);vtkSmartPointer<vtkActor> originalActor =vtkSmartPointer<vtkActor>::New();originalActor->SetMapper(originalMapper);originalActor->SetBackfaceProperty(backfaceProp);originalActor->GetProperty()->SetDiffuseColor(1.0, 0.3882, 0.2784);vtkSmartPointer<vtkPolyDataMapper> edgeMapper =vtkSmartPointer<vtkPolyDataMapper>::New();edgeMapper->SetInputData(featureEdges->GetOutput());vtkSmartPointer<vtkActor> edgeActor =vtkSmartPointer<vtkActor>::New();edgeActor->SetMapper(edgeMapper);edgeActor->GetProperty()->SetEdgeColor(0., 0., 1.0);edgeActor->GetProperty()->SetEdgeVisibility(1);edgeActor->GetProperty()->SetLineWidth(5);vtkSmartPointer<vtkPolyDataMapper> filledMapper =vtkSmartPointer<vtkPolyDataMapper>::New();filledMapper->SetInputData(normals->GetOutput());vtkSmartPointer<vtkActor> filledActor =vtkSmartPointer<vtkActor>::New();filledActor->SetMapper(filledMapper);filledActor->GetProperty()->SetDiffuseColor(1.0, 0.3882, 0.2784);vtkSmartPointer<vtkRenderer> leftRenderer =vtkSmartPointer<vtkRenderer>::New();leftRenderer->SetViewport(leftViewport);leftRenderer->AddActor(originalActor);leftRenderer->AddActor(edgeActor);leftRenderer->SetBackground(1.0, 1.0, 1.0);vtkSmartPointer<vtkRenderer> rightRenderer =vtkSmartPointer<vtkRenderer>::New();rightRenderer->SetViewport(rightViewport);rightRenderer->AddActor(filledActor);rightRenderer->SetBackground(1.0, 1.0, 1.0);vtkSmartPointer<vtkRenderWindow> renderWindow =vtkSmartPointer<vtkRenderWindow>::New();renderWindow->AddRenderer(leftRenderer);renderWindow->AddRenderer(rightRenderer);renderWindow->SetSize(640, 320);renderWindow->Render();renderWindow->SetWindowName("PolyDataClosed");vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =vtkSmartPointer<vtkRenderWindowInteractor>::New();renderWindowInteractor->SetRenderWindow(renderWindow);leftRenderer->GetActiveCamera()->SetPosition(0, -1, 0);leftRenderer->GetActiveCamera()->SetFocalPoint(0, 0, 0);leftRenderer->GetActiveCamera()->SetViewUp(0, 0, 1);leftRenderer->GetActiveCamera()->Azimuth(30);leftRenderer->GetActiveCamera()->Elevation(30);leftRenderer->ResetCamera();rightRenderer->SetActiveCamera(leftRenderer->GetActiveCamera());renderWindow->Render();renderWindowInteractor->Start();
}

由于没有查找一个AddNode函数所以这里是c++的代码实现。

3）效果

左图为原始模型，并将检测的边界边用蓝色显示，右图为漏洞填补后的结果。

4）说明

        这里建立了一个球面网格，并去除了其中两个三角面片(单元)。将该数据作为 vtkFeatureEdges 的输入，BoundaryEdgesOn()函数设置提取边界边，而其他三种类型的边则不予考虑。执行完毕，其输出 GetOutput()为一个包含边信息的 vkPolyData 数据。可以通过判断边界边的数目来确定网格是否封闭:
        int numberOfOpenEdges = featureEdges->GetOutput()->GetNumberOfCellsO;
        当网格为非封闭时，可以为检测结果建立相应的vtkPolyDataMapper和 vtkActor 对象，将边界边与原网格同时显示，以观察检测结果是否正确。

仅仅检测出是否封闭是不够的，很多情况下，还需将这些漏洞填补起来。VTK中有现成的类来完成这个功能--vtkFilHolesFilter。其内部执行过程是首先检测出网格中的所有边界边，然后找出这些边界边中的每一个闭合回路，最后将这些闭合回路进行三角化(即生成三角网格)以实现填补的目的。该类的使用非常简单，只需要设置需要填补的网格数据即可:

vtkSmartPointer<vtkFillHolesFilter> fillHolesFilter =vtkSmartPointer<vtkFillHolesFilter>::New();fillHolesFilter->SetInputData(input);fillHolesFilter->Update();vtkSmartPointer<vtkPolyDataNormals> normals =vtkSmartPointer<vtkPolyDataNormals>::New();normals->SetInputConnection(fillHolesFilter->GetOutputPort());normals->ConsistencyOn();normals->SplittingOff();normals->Update();

需要注意的是，有些边界的闭合回路是不需要三角化的，例如一个平面网格，若填补其四周的边界边，则会与原网格产生覆盖。vtkFillHolesFilter 中的 SetHoleSize()函数可用于控制需要修补的漏洞的面积最大值，大于该值的漏洞则不需要填补处理。

vtkPolyDataNormals这个类，主要是用来计算法向量的，而法向量则与光照和阴影计算相关。单元的法向量朝向则与单元的点顺序相关，只有保持所有单元的点顺序一致才能得到正确的法向量，否则在网格模型显示时会得到意外的结果。由于经过漏洞填充，模型的所有单元的点顺序并不一致，因此使用vtkPolyDataNormal::ConsistencyOn()进行调整。

3、连通区域分析

1）概述

许多图形数据中，并非只包含一个对象(连通区域)。而在处理这些图形数据时，有时需要对每一个对象单独处理或者让其单独显示。比如利用MarchingCube 方法提取三维图像中的等值面，得到的结果往往是存在多个连通的对象区域，这时就需要对图形数据做连通区域分析，提取每个连通区域并计算其属性信息，以此来得到需要的连通区域。

2）代码

 private void TestCompExtract(){//球面数据vtkSphereSource sphereSource = vtkSphereSource.New();sphereSource.SetRadius(10);sphereSource.SetThetaResolution(10);sphereSource.SetPhiResolution(10);sphereSource.Update();//锥体数据vtkConeSource coneSource = vtkConeSource.New();coneSource.SetRadius(5);coneSource.SetHeight(10);coneSource.SetCenter(25, 0, 0);  //锥体中心设置，防止两个数据之间存在覆盖coneSource.Update();//构造一个含有多个连通区域的模型数据vtkAppendPolyData appendFilter = vtkAppendPolyData.New();appendFilter.AddInputData(sphereSource.GetOutput());appendFilter.AddInputData(coneSource.GetOutput());appendFilter.Update();//实现连通区域分析vtkPolyDataConnectivityFilter connectivityFilter = vtkPolyDataConnectivityFilter.New();connectivityFilter.SetInputData(appendFilter.GetOutput());connectivityFilter.SetExtractionModeToLargestRegion();   //提取具有多点的连通区域connectivityFilter.AddSeed(100);connectivityFilter.Update();ShowImageBy2(appendFilter.GetOutput(), connectivityFilter.GetOutput());}

3）效果

4）说明

示例构造一个含有多个连通区域的模型数据。vtkAppendPolyData可以实现vtkPolyData的合并，使用该类可以方便地构造含有多个连通区域的数据，该类接收两个或者多个vtkPolyData 数据输入，合并结果包含输入数据的所有几何和拓扑数据。若输入的两个或者多个数据都含有点属性数据，则将其存储至输出结果中;对于单元属性数据也是如此。

示例先定义了-个球面和一-个锥体数据。为了防止两个数据之间存在覆盖，这里将锥体中心设置到(25,0,0);然后将两个数据设置到 vtkAppendPolyData 的输入，即可生成一个包含两个对象的 vkPolyData 数据,每个输入模型为一个连通区域。VTK 中的 vtkPolyDataConne-ctivityFilter 类可用于实现连通区域分析，该类接受 vtkPolyData数据作为输入，其使用如下:

  //实现连通区域分析vtkPolyDataConnectivityFilter connectivityFilter = vtkPolyDataConnectivityFilter.New();connectivityFilter.SetInputData(appendFilter.GetOutput());connectivityFilter.SetExtractionModeToLargestRegion();   //提取具有多点的连通区域connectivityFilter.AddSeed(100);connectivityFilter.Update();

SetExtractionModeToLargestRegion()函数用于提取具有最多点的连通区域，因此在该例中
得到的结果为球面数据。

除了能够提取最大连通区域，vtkPolyDataConnectivityFilter 类还支持以下5种模式

1)SetExtractionModeToAllRegions()。该模式主要用于连通区域标记，配合函数 ColorRegionsOn()使用，在连通区域检测的同时，生成一个名为Regionld 的点属性数据。

2)SetExtractionModeToSpecifiedRegions()。该模式用于提取一个或者多个连通区域。在该模式下，需要通过 AddSpecifiedRegion()来添加需要提取的区域号，区域号从0开始。

3)SetExtractionModeToClosestPointRegion()。该模式下需要使用 SetClosestPoint()函数设置一个空间点坐标，其执行结果为离该点最近的连通区域。

4)SetExtractionModeToPointSeededRegions()。该模式下需要使用 AddSeed()函数添加种子点，提取种子点所在的区域。

5)SetExtractionModeToCellSeededRegions()。该模式下需要使用 AddSeedO函数添加种子单元，提取种子单元所在的区域。

5）其它

许多情况下，事先并不知道整个数据中存在多少个区域，因此需要先获取连通区域的数目。

private void TestConnectedAllCompExtract()
{//球面数据vtkSphereSource sphereSource = vtkSphereSource.New();sphereSource.SetRadius(10);sphereSource.SetThetaResolution(10);sphereSource.SetPhiResolution(10);sphereSource.Update();//锥体数据vtkConeSource coneSource = vtkConeSource.New();coneSource.SetRadius(5);coneSource.SetHeight(10);coneSource.SetCenter(25, 0, 0);  //锥体中心设置，防止两个数据之间存在覆盖coneSource.Update();//构造一个含有多个连通区域的模型数据vtkAppendPolyData appendFilter = vtkAppendPolyData.New();appendFilter.AddInputData(sphereSource.GetOutput());appendFilter.AddInputData(coneSource.GetOutput());appendFilter.Update();//实现连通区域分析vtkPolyDataConnectivityFilter connectivityFilter = vtkPolyDataConnectivityFilter.New();connectivityFilter.SetInputData(appendFilter.GetOutput());//该模式用于提取一个或者多个连通区域。//在该模式下，需要通过 AddSpecifiedRegion()来添加需要提取的区域号，区域号从0开始。connectivityFilter.SetExtractionModeToLargestRegion();   connectivityFilter.Update();int regionNum = connectivityFilter.GetNumberOfExtractedRegions();for(int i=0;i<regionNum; i++){vtkPolyDataConnectivityFilter connectivityFilter2 = vtkPolyDataConnectivityFilter.New();connectivityFilter2.SetInputData(appendFilter.GetOutput());connectivityFilter2.InitializeSpecifiedRegionList();connectivityFilter2.SetExtractionModeToSpecifiedRegions();connectivityFilter2.AddSpecifiedRegion(i);connectivityFilter2.Update();vtkPolyDataWriter writer = vtkPolyDataWriter.New();writer.SetFileName($"{i:000}.vtk");writer.SetInputData(connectivityFilter2.GetOutput());writer.Update();}ShowImageBy2(appendFilter.GetOutput(), connectivityFilter.GetOutput());
}

示例在 SetExtractionModeToAllRegions()模式下进行连通区域分析，并通过 GetNumberOfExtractedRegions()函数获取连通区域的数目;然后通过一个循环，在SetExtractionModeToSpecifiedRegions()模式下,每次循环通过 AddSpecifedRegion()函数设置提取的连通区域号并单独保存为一个.vtk 文件。InitializeSpecifiedRegionList()函数用于清空要提取的连通区域号的列表。