6月的太阳炙烤大地，在办公室编写代码必须有空调的辅助。旁边的UI同事还买了一杯冰镇饮料，羡慕得我赶紧去接了一杯温开水解解渴。
几天前，我接到新的任务。用户有时候会在我们的软件中创建新的牙齿模型替换旧的模型，新的牙齿包含凸起，这凸起正对应着旧牙齿上的附件。我们需要磨平凸起。
为此，设想的大体步骤是这样的：先找到新牙和附件的交线，然后根据交线找出线内的区域，再去除区域内的cell，最后缝补造成的“洞”。
求交线并不是一件复杂的事情，vtkIntersectionPolyDataFilter 可以帮助我们完成。获取交线之后，还需要把交线内的牙齿区域找出来，这里用到了vtkSelectPolyData。它可以将交线内的points和交线外的points通过标记不同的scalar区分开来。scalar的值小于0的话，则表明属于交线内的点，也就是我们要找的凸起。
demo如下：
#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 <vtkSTLReader.h>
#include <vtkXMLPolyDataReader.h>
#include <vtkProperty.h>
#include <vtkSelectPolyData.h>
using namespace std;
#define vtkSPtr vtkSmartPointer
#define vtkSPtrNew(Var, Type) vtkSPtr<Type> Var = vtkSPtr<Type>::New();
int main()
{
vtkSPtrNew( toothReader, vtkSTLReader );
toothReader->SetFileName( "../newData.STL" );
toothReader->Update();
vtkSPtrNew( circleReader, vtkXMLPolyDataReader );
circleReader->SetFileName( "../circle.vtp" );
circleReader->Update();
vtkSPtrNew( circleMapper, vtkPolyDataMapper );
circleMapper->SetInputData( circleReader->GetOutput() );
vtkSPtrNew( circleActor, vtkActor );
circleActor->SetMapper( circleMapper );
circleActor->GetProperty()->SetColor( 1, 0, 0 );
vtkSmartPointer<vtkSelectPolyData> loop =
vtkSmartPointer<vtkSelectPolyData>::New();
loop->GlobalWarningDisplayOn();
loop->DebugOn();
loop->SetInputData( toothReader->GetOutput() );
loop->SetLoop( circleReader->GetOutput()->GetPoints() );
loop->GenerateSelectionScalarsOn();
loop->SetSelectionModeToSmallestRegion(); //negative scalars inside
loop->Update();
vtkSPtrNew( toothMapper, vtkPolyDataMapper );
toothMapper->SetInputData( loop->GetOutput() );
vtkSPtrNew( toothActor, vtkActor );
toothActor->SetMapper( toothMapper );
vtkSPtrNew( renderer, vtkRenderer );
renderer->AddActor( toothActor );
renderer->AddActor( circleActor );
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;
}
通过阅读源码我们知道，可以通过如下方式获取points的scalar信息。
有了points的scalar信息，就可以给不同的cell赋予不同的scalar值，这对于后面的删除cell和补洞操作意义重大。没有vtkSelectPolyData计算出的points scalar就没有后面的一切。
auto *ptData = loop->GetOutput()->GetPointData();
vtkDataArray *dataArray = ptData->GetAttribute( vtkDataSetAttributes::SCALARS ); 接着是删除凸起，我们可以不更改points，移除相关的cell即可。
类似于这样：
for( int i = 0; i < cellScalars->GetNumberOfTuples(); ++i )
{
if( fabs( TARGET_TYPE - cellScalars->GetTuple1( i ) ) < 1e-6 )
{
newData->DeleteCell( i );
}
}
newData->RemoveDeletedCells();最后是缝补开口，这一部分内容比较复杂。可以简单地认为是这样完成的：
我们寻找“洞”的边界，然后查找夹角最小的角，这个角对应的顶点和另外两个点构造新的三角形，一直这样迭代下去，直到世界的尽头，哦不，直到边界上的点没了，完成缝补。
我以为事情就这样结束了，然后顺利的提交了代码，来一杯枸杞，享受一下难得的成就感。
然而，测试同事反馈的问题颠覆了我的三观。
在完成一次补洞之后，如果想要再来delete cell（假设第一次没有删除干净），vtkSelectPolyData返回的数据是0个point，0个cell。第二次计算的交线是没有问题的，它就在新的模型上。
我打开了vtk的debug窗口。
loop->GlobalWarningDisplayOn();
loop->DebugOn();发现，第二次工作的时候，总会有错误出现：
ERROR: In D:\UlabVTK\VTK-8.2.0\Filters\Modeling\vtkSelectPolyData.cxx, line 263
vtkSelectPolyData (000001A5FEC7FD80): Can't follow edge看了对应的源码，好神奇。在使用vtkSelectPolyData之前，我们已经完成了links和cells的build，我无法理解为什么会找不到point的邻间关系。
然后，我在第二次尝试计算交线内的cell时提前将牙齿和交线分别保存成STL和VTP文件，我想用上面的demo去查看细节，比如我们可以很方便的查看模型的网状结构，或者用vtkActor2D显示一些内部的points信息。
瘆人的结果出现了，在demo程序里居然没有出现任何错误，交线外的大部分points对应的scalar是正数，交线内的points scalar是负数。
为什么同样的数据会有不同的结果？两个程序使用VTK版本以及QT环境一模一样，就连文本编辑器、IDE都是一样的呀。
难道我遇到了量子物理的观察者效应？别开玩笑了，快code freeze了，波尔大神，我玩不起啊。
在不断怀疑的过程中，我将注意力集中在了这个问题上：“这两个程序使用的数据真的是一样的吗？” 通过对比，我发现demo使用的数据是从文件中读出来的，而工程中使用的数据是内存中的，从逻辑上讲，文件中记录的信息就是内存中的信息，两者没有任何区别。但是在这里，我不得不猜测，它们就是不一样的！
于是，我在工程中也做了类似的事情，在每次使用vtkSelectPolyData处理数据之前，我都先将牙齿模型存成STL然后从文件中读取出来，然后再将读取出来的信息作为输入数据，这看起来就是画蛇添足。
结果出来了，工程中的vtkSelectPolyData也开始正常工作！这是为什么，因为STL的读写操作影响了polydata吗？单纯看“Can’t follow edge”的错误提示，我没能得到有价值的线索。不管怎么说，在使用该对象处理数据之前居然需要写读一遍信息本身就显得很诡异，或者，这就是VTK8.2上的一个bug？也许，当我升级VTK之后才能找到答案。
相关文章：
VTK – Sort The Points On The Circle
VTK – Expand Sample Area From A Cell
[…] can help us to avoid error Can't follow edge in vtkSelectPolyData as described in the article Is there a bug in vtkSelectPolyData?. The data after I/O operations in STL file format, it can be handled correctly by […]
I have issues with “Can’t follow edge” too.
Could it be that you wrote the PolyData to an STL and read it back? This way the mesh get’s triangulated, as STL can only store triangles.
Maybe your hole filling algorithm could have left n-gons which makes it difficult for the select to find a path between points of your PolyLine.
Yes, I handled the mesh before using vtkSelectPolyData, after that there are a few edge points that we can’t find a valid edge for them.
I also found that my old PolyData has many invalid points that are not used in cells, but vtkSelectPolyData takes all points in the data. The frequency of the error occurrence became smaller after I removed the points which are not used in cells. Relevant link: https://www.weiy.city/2020/07/handle-invalid-points-in-vtkstlwriter-and-vtkxmlpolydata