基于新型相关函数约束的cuda架构并行优化立体变形测量方法
技术领域
本发明属于基于视觉技术的立体变形非接触测量
技术领域
，尤其是涉及一种基于新型相关函数约束的CUDA架构并行优化立体变形测量方法及系统。背景技术
基于视觉测量模式的立体变形信息测量为被测目标的性能分析、几何形变监测、承载能力评估等提供了数据参考与支撑。然而在立体变形测量过程中运算数据规模较大，例如为了提高变形测量的准确性，通常需要大量的散斑点参与配准运算，因此配准的数据量规模较大，同时随着图像帧数与图像尺寸的增加，散斑子区的相关运算量以及亚像素的插值计算量都会大幅度增大，巨大的耗时降低了三维变形测量效率。导致当前要实现基于视觉测量模式的试件三维变形高精度测量，均需要以时间和较大计算量为代价，仅在算法上的优化和改进难以使解算速度获得跨数量级的提升。因此需要对变形测量的理论算法进行改进，同时需要结合改进后算法的特点与复杂度，搭建最优的并行耦合运算机制，提高基于视觉测量模式的立体变形测量效率。
当前在基于数字图像相关方法的变形测量中，主要存在两个关键问题。问题一是不同相机对间的散斑立体配准问题，如果依旧采用传统相关函数作为不同相机间散斑图像子区相似度的评判标准找到最大相关系数所对应的像点坐标，那么在立体匹配过程中需要搜索整幅目标图像，不仅降低了搜索效率，而且由于视差的存在，直接通过二维相关计算得到的同名散斑点亚像素坐标会存在较大的偏差，将会严重影响后续散斑点空间三维坐标的重构精度。问题二是关于散斑立体配准的硬件加速机理问题，目前基于CUDA架构进行散斑子区配准算法分解与并行优化的研究依然鲜少，这些方法中一般都是采用硬件组合的方式进行处理层的加速。例如，国外部分研究人员通过分析数字图像相关变形测量中相关系数的并行运算问题，采用GPU和FPGA相结合的平台模式，解决散斑匹配受限于运算资源的问题；还有一些学者通过多图形处理器寻求相关函数匹配的峰值分布，减少了总的运行时间。此外，一些方法采用CUDA Jacket模型进行散斑匹配的并行加速计算，但没有考虑算法的复杂度和适应性，导致容易存在循环拆分不当导致内存溢出等局限性问题。国内对于这方面的专业研究更为稀缺，经调研发现少量研究人员通过多线程实现变形测量中散斑时序匹配的加速研究，并取得了一定的加速效果。然而以上方法都是针对视觉测量模式中散斑子区二维匹配的加速，并没有考虑和进行图像对间散斑子区立体配准的并行耦合运算，同时也没有给出基于CUDA架构散斑立体配准与时序配准并行计算的具体实现方式和优化方案。本发明因此而来。
发明内容
1、本发明的目的
针对上述基于视觉测量模式的立体变形信息测量中，散斑点配准的数据量规模较大，同时随着图像帧数与图像尺寸的增加，散斑子区的相关运算量以及亚像素的插值计算量都会大幅度增大，巨大的耗时降低了三维变形测量效率的技术问题，提出了一种基于新型相关函数约束的CUDA架构并行优化立体变形测量方法及系统，通过多相机间的联合约束关系建立新型相关函数，能够将图像对间立体配准点的搜索限制在对极线附近区域，而不是整幅图像上，从而缩小了搜索空间，大大提高了三维变形测量效率和精度。
2、本发明所采用的技术方案
一种基于新型相关函数约束的CUDA架构并行优化立体变形测量方法，包括以下步骤：
S01：获得多相机间的相对位姿参数关系；
S02：基于相机间的联合约束关系建立散斑立体配准的新型相关函数，所述新型相关函数将图像对间立体配准点的搜索限制在对极线附近区域，对不同序列图像中采样散斑点进行时序配准和立体配准；
S03：对散斑点的空间坐标进行三维重构。
优选的技术方案中，所述步骤S01中获得多相机间的相对位姿参数关系的方法包括：
S11：设定测量系统共有η+1个相机组成，主相机为ξ，其余η个均为从相机，则根据透视投影原理建立基于空间点P的主相机共线性方程：
λξ[uξ vξ 1]＝Mξ[Rξ tξ][Xw Yw Zw 1]T (3)
其中，λξ是主相机投影尺度系数；(uξ,vξ)为主相机采集图片中特征点的像素坐标；[Xw Yw Zw 1]T是空间三维点齐次坐标；Mξ为主相机的内部参数；(Rξ,tξ)为主相机的外参矩阵；
S12：通过式(3)对相机外参和空间点三维坐标分量分别求一阶偏导，得到共线性误差方程：
其中，Ξξ是主相机共线性误差的偏差值；Bξ,Cξ分别是(uξ,vξ)T对相机外部参数和空间点三维坐标分量的一阶偏导数；是旋转矩阵Rξ对应的欧拉角；txξ,tyξ,tzξ是主相机平移矩阵的元素，[δξc δξp]T为主相机的外参修正值与三维坐标修正值，即每次迭代的步长：
式(5)中，分别是相机外部参数一阶偏导数的增量；ΔX,ΔY,ΔZ分别是三维坐标分量的一阶偏导数；
S13：矩阵γξ表示特征点的实际坐标与用共线性方程计算得到的重投影坐标之间的差值：
S14：根据主从相机间的相对外参关系，第η个相机的共线性方程表示为：
λη[uη vη 1]T＝Mη[Rξη Rξ tξη+Rξηtξ][Xw Yw Zw 1]T (7)
式(7)中，λ是从相机投影尺度系数；(uη,vη)为第η个相机下图像特征点像素坐标；Mη为第η个相机的内参矩阵；(Rξη，tξη)为从相机的外参矩阵；
S15：同理得到第η个从相机的共线性误差方程：
式(8)中，Ξη是从相机共线性误差的偏差值；Bη,Cη分别是从相机图像坐标对外部参数和空间点坐标的一阶偏导数；[δηc δηp]T分别是从相机的外参修正值与三维坐标修正值；矩阵γη表示特征点的实际坐标与用共线性方程计算得到的重投影坐标之间的差值；
S16：记B＝[Bξ B1 B2 … Bη]T,C＝[Cξ C1 C2 … Cη]T，则Jacobi矩阵J＝[B C]，故法化矩阵表示为：
S17：运用Levenberg-Marquardt算法进行迭代求解，建立的法方程表示为：
其中，γ＝[γξ γ1 γ2 … γη]T；I是单位矩阵；
得到多相机网络外参相对关系矩阵。
优选的技术方案中，所述步骤S02中建立散斑立体配准的新型相关函数的方法包括：
S21：对选定的相机对进行散斑图像配准，右相机成像面的散斑搜索范围表示为：
Fp(x,e)＝k0x+e(emin≤e≤emax) (11)
式(11)中，x是左相机目标图像中任意一点的横坐标，k0是极线方程的斜率，e为极线与右相机图像中坐标系纵轴的交点坐标；emin和emax分别表示在右相机图像中极线附近搜索对应散斑点亚像素的下限和上限；
S22：在极线方程的上下界内浮动，寻求最优的相关系数峰值，接着将基于多相机联合约束优化得到的极线方程代入到零均值归一化最小平方距离相关函数中，则生成的新型相关函数表达式改写成关于xi′和e的约束形式：
式(12)中，f(xi,yj)是变形前参考图像在点(xi,yj)处的灰度值；k的含义，g(x′i,y′j)是变形后目标图像中对应的同名点(x′i,y′j)处的灰度值，fm是变形前子区平均灰度值；gm是变形后子区平均灰度值；M是配准子区中心点到窗口边界距离像素尺寸。
优选的技术方案中，所述步骤S02中对不同序列图像中采样散斑点进行时序配准和立体配准时采用基于CUDA架构实施并行加速，并构建新型相关函数的Kernel函数接口。
优选的技术方案中，所述并行加速方法包括：
S23：在启动并行程序前将散斑子区亚像素匹配数据导入到常量内存，图像数据导入到纹理内存；待散斑配准结束后，再将计算结果回传至主机端；主机端与设备端之间相互连接的并行连接点是Kernel函数；
S24：根据预先配置的线程块与线程数量启动Kernel函数进行散斑图像的整像素和亚像素搜索任务，同时返回匹配迭代收敛的散斑点的标志位；完成所有的线程计算后，将散斑匹配得到的数据回传至CPU内存区域，散斑点的亚像素搜索过程结束；数据传输与访问时，通过线程同步方式进行存储器的累加操作与存取。
优选的技术方案中，所述Kernel函数的处理步骤包括：
输入线程数、左右相机参考图像和目标图像序列、定义配准矩阵空间、校准矩阵空间；
下载预先完成计算的数据包，k0,e,emin,emax；
定义Kernel函数线程索引，初始化并行计算参数；
对新型相关函数进行分解，包括：创建Hessian矩阵的逆矩阵空间、创建Hessian矩阵空间、创建Jacobi矩阵空间、创建形函数参数空间；
对新型相关函数的散斑立体配准迭代求解，得到形函数向量和对应的目标图像散斑点坐标(x′i,y′j)；
阻塞线程，保证线程同步，释放内存空间；
输出时序配准数据和立体配准数据。
优选的技术方案中，所述步骤S24中对亚像素插值进行并行计算，通过创建2M+1个线程分别对2M+1列的散斑点进行亚像素插值并行计算，其中，2M+1是散斑配准子区窗口的尺寸。
本发明还公开了一种基于新型相关函数约束的CUDA架构并行优化立体变形测量系统，包括：
多相机网络外参相对关系矩阵获取模块，获得多相机间的相对位姿参数关系；
相关函数构建及配准模块，基于相机间的联合约束关系建立散斑立体配准的新型相关函数，所述新型相关函数将图像对间立体配准点的搜索限制在对极线附近区域，对不同序列图像中采样散斑点进行时序配准和立体配准；
重构模块，对散斑点的空间坐标进行三维重构。
本发明又公开了一种基于新型相关函数约束的CUDA架构并行优化立体变形测量设备，包括上述的基于新型相关函数约束的CUDA架构并行优化立体变形测量方法。
本发明又公开了一种存储介质，所述存储介质存储有上述的基于新型相关函数约束的CUDA架构并行优化立体变形测量方法。
3、本发明所采用的有益效果
(1)针对立体变形测量中散斑配准耗时较长问题，提出基于新型相关函数约束的CUDA异构并行优化立体变形测量方法。首先通过多相机间的联合约束关系建立新型相关函数，能够将图像对间立体配准点的搜索限制在对极线附近区域，而不是整幅图像上，从而缩小了搜索空间。然后采用NVCC编译CUDA源程序的GPU并行机制，给出基于新型相关函数立体配准的并行运算程序设计原则和散斑子区配准的优化方案。本发明中的技术耦合Kernel函数进行新型相关函数的分解与并行算法的设计，不受重载函数的限制。从散斑子区立体配准算法优化以及与硬件的耦合加速两方面，使立体变形测量方法运算性能达到一个较优状态。实验结果表明，基于新型相关函数约束的CUDA异构并行优化算法在保证测量精度和稳定性的情况下，使散斑时序配准和立体配准的计算加速比分别达到了20.39倍和17.87倍，为试件三维变形的实时测量与输出提供了技术铺垫与参考。
(2)本发明从理论算法上对基于视觉模式的变形测量方法进行了改进，减小了算法的复杂度，提高了散斑立体配准的精度和稳定性，这是基于图像纹理匹配中的一项关键技术。然后考虑改进后算法与硬件加速的耦合因素，采用NVCC编译CUDA源程序的GPU并行编程方式，解决了Mex脚本和不同编程语言交互时的障碍问题，同时不受重载函数的限制，实现了三维变形信息测量中散斑图像子区时序匹配和立体匹配的高效运算。此外，考虑到亚像素插值计算的独立性，实现了散斑图像亚像素行、列插值的并行计算模式，进一步提高了散斑匹配的亚像素搜索速度。由此，即便选用硬件配置不是较高的计算机系统，基于本发明中的技术依然可以获得很好的加速效果，得到准确、稳定的散斑子区配准数据和三维重建云图，减小了基于视觉变形快速测量系统的硬件造价和研制成本，对实现立体变形的快速测量与结果再现具有实际的工程应用价值。
附图说明
图1为本发明基于新型相关函数约束的CUDA架构并行优化立体变形测量方法的流程图；
图2为本实施例基于多相机网络联合约束优化的散斑立体配准搜索区域示意图；
图3为本发明基于新型相关函数约束的CUDA架构并行优化立体变形测量系统的原理框图；
图4为本发明图像子区匹配并行运算流程图；
图5为本发明时序匹配耗时对比示意图；
图6为本发明立体匹配耗时对比示意图；
图7为本发明不同时刻测试散斑点的位移曲线；
图8为本发明t＝0.8s时Z向位移云图；
图9为本发明t＝1.5s时Z向位移云图；
图10为本发明t＝2.0s时Z向位移云图。
具体实施方式
下面结合本发明实例中的附图，对本发明实例中的技术方案进行清楚、完整地描述。显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域技术人员在没有做创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。
下面将结合附图对本发明实例作进一步地详细描述。
实施例1
如图1所示，一种基于新型相关函数约束的CUDA架构并行优化立体变形测量方法，包括以下步骤：
S01：获得多相机间的相对位姿参数关系；
S02：基于相机间的联合约束关系建立散斑立体配准的新型相关函数，所述新型相关函数将图像对间立体配准点的搜索限制在对极线附近区域，对不同序列图像中采样散斑点进行时序配准和立体配准；
S03：对散斑点的空间坐标进行三维重构。
一较佳的实施例中，步骤S01中获得多相机间的相对位姿参数关系的方法包括：
S11：设定测量系统共有η+1个相机组成，主相机为ξ，其余η个均为从相机，则根据透视投影原理建立基于空间点P的主相机共线性方程：
λξ[uξ vξ 1]＝Mξ[Rξ tξ][Xw Yw Zw 1]T (3)
其中，λξ是主相机投影尺度系数；(uξ,vξ)为主相机采集图片中特征点的像素坐标；[Xw Yw Zw 1]T是空间三维点齐次坐标；Mξ为主相机的内部参数；(Rξ,tξ)为主相机的外参矩阵；
S12：通过式(3)对相机外参和空间点三维坐标分量分别求一阶偏导，得到共线性误差方程：
其中，Ξξ是主相机共线性误差的偏差值；Bξ,Cξ分别是(uξ,vξ)T对相机外部参数和空间点三维坐标分量的一阶偏导数；是旋转矩阵Rξ对应的欧拉角；txξ,tyξ,tzξ是主相机平移矩阵的元素，[δξc δξp]T为主相机的外参修正值与三维坐标修正值，即每次迭代的步长：
式(5)中，分别是相机外部参数一阶偏导数的增量；ΔX,ΔY,ΔZ分别是三维坐标分量的一阶偏导数；
S13：矩阵γ表示特征点的实际坐标与用共线性方程计算得到的重投影坐标之间的差值：
S14：根据主从相机间的相对外参关系，第η个相机的共线性方程表示为：
λη[uη vη 1]T＝Mη[RξηRξ tξη+Rξηtξ][Xw Yw Zw 1]T (7)
式(7)中，λη是从相机投影尺度系数；(uη,vη)为第η个相机下图像特征点像素坐标；Mη为第η个相机的内参矩阵；(Rξη，tξη)为从相机的外参矩阵；
S15：同理得到第η个从相机的共线性误差方程：
式(8)中，Ξη是从相机共线性误差的偏差值；Bη,Cη分别是从相机图像坐标对外部参数和空间点坐标的一阶偏导数；[δηc δηp]T分别是从相机的外参修正值与三维坐标修正值；矩阵γη表示特征点的实际坐标与用共线性方程计算得到的重投影坐标之间的差值；
S16：记B＝[Bξ B1 B2 … Bη]T,C＝[Cξ C1 C2 … Cη]T，则Jacobi矩阵J＝[B C]，故法化矩阵表示为：
S17：运用Levenberg-Marquardt算法进行迭代求解，建立的法方程表示为：
其中，γ＝[γξ γ1 γ2 … γη]T；I是单位矩阵；
得到多相机网络外参相对关系矩阵。
一较佳的实施例中，步骤S02中建立散斑立体配准的新型相关函数的方法包括：
S21：对选定的相机对进行散斑图像配准，右相机成像面的散斑搜索范围表示为：
Fp(x,e)＝k0x+e(emin≤e≤emax) (11)
式(11)中，x是左相机目标图像中任意一点的横坐标，k0是极线方程的斜率，e为极线与右相机图像中坐标系纵轴的交点坐标；emin和emax分别表示在右相机图像中极线附近搜索对应散斑点亚像素的下限和上限；所限定的搜索区域如图2所示。
S22：图2中通过[emin,emax]设定了灰色搜索区域，在极线方程的上下界内浮动，寻求最优的相关系数峰值，接着将基于多相机联合约束优化得到的极线方程代入到零均值归一化最小平方距离相关函数中，则生成的新型相关函数表达式改写成关于x′i和e的约束形式：
式(12)中，f(xi,yj)是变形前参考图像在点(xi,yj)处的灰度值；k的含义，g(x′i,y′j)是变形后目标图像中对应的同名点(x′i,y′j)处的灰度值，fm是变形前子区平均灰度值；gm是变形后子区平均灰度值；M是配准子区中心点到窗口边界距离像素尺寸。
一较佳的实施例中，步骤S02中对不同序列图像中采样散斑点进行时序配准和立体配准时采用基于CUDA架构实施并行加速，并构建新型相关函数的Kernel函数接口。
并行加速方法的一个具体实现方案可以为包括以下步骤：
S23：在启动并行程序前将散斑子区亚像素匹配数据导入到常量内存，图像数据导入到纹理内存；待散斑配准结束后，再将计算结果回传至主机端；主机端与设备端之间相互连接的并行连接点是Kernel函数；
S24：根据预先配置的线程块与线程数量启动Kernel函数进行散斑图像的整像素和亚像素搜索任务，同时返回匹配迭代收敛的散斑点的标志位；完成所有的线程计算后，将散斑匹配得到的数据回传至CPU内存区域，散斑点的亚像素搜索过程结束；数据传输与访问时，通过线程同步方式进行存储器的累加操作与存取。
一较佳的实施例中，Kernel函数的处理步骤包括：
输入线程数、左右相机参考图像和目标图像序列、定义配准矩阵空间、校准矩阵空间；
下载预先完成计算的数据包，k0,e,emin,emax；
定义Kernel函数线程索引，初始化并行计算参数；
对新型相关函数进行分解，包括：创建Hessian矩阵的逆矩阵空间、创建Hessian矩阵空间、创建Jacobi矩阵空间、创建形函数参数空间；
对新型相关函数的散斑立体配准迭代求解，得到形函数向量和对应的目标图像散斑点坐标(x′i,y′j)；
阻塞线程，保证线程同步，释放内存空间；
输出时序配准数据和立体配准数据。
一较佳的实施例中，步骤S24中对亚像素插值进行并行计算，通过创建2M+1个线程分别对2M+1列的散斑点进行亚像素插值并行计算，其中，2M+1是散斑配准子区窗口的尺寸。
如图3所示，本发明还公开了一种基于新型相关函数约束的CUDA架构并行优化立体变形测量系统，包括：
多相机网络外参相对关系矩阵获取模块10，获得多相机间的相对位姿参数关系；
相关函数构建及配准模块20，基于相机间的联合约束关系建立散斑立体配准的新型相关函数，所述新型相关函数将图像对间立体配准点的搜索限制在对极线附近区域，对不同序列图像中采样散斑点进行时序配准和立体配准；
重构模块30，对散斑点的空间坐标进行三维重构并拟合变形前后云图。
下面通过完整的实施例，对本发明的测量系统的具体技术方案作进一步具体的说明。测量系统的工作流程包括以下步骤：
步骤一：基于相机间的联合约束关系建立散斑立体配准的新型相关函数。
步骤11：设定测量系统共有η+1个相机组成，主相机记为ξ，其余η个均为从相机，则根据透视投影原理可以建立基于空间点P的主相机共线性方程：
λξ[uξ vξ 1]＝Mξ[Rξ tξ][Xw Yw Zw 1]T (3)
其中，λξ是主相机投影尺度系数；(uξ,vξ)为主相机采集图片中特征点的像素坐标；[Xw Yw Zw 1]T是空间三维点齐次坐标；Mξ为主相机的内部参数；(Rξ,tξ)为主相机的外参矩阵。
步骤12：通过式(3)对相机外参和空间点三维坐标分量分别求一阶偏导，得到共线性误差方程：
其中，
此外，Ξξ是主相机共线性误差的偏差值；Bξ,Cξ分别是(uξ,vξ)T对相机外部参数和空间点三维坐标分量的一阶偏导数；是旋转矩阵Rξ对应的欧拉角；txξ,tyξ,tzξ是主相机平移矩阵的元素。[δξc δξp]T为主相机的外参修正值与三维坐标修正值，即每次迭代的步长：
式(5)中，分别是相机外部参数一阶偏导数的增量；ΔX,ΔY,ΔZ分别是三维坐标分量的一阶偏导数。
步骤13：矩阵γξ表示特征点的实际坐标与用共线性方程计算得到的重投影坐标之间的差值：
步骤14：根据主从相机间的相对外参关系，第η个相机的共线性方程可以表示为：
λη[uη vη 1]T＝Mη[RξηRξ tξη+Rξηtξ][Xw Yw Zw 1]T (7)
式(7)中，λη是从相机投影尺度系数；(uη,vη)为第η个相机下图像特征点像素坐标；Mη为第η个相机的内参矩阵；(Rξη，tξη)为从相机的外参矩阵。
步骤15：同理可以得到第η个从相机的共线性误差方程：
式(8)中，Ξη是从相机共线性误差的偏差值；Bη,Cη分别是从相机图像坐标对外部参数和空间点坐标的一阶偏导数；[δηc δηp]T分别是从相机的外参修正值与三维坐标修正值；矩阵γη表示特征点的实际坐标与用共线性方程计算得到的重投影坐标之间的差值。
步骤16：记B＝[Bξ B1 B2 … Bη]T,C＝[Cξ C1 C2 … Cη]T，则Jacobi矩阵J＝[B C]，故法化矩阵可以表示为：
步骤17：运用Levenberg-Marquardt算法进行迭代求解，建立的法方程可以表示为：
其中，γ＝[γξ γ1 γ2 … γη]T；I是单位矩阵。
步骤18：完成对式(10)法方程的求解之后，可以得到多相机网络外参相对关系的精确矩阵，结合对极几何原理，使得不同相机成像的图像对间同名散斑点的立体配准搜索区域被限制在极线约束的邻近范围内。对于选定的相机对进行散斑图像配准，右相机成像面的散斑搜索范围可表示为：
Fp(x,e)＝k0x+e(emin≤e≤emax) (11)
式(11)中，x是左相机目标图像中任意一点的横坐标，k0是极线方程的斜率，e为极线与右相机图像中坐标系纵轴的交点坐标；emin和emax分别表示在右相机图像中极线附近搜索对应散斑点亚像素的下限和上限，所限定的搜索区域如图2所示。
步骤19：图2中通过[emin,emax]设定了黄色搜索区域，在极线方程的上下界内浮动，寻求最优的相关系数峰值。接着将基于多相机联合约束优化得到的极线方程代入到零均值归一化最小平方距离相关函数中，则生成的新型相关函数表达式可以改写成关于x′i和e的约束形式：
式(12)中，f(xi,yj)是变形前参考图像在点(xi,yj)处的灰度值；k的含义g(x′i,y′j)是变形后目标图像中对应的同名点(x′i,y′j)处的灰度值。fm是变形前子区平均灰度值；gm是变形后子区平均灰度值；M是配准子区中心点到窗口边界距离像素尺寸。所建立的新型相关函数能够缩小相机对间散斑立体配准的搜索区域，提高了搜索效率，在一定程度上提高了三维变形测量的整体速度。
步骤二：设计采用CUDA架构并基于新型相关函数约束的散斑子区立体配准的并行实现与优化方案，同时构建新型相关函数的Kernel函数接口。
步骤21：在三维变形测量中，相关系数的计算数据占据了总体计算量的绝大部分，利用基于CUDA架构实施的并行加速机制，其目的就是要实现散斑点匹配的高速计算。基于CUDA架构的散斑并行匹配流程图如图4所示，散斑图像子区的亚像素搜索流程包括串行处理和并行计算两大部分。
步骤22：由于CPU端与显存模块之间的数据传递会降低GPU加速的整体效率，所以在启动并行程序前将散斑子区亚像素匹配数据导入到常量内存，图像数据导入到纹理内存；待散斑配准结束后，再将计算结果回传至主机端。
步骤23：主机端与设备端之间相互连接的并行连接点是Kernel函数，本发明设计的基于新型相关函数散斑立体匹配的Kernel函数实现如表1所示。
表1基于新型相关函数的散斑立体匹配Kernel函数实现
步骤24：最后根据预先配置的线程块与线程数量启动Kernel函数进行散斑图像的整像素和亚像素搜索任务，同时返回匹配迭代收敛的散斑点的标志位。完成所有的线程计算后，将散斑匹配得到的数据回传至CPU内存区域，散斑点的亚像素搜索过程结束。在调用mexw64脚本文件的程序中通过计时语句来记录并行计算部分的运行时间。
步骤25：对亚像素插值也进行并行计算。在目标图像的散斑子区匹配过程中，各个散斑点间亚像素位置处的灰度值插值运算互不影响，可以通过创建2M+1个线程分别对2M+1列的散斑点进行亚像素插值并行计算。其中，2M+1是散斑配准子区窗口的尺寸。这样在并行运算模式下通过2M+1次计算便可达到串行模式下(2M+1)2次的计算量。
步骤26：基于CUDA架构的并行程序设计和执行流程对运算速率有很大的影响，因此在图像散斑子区匹配优化方案的设计中，涉及数据传输与访问时，通过线程同步方式来解决存储器的累加操作与存取的先后性问题：尤其当并行程序需要进行频繁的局部空间操作时，合理的使用纹理内存可以增强程序的整体性能。
另一实施例中，本发明又公开了一种基于新型相关函数约束的CUDA架构并行优化立体变形测量设备，包括上述的基于新型相关函数约束的CUDA架构并行优化立体变形测量方法。
另一实施例中，本发明又公开了一种存储介质，所述存储介质存储有上述的基于新型相关函数约束的CUDA架构并行优化立体变形测量方法。
性能测试实验结果及分析如下：
用于测试的计算机配有IntelXeon E5-2620处理器，主频为2.10GHz，运行内存是32GB，安装CUDAToolkitV8.0。GPU是NVIDIAQuadro P2000系列显卡，内置1024个流处理器，显存位宽是128bit，具有5GB的显存容量。首先需要通过NVIDIAGPU的NVCC编译器在终端或通过函数调用，进行CUDA文件代码和文件目标的编译，供后续的Mex或者Mexcuda调用。在后缀为cu的文件中完成并行计算代码编写后，通过Mexcuda指令将其编译生成一个Mex-Function文件，供后续调用实现GPU并行计算。可以在Visual Studio工程中创建CU文件，并在Build Customizaiton中选择CUDA文件。在实际测量过程中，所采用的双三次样条插值法是以待插值区域中心的7×7窗口为计算单元，对散斑图像行与列的亚像素位置处的灰度值进行求解。
通过基于CPU运行的数字图像相关匹配算法与CUDA异构并行匹配运算分别对相同的散斑ROI区域进行测试，进行计算效率上的对比。依次选取5、20、50、100、150、200、250、300个测试点分别进行时序亚像素搜索和立体亚像素搜索，匹配运行时间如图5和图6所示。
图5和图6中的结果显示，当参与匹配的测试点数量较少时，基于CPU的串行匹配方式与CUDA异构并行运算耗时差别不大，并行加速比较小。这是因为参与匹配运算的数据量较小，还没有完全处于CPU计算能力的饱和状态；而启动GPU端核函数进行并行运算时，计算数据需要先从主机端传输到设备端，运算完成后还需将结果传回主机端，当数据规模较小时，并行计算的效率不足以掩盖访存和数据传输的延迟，因此，并行加速效果不明显。但随着匹配点个数的不断增加，通过传统串行运算模式的耗时增长趋势较大，而并行运算耗时的增长却较为缓慢。耗费在GPU上的非计算时间占算法执行总时间的比例越来越小，对于时序匹配耗时和立体匹配耗时，基于CUDA的异构并行计算速度要远远高于传统串行计算方式。而且随着匹配点数规模不断增大，加速效果越明显。当匹配点数为300时，通过本发明技术使散斑子区时序匹配加速比达到20.39倍，立体匹配加速比达到17.87倍，加速效果显著。
接着验证并行加速优化方案是否对变形测量精度造成影响，搭建测量系统硬件结构，试件变形是通过外力来实现，将试件的左侧通过螺栓刚性连接在目标模拟器滑块上，右侧固定在台体上，这样通过目标模拟器的运动来实现模拟翼面的动态变形。相机的采样频率为5Hz，曝光时间为7.07ms。目标模拟器的运动速度为2cm/s，运动时间为2s。
通过等间距选取一定线条上的散斑点进行空间三维重构，由于目标模拟器进行的是水平方向上的运动，因此针对性的解算出深度方向上的变形量，得到Z方向上的曲线如图7所示。此外，通过全站仪对试件变形量进行采样验证，得到的校准数据与图7中的测量结果是一致的，由此说明了本发明中的技术在实现速度变形测量的基础上，能够保证测量结果的精度和稳定性。
为了方便分析和重现实际效果，将动态采集图像中t＝0.8s，t＝1.5s和t＝2.0s时刻的散斑图在Z方向上的变形量进行拟合，得到不同时刻的拟合结果云图如图8、图9、图10所示。从图中可以看出，Z方向上的位移云图与实际的变形情况是相吻合的。
综上所述，本发明针对立体变形测量过程中数据规模较大，导致实现高精度变形测量需要以时间和较大计算量为代价，仅在算法上的优化和改进难以使解算速度获得跨数量级的提升等问题，通过从散斑子区立体配准算法优化以及与硬件的耦合加速两方面，使立体变形测量方法运算性能达到一个较优状态，实现对三维变形测量中散斑图像子区时序匹配和立体匹配的高效并行运算，进而实现了基于视觉非接触测量模式的立体变形高精度、快速测量功能。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。
