基于GPU的Zernike矩快速计算方法

摘要

基于GPU的Zernike矩快速计算方法属正交矩加速技术领域，本发明提出存径向多项式系数结合图像重布局的八卦限对称性算法，该混合算法能显著地加快任意单个矩的计算速度。提出在一个内核Kernel中组合Block的组包方案，加快一族或一组矩的计算。而合并内核的组包方案能进一步缩短一组或者一族图像的矩的计算时间，特别是能克服小尺寸图像中由于计算量小，不能有效实施混合算法的瓶颈。本发明能推动Zernike矩在实时图像处理和模式识别等领域中的发展，可应用于视频图像水印、镜头分割、光学系统中波前重建、机器视觉系统中目标定位和识别检测，其它正交矩的快速计算提供了非常有价值的参考。

说明书

技术领域

本发明属正交矩加速技术领域，具体涉及Zernike矩理论，以及基于GPU的快速计算方法。

背景技术

图像是非常有用的信息媒介和交流工具，它能够以一种紧凑而高效的的方式来表达和处理复杂场景。矩不变理论作为描述图像的特征广泛应用于计算机视觉、遥感处理、医学成像、模式识别及图像处理等领域。矩理论源自数学、物理学和统计学中，是一个用于表征函数及捕捉函数重要特性的标量。矩描述子是基于区域的图像特征，因为它们使用图像的所有信息，即图像轮廓和它的内容。不同于仅基于轮廓的描述符，如使用图像边界信息的傅立叶描述符。

总结来说，现有的图像矩主要有非正交矩与正交矩两种。非正交(几何矩与复数矩)矩将图像投影到一组非正交的函数多项式上。正交矩是将图像投影到一组正交多项式上。现有的正交矩分为基于矩形正交矩和基于圆的正交矩两种。基于矩形正交矩，如Legendre矩、Tchebichef矩、小波矩、Krawtchouk矩定义在笛卡尔坐标系下，其几何不变性，特别是旋转不变性并不成立。基于圆的正交矩，如Zernike矩、伪Zernike矩、Fourier-Mellin矩定义在单位圆上，其幅值具有本质上的旋转不变性。该特性使Zernike矩具有稳定的数值性质和良好的重建能力，因而在实际中广泛应用。

但是，Zernike矩的定义复杂。其中，Zernike径向多项式是Zernike矩的核心，它是定义在单位圆内完备的正交集，含有复杂的阶乘和幂函数运算。此外，还需要对图像的每一个像素进行多项式的映射。因此，对它们的计算是非常耗时的。同时，在实际应用中通常需要一组或一族矩描述图像的特征。如果待处理的图像较大时，如一些遥感图像；或者需要计算高阶矩，如一些医学图像；或者需要高精度的矩，如对图像水印技术或重建，矩的计算时间更长。另一方面，在工程应用中需要快速计算矩以适应实时系统的应用。例如：视频图像水印、视频监控中人脸识别、在线产品缺陷检测等。为消除由于计算时间而引起的Zernike矩的应用的限制，需要对Zernike矩的快速计算方法进一步研究。

当前对Zernike矩的加速主要集中在CPU中进行算法的改进，和直接定义法比较其加速比最大到几十倍。GPU(GraphicsProcessingUnits)的架构由大量简单的处理单元组成，由于GPU对并行数据加速的优势，将GPU作为加速硬件的需求日俱增加。图像处理算法通常具有数据量大以及计算访存密集的特点，因此，GPU被广泛用于图像处理及模式识别中。虽然GPU提供很高的计算能力，使用定义加速Zernike矩相对于CPU的加速比能达到几百甚至上千。但相对于GPU程序实现，优化过程却复杂得多，因为不仅要考虑算法特征，还要了解底层硬件架构的特征，以获得这两种特征的高效映射。因此，需要研究Zernike矩的算法设计、GPU的优化程序方法、存储器使用和程序指令等，进一步提高Zernike矩计算的加速比。Zernike矩的快速计算方法研究对于图像处理的发展有着推动的作用，其取得的成果可直接应用于图像水印、镜头分割、光学系统中波前重建、机器视觉系统中目标定位和识别检测。

发明内容

本发明的目的在于提供一种能加快Zernike矩运算速度的计算方法。有效消除了阶乘计算的精度限制，提高了Zernike矩的计算阶数。在GPU的算法设计中避免了大量条件语句引起的线程分歧，并消除小尺寸图像由于占用率低引起的计算瓶颈。

本发明的基于GPU的Zernike矩快速计算方法包括两部分内容:基于GPU的存径向多项式系数结合图像重布局的八卦限对称性算法的混合算法和基于GPU的内核合并的组包方案。

1.所述的基于GPU的存径向多项式系数结合图像重布局的八卦限对称性算法的混合算法包括下列步骤：

1.1将笛卡尔坐标下N×N图像的坐标[i,j]，映射到单位内切圆[x,y]的映射变换，归一化以后的坐标为：

其中，i＝tid％N，k＝tid/N。tid是线程索引，其值从0到(N×N-1)。x轴、y轴、直线y＝x和直线y＝-x将单位圆分为8个象限，称为八卦限。用8个一维数组h₁,h₂,h₃…h₈表示八卦限中像素的值。利用第1卦限中的地址索引值计算对应的极径ρ和极角θ的值，即0≤ρ≤1，将计算结果放入2个一维数组，存储在全局内存中。

1.2将步骤1.1中的所有一维数组中对角线上重复的像素值置0。本算法中将h₂,h₃,h₆,h₇中对角线上的像素置0，这样避免了在对角线上重复的像素叠加导致错误的结果。

1.3计算径向多项式系数将计算结果存储到常量缓存中。其中，n为阶数(order)，m为角频(repetitions)，n是非负正整数，m是整数，且满足n-|m|为偶数，|m|＜n，s为[0-(n-|m|)/2]的整数。

1.4通过4个Kernel在GPU中计算n阶m角频的Zernike矩，计算过程分为下列4个步骤:

1.4.1在GPU的Kernel1中，将八卦限的灰度值按照相同的索引地址是对称点的规则，重新载入到8个一维的数组h₁,h₂,h₃…h₈。此步骤和步骤1.2构成重布局后的图像信息。不同卦限的灰度值被重新布局在8个一维的数组，确保了8个数组中相同的索引地址自然地是对称点。因此，在后续的计算中，不需要任何条件语句去获取对称点，进而避免产生线程分歧。其次，所有卦限数据转为一维数组后，确保了GPU中图像数据的访问是对齐的和连续的。

1.4.2在GPU的Kernel2中，将步骤1.1中的第一卦限的极角θ通过查找表取出，计算八卦限映射和

1.4.3在GPU的Kernel3中，将步骤1.1的第1卦限中的极径ρ和步骤1.3的径向多项式系数通过查找表的方式取出，代入径向多项式：

1.4.4在GPU的Kernel3中，将步骤1.4.3径向多项式计算结果和步骤1.4.2八卦限映射与相乘求Zernike矩的映射；

1.4.5在GPU的Kernel4中将步骤1.4.4的Zernike矩的映射使用并行归约求和，并将最终结果输出到CPU。

2.所述的基于GPU的内核合并的组包方案，采用了顺序执行模式，计算第n阶m个角频的一族Zernike矩时，包括下列步骤：

2.1重复步骤1.1至1.3中的计算；

2.2重复步骤1.4.1的Kernel1计算步骤；

2.3重复步骤1.4.2的Kernel2计算步骤；

2.4在Kernel3中，采用顺序执行模式，保留原线程块Block和线程束Warp调用方式，通过线程块Block组包，将n阶{floor(n/2)+1}个Zernike矩的八卦限映射与径向多项式相乘的计算合并到一起，由图像尺寸和阶数n决定组包线程块数量；

2.5在Kernel4中，同样用线程块Block组包方案，使用并行归约将n阶的{floor(n/2)+1}个Zernike矩的映射求和，结果输出到CPU。

本发明通过对Zernike矩快速计算的研究，提出存径向多项式系数结合图像重布局的八卦限对称性算法，该混合算法能加快任意单个矩的计算速度。提出在一个内核Kernel中组合Block的组包方案，加快一族或一组矩的计算。实验证明在GPU中混合算法对大尺寸图像的Zernike矩计算加速效果显著。而合并内核的组包方案能进一步缩短一组或者一族图像的矩的计算时间，特别是能克服小尺寸图像中由于计算量小，不能有效实施混合算法的瓶颈。应用本发明提出的算法最终达到了显著的加速效果，能推动Zernike矩在图像处理和模式识别等领域中的发展，同时为其它正交矩的快速计算提供了非常有价值的参考。

附图说明

图1为含有图像索引的内切圆映射示意图

图2为8×8图像的重布局和第1卦限的极坐标排列

图3为CPU执行switch控制指令启动GPUKernel计算八卦限映射的伪代码

图4为Zernike矩映射zr-map和zi-map的伪代码

图5为Zernike矩计算流程图

图6为顺序执行模式

图7为64×64图像阶数n＝8的一族Zernike矩的组包方案示意图

图8为组包方案伪代码

图9为小尺寸图像实验对比

图10为大尺寸图像实验对比

具体实施方式

下面结合附图对本发明进行描述。

1.存径向多项式系数结合图像重布局的八卦限对称性算法的混合算法

1.1存储第一卦限的极径和极角值

图1说明了笛卡尔坐标下一个N×N图像坐标[i，j]映射到内切单位圆[x，y]的映射变换。每个方格内的数字代表了像素的地址索引。归一化以后的坐标为：

其中，i＝tid％8，k＝tid/8,tid是线程索引，其值从0到63。按照图2(b)所示的读取顺序载入索引值，利用图2(a)中第1卦限中的地址索引值计算对应的极径ρ和极角θ的值。其计算结果如图2(d)所示，将它们存储在全局内存中。

1.2存径向多项式系数

存储阶乘方法计算径向多项式虽然能节省阶乘计算时间，但存储的阶乘数的大小受控于计算的精度，单精度的只能算到18，而双精度只能算到42。一个更好的算法是避开存储阶乘所受到的精度限制，同时，还能进一步提高计算速度。该算法是将径向多项式系数

提前计算，然后存到常量存储器。由于预计算径向多项式系数时，分子分母阶乘项在提前计算时被约掉，有效地避开了大数量的阶乘运算，计算的结果适用于存储。同时，可以得到更高阶的Zernike矩。

1.3在GPU的Kernel1中，图像重布局；

在如图2(a)所示8×8图像中，第1卦限按照图2(b)的顺序载入像素信息到一维数组h₁。数组内对应的[h₁(0)...h₁(7)]是对应了地址(28,29,30,31,21,22,23,14)的像素值。为了确保载入第2卦限的像素和第1卦限对称，它的载入顺序必须遵循地址(28,20,12,4,21,13,5,14)。其余的6个卦限的载入顺序见图2(c)中的地址。观察图2(c)中每一列，会发现一维数组中相同地址索引所包含的像素信息在原图中恰好是按照八卦限对称。这样就避免了使用条件语句判断对称点而引起的线程分歧，8个卦限的像素灰度值按照这样的方式被重新布局到8个一维数组h₁,h₂,h₃…h₈中。同时，需要提前设置对角线像素中重复的像素值为0。如图2(c)将h₂,h₃,h₆,h₇中重复的像素值置0。这样对角线像素相当于只载入一次，避免了像素的重复叠加。

1.4在GPU的Kernel2中，计算八卦限映射和

在GPUKernel2中，用重布局的图像数据和步骤1.1中计算的第1卦限极角数组按照下面公式计算和为了避免GPU中使用条件语句判断4种情况产生的线程分歧，此处采用在CPU中执行switch控制指令，然后在CPU中启动GPUKernel。伪代码如图3所示，在CPU中当m除以4的余数分别是0，1，2，3时，顺次对应gr_gi_0,gr_gi_1,gr_gi_2和gr_gi_3，即下面公式在GPU执行的4种计算。与分别用gr和gi简化表示。

1.5计算径向多项式

在GPU的Kernel 3中，将步骤1.1中的第1卦限中的极径和步骤1.2的径向多项式系数通过查找表的方式取出，计算径向多项式

1.6在GPU的Kernel3中，将步骤1.5径向多项式和步骤1.4八卦限映射和相乘求Zernike矩的映射，伪代码如图4所示。

1.7在GPU的Kernel4中将步骤1.6的Zernike矩的映射使用并行归约求和，并将最终结果输出到CPU，上述的混合算法流程如图5所示。实验中GPU选择了NVIDIATeslaK40，它是当前比较流行的开普勒架构。在K40中使用混合算法对单个Zernike矩的计算时间及与直接法对比的加速比如表1所示。在K40中使用混合算法对Zernike一族矩的计算时间及和直接法对比的加速比如表2所示。在大尺寸图像中，所提混合算法有显著的加速比。

表1双精度K40中不同尺寸图像单个Zernike矩的计算时间及加速比(μs)

表2双精度K40中不同尺寸图像一族Zernike矩的计算时间及加速比(ms)

2.使用提出的基于GPU的内核合并的组包方案，进一步加速计算一组或者一族Zernike矩。

2.1如图7所示，首先重复步骤1.1、1.2、1.3、1.4和1.5。

2.2在Kernel3中，采用图6顺序执行模式，保留原线程块Block和线程束Warp调用方式。通过线程块Block组包，将n阶{floor(n/2)+1}个Zernike矩的八卦限映射与径向多项式相乘的计算合并到一起。图像尺寸和阶数n决定组包线程块数量。

2.3在Kernel4中，同样适用线程块Block组包，使用并行归约将n阶的{floor(n/2)+1}个Zernike矩的映射求和，结果输出到CPU。

图7中示意了阶数n为8的64×64图像Zernike矩计算中Kernel3和Kernel4的内核合并过程。此时共有5个计算任务，分别是Z₈₀，Z₈₂，Z₈₄，Z₈₆和Z₈₈。在内核Kernel中每个线程块Block含有128线程，共20个线程块，平均分配给每个计算任务4个Block线程块。其组包方案的伪代码如图8所示。将不同计算任务Block，组合在一个Kernel中，原来未利用的资源空间被安排一些线程块，执行不同计算任务，从而提高了并发度、资源利用率。

图像尺寸和阶数决定组包线程块数量。表3示意了一些尺寸下的线程块数。图9和图10可视化地展示了在不同尺寸下提出的所有算法的40阶以内所有矩的双精度累计计算时间，这些算法包括直接法，提出的混合算法，组包方案。可以看到，直接法的计算时间随着输入图像尺寸的增加而急剧增加，所提出的混合算法运算时间显著地减小了运算时间。而使用提出的组包方案，解决了小尺寸图像引起的GPU流处理器占用率低、资源闲置的情况的问题，特别对图像尺寸小于256小尺寸图像，其加速效果显著。

表3组包方案对于不同尺寸图像和不同阶数分配Block的块数