VTK tips

记录一些VTK filter使用以及基于VTK进行二次开发过程中的tricks以及容易混淆的地方和注意事项。

Clip filter 性能问题

在一次实验中发现vtk原本的clip filter速度特变慢，但是换成了vtkTableClipfilter之后速度就变得很快了，差距有10多倍的样子。这里记录一下它们的区别。

大致上它们的区别在与tableClip采用了空间换时间的思路，通过查表的方式提升了一些速度。本质上它们的时间复杂度都是O(N)的。原本的Clip filter中有一个点合并的操作，这个操作是相对比较费时间的。

clip通常的使用是把一个数据切开，然后看里面的结果，但最近发现，它在分割数据的时候还是有一些妙用，特别是对于基于集合坐标的划分上。比如一个压气机叶片模拟的结果，设计的时候动叶的区域与静叶的区域应该是不同的zone，旋转复制出来一周看效果的话，动叶可能是32份然后静叶可能是16份。但是这里的问题是模拟的数据可能把一个流到合在一起成一个zone计算了，这样就不好分成不同的区域进行旋转了。一个比较好的思路就是实用clip从这两个zone分界的地方把数据给切开，因为数据是比较规则的，通过几何定义能很轻松的找到分界线的区域，这样就把一个整体的数据分成了多个zone，然后可以单独对每个zone进行不同的操作。

Slice丢cell data 的问题

有一个项目中需要对存在cell value的值进行切片，比如有一个划分好的网格，每个网格上要有一个质量监测的参数，这个时候如果用slice的话，那新的cell和旧的cell其实就对应了，新的cell是新生成的，这个时候用vtk extracegeometry是比较好的方式

        plane = vtk.vtkPlane()
        plane.SetOrigin(origin)
        plane.SetNormal(normal)

        extract = vtk.vtkExtractGeometry()
        extract.SetImplicitFunction(plane)

        extract.ExtractBoundaryCellsOn()
        extract.ExtractInsideOn()

        extract.SetInputData(merged_block_ds)
        extract.Update()

这个时候会把穿过结合平面的所有cell都保存下来，结果了ExtrctInsideOn之后的效果有点类似于clip，只不过平面穿过的地方有点凹凸不平的感觉。

三角化的操作

直接说三角化的时候容易混淆triangulate以及d3两个操作的区别，vtkTriangulate 只处理曲面多边形，把它们拆成三角面；不填充体、不移动点；输出 PolyData。其主要目的是网格单元类型的转换。将非三角形单元（如四边形、多边形）转换为三角形单元，保持原有几何形状不变。

而Delaunay 2D、3D会对所有的点进行二维或者三维的三角剖分，会产生新的单元。整个网格数据的形状会发生改变。主要使用的场景是从点云生成网格，从一组无序的二维或投影到平面的三维点生成一个优化的三角网格。

Cell array的使用场景

好多时候都用的point array 但是 cell array 似乎不常用。一个使用场景是求叶轮机械的子屋面的之后，会沿着轴向对左右的值进行平均，这个时候会对新的除去了theta之后的平面重新划分网格，之后mapping的时候把所有的值映射到对应的grid中，之后再对grid中的值求平均，然后就是cell value，在绘制的时候再把cell value 映射到point value。

求网格划分的质量的时候也是一个使用cell array的场景，每个网格有一个独立的值。

structured gird, unstructured grid, image 区别

对于structured grid 理解的时候，这里所谓的结构化，是说每个cell可以按照 i j k 来进行索引，不存connectiveity data。具体cell的形状可以是正方形的（比如2d场景下），或者是平行四边形的，或者是曲线的（structured curvlinear grid）这里说的structured，也就是结构化或者是规则网格，主要是说的是否为规则的拓扑结构，而不是说的横平竖直的那种情况。

有一个误区就是，如果所有的cell都是一样的，比如全是quad或者是hexahedral （六面体）就能使用structured grid，实际并非如此，比如一个弯曲的空间形状，上面2个cell，下面一个cell，这样即使都是一样的cell，也需要把所有connectiveity保存下来，这就不是structured grid了，而是unstructured grid，因为这种结构没法用i j k 去对每个cell进行索引。

vtk 里面的image data 可以理解为横平竖直的 structured grid。
structured grid 与unstructured grid的本质区别是是否必须要存储connectivity，也就是节点之间的连接关系。如果直接可以通过存储point(i,j,k)就表示网格的形状了，那就是结构化网格，对应的assumpition就是每个cell都是按照一定的规则或者顺序进行连接的，反之如果必须得存储connectivity，这样的就是非结构的网格了。

structured gird的 i j k 索引

structure grid构建的时候就是按照 i j k 的顺序插入很多的点，然后网格会自动根据点的顺序生成。具体的cell可以是扭曲的那种长方形，这是一个相关的例子，比如在求子午面然后画网格的时候，只要计算出来每个关键节点的坐标，剩下的cell就能自动生成了。

        sg = vtk.vtkStructuredGrid()
        sg.SetDimensions(Ns,Nr,1)
        points = vtk.vtkPoints()
        # x 先变化 然后 y 再变化
        for j in range(Nr): 
            for i in range(Ns):      
                x = float(grid_points[i,j][0])   # X 用 z 坐标
                y = float(grid_points[i,j][1])   # Y 用 r 坐标
                z = 0.0                    # 平面网格，Z=0
                points.InsertNextPoint(x, y, z)
        sg.SetPoints(points)

这种网格便利的地方就是按照i j k 进行切片，比如对于tecplot来说，它有按照 i j k 进行切片的方式，如果压气机设计中的一个多级的流道，每个流道是按照 i j k 进行结构化网格索引的，后续的分析操作就比较方便了。这样的好处就是简化了几何形状，如果加上叶片的话，肯定只能是用无结构网格了，比如对于压气机流到模拟的话，这样的方式就不加上叶片了，之对流道进行划分，这样才能保证空间是封闭的，然后按照结构化网格的方式进行划分

几何形变

最近遇到一些关于几何形变的方法的使用，比如warpbyscalar，warpbyvector感觉还是挺奇妙的。

warpByScalar可以适用于二维或者三维数据，把几何形状按照scalar的值进行一个拉伸的操作，具体scale的方向可以手动进行指定，也可以自动计算每个point的normal的方向然后进行扩展。类似地还有warpByVector的操作，这要求数据里面包含有矢量，之后scale的时候需要按照这个vector的方向进行形变。

一个实际的应用是那种算应力形变的场景。比如一个叶片上的应力或者应变，在采用了warpByScalar之后就可以通过形变很明显地观察到哪里的受力大一些，很生动形象。或者是那种二维升三维的场景，比如海平面的一个标量值升到三维，升维之后观察起来会非常直观。

二次开发的思路

在做一些二次开发的时候如果把自己的工作定位成一个写vtk算法的层面的话就会很痛苦，相当于是重新造轮子，而且还可能造的还没有原来的好。这样既没有成就感也没有做真正有用的事情。以vtk或者paraview的这种成熟的程度，基础可视化算法相对于vtk就等价于数组运算对于numpy以及矩阵运算对于eigen。

所以合理的定位还是产品或者解决方案的角度，拿vtk里面的算法去解决某个实际的问题，形成一个具体领域的产品。从这个角度看，自己到底真正要实现什么或者项目书里怎么写就比较明确了。

一些可能的点应该就是ensmeble data 的vis，以及基于workflow或者llm的vis等等。

渲染pipeline

在vtk rendering pipeline的部分几个经常用到的对象包括 render window, vktRender, actor 以及 mapper。其中data set 是mapper的输入，然后mapper又是actor的输入，actor会被添加到render中。

通常render window以及 vtkRender 这些只需要在整个渲染管线初始化的时候设置好就行了。一旦data set 也确定下来，mapper以及actor也就确定下来了，不用再重新设置了，最基本的意识是一些操作是要进行缓存，比如actor，render一般数据不变的话就没必要重新创建了。

之后再不同时刻的时候唯一变换的就是camera的parameter，这个时候渲染的整个操作都是在GPU上计算的了，速度就相当快。在测试中几百兆的数据离线渲染基本上都是几十ms就好了，第一次render的时候需要一些额外的开销，可能需要1s多。

比较需要注意的地方就是在每次render操作执行的时候，如果设置不太好可能会导致mapper被重新执行一次，造成额外的开销，比如我最近的一个程序里，在每次render的时候，不小心执行了vtk data set 的set activate scalar 操作，这样vtk就认为每次data set 都是新创建的（虽然实际上data set就一直没变过），这样每次改变了camera，mapper还是被重新执行了一遍，基本上是O(cell num)的时间复杂度，这就引入了许多额外的开销。

后来是通过查看modified time发现这个原因的，vtk中的这个设计还是挺不错的，因为渲染管线本身就比较复杂，通过modified time查看哪些对象被修改，然后根据修改的情况执行对应的渲染更新操作就ok了。

offscreen时候 colorbar 以及 axis 设定

这个code是AI生成的，但是感觉可能常用到，就放在这里了。主要是offscreen渲染的时候，如何把axis坐标轴也加上去。这里还是基于vtk的，具体实现的逻辑是用另外一层layer作为 ui layer，把axis加上去，axis要和主相机进行同步，具体同步的时候需要一些计算。主要是只计算旋转的方向，focal point以及distance不发生改变，需要提前把一些关键的参数归一化。

    def _setup_axis_offscreen(self):
        if hasattr(self, "_orientation_marker_set") and self._orientation_marker_set:
            return

        # 1. 创建坐标轴 Actor
        axes = vtk.vtkAxesActor()
        axes.SetTotalLength(1.0, 1.0, 1.0) # 轴本身的几何长度
        axes.AxisLabelsOn()
        
        # 2. 创建独立的轴渲染器
        self.axis_renderer = vtk.vtkRenderer()
        # Viewport 分别是 [x_min, y_min, x_max, y_max]，左下角占 15% 宽度
        self.axis_renderer.SetViewport(0.0, 0.0, 0.15, 0.15) 
        self.axis_renderer.AddActor(axes)
        
        # 3. 设置层级：主场景在 0，坐标轴在 1 (最上层)
        self.render_window.SetNumberOfLayers(2)
        self.renderer.SetLayer(0)
        self.axis_renderer.SetLayer(1)
        self.axis_renderer.SetBackgroundAlpha(0.0) # 透明背景
        self.axis_renderer.InteractiveOff()

        # 4. 给轴渲染器配一个独立的相机，不要 SetActiveCamera(主相机)
        axis_cam = vtk.vtkCamera()
        self.axis_renderer.SetActiveCamera(axis_cam)

        self.render_window.AddRenderer(self.axis_renderer)
        self._orientation_marker_set = True

    def _sync_axis_camera(self):
        """
        同步主相机的旋转角度到坐标轴相机，但保持坐标轴相机原地不动
        """
        if not hasattr(self, "axis_renderer"):
            return

        main_cam = self.renderer.GetActiveCamera()
        axis_cam = self.axis_renderer.GetActiveCamera()

        # A. 获取主相机的方向向量
        # 方向 = 位置 - 焦点
        p = np.array(main_cam.GetPosition())
        f = np.array(main_cam.GetFocalPoint())
        v_up = main_cam.GetViewUp()
        
        direction = p - f
        # 归一化，我们只关心旋转角度
        norm = np.linalg.norm(direction)
        if norm == 0: return
        direction = direction / norm

        # B. 重新计算轴相机的位置
        # 我们让轴相机始终距离原点 (0,0,0) 固定 5 个单位长度
        # 这样无论主模型缩放多大，坐标轴大小永远不变
        distance = 5.0 
        new_pos = direction * distance

        axis_cam.SetPosition(new_pos[0], new_pos[1], new_pos[2])
        axis_cam.SetFocalPoint(0, 0, 0) # 永远盯着原点转
        axis_cam.SetViewUp(v_up)
        
        # 自动重置剪切平面，防止坐标轴旋转时被切掉
        self.axis_renderer.ResetCameraClippingRange()

自己编译VTK时候的注意事项

conda 自带的vtk有时候不带MPI以及python，这个时候需要自己编译一下，cmake相关的知识就不详细介绍了，需要知道的就是把USE_MPI打开，然后python 相关的参数包括下面这些，如下是一个例子

Python3_EXECUTABLE=D:/miniforge/envs/ISP_py312/python.exe
Python3_INCLUDE_DIR=D:/miniforge/envs/ISP_py312/include
Python3_LIBRARY=D:\miniforge\envs\ISP_py312\libs\python312.lib

之后正常 configure 以及 generate 对应的配置

build的时候在windows环境下一个是可以通过sln文件直接打开visual studio 然后进行编译，还可以在终端通过命令行的方式，比如

building command:
cmake --build . --config Release --target ALL_BUILD -- /m:16

installing command (this will install eveything into the install dir):
cmake --build . --config Release --target INSTALL

需要注意的一个地方是，在额外使用build好的python wrapper的时候，即使是conda 环境或者是 python env 也不用单独把vtk的包再给装到当前的环境里面了，最稳妥的方式是直接在虚拟环境中设置python path，比如在windows环境下直接

set PYTHONPATH=D:\working_dir\cworkspac
e\install_vtk_python312_mpi\bin\Lib\site-packages;%PYTHONPATH%

这样能保证路径不会出问题。对于conda环境来说，它的目录组织方式有点像一个小型的linux系统，很多包都放在根libs下了，要是不用conda然后自己去设置各种依赖的话，很容易出问题。

issue of loading time

今天遇到一个有意思的问题，两个同样结构，同样大小的vtk数据在loading的时候时间差了好几倍，一直是以为代码有问题，后来一问LLM才发现，原来根本的问题是两个vtk数据的格式不一样，一个是binary的格式，一个是ASCII的格式，在读取的时候，如果是VTK binary的文件格式，数据的加载速度会快上很多。

surface extraction 还是包含内部信息

surface extraction就是提取表面的意思，但是有好多时候cfd网格划分产生的数据不是那么理想，就容易导致surface extraction的过程失败，最常见的就是一个点存了两份，这个时候网格可能看起来是相邻的，或者说是在物理上是相邻的，但是在实际的数据结构上是不相邻的。这个时候常见的操作就是clean to grid 把数据给转成vtu的。

最近还见到了一个情况就是命名是四面体网格，结果存的时候就按照六面体网格的格式给存了，这样在surface extraction判定的时候计算normal方向的时候还是按照六面体的公式去计算，就容易出现问题了。这种时候就需要手动去过滤一遍，过滤每个cell，然后弄一个ordered sort，把如果点有重复的就通过ordered sort给过滤出去，就是把整个的mesh给重建了一遍。然后点还是原先的点。

目前就遇到了这两种情况。