如何从VTK文件中提取六面体网格的外表面？

问题描述 #

我有一个由六面体组成的立方体模型，希望获取所有六面体的外表面（而非整体大立方体的表面），但未能实现。我认为可行方法是遍历所有面，筛选出无相邻面的面，但不确定是否正确。我的代码基于Python，同时欢迎C/C++实现方案。

初始尝试代码 #

# Read the source file.
reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName(file_name)
reader.Update()

output = reader.GetOutput()
points = output.GetPoints()

# get faces
faces = []
for i in range(output.GetNumberOfCells()):
    cell = output.GetCell(i)
    #print(type(cell)) hexaedron
    cell_points_ids = []
    for face_index in range(cell.GetNumberOfFaces()):
        a=vtk.reference([0,0,0,0])
        cell.GetFaceToAdjacentFaces(face_index,a)

更新尝试1 #

gf = vtk.vtkGeometryFilter()
gf.SetInputConnection(reader.GetOutputPort())
gf.Update()
polydata_output = gf.GetOutput()

print(polydata_output.GetNumberOfCells()) # 291024
print(output.GetNumberOfCells()) #48504 

更新尝试2 #

sc = vtk.vtkStaticCleanUnstructuredGrid()
sc.Update()
sc.SetInputData(output)

sc.Update()
gf = vtk.vtkGeometryFilter()
sc_output = sc.GetOutput()
gf.SetInputData(sc_output)
#gf.SetInputConnection(reader.GetOutputPort())
gf.Update()
polydata_output = gf.GetOutput()

print(polydata_output.GetNumberOfCells()) # 19120
print(output.GetNumberOfCells()) #48504 

解决方案 #

你的思路是正确的：遍历每个六面体的所有面，筛选出没有相邻六面体的面，这些就是单个六面体的外表面。vtkGeometryFilter默认提取的是整个网格的外表面（仅暴露在整个模型外部的面），无法满足需求，所以需要手动判断每个面的相邻情况。

Python 实现

import vtk

def extract_hexahedron_surfaces(file_name):
    # 读取网格数据
    reader = vtk.vtkXMLUnstructuredGridReader()
    reader.SetFileName(file_name)
    reader.Update()
    grid = reader.GetOutput()

    # 创建存储结果的PolyData
    result_polydata = vtk.vtkPolyData()
    points = vtk.vtkPoints()
    polygons = vtk.vtkCellArray()

    # 遍历每个六面体单元
    for cell_idx in range(grid.GetNumberOfCells()):
        cell = grid.GetCell(cell_idx)
        if cell.GetCellType() != vtk.VTK_HEXAHEDRON:
            continue  # 只处理六面体

        # 遍历六面体的每个面（共6个）
        for face_idx in range(cell.GetNumberOfFaces()):
            # 获取当前面的相邻单元信息
            adjacent_cells = vtk.reference([-1, -1, -1, -1])
            cell.GetFaceToAdjacentFaces(face_idx, adjacent_cells)
            
            # 判断是否有相邻单元：如果相邻单元ID都是-1，说明该面没有共享
            has_adjacent = False
            for adj_cell in adjacent_cells:
                if adj_cell != -1:
                    has_adjacent = True
                    break
            
            if not has_adjacent:
                # 获取当前面的点ID
                face = cell.GetFace(face_idx)
                point_ids = []
                for pt_idx in range(face.GetNumberOfPoints()):
                    global_pt_id = cell.GetPointId(face.GetPointId(pt_idx))
                    point_ids.append(global_pt_id)
                    # 添加点到结果
                    pt = grid.GetPoint(global_pt_id)
                    points.InsertNextPoint(pt)
                
                # 添加面到结果
                polygons.InsertNextCell(len(point_ids), point_ids)

    # 设置PolyData的点和面
    result_polydata.SetPoints(points)
    result_polydata.SetPolys(polygons)

    # 清理重复的点和面
    cleaner = vtk.vtkCleanPolyData()
    cleaner.SetInputData(result_polydata)
    cleaner.Update()
    cleaned_polydata = cleaner.GetOutput()

    # 保存结果
    writer = vtk.vtkXMLPolyDataWriter()
    writer.SetFileName("hexahedron_surfaces.vtp")
    writer.SetInputData(cleaned_polydata)
    writer.Write()

    return cleaned_polydata

# 调用示例
extract_hexahedron_surfaces("your_input_file.vtu")

C++ 实现

#include <vtkXMLUnstructuredGridReader.h>
#include <vtkPolyData.h>
#include <vtkPoints.h>
#include <vtkCellArray.h>
#include <vtkHexahedron.h>
#include <vtkCleanPolyData.h>
#include <vtkXMLPolyDataWriter.h>
#include <vtkSmartPointer.h>

vtkSmartPointer<vtkPolyData> ExtractHexahedronSurfaces(const std::string& fileName)
{
    // 读取网格数据
    auto reader = vtkSmartPointer<vtkXMLUnstructuredGridReader>::New();
    reader->SetFileName(fileName.c_str());
    reader->Update();
    auto grid = reader->GetOutput();

    // 创建结果PolyData
    auto resultPolydata = vtkSmartPointer<vtkPolyData>::New();
    auto points = vtkSmartPointer<vtkPoints>::New();
    auto polygons = vtkSmartPointer<vtkCellArray>::New();

    // 遍历每个六面体单元
    for (vtkIdType cellIdx = 0; cellIdx < grid->GetNumberOfCells(); ++cellIdx)
    {
        auto cell = vtkHexahedron::SafeDownCast(grid->GetCell(cellIdx));
        if (!cell)
            continue; // 跳过非六面体单元

        // 遍历六面体的每个面
        for (int faceIdx = 0; faceIdx < cell->GetNumberOfFaces(); ++faceIdx)
        {
            // 获取相邻单元ID
            vtkIdType adjacentCells[4];
            cell->GetFaceToAdjacentFaces(faceIdx, adjacentCells);

            // 判断是否有相邻单元
            bool hasAdjacent = false;
            for (int i = 0; i < 4; ++i)
            {
                if (adjacentCells[i] != -1)
                {
                    hasAdjacent = true;
                    break;
                }
            }

            if (!hasAdjacent)
            {
                // 获取当前面的点
                auto face = cell->GetFace(faceIdx);
                vtkIdType pointIds[4];
                face->GetPointIds()->GetIds(pointIds);

                // 转换为全局点ID并添加到点集合
                vtkIdType globalPointIds[4];
                for (int i = 0; i < 4; ++i)
                {
                    globalPointIds[i] = cell->GetPointId(pointIds[i]);
                    double pt[3];
                    grid->GetPoint(globalPointIds[i], pt);
                    points->InsertNextPoint(pt);
                }

                // 添加面到单元集合
                polygons->InsertNextCell(4, globalPointIds);
            }
        }
    }

    // 设置PolyData的数据
    resultPolydata->SetPoints(points);
    resultPolydata->SetPolys(polygons);

    // 清理重复点和面
    auto cleaner = vtkSmartPointer<vtkCleanPolyData>::New();
    cleaner->SetInputData(resultPolydata);
    cleaner->Update();

    // 保存结果
    auto writer = vtkSmartPointer<vtkXMLPolyDataWriter>::New();
    writer->SetFileName("hexahedron_surfaces.vtp");
    writer->SetInputData(cleaner->GetOutput());
    writer->Write();

    return cleaner->GetOutput();
}

// 调用示例
int main()
{
    ExtractHexahedronSurfaces("your_input_file.vtu");
    return 0;
}

关键说明 #

1. 相邻面判断：通过GetFaceToAdjacentFaces获取当前面的相邻单元ID，若所有ID都是-1，说明该面未与其他六面体共享，属于单个六面体的外表面。
2. 数据去重：使用vtkCleanPolyData清理重复的点和面，避免结果中出现冗余数据。
3. 与vtkGeometryFilter的区别：vtkGeometryFilter仅提取整个网格的外表面，而上述代码提取的是每个六面体自身未被共享的面，包括模型内部六面体的暴露面。

内容的提问来源于stack exchange，提问作者manuersuper