CN104318591A - 一种带边界平面流场的动态绘制方法 - Google Patents

Abstract

本发明公开一种基于图像融合技术的带边界平面流场的动态绘制方法。该方法首先确定平面流场的可视化范围并对该范围进行栅格化处理，根据给定的边界曲线确定流场的绘制区域；然后确定流场显示图像的大小并构造背景噪声图像集；根据流场和边界计算当前帧各栅格顶点速度的矢量，运用图像融合技术和纹理映射技术构造当前帧流场图像。通过逐帧构造流场图像，实时绘制流场的动态运行过程。相对于现有技术，本发明方法经济、视觉效果好、适应性强，适用于各种静态与动态流场的可视化。

Description

一种带边界平面流场的动态绘制方法

技术领域

本发明提供一种基于图像融合技术的带边界平面流场的动态绘制方法，其属于可视化技术领域，内容涉及图像融合技术及流场可视化技术，尤其涉及一种带边界平面流场的动态绘制方法。

背景技术

流场在科学计算和工程分析中扮演非常重要的角色，流场的可视化是理解流场运动变化规律的有力工具，在科学计算、工程分析、教学科普等领域中发挥着重要作用，具有重要的实际意义和研究价值。

点图标法(参见：王成恩,姜传彦,武新宇.基于几何图形的矢量场可视化技术.东北大学学报(自然科学版),2011,32(3),327-330)和矢量线法(参见：胡星,杨光.流线可视化技术研究与进展.计算机应用研究,2002,(05):8-11)是传统的流场可视化方法的代表，这两种方法都要求有较高的采样频率，否则可能会丢失流场的一些重要细节，而采样点过多又会导致图像混乱，并且这两种方法都不适用于对动态流场进行可视化。

基于纹理的方法中，点噪声方法(参见：Van Wijk J J.Spot noise:texture synthesis for datavisualization.Computer Graphics,1991,25(4):309-318)是最早的基于纹理的流场可视化方法，该方法沿流场方向对噪声滤波来生成流场图像，依靠改变点的属性，可整体或局部地控制所产生的纹理模式。该方法不适合绘制流场中流线方向变化剧烈的局部区域。LIC方法(参见：CabraL B,Leedo M C.Imaging vector fields using line integral convolution.ACM Computer Graphics,1993,27(4):263-270)用一维卷积核函数沿流线方向对噪声纹理进行卷积计算，生成可以表达流场模式的静态图像，该方法计算耗时，容易引入高频噪声，且很难反映流场动态强度。

刘庆平等申请的专利“流场可视化方法及装置”(申请号：200710055770.6)使用照相机、光源、遮光板、聚光镜和盛粒子器皿等装置实现流场观测与分析。该方法需要制作特别的硬件装置，其实现成本较高。

发明内容

本发明针对现有技术的不足，给出一种带边界平面流场的动态绘制方法，是一种基于图像融合技术的流场绘制方法。本方法即可以产生静态流场图像，也可以模拟流场动态运行过程，适用于各种类型的流场可视化绘制。

为了叙述方便，对相关的概念、符号做如下说明：

1.Ω＝[a,b]×[c,d]：即平面可视化区域，其宽为a～b，高为c～d，围成此区域的四条线段称为Ω的边界。

2.F(x,y；t)：表示流场，(x,y)为平面点坐标，t为时间。F(x,y；t)的值为点(x,y)处在t时刻的场值。当场值与t无关时，F(x,y；t)为稳定场；否则，为非稳定场。

3.C＝{c_k；k＝0,1,2,...}：区域Ω中的曲线集合，k为曲线标号。C中各曲线可以嵌套但不相交。若c_k为不封闭曲线，则其与Ω的局部边界构成封闭曲线。附图2为曲线示意图。

本发明采用的技术方案如下：

一种带边界平面流场的动态绘制方法，具体包括以下步骤：

A.确定可视化范围并进行栅格化处理：确定可视化范围Ω＝[a,b]×[c,d]中a,b,c,d的值，将Ω划分为N_Y行N_X列的四边形网格；计算栅格单元的宽DX和高DY，以及各栅格顶点坐标(x,y)；

B.确定流场的可视化边界：

预处理工作：将C中边界曲线c_k转化成由有序点列构成的折线；形成的折线表示为其中n+1为点列中点的个数，表示点列中第m个点；

确定与每条折线l_k对应的流场可视化边界：计算折线l_k与栅格线的交点，记为其中i为交点标号，i＝0,1,...；各交点形成沿折线l_k的有序点列，顺次连接这些交点，形成一条折线，将形成的折线设定为与l_k对应的流场可视化边界；

C.确定流场图像大小并构造背景噪声图像集，技术方案如下：

1)设定背景噪声图像为灰度图像；确定流场图像和背景噪声图像的宽、高；设定噪声图像集中图像幅数N；产生L个随机的图像灰度值，并存入环形有序序列A中；将A中元素的序号称为背景噪声图像像素的相位，计算变量t(t＜L)，称t为相位差；

2)构造第一幅背景噪声图像：对于各像素e，随机取A中的元素，将此元素值作为e的灰度值，并将此元素的序号作为e的相位；

3)构造其余背景噪声图像：对于第k(k＝2,3..,N)幅噪声图像中的像素e，记第k-1幅图像中e的相位为ph，取A中ph后第t个元素，将该元素值作为e的灰度值，并将该元素的序号作为e的相位；

D.计算各栅格顶点的速度矢量，技术方案如下：

1)由F(x,y；t)，计算各栅格顶点的速度(vx,vy)；

2)调整与l_k对应的流场可视化边界的速度矢量：记q_k为l_k与栅格线的任一交点，调整该点速度方向为向量q_kq_k+1的方向；记q_k所在的栅格线段为p₀p₁，将p₀,p₁中距q_k较近的点作为边界栅格顶点，并调整其速度方向与q_k速度方向相同；

3)以边界栅格顶点为界，逐列(或行)对栅格顶点分类：以逐列方式为例，对于某列栅格线上的各边界栅格顶点p_k(k＝0,1,...,n)，设定p_2k～p_2k+1间的栅格顶点为“流场边界外栅格顶点”，p_2k+1～p_2k+2间的栅格顶点为“流场边界内栅格顶点”；

4)调整不同类型栅格顶点速度：首先，逐行(或列)调整流场边界内栅格顶点速度，记相邻两边界栅格顶点p₀,p₁速度矢量与X轴正向夹角为α₀,α₁，计算p₀p₁间栅格单元个数n及角度变量Δα＝(α₁-α₀)/n；以Δα为递增量，依次重新设定p₀p₁间各栅格顶点的速度矢量与X轴正向的夹角α，并按此角度调整其速度方向；然后，置流场边界外栅格顶点的速度为0；

E.构造单帧流场图像，将第一幅背景噪声图像作为第一帧流场图像，构造后续单帧流场图像的技术方案如下：

1)构造临时流场图像：计算前一帧时刻各栅格顶点的速度矢量，前一帧流场图像各像素按前一帧时刻对应位置处栅格顶点的速度移动，形成新图像LI，称为当前帧临时流场图像；

2)图像融合：取背景噪声图像集中与当前帧对应的图像g，将图像LI和图像g融合，融合后的结果即为当前帧流场图像；

F.绘制流场动态运行过程：用上述方法构造单帧流场图像，按时间顺序，每构造一帧图像，将其绘制出来，如此，显示流场的动态运行过程。

本发明方法相对于现有技术，具有如下优点：

1.本发明方法可以产生各时刻静态流场图像，也可以通过绘制连续帧流场图像显示流场的动态运行过程。

2.本发明方法适用于带一条或多条(不交叉)任意边界线的平面流场，同时也适用于不带边界的平面流场。

3.稳定和非稳定流场都可以用本发明方法可视化。

附图说明

图1本发明方法的流程图；

图2可视化范围内曲线示意图；

图3参数曲线转化为折线示意图；

图4隐式曲线转化为折线示意图；

图5可视化范围示意图；

图6可视化范围栅格化示意图；

图7求折线相邻两点所成线段与栅格线交点示意图；

图8对折线相邻两点所成线段与栅格线交点排序示意图，(a)为折线相邻点所成的线段与栅格线的交点示意图，(b)为对折线相邻点所成的线段与栅格线的交点进行排序的示意图；

图9调整流场边界栅格顶点速度矢量示意图；

图10构造背景噪声图像集示意图；

图11计算栅格顶点速度矢量示意图；

图12调整流场边界栅格顶点速度矢量示意图；

图13确定流场边界范围内、外区域示意图；

图14调整流场边界范围内栅格顶点速度矢量示意图；

图15构造单帧流场图像示意图。

具体实施方式

结合相关概念和附图，对本发明进行详细说明。如附图1所示，一种带边界平面流场的动态绘制方法，其具体实施步骤如下：

A.确定可视化范围并进行栅格化处理，步骤如下：

1)设定a,b,c,d的值，确定Ω。本实例取a＝-b＝-1，c＝-d＝1(参见附图5)。

2)将Ω区域栅格化为N_Y行N_X列的四边形网格，称栅格线交点为栅格顶点，本实例取N_Y＝N_X＝100。根据Ω的大小及N_X与N_Y，计算栅格单元的宽和高DX,DY：

$\left\{\begin{array}{c}\mathrm{DX}=2\·b/N_X\\ \mathrm{DY}=2\·d/N_{\text{Y}}\end{array}\right.---\left(1\right)$

进而计算各栅格顶点的坐标P_ij(x_ij,y_ij)(0≤i≤N_Y,0≤j≤N_X)(i,j为栅格线的行列号)：

$\left\{\begin{array}{c}{x}_{\mathrm{ij}}=-b+i\·\mathrm{DX}\\ y_{\mathrm{ij}}=-d+j\·\mathrm{DY}\end{array}\right.---\left(2\right)$

附图6是大小为[a₁,b₁]×[c₁,d₁]，6行8列的平面栅格。

B.确定流场的可视化边界：

预处理工作，即将C中边界曲线转化成由有序点列构成的折线l_k：

Ⅰ.若c_k是参数曲线，参数为t，参数取值范围为a≤t≤b，则c_k表示为：

$\left\{\begin{array}{c}x=f_x\left(t\right)\\ y=f_y\left(t\right)\end{array}\right.(a\≤t\≤b)---\left(3\right)$

t在参数范围内取不同的值，得到一组有序点列，连接有序点列形成与c_k对应的折线，表示为其中n+1为点列中点的个数，表示点列中第m个点(参见附图3)。

Ⅱ.若c_k是隐式曲线，可表示为f(x,y)＝c，其中x_a≤x≤x_b,y_a≤y≤y_b，c为常数。参见附图4，将曲线c_k转化为折线l_k的步骤为：

1)设定变量len，其值根据曲线中x的范围确定，实例中取len＝|(x_b-x_a)/20|。

2)从x＝x_a开始，以len为步长，查找第一个满足f(x,y)＝c的点，记为将其放入点列。

3)依次查找其余点，并将其放入点列：

为表示l_k的点列中的当前点，计算下一点时有两种情况：

a)若曲线c_k在点处可导

首先，计算c_k在处的导数，规定该点处切向量长度为len，进而计算该点处的切向量V_m，并记其与X轴正向夹角为α，由此计算x_m+1的值x_m+1＝x_m+len·cos(α)。

然后，将x_m+1代入f(x,y)＝c，得到满足曲线的各点p_k(x_m+1,y_k)(k＝0,1,...)，并计算各纵坐标差值Δy_k＝y_k-y_m。判断α，若0≤α≤π，取Δy_k＞0中Δy_k最小的点，即为若π＜α＜2π，取Δy_k＜0中-Δy_k最小的点，即为

b)若c_k在点(附图4用表示)处不可导

设定变量ε＜len，实例取ε＝len/2，计算横坐标x_m+1＝x_m+ε，随即将x_m+1代入f(x,y)＝c，得到纵坐标y_m+1，即得到点(附图4用表示)。

按如上步骤，计算出表示折线l_k的点列中的所有点。

确定与l_k对应的流场可视化边界的步骤为：

1)求线段与栅格线的交点，步骤为：

1.1)计算各点所在的栅格单元M_IJ：记点所在的栅格单元行列号分别为I,J，用表示不大于Ζ的最大整数，计算I,J值的公式为：

1.2)计算交点：记点和点所在的栅格单元行列号分别为i₀,j₀和i₁,j₁。线段与横向栅格线交点记为p_i(x_i,y_i)，则y_i＝-d+i·DY(i为栅格线行号，i₀≤i≤i₁)；与纵向栅格线交点记为q_j(x_j,y_j)，则x_j＝-b+j·DX(j为栅格线列号，j₀≤j≤j₁)(参见附图7)。y_i,x_j的计算公式为：

$\left\{\begin{array}{cc}{x}_i=x_0-k_1\·(y_i-y_0)& (k_1=(y_1-y_0)/(x_1-x_0))\\ y_j=y_0-k_2\·(x_j-x_0)& (k_2=(x_1-x_0)/(y_1-y_0))\end{array}\right.---\left(5\right)$

将求得的交点依次存入有序序列List。

2)对交点排序：若点p(x_p,y_p),q(x_q,y_q)满足|x_p-x_q|≤DX且|y_p-y_q|≤DY，称点p,q为临近点。对List中的点排序的方法为：首先，记点为P，从List中找到P的临近点，将其与List中第一个元素交换，并记该点为P；然后，从List其余元素中找到P的临近点，与P后面的第一个元素交换，并记该点为P，如此，循环执行此操作，直至将所有点排序(参见附图8)。

3)确定流场的可视化边界：按如上方法计算点列中每相邻两点所成线段与栅格线的交点，得到l_k与栅格线的所有交点，顺次连接这些交点，形成的折线即为与c_k对应的流场的可视化边界(参见附图9左图)。

C.确定流场图像大小并构造背景噪声图像集：设定背景噪声图像为灰度图像。流场图像的宽高记为W,H，背景噪声图像的宽高记为BW,BH，背景噪声图像像素所在的行列号用r,c表示，则0≤r＜BH,0≤c＜BW。

将各背景噪声图像编号为k＝0,1,...N-1，产生L个(实例取L＝256)0～255间的随机整数，并存入环形有序序列A中(参见附图10)；将A中元素的序号称为背景噪声图像像素的相位，A_m表示A中第m个元素。构造背景噪声图像集的步骤如下：

1)设定流场图像和背景噪声图像大小。分别设定W,H,BW,BH的值，须满足BW:BH＝W:H，记背景噪声图像总幅数为N。噪声图像各灰度值取自A，记相邻两噪声图像相同位置处像素的相位差为t。本实例取W＝H＝512，BW＝BH＝64，N＝32，计算

2)构造第一幅背景噪声图像：首先，设定r行c列处像素的相位(0为第一幅背景噪声图像编号)为0～L-1间的随机整数；然后，设定此处像素的灰度值为序列A中序号为的元素的值(参见附图10)。

3)构造其余背景噪声图像：依次构造编号为k(1,2..,N-1)的背景噪声图像，步骤为：首先，设定r行c列处像素的相位随后，将元素的值作为该像素的灰度值(参见附图10)。

D.计算各栅格顶点的速度矢量，步骤如下：

1)将栅格顶点坐标P_ij(x_ij,y_ij)代入流场F(x,y；t)，求得顶点速度矢量V_ij(vx_ij,vy_ij)，并计算矢量长度v_ij及矢量与X轴正向的夹角α_ij：

本实例设定流场为稳定流场F(x,y；t){vx＝vy＝0.1,t≥0}，则对于各栅格顶点α_ij＝π/4(参见附图11)。

2)调整流场可视化边界的速度矢量：按如下步骤调整l_k与栅格线各交点的速度方向：如附图9所示，记q_k为l_k与栅格线交点中的任一点，线段p₀p₁为q_k所在的栅格线段；记q_k,q_k+1的坐标为q_k(x_k,y_k)，q_k+1(x_k+1,y_k+1)；记Δx＝x_k+1-x_k，Δy＝y_k+1-y_k；首先，计算向量q_kq_k+1与X轴正向的夹角α：

然后，计算q_k与点p₀,p₁的距离a,b，并判断：若a≤b，将点p₀作为流场边界栅格顶点，若a＞b，将点p₁作为流场边界栅格顶点。接着，记流场边界栅格顶点所在的行列号分别为I,J，令α_IJ＝α，由此计算该顶点速度矢量V_IJ(vx_IJ,vy_IJ)：

$\left\{\begin{array}{c}{\mathrm{vx}}_{\mathrm{IJ}}=v_{\mathrm{IJ}}\·\cos \left({\mathrm{\α}}_{\mathrm{IJ}}\right)\\ {\mathrm{vy}}_{\mathrm{IJ}}=v_{\mathrm{IJ}}\·\sin \left({\mathrm{\α}}_{\mathrm{IJ}}\right)\end{array}\right.---\left(8\right)$

附图12为调整流场边界栅格顶点速度矢量示意图。

3)逐列对栅格顶点分类：记某列流场边界栅格顶点依次为p_k(k＝0,1,...,n)，其中p₀,p_n为该列栅格线与Ω边界的交点，设定p_2k～p_2k+1间的栅格顶点为“流场边界外栅格顶点”，p_2k+1～p_2k+2间的栅格顶点为“流场边界内栅格顶点”(参见附图13右图)。附图13中阴影部分为流场边界范围外区域，反之，为流场边界范围内区域。

4)调整不同类型栅格顶点速度，步骤为：

4.1)逐行调整流场边界内栅格顶点速度方向：如附图14所示，记某行栅格线上相邻两流场边界栅格顶点为p₀,p₁，其速度与X轴正向夹角为α₀,α₁，依次将p₀到p₁(包括p₀,p₁)间的栅格顶点标记为k＝0,1,...，并计算p₀到p₁间的栅格单元个数n，计算角度插值变量Δα＝(α₁-α₀)/n，随后，依次计算这些栅格顶点p_ij的速度矢量与X轴正向的夹角α_ij＝α₀+k·Δα，按此夹角调整点p_ij的速度方向，大小不变；计算点p_ij的速度矢量V_ij(vx_ij,vy_ij)的公式为：

$\left\{\begin{array}{c}{\mathrm{vx}}_{\mathrm{ij}}=v_{\mathrm{ij}}\·\cos \left({\mathrm{\α}}_{\mathrm{ij}}\right)\\ {\mathrm{vy}}_{\mathrm{ij}}=v_{\mathrm{ij}}\·\sin \left({\mathrm{\α}}_{\mathrm{ij}}\right)\end{array}\right.---\left(9\right)$

4.2)处理流场边界外栅格顶点速度，置vx_ij＝0,vy_ij＝0,v_ij＝0,α_ij＝0。

E.构造单帧流场图像。将第一幅背景噪声图像作为第一帧流场图像。记当前时刻t对应的流场为第k(k＝2,3,...)帧流场；记相邻两帧间的时间间隔为Δt，t-Δt时刻的流场记为第k-1帧流场，实例取Δt＝1。参见附图15，构造第k帧流场图像的具体步骤如下：

1)构造临时流场图像：

1.1)保存t-Δt时刻各栅格顶点坐标记为则：

$\left\{\begin{array}{c}{\mathrm{tx}}_{\mathrm{ij}}^k=x_{\mathrm{ij}}^{k-1}\\ {\mathrm{ty}}_{\mathrm{ij}}^k=y_{\mathrm{ij}}^{k-1}\end{array}\right.---\left(10\right)$

1.2)计算t-Δt时刻各栅格顶点速度矢量，记为并计算t时刻各栅格顶点坐标，记为则：

$\left\{\begin{array}{c}{x}_{\mathrm{ij}}^k=x_{\mathrm{ij}}^{k-1}+{\mathrm{vx}}_{\mathrm{ij}}^{k-1}\·\mathrm{\Δt}\\ y_{\mathrm{ij}}^k=y_{\mathrm{ij}}^{k-1}+{\mathrm{vy}}_{\mathrm{ij}}^{k-1}\·\mathrm{\Δt}\end{array}\right.---\left(11\right)$

将移动后的各栅格顶点形成的栅格记为M。

1.3)将第k-1帧流场图像作为纹理，以为纹理坐标，用纹理映射技术将其映射到栅格M，形成的图像记为第k帧临时流场图像LI_k。

2)图像融合：取背景噪声图像集中第k％N幅图像g_k，设定变量β(0＜β＜1.0)，实例取β＝0.5，将图像LI_k和图像g_k按β:(1-β)融合，融合后的结果作为第k帧流场图像I_k＝β·LI_k+(1-β)·g_k。

附图15中，a为第k-1帧流场图像，b为第k帧栅格，c为第k帧临时流场图像，d为与k帧流场图像对应的背景噪声图像，e为第k帧流场图像。

F.绘制流场动态运行过程：通过流场F(x,y；t)，按照上述方法，每构造一帧流场图像，将其绘制出来。如此，通过实时生成并绘制有序帧流场图像的方式，呈现流场的动态运行过程。

Claims (1)

1.一种带边界平面流场的动态绘制方法，其特征在于，具体步骤是：

A.确定可视化范围并进行栅格化处理：确定平面流场可视化区域的范围，并将其栅格化为N_Y行N_X列的四边形网格，计算各栅格顶点的坐标；

B.确定流场的可视化边界：对于给定的每条边界曲线，求出该曲线与栅格线的交点；用边界曲线上的交点序列定义的折线，表示相应的可视化边界曲线；

C.确定流场图像大小并构造背景噪声图像集：设定流场图像和背景噪声图像的分辨率，以及背景噪声图像的幅数N；产生L个随机的图像灰度值，并存入环形有序序列A中；构造N幅有序的背景噪声图像，其中第一幅图像各像素的灰度值均随机取自序列A，第一幅图像与其余图像位于相同位置处像素的灰度值在序列A中呈周期性变化，且周期相同；

D.计算各栅格顶点的速度矢量：由给定的流场计算各栅格顶点速度矢量；沿可视化边界曲线方向适当调整曲线与栅格线交点的速度方向；在与交点相邻的两栅格顶点中，将离交点较近的顶点设定为流场的边界栅格顶点，并调整此顶点速度方向与交点速度方向相同；以边界栅格顶点为界，将栅格顶点划分为“流场边界内栅格顶点”和“流场边界外栅格顶点”两类，并适当调整这些栅格顶点的速度；

E.构造单帧流场图像：将第一幅背景噪声图像作为第一帧流场图像；后续每一帧流场图像，是由前一帧流场图像和前一帧流场混合后，再与对应当前帧的背景噪声图像融合而形成的；

F.绘制流场动态运行过程：每构造一帧流场图像，将其绘制出来，如此，用逐帧构造并实时绘制流场图像的方式，显示流场的动态运行过程。