上一节课内容回顾 第八讲 三维标量场体可视化 体绘制与面绘制的比较 体光照模型 体光线跟踪法 体单元投影法 三维扫描变换 VTK中的体绘制技术

8．1 体绘制与面绘制的比较 在自然环境和计算模型中，许多对象和现象只能用三维体数据场表示，和传统的计算机图形学相比，对象体不再是用几何曲面或曲线表示的三维实体，而是以体素(voxel)作为基本造型单元．对于每一体素，不仅其表面而且其内部都包含了对象信息，这是仅用曲线或曲面等几何造型方法所无法表示的． 体绘制的目的就在于提供一种基于体素的绘制技术，它有别于传统的基于面的绘制，能显示出对象体的丰富的内部细节． 体绘制技术包括体数据的表示、操作和绘制三部分内容．

图8．1 两种绘制方法的各自过程

8．2 体光照模型 体光照模型的研究是进行体绘制的基础．当光线穿过体素与光线遇到一曲面时，会发生不同的光学现象．前者如光线穿过云层会发生吸收、散射等现象，后者如光线射到桌面上，有漫射、反射、透射等现象．不同的物理背景决定了体光照强度的计算与面光照强度的计算有着不同的模型和方法． 体光照模型就是研究光线穿过体素时的光强变化，将光线穿过体素时的物理现象用数学模型来描述．在目前的体绘制中，采用得较多的有源—衰减模型、变密度发射模型和材料分类及组合模型．

体光照模型提供了体数据中各数据点光照强度的计算方法，体绘制方法提供的是二维结果图像的生成方法． 8．3 体绘制方法分类和比较 体光照模型提供了体数据中各数据点光照强度的计算方法，体绘制方法提供的是二维结果图像的生成方法． 首先根据数据点值对每一数据点赋以不透明度值和颜色值；再根据各数据点所在点的梯度及光照模型计算出各数据点的光照强度，然后将投射到图像平面中同一象素点的各数据点的半透明度和颜色值综合在一起，形成最终的结果图像． 根据不同的绘制次序，体绘制方法目前主要分两类： 以图像空间为序的体绘制算法 以对象空间为序的体绘制方法．

（1）以图像空间为序的体绘制算法——光线跟踪 以图像空间为序的绘制算法是从屏幕上的每一象素点出发，根据设定的视点方向，发出一条射线，这条射线穿过三维数据场的体素矩阵，沿这条射线选择K个等距采样点，由距离某一采样点最近的8个体素的颜色值及不透明度值做三维线性插值，求出该采样点的不透明度值及颜色值．在求出该条射线上所有采样点的颜色值及不透明度值以后，可以采用由后到前或由前到后的两种不同的方法将每一采样点的颜色及不透明度进行组合，从而计算出屏幕上该象素点处的颜色值． 光线跟踪的主要步骤： For 每条光线 Do For 每个与光线相交的体素 Do 计算该体素对图像空间对应象素的贡献

（2）以对象空间为序的体绘制算法——单元投影法 该类算法首先根据每个数据点的函数值计算该点的不透明度及颜色值，然后根据给定的视平面和观察方向，将每个数据点的坐标由对象空间变换到图像空间．再根据选定的光照模型，计算出每个数据点处的光照强度．然后根据选定的重构核函数计算出从三维数据点光照强度到二维图像空间的映射关系，得出每个数据点所影响的二维象素的范围及对其中每个象素点的光照强度的贡献．最后将不同的数据点对同一象素点的贡献加以合成． 单元投影法的主要步骤： For 每一体素或单元 Do For 该体素在视平面投影区域内的每一象素 Do 计算象素点获得的光照强度

图8．15 光线投射法

图8．16 光线投射流程

图8．20 体单元投影

8．7 VTK中的体绘制技术: 1 体元道具:vtkVolume ----几何表面可视化中vtkActor volume = vtkVolume() volume.SetMapper(volumeMapper) volume.SetProperty(volumeProperty) 2 体元特性:vtkVolumeProperty (颜色和透明度) 灰度: SetColor (vtkPiecewiseFunction ) 彩色: SetColor (vtkColorTransferFunction) 标量的透明度:SetScalarOpacity (vtkPiecewiseFunction ) 梯度的透明度: SetGradientOpacity (vtkPiecewiseFunction ) volumeProperty = vtkVolumeProperty() volumeProperty.SetColor(colorTransferFunction) volumeProperty.SetScalarOpacity(opacityTransferFunction) opacityTransferFunction = vtkPiecewiseFunction() opacityTransferFunction.AddPoint(20,0.0) opacityTransferFunction.AddPoint(255,0.2)

3 体元Mapper: vtkVolumeMapper vtkVolumeRayCastMapper: colorTransferFunction = vtkColorTransferFunction() colorTransferFunction.AddRGBPoint(0.0,0.0,0.0,0.0) colorTransferFunction.AddRGBPoint(64.0,1.0,0.0,0.0) colorTransferFunction.AddRGBPoint(128.0,0.0,0.0,1.0) colorTransferFunction.AddRGBPoint(192.0,0.0,1.0,0.0) colorTransferFunction.AddRGBPoint(255.0,0.0,0.2,0.0) 3 体元Mapper: vtkVolumeMapper vtkVolumeRayCastMapper: 设置射线合成函数: SetVolumeRayCastFunction(vtkVolumeRayCastFunction) 设置梯度计算函数: SetGradientEstimator(vtkEncodedGradientEstimator) volumeMapper = vtkVolumeRayCastMapper() compositeFunction = vtkVolumeRayCastCompositeFunction() volumeMapper.SetVolumeRayCastFunction(compositeFunction) volumeMapper.SetInput(reader.GetOutput()) (2)vtkVolumeProMapper用于VolumePro卡

本节课内容 第九讲 矢量场可视化 概述 数据空间及转换 基本点图标 矢量场线和面的生成 质点跟踪 矢量场拓扑 VTK中矢量场可视化 EVS中矢量场可视化

9．1 概 述 可视化的过程是一个把实验数据或计算数据转换成人类视觉能理解的表示形式的过程，根据不同的原始数据和不同的显示要求，可视化有许多具体的方法，但整个可视化过程可以分成三个阶段，即数据处理、可视化映射和绘制。 与标量场相比，矢量场的最大不同点在于每一物理量不仅具有大小而且具有方向，这种方向性的可视化要求决定了它是与标量场完全不同的可视化映射方法。 可视化映射主要是为数据设计一种能充分表达相关信息的表示方法，这种具体的表示被称为图标(Icon)．具体地讲，映射中的图标被定义为一几何对象，映射的过程就是一个根据给定点的数据，通过改变图标的几何特点，如长度、角度等，或改变其属性，如颜色、不透明度等，对给定点进行编码的过程．

对于标量场，许多图形学中现有的技术都能直接应用．主要是因为标量数据可直接映射到许多现有的图形参数上，如温度场可直接映射到物体颜色上，显示出连续的温度场变化． 而对于矢量场，可直接应用的技术就相对少多了．虽然可以用箭头对矢量场进行映射，但结果往往不理想，过多的箭头往往导致图像杂乱无章，无法表示矢量场的连续变化．问题的实质并不在于矢量场的大小，而在于缺少行之有效的矢量映射方法．对于张量场，问题就更突出了．因而与目前标量场的研究方向不同，矢量场的研究主要是集中于矢量的映射表示上，希望寻找一种既能反映矢量大小方向又不易引起混乱的映射图标，具体包括进行逐点映射的基础图标和通过特征抽取表示矢量场整体信息的局部图标和全局图标． 局部图标和全局图标是在不牺牲数据基本特征的前提下抽取数据特征的拓扑结构，它减少了数据集所包含的信息量，保留了局部和全局的关键特征，被认为是目前矢量场和张量场可视化中一个可行的方法。

从符号学角度出发，图标其实是一种逻辑符号，它表示了数据与其解释之间的一种关系．符号与数据之间有三种相关方式：图标(icon)、索引(index)和符号(symbol)． 图标是基于数据与表示之间的一种相似性设计的，如表示分子结构的化学分子式，它受到数据的影响，与数据有一种动态的联系． 索引是根据数据的真实影响来表示对象，如表示时间的钟是一种索引表示．与图标相比，索引是一种因果对应关系，而图标主要是一种模拟表示． 符号与数据对象之间根据任意变换都能对应，如英文字母不论如何变换其读音是不变的． 事实上，在可视化映射中所指的图标包括了符号学中的图标和索引．最简单的矢量图标是箭头，箭头本身是一个图标，它定性地描述了矢量的特征：方向和大小．而箭头中表示方向和长度的是一种索引，与数据值之间是一种因果对应关系．

9．1．1 矢量场可视化 图标及其相关属性是映射讨论的主要问题．根据所表示的数据属性，图标分成基本图标、局部图标和全局图标三大类．基本图标是一种典型的逐点模拟矢量场的技术，而局部图标和全局图标是一种更抽象的表示．图9．1总结了现有矢量场中的图标映射技术． 图9．1 矢量场图标

图标根据其表示的具体形式又能分为点、线、面、体图标．最直接的矢量图标是逐点映射的基本图标，被称为基本点图标．对于每一采样点，用具有大小和方向的图标映射矢量的大小和方向．常见的有箭头、锥体、有向线段等等，这类图标有时被统称为刺状体(hedgehogs)图标(见图9．2)． ． 图9．2 常见基本点图标

对于大型数据集，由基本点图标逐点映射出所有矢量会导致显示图像杂乱无章，显示太少又不能准确地把握矢量场的变化情况，而且点图标缺乏连续性，无法表示出矢量场中点数据之间的连续分布情况．如流场中涡流点附近流速变化较大，在这些区域内很难用箭头表示出涡流的具体特征．从中也可以看到，同是用点图标，表示标量场就要简单容易得多了．选择正确合适的插值方式，可以显示出标量场内每点的标量大小，如要显示某一物体的温度场，将温度映射到颜色，过密的点表示并不会导致杂乱无章，反而会使整个图像更细致，这就是矢量场映射与标量场映射的不同点．也是矢量场映射的困难之处． 为了表示出矢量场内矢量分布的连续性，从许多物理实验中得到启发，引入了线图标的表示方法，例如有质点轨迹、纹线、流线等表示．准确地讲，流线应称为场线，它显示了某一时刻，与经过点矢量相切的一条有向线段，在流场中被称作流线，在磁场中被称作磁力线，在电场中又被称作电力线．

质点轨迹显示的是在一定的时间间隔内，某一质点(集)经过矢量场后留下的运动轨迹．纹线定义的是某一时刻矢量场中所有单元质点的运动轨迹．在稳定场中，这三者的显示图像是一样的，只有在不稳定场中才会不一样． 线图标映射的有效性在很大程度上依赖于初始点的选取，同一矢量场，初始点位置选择不同，得到的场线也肯定不同，表示出的矢量场结构也不一样，要准确地揭示出矢量场中一些关键区域的矢量分布情况，选择的初始点必须经过这些区域．有人经过统计，对于一个有1M(1024×1024×l024)数据点的不稳定矢量场，大约需通过l0万个初始点的跟踪才能显示出其准确的结构．如此众多的初始点，其得到的结果同样是杂乱无章的一堆线段，难以揭示出准确的结构信息．因此又导致出现了面图标． 基本面图标主要是场面，在流场中称之为流面。它是场线的一个直接扩充，也就是与矢量场中经过的每一点的矢量都相切的一张曲面．场线是初始点跟踪的结果，相应地场面就是初始线段跟踪的结果．

由以上介绍可知基本图标不论是点图标还是线面标、面图标，其映射是一种根据给定点矢量值直接映射的方法，反映出的是给定点矢量的大小和方向，因而被认为是一种逐点映射的方法，无法反映出该点附近或矢量场整体的结构信息． 局部图标表示的是矢量场的梯度，现有的局部矢量图标包括局部点图标，如临界点，和局部线图标，如流带、流管等表示． 所谓临界点是矢量大小为零的点．在矢量场可视化中临界点对于了解矢量场内部拓扑结构有着非常关键的作用，可以根据终止于临界点的流线集合来判定这些流线附近的流线走向，因为流线之间除临界点外是不会相交的．在实际计算中，是根据临界点矢量的Jacobian矩阵的特征值对该点附近局部区域进行分类，分出交点、聚点、马鞍点等情况。 由此可以判定出局部区域内的矢量分布．

在不稳定场中，由于矢量梯度不同，流体单元在沿流线运动时会发生两种局部的形变：弹性形变和刚体旋转。要显示出这两种局部形变导致出现了局部的线图标表示，即流带和流管．流带是由两相邻流线定义的一条狭窄的面，它们的宽度可以反映出矢量场的散度，而扭曲的程度反映出沿流带的旋度．将N条流带连接起来而形成的流管，也就是一N边多边形沿流线的运动轨迹．在流管中，边的旋转反映了流线的旋度和横向形变，其断面显示了流场的散度． 局部图标表示了某些特定点附近矢量的分布特征，它是以一种简化抽象的表示来描述矢量场的整体结构．显然用一个图标来描述出整个矢量场的结构是件非常复杂和困难的事．

全局图标的基本方法是将整个区域以交于临界点的分界线(2D)或面(3D)分割成若干个等同于稳定流的子区域，再加上接近物体的切向矢量，将这些点和线连接起来就形成了矢量场拓扑结构的框架表示，沿着分割线上的点进行跟踪能生成表示流体分割的拓扑表示． 要生成整体图标表示关键是决定分割的拓扑形式，现有的技术还无法从原数据集中自动推导并构造这些逻辑结构，还需要在用户辅助下构造整体拓扑结构，特别是不同的局部结构之间的连接，这也是全局图标难以实现自动化的主要困难． ．

9．1．2 张量场可视化 在许多计算或分析中，张量场数据集是另一类非常重要的数据类型，它的可视比技术由于其成分量更多就更为复杂和困难了，目前还缺少强有力的表示图标来映射张量场数据，张量场的可视化可以说是刚刚起步．一个二阶张量(3×3)就有9个分量，目 前的工作是面向最简单的张量：二阶实数对称张量．图9．3显示出现有的几种张量场图标表示． 图9.3张量场图标

二阶实对称张量意味着该张量有三个实特征值和特征矢量，这样对一个张量场的映射可等同于同时对三个正交的特征矢量场的可视化映射．问题转化为在同一图像中同时显示出每点三个正交矢量的分布情况． 要显示每点三个特征值、三个特征矢量这样一种数据分布类型，最直接的基本点映射是采用椭球张量图标(见图9.4(a)) ，也就是用三个相互正交的特征矢量表示椭球的三个轴的方向，相应的特征值表示为三个轴的长度．对椭球张量图标的另一改进的表示方法是圆柱十椭圆图标表示(见图9.4(b)) ，用圆柱体的朗向来表示出最大特征值对应的主持征方向，用与圆柱垂直的一椭圆的长短轴表示次特征值和最小特征值的大小及特征矢量的方向．这种表示较之椭球图标更突出了主特征方向的意义。 图9.4 张量基本点图标

张量基本点图标比矢量基本点图标更复杂，当数据量增大时，分布不连续和杂乱天章的问题就更突出了．为表示出张量分布连续性同样也出现了张量基本线图标：超流线，当具有一定大小的几何单元沿某一特征矢量扫过张量场时，得到一条轨迹，在每一 经过点，该几何单元受到另两个特征矢量的变换操作，沿轨迹将变换后的几何单元连接生成一曲面，该曲面被称为超流线．超流线的面可以用一用户定义的函数对其颜色编码。 一般是用特征值的大小直接编码．超流线的颜色和轨迹完全表示了某一特征矢量场的分布，而其断面表示了另两个特征矢量场，这样整个张量场就可用超流线图标表示，同时也可以看到，在超流线图标映射时的两个关键问题：超流线轨迹生成和断面类型构造．

超流线的轨迹生成与流线相仿，断面上不同几何单元的选择将影响着超流线断面上另两个特征矢量的编码，常见的几何单元有圆和十字叉．由圆作为几何单元生成的超流线被称为管(tube)，由十字叉生成的超流线被称为螺旋线(helix)。管和螺旋线各有特点，前者适合于表示在断面上的张量变化，而后者又适合于表示主方向的改变． 与矢量场类似，张量场的全局图标是对许多超流线综合性质的映射，也就是主要特征的抽取映射．从轨迹临界点、断面奇点等点中显示出张量场的全局结构，当然比矢量的全局图标更为复杂．与标量场相比，由于在显示空问方向上的难度，矢量场和张量场的许多问题都有待于进一步研究，关键问题仍是难于表示出空间大量矢量数据或张量数据的方向性，张量场的难度就更大了，从使用点图标到线图标和面国标，映射方法的改变改进了基本信息的显示手段，用于显示数据场梯度的局部图标和揭示数据场整体结构的全局图标的研究提供了一种减少矢量场显示的数据信息量，把握关键特征的拓扑结构显示方法．由于拓扑结构一般是用曲面(线)表示，因而拓扑结构比原矢量场更易显示和理解．

9．2 数据空间及转换 尽管矢量场或张量场数据既能来自于实验也能来自于计算，但就目前可视化的研究现状而言，可视化主要还是作为计算仿真的视觉再现，也就是说，其数据主要还是来源于计算．因此，有必要从计算角度出发讨论数据空间以及对可视化技术的影响，特别是在矢量场中． 从计算仿真到图像显示，整个可视化过程涉及到三个数据空间，即物理空间、计算空间和图形空间． 物理空间(P空间)：是现实世界对象的存在空间，一般也是物理方程直接定义的空间．对其中对象进行网格化离散处理后，通常其网格是曲线网格，网格与具体对象边界一致，物理量可以在网格点上直接计算． 计算空间(C空间)：是为了简化数值计算的需要而提出的空间，其中的网格通常是正交网格。正是由于其网格简单，物理方程在计算空间就更易迭代求解．P空间与C空间之间的网格单元是相互对应的，通过转换两个空间是相互等价的．

图形空间(G空间)：是进行可视化图形处理的空间，经过映射，数据场被转换成图形空间中的几何表示．G空间的几何表示有曲面、实体、层次等表示形式，通过各种造型手段，G空间所要表示的就是现实世界的物理对象，与P空间是同一对象．而P空间网格离散化一般已拥有较高的精度，没有必要为了图形绘制再对P空间中的对象重新离散化．从可视化映射的角度出发，要重构必然要对物理量进行再采样。结果往往会导致欠采样或过采样，从而使可视化结果的精度降低，所以G空间中的对象一般就直接采用了P空间中的表示形式，两者是一致的．当然最终G空间被表示成象素，即可视化结果图像． 这些不同的空间之间是可以相互转换的，围绕着同一物理实体，其表示能在不同的空间之间进行转换．由于一般G空间与P空间一致，因而这种转换主要是讨论P空间与C空间之间的转换(见图9．5)

图9．5 物理空间与计算空间的相互转换

作为数值计算预处理的第一步是网格离散化，能生成许多类型的网格：结构或非结构，直线或曲线，或几种组合．本讲所讨论的矢量场映射限定在结构化网格，网格单元的拓扑是一正规六面体，每一网格单元以下标(i，j，k)标识，它是网格单元中与网格原点最近的顶点．P空间中的点一般写成Xp(x，y，z)，其矢量写成Vp(i，j，k)＝(u，v，w)． 假设计算空间C是正交的笛卡儿网格，与P空间网格单元具有相同单元下标(i，j，k)和点Xc(ξ,η,ζ)．在C空间网格点的矢量可写成 Vc(i，j，k)＝(u’，v’，w’) 一般而言，P空间与C空间之间：不存在全局的转换函数关系，只存在网格单元(i，j，k)及其邻域在P空间与C空间的局部转换关系．非网格点的转换可用转换函数插值求取． 设L是网格单元(i，j，k)从C到P的局部转换，该转换就是有限元中常讲的形函数(见图9．6)。类似地，L-1就是从P空间到C空间的转换(见图9．5)

图9．6 单元间转换

实现这种相互转换是很必要的．这是因为不同的矢量场映射会在不同的空间中进行，箭头等基本点图标的映射一般是在物理空间中直接进行，而流线的计算一般是在计算空间中进行，这主要是从映射计算的简便出发的． 在可视化映射的计算中经常遇到的一个问题是网格空间中的点定位问题．给定一点坐标，找出该点所在网格单元下标(i，j，k)：和偏移量(α,β,γ)．如果点是在C空间中一点Xc(ξ,η,ζ) ，且C空间是正交网格，则找出(i，j，k)是很方便的．困难的是从P空间中一点Xp(x，y，z),找出C空间中对应单元(i，j，k)和偏移量(α,β,γ)．通常采用的是循环逼近的方法。 这种数据空间中的定位主要用于矢量场中质点轨迹的跟踪，但在具体实现中，前一跟踪点可作为初始估计点，这样估计的命中率大为提高．在非结构化网格单元中，初始估计点的准确性就更重要了．

9．3 基本点图标 矢量场的基本点图标是最简单的映射图标，其中最常见用得最多的是箭头，其它还有锥体、有向线段等多种表示，所有这些点图标也被称作刺状体(hedgehog)． 9．3．1 箭头 二维中的箭头可直接由线段组成，尽管简单，但对于表示二维矢量场仍不失为一种有效的表示，一般要求，按大小放大后的箭头不要相互覆盖，这样一幅二维矢量场图还是清晰的． 对于三维矢量场，箭头的绘制就要复杂多了．在没有有向光反射、光照模型等提供三维深度信息时只有通过透视变换才是唯一方法，因而三维线箭头在方向性的解释上往往会出现模糊，观察者很难从图9．7中区分哪一解释是正确的．

图9．7 三维箭头方向的模糊性

箭头的大小由三个因子决定: 矢量的大小、相对与图像平面的矢量方向和透视变换．为了提高方向性，减少杂乱的程度，对箭头的映射和显示提出了许多方法． 考虑一个由边界面包围的区域(见图9．8)，区域R内每一点有一矢量值V，对其中断面上的每点P采样有(Vx，Vy，Vz)三个分量，若干个断面上的矢量场组成了R中整个矢量场． 采用实体箭头表示R中的矢量场，其映射有下列属性： 1)由矢量的大小决定箭头的组成：头和尾的大小和长度； 2)箭头不透明； 3)箭头能旋转交互．

图9．8 定义矢量场的区域R和其上的采样平面

整个箭头由头和尾组成(见图9.9)，头是一矩形金字塔，尾是一条长为l的线段． 令L是箭头的总长，即矢量的大小，h＝l＝L/2，头的底部矩形大小为 a＝pL，b＝qL (0＜p，q＜1) 参数p，q可控制箭头的形状，常取p＝1/2，q＝1/8． 角β是矢量与断面的夹角，角α是矢量在断面上的投影与坐标轴的夹角，角αβ定义了箭头的方向，同时也决定了箭头中其它线段的可见性．

图9．9 三维箭头表示

在箭头绘制时，还应注意以下，几个问题： 1)箭头的长度应小于网格单元的大小，以免引起混乱．如果箭头方向在邻域内有反转的情况．其长度就应小于网格大小的1／2； 2)为了加强箭头的三维特性，增加旋转交互。通过旋转了解矢量场整体结构，避免因脱化、遮挡等引起的模糊性； 3)除了用旋转等交互手段后，表示深度信息还可用阴影，在光强计算中绘制出阴影效果. 4)可采用正交投影，而不用透视变换．这是因为在透视变换下不同位置相同大小的矢量，看上去会不一致，易引起二义性； 5)箭头也可用其它各种形式，如图9．11所示．为了区别朝向，在线段的一端上要有一些特定的标志．据研究，拥有可变宽度的线段更易为人所感知．

图9．11 各种箭头表示形式

9．4 矢量场线和面的生成 9．4．1 矢量线的生长算法 在许多研究上，经常用场线来提高对矢量场数值解的理解．所谓场线，是一条虚拟曲线，其上任一点的方向与矢量场在该点的方向一致．场线在流体速度场中称为流线，在电场中称为电力线等等。 考虑一个三维矢量场，P是其中一点，其位置为r．通过P的场线是一条空间曲线，r作为t的函数，t可以是时间、弧长等多种参量． 按定义场线为(见图9．12) 求解该方程，能构造出 某一瞬时的一条场线 图9．12 场线定义

设场线为 是方程的解，则在物理空间场线的积分表达式为 上式表明场线是从一初始位置开始，以小步长t生成，r(0)是 初始条件．

算法9．1 场线生成算法 令rijk是矢量fijk的单元角点，其中i＝1，…，Ni，j＝1,…,Nj，k=1,…,Nk，Cijk是单元号(rijk是Cijk的原点)，计算空间坐标为(ξ，η，ζ)． 1)计算出场线初始点r0(ξ0，η0，ζ0)； 2)计算出网格映射形函数Lijk的参数； 3) 计算出八个单元顶点在计算空间的矢量场gs=Js-1*fs； 4)用gs的三线性插值找出采样点矢量值，对方程式(9．4)积分，用Runge—Kutta法解移动质点，求出场线下一新的位置Lijk (ξ，η，ζ) ，连接新老位置，直到场线到达Cijk的某一面，如 ξ= ξ* =1，η= η* ，ζ= ζ* 5)如果面是边界面(如i=Ni-1)，则停止, 否则 i＝i+1， ξ = -1， η= η* ，ζ= ζ* go to 1)继续算法

用数值方法生成流线存在着以下几个误差源： 1)跟踪算法所无法控制的误差，包括收敛性误差，截断误差，高度扭曲网格的误差； 2)与跟踪算法有关的误差，如矢量场插值，方程数值积分的阶数，步长选择等导致的误差．特别是在高曲率网格的情况下，最常用的三线性插值的误差较大．其主要原因是在三线性插值模式下，流体质量不守恒，产生的误差靠改进积分方法也难以消除，特别是在曲率较大时．Yeung使用了基于三次样条维持流体质量守恒的高阶插值模式，精度提高了，但计算更费时，难以达到交互生成的程度．

8．4．2 流函数构造矢量线 为了解决场线生长模型的误差，提高场线的精度，Kenwright提出了一高精度的求连续场线的方法．该方法以流场为基础，将一单元内离散定义的流体转换成由二个三维流函数表示的流体，流线就是该两流函数的交线．跟踪每一组流函数的一个常数解就能生成一条流线．这一跟踪过程是质量守恒的，无需按步长积分，避免了效值积分导致的误差. 图9．14 流面和流线

首先，流线被定义为两独立面f和g算的一条交线(见图9．14)，设 f＝ f(x，y，z) g＝ g(x，y，z) 对于f或g的任一常数解，上两函数分别表示了一张曲面，如果 这些函数是流线方程的积分解，其定义的曲面就称为流面．在流面上没有流体穿过，总是与流场相切，这两个面的交线也总是与流场相切，也就是流线. 与前节中的流线生长算法相比，流函数方法更有优势．首先速度加快，免去了生长模型中解微分方程的Runge—Kutta数值方法；其次精度提高了；再者流函数本身是一个很好的流面构造工具，可用等值面方法生成．更重要的是提供了一种思路，将矢量场转化为标量场处理的方法．

8．4．3 矢量面的生长算法 在二维场景中，场线技术对于给出矢量场的走势是很有效的，对于矢量场的大小也能给出间接的解释，这是因为场线间的距离与矢量场的大小是成反比的． 但是在三维场景中，场线技术就不如在二维场景中那么有效了．主要是需要深度信息来确定场线的空间位置，但场线往往较细，如果直接计算光照强度，效果并不十分理想．因此，Hultquist提出若干条场线插值连接生成矢量面的构造技术，在流场中被称为流面(stream surface)．若只使用两条相邻场线，生成的面被称为流带(stream ribbon)．

矢量面是一无限可伸展线段经过矢量场所形成的二维轨迹．一般的矢量场较复杂，分叉很多，嵌入的曲面往往变形严重(扭曲、重叠)，所以构造矢量面比构造等值面要更困难． 由场中的矢量线组成的面称为矢量面，对应于矢量面上每一点(X)，该点的矢量V(X)必在过该点的矢量面的切面上．例如在所考虑的矢量场内任取一异于矢量线的曲线c，并过c的每一点引出一矢量线，则这些矢量线的几何轨迹给出了一个矢量面．当曲线c为封闭时，所得的矢量面是管状的，称为矢量管，在流场中被称为流管(stream pipe)． 因此，基于矢量线的构造，矢量面的构造是类似的．在矢量场中选定一初始质点集合R组成曲线c上的点，由R中的每一初始点生成一条矢量线，将相邻两矢量线间的空间用三角片连接生成矢量面．

图9．17 矢量场的(s，t)网格定义 图9．18 矢量面连接

矢量面可用(s，t)参数空间上的st矩形网格构造(见图9．17)，在每一网格上有(si，ti) ，其中si是矢量线，是初始点集R在矢量场中的矢量线. ti是矢量线的参数，可取为时间,表示了在ti时刻R中所有质点所在的位置． 由于矢量面有时卷曲很严重，R中质点的采样分布应根据曲率的变化再采样以作调整． 根据两条相邻矢量线上的采样点生成三角片时，可采用类似等值面的连接方法． 如L，R是两条流线，如图9.18所示，L0，R0是在t0时刻两质点位置，L1，R1是在t1时刻两质点位置，可采用连接较短者的方法．即比较L1R0和L0R1的距离，连接其中较短的一条对角线．这样每一条带上包含了左、右两条流线上的点．除了s0，s1外，其它矢量线都被两条带所共享．

9．5．1 基本质点跟踪 9．5 质点跟踪 质点跟踪技术来源于流场的实验，在流场中将烟雾、染料等注入流场，观察其运动轨迹的分布，以了解流场的内部结构．该技术被广泛地应用于流场实验．质点跟踪技术是数值矢量场中的模拟技术，质点被认为是一点状发光的运动质点，质点轨迹是一条曲线，是由该质点在矢量场中一系列连续的位置组成． 质点的运动方程是:

其中X是在t时刻质点的位置，V(X)是在矢量场中X处的矢量值，主要是速度场． 质点跟踪的过程就是确定质点在物理空间中一系列具体位置的过程．从初始位置X(x，y，z)开始，在计算空间中跟踪．具体步骤为 1) 在计算空间中为位于物理空间X的质点定位，求出单元(i，j，k)和偏差(α，β，γ)； 2)在网格单元内，通过插值求出该点的矢量值； 3)解微分方程，找出下一质点位置．

质点的形状最常见的是单独的点，采用点可显示出矢量场整体的结构．但质点还可采用其它一些特殊类型的几何单元，如二维中的圆和三维中的球，但这些单元和点一样无法显示出朝向，即使使用有向点光源，球也只能显示出其位胃形状，无法决定其空间 的朝向．而作为矢量场，显尔其中的朝问是很必要的．因而为f提供空间的朝向，可选择平面的或拉伸形状的对象作为质点形状． 为了改进静态图像中的方向性效果，每个质点可加上一个尾，尾显示了前一时间片质点的运动轨迹，其长度表明速度值．另一方法是采用流箭头(stream arrow),箭头的头指向质点位置，表明了运动方向，其尾部是运动轨迹的一部分(见图9．22)． 图9．22 流箭头

9．6 矢量场拓扑 矢量线和矢量面的构造都存在着一个问题，即初始点位置的选择，准确的选择对于了解矢量场内部结构是非常重要的．但这种选择却是一件非常困难的事情，尤其是在不知道矢量场内部结构的前提下． 从1987年起由Helman和Hesselink开始的流体拓扑结构的分析和可视化研究为矢量场的可视化提供了一种从全局了解矢量场结构的新技术，它是一种全局图标的构造和映射方法．该技术可以从流场扩展到其它矢量场中，如涡流场、压强梯度场等． 矢量场的分析和可视化由以下几步组成： 1)临界点位置的计算； 2)对临界点进行分类； 3)计算积分曲线或曲面．

9．6．1 临界点及其分类 临界点是矢量三个成分同时为零的点．在矩形网格中其准确位置可以由插值计算出来．如果是曲线网络，临界点的位置可以通过递归分割网格单元计算出来，或通过Newton下降法等数值方法计算出来． 由于在临界点速度矢量值为零，即V＝0，因此临界点附近的矢量场完全由偏导 所决定．临界点就是根据 的特征值和特征向量决定分类。 求解在临界点位置矢量的Jacobian矩阵：

Jacobian矩阵特征值实部的正或负分别表示了相吸和相斥的特征．具体地讲，正的特征值表示V从临界点发散(出)，负的特征值表示V从临界点聚扰(入)，共扼复数对表明了V是螺旋入或出(由实部符号决定其入或出)． 因此，根据 待征值，临界点被分成交点、聚点和马鞍点三类．交点和聚点又可进一步分成吸引的和排斥的．

1 图9.25 二维临界点分类表

在完成临界点分类后，要通过积分曲线(面)将临界点连接起来，构造矢量场的拓扑结构．这些积分曲线(面)以临界点作为起始点和终止点，以临界点的特征向量作为起始点的方向．在包含对象物体的流场中，在无滑动边界上还存在另一类V＝0的点，这类点称为壁点．在具体构造积分曲线时，壁点加上临界点成为起始点和终止点的所有点集．为了便于控制，要求起始点只有有限条线段交于该点．这样由马鞍点、壁点组成了起始点，而由其它点组成终止点．积分曲线从起始点开始，马鞍点共有4条，而壁点有一条，由这些曲线描述了流场结构．

小结: 在矢量场可视化上，介绍了传统的基本点图标的矢量场映射方法，包括刺状体(箭头等有向线段)，场线，场面(流带、流管)等生成算法，以及质点跟踪技术。 在线图标和面图标的映射上，介绍了二维和三维矢量场拓扑结构的分析和抽取算法．

（1） 矢量场的流线表示 vtkStreamLine sl = vtkStreamLine() sl.SetIntegrator(integ) 设置积分函数 sl.SetInput(aPolyDataV) 设置矢量数据 sl.SetSource(seeds.GetOutput()) 设置流线种子点 sl.SetMaximumPropagationTime(1000.0) sl.SetIntegrationStepLength(0.1) sl.SetIntegrationDirectionToIntegrateBothDirections() sl.SpeedScalarsOn() sl.SetStepLength(0.1) 生成线段vtkLineSource seeds = vtkLineSource() seeds.SetPoint1(1430, 14200.0, 25000.0) seeds.SetPoint2(1430, 8200.0, 25000.0) seeds.SetResolution(30) 积分函数vtkRungeKutta4 integ = vtkRungeKutta4()

pl3d = vtkPLOT3DReader() pl3d.SetXYZFileName("./combxyz.bin") pl3d.SetQFileName("./combq.bin") pl3d.SetScalarFunctionNumber(100) pl3d.SetVectorFunctionNumber(202) pl3d.Update() rake = vtkLineSource() rake.SetPoint1(15,-5,32) rake.SetPoint2(15,5,32) rake.SetResolution(21) integ = vtkRungeKutta4() sl = vtkStreamLine() sl.SetInput(pl3d.GetOutput()) sl.SetSource(rake.GetOutput()) sl.SetIntegrator(integ) sl.SetMaximumPropagationTime(0.1) sl.SetIntegrationStepLength(0.1) sl.SetIntegrationDirectionToBackward() sl.SetStepLength(0.001) mapper = vtkPolyDataMapper() mapper.SetInput(sl.GetOutput()) actor = vtkActor() actor.SetMapper(mapper) 演示

（2） 矢量场的流点表示 vtkStreamPoints sl = vtkStreamPoints() 流点生成对象 integ = vtkRungeKutta4() sl.SetInput(pl3d.GetOutput()) sl.SetSource(rake.GetOutput()) 设置矢量数据 sl.SetIntegrator(integ) sl.SetMaximumPropagationTime(0.1) sl.SetIntegrationStepLength(0.1) sl.SetIntegrationDirectionToBackward() sl.SetTimeIncrement(0.001) cone = vtkConeSource() 流点表示类型 cone.SetResolution(8) cones = vtkGlyph3D() 装配对象 cones.SetInput(sl.GetOutput()) cones.SetSource(cone.GetOutput()) cones.SetScaleFactor(0.001) cones.SetScaleModeToScaleByVector()

integ = vtkRungeKutta4() sl = vtkStreamPoints() sl.SetInput(pl3d.GetOutput()) sl.SetSource(rake.GetOutput()) sl.SetIntegrator(integ) sl.SetMaximumPropagationTime(0.1) sl.SetIntegrationStepLength(0.1) sl.SetIntegrationDirectionToBackward() sl.SetTimeIncrement(0.001) cone = vtkConeSource() cone.SetResolution(8) cones = vtkGlyph3D() cones.SetInput(sl.GetOutput()) cones.SetSource(cone.GetOutput()) cones.SetScaleFactor(0.001) cones.SetScaleModeToScaleByVector() mapper = vtkPolyDataMapper() 演示 mapper.SetInput(cones.GetOutput())

（3） 矢量场的流面表示 vtkRuledSurfaceFilter 首先生成流线， 然后用vtkRuledSurfaceFilter生成流面（两条流线生成一流面）： scalarSurface = vtkRuledSurfaceFilter() scalarSurface.SetInput(sl.GetOutput()) 设置流线 scalarSurface.SetOffset(0) 设置显示的第一条流面 scalarSurface.SetOnRatio(2) 设置不显示的流面比率 scalarSurface.PassLinesOn() 是否显示流线

integ = vtkRungeKutta4() sl = vtkStreamLine() sl.SetInput(pl3d.GetOutput()) sl.SetSource(rake.GetOutput()) sl.SetIntegrator(integ) sl.SetMaximumPropagationTime(0.1) sl.SetIntegrationStepLength(0.1) sl.SetIntegrationDirectionToBackward() sl.SetStepLength(0.001) scalarSurface = vtkRuledSurfaceFilter() scalarSurface.SetInput(sl.GetOutput()) scalarSurface.SetOffset(0) scalarSurface.SetOnRatio(2) scalarSurface.PassLinesOn() scalarSurface.SetRuledModeToPointWalk() scalarSurface.SetDistanceFactor(30) mapper = vtkPolyDataMapper() mapper.SetInput(scalarSurface.GetOutput()) mapper.SetScalarRange(pl3d.GetOutput().GetScalarRange()) 演示

EVS中矢量场可视化

第三套 作业： 作业要求： 一 创建交互器和交互窗口 a 创建1个交互器， b 创建1个绘制窗口， c 在上述绘制窗口中，创建3个绘制器。 二 读入第7讲例子中所附的数据RectGrid.vtk 创建vtkContourFilter对象，在第一个绘制器中以等值面方式显示RectGrid.vtk数据。 要求： 只显示RectGrid.vtk的等值值是0.5, 1.0, 2.5，3.0的等值面. 三 创建vtkStreamLine对象，在第二个绘制器中绘制数据v252000.dat.vtk的流线 三 创建vtkStreamPoints对象，在第三个绘制器中绘制数据v252000.dat.vtk的流点 作业参考： 参考第7讲例子中的iso-plane2.py streamline.py streampoints.py kkmodel.py

下一节课内容： 第十讲 三维交互技术 三维交互工具 三维交互算法 立体图绘制 交互视算