CN1873640A - 扩展的矢量场线卷积积分方法 - Google Patents

Abstract

扩展的矢量场线卷积积分LIC方法，(1)从矢量场中寻找下一个未被标记的采样点；(2)以上述找到的采样点为起点，采用基于数值计算的DDA曲线形成方法形成DDA曲线，首先从起点出发进行曲线的前向预测，然后还是从起点出发进行曲线的后向预测，最后将前向预测得到的曲线和后向预测得到的曲线合并，形成最终的DDA曲线；(3)对此轮循环形成的DDA曲线进行有色噪声序列直接赋值，计算DDA曲线上各点的LIC结果；(4)依据此轮循环形成的DDA曲线对整个矢量场进行区域标记，依据该曲线上的矢量场采样点的幅度大小确定覆盖区域，标记该曲线区域中的所有矢量场采样数据点；(5)对矢量场的LIC值进行归一化处理。本发明该方法既可以面向非结构矢量场，又可面向结构矢量场，增强了矢量幅度及矢量场结构的表现力，提高了预处理计算速度。

Description

扩展的矢量场线卷积积分方法

所属技术领域

本发明涉及一种扩展的矢量场线卷积积分LIC方法，属于科学计算可视化领域。

背景技术

在科学计算可视化研究领域，矢量场可视化的研究目的在于运用计算机提供的图形、图像、视频等多种手段将科学实验、仿真、测量得到的大量矢量场数据准确、有效地展现出来，并支持人机间的多种交互处理。矢量场可视化的方法首先必须保证矢量场信息描述的准确性和完整性；其次还应具备较高的计算效率和较强的针对矢量场结构及幅度的表现力；矢量场可视化方法中的交互处理主要指多角度及多分辨处理。为了达成上述目标，矢量场可视化的发展大致经历了三个阶段，第一个阶段是基于图标(例如小的箭头)的方法，此类方法的主要缺陷在于无法有效地表现数据密集的矢量场；第二个阶段是基于流线，流面和流体的方法，此类方法的主要缺陷在于流线、流面、流体的形成预先选定的种子点关系密切，有可能忽略矢量场中重要的特征信息；第三个阶段是基于纹理映射的方法，主要有点噪声纹理合成的方法和线卷积积分方法LIC。基于纹理映射的可视化方法尤其是基于LIC的方法，由于能够较好的克服前两类方法的主要缺陷，因而得到了快速的发展，其应用的领域十分广泛。

目前，市面上有一些针对可视化的商业或开源软件，其中有一部分提供了矢量场可视化功能，但是大多采用基于几何形状或基于颜色、光学特征的矢量场映射方法，而且针对的都是结构化矢量场。其中，AVS/Express使用流线、矢量、图形符来表示矢量；IDL使用箭头、颜色等表示矢量；OpenDX也使用颜色和箭头来表示矢量；VTK使用流线等几何形状表示矢量，因属于开源项目，目前已有人为VTK开发了基于LIC方法的矢量场可视化方法，但针对的是结构化矢量场；ParaView本身是基于VTK的，因此其映射方法与VTK相同。

非结构矢量场的最大特点是采样点位置随机，将传统的结构化的线卷积积分方法LIC直接应用于非结构矢量场的LIC处理存在的问题主要有两个方面：

(1)结构化的LIC方法会因分辨率等因素导致矢量场信息的失真。以图2所示的非结构矢量场为例，有两种方法可以将其转换为结构化矢量场，一种方法是在网格内对随机采样点上的矢量进行加权处理，获得一个位于网格中心的矢量(如图4所示)，通过对比图4b和图2b可以看出，这种加权处理由于降低了采样点的分辨率，导致矢量场信息失真。另一种方法是通过加密网格使非结构矢量场的随机采样点均落在单一网格的中心，从而在保证分辨率的前提下将非结构矢量场转化成为结构化矢量场，但由于随机采样点的位置有可能非常接近，这样的处理方法人为提高了网格密度，导致计算效率低下。此外，结构化LIC的DDA曲线(矢量线，下同)形成过程是以网格单元为基本单位的，这就意味着结构化LIC处理的最高分辨率为网格单元的大小，分辨率的限制也会导致矢量场信息的失真。

(2)结构化LIC方法的矢量幅度表现力不强。在矢量场的研究领域，人们已经习惯于用矢量线的疏密来表示矢量场幅度的大小(如图5所示)，而传统的结构化LIC方法并不支持这种幅度表现方式。

本发明经过国内外专利文献和非专利文献的检索尚未发现针对非结构矢量场的线卷积积分方法。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供一种扩展的矢量场线卷积积分LIC方法，该方法既可以面向非结构矢量场，又可面向结构矢量，增强矢量幅度的表现力，提高了预处理计算速度。

本发明的技术解决方案：扩展的矢量场线卷积积分LIC方法，其特点在于通过以下步骤实现：

(1)从矢量场中寻找下一个未被标记的采样点，如果没有找到，则对矢量场的LIC值进行归一化处理，进入步骤(5)；如果找到，则进入下个步骤(2)；

(2)以上述步骤(1)找到的采样点为起点，采用基于数值计算的DDA曲线形成方法形成DDA曲线，该方法是首先从起点出发进行曲线的前向预测，然后还是从点出发进行曲线的后向预测，最后将前向预测得到的曲线和后向预测得到的曲线合并，形成最终的DDA曲线；

(3)对此轮循环形成的DDA曲线进行有色噪声序列直接赋值，计算DDA曲线上各点的LIC结果；

(4)依据此轮循环形成的DDA曲线对整个矢量场进行区域标记，依据该曲线上的矢量场采样点的幅度大小确定覆盖区域，标记该曲线区域中的所有矢量场采样数据点，并返回步骤(1)；

(5)对矢量场的LIC值进行归一化处理。

本发明与现有技术相比的有益效果是：

(1)现有的LIC方法不适用于非结构矢量场的LIC处理，而本发明给出的方法和流程即能处理各类结构化矢量场，也能处理各类非结构矢量场，这就为各类科学计算、仿真、及实测获得的矢量数据(如医学检测数据、地质勘探测量数据、气象数据、计算流体力学数据等)的可视化提供了一个更为方便和实用的手段。

(2)现有的LIC方法一方面不能很好地表现矢量场的幅度信息，另一方面会导致三维情况下处理结果因密集而互相遮挡，影响矢量场场结构的表现效果。而本发明采用的是区域标记技术，该技术在矢量场幅度方面具有更强的表现力，体现为纹理稀疏的位置幅度小，纹理密集的位置幅度大，这一方法合理地实现了LIC结果的稀疏化，因而在矢量场结构信息的方面也具有更好的表现力，有利于各领域的分析人员更好地把握矢量场的结构，从而找出最有价值的数据信息。

(3)现有的LIC方法采用随机纹理累加计算的方式进行曲线纹理映射，而本发明提出可以通过有色噪声序列直接进行曲线纹理映射，一方面适用于结构化的数据场，另一方面解决了非结构矢量场中无结构化输入纹理情况下的纹理映射问题，同时可以通过预处理提高计算速度。

附图说明

图1为本发明流程图；

图2为本发明DDA曲线前向预测示意图，其中P₀为当前DDA曲线点的位置，A₀为P₀点处的矢量，P为前向预测的下一点位置，A为P点处的矢量，P_i，i＝1，2，3，4为预测区域内的矢量场采样点位置，A_i，i＝1，2，3，4为P_i处的矢量；

图3为本发明DDA曲线后向预测示意图，其中P₀为当前DDA曲线点的位置，A₀为P₀点处的矢量，P为后向预测的下一点位置，A为P点处的矢量，P_i，i＝1，2，3，4为预测区域内的矢量场采样点位置，A_i，i＝1，2，3，4为P_i处的矢量；

图4为本发明区域标记示意图，其中，c为DDA曲线上一点矢量幅度的大小，w为该点相应的区域宽度；

图5-7为传统LIC方法与本发明方法获得的处理结果的对比图，其中

图5为放射形矢量场的LIC结果，左图为传统LIC的结果，右图为本发明方法的结果，图像分辨率均为210*210象素；

图6为圆形矢量场的LIC结果，左图为传统LIC的结果，右图为本发明方法的结果，图像分辨率均为210*210象素；

图7为雷达矢量场的LIC结果，左图为传统LIC的结果，右图为本发明方法的结果，图像分辨率均为210*210象素。

具体实施方式

本发明的具体实现方法如下：

(1)从矢量场中寻找下一个未被标记的采样点，如果没有找到，则对矢量场的LIC值进行归一化处理，结束；如果找到，则进入步骤(2)；

(2)以步骤(1)找到的采样点为起点，采用基于数值计算的DDA曲线形成方法形成DDA曲线，首先从起点出发进行曲线的前向预测，然后还是从点出发进行曲线的后向预测，最后将前向预测得到的曲线和后向预测得到的曲线合并，形成最终的DDA曲线。

如图2所示，DDA曲线前向预测方法如下：

前向预测体现为一个迭代循环，如为第一次循环，则以步骤(1)获得的采样点为当前的DDA曲线点，以该采样点处的矢量作为当前DDA曲线点处的矢量；否则以上一次前向预测得到的点作为当前DDA曲线点，以该点处的矢量作为当前DDA曲线点处的矢量。

每一次循环的处理方法如下：

为表述方便，设当前DDA曲线点的位置为P₀＝[P_0，x，P_0，y，P_0，z]，P₀点处的矢量为A₀＝[A_0，x，A_0，y，A_0，z]，以P₀为原点，以A₀为基准，根据事先设定的半径R(例如取10.0)和前向角θ(例如取180°)，可以得到一个扇形区域(三维情况下为一个圆锥体)，则矢量场在该区域中存在一组随机采样点P_i＝[P_i，x，P_i，y，P_i，z]，i＝1，2，…N，P_i点上的矢量为A_i＝[A_i，x，A_i，y，A_i，z]，P_i和A_i为给定的矢量场中的已知值。

沿前向计算两邻近点位置：事先设定一个Δk(例如取0.1)，可以计算出两点位置P_Δk＝[P_Δk，x，P_Δk，y，P_Δk，z]和P_2Δk＝[P_2Δk，x，P_2Δk，y，P_2Δk，z]：

$P_{\mathrm{\Δk},x}=P_{0,x}+\mathrm{\Δk}\frac{A_{0,x}}{\left|A_0\right|},P_{\mathrm{\Δk},y}=P_{0,y}+\mathrm{\Δk}\frac{A_{0,y}}{\left|A_0\right|},P_{\mathrm{\Δk},z}=P_{0,z}+\mathrm{\Δk}\frac{A_{0,z}}{\left|A_0\right|}$

$P_{2\mathrm{\Δk},x}=P_{0,x}+2\mathrm{\Δk}\frac{A_{0,x}}{\left|A_0\right|},P_{2\mathrm{\Δk},y}=P_{0,y}+2\mathrm{\Δk}\frac{A_{0,y}}{\left|A_0\right|},P_{2\mathrm{\Δk},z}=P_{0,z}+2\mathrm{\Δk}\frac{A_{0,z}}{\left|A_0\right|}$

根据两邻近点位置计算两邻点处的矢量：获得P_Δk和P_2Δk两点位置后，可以计算两点P_Δk和P_2Δk处的矢量A_Δk＝[A_Δk，x，A_Δk，y，A_Δk，z]，A_2Δk＝[A_2Δk，x，A_2Δk，y，A_2Δk，z]：

$A_{\mathrm{\Δk},x}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\exp \{-\sqrt{{(P_{\mathrm{\Δk},x}-P_{i,x})}^2+{(P_{\mathrm{\Δk},y}-P_{i,y})}^2+{(P_{\mathrm{\Δk},z}-P_{i,z})}^2}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp \{-\sqrt{{(P_{\mathrm{\Δk},x}-P_{j,x})}^2+{(P_{\mathrm{\Δk},y}-P_{j,y})}^2+{(P_{\mathrm{\Δk},z}-P_{j,z})}^2}}{A}_{i,x}$

$A_{\mathrm{\Δk},y}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\exp \{-\sqrt{{(P_{\mathrm{\Δk},x}-P_{i,x})}^2+{(P_{\mathrm{\Δk},y}-P_{i,y})}^2+{(P_{\mathrm{\Δk},z}-P_{i,z})}^2}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp \{-\sqrt{{(P_{\mathrm{\Δk},x}-P_{j,x})}^2+{(P_{\mathrm{\Δk},y}-P_{j,y})}^2+{(P_{\mathrm{\Δk},z}-P_{j,z})}^2}}{A}_{i,y}$

$A_{\mathrm{\Δk},z}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\exp \{-\sqrt{{(P_{\mathrm{\Δk},x}-P_{i,x})}^2+{(P_{\mathrm{\Δk},y}-P_{i,y})}^2+{(P_{\mathrm{\Δk},z}-P_{i,z})}^2}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp \{-\sqrt{{(P_{\mathrm{\Δk},x}-P_{j,x})}^2+{(P_{\mathrm{\Δk},y}-P_{j,y})}^2+{(P_{\mathrm{\Δk},z}-P_{j,z})}^2}}{A}_{i,z}$

以及

$A_{2\mathrm{\Δk},x}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\exp \{-\sqrt{{(P_{2\mathrm{\Δk},x}-P_{i,x})}^2+{(P_{2\mathrm{\Δk},y}-P_{i,y})}^2+{(P_{2\mathrm{\Δk},z}-P_{i,z})}^2}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp \{-\sqrt{{(P_{2\mathrm{\Δk},x}-P_{j,x})}^2+{(P_{2\mathrm{\Δk},y}-P_{j,y})}^2+{(P_{2\mathrm{\Δk},z}-P_{j,z})}^2}}{A}_{i,x}$

$A_{2\mathrm{\Δk},y}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\exp \{-\sqrt{{(P_{2\mathrm{\Δk},x}-P_{i,x})}^2+{(P_{2\mathrm{\Δk},y}-P_{i,y})}^2+{(P_{2\mathrm{\Δk},z}-P_{i,z})}^2}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp \{-\sqrt{{(P_{2\mathrm{\Δk},x}-P_{j,x})}^2+{(P_{2\mathrm{\Δk},y}-P_{j,y})}^2+{(P_{2\mathrm{\Δk},z}-P_{j,z})}^2}}{A}_{i,y}$

$A_{2\mathrm{\Δk},z}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\exp \{-\sqrt{{(P_{2\mathrm{\Δk},x}-P_{i,x})}^2+{(P_{2\mathrm{\Δk},y}-P_{i,y})}^2+{(P_{2\mathrm{\Δk},z}-P_{i,z})}^2}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp \{-\sqrt{{(P_{2\mathrm{\Δk},x}-P_{j,x})}^2+{(P_{2\mathrm{\Δk},y}-P_{j,y})}^2+{(P_{2\mathrm{\Δk},z}-P_{j,z})}^2}}{A}_{i,z}$

根据两邻点处的矢量计算误差的一阶导数：获得P_Δk和P_2Δk的位置及相应的矢量后，可以计算出误差的一阶导数

$e_1^{\′}=\frac{\sqrt{{(A_{\mathrm{\Δk},x}-A_{0,x})}^2+{(A_{\mathrm{\Δk},y}-A_{0,y})}^2+{(A_{\mathrm{\Δk},z}-A_{0,z})}^2}}{\mathrm{\Δk}}$

$e_2^{\′}=\frac{\sqrt{{(A_{\mathrm{\Δk},x}-A_{0,x})}^2+{(A_{\mathrm{\Δk},y}-A_{0,y})}^2+{(A_{\mathrm{\Δk},z}-A_{0,z})}^2}}{\mathrm{\Δk}}$

根据误差的一阶导数计算误差的二阶导数：获得误差的一阶导数后，可以计算误差的二阶导数

$e^{\′\′}=\frac{e_2^{\′}-e_1^{\′}}{\mathrm{\Δk}}$

根据误差的一阶及二阶导数计算预测步长：获得误差的一阶及二阶导数后，可以计算出预测步长 $k=\frac{{-e}_1^{\′}+\sqrt{{\left(e_1^{\′}\right)}^2-{2e}^{\′\′}\left|A_0\right|\mathrm{\μ}}}{e^{\′\′}},$ 其中μ为事先设定的允许误差比例，例如取0.1。将k与最大允许步长v(v为事先设定的值，例如取1.0)加以比较，如果k大于v则令k＝v，如果k小于等于v，则k不变。

根据预测步长计算此次循环的下一点预测位置：获得此轮循环的预测步长k后，可以计算出前向预测的下一点位置P＝[P_x，P_y，P_z]：

$P_x=P_{0,x}+k\frac{A_{0,x}}{\left|A_0\right|},P_y=P_{0,y}+k\frac{A_{0,y}}{\left|A_0\right|},P_z=P_{0,z}+k\frac{A_{0,z}}{\left|A_0\right|}$

根据下一点预测位置计算该点处的矢量：获得P点位置后，可以计算P点处的矢量A＝[A_x，A_y，A_z]：

$A_x=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\exp \{-\sqrt{{(P_x-P_{i,x})}^2+{(P_y-P_{i,y})}^2+{(P_z-P_{i,z})}^2}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp \{-\sqrt{{(P_x-P_{j,x})}^2+{(P_y-P_{j,y})}^2+{(P_z-P_{j,z})}^2}}{A}_{i,x}$

$A_y=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\exp \{-\sqrt{{(P_x-P_{i,x})}^2+{(P_y-P_{i,y})}^2+{(P_z-P_{i,z})}^2}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp \{-\sqrt{{(P_x-P_{j,x})}^2+{(P_y-P_{j,y})}^2+{(P_z-P_{j,z})}^2}}{A}_{i,y}$

$A_z=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\exp \{-\sqrt{{(P_x-P_{i,x})}^2+{(P_y-P_{i,y})}^2+{(P_z-P_{i,z})}^2}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp \{-\sqrt{{(P_x-P_{j,x})}^2+{(P_y-P_{j,y})}^2+{(P_z-P_{j,z})}^2}}{A}_{i,z}$

至此完成了前向预测一次循环过程。

不断以预测得到的点P及该点矢量A作为新的当前DDA曲线点及相应的矢量，并按上述步骤迭代循环，直到P超出矢量场的空间范围或矢量A的模为0，则结束前向预测，此时可以得到一串点构成的前向预测曲线及各点相应的矢量。

如图3所示，DDA曲线后向预测方法如下：

与前向预测类似，后向预测也体现为一个迭代循环，如为第一次循环，则以步骤(1)获得的采样点为当前的DDA曲线点，以该采样点处的矢量作为当前DDA曲线点处的矢量；否则以上一次后向预测得到的点作为当前DDA曲线点，以该点处的矢量作为当前DDA曲线点处的矢量。

每一次循环的处理方法如下：

为表述方便，设当前DDA曲线点的位置为P₀＝[P_0，x，P_0，y，P_0，z]，P₀点处的矢量为A₀＝[A_0，x，A_0，y，A_0，z]，以P₀为原点，以-A₀为基准，根据事先设定的半径R(例如取10.0)和后向角θ(例如取180°)，可以得到一个扇形区域(三维情况下为一个圆锥体)，则矢量场在该区域中存在一组随机采样点P_i＝[P_i，x，P_i，y，P_i，z]，i＝1，2，…N，P_i点上的矢量为A_i＝[A_i，x，A_i，y，A_i，z]，P_i和A_i为给定的矢量场中的已知值。

沿后向计算两邻近点位置：事先设定一个Δk(例如取0.1)，可以计算出两点位置P_Δk＝[P_Δk，x，P_Δk，y，P_Δkz]和P_2Δk＝[P_2Δk，x，P_2Δk，y，P_2Δk，z]：

$P_{\mathrm{\Δk},x}=P_{0,x}-\mathrm{\Δk}\frac{A_{0,x}}{\left|A_0\right|},P_{\mathrm{\Δk},y}=P_{0,y}-\mathrm{\Δk}\frac{A_{0,y}}{\left|A_0\right|},P_{\mathrm{\Δk},z}=P_{0,z}-\mathrm{\Δk}\frac{A_{0,z}}{\left|A_0\right|}$

$P_{2\mathrm{\Δk},x}=P_{0,x}-2\mathrm{\Δk}\frac{A_{0,x}}{\left|A_0\right|},P_{2\mathrm{\Δk},y}=P_{0,y}-2\mathrm{\Δk}\frac{A_{0,y}}{\left|A_0\right|},P_{2\mathrm{\Δk},z}=P_{0,z}-2\mathrm{\Δk}\frac{A_{0,z}}{\left|A_0\right|}$

根据两邻近点位置计算两邻点处的矢量：获得P_Δk和P_2Δk两点位置后，可以计算两点P_Δk和P_2Δk处的矢量A_Δk＝[A_Δk，x，A_Δk，y，A_Δk，z]，A_2Δk＝[A_2Δk，x，A_2Δk，y，A_2Δk，z]：

$A_{\mathrm{\Δk},x}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\exp \{-\sqrt{{(P_{\mathrm{\Δk},x}-P_{i,x})}^2+{(P_{\mathrm{\Δk},y}-P_{i,y})}^2+{(P_{\mathrm{\Δk},z}-P_{i,z})}^2}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp \{-\sqrt{{(P_{\mathrm{\Δk},x}-P_{j,x})}^2+{(P_{\mathrm{\Δk},y}-P_{j,y})}^2+{(P_{\mathrm{\Δk},z}-P_{j,z})}^2}}{A}_{i,x}$

$A_{\mathrm{\Δk},y}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\exp \{-\sqrt{{(P_{\mathrm{\Δk},x}-P_{i,x})}^2+{(P_{\mathrm{\Δk},y}-P_{i,y})}^2+{(P_{\mathrm{\Δk},z}-P_{i,z})}^2}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp \{-\sqrt{{(P_{\mathrm{\Δk},x}-P_{j,x})}^2+{(P_{\mathrm{\Δk},y}-P_{j,y})}^2+{(P_{\mathrm{\Δk},z}-P_{j,z})}^2}}{A}_{i,y}$

$A_{\mathrm{\Δk},z}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\exp \{-\sqrt{{(P_{\mathrm{\Δk},x}-P_{i,x})}^2+{(P_{\mathrm{\Δk},y}-P_{i,y})}^2+{(P_{\mathrm{\Δk},z}-P_{i,z})}^2}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp \{-\sqrt{{(P_{\mathrm{\Δk},x}-P_{j,x})}^2+{(P_{\mathrm{\Δk},y}-P_{j,y})}^2+{(P_{\mathrm{\Δk},z}-P_{j,z})}^2}}{A}_{i,z}$

以及

$A_{2\mathrm{\Δk},x}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\exp \{-\sqrt{{(P_{2\mathrm{\Δk},x}-P_{i,x})}^2+{(P_{2\mathrm{\Δk},y}-P_{i,y})}^2+{(P_{2\mathrm{\Δk},z}-P_{i,z})}^2}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp \{-\sqrt{{(P_{2\mathrm{\Δk},x}-P_{j,x})}^2+{(P_{2\mathrm{\Δk},y}-P_{j,y})}^2+{(P_{2\mathrm{\Δk},z}-P_{j,z})}^2}}{A}_{i,x}$

$A_{2\mathrm{\Δk},y}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\exp \{-\sqrt{{(P_{2\mathrm{\Δk},x}-P_{i,x})}^2+{(P_{2\mathrm{\Δk},y}-P_{i,y})}^2+{(P_{2\mathrm{\Δk},z}-P_{i,z})}^2}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp \{-\sqrt{{(P_{2\mathrm{\Δk},x}-P_{j,x})}^2+{(P_{2\mathrm{\Δk},y}-P_{j,y})}^2+{(P_{2\mathrm{\Δk},z}-P_{j,z})}^2}}{A}_{i,y}$

$A_{2\mathrm{\Δk},z}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\exp \{-\sqrt{{(P_{2\mathrm{\Δk},x}-P_{i,x})}^2+{(P_{2\mathrm{\Δk},y}-P_{i,y})}^2+{(P_{2\mathrm{\Δk},z}-P_{i,z})}^2}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp \{-\sqrt{{(P_{2\mathrm{\Δk},x}-P_{j,x})}^2+{(P_{2\mathrm{\Δk},y}-P_{j,y})}^2+{(P_{2\mathrm{\Δk},z}-P_{j,z})}^2}}{A}_{i,z}$

根据两邻点处的矢量计算误差的一阶导数：获得P_Δk和P_2Δk的位置及相应的矢量后，可以计算出误差的一阶导数

$e_1^{\′}=\frac{\sqrt{{(A_{\mathrm{\Δk},x}-A_{0,x})}^2+{(A_{\mathrm{\Δk},y}-A_{0,y})}^2+{(A_{\mathrm{\Δk},z}-A_{0,z})}^2}}{-\mathrm{\Δk}}$

$e_2^{\′}=\frac{\sqrt{{(A_{\mathrm{\Δk},x}-A_{0,x})}^2+{(A_{\mathrm{\Δk},y}-A_{0,y})}^2+{(A_{\mathrm{\Δk},z}-A_{0,z})}^2}}{-\mathrm{\Δk}}$

根据误差的一阶导数计算误差的二阶导数：获得误差的一阶导数后，可以计算误差的二阶导数

$e^{\′\′}=\frac{e_2^{\′}-e_1^{\′}}{-\mathrm{\Δk}}$

根据误差的一阶及二阶导数计算预测步长：获得误差的一阶及二阶导数后，可以计算出预测步长 $k=\frac{{-e}_1^{\′}+\sqrt{{\left(e_1^{\′}\right)}^2-{2e}^{\′\′}\left|A_0\right|\mathrm{\μ}}}{e^{\′\′}},$ 其中μ为事先设定的允许误差比例，例如取0.1。将k与最大允许步长v(v为事先设定的值，例如取1.0)加以比较，如果k大于v则令k＝v，如果k小于等于v，则k不变

根据预测步长计算此次循环的下一点预测位置：获得此轮循环的预测步长k后，可以计算出前向预测的下一点位置P＝[P_x，P_y，P_z]：

$P_x=P_{0,x}-k\frac{A_{0,x}}{\left|A_0\right|},P_y=P_{0,y}-k\frac{A_{0,y}}{\left|A_0\right|},P_z=P_{0,z}-k\frac{A_{0,z}}{\left|A_0\right|}$

根据下一点预测位置计算该点处的矢量：获得P点位置后，可以计算P点处的矢量A＝[A_x，A_y，A_z]：

$A_x=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\exp \{-\sqrt{{(P_x-P_{i,x})}^2+{(P_y-P_{i,y})}^2+{(P_z-P_{i,z})}^2}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp \{-\sqrt{{(P_x-P_{j,x})}^2+{(P_y-P_{j,y})}^2+{(P_z-P_{j,z})}^2}}{A}_{i,x}$

$A_y=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\exp \{-\sqrt{{(P_x-P_{i,x})}^2+{(P_y-P_{i,y})}^2+{(P_z-P_{i,z})}^2}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp \{-\sqrt{{(P_x-P_{j,x})}^2+{(P_y-P_{j,y})}^2+{(P_z-P_{j,z})}^2}}{A}_{i,y}$

$A_z=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\exp \{-\sqrt{{(P_x-P_{i,x})}^2+{(P_y-P_{i,y})}^2+{(P_z-P_{i,z})}^2}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp \{-\sqrt{{(P_x-P_{j,x})}^2+{(P_y-P_{j,y})}^2+{(P_z-P_{j,z})}^2}}{A}_{i,z}$

至此，完成了后向预测一次循环过程，不断以预测得到的点P及该点矢量A作为新的当前DDA曲线点及相应的矢量，并按上述步骤迭代循环，直到P超出矢量场的空间范围或矢量A的模为0，则结束后向预测，此时可以得到一串点构成的后向预测曲线及各点相应的矢量。

将前向预测得到的曲线和后向预测得到的曲线合并的方法如下：从后向预测的最后一点开始，按逆序连接后向预测得到的曲线上各点，再按顺序连接前向预测得到的曲线上各点，即可形成一条完整的DDA曲线。

(3)对此轮循环形成的DDA曲线进行有色噪声序列直接赋值，计算DDA曲线上各点的LIC结果。

为表述方便，设此轮循环生成的DDA曲线有M个点，记为C_i＝[C_i，x，C_i，y，C_i，z]，i＝1，2，…，M，并设到此轮循环开始前已经获得的结果点集为R_j＝[R_j，x，R_j，y，R_j，z]，j＝1，2，…，N，各结果点相应的LIC值为L_j，j＝1，2，…，N，各结果点相应的相关点数为n_j，j＝1，2，…，N。

具体的计算步骤如下：

生成有色噪声序列并为DDA曲线点赋LIC值：具体方法为，先生成一个随机白噪声序列S_k，k＝1，2，…，M+K，K为附加长度(例如可以取10)，再由白噪声序列构造有色噪声序列Sc_i，i＝1，2，…，M， ${\mathrm{Sc}}_i=\underset{k=i}{\overset{i+K}{\mathrm{\Σ}}}{S}_k.$ Sc_i即为DDA曲线点C_i对应的LIC值。

判断DDA曲线上各点与已获得的结果点集中点的关系，并做相应处理：将DDA曲线各点C_i，i＝1，2，…，M分别与以往循环获得的结果集中的点R_j＝[R_j，x，R_j，y，R_j，z]，j＝1，2，…，N加以比较，如果C_i与所有R_j，j＝1，2，…，N的距离 $\sqrt{{(C_{i,x}-R_{j,x})}^2+{(C_{i,y}-R_{j,y})}^2+{(C_{i,z}-R_{j,z})}^2}>G$ (G为给定门限值，例如取1.0)，表明C_i与以往循环获得的结果集中的任何一点均不相同同，则将C_i做为新的结果加入到结果集合中，并以C_i对应的有色噪声序列中的值Sc_i作为新结果的LIC值，同时将该新结果的相关次数置1；如果C_i与以往循环获得的结果集中某R_j的距离 $\sqrt{{(C_{i,x}-R_{j,x})}^2+{(C_{i,y}-R_{j,y})}^2+{(C_{i,z}-R_{j,z})}^2}\≤G$ (G为给定门限值，例如取1.0)，表明C_i与R_j相同，则不生成新的结果，而将C_i对应的有色噪声序列中的值Sc_i加到R_j的LIC值上，即新的L_j等于原L_j+Sc_i，并将R_j的相关次数加1，即新的n_j等于原n_j+1。

(4)依据此轮循环形成的DDA曲线对整个矢量场进行区域标记，依据该曲线上的矢量场采样点的幅度大小确定覆盖区域，标记该曲线区域中的所有矢量场采样数据点。具体的计算步骤如下：

如图4所示，以此轮循环获得的DDA曲线为中心线，根据曲线上各点的矢量幅度c的大小决定相应的区域宽度W，c越大则W越小(可令W＝(c_max-c)/K，K为某预先设定的常数，c_max为矢量场中采样点矢量幅度的最大值)，从而形成一个标记区域(在三维空间，这样的区域是一个旋转体围成的柱空间)，对落入该标记区域的矢量场中的随机采样点进行标记，被标记的点不在作为以后各轮循环的DDA曲线形成的起点。

(5)对矢量场的LIC值进行归一化处理。具体的计算步骤如下：

对最终形成的LIC结果集{R_j，j＝1，2，…，N}，其中每一个R_j均对应有一个LIC值L_j和一个相应的相关点数n_j，归一化即将R_j的LIC值L_j除以相关点数n_j作为R_j的最终LIC值。

矢量场数据预处理得到的最终结果为结果点集的位置序列{R_j，j＝1，2，…，N}，其中R_j＝[R_j，x，R_j，y，R_j，z]，及相应的LIC值序列{L_j，j＝1，2，…，N}

如图5-7为传统LIC方法与本发明方法获得的处理结果的对比，左图为传统LIC的结果，右图为本发明方法的结果，图像分辨率均为210*210象素、以三种典型的矢量场为例进行对比实验，在每组实验中，随机选择40000个采样点的位置，并给出采样点上位置上的矢量，从而构造出相应的非结构矢量场。再以非结构化的矢量场为基础，采用矢量插值的方法构造出结构化的矢量场。分别运用传统的LIC方法和本文的LIC方法进行LIC处理。从两种方法的结果对比可以看出，传统的LIC方法因分辨率的原因不能准确地给出矢量线的走向，而本发明的方法由于采用了基于数值的DDA曲线形成方法，因而能够更准确地描绘出矢量线走向的细节信息。此外，和传统的结构化LIC方法相比，本发明方法采用的区域标记方法在矢量场幅度方面具有更强的表现力，体现为纹理稀疏的位置幅度小，纹理密集的位置幅度大，这一方法合理地实现了LIC结果的稀疏化，因而在矢量场结构信息方面也具有更好的表现力。

Claims (6)

1、扩展的矢量场线卷积积分方法，其特征在于通过以下步骤实现：

(1)从矢量场中寻找下一个未被标记的采样点，如果没有找到，则对矢量场的LIC值进行归一化处理，结束；如果找到，则进入步骤(2)；

(2)以步骤(1)找到的采样点为起点，采用基于数值计算的DDA曲线形成方法形成DDA曲线，即首先从起点出发进行曲线的前向预测，然后还是从起点出发进行曲线的后向预测，最后起将前向预测得到的曲线和后向预测得到的曲线合并，形成最终的DDA曲线；

(3)对此轮循环形成的DDA曲线进行有色噪声序列直接赋值，计算DDA曲线上各点的LIC结果；

(4)依据此轮循环形成的DDA曲线对整个矢量场进行区域标记，依据该曲线上的矢量场采样点的幅度大小确定覆盖区域，标记该曲线区域中的所有矢量场采样数据点；

(5)对矢量场的LIC值进行归一化处理。

2、根据权利要求1所述的扩展的矢量场线卷积积分方法，其特征在于：所述步骤(2)中的DDA曲线前向预测方法如下：

前向预测体现为一个迭代循环，如为第一次循环，则以步骤(1)获得的采样点为当前的DDA曲线点，以该采样点处的矢量作为当前DDA曲线点处的矢量；否则以上一次前向预测得到的点作为当前DDA曲线点，以该点处的矢量作为当前DDA曲线点处的矢量，

每一次循环的处理方法如下：

(1)沿前向计算两邻近点位置

(2)根据两邻近点位置计算两邻点处的矢量

(3)根据两邻点处的矢量计算误差的一阶导数

(4)根据误差的一阶导数计算误差的二阶导数

(5)根据误差的一阶及二阶导数计算预测步长

(6)根据预测步长计算此次循环的下一点预测位置

(7)根据下一点预测位置计算该点处的矢量

(8)将此次循环获得的预测位置及矢量作为新的循环的当前DDA曲线点及矢量迭代循环，直至预测位置超出矢量场空间范围或预测位置处的矢量为0矢量，终止前向预测，从而获得由一串点构成的前向曲线及相应点处的矢量。

3、根据权利要求1所述的扩展的矢量场线卷积积分方法，其特征在于：所述步骤(2)中的DDA曲线后向预测方法如下：

后向预测也体现为一个迭代循环，如为第一次循环，则以步骤(1)获得的采样点为当前的DDA曲线点，以该采样点处的矢量作为当前DDA曲线点处的矢量；否则以上一次后向预测得到的点作为当前DDA曲线点，以该点处的矢量作为当前DDA曲线点处的矢量，

每一次循环的处理方法如下：

(1)沿后向计算两邻近点位置

(2)根据两邻近点位置计算两邻点处的矢量

(3)根据两邻点处的矢量计算误差的一阶导数

(4)根据误差的一阶导数计算误差的二阶导数

(5)根据误差的一阶及二阶导数计算预测步长

(6)根据预测步长计算此次循环的下一点预测位置

(7)根据下一点预测位置计算该点处的矢量

(8)将此次循环获得的预测位置及矢量作为新的循环的当前DDA曲线点及矢量迭代循环，直至预测位置超出矢量场空间范围或预测位置处的矢量为0矢量，终止后向预测，从而获得由一串点构成的后向曲线及相应点处的矢量。

4、根据权利要求1所述的扩展的矢量场线卷积积分方法，其特征在于：所述的步骤(2)中将前向预测得到的曲线和后向预测得到的曲线合并的方法如下：从后向预测的最后一点开始，按逆序连接后向预测得到的曲线上各点，再按顺序连接前向预测得到的曲线上各点，即可形成一条完整的DDA曲线。

5、根据权利要求1所述的扩展的矢量场线卷积积分方法，其特征在于：所述步骤(3)的DDA曲线进行有色噪声序列直接赋值，计算DDA曲线上各点的LIC结果的方法如下：

(1)生成有色噪声序列并为DDA曲线点赋LIC值；

(2)判断DDA曲线上各点与已获得的结果点集中点的关系，并做相应处理，如果DDA曲线上某点不与结果点集中的任何点相同，则将该DDA曲线点作为新的结果点加入结果点集；如果DDA曲线上某点与结果点集中的一点相同，则将DDA曲线上该点的LIC值叠加到结果点集中与之相同的点的LIC值上，并对该结果点集中点的相关点数做加1处理。

6、根据权利要求1所述的扩展的矢量场线卷积积分方法，其特征在于：所述步骤(4)的区域标记方法如下：以DDA曲线为中心线，根据曲线上各点的矢量幅度c的大小决定相应的区域宽度w，c越大则w越小，从而形成一个标记区域，对落入该标记区域的所有矢量场中的随机采样点进行标记。

Images