一种基于GPU加速的三维粒子场及速度场重建方法与流程

一种基于gpu加速的三维粒子场及速度场重建方法
技术领域
1.本发明涉及粒子图像测速技术领域，尤其涉及了一种基于gpu加速的三维粒子场及速度场重建方法。


背景技术：

2.粒子图像测速(particle image velocimetry，piv)技术是一种非接触无干扰的空间流场定量测量手段。piv技术利用均匀散布的微小粒径示踪粒子跟随流体的运动，通过连续“打光”和“拍照”方式记录示踪粒子的运动，基于图像处理算法从连续示踪粒子图像提取出粒子的运动信息，获得流体空间运动速度的近似。层析粒子图像测速(tomographic particle image velocimetry，tomo
‑
piv)技术将piv技术和计算机断层诊断技术(computed tomography，ct)相结合，使用多视角(至少三个相机)同时拍摄并记录被照亮的粒子场，并通过光学层析方法重建光强的三维分布，获得示踪粒子的三维空间分布；然后根据粒子场通过互相关算法重建获得流场的三维速度场分布。
3.相比传统二维piv技术，tomo
‑
piv技术能够实现瞬时的三维流场速度测量和空间流场的全场定量测量，对于具有复杂流场结构(如湍流和三维涡结构的流场)的研究能够提供很大的便利。tomo
‑
piv技术具有操作简单，测量精度高，设备简易，应用范围广及获得信息数量更多等优势；但也存在需要更高的体积照明能量、庞大的数据存储量及计算效率低下等局限性。
4.tomo
‑
piv技术中获取三维粒子场的核心步骤为基于图像的粒子场层析重建算法，计算过程中权重矩阵的计算需要12h以上，且需要占用200
‑
500gb的存储空间，耗费极大的计算成本，成为piv技术测量流场的主要技术瓶颈。其次，获得粒子分布之后的速度场重建也是tomo
‑
piv技术的核心步骤，由于需要详细对比两帧粒子场分布的每一个粒子(每一个粒子近似为立方体网格，称为体素)，而往往层析重建得到的体素数量也极其庞大，因此速度场重建的效率问题也是tomo
‑
piv技术的应用阻碍之一。
5.近年来国内外关于粒子场层析重建的研究重点集中于提升tomo
‑
piv技术的计算精度，提出了包括乘法代数重建算法(multiplicative algebraic reconstruction techniques，mart)、体自标定技术、multiplicative first guess(mfg)算法、mlos
‑
mart算法及inte
‑
mart算法等在内的众多重建算法，提高了三维粒子场重建的质量和三维速度场重建的精度。但由于这些算法运行在cpu上，效率的提升并不明显。而piv技术的速度场重建中，国内外学者已经开始使用gpu(graphics processing unit，图像处理器)加速传统二维piv的或者三维piv的互相关算法，重建速度场；但这些研究没有更为深入的优化gpu线程设计，仅仅依赖gpu函数库中的cufft模块替换了互相关算法中的快速傅里叶变换，对于多线程并发模式没有进行深入研究，无法使速度场重建模块充分发挥gpu的优势。
6.总之，粒子场重建和速度场重建的效率问题严重阻碍了tomo
‑
piv技术的大规模应用，需要对现有技术进行改进。


技术实现要素：

7.本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种有效减少计算时间开销、提升效率的基于gpu加速的三维粒子场及速度场重建方法。
8.本发明的目的可以通过以下技术方案来实现：
9.一种基于gpu加速的三维粒子场及速度场重建方法，包括以下步骤：
10.获取两组具有一定时间间隔的多相机图像数据，以及每个相机的权重矩阵查找曲线；
11.采用gpu，基于所述多相机图像数据和权重矩阵查找曲线进行三维粒子场层析重建，获得重建的三维粒子场，在所述层析重建过程中，使用与体素数量等同的线程数量并行计算所有体素的光强重建；
12.对所述重建的三维粒子场进行分块，提取每个数据块的两个诊断窗口数据，采用gpu并发进程进行三维互相关计算，获得所有窗口的平均速度，重建获得三维速度场。
13.进一步地，所述权重矩阵查找曲线通过以下步骤获得：
14.进行多相机的标定，得到相机标定矩阵；
15.基于所述相机标定矩阵，将实验图像坐标映射到空间三维真实坐标系中；
16.根据映射得到的空间三维真实坐标，将空间三维真实坐标系下的每个图像像素点视线近似为圆柱，三维真实坐标系下的每个体素近似为球体，每个体素对应图像像素点的贡献值近似为圆柱与体素相交的体积与整个体素体积之比，拟合出每个体素对应每个像素的贡献值，即权重矩阵查找曲线。
17.进一步地，所述相机标定矩阵基于控制常数和标定数据获得，所述控制常数包括粒子场层析重建参数和速度场重建参数。
18.进一步地，所述粒子场层析重建参数包括图像数据编号、图像数据类型、图像数据尺寸、图像数据存放路径、相机数量、测量区域大小和每毫米体素数量，
19.所述速度场重建参数包括诊断窗口尺寸和窗口重叠率。
20.进一步地，所述三维粒子场层析重建中，使用基于gpu并行的mart算法层析重建得到三维粒子场的空间分布。
21.进一步地，所述mart算法的计算公式如下：
[0022][0023]
式中，k为迭代次数，μ为松弛参数，e(x
j
,y
j
,z
j
)为体素灰度值，i(x
i
,y
j
)为相机图像像素值，w
i,j
为每个体素j对每个像素i灰度值的贡献权重矩阵，该贡献权重矩阵根据所述多相机图像数据和权重矩阵查找曲线获得，n
i
为体素总数。
[0024]
进一步地，基于所述诊断窗口尺寸和窗口重叠率对所述重建的三维粒子场进行分块，并通过自适应选窗技术提取每个数据块的两个诊断窗口数据。
[0025]
进一步地，所述三维互相关计算具体为：
[0026]
对所述两个诊断窗口数据，分别基于gpu的cufft快速傅里叶变换进行三维傅里叶变换，得到两个诊断窗口的频域变换结果；
[0027]
基于gpu的cufft卷积计算，得到两个诊断窗口粒子场的能谱；
[0028]
基于gpu的cufft快速傅里叶逆变换进行反向傅里叶变换，得到两个诊断窗口粒子场的互相关函数；
[0029]
求出互相关函数的峰值大小以及峰值所在的三维数组索引，计算得到诊断窗口中示踪粒子的平均速度。
[0030]
进一步地，在获得所有窗口的平均速度后，查找其中是否存在错误速度矢量，若是，则进行替换处理。
[0031]
进一步地，所述时间间隔为0.001s。
[0032]
与现有技术相比，本发明具有以下有益效果：
[0033]
1、本发明运用新型gpu硬件加速tomo
‑
piv技术中的粒子场层析重建算法和速度场重建过程的三维互相关算法，能够大幅度提升tomo
‑
piv技术的计算效率，减少使用者在使用tomo
‑
piv技术相关硬件设备时的繁琐数据处理工作，以节省大量的时间，加快实验流体力学研究速度，将对tomo
‑
piv技术在流场重建的广泛应用提供基础，也对于加速实验流体研究具有重要意义。
[0034]
2、本发明基于gpu并行的层析重建能够极大地提高粒子场重建的效率，降低层析粒子图像测速技术中绝大部分的时间开销。
[0035]
3、本发明基于多任务并发的思维，将粒子场大量诊断窗口的三维互相关计算分解为多进程并发计算，极大地降低了速度场重建的计算时间开销。
附图说明
[0036]
图1为本发明方法的流程图；
[0037]
图2为本发明三维粒子场层析重建流程图；
[0038]
图3为本发明三维速度场重建流程图；
[0039]
图4为本发明实施例中计算得到的权重矩阵查找曲线图；
[0040]
图5为本发明实施例中读取的8张图像；
[0041]
图6为本发明实施例中基于gpu加速重建获得的两个时刻的三维粒子场；
[0042]
图7为本发明实施例中基于三维粒子场提取的两个时刻诊断窗口数据；
[0043]
图8为本发明实施例中并发4进程，基于三维粒子场提取的4组两个时刻诊断窗口数据；
[0044]
图9为本发明实施例中基于gpu加速三维互相关获得的互相关函数分布；
[0045]
图10为本发明实施例中三维速度场结果。
具体实施方式
[0046]
下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。
[0047]
本发明提供一种基于gpu加速的三维粒子场及速度场重建方法，该方法采用gpu，基于多相机图像数据和权重矩阵查找曲线进行三维粒子场层析重建，获得重建的三维粒子场，在层析重建过程中，使用与体素数量等同的线程数量并行计算所有体素的光强重建；对重建的三维粒子场进行分块，提取每个数据块的两个诊断窗口数据，采用gpu并发进程进行
三维互相关计算，获得所有窗口的平均速度，重建获得三维速度场。上述方法基于gpu并行的层析重建能够极大地提高粒子场重建的效率，降低层析粒子图像测速技术中绝大部分的时间开销；同时基于多任务并发的思维，将粒子场大量诊断窗口的三维互相关计算分解为多进程并发计算，极大地降低了速度场重建的计算时间开销。
[0048]
参照图1
‑
图3，为本发明的基于gpu加速层析粒子图像测速的三维粒子场及速度场重建方法的流程图，具体包括以下步骤：
[0049]
1)从输入文件中读取计算的控制常数以及标定数据，进行相机的标定得到相机标定矩阵。
[0050]
11)读取本实施例的控制常数，包括层析重构所需的图像数据编号、图像数据类型、图像数据尺寸、图像数据存放路径、相机数量、测量区域实际尺寸和每毫米体素数量等；包括速度场重建所需的诊断窗口尺寸、窗口重叠率等参数如表1所示。
[0051]
表1控制参数
[0052][0053]
12)读取标定数据，包括将进行相机标定和计算标定矩阵a所需的全部标定数据。
[0054]
13)相机的标定，根据步骤11)和步骤12)的控制常数数据和标定数据，求解得到相机标定矩阵a中的12个未知数a
11
，a
12
，a
13
，a
14
，a
21
，a
22
，a
23
，a
24
，a
31
，a
32
，a
33
，a
34
，构建的空间三维真实坐标系和图像平面像素坐标系的标定函数关系为：
[0055][0056]
其中，λ为尺度因子，(u,v)为图像点p的像素坐标，(x,y,z)为三维真实空间体素点p的三维真实坐标；a为相机投影变换矩阵。
[0057]
2)根据步骤1)计算得到的相机标定矩阵，通过标定函数关系将实验图像坐标映射到空间三维真实坐标，建立权重矩阵查找曲线。
[0058]
21)根据步骤1)得到的相机标定矩阵a，读取两张实验图像，将实验图像通过标定函数关系将实验图像坐标映射到空间三维真实坐标系中。
[0059]
22)根据步骤21)得到的空间三维真实坐标结果，将空间三维真实坐标系下的每个图像像素点视线近似为圆柱，三维真实坐标系下的每个体素近似为球体，每个体素对应图像像素点的贡献值近似为视线圆柱与体素相交的体积与整个体素体积之比，拟合出每个体素对应每个像素的贡献值，即权重矩阵查找曲线。
[0060]
本实施例中，相机数量为4，每个相机均需要计算一条权重矩阵查找曲线，计算得到4条权重矩阵查找曲线，每张表数据长度为18000个点，4条权重矩阵查找曲线结果如图4所示。
[0061]
3)基于gpu的三维粒子场层析重建。
[0062]
如图2所示，三维粒子场层析重建具体包括：
[0063]
31)读取时间间隔为0.001s的两组多相机图像数据，包括层析重建过程所需的相
机数量、每个相机记录的图像以及每张图像每个像素的坐标及灰度值，并将其传递到gpu设备内存当中。
[0064]
本实施例中，相机数量为4，拍摄时间间隔为0.001s，两个时刻四个相机各拍摄一张图像，图像总数为8，图像大小为1000
×
2000pixel，将图像编号为1～4和5～8，1～4为第一个时刻的4张图像，5～8为第二个时刻的4张图像，读取的所有图像如图5所示。
[0065]
32)将步骤2)中得到的每个相机的权重矩阵查找曲线拷贝传递给gpu内存。
[0066]
33)基于三维空间中所有体素中心坐标的体素灰度值e(x
j
,y
j
,z
j
)与每台相机拍摄图像所有像素i(x
i
,y
i
)的投影关系，根据步骤31)传入gpu的相机图像数据，和步骤32)传入gpu的权重矩阵查找曲线计算得到每个体素j对每个相机图像的每个像素i灰度值的贡献权重矩阵w
i,j
。
[0067]
本实施例中，测量区域尺寸为50mm
×
100mm
×
100mm，每毫米的体素数量为11voxels。使用gpu多线程同时计算两个时刻所有体素j对图像编号1～8的所有图像像素i的贡献权重矩阵w
i,j
。
[0068]
34)根据步骤33)计算得到的贡献权重矩阵w
i,j
，使用基于gpu并行的mart算法层析重建得到三维粒子场的灰度值空间分布。基于gpu的层析重建计算，将每个体素j与所有相机像素的光强重建过程采用一个gpu线程来进行计算，利用gpu设备中线程数量庞大的特点，使用与体素数量等同的线程数量并行计算所有体素的光强重建过程。加速之前基于cpu的层析重建计算流程，以较高的计算效率重建得到三维粒子场。
[0069]
所述mart算法的计算公式如下：
[0070][0071]
式中，k为迭代次数，μ为松弛参数，e(x
j
,y
j
,z
j
)为体素灰度值，i(x
i
,y
j
)为相机图像像素值，w
i,j
为每个体素j对每个像素i灰度值的贡献权重矩阵，该贡献权重矩阵根据所述多相机图像数据和权重矩阵查找曲线获得，n
i
为体素总数。
[0072]
本实施例中，相机数量为4，拍摄时间间隔为0.001s，两个时刻四个相机各拍摄一张图像，图像总数为8，图像大小为1000
×
2000pixel，将图像编号为1～4和5～8，1～4为第一个时刻的4张图像，5～8为第二个时刻的4张图像；测量区域尺寸为50mm
×
100mm
×
100mm，每毫米的体素数量为11voxels。使用gpu同时多线程同时计算两个时刻所有粒子(体素)的灰度值，两个时刻的三维粒子场如图6所示。
[0073]
4)基于gpu的三维速度场重建。
[0074]
如图2所示，三维速度场重建具体包括：
[0075]
41)根据步骤3)得到的两个时刻的三维粒子场数据，基于步骤1)得到的诊断窗口大小，对粒子场数据进行分块，通过自适应选窗技术提取每个数据块中用于三维互相关计算的两个诊断窗口数据，将其存入gpu内存。
[0076]
本实施例中，诊断窗口大小32voxel
×
32voxel
×
32voxel，窗口重叠率为50％，测量区域三个维度的诊断窗口总数为32
×
64
×
32；提取第一个时刻，体素(粒子)坐标中心为(17voxel,17voxel,17voxel)的诊断窗口数据a，提取第二个时刻，体素坐标中心为(17voxel,17voxel,17voxel)的诊断窗口数据b，如图7所示，将其拷贝存入gpu内存。
[0077]
42)根据计算设备gpu内存的大小，自适应设计并发进程数量，每条进程均可以完成步骤41)的诊断窗口数据选择，并独立处理一个三维互相关计算模块，最大并发进程同时处理所有的诊断窗口三维互相关计算。
[0078]
本实施例中，根据计算设备gpu内存的大小，设置并发进程数量为4；4条进程分别提取第一个时刻的诊断窗口数据a1，a2，a3，a4；和第二个时刻的诊断窗口数据b1，b2，b3，b4，如图8所示，将其拷贝存入gpu内存。
[0079]
43)对步骤42)的每条进程，根据步骤41)提取的两个诊断窗口数据，分别基于gpu的cufft快速傅里叶变换模块进行三维傅里叶变换，得到两个诊断窗口的频域变换结果。
[0080]
44)对步骤43)计算获得的两个诊断窗口的频域变换结果进行基于gpu的cufft卷积计算，得到两个诊断窗口粒子场的能谱；然后基于gpu的cufft快速傅里叶逆变换模块进行一次反向傅里叶变换，获得4条进程两个诊断窗口粒子场的互相关函数如图9所示。
[0081]
45)根据步骤44)得到的两个诊断窗口的互相关函数，求出互相关函数的峰值大小，以及峰值所在的三维数组索引，计算得到诊断窗口中示踪粒子的平均速度。
[0082]
本实施例中，根据图9的四条进程所示的互相关函数分布可得峰值坐标、峰值及三个维度的平均速度如表2所示。
[0083]
表2四条进程计算所得互相关函数峰值坐标、峰值及三个维度的平均速度
[0084][0085]
45)重复步骤42)，步骤43)，步骤44)，直到完成所有诊断窗口的平均速度计算。
[0086]
本实施例中，诊断窗口大小32voxel
×
32voxel
×
32voxel，窗口重叠率为50％，测量区域三个维度的诊断窗口总数为32
×
64
×
32；一共执行32
×
64
×
32＝65536次互相关循环。
[0087]
46)对步骤45)得到的所有诊断窗口中示踪粒子的平均速度进行分析，替换其中错误的速度矢量。
[0088]
5)对步骤2)、步骤3)、步骤4)计算得到的三维粒子场和三维速度场结果进行输出保存，并结束计算。
[0089]
本实施例的三维速度场结果如图10所示。
[0090]
本实施例使用cpu和不同gpu的性能提升如表3所示。
[0091]
表3性能分析
[0092][0093][0094]
以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的
技术方案，皆应在由权利要求书所确定的保护范围内。