CN107480102B - 一种基于并行lu分解的psf反变换方法 - Google Patents

Abstract

本发明公开了一种基于并行LU分解的PSF反变换方法，包括：A1：输入二维化的PSF矩阵H与向量化的像平面图像P，分别对矩阵H和向量P进行分块；A2：初始化s＝0；A3：令s＝s+1，对矩阵H的子矩阵H_s‑s与n‑s个子矩阵H_i‑s进行第一变换运算；A4：对矩阵H的子矩阵H_i‑j进行第二变换运算；A5：判断s是否等于n‑1，如是则进入步骤A6，如否则返回步骤A3；A6：将矩阵H中处于主对角线以下的子矩阵集合作为方程组系数矩阵，将向量P作为方程组常数项，计算得到作为中间变量的分块向量Q；A7：将矩阵H中处于主对角线及主对角线以上的子矩阵集合作为方程组系数矩阵，将分块向量Q作为方程组常数项，计算得到物平面图像向量B。本发明能有效地减少PSF反变换过程所消耗的时间。

Description

一种基于并行LU分解的PSF反变换方法

技术领域

本发明涉及计算机视觉与数字图像处理领域，尤其涉及一种基于并行LU分解的PSF反变换方法。

背景技术

近年来，光场相机所捕获的光场图像在立体显示、三维重建、虚拟现实等计算机视觉领域引发了人们的广泛关注，对光场相机的建模分析也成为了研究重点，根据菲涅耳衍射定理，光场相机的理论模型往往都可以对应成一个线性方程组，其中方程组的系数矩阵即为二维化的联合PSF(点扩散函数)矩阵，方程组的未知量即为向量化的物平面图像信息，方程组的常数项即为向量化的成像平面图像信息。PSF反变换的目的就是通过求解这个线性方程组，从PSF矩阵与相机所捕获的图像，反推得到物平面上的真实图像信息。但是随着成像平面像素值的增加，方程组的维数也会急剧增加，这就导致利用传统求解方程组的方法很难求解出结果，或者求解的过程十分复杂，时间消耗非常严重。

在目前已有的求解PSF反变化的方法中：一类是直接求解联合PSF矩阵的逆矩阵，将逆矩阵左乘常数项，得出最终结果；该方案理论上比较直观，步骤上比较简洁，但是需要消耗的时间非常多，而且求解过程中对计算机的内存、CPU要求也较高；另一类则基于MATLAB自带的左除功能，该方案求解时间短，理论上较直观，但是适用平台较为局限，而且代码封装后可拓展性较差，求解速度受到限制，运算效率较差。

以上背景技术内容的公开仅用于辅助理解本发明的构思及技术方案，其并不必然属于本专利申请的现有技术，在没有明确的证据表明上述内容在本专利申请的申请日已经公开的情况下，上述背景技术不应当用于评价本申请的新颖性和创造性。

发明内容

为解决上述技术问题，本发明的主要目的在于提供一种基于并行LU分解的PSF反变换方法，能有效地减少PSF反变换过程所消耗的时间，提高运算效率。

为了达到上述目的，本发明采用以下技术方案：

本发明公开了一种基于并行LU分解的PSF反变换方法，包括以下步骤：

A1：输入二维化的PSF矩阵H与向量化的像平面图像P，分别对矩阵H和向量P进行分块，将矩阵H分为n×n个k阶子矩阵，将向量P分成n个含k个元素的子向量，其中n和k均为正整数；

A2：初始化s＝0；

A3：令s＝s+1，对矩阵H的子矩阵H_s-s与n-s个子矩阵H_i-s进行第一变换运算，其中s+1≤i≤n；

A4：对矩阵H的子矩阵H_i-j进行第二变换运算，其中s+1≤i≤n，s+1≤j≤n；

A5：判断s是否等于n-1，如果是，则进入步骤A6，如果否，则返回步骤A3；

A6：将矩阵H中处于主对角线以下的子矩阵集合作为方程组的系数矩阵，将向量P作为方程组的常数项，计算得到作为中间变量的分块向量Q；

A7：将矩阵H中处于主对角线及主对角线以上的子矩阵集合作为方程组的系数矩阵，将分块向量Q作为方程组的常数项，计算得到物平面图像向量B。

优选地，步骤A1中n与k的取值满足下述两式：

m＝n×k (1)

其中，m为分块前的矩阵H的阶数，n为分块后的矩阵H的阶数，k为分块后的矩阵H中每个子矩阵的阶数，s为用于并行计算的处理器个数，O₁为计算所采用的矩阵除法的时间复杂度，O₂为计算所采用的矩阵乘法的时间复杂度。

优选地，步骤A3中对矩阵H的子矩阵H_s-s与n-s个子矩阵H_i-s进行第一变换运算具体为：H_i-s＝H_i-s/H_s-s。

优选地，步骤A4中对矩阵H的子矩阵H_i-j进行第二变换运算具体为：H_i-j＝H_i-j-H_i-s×H_s-j。

优选地，步骤A6中具体包括：将矩阵H中处于主对角线以下的子矩阵集合补全为下三角矩阵L，并将该矩阵L作为方程组的系数矩阵，将向量P作为方程组的常数项，采用并行的高斯消元法计算得到作为中间变量的分块向量Q。

优选地，步骤A6中将矩阵H中处于主对角线以下的子矩阵集合补全为下三角矩阵L具体为：L＝H_d+E，其中，H_d为矩阵H中主对角线以下的子矩阵集合，E为n×k阶的单位矩阵，补全后的矩阵L为n×n个k阶子矩阵。

优选地，步骤A6中采用并行的高斯消元法计算得到分块向量Q具体包括：对子矩阵H_i-j与向量P作如下运算，其中j+1≤i≤n，并初始化j＝1；

H_i-j＝H_i-j×P_j (5)

P_i＝P_i-H_i-j (6)

j＝j+1 (7)

其中式(5)和式(6)的计算包括i取值为j+1≤i≤n的所有运算；

重复式(5)～(7)，直至j＝n为止，此时计算得到分块向量Q＝P。

优选地，步骤A7中具体包括：将矩阵H中处于主对角线及主对角线以上的子矩阵集合补全为上三角矩阵U，并将该矩阵U作为方程组的系数矩阵，将分块向量Q作为方程组的常数项，采用并行的高斯消元法计算得到物平面图像向量B。

优选地，步骤A7中将矩阵H中处于主对角线及主对角线以上的子矩阵集合补全为上三角矩阵U具体为：U＝H_u+Z，其中，H_u为矩阵H中主对角线及主对角线以上的部分，Z为n×k阶的零矩阵，补全后的矩阵U为n×n个k阶子矩阵。

优选地，步骤A7中采用并行的高斯消元法计算得到物平面图像向量B具体包括：对子矩阵H_i-j与分块向量Q作如下运算，其中1≤i≤j-1，并初始化j＝n；

Q_j＝Q_j/H_i-j (8)

H_i-j＝H_i-j×Q_j (9)

Q_i＝Q_i-H_i-j (10)

j＝j-1 (11)

其中式(9)和式(10)的计算包括i取值为1≤i≤j-1的所有运算；

重复式(8)～(11)，直至j＝1为止，此时计算得到物平面图像向量B＝Q。

与现有技术相比，本发明的有益效果在于：本发明的基于并行LU分解的PSF反变换方法利用并行化的LU分解法，通过并行的方法对PSF矩阵进行分块与变换，以将其改造成LU分解中两个系数矩阵的组合的形式，再通过两次计算得到相机物平面的有效图像信息，该方法充分利用了LU分解适用于并行化的特点，并行化程度高、内存复用率高，有效地减少了PSF反变换所消耗的时间，提高运算效率。

在进一步的方案中，进一步利用可以并行化的高斯消元法，从而更进一步减少PSF反变换中所消耗的时间，提高运算效率。

附图说明

图1是本发明优选实施例的基于并行LU分解的PSF反变换方法的流程示意图；

图2是本发明优选实施例的对矩阵H和向量P的分块示意图；

图3是步骤A3中第s层的计算过程示意图；

图4是图3中s＝1的计算过程示意图；

图5是步骤A4中第s层的计算过程示意图；

图6是图5中s＝1的计算过程示意图；

图7是本发明优选实施例的对矩阵H进行LU分解的运算过程示意图。

具体实施方式

下面对照附图并结合优选的实施方式对本发明作进一步说明。

为了提升PSF反变换的效率，尽量保持算法的可拓展性，如图1所示，本发明的优选实施例公开了一种基于并行LU分解的PSF反变换方法，本发明的思路是基于LU分解在求解线性方程组中的便捷性，利用矩阵分块的方法，使得LU分解中使用并行处理成为可能，并利用并行化的LU分解技术和高斯消元法，有效地减少了PSF反变换所消耗的时间；具体地包括以下步骤：

A1：输入二维化的联合PSF矩阵H与向量化的像平面图像P，分别对矩阵H和向量P进行分块，将矩阵H分为n×n个k阶子矩阵，将向量P分成n个含k个元素的子向量；

具体地，选择一个合适的分块尺寸k，分别对矩阵H和向量P进行分块，其中，本发明优选实施例针对的矩阵H为方阵，在PSF矩阵中，向量P的维度与矩阵H的列数相同，可以记矩阵H的行数、列数以及向量P的维度均为m，为保证矩阵H与向量P分块后，所有子矩阵阶数、自向量维数大小一致，分块尺寸需要是矩阵H行数、列数的公约数，即有m＝k×n，其中n为正整数。利用满足以上两点的分块尺寸k，将矩阵H分为n×n个k阶子矩阵，将向量P分成n个含k个元素的子向量，其中n和k均为正整数，分块结果如图2所示。

其中，分块尺寸k的选取，需考虑具体计算机的性能，分块尺寸越大，单次计算量越大，但循环次数、每层运算中的计算次数越少，在本实施例中，按照下述两式来对n和k进行取值较为合理，此时单次的计算量不算太大，循环次数、每层运算的运算次数也不太多：

m＝n×k (1)

其中，m为分块前的矩阵H的阶数，n为分块后的矩阵H的阶数，k为分块后的矩阵H中每个子矩阵的阶数，s为用于并行计算的处理器个数，O₁为计算所采用的矩阵除法的时间复杂度，O₂为计算所采用的矩阵乘法的时间复杂度。

A2：初始化s＝0；

A3：对矩阵H进行LU分解运算，令s＝s+1，在第s层运算中，对矩阵H的子矩阵H_s-s与n-s个子矩阵H_i-s进行第一变换运算，其中s+1≤i≤n；

具体地，基于矩阵除法与并行运算，对矩阵H的指定一列进行改造。为了将矩阵H改造成满足LU分解的形式，分层次地对矩阵H内子矩阵进行变换，在每一层计算中，需要对该层计算所在的列进行改造，在第s层计算中，选取主对角线上的H_s-s与n-s个子矩阵H_i-s(s+1≤i≤n)进行如下第一变换运算(并行除法计算)：

H_i-s＝H_i-s/H_s-s (3)。

因为在同一层运算中，不同i值对应的运算相互无关，所以可以使用并行计算来加速。计算过程如图3所示，其中s＝1时的计算过程如图4所示，其中斜线块代表H_s-s，十字线块代表H_i-s。

A4：对矩阵H的子矩阵H_i-j进行第二变换运算，其中s+1≤i≤n，s+1≤j≤n；

在第s层的运算中，并行地作如下第二变换运算：

H_i-j＝H_i-j-H_i-s×H_s-j (4)

其中，s+1≤i≤n，s+1≤j≤n；

在步骤A3中所提到的第s层运算中，在对特定一列进行改造后，还需要对部分其他子矩阵进行变换，具体如上式。

因为在同一层运算中，不同i、j值对应的运算互不无关，因此可以采用并行方法来加速计算。计算过程如图5所示，其中s＝1时的计算过程如图6所示，其中圆点块代表H_i-j，竖线块代表H_s-j，十字线块代表H_i-s。

A5：判断s是否等于n-1，如果是，则进入步骤A6，如果否，则返回步骤A3；

在本实施例中，对矩阵H进行LU分解运算总共分为n-1层，即重复步骤A3和A4，直至n-1层计算结束，循环过程如图7所示。

A6：将矩阵H中处于主对角线以下的子矩阵集合作为方程组的系数矩阵，将向量P作为方程组的常数项，计算得到作为中间变量的分块向量Q；

具体地，对矩阵H中处于主对角线以下的子矩阵集合，即可视为LU分解中的L矩阵内处于主对角线以下的部分，将其补全为下三角矩阵L，并作为方程组的系数矩阵，将向量P作为方程组的常数项，采用并行的高斯消元法，计算得到分块向量Q。

其中将矩阵H中处于主对角线以下的子矩阵集合补全为下三角矩阵L具体为：L＝H_d+E，其中，H_d为矩阵H中主对角线以下的子矩阵集合，E为m阶的单位矩阵，补全后的矩阵L为n×n个k阶子矩阵。

采用并行的高斯消元法具体而言，是对子矩阵H_i-j与向量P(j＝1，j+1≤i≤n)，作如下运算：

H_i-j＝H_i-j×P_j (5)

P_i＝P_i-H_i-j (6)

j＝j+1 (7)

其中式(5)和式(6)的计算包括i取值为j+1≤i≤n的所有运算；

重复式(5)～(7)，直至j＝n为止，此时计算得到分块向量Q＝P。

因为对于同一个j而言，式(5)、(6)在不同i之间的运算互不相关，因此可以采用并行的方法加速计算过程。

A7：将矩阵H中处于主对角线及主对角线以上的子矩阵集合作为方程组的系数矩阵，将分块向量Q作为方程组的常数项，计算得到物平面图像向量B。

具体地，对矩阵H中处于主对角线及主对角线以上的子矩阵集合，即可视为LU分解中的U矩阵内处于主对角线及主对角线以上的部分，将其补全为上三角矩阵U，并作为方程组的系数矩阵，将分块向量Q作为方程组的常数项，采用并行的高斯消元法，得到物平面图像向量B。

其中将矩阵H中处于主对角线及主对角线以上的子矩阵集合补全为上三角矩阵U具体为：U＝H_u+Z，其中，H_u为矩阵H中主对角线及主对角线以上的部分，Z为m阶的零矩阵，补全后的矩阵U为n×n个k阶子矩阵。

采用并行的高斯消元法具体而言，是对子矩阵H_i-j与分块向量Q(j＝n，1≤i≤j-1)，作如下运算：

Q_j＝Q_j/H_i-j (8)

H_i-j＝H_i-j×Q_j (9)

Q_i＝Q_i-H_i-j (10)

j＝j-1 (11)

其中式(9)和式(10)的计算包括i取值为1≤i≤j-1的所有运算；

重复式(8)～(11)，直至j＝1为止，此时计算得到物平面图像向量B＝Q。

因为对于同一个j而言，式(9)、(10)在不同i之间的运算互不相关，因此可以采用并行的方法加速计算过程。

本发明的优选实施例公开的是一种针对光场建模的基于并行LU分解的PSF反变换方法，首先将输入的二维化的PSF矩阵H与向量化的像平面图像P进行分块，通过分块可以减少后续并行运算中的循环次数，提高并行效率；然后引入并行LU分解方法，在PSF反变换中，将线性方程组的系数矩阵改造为传统LU分解中的L矩阵与U矩阵结合的形式；最后，在求解计算过程中，使用两次并行化的高斯消元法，得到最终的物平面图像向量；本发明充分利用了LU分解适用于并行化的特点，在运算过程中各个子矩阵的运算互不相关，能够采用并行运算来进行加速，因此，本发明的方法并行化程度高、内存复用率高，可以有效地减少PSF反变换中所消耗的时间，提高运算效率。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干等同替代或明显变型，而且性能或用途相同，都应当视为属于本发明的保护范围。

Claims (6)

1.一种基于并行LU分解的PSF反变换方法，其特征在于，包括以下步骤：

A1：输入二维化的PSF矩阵H与向量化的像平面图像P，分别对矩阵H和向量P进行分块，将矩阵H分为n×n个k阶子矩阵，将向量P分成n个含k个元素的子向量，其中n和k均为正整数，m＝n×k，m为分块前的矩阵H的阶数，n为分块后的矩阵H的阶数，k为分块后的矩阵H中每个子矩阵的阶数；

A2：初始化s＝0；

A3：令s＝s+1，对矩阵H的子矩阵H_s-s与n-s个子矩阵H_i-s进行第一变换运算：H_i-s＝H_i-s/H_s-s，其中s+1≤i≤n；

A4：对矩阵H的子矩阵H_i-j进行第二变换运算：H_i-j＝H_i-j-H_i-s×H_s-j，其中s+1≤i≤n，s+1≤j≤n；

A5：判断s是否等于n-1，如果是，则进入步骤A6，如果否，则返回步骤A3；

A6：将矩阵H中处于主对角线以下的子矩阵集合补全为下三角矩阵L，并将该矩阵L作为方程组的系数矩阵，将向量P作为方程组的常数项，采用并行的高斯消元法计算得到作为中间变量的分块向量Q；

A7：将矩阵H中处于主对角线及主对角线以上的子矩阵集合补全为上三角矩阵U，并将该矩阵U作为方程组的系数矩阵，将分块向量Q作为方程组的常数项，采用并行的高斯消元法计算得到物平面图像向量B。

2.根据权利要求1所述的PSF反变换方法，其特征在于，步骤A1中n与k的取值还满足下式：

其中，s为用于并行计算的处理器个数，O₁为计算所采用的矩阵除法的时间复杂度，O₂为计算所采用的矩阵乘法的时间复杂度。

3.根据权利要求1所述的PSF反变换方法，其特征在于，步骤A6中将矩阵H中处于主对角线以下的子矩阵集合补全为下三角矩阵L具体为：L＝H_d+E，其中，H_d为矩阵H中主对角线以下的子矩阵集合，E为n×k阶的单位矩阵，补全后的矩阵L为n×n个k阶子矩阵。

4.根据权利要求1或3所述的PSF反变换方法，其特征在于，步骤A6中采用并行的高斯消元法计算得到分块向量Q具体包括：

对子矩阵H_i-j与向量P作如下运算，其中j+1≤i≤n，并初始化j＝1；

H_i-j＝H_i-j×P_j (5)

P_i＝P_i-H_i-j (6)

j＝j+1 (7)

其中式(5)和式(6)的计算包括i取值为j+1≤i≤n的所有运算；

重复式(5)～(7)，直至j＝n为止，此时计算得到分块向量Q＝P。

5.根据权利要求1所述的PSF反变换方法，其特征在于，步骤A7中将矩阵H中处于主对角线及主对角线以上的子矩阵集合补全为上三角矩阵U具体为：U＝H_u+Z，其中，H_u为矩阵H中主对角线及主对角线以上的部分，Z为n×k阶的零矩阵，补全后的矩阵U为n×n个k阶子矩阵。

6.根据权利要求1或5所述的PSF反变换方法，其特征在于，步骤A7中采用并行的高斯消元法计算得到物平面图像向量B具体包括：

对子矩阵H_i-j与分块向量Q作如下运算，其中1≤i≤j-1，并初始化j＝n；

Q_j＝Q_j/H_i-j (8)

H_i-j＝H_i-j×Q_j (9)

Q_i＝Q_i-H_i-j (10)

j＝j-1 (11)

其中式(9)和式(10)的计算包括i取值为1≤i≤j-1的所有运算；

重复式(8)～(11)，直至j＝1为止，此时计算得到物平面图像向量B＝Q。

Images