第五节《VTK 点与网格模型求交处理技巧》

116 阅读 0 评论 77 点赞

我是靠谱客的博主现代航空，这篇文章主要介绍第五节《VTK 点与网格模型求交处理技巧》，现在分享给大家，希望可以做个参考。

使用VTK做项目常常需要用到单点或点集与网格模型求交计算，常用手段是使用vtkOBBTree进行求交计算。我基于VTK提供的算法总结了3种点与网格求交计算的方法，第一种vtkOBBTree直线与网格求交计算；第二种局部网格直线求交计算；第三种基于深度缓存求交计算；

第一种vtkOBBTree直线与网格求交计算，vtkOBBTree是用于生成模型OBB树的对象数据结构，它基于模型创建一个有向包围盒，对包围盒空间划分了更多层次的空间，直线与模型相交能够快速锁定某个区域，进行求交计算，用法如下：

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
// polyData 待求交计算的网格模型
vtkSmartPointer<vtkOBBTree> obbTree = vtkSmartPointer<vtkOBBTree>::New();
obbTree->SetDataSet(polyData);
obbTree->BuildLocator();
double p1[3];
double p2[3];
vtkSmartPointer<vtkPoints> intersectPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkIdList> intersectCells = vtkSmartPointer<vtkIdList>::New();
obbTree->IntersectWithLine(p1, p2, intersectPoints, intersectCells);
参数p1 线段起点
参数p2 线段终点
参数intersectPoints 输出交点
参数intersectCells 输出交到的三角形

vtk提供的接口是一条线段与模型求交，输出结果为线段内所有的交点和所有的网格单元，实际应用要比这个复杂一些，大部分应用我们希望输入是一个点p和一个向量n来表示一条直线，输出有多种情况，1.输出距离p点最近的交点；2.输出投影点前方并且距离p点最近的交点；3.输出投影点前方并且距离p点最近，并且交到的网格单元正面；4.输出投影点前方并且距离p点最近，并且交到的网格单元背面。下面我展示一个函数来完成上述的过程：此函数所在位置是FreePolyData类里面，使用前需要执行buildOBBTree（）；

复制代码

enum MAP_TYPE
{
DIST_MIN,
// 距离最近的交点
DIST_MIN_LINEFRONT,
// 投影点前方，距离最近的点
DIST_MIN_LINEFRONT_CELLFRONT,
// 投影点前方，距离最近的点，交到三角形正面
DIST_MIN_LINEFRONT_CELLBACK
// 投影点前方，距离最近的点，交到三角形背面
};
int FreePolyData::pointMappingToPolyData(const FreeVector &p, const FreeVector &n, MAP_TYPE mapType, FreeVector &outPut)
{
FreeVector p1 = p;
FreeVector p2 = p;
p1.move(n, -100000);
p2.move(n, 100000);
vtkSmartPointer<vtkPoints> intersectPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkIdList> intersectCells = vtkSmartPointer<vtkIdList>::New();
if(obbTree == NULL)
{
buildOBBTree();
}
obbTree->IntersectWithLine(p1.GetData(), p2.GetData(), intersectPoints, intersectCells);
int outPutID = -1;
double dist = 10000000;
for(int i = 0; i < intersectPoints->GetNumberOfPoints(); i++)
{
double resP[3];
intersectPoints->GetPoint(i, resP);
int cellID = intersectCells->GetId(i);
FreeVector theP(resP);
if(mapType == DIST_MIN) // 所有交点
{
double d = theP.point3Distance(p);
if(d < dist)
{
dist = d;
outPut = theP;
outPutID = cellID;
}
} else if(mapType == DIST_MIN_LINEFRONT) { // 投影点前方的
FreeVector n1 = p2 - p;
FreeVector n2 = theP - p;
if(n1*n2 > 0)
{
double d = theP.point3Distance(p);
if(d < dist)
{
dist = d;
outPut = theP;
outPutID = cellID;
}
}
} else if(mapType == DIST_MIN_LINEFRONT_CELLFRONT){
FreeVector n1 = p2 - p;
FreeVector n2 = theP - p;
FreeVector cellNormal = getCellNormal(cellID);
if(n1*n2 > 0 && cellNormal*n2 < 0)
{
double d = theP.point3Distance(p);
if(d < dist)
{
dist = d;
outPut = theP;
outPutID = cellID;
}
}
} else if(mapType == DIST_MIN_LINEFRONT_CELLFRONT) {
FreeVector n1 = p2 - p;
FreeVector n2 = theP - p;
FreeVector cellNormal = getCellNormal(cellID);
if(n1*n2 > 0 && cellNormal*n2 > 0)
{
double d = theP.point3Distance(p);
if(d < dist)
{
dist = d;
outPut = theP;
outPutID = cellID;
}
}
}
}
return outPutID;
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88

enum MAP_TYPE
{
DIST_MIN,
// 距离最近的交点
DIST_MIN_LINEFRONT,
// 投影点前方，距离最近的点
DIST_MIN_LINEFRONT_CELLFRONT,
// 投影点前方，距离最近的点，交到三角形正面
DIST_MIN_LINEFRONT_CELLBACK
// 投影点前方，距离最近的点，交到三角形背面
};
int FreePolyData::pointMappingToPolyData(const FreeVector &p, const FreeVector &n, MAP_TYPE mapType, FreeVector &outPut)
{
FreeVector p1 = p;
FreeVector p2 = p;
p1.move(n, -100000);
p2.move(n, 100000);
vtkSmartPointer<vtkPoints> intersectPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkIdList> intersectCells = vtkSmartPointer<vtkIdList>::New();
if(obbTree == NULL)
{
buildOBBTree();
}
obbTree->IntersectWithLine(p1.GetData(), p2.GetData(), intersectPoints, intersectCells);
int outPutID = -1;
double dist = 10000000;
for(int i = 0; i < intersectPoints->GetNumberOfPoints(); i++)
{
double resP[3];
intersectPoints->GetPoint(i, resP);
int cellID = intersectCells->GetId(i);
FreeVector theP(resP);
if(mapType == DIST_MIN) // 所有交点
{
double d = theP.point3Distance(p);
if(d < dist)
{
dist = d;
outPut = theP;
outPutID = cellID;
}
} else if(mapType == DIST_MIN_LINEFRONT) { // 投影点前方的
FreeVector n1 = p2 - p;
FreeVector n2 = theP - p;
if(n1*n2 > 0)
{
double d = theP.point3Distance(p);
if(d < dist)
{
dist = d;
outPut = theP;
outPutID = cellID;
}
}
} else if(mapType == DIST_MIN_LINEFRONT_CELLFRONT){
FreeVector n1 = p2 - p;
FreeVector n2 = theP - p;
FreeVector cellNormal = getCellNormal(cellID);
if(n1*n2 > 0 && cellNormal*n2 < 0)
{
double d = theP.point3Distance(p);
if(d < dist)
{
dist = d;
outPut = theP;
outPutID = cellID;
}
}
} else if(mapType == DIST_MIN_LINEFRONT_CELLFRONT) {
FreeVector n1 = p2 - p;
FreeVector n2 = theP - p;
FreeVector cellNormal = getCellNormal(cellID);
if(n1*n2 > 0 && cellNormal*n2 > 0)
{
double d = theP.point3Distance(p);
if(d < dist)
{
dist = d;
outPut = theP;
outPutID = cellID;
}
}
}
}
return outPutID;
}

第二种局部网格直线求交计算：上述采用vtkOBBTree进行求交计算，vtkOBBTree构建了层级结构，加快了求交计算方法，但有些时候我们的数据量过大时，采用vtkOBBTree计算依然不能满足计算时间，并且我们待计算的点可能在网格附近，此时采用第二种求交计算策略恰当好处。局部网格求交计算讲述的是，网格模型周围的点或点集向网格模型上的投影；

适用场景，例如我们想在模型表面显示一条参数化曲线，曲线控制点在模型表面，当进行拟合计算后，曲线的拟合点会脱离模型表面，此时需要将曲线的拟合点向模型表面进行投影，这样才能保证显示的曲线在模型表面，下面请看局部点或点集投影算法：

复制代码

/**
* @brief FreePolyData::pointMappingToPolyData 点按照n方向与模型的投影
* @param p 点
* @param n 投影方向
* @param range 投影范围
* @param outPut 输出交点
* @return 输出投影到的Cell ID
*/
int FreePolyData::pointMappingToPolyData(const FreeVector &p, const FreeVector &n, double range, FreeVector &outPut)
{
if(kdTree == NULL)
{
buildKdtree();
}
FreeVector p1 = p;
FreeVector p2 = p;
p1.move(n, -100);
p2.move(n, 100);
vtkSmartPointer<vtkIdList> result = vtkSmartPointer<vtkIdList>::New();
kdTree->FindPointsWithinRadius(range, p.GetData(), result);
std::vector<bool> pointMark(polyData->GetNumberOfPoints(), false);
for(int i = 0; i < result->GetNumberOfIds(); i++)
{
pointMark[result->GetId(i)] = true;
}
for(int i = 0; i < polyData->GetNumberOfCells(); i++)
{
vtkCell* cell = polyData->GetCell(i);
int v1 = cell->GetPointId(0);
int v2 = cell->GetPointId(1);
int v3 = cell->GetPointId(2);
if(pointMark[v1]||pointMark[v2]||pointMark[v3])
{
double t;
double intersectionCoordinates[3];
double parametricCoordinates[3];
int subId;
int res = cell->IntersectWithLine(p1.GetData(), p2.GetData(), 0.0001, t, intersectionCoordinates, parametricCoordinates, subId);
if(res == 1)
{
outPut = FreeVector(intersectionCoordinates[0], intersectionCoordinates[1], intersectionCoordinates[2]);
return i;
}
}
}
return -1;
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
/**
* @brief FreePolyData::pointMappingToPolyData 点按照n方向与模型的投影
* @param p 点
* @param n 投影方向
* @param range 投影范围
* @param outPut 输出交点
* @return 输出投影到的Cell ID
*/
int FreePolyData::pointMappingToPolyData(const FreeVector &p, const FreeVector &n, double range, FreeVector &outPut)
{
if(kdTree == NULL)
{
buildKdtree();
}
FreeVector p1 = p;
FreeVector p2 = p;
p1.move(n, -100);
p2.move(n, 100);
vtkSmartPointer<vtkIdList> result = vtkSmartPointer<vtkIdList>::New();
kdTree->FindPointsWithinRadius(range, p.GetData(), result);
std::vector<bool> pointMark(polyData->GetNumberOfPoints(), false);
for(int i = 0; i < result->GetNumberOfIds(); i++)
{
pointMark[result->GetId(i)] = true;
}
for(int i = 0; i < polyData->GetNumberOfCells(); i++)
{
vtkCell* cell = polyData->GetCell(i);
int v1 = cell->GetPointId(0);
int v2 = cell->GetPointId(1);
int v3 = cell->GetPointId(2);
if(pointMark[v1]||pointMark[v2]||pointMark[v3])
{
double t;
double intersectionCoordinates[3];
double parametricCoordinates[3];
int subId;
int res = cell->IntersectWithLine(p1.GetData(), p2.GetData(), 0.0001, t, intersectionCoordinates, parametricCoordinates, subId);
if(res == 1)
{
outPut = FreeVector(intersectionCoordinates[0], intersectionCoordinates[1], intersectionCoordinates[2]);
return i;
}
}
}
return -1;
}

需要使用vtkKdtree寻找附近点集，然后与找到的点集所在的Cell进行求交计算，范围越大计算越慢，范围越小计算越快。此函数所在位置是FreePolyData类里面，使用前需要执行buildKdtree（）；

第三种基于深度缓存求交计算：上述两种方法求交计算还是比较耗时的，不适用大量点云计算的，如果有大量点云需要求交计算应该怎么处理呢，我的想法是想利用vtk中的相机投影矩阵进行计算，将相机视口2D坐标转换到3D坐标，获取深度信息，即可完成投影计算，条件是点集所有点的投影方向是一致的。

1.构建局部坐标系：将投影方向当做Z轴构建坐标系

2.将点和待求交的模型变换到局部坐标系内

3.计算模型包围盒，计算窗口尺寸

4.创建相机视口，并将模型数据天骄到视口中

5.创建基于深度缓存拾取器，计算每个点的交点

6.将结果点集变回原来坐标系，没投上的点至0,0,0

这里函数输出点集与输出点集个数相等，顺序相同，没投上的点是（0,0,0）；

复制代码

/**
* @brief pointsDepthMapping 点集与模型求交
* @param inputPoints 输入点集
* @param mapNormal 投影方向
* @param polyData 模型
* @param outPoints 输出点集，个数与输入对应，没投上是0，0，0
*/
void pointsDepthMapping(std::vector<FreeVector> &inputPoints, FreeVector mapNormal, vtkPolyData *polyData, std::vector<FreeVector> &outPoints)
{
vtkSmartPointer<vtkPolyData> tPolyData = vtkSmartPointer<vtkPolyData>::New();
tPolyData->DeepCopy(polyData);
// 1.构建局部坐标系：将投影方向当做Z轴构建坐标系
FreeVector zDir, xDir, yDir;
zDir = -mapNormal;
zDir.Normalize();
constructCoordinateByZAxle(zDir, xDir, yDir); // 利用Z轴构建一个坐标系
vtkSmartPointer<vtkMatrix4x4> m = vtkSmartPointer<vtkMatrix4x4>::New();
m->SetElement(0, 0, xDir[0]);
m->SetElement(1, 0, xDir[1]);
m->SetElement(2, 0, xDir[2]);
m->SetElement(0, 1, yDir[0]);
m->SetElement(1, 1, yDir[1]);
m->SetElement(2, 1, yDir[2]);
m->SetElement(0, 2, zDir[0]);
m->SetElement(1, 2, zDir[1]);
m->SetElement(2, 2, zDir[2]);
vtkSmartPointer<vtkMatrix4x4> m2 = vtkSmartPointer<vtkMatrix4x4>::New();
m2->DeepCopy(m);
m2->Invert();
// 2.将点和待求交的模型变换到局部坐标系内
transformationCurveData(m2, inputPoints);
transformationObjectData(m2, tPolyData);
// 3.计算模型包围盒,计算窗口尺寸
tPolyData->ComputeBounds();
double* bound = tPolyData->GetBounds();
double modelW = bound[1] - bound[0];
double modelH = bound[3] - bound[2];
int xres = 0;
int yres = 0;
double pix = 2000;
if(modelW > modelH)
{
xres = pix;
yres = pix*(modelH/modelW);
} else {
xres = pix*(modelW/modelH);
yres = pix;
}
// 4.创建相机视口
vtkSmartPointer<vtkRenderWindow> render_win = vtkSmartPointer<vtkRenderWindow>::New ();
vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New ();
render_win->SetPosition(100000, 100000); // 起始位置移除电脑屏幕外，防止界面看到
render_win->AddRenderer (renderer);
render_win->SetSize (xres, yres);
renderer->SetBackground (1.0, 0, 0);
//render view
vtkSmartPointer<vtkPolyDataMapper> mapper = vtkSmartPointer<vtkPolyDataMapper>::New ();
mapper->SetInputData(tPolyData);
mapper->Update ();
vtkSmartPointer<vtkActor> actor_view = vtkSmartPointer<vtkActor>::New ();
actor_view->SetMapper (mapper);
//
renderer->SetActiveCamera (cam);
renderer->AddActor (actor_view);
renderer->ResetCamera();
renderer->ResetCameraClippingRange();
renderer->GetActiveCamera()->ParallelProjectionOn();
render_win->Render ();
// 5.创建基于深度拾取器
vtkSmartPointer<vtkWorldPointPicker> worldPicker = vtkSmartPointer<vtkWorldPointPicker>::New ();
std::vector<bool> mark(inputPoints.size(), false);
for(int i = 0; i < inputPoints.size(); i++)
{
FreeVector p = inputPoints[i];
double worldCoord[3] = {p[0], p[1], p[2]};
vtkSmartPointer<vtkCoordinate> pCoorPress = vtkSmartPointer<vtkCoordinate>::New();
pCoorPress->SetCoordinateSystemToWorld();
pCoorPress->SetValue(worldCoord);
int *dispCoord = pCoorPress->GetComputedDisplayValue(renderer);
int x = dispCoord[0];
int y = dispCoord[1];
float value = render_win->GetZbufferDataAtPoint(x, y);
if(value == 1)
{
outPoints.push_back(FreeVector(0, 0, 0));
mark[i] = true;
continue;
}
worldPicker->Pick (x, y, value, renderer); // 计算交点
double coords[3] = {0};
worldPicker->GetPickPosition (coords);
p[2] = coords[2];
outPoints.push_back(p);
}
// 开始将点变换回去
// 旋转回去
transformationCurveData(m, outPoints);
for(int i = 0; i < outPoints.size(); i++) // 没投上的点置零
{
FreeVector& p = outPoints[i];
if(mark[i])
{
p = FreeVector(0, 0, 0);
}
}
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
/**
* @brief pointsDepthMapping 点集与模型求交
* @param inputPoints 输入点集
* @param mapNormal 投影方向
* @param polyData 模型
* @param outPoints 输出点集，个数与输入对应，没投上是0，0，0
*/
void pointsDepthMapping(std::vector<FreeVector> &inputPoints, FreeVector mapNormal, vtkPolyData *polyData, std::vector<FreeVector> &outPoints)
{
vtkSmartPointer<vtkPolyData> tPolyData = vtkSmartPointer<vtkPolyData>::New();
tPolyData->DeepCopy(polyData);
// 1.构建局部坐标系：将投影方向当做Z轴构建坐标系
FreeVector zDir, xDir, yDir;
zDir = -mapNormal;
zDir.Normalize();
constructCoordinateByZAxle(zDir, xDir, yDir); // 利用Z轴构建一个坐标系
vtkSmartPointer<vtkMatrix4x4> m = vtkSmartPointer<vtkMatrix4x4>::New();
m->SetElement(0, 0, xDir[0]);
m->SetElement(1, 0, xDir[1]);
m->SetElement(2, 0, xDir[2]);
m->SetElement(0, 1, yDir[0]);
m->SetElement(1, 1, yDir[1]);
m->SetElement(2, 1, yDir[2]);
m->SetElement(0, 2, zDir[0]);
m->SetElement(1, 2, zDir[1]);
m->SetElement(2, 2, zDir[2]);
vtkSmartPointer<vtkMatrix4x4> m2 = vtkSmartPointer<vtkMatrix4x4>::New();
m2->DeepCopy(m);
m2->Invert();
// 2.将点和待求交的模型变换到局部坐标系内
transformationCurveData(m2, inputPoints);
transformationObjectData(m2, tPolyData);
// 3.计算模型包围盒,计算窗口尺寸
tPolyData->ComputeBounds();
double* bound = tPolyData->GetBounds();
double modelW = bound[1] - bound[0];
double modelH = bound[3] - bound[2];
int xres = 0;
int yres = 0;
double pix = 2000;
if(modelW > modelH)
{
xres = pix;
yres = pix*(modelH/modelW);
} else {
xres = pix*(modelW/modelH);
yres = pix;
}
// 4.创建相机视口
vtkSmartPointer<vtkRenderWindow> render_win = vtkSmartPointer<vtkRenderWindow>::New ();
vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New ();
render_win->SetPosition(100000, 100000); // 起始位置移除电脑屏幕外，防止界面看到
render_win->AddRenderer (renderer);
render_win->SetSize (xres, yres);
renderer->SetBackground (1.0, 0, 0);
//render view
vtkSmartPointer<vtkPolyDataMapper> mapper = vtkSmartPointer<vtkPolyDataMapper>::New ();
mapper->SetInputData(tPolyData);
mapper->Update ();
vtkSmartPointer<vtkActor> actor_view = vtkSmartPointer<vtkActor>::New ();
actor_view->SetMapper (mapper);
//
renderer->SetActiveCamera (cam);
renderer->AddActor (actor_view);
renderer->ResetCamera();
renderer->ResetCameraClippingRange();
renderer->GetActiveCamera()->ParallelProjectionOn();
render_win->Render ();
// 5.创建基于深度拾取器
vtkSmartPointer<vtkWorldPointPicker> worldPicker = vtkSmartPointer<vtkWorldPointPicker>::New ();
std::vector<bool> mark(inputPoints.size(), false);
for(int i = 0; i < inputPoints.size(); i++)
{
FreeVector p = inputPoints[i];
double worldCoord[3] = {p[0], p[1], p[2]};
vtkSmartPointer<vtkCoordinate> pCoorPress = vtkSmartPointer<vtkCoordinate>::New();
pCoorPress->SetCoordinateSystemToWorld();
pCoorPress->SetValue(worldCoord);
int *dispCoord = pCoorPress->GetComputedDisplayValue(renderer);
int x = dispCoord[0];
int y = dispCoord[1];
float value = render_win->GetZbufferDataAtPoint(x, y);
if(value == 1)
{
outPoints.push_back(FreeVector(0, 0, 0));
mark[i] = true;
continue;
}
worldPicker->Pick (x, y, value, renderer); // 计算交点
double coords[3] = {0};
worldPicker->GetPickPosition (coords);
p[2] = coords[2];
outPoints.push_back(p);
}
// 开始将点变换回去
// 旋转回去
transformationCurveData(m, outPoints);
for(int i = 0; i < outPoints.size(); i++) // 没投上的点置零
{
FreeVector& p = outPoints[i];
if(mark[i])
{
p = FreeVector(0, 0, 0);
}
}
}