一种基于流体运动矢量场的并行分析方法
【专利摘要】本发明提供了一种基于流体运动矢量场的并行分析方法，具体步骤包括：将图像序列的各图像分解为同等大小的若干块小区域；构建灰度分布矩阵；并行化获取图像序列的初始运动矢量场；通过新的平滑性约束项对初始运动矢量场进行并行化处理去噪；通过新的能量约束函数对去噪后的运动矢量场进行并行化处理优化；获取包含流体运动图像序列的主方向运动矢量。本发明所述的流体运动矢量分析方法结合能量优化技术并利用图形处理器GPU的单指令多线程SIMT特性，使流体分析从串行转换为并行，实现对图像的并行化处理，提高了图像处理的性能和效率，大大提高了图像处理的精度，也降低系统处理图像序列的时间。
【专利说明】
一种基于流体运动矢量场的并行分析方法
技术领域
[0001] 本发明涉及计算机图形学领域，尤其涉及一种基于流体运动矢量场的并行分析方 法。
【背景技术】
[0002] 在进行图像运动分析过程中，一般存在两种分析方法:光流分析方法和运动矢量 分析方法。光流分析方法主要表示几何变化情况，利用图像序列之间的像素强度值的运动 变化，来确定图像中每一个像素的运动变化情况，构建出图像中每个像素点的运动矢量;运 动矢量分析方法主要针对相邻两帧图像之间统一对应点的变换情况，来获得每个像素点的 运动矢量，从而构建一个运行矢量场。
[0003] 传统的图像序列处理都是利用运动矢量分析方法或者光流分析方法进行形态矩 阵初始化表示，构建出形矩阵中的基本数据，将形态矩阵表示的图像数据进行去噪处理，得 到了比较好的数据表示，最后利用能量约束算法进行数据的深入处理，获得到实验的结果。
[0004] 流体运动的变化时而缓慢，时而剧烈，造就了流体运动相邻之间图像灰度值千变 万化。在流体运动变化相对缓慢的状态下，形态矩阵能够精确地表示图像的变换情况;但在 流体运动变化相对剧烈的状态下，包含流体运动的图像序列的相邻两帧图像序列之间变化 剧烈，其灰度值变化比较频繁，形态矩阵难以精确地表示图像的变换情况。
[0005] 光流分析方法处理包含流体运动的图像序列时，难以处理相邻两帧图像之间存在 剧烈变化的灰度值的图像序列，其处理结果存在巨大误差;而传统的运动矢量分析方法对 包含流体运动的图像序列进行分析处理的过程中也存在以下几个方面的不足:其一，形态 矩阵在表示图像序列的数据时，无法对流体运动中频繁变化的灰度值进行表示;其二，传统 的去噪算法不能处理包含流体运动的图像序列中剧烈变化的灰度，如连续性方程和灰度连 续性约束算法在处理剧烈变化的灰度值图像时，会出现很大的误差;其三，传统的能量约束 函数在处理较大灰度变化的图像序列的时候，会因为参数的取值不同，而导致在处理图像 的时候，出现较大的偏差。
【发明内容】
[0006] 为了克服现有的图像序列处理方法在处理流体运动图像序列过程中存在的处理 效率低和精度不高的问题，本发明提出了一种基于流体运动矢量场的并行分析方法，利用 已有的形态矩阵和图像序列灰度值构建灰度分布矩阵(Gray Distribution Matrix，简称 GDM)，提高表示图像基本信息的准确性;改进去噪算法，对灰度分布矩阵进行优化操作；改 进能量约束函数，实现新的能量约束优化;利用图形处理器(Graphic Processing Unit，简 称GPU)的单指令多线程(Single Instruction Multiple Threads，简称SIMT)的特性，将对 包含流体运动图像的分析从串行处理转换为并行处理，提高了图像处理的性能、效率和精 度，也降低系统处理图像序列的时间。
[0007] 本发明提出了一种基于流体运动矢量场的并行分析方法，包括步骤如下：
[0008]将图像序列的各图像分解为同等大小的aXa块小区域；
[0009]构建灰度分布矩阵；
[0010]并行化获取图像序列的初始运动矢量场；
[0011]通过新的平滑性约束项对初始运动矢量场进行并行化处理去噪；
[0012]通过新的能量约束函数对去噪后的运动矢量场进行并行化处理优化；
[0013 ]获取包含流体运动图像序列的主方向运动矢量。
[0014] 其中，所述构建灰度分布矩阵的具体步骤如下：图像分解为aXa块，构建一个对应 大小为a X a的区域下面的灰度分布矩阵，故灰度分布矩阵描述如下：
[0015] GDM(x,y,t) = (S°(x,y,t),S1(x,y,t)^. . .SL(x,y,t)) ?…V
[0016] 其中，灰度级信息矩阵S'(.、:，.w)= : ···：; p. "， p， 、^??-Ι J
[0017] 式中，L表示为灰度级;每个灰度级表示某个区间段的灰度值，灰度值的范围就是0 ~255。在一种优选实施例中，灰度级为0表不[0,10]区间段的灰度值。ri表不的是图像的第 (c+1)块区域，c = 0，l, . . .，aXa-l;<表示的是r。区域内所有的灰度级L的值相加之和。
[0018] 在一种优选实施例中，令a = 4，即将图像分解为4 X 4块，构建一个大小为4 X 4的区 域下面的灰度分布矩阵，L的取值范围为0~25,令灰度级L = 25故灰度分布矩阵描述如下：
[0019] GDM(x,y,t) = {S°(x,y,t),S1(x,y,t)^ · . .S25(x,y,t)}；
[K …
[0020] 其中，灰度级信息矩阵S25(xj,i)= ? ; U,…々J
[0021] 式中，r。表示的是图像的第(c+1)块区域，c = 0，l，. . .，15;及表示的是r。区域内所 有的灰度级的值相加之和。
[0022] 其中，所述构建初始运动矢量场的具体步骤包括:将图像分解为a X a块，每一块区 域大小相同，利用图形处理器GPU的单指令多线程SMT特性对各块区域进行并行化处理， GPU中每一个线程处理一个小区域的搜索窗口；每一个线程在其处理的小区域内进行循环 遍历操作，针对像素点(x，y)计算原像素和当前搜索窗口里面的每个像素之间的距离值，取 距离最小值作为对应的像素信息，可得到运动场中点(x，y)对应的运动矢量，循环迭代获得 该小区域内每一个像素点的运动矢量，并行化处理获得各个小区域的运动矢量场，由此得 到包含流体运动的图像序列对应的运动矢量场。
[0023] 在一种优选实施例中，所述对初始运动矢量场进行去噪的具体步骤包括:通过将 图像序列分成若干同等大小的模块，在每个子区域模块中，统计分析出大多数像素点的方 向：上、下、左和右四个方向；针对每一个子区域，取其中像素点最多的方向作为该子区域运 动矢量方向，移除该子区域内不在矢量方向的像素信息，从而得到该子区域去噪后的运动 矢量场，将所有去噪后的子区域的运动矢量场合成一个运动矢量场，即为包含流体运动的 图像序列去噪后的运动矢量场。
[0024] 其中，所述新的平滑约束项的获取方式包括:根据Snake模型，结合去噪后的运动 矢量场数据，得到新的平滑约束项。
[0025] Snake模型由Kass等人在1987年提出，基本思想为：以一些图像中的一些控制点为 轮廓线，通过模板自身的弹性形变，结合图像局部特征因素相匹配，最终达到调和，从而完 成图像的分割。
[0026] 其中，所述新的平滑性约束项如下： xfx y+n
[0027] U"，'，)= Σ Σ - v".): i-χ-η j- v-n
[0028] 式中，Fimage(u，v)表示大小为u*v区域的平滑约束函数，Uij和 Vij分别为运动矢量场 中位于(X，y)在X轴和y轴方向的运动矢量分量，i表示当前正在计算像素点U方向的坐标信 息，j表示当前正在计算像素点v方向的坐标信息，η表示我们整个能量函数的搜索状态空间 大小。
[0029] 在一种优选实施例中，所述对去噪后的运动矢量场进行新的能量约束优化的具体 步骤包括:通过新的能量约束函数计算出图像分割的轮廓线和图像序列中每个点位置上的 运动矢量，根据去噪后的矢量场数据信息，遍历每个子区域，通过能量约束函数迭代操作， 对能量函数进行最小化，寻找能量函数的最小值，即找到去噪后的运动矢量场中缺失位置 的原像素对应的点，也就得到了更加详细的灰度分布数据信息，即运动矢量场中缺失的位 置上像素信息，实现处理包含较大灰度变化的图像序列。
[0030] 其中，所述新的能量约束函数的获取方式包括:结合新的平滑约束条件和经过平 滑约束优化后的灰度分布矩阵，得到新的能量约束函数。
[0031] 其中，所述新的能量约束函数F(u，v)为：
[0032] F(u, v) =Fimage(u, v)+aFcon(u, v)
[0033] 其中，^。"表示的是区域以外的能量的约束与小区域内的图像信息约束对应关系， α表示的是权衡因子。
[0034]其中，对流体运动图像序列进行并行处理的方式为:基于图形处理器GPU的SMT特 性对流体运动图像序列进行并行化处理。
[0035] 在一种优选实施例中，对流体运动图像序列进行并行化处理即在能量约束优化阶 段，充分利用能量函数最小化过程是可并行化的特性，提高程序的执行效率，缩短能量约束 优化过程的处理时长。
[0036] 本发明提出了一种基于流体运动矢量场的并行分析方法，通过基于形态矩阵的新 型灰度分布矩阵表示图像序列的灰度值信息，从而更精确地表示图像序列的运动矢量信 息，同时利用GPU的SMT并行化特性进行数据处理，实现对图像的并行化处理，提高了图像 处理的性能和效率，同时也大大提高了图像处理的精度。
【附图说明】
[0037] 图1是本发明的流程图；
[0038] 图2为实施例中的图像序列不意图；
[0039] 图3为实施例中的灰度分布矩阵求解过程伪代码图；
[0040] 图4为实施例中的去噪算法实施细节伪代码图；
[0041] 图5为实施例中的能量约束实施细节伪代码图；
[0042] 图6为本发明与现有技术对本实施例处理得到的效果对比图。
【具体实施方式】
[0043] 结合以下具体实施例和附图，对本发明作进一步的详细说明。实施本发明的过程、 条件、试验方法等，除以下专门提及的内容之外，均为本领域的普遍知识和公知常识，本发 明没有特别限制内容。
[0044] 实施例
[0045] 如图1所示，为本实施例的分析过程流程示意图。
[0046] 如图2所示，为本实施例中的待分析的包含流体运动的图像序列，将图像序列的各 图像分解为同等大小的4 X 4块小区域。
[0047] 根据分解的同等大小的4 X 4块小区域，构建一个对应的4 X 4灰度分布矩阵，令L = 25,即包含26个灰度级，每个灰度级表示某个区间段的灰度值。利用opencv库中相应的函 数，获得图像的灰度值。根据图像中每个像素的RGB值可以得到图像的灰度值，灰度值的范 围就是0~255。本实施例中，26个灰度级与其对应的灰度值的区间段如下表所示：
[0050] 灰度分布矩阵描述如下：
[0051] GDM(x,y,t) = {S°(x,y,t),S1(x,y,t)^ · . .S25(x,y,t)}；
[0052] 当L = 25,灰度级信息矩阵、
[0053] 式中，r。表示的是图像的第(c+1)块区域，c = 0，l，. . .，15;&表示的是r。区域内所 有像素点的灰度级的值相加之和，循环迭代整个图像序列，得到完整的灰度分布矩阵。
[0054] 为流体运动矢量计算并行化过程中灰度分布矩阵求解的过程。通过灰度值和灰度 级之间的对应关系，在进行计算灰度分布矩阵之前，首先将图像中的灰度值信息进行转化， 将其转化为灰度级信息，从而在进行灰度分布矩阵计算的时候，识别灰度级信息。
[0055] 利用本发明提出的基于形态矩阵以及图像灰度值构建的灰度分布矩阵可以极大 地降低数据表示的不准确性和误差;表现的数据信息更加精确，相比较传统的数据表示方 法，灰度分布矩阵中的噪声数据更加少。
[0056]并行化处理图像序列，获取初始运动矢量场。利用GPU的SMT特性对各小区域进行 并行化处理，将待搜索的图像的每一块小区域看作一个搜素窗口，每一个线程处理一块小 区域的搜索窗口；获取每个小区域的运动矢量场，从而获得图像序列的初始运动矢量场。如 图3所示，为每块小区域的处理流程示意图，具体过程包括：
[0057] 在一个小区域内进行遍历操作，针对图像中的一个像素点(x，y)，在搜索窗口内寻 找与该像素点的灰度值接近的像素点(χ+i，y+j);
[0058] 利用compare⑶Μ函数计算出整个图像中像素点（x，y)和当前对应的搜索窗口中像 素点(x+i，y+j)之间的距离值；
[0059]取距离最小值作为对应的像素信息，则就可以求出运动场中点（x，y)对应的初始 运动矢量；
[0060]循环迭代获得该小区域内每一个像素点的运动矢量，通过并行化处理获得各个小 区域的初始运动矢量场；
[0061 ]由此得到包含流体运动的图像序列对应的初始运动矢量场。
[0062] 通过新的平滑约束项对初始运动矢量场进行去噪实现噪声元素信息的初步优化， 进一步缩减灰度分布矩阵中的噪声数据。
[0063] 如图4所示，为初始运动矢量场去噪流程图，具体包括步骤：
[0064] 将初始运动矢量场分解为同等大小4X4块小区域，进行并行化分析处理；
[0065] 针对每一个小区域，寻找表示该区域内的像素的运动矢量；
[0066] 图像序列的主方向表示的是图像运行中运行矢量占绝大多数的所指的方向；
[0067] 判断当前正在搜索的区域内的运动矢量场与图像序列的主方向运动矢量场是否 一致，若不一致，这表示该区域没有对应的运动矢量场;反之，该小区域的运动矢量场可进 行保留和填充。
[0068] 通过新的能量约束函数对去噪后的运动矢量场进行并行化能量约束优化，获得到 更加精确的图像序列的运动矢量场，如图5所示，能量约束优化流程图，具体步骤包括：
[0069] 将去噪优化之后的运动矢量场进行切分，切分成同等大小的4X4块小区域，针对 每一块小区域进行并行化处理；
[0070] 对于每一个小区域中的每一个位置上的点，我们通过计算寻找到当前搜索窗口中 能量最小的点，最终确定搜索点的运动矢量；
[0071] 循环迭代计算过程，可以弥补每个小区域内中去噪后空白的位置信息元素的运动 矢量。
[0072] 对同一组运动的图像序列采用现有技术和本发明提出的技术方案分别对其进行 图像运动分，取得的效果对比，如图6所示，最左边的图片为原始图片，中间两张图片是现有 的两种技术方法经过处理以后得到的运动矢量数据信息，最右边的图片是经过本发明方法 处理得到的实验数据，通过两者的对比，我们可以发现经过本发明提出的技术方案处理后 的数据，能大大提高图像处理的精度，选择其中的一种现有方法与本发明的方法的运行时 间进行比较，如下表所示：
[0074] 由表可知，本发明的方法比传统的方法提高了时间效率。
[0075] 虽然针对于大多数图像序列而言，一块小区域中的灰度值变化是微乎其微的，但 在包含有流体运动的图像序列中，图像序列之间的灰度值变换相对比较频繁。因而，利用揭 示了流体密度和流速之间的重要关系的连续性方程式
·其中V表示3D的运 动矢量场，P是流体的密度，r是某一个时间点;会影响整个对应点的搜索过程。并且，由于在 高度方向的流体变化，以及在进行数据采集的时候，摄像机的角度变换等因素，造成了图像 即是具有相同的二维运动矢量数据信息，也会存在不同的灰度值信息。所以，基于连续性方 程的约束，在计算剧烈的流体运动时，存在较大的误差。故而，需要进行改进，本发明就是在 原有的理论基础之上，进行较大的优化，得到了一种基于灰度分布矩阵数据保存处理方式。
[0076] 本发明提出了一种基于流体运动矢量场的并行分析方法，通过基于形态矩阵的新 型灰度分布矩阵表示图像序列的灰度值信息，从而更精确地表示图像序列的运动矢量信 息，同时利用GPU的SMT并行化特性进行数据处理，实现对图像的并行化处理，提高了图像 处理的性能和效率，同时也大大提高了图像处理的精度。
[0077] 以上对本发明的具体实施例进行了详细描述，但其只是作为范例，本发明并不限 制于以上描述的具体实施例。对于本领域技术人员而言，任何对本发明进行的等同修改和 替代也都在本发明的范畴之中。因此，在不脱离本发明的精神和范围下所作的均等变换和 修改，都应涵盖在本发明的范围内。
【主权项】
1. 一种基于流体运动矢量场的并行分析方法，其特征在于，包括步骤如下： 将图像序列的各图像分解为同等大小的aXa块小区域； 构建灰度分布矩阵； 并行化获取图像序列的初始运动矢量场； 通过新的平滑性约束项对初始运动矢量场进行并行化处理去噪； 通过新的能量约束函数对去噪后的运动矢量场进行并行化处理优化； 获取包含流体运动图像序列的主方向运动矢量。2. 根据权利要求1所述的一种方法，其特征在于，所述构建灰度分布矩阵的具体步骤如 下：图像分解为a X a块，构建一个对应大小为a X a的灰度分布矩阵，灰度分布矩阵描述如 下： GDM(x，y，t) = {S〇(x，y，t)，Si(x，y，t)，...SL(x,y,t)}; 其中，灰度级信息矩時式中，L表示为灰度级;每个灰度级表示某个区间段的灰度值;ri表示的是图像的第（i+ 1)块区域，i = 〇，l，. . .，aXa-l;巧表示的是ri区域内所有的灰度级L的值相加之和。3. 根据权利要求1所述的一种方法，其特征在于，所述获取初始运动矢量场的具体步骤 包括： 利用图形处理器GPU的单指令多线程SIMT特性对分解的图像序列的aXa块小区域进行 并行化处理，每块小区域作为一个捜索窗口，每一个线程处理一块小区域的捜索窗口；各个 线程在其处理的小区域内进行遍历操作，针对像素点（x，y)计算原像素和当前捜索窗口里 面的像素相差的距离值，取距离最小值作为对应的像素信息，可得到运动场中点(x，y)对应 的运动矢量，循环迭代获得该小区域内每一个像素点的运动矢量，并行化处理获得各个小 区域的运动矢量场，由此得到包含流体运动的图像序列的初始运动矢量场。4. 根据权利要求1所述的一种方法，其特征在于，所述对初始运动矢量场进行去噪的具 体步骤包括:通过将图像序列分成若干同等大小的模块，在每个子区域模块中，统计分析出 大多数像素点的方向：上、下、左和右四个方向；针对每一个子区域，取其中像素点最多的方 向作为该子区域运动矢量方向，移除该子区域内不在矢量方向的像素信息，从而得到该子 区域去噪后的运动矢量场，将所有去噪后的子区域的运动矢量场合成一个运动矢量场，即 为包含流体运动的图像序列去噪后的运动矢量场。5. 根据权利要求1所述的一种方法，其特征在于，所述新的平滑约束项的获取方式包 括:根据Snake模型，结合去噪后的运动矢量场数据，得到新的平滑约束项。6. 根据权利要求1或5所述的一种方法，其特征在于，所述新的平滑性约束项如下：式中，Fimage(U，V)表示大小为U*V区域的平滑约束函数，叫和^分别为运动矢量场中位 于(x，y)在X轴和y轴方向的运动矢量分量，i表示当前正在计算像素点U方向的坐标信息，j 表示当前正在计算像素点V方向的坐标信息，η表示我们整个能量函数的捜索状态空间大 小。7. 根据权利要求1所述的一种方法，其特征在于，所述对去噪后的运动矢量场进行新的 能量约束优化的具体步骤包括:通过新的能量约束函数计算出图像分割的轮廓线和图像序 列中每个点位置上的运动矢量，根据去噪后的矢量场数据信息，遍历每个子区域，通过能量 约束函数迭代操作，对能量函数进行最小化，寻找能量函数的最小值，即找到去噪后的运动 矢量场中缺失位置的原像素对应的点，实现处理包含较大灰度变化的图像序列。8. 根据权利要求1所述的一种方法，其特征在于，所述新的能量约束函数的获取方式包 括:结合新的平滑约束条件和经过平滑约束优化后的灰度分布矩阵，得到新的能量约束函 数。9. 根据权利要求1或8所述的一种方法，其特征在于，所述新的能量约束函数F(u，v)为： F(u,v)=Fimage(u,v)+aFcon(u,v) 其中，Fc。。表示的是区域W外的能量的约束与小区域内的图像信息约束对应关系，α表 示的是权衡因子。10. 根据权利要求1所述的一种方法，其特征在于，对包含流体运动图像序列进行并行 化处理的方式为:基于图形处理器GPU的单指令多线程SIMT特性对包含流体运动图像序列 每一个步骤都进行并行化处理。
【文档编号】G06T1/20GK106097396SQ201610490479
【公开日】2016年11月9日
【申请日】2016年6月29日
【发明人】陈铭松, 王红祥, 庄涵, 徐思远, 宋进忠
【申请人】华东师范大学