当前位置：首页 > news >正文

【实战演练】基于VTK的散点凹包计算实战：从代码逻辑到实现思路

news 2025/10/15 6:30:37

实战问题背景：

在地理信息系统、计算机视觉、机器人路径规划等领域，凹包（Concave Hull）是描述点集边界的重要形态，相比凸包能更精确地贴合非凸分布的点集。本文将基于VTK（Visualization Toolkit）库，解析一段散点凹包计算代码的设计思路与实现逻辑，帮助开发者理解如何利用VTK的几何处理能力快速构建凹包算法。

凹包计算的核心需求与VTK优势

凹包的核心需求是：对于给定的散点集，生成一个能包裹所有点且尽可能贴合点集分布的非凸多边形边界。与凸包（Convex Hull）不同，凹包允许边界向内凹陷，更适合描述“有凹痕”的点集形态（如湖泊轮廓、复杂区域边界等）。

VTK作为开源可视化库，提供了丰富的几何处理工具（如三角化、边界提取、拓扑操作等），无需从零实现基础几何算法，可大幅简化凹包计算的开发流程。本文解析的代码正是基于VTK的这些特性，实现了从3D散点到凹包边界的自动化计算。

代码整体设计思路

代码核心函数为CalcConvexhull（注：函数名可能存在笔误，实际功能为凹包计算），输入为3D点集（vtkPoints*）、输出结果（vtkPolyData*）和容差参数（tol），整体流程可分为5个关键步骤：

3D点集降维到2D平面：将输入的3D点投影到指定平面（此处为Z=0平面），因为凹包计算通常在2D空间中进行；
Delaunay三角化构建初始网格：通过三角化将散点连接为三角形网格，为后续边界提取提供基础；
边界边筛选与标记：从三角化网格中提取边界边，并基于容差参数筛选出“有效边界边”；
迭代剔除内部三角形：通过删除与有效边界边关联的内部三角形，逐步暴露凹包边界；
边界整合与输出：清理冗余数据，连接连续边界线段，生成最终的凹包多边形。

代码逻辑分步解析

步骤1：3D点集投影到2D平面

double normz[] = { 0,0,1 }; // Z轴方向法向量（对应Z=0平面）
vtkNew<vtkTransform> trans;
RotateToPlane(pts, trans, normz); // 将点集旋转到Z=0平面

VTK中凹包计算依赖2D几何处理，因此需先将3D点集投影到平面。这里通过RotateToPlane函数（未展示实现，推测为坐标变换）将点集旋转至Z=0平面，并用vtkTransform记录变换矩阵，后续处理将基于该平面进行。

步骤2：Delaunay三角化生成初始网格

vtkNew<vtkPolyData> ptsData;
ptsData->SetPoints(pts); // 包装点集为PolyData
vtkNew<vtkDelaunay2D> del;
del->SetTransform(trans); // 应用平面变换
del->SetInputData(ptsData);
del->Update(); // 执行Delaunay三角化

Delaunay三角化是将散点连接为无重叠、无交叉三角形网格的经典算法，其生成的网格边界是点集的凸包候选。这里通过vtkDelaunay2D对投影后的2D点集进行三角化，得到包含所有点的三角形网格，作为后续处理的“原始画布”。

步骤3：提取边界边并筛选有效边

// 为点添加原始ID标记，便于后续追踪
vtkNew<vtkGenerateIds> idsFilter;
idsFilter->SetInputData(del->GetOutput());
idsFilter->SetPointIdsArrayName("origin");
idsFilter->Update();
auto mesh = vtkPolyData::SafeDownCast(idsFilter->GetOutput());// 提取网格的边界边（非共享边）
vtkNew<vtkFeatureEdges> edgeFilter;
edgeFilter->SetInputData(mesh);
edgeFilter->SetBoundaryEdges(1); // 只提取边界边
edgeFilter->SetFeatureEdges(0);
edgeFilter->SetManifoldEdges(0);
edgeFilter->SetNonManifoldEdges(0);
edgeFilter->Update();
auto edges = edgeFilter->GetOutput();

三角化网格的“边界边”是指仅属于一个三角形的边（凸包的初始边界）。通过vtkFeatureEdges可快速提取这类边。同时，vtkGenerateIds为点添加原始ID，避免后续处理中坐标变换导致的ID混乱。

// 基于容差筛选有效边界边（长度≥tol的边）
vtkNew<vtkCellArray> InvalidEdgeLines;
auto idArr = edges->GetPointData()->GetArray("origin");
for (int i = 0; i < nbound; i++) {vtkNew<vtkIdList> eIds;edges->GetCellPoints(i, eIds); // 获取边的两个端点ID// 计算两点距离double pt0[3], pt1[3];edges->GetPoint(pid0, pt0);edges->GetPoint(pid1, pt1);if (sqrt(vtkMath::Distance2BetweenPoints(pt0, pt1)) >= tol) {// 转换为原始ID并保存for (int e = 0; e < eIds->GetNumberOfIds(); e++) {eIds->SetId(e, idArr->GetVariantValue(eIds->GetId(e)).ToInt());}InvalidEdgeLines->InsertNextCell(eIds);}
}

容差参数tol用于过滤短边界边（可能是噪声或局部凸起导致），仅保留长度足够的边作为“有效边界边”，存入InvalidEdgeLines（命名可能为“待处理边界边”）。

步骤4：迭代剔除内部三角形，暴露凹包边界

// 循环处理有效边界边，删除关联的内部三角形
while (InvalidEdgeLines->GetNumberOfCells() > 0) {vtkNew<vtkCellArray> tempEdges;InvalidEdgeLines->InitTraversal();vtkNew<vtkIdList> eIds;while (InvalidEdgeLines->GetNextCell(eIds)) {// 找到边的两个端点共同所属的三角形（唯一）int pid0 = eIds->GetId(0);int pid1 = eIds->GetId(1);vtkNew<vtkIdList> triangleIds0, triangleIds1;mesh->GetPointCells(pid0, triangleIds0); // 包含pid0的三角形mesh->GetPointCells(pid1, triangleIds1); // 包含pid1的三角形triangleIds0->IntersectWith(triangleIds1); // 求交集（共同三角形）int delCellId = triangleIds0->GetId(0); // 待删除的内部三角形// 提取该三角形的其他边，若长度≥tol则加入下一轮处理vtkCell* delCell = mesh->GetCell(delCellId);for (int e = 0; e < delCell->GetNumberOfEdges(); e++) {auto edCell = delCell->GetEdge(e);vtkIdList* edgeIds = edCell->GetPointIds();if (edgeIds->IsId(pid0) != -1 && edgeIds->IsId(pid1) != -1) {continue; // 跳过当前边}// 计算边长度，符合条件则加入tempEdgesdouble p0[3], p1[3];mesh->GetPoint(edgeIds->GetId(0), p0);mesh->GetPoint(edgeIds->GetId(1), p1);if (sqrt(vtkMath::Distance2BetweenPoints(p0, p1)) >= tol) {tempEdges->InsertNextCell(edgeIds);}}mesh->DeleteCell(delCellId); // 删除内部三角形}InvalidEdgeLines->ShallowCopy(tempEdges); // 更新待处理边mesh->RemoveDeletedCells(); // 清理已删除的三角形
}

这是凹包计算的核心步骤：通过迭代删除与有效边界边关联的三角形，逐步“剥离”凸包内部的三角形，暴露凹进去的边界。每次删除一个三角形后，其未被处理的边（长度≥tol）将作为新的待处理边，直到没有更多边需要处理，此时剩余的网格边界即为凹包。

步骤5：边界整合与输出

// 清理冗余数据（重复点、无效单元）
vtkNew<vtkCleanPolyData> clean;
clean->SetInputData(mesh);
clean->Update();// 提取最终边界并连接为连续线段
vtkNew<vtkFeatureEdges> edgeFilter;
edgeFilter->SetInputData(clean->GetOutput());
edgeFilter->Update();
vtkNew<vtkStripper> strip;
strip->SetJoinContiguousSegments(true); // 连接连续线段
strip->SetInputData(edgeFilter->GetOutput());
strip->Update();// 筛选最长连续线段作为凹包边界
auto temp = strip->GetOutput();
int nmaxCell = temp->GetMaxCellSize();
for (int pc = 0; pc < temp->GetNumberOfCells(); pc++) {if (temp->GetCellSize(pc) != nmaxCell)temp->DeleteCell(pc);
}
temp->RemoveDeletedCells();// 输出结果
clean->SetInputData(temp);
clean->Update();
out->ShallowCopy(clean->GetOutput());

最后通过vtkCleanPolyData清理冗余点和单元，vtkStripper将分散的边界线段连接为连续多边形，再筛选出最长的连续线段（即凹包的完整边界），最终输出到out。

完整代码

void CalcConvexhull(vtkPoints* pts, vtkPolyData* out, double tol)
{double normz[] = { 0,0,1 };vtkNew<vtkTransform> trans;RotateToPlane(pts, trans, normz);//计算变化vtkNew<vtkPolyData> ptsData;ptsData->SetPoints(pts);vtkNew<vtkDelaunay2D> del;del->SetTransform(trans);del->SetInputData(ptsData);del->Update();//out->ShallowCopy(del->GetOutput());//return;vtkNew<vtkGenerateIds> idsFilter;idsFilter->SetInputData(del->GetOutput());idsFilter->SetCellIds(0);idsFilter->SetFieldData(0);idsFilter->SetPointIdsArrayName("origin");idsFilter->Update();auto mesh = vtkPolyData::SafeDownCast(idsFilter->GetOutput());vtkNew<vtkFeatureEdges> edgeFilter;edgeFilter->SetInputData(mesh);edgeFilter->SetBoundaryEdges(1);edgeFilter->SetFeatureEdges(0);edgeFilter->SetManifoldEdges(0);edgeFilter->SetNonManifoldEdges(0);edgeFilter->Update();auto edges = edgeFilter->GetOutput();int nbound = edges->GetNumberOfLines();vtkNew<vtkCellArray> InvalidEdgeLines;auto idArr = edges->GetPointData()->GetArray("origin");//第一次计算所有无效边界边for (int i = 0; i < nbound; i++) {vtkNew<vtkIdList> eIds;edges->GetCellPoints(i, eIds);assert(eIds->GetNumberOfIds() == 2);int pid0 = eIds->GetId(0);int pid1 = eIds->GetId(1);double pt0[3], pt1[3];edges->GetPoint(pid0, pt0);edges->GetPoint(pid1, pt1);//如果边界边长度超过阈值，加入无效集合if (sqrt(vtkMath::Distance2BetweenPoints(pt0, pt1)) >= tol) {//转换idfor (int e = 0; e < eIds->GetNumberOfIds(); e++) {eIds->SetId(e, idArr->GetVariantValue(eIds->GetId(e)).ToInt());}InvalidEdgeLines->InsertNextCell(eIds);}}//遍历所有无效边边界while (InvalidEdgeLines->GetNumberOfCells() > 0) {vtkNew<vtkCellArray> tempEdges;InvalidEdgeLines->InitTraversal();vtkNew<vtkIdList> eIds;while (InvalidEdgeLines->GetNextCell(eIds)) {//查找edge关联的单元delCellIdassert(eIds->GetNumberOfIds() == 2);int pid0 = eIds->GetId(0);int pid1 = eIds->GetId(1);vtkNew<vtkIdList> triangleIds0, triangleIds1;mesh->GetPointCells(pid0, triangleIds0);mesh->GetPointCells(pid1, triangleIds1);triangleIds0->IntersectWith(triangleIds1);int nlinkCell = triangleIds0->GetNumberOfIds();if (nlinkCell == 0)continue;assert(nlinkCell == 1);int delCellId = triangleIds0->GetId(0);vtkCell* delCell = mesh->GetCell(delCellId);int nedge = delCell->GetNumberOfEdges();for (int e = 0; e < nedge; e++) {auto edCell = delCell->GetEdge(e);vtkIdList* edgeIds = edCell->GetPointIds();if (edgeIds->IsId(pid0) != -1 && edgeIds->IsId(pid1) != -1) {continue;}double p0[3], p1[3];mesh->GetPoint(edgeIds->GetId(0), p0);mesh->GetPoint(edgeIds->GetId(1), p1);//如果边界边长度不超过阈值，继续下一条if (sqrt(vtkMath::Distance2BetweenPoints(p0, p1)) < tol) {continue;}//将新的无效边加入tempEdges->InsertNextCell(edgeIds);}mesh->DeleteCell(delCellId);}InvalidEdgeLines->ShallowCopy(tempEdges);//更新mesh的拓扑mesh->RemoveDeletedCells();}vtkNew<vtkCleanPolyData> clean;clean->SetInputData(mesh);clean->Update();//重新提取边界edgeFilter->SetInputData(clean->GetOutput());edgeFilter->Update();vtkNew<vtkStripper> strip;strip->SetJoinContiguousSegments(true);strip->SetInputData(edgeFilter->GetOutput());strip->Update();auto temp = (strip->GetOutput());//抽取最大Cellint nmaxCell = temp->GetMaxCellSize();for (int pc = 0; pc < temp->GetNumberOfCells(); pc++) {if (temp->GetCellSize(pc) != nmaxCell)temp->DeleteCell(pc);}temp->RemoveDeletedCells();//清理clean->SetInputData(temp);clean->Update();out->ShallowCopy(clean->GetOutput());vtkNew<vtkIdList> lids;out->GetLines()->GetCell(0, lids);int npt = out->GetNumberOfPoints();//保证点不闭合if (npt < lids->GetNumberOfIds())lids->Resize(npt);out->GetLines()->Initialize();out->GetLines()->InsertNextCell(lids);out->Modified();return;
}

关键参数与注意事项

容差tol的设置：tol决定了边界边的最小长度，过小会保留过多噪声导致边界杂乱，过大会丢失细节，需根据点集密度调整（建议设为点平均间距的1-2倍）。
3D到2D的投影平面：代码中固定投影到Z=0平面，实际应用中可根据点集分布动态选择投影平面（如通过主成分分析确定最佳投影方向）。
效率考量：Delaunay三角化和迭代删除步骤的时间复杂度与点数量正相关，对于大规模点集（10万+），建议先降采样再计算。