一种基于流线重心Voronoi图的流场可视化方法

摘要

本发明公开了一种基于流线重心Voronoi图的流场可视化方法。本发明将点的重心Voronoi图的概念推广到流线上，以流线重心Voronoi图为基础给出一种有效的基于流线的流场可视化方法。本发明既可以实现更佳效果的流线布置，又能同时兼顾到流线的覆盖性、均匀性、连续性，并且可以迭代的完成流线布置。本发明的具体步骤为：1）随机生成流线的种子点集合；2）由种子点出发，以数值积分的方法追踪得到长流线；3）计算长流线的Voronoi区域；4）确定并极小化长流线的目标函数，优化目标函数，进而得到新的种子点集合，返回步骤2），直至流线在流场中均匀分布。

说明书

技术领域

本发明涉及一种流场可视化方法，具体涉及一种基于流线重心Voronoi图的流场可视化方法。 

背景技术

流场可视化是科学可视化中非常重要的一部分，特别是在计算力学，流体力学和空气动力学、天气预报和地球物理模型等领域有着广阔的应用前景。流场按照流速是否随着时间变化分为定常流场和非定常流场，按照数据采集维度分为二维流场、三维表面流场、三维空间流场三类。一般的，不同的流场类型需要不同的流场可视化算法。归纳地，流场的可视化技术可以分为五类： 

（1）直接可视化：数据不需要提前处理，直接进行可视化，但是由于能力有限，该类方法往往造成视觉上的混乱。 

（2）基于纹理的可视化：该类方法一般用于二维流场和三维平面流场，但是他们存在关键点定位问题。 

（3）基于特征的可视化：该类方法更适合大规模数据流场，因为它更关注特征点信息。 

（4）基于空间剖分的可视化：基于空间剖分的方法采用空间或积分线的相似性聚类等技术对可视区域进行剖分，可能存在结果的可信性验证等问题。 

（5）基于几何的可视化：该类方法利用几何信息进行流场的可视化，具有表示清晰、简洁等优点。 

流场可视化是科学可视化中非常重要的组成部分，在流体力学、空气动力学、天气预报和地球物理模拟等领域有着广阔的应用前景。基于几何的流场可视化因其能够反映流场几何特性的优点是一类重要分支，其中，流线方法是基于几何的可视化方法中最常用的技术。 

流线是基于几何的可视化方法中最常用的技术，流线即是处处与流场相切的曲线，能有效地表现流场的局部特征。恰当的流线放置可以清晰，简洁的反映流场的特点。评价流线布置的质量，一般根据以下原则： 

（1）覆盖性：流线应涵盖整个域，不会错过任何重要的特征。 

（2）均匀性：流线应均匀地分布在流场中。 

（3）连续性：较少的长流线效果要优于众多的短流线。 

目前，基于流线的算法可以分为两类： 

基于图像的算法：将流线布局看成一幅二值图像，通过比较当前布局与一幅参考图像的差别来优化流线位置。如Turk和Banks定义了一个基于当前布局的低通滤波与一副均匀灰度图像之差的能量函数，通过一系列操作如组合、删除、生成、增长、缩短等使得函数能量下降进而优化得到一组较短的流线（G.Turk and D.Banks,“Image-guided streamline place-ment,”SIGGRAPH’96:Proceedings of the 23rd annual conference on Computer graphics and interactive techniques，“基于图像的流线放置方法”，第23届计算机图形和交互技术会议，1996，453-460）。然而，由于每一步都可能需要大量的操作，这种方法会出现指数的运行时间。采样控制算法：该方法通过不同的采样原则直接控制种子点的生成。如Jobard和Lefer基于一个叠加的正规网格计算种子点与流线的距离，从而选择种子点的位置（B.Jobard and W.Lefer,“Creating evenly-spaced stream-lines of arbitrary density,”in In Proc.of 8thEurographics Workshop on Visualization in Scientific Computing，“基于任意密度的流线均与分布的流场可视化”，第8届科学计算可视化欧洲研讨会；“Unsteady flow visualization by animating evenly-spaced streamlines,”Computer Graphics Forum,vol.19,no.3,pp.31–39,2000.[Online].，“基于均与间隔的非定常流场可视化”计算机图形学论坛，第二卷，2000，31–39）。这个方法可以很好保证流场的均匀性，但是由于使用短流线而忽略了流线的连续性。 

Rosanwo等人利用与流场相切以及正交的两组对偶流线进行可视化（O.Rosanwo,C.Petz,S.Prohaska,H.-C.Hege,and I.Hotz,“Dual streamline seeding”in Visualization Sym-posium,2009.Pacific Vis’09.IEEE Pacific，“双流线算法”，可视化研讨会，2009。PacificVis’09.IEEEPacific,April,pp.9–16），需要首先计算流场的拓扑骨架以获得初始对偶流线。 

以上方法都可以满足覆盖性、均匀性、连续性的部分要求，但是几乎所有的算法都是通过一次迭代确定性的完成流线布置，无法进一步优化流线的位置。Du等人利用点的重心Voronoi图给出了一个流场可视化的算法（Q.Du and X.Wang,“Centroidal voronoi tessellation based algorithms for vector fields visualization and segmentation”“基于CVT的流场可视化和分割算法”，2004，可视化，IEEE，10月，第43-50页），但是以箭头为可视元素，效果不甚理想。 

发明内容

本发明的目的是为克服上述现有技术的不足，提供一种基于流线重心Voronoi图的流场可视化方法。本发明将点的重心Voronoi图的概念推广到流线上，以流线重心Voronoi图为基础给出一种有效的基于流线的流场可视化方法。本发明既可以实现更佳效果的流线布置，又能同时兼顾到流线的覆盖性、均匀性、连续性，并且可以迭代的完成流线布置。 

为实现上述目的，本发明采用下述技术方案： 

一种基于流线重心Voronoi图的流场可视化方法，具体步骤为： 

1）随机生成流线的种子点集合； 

2）由种子点出发，以数值积分的方法追踪得到长流线； 

3）计算长流线的Voronoi区域； 

4）确定并极小化长流线的目标函数，优化目标函数，进而得到新的种子点集合，返回步骤2），直至流线在流场中均匀分布。 

所述步骤1）中的种子点集合为其中n为大于1的整数。 

所述步骤2）的具体步骤为： 

21）选取种子点集合中的第一个种子点x_i，并获取其坐标p_j及速度V(p_j)； 

22）种子点x_i根据欧拉公式p_k+1＝p_k+δV(p_k)同时向前向后积分，从而获得长流线上的所有采样点s[x_i]={p₁,…,p_j，…,p_m},其中δ表示积分步长，1≤k≤m，且j、k、m均为正整数。 

所述步骤3）的具体步骤为： 

31）根据点的Voronoi区域计算公式 

$\mathrm{Vor}\left(x_p\right)=\{x\in R^N|\left|\right|x-x_p\left|\right|\le \left|\right|x-x_q\left|\right|,\forall q\ne p\}$计算长流线上每个采样点的Voronoi区域； 

32）用长流线上的点{p₁,…,p_m}代替每条长流线s[x_i]，其中x_i为s[x_i]的种子点，用p_j＝x_i来表示它，即将长流线的Voronoi区域表示为 

$\mathrm{Vor}\left(s\right[x_i\left]\right)\cong \underset{k}{\cup }\mathrm{Vor}\left(p_k\right),1\le k\le m.$

所述步骤4）的具体步骤为： 

41）目标函数为$\tilde{F}\left(X\right)=\underset{i\in S\left[X_i\right]}{\Sigma }\underset{k}{\Sigma }f\left(p_k\right);$

42）利用拟牛顿法极小化目标函数，具体步骤如下： 

当目标函数的梯度时，其中，ε为收敛判断条件，且ε>0， 

421)求解梯度下降方向d，使之满足

422)确定步长α以满足值下降； 

423）对流线的种子点进行一个αd的位移，即X←X+αd； 

其中和分别表示目标函数的梯度和Hessian矩阵； 

43）利用L-BFGS算法框架（D.C.Liu and J.Nocedal,“On the limited memory BFGS method for large scale optimization,”Mathematical Pro-gramming“有限内存大规模优化的BFGS框架”，数学专业编程：系列A和B，第一卷。第503-528页）优化求解目标函数。 

所述步骤43）中，L-BFGS算法框架的具体步骤是，对于每个牛顿迭代： 

431)对于每条长流线s[x_i]，计算其采样点

432)对于每个采样点p_k，计算其裁剪Voronoi区域，即Vor(p_k)∩Ω，其中Ω表示给定区域； 

433)将每个采样点p_k的贡献加到和中。 

本发明的有益效果： 

1.本发明利用长流线对流场进行可视化，可以迭代的完成流线的布置，在初始化流线之后，可以进一步优化流线的位置； 

2.本发明使用的长流线可以很好地满足流线的连续性特点； 

3.由于流线重心Voronoi图本身的优良特性，本发明能够很好的满足流线的均匀性和覆盖性两个特点。 

附图说明

图1a为7个点的Voronoi图； 

图1b为7个点的重心Voronoi图； 

图2a为流场数据； 

图2b为图2a流场中的两条流线的Voronoi图； 

图2c示出了图2b中的流线Voronoi图由流线上采样点近似计算而成； 

图2d为图2c的局部放大图； 

图3a为流场数据； 

图3b为本发明流线算法初始化效果图； 

图3c为图3b中流线算法迭代2次之后效果图； 

图3d为图3b中流线算法迭代6次之后效果图； 

图4为本发明的流程框架图。 

具体实施方式

下面结合附图和实施例对本发明进行进一步的阐述，应该说明的是，下述说明仅是为了解释本发明，并不对其内容进行限定。 

本发明提供了一种基于流线重心Voronoi图的流场可视化方法，图4为根据本发明的实施例的流场可视方法的流程框架图，步骤如下： 

1）随机生成流线的种子点集合，如图中401-402； 

2）由种子点出发，以数值积分的方法追踪得到长流线，如图中403； 

3）计算长流线的Voronoi区域，如图中404； 

4）确定并极小化长流线的目标函数，优化目标函数，进而得到新的种子点集合，返回步骤2），直至流线在流场中均匀分布，如图中405。 

所述步骤1）中的种子点集合为其中n为大于1的整数。 

所述步骤2）的具体步骤为： 

21）选取种子点集合中的第一个种子点x_i，并获取其坐标p_j及速度V(p_j)； 

22）种子点x_i根据欧拉公式p_k+1＝p_k+δV(p_k)同时向前向后积分，从而获得长流线上的所有采样点s[x_i]={p₁,…,p_j,…,p_m},其中δ表示积分步长，1≤k≤m，且j、k、m均为正整数。 

流线是某一时刻在流场中画出的一条空间曲线，在该时刻，曲线上所有质点的速度矢量均与这条曲线相切。流线是基于欧拉方法描述流场的一种技术手段，在实际使用中，常用多线段（polyline）的形式来表示，即由种子点出发，用数值积分的方法跟踪得到流线的点序列。但是如何恰当的选择积分步数并不是一件简单的事情。由于长流线可以更好的表现出流线的连续性，本发明中的算法选择使用长流线作为表现流场的几何元素。算法中的长流线当遇到流场边界或速度为0的点时，会停止积分。 

以向前积分为例，欧拉积分法表示为：p_k+1＝p_k+δV(p_k)，其中，1≤k≤m，且k、m均为正整数。V(p_k)代表点p_k处的速度，δ代表积分步长。为了算法后半部分方便计算梯度，将流线的种子点x_i以及其流场方向V_i和正交向量作为参考标架，将采样点表示为该标架中的向量，即其中x_i表示种子点x_i的位置信息，α_k、β_k是点p_k的标架参数，$V_i=V\left(x_i\right)=\left(\begin{array}{c}{V}_{\mathrm{ix}}\\ V_{\mathrm{iy}}\end{array}\right),$$V_i^{\perp }=\left(\begin{array}{c}-V_{\mathrm{iy}}\\ V_{\mathrm{ix}}\end{array}\right).$

所述步骤3）的具体步骤为： 

31）根据点的Voronoi区域计算公式 

$\mathrm{Vor}\left(x_p\right)=\{x\in R^N|\left|\right|x-x_p\left|\right|\le \left|\right|x-x_q\left|\right|,\forall q\ne p\}$计算长流线上每个采样点的Voronoi区域； 

32）用长流线上的点{p₁,…,p_m}代替每条长流线s[x_i]，其中x_i为s[x_i]的种子点点，用p_j＝x_i来表示它，即将长流线的Voronoi区域表示为 

$\mathrm{Vor}\left(s\right[x_i\left]\right)\cong \underset{k}{\cup }\mathrm{Vor}\left(p_k\right),1\le k\le m.$

因为长流线用多线段表示，所以多线段s[x_i]的Voronoi区域可以表示为 $\mathrm{Vor}\left(s\right[x_i\left]\right)=\{z\in R^N|d(z,s[x_i\left]\right)\le d(z,s[x_j\left]\right),\forall j\ne i\in S\},$d(z,s[x_i])表示点距离多线段的欧拉距离。但是一般的，采用近似计算的方式计算长流线的Voronoi图，原因如下： 

1)多线段Voronoi图的计算复杂性：两个线段的Voronoi中分线是二次曲线，生成多线段的Voronoi图涉及二次曲线求交等计算，代价极高。 

2)即使得到了多线段的Voronoi图计算，对目标函数进行优化时在Voronoi单元边界上进行二次函数的积分也非常困难。 

因此用流线上的点{p₁,…,p_m}来代替每条流线s[x_i]。这里x_i是s[x_i]的种子点，用p_j＝x_i来表示它。因为对给定包含n个种子点的点集x_i的Voronoi区域定义为：$\mathrm{Vor}\left(x_i\right)=\{x\in R^N|\left|\right|x-x_i\left|\right|\le \left|\right|x-x_j\left|\right|,\forall j\ne i\}.$所以，流线的Voronoi单元可以用生成该流线的所有积分点序列的Voronoi单元的并集近似表示 $\mathrm{Vor}\left(s\right[x_i\left]\right)\cong \underset{k}{\cup }\mathrm{Vor}\left(p_k\right).$

图2给出了一个流线重心Voronoi图的近似计算示例。其中，对于(a)中的流场，(b)中显示的是两条流线（曲线）的Voronoi图（中间直线段），(c)中给出了利用采样点近似计算得到的各条流线的Voronoi图，(d)中显示的是(c)中图的局部放大效果。 

所述步骤4）的具体步骤为： 

41）目标函数为$\tilde{F}\left(X\right)=\underset{i\in S\left[X_i\right]}{\Sigma }\underset{k}{\Sigma }f\left(p_k\right);$

42）利用拟牛顿法极小化目标函数，具体步骤如下： 

当目标函数的梯度时（ε为收敛判断条件，ε>0）， 

421)求解梯度下降方向d，使之满足

422)确定步长α以满足值下降； 

423）对流线的种子点进行一个αd的位移，即X←X+αd； 

其中和分别表示目标函数的梯度和Hessian矩阵； 

43）利用L-BFGS算法框架优化求解目标函数。 

所述步骤43）中，L-BFGS算法框架的具体步骤是，对于每个牛顿迭代： 

431)对于每条长流线s[x_i]，计算其采样点

432)对于每个采样点p_k，计算其裁剪Voronoi区域，即Vor(p_k)∩Ω，其中Ω表示给定区域； 

433)将每个采样点p_k的贡献加到和中。 

给定包含n个种子点的点集x_i的Voronoi区域定义为 $\mathrm{Vor}\left(x_i\right)=\{x\in R^N|\left|\right|x-x_i\left|\right|\le \left|\right|x-x_j\left|\right|,\forall j\ne i\};$CVT可以通过变分的形式定义，其对应能量函数： 

$F\left(X\right)=\underset{i=1}{\overset{n}{\Sigma }}f\left(x_i\right);f\left(x_i\right)={\int }_{\mathrm{Vor}\left(x_i\right)\cap \Omega }{\left|\right|x-x_i\left|\right|}^2\mathrm{dx}---\left(1\right)$

的极小解。这里，我们假设Ω的密度分布为常数。 

如图1展示了点的普通Voronoi图和重心Voronoi图的区别，图1中空心点代表生成点，实心点代表对应Voronoi区域的重心。 

因为流线的Voronoi图可以由所有流线生成点的Voronoi图近似表示。所以流线重心Voronoi图能量函数是目标函数： 

$\tilde{F}\left(X\right)=\underset{i\in S}{\Sigma }\underset{k}{\Sigma }f\left(p_k\right)---\left(2\right)$

的极小解。 

通过极小化流线重心Voronoi图能量，从而达到优化流线布置，直至在流场中均匀分布。为了极小化目标函数（请见公式（2）），利用拟牛顿法求解。由于计算Hessian矩阵以及其求逆相当耗时，采用L-BFGS算法，即不计算二阶导数，而是通过利用以前的梯度和修改值去近似Hessian矩阵的逆矩阵。L-BFGS算法具有近似的算法结构，不同之处在于在步骤1）中Hessian矩阵被替换为由以前的梯度值构造的简单矩阵。在实际中，我们可以采用优 化工具包TAO（T.Munson,J.Sarich,S.Wild,S.Benson,and L.C.McInnes,“TAO 2.0 users manual,” Mathematics and Computer Science Division, Argonne National Laboratory, Tech. Rep.ANL/MCS-TM-322，“TAO 2.0用户手册”，数学集成电路和计算机科学部，美国Argonne国家实验室）或科学计算工具包PETSc（S.Balay,J.Brown,,K.Buschelman,V.Eijkhout,W.D.Gropp,D.Kaushik,M.G.Knepley,L.C.McInnes,B.F.Smith,and H.Zhang,“PETSc users manual,” “PETSc用户手册”，美国Argonne国家实验室，Tech.Rep.ANL-95/11–Revision 3.2,2011）等提供的L-BFGS算法实现。 

和的计算过程如下： 

对于一个三角形T＝(p_k,q₁,q₂)，其中，p_k,q₁,q₂为T的三个顶点，它所关联的重心Voronoi函数能量可以通过下式计算得到：其中U_i＝q_i-p_k，|T|表示三角形T的带符号面积。重心Voronoi函数的梯度可以通过下式计算得到 其中，m_k和c_k分别表示三角形T的质量和质心，有 $m_k={\int }_{\mathrm{Vor}\left(p_k\right)\cap \Omega }\mathrm{dx},$$c_k=\frac{1}{m_k}{\int }_{\mathrm{Vor}\left(p_k\right)\cap \Omega }\mathrm{xdx}.$基于链式法则，可以通过下式计算： 

$\frac{\partial \tilde{F}}{\partial x_i}=\underset{k}{\Sigma }\frac{\partial f\left(p_k\right)}{\partial x_i}=\underset{k}{\Sigma }\frac{\partial f\left(p_k\right)}{\partial p_k}\frac{\partial p_k}{\partial x_i}=2\underset{k}{\Sigma }{m}_k\frac{\partial p_k}{\partial x_i}(p_k-c_k)$

流线上的采样点相对于种子点的梯度如下： 

$\frac{\partial p_k}{\partial x_i}=1+{\alpha }_k\left(\begin{array}{cc}\frac{\partial V_x}{\partial x}& \frac{\partial V_y}{\partial x}\\ \frac{\partial V_x}{\partial y}& \frac{\partial V_y}{\partial y}\end{array}\right)+{\beta }_k\left(\begin{array}{cc}\frac{-\partial V_y}{\partial x}& \frac{\partial V_x}{\partial x}\\ \frac{-\partial V_y}{\partial y}& \frac{\partial V_x}{\partial y}\end{array}\right),$其中α_k、β_k是点p_k的标架参数。 

图3展示了本发明算法的迭代过程。图3(a)为流场数据信息；图3(b)为流线初始化后的结果；图3(c)为算法迭代执行2次之后的结果；图3(d)为算法迭代执行6次之后的结果。实验证明，我们的算法在10次迭代之内就可以达到稳定的状态。 

而且由于流线重心Voronoi图本身的特点，本发明同时兼顾到了均匀性和覆盖性。 

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。