台南市网站建设_网站建设公司_博客网站_seo优化

一：主要的知识点

1、说明

本文只是教程内容的一小段，因博客字数限制，故进行拆分。主教程链接：vtk教程——逐行解析官网所有Python示例-CSDN博客

2、知识点纪要

本段代码主要涉及的有①如何根据顶点索引找到邻接顶点的索引，②特征边缘，③曲率的计算

二：代码及注释

import vtkmodules.vtkRenderingOpenGL2 from vtkmodules.vtkIOXML import vtkXMLPolyDataReader from vtkmodules.vtkFiltersGeneral import vtkCurvatures from vtkmodules.vtkCommonDataModel import vtkPolyData from vtkmodules.numpy_interface import dataset_adapter as dsa from vtkmodules.vtkFiltersCore import vtkIdFilter, vtkFeatureEdges from vtkmodules.vtkCommonCore import vtkIdList, VTK_DOUBLE import numpy as np from vtkmodules.vtkRenderingAnnotation import vtkScalarBarActor from vtkmodules.util import numpy_support from vtkmodules.vtkCommonColor import vtkColorSeries, vtkNamedColors from vtkmodules.vtkRenderingCore import ( vtkActor, vtkColorTransferFunction, vtkPolyDataMapper, vtkRenderWindow, vtkRenderWindowInteractor, vtkRenderer ) from vtkmodules.vtkInteractionWidgets import vtkCameraOrientationWidget def point_neighbourhood(source: vtkPolyData, pt_id: int): """ 找到对应索引顶点的拓扑邻接顶点索引 """ cell_ids = vtkIdList() """ GetPointCells 拿到模型中与点p_id相连的所有面的ID """ source.GetPointCells(pt_id, cell_ids) neighbour = set() for cell_idx in range(0, cell_ids.GetNumberOfIds()): cell_id = cell_ids.GetId(cell_idx) cell_point_ids = vtkIdList() """ GetCellPoints 拿到指定单元(cell)的所有顶点的ID """ source.GetCellPoints(cell_id, cell_point_ids) for cell_pt_idx in range(0, cell_point_ids.GetNumberOfIds()): neighbour.add(cell_point_ids.GetId(cell_pt_idx)) return neighbour def compute_distance(source: vtkPolyData, pt_id_a: int, pt_id_b: int) -> float: pt_a = np.array(source.GetPoint(pt_id_a)) pt_b = np.array(source.GetPoint(pt_id_b)) return np.linalg.norm(pt_a - pt_b) def adjust_edge_curvatures(source: vtkPolyData, curvature_name: str, epsilon=1.0e-08) -> None: """ source.GetPointData().SetActiveScalars(cruvature_name) 这段代码的含义是 将模型的指定曲率数组设置为活动标量（active scalars） 一个 vtkPolyData 对象可以包含多个不同的数据数组，比如法线，TextureCoordinates 纹理坐标，Gauss_Curvature 等 source: 指的是你的 vtkPolyData 对象，也就是牛头模型的数据。 GetPointData(): 获取与模型每个点相关联的所有数据数组的集合（vtkPointData）。这包括坐标、法线、纹理坐标以及像曲率这样的自定义数据。 SetActiveScalars(curvature_name): 这是最关键的方法。它告诉 VTK 的管线，在接下来的操作中（例如映射和渲染），应该使用哪个数据数组作为主要的标量数据。 在这里，curvature_name 可以是 'Gauss_Curvature' 或 'Mean_Curvature'。通过这行代码，你明确指定了： “请把这个模型上的 Gauss_Curvature（或 Mean_Curvature）数组作为主要颜色值来源，用于渲染。 """ source.GetPointData().SetActiveScalars(curvature_name) """ dataset_adapter 了把 VTK 的数据对象（vtkDataSet, vtkArray, …）和 NumPy 数组 之间建立一个桥梁 """ np_source = dsa.WrapDataObject(source) """ WrapDataObject 把一个 VTK 数据对象包一层壳（adapter），让它可以直接用 NumPy 的方式访问 np_source.Points → 得到点坐标的 numpy.ndarray np_source.PointData -> 得到点数据（标量场/向量场）的数组字典 np_source.CellData -> 得到单元数据 """ curvatures = np_source.PointData[curvature_name] # 得到了各个点的高斯曲率大小，长度与点数保持一致 # 下面的一段代码是给数据source的点生产成一个ID数组，并存储在指定名字的数组里 array_name = 'ids' """ vtkIdFilter 把隐含的索引（点/单元 ID）转换成显式的数组，在数据中生成点 ID 和/或单元 ID，并把它们存储为数组，方便后续访问、映射或可视化 在 VTK 里，几乎所有的数据集都有 点 (Points) 和 单元 (Cells)： 点 ID：每个点的索引（从 0 开始编号） 单元 ID：每个单元的索引（也是从 0 开始编号） 平时这些 ID 是隐含存在的（你用 GetPoint(i)、GetCell(i) 可以拿到），但并不是数据数组的一部分。 如果你想把这些 ID 当作数组存到 PointData 或 CellData 中，就可以用 vtkIdFilter """ id_filter = vtkIdFilter() id_filter.SetInputData(source) id_filter.SetPointIds(True) # 是否生成点ID id_filter.SetCellIds(False) # 是否生成单元的ID id_filter.SetPointIdsArrayName(array_name) # 指定生成的ID数组的名字 # id_filter.SetCellIdsArrayName(array_name) 没有生成单元的ID，没必要 id_filter.Update() """ vtkFeatureEdges 从一个曲面（PolyData）中提取“特征边缘”，并把这些边缘输出为一条条 polyline（折线） 什么是“特征边缘”？ 在 VTK 里，三角形网格的边大体可以分成几类： Boundary Edges（边界边） 只属于一个多边形的边（即网格开口处的边）。 打开：BoundaryEdgesOn() Feature Edges（特征边） 被两个多边形共享，但这两个面的夹角大于某个阈值（默认 30°）。 打开：FeatureEdgesOn() 用 SetFeatureAngle(angle) 设置角度阈值。 Manifold Edges（流形边） 被恰好两个多边形共享，且夹角小于阈值（即正常相邻面）。 打开：ManifoldEdgesOn() Non-Manifold Edges（非流形边） 被两个以上的多边形共享。 打开：NonManifoldEdgesOn() """ edges = vtkFeatureEdges() edges.SetInputConnection(id_filter.GetOutputPort()) edges.BoundaryEdgesOn() # 检测边缘边 edges.ManifoldEdgesOff() # 不检测流行边 edges.NonManifoldEdgesOff() # 不检测非流行边 edges.FeatureEdgesOff() # 不检测特征边 edges.Update() edge_array = edges.GetOutput().GetPointData().GetArray(array_name) boundary_ids = [] # mesh边缘顶点的索引 for i in range(edges.GetOutput().GetNumberOfPoints()): boundary_ids.append(edge_array.GetValue(i)) # 删除可能存在的边缘索引 p_ids_set = set(boundary_ids) count_invalid = 0 for p_id in boundary_ids: # 获取到p_id这个索引定点的邻接点的索引 p_ids_neighbours = point_neighbourhood(source, p_id) # 这里减去了边缘索引，说明只想保留非边缘顶点的索引 p_ids_neighbours -= p_ids_set curvs = [curvatures[p_id_n] for p_id_n in p_ids_neighbours] # 计算两个顶点的欧式距离 dists = [compute_distance(source, p_id_n, p_id) for p_id_n in p_ids_neighbours] curvs = np.array(curvs) dists = np.array(dists) # 避免因计算精度原因导致的出现距离为0的现象，或者避免比较自身点 curvs = curvs[dists > 0] dists = dists[dists > 0] if len(curvs) > 0: # 距离的大小与权重成反比 weights = 1 / np.array(dists) weights /= weights.sum() new_curv = np.dot(curvs, weights) else: count_invalid += 1 new_curv = 0.0 curvatures[p_id] = new_curv if epsilon != 0: """ 下面这段代码的含义是指：如果某个曲率值的绝对值小于阈值 epsilon，就把它替换成 0，否则保持原值 """ curvatures = np.where(abs(curvatures) < epsilon, 0, curvatures) """ 把 NumPy 数组转换成 VTK 的 vtkDataArray deep=True 表示复制一份数据，而不是引用同一块内存 """ curv = numpy_support.numpy_to_vtk(num_array=curvatures.ravel(), deep=True, array_type=VTK_DOUBLE) curv.SetName(curvature_name) # 给这个vtk数组命名 source.GetPointData().RemoveArray(curvature_name) # 将之前同名的数组删掉，避免冲突 source.GetPointData().AddArray(curv) source.GetPointData().SetActiveScalars(curvature_name) def main(): file_name = "Data/cowHead.vtp" color_map_idx = 16 gaussian_curvature = True if gaussian_curvature: curvature = "Gauss_Curvature" # 高斯曲率 else: curvature = "Mean_Curvature" # 平均曲率 reader = vtkXMLPolyDataReader() reader.SetFileName(file_name) reader.Update() source = reader.GetOutput() """ vtkCurvatures 计算曲面曲率的一个滤波器 输入：vtkPolyData（通常是一个三角网格曲面）。 输出：在网格的 点数据（point data） 中生成一个标量数组，表示每个点的曲率值。 这个过滤器会将计算出的曲率值作为新的点数据数组添加到模型中 """ cc = vtkCurvatures() cc.SetInputConnection(reader.GetOutputPort()) if gaussian_curvature: cc.SetCurvatureTypeToGaussian() # 计算高斯曲率 cc.Update() else: cc.SetCurvatureTypeToMean() # 计算平均曲率 cc.Update() """ adjust_edge_curvatures 这个函数通过将曲面边缘的曲率值替换为其邻域内点的曲率平均值，来调整边缘的曲率 该函数首先通过 vtkFeatureEdges 过滤器识别出模型的边界点。 然后，它遍历每个边界点，并计算其邻近内部点的曲率值的加权平均值。 最后，它用这个加权平均值替换掉原始的边界曲率值，从而平滑了边界上的曲率，使得渲染结果更加合理 """ adjust_edge_curvatures(cc.GetOutput(), curvature) """ GetPointData().GetAbstractArray 和 GetPointData().GetArray的区别 GetArray(name_or_index) 返回类型：vtkDataArray（或其子类） vtkDataArray 是数值数组，专门存 数值型标量/向量/张量，比如 float、double、int 也就是说，GetArray() 只适用于 数值数组，对非数值的数组（例如字符串）会失败或者返回 None GetAbstractArray(name_or_index) 返回类型：vtkAbstractArray vtkAbstractArray 是 vtkDataArray 的基类，它更泛化，可以表示 任意类型的数组，不仅仅是数值型 数值数组：vtkDoubleArray、vtkFloatArray … 字符串数组：vtkStringArray 稀疏数组等其他类型 适合在你 不确定数组类型 时使用，比如不知道是数值还是字符串 """ source.GetPointData().AddArray(cc.GetOutput().GetPointData().GetAbstractArray(curvature)) """ GetScalars 返回 点数据 (PointData) 或单元数据 (CellData) 中的 “活动 (active)”标量数组 下面是一种类似语法糖的写法，当正确设置活动标量后 应该这样写 source.GetPointData().GetScalars().GetRange() GetRange() 获取这个标量值的最大和最小值 """ scalar_range = source.GetPointData().GetScalars(curvature).GetRange() # vtkColorSeries 快速生成一组颜色用于着色 color_series = vtkColorSeries() # SetColorScheme 选择某一种预定义的配色方案 color_series.SetColorScheme(color_map_idx) """ vtkColorTransferFunction 专门用来定义标量数值到颜色之间映射关系 的类 输入：一个 标量值（通常是数据里的 scalar，比如密度、强度、温度等）。 输出：一个 颜色（RGB 值）。 """ lut = vtkColorTransferFunction() lut.SetColorSpaceToHSV() for i in range(0, color_series.GetNumberOfColors()): color = color_series.GetColor(i) # rgb三色，范围0~255 double_color = list(map(lambda x: x / 255.0, color)) # 映射到0~1 """ 下面这段代码，是将标量范围线性映射到颜色索引 把整个标量范围 等分 为 N-1 段（N = 颜色数），第 i 个颜色点对应的标量位置就是 t """ t = scalar_range[0] + (scalar_range[1] - scalar_range[0]) / (color_series.GetNumberOfColors() - 1) * i """ AddRGBPoint 往lut添加一个颜色控制点 t： 这个控制点对应的标量值 double_color[0], double_color[1], double_color[2] → 这个标量值对应的 RGB 颜色 这样，渲染时数据里的数值如果等于 t，就会显示为这个 RGB； 如果数值在两个控制点之间，VTK 会自动插值 """ lut.AddRGBPoint(t, double_color[0], double_color[1], double_color[2]) colors = vtkNamedColors() mapper = vtkPolyDataMapper() mapper.SetInputData(source) mapper.SetScalarModeToUsePointFieldData() # 使用点数据上色 mapper.SelectColorArray(curvature) # 选择具体的标量数组进行上色 mapper.SetScalarRange(scalar_range) # 设置标量值映射到颜色查找表的范围，超过的会被裁剪到边界颜色 mapper.SetLookupTable(lut) actor = vtkActor() actor.SetMapper(mapper) window_width = 800 window_height = 800 """ vtkScalarBarActor 专门用于创建和显示一个颜色标尺（color bar 或 legend） 它的主要作用是为你的可视化结果添加一个解释性的图例。 当你使用颜色来映射数据（比如温度、压力或曲率）时， 这个颜色标尺会告诉你每种颜色代表的数值范围，从而让你的可视化结果更具可读性和科学性 """ scalar_bar = vtkScalarBarActor() scalar_bar.SetLookupTable(mapper.GetLookupTable()) scalar_bar.SetTitle(curvature.replace('_', '\n')) """ UnconstrainedFontSizeOn 让颜色标尺上的字体大小不受限制，从而能够根据标尺的大小自动调整 """ scalar_bar.UnconstrainedFontSizeOn() scalar_bar.SetNumberOfLabels(5) scalar_bar.SetMaximumWidthInPixels(window_width // 8) scalar_bar.SetMaximumHeightInPixels(window_height // 3) renderer = vtkRenderer() ren_win = vtkRenderWindow() ren_win.AddRenderer(renderer) ren_win.SetSize(window_width, window_height) ren_win.SetWindowName('Curvatures') iren = vtkRenderWindowInteractor() iren.SetRenderWindow(ren_win) # Important: The interactor must be set prior to enabling the widget. iren.SetRenderWindow(ren_win) """ vtkCameraOrientationWidget 控制相机视角的组件 """ cam_orient_manipulator = vtkCameraOrientationWidget() cam_orient_manipulator.SetParentRenderer(renderer) # Enable the widget. cam_orient_manipulator.On() # Add the actors to the scene renderer.AddActor(actor) renderer.AddViewProp(scalar_bar) renderer.SetBackground(colors.GetColor3d('DarkSlateGray')) # Render and interact ren_win.Render() iren.Start() print("") if __name__ == '__main__': main()