Replicated边界条件下图像模糊矩阵与图像矢量乘积的一种替代计算方法与流程

本发明涉及图像处理，特别地，本发明涉及到Replicated边界条件下的线性图像恢复中，对图像模糊矩阵及其转置与图像矢量乘积的一种替代计算。

背景技术：

图像恢复是对线性位移不变引起的下述退化问题进行求解

f＝Hx+η (1)

其中：图像模糊矩阵H∈R^mn×mn由点扩展函数(Point Spread Function,PSF)h∈R^p×q的元素构成，其结构唯一地由所采用的边界条件决定；f∈R^mn×1由退化图像F∈R^m×n按行优先排列得到的矢量；η为加性噪声的重排矢量。

在对(1)的求解中，主要的计算任务,是计算H与f乘积(乘积一)，即g1＝Hf，以及H的伴随矩阵H^*(在实数域中，即为H的转置矩阵H^T)与f的乘积(乘积二)，即g2＝H^Tf。

一般地，H和H^T皆为严重稀疏的超大型矩阵，使得乘积一和乘积二难于直接计算。对此，针对一些边界条件下的乘积一和乘积二，目前已有两种替代计算方式：一是基于Kronecker乘积的替代计算方式，一种是基于对图像模糊矩阵进行对角化的替代计算方式。在确定性(Determinate)边界条件下，乘积一和乘积二能否被替代计算，以及能的话，被上述哪一种方式替代计算更为有效，这两个问题都取决于采用的是何种边界条件，以及点扩展函数是否对称。

目前，确定性的图像边界条件有Zero边界条件、Periodic边界条件、Neumann边界条件、Anti-Reflective边界条件、Mean边界条件和Replicated边界条件，在这些条件下，图像模糊矩阵具有与退化图像不相关的确定性结构，且结构只唯一地由边界条件决定，这使得我们可以去分析这些结构，并将其应用到乘积一和乘积的计算,和对(1)式的预置器的设计中。

Replicated边界条件是一种最新的边界条件，由我国学者Xu zhou等于2014年在文献”A boundary condition based deconvolution framework for image deblurring”中提出，Xu zhou等认为该边界条件考虑了图像像素符合“重尾”分布的情况(即图像中大多像素和它的邻域像素值相近)，因此是一种较为符合图像实际的边界条件，可应用在对具有相应特征的图像的恢复和处理中。

Xu zhou等在上述同一文献中还提出了基于Replicated边界的图像去模糊框架：1)对H进行乘法分解，即将H转化为非方的有偏模糊矩阵K∈R^{(m+2p1)(n+2q1)×mn}与非方的边界条件矩阵P∈R^{mn×(m+2p1)(n+2q1)}的乘积，即H＝KP，从而将图像退化的问题f＝Hx+ζ转化为f＝KPx+ζ；2)采用共轭梯度算法求解线性方程P^*K^*y＝P^*K^*KPx，来恢复未知图像x，并在对应的迭代求解中，利用两次快速傅立叶变换(Fast Fourier Transformation，FFT)和一次快速傅立叶逆变换来计算乘积一和乘积二的矩阵形式；3)确定了Replicated边界条件下P和P^*的结构，使得Px和P^*x可被计算，从而提出了Replicated边界条件下乘积一和乘积二的一种替代计算方法，也是该边界条件下唯一的一种替代计算方法。Xu Zhou等对乘积一和乘积二的替代计算方法算法复杂度为O(m×n×Nnz)，其中Nnz为点扩展函数中非零元素的个数。在其他边界条件下，该框架是否可行，则取决于边界矩阵P和P^*的结构是否已确定。

技术实现要素：

本发明所要解决的技术问题是：设计Replicated边界条件下乘积一和乘积二的一种新的间接计算方法，该方法不对模糊矩阵进行乘法分解，而是进行加法分解，使对乘积一和乘积二的计算精度更高，并且在点扩展函数为常规大小时，时间耗费与Xu Zhou的方法相当。

本发明针对上述问题，给出一种Replicated边界条件下乘积一和乘积二的间接计算方法，思路是：在Replicated边界条件下，根据图像模糊矩阵的结构，先将乘积一和乘积二，分别分解为多个矩阵与矢量乘积的和，保证所产生的分解矩阵带可利用的分块结构；利用各分块矩阵的结构，构造其卷积核；用卷积替代地计算各分解矩阵与图像矢量的乘积的矩阵形式；综合各计算结果，得到乘积一和乘积二。

为支持上述思路，本发明中：实施例1给出Replicated边界条件下的乘积一和乘积二的替代计算方法的具体步骤；实施例2将Replicated边界条件下的乘积一和乘积二，分别转化为多个模糊矩阵与图像矢量的乘积之和；实施例3针对乘积一与乘积二的分解式中的各分解矩阵，构造对应的点扩展；实施例4基于构造的点扩展函数,采用卷积，计算乘积一和乘积二各个分解矩阵与图像矢量乘积的矩阵形式；实施例5通过实验结果，对本发明给出的卷积替代方法的有效性进行了说明。

本发明所给出的计算方法，提供了Replicated边界条件下，对应大型乘积一和乘积二的替代计算，可集成到Replicated边界条件下的图像滤波和图像恢复中，具有应用价值。本发明所给出的乘积一和乘积二的计算精度较高。本算法的时间复杂度，主要来自于对零边界条件下的乘积一和乘积二的卷积运算复杂度，即O(m×n×p×q)，因此当点扩展函数为普通情形，即普通大小且元素值随机时，本方法与Xu Zhou的方法时间耗费相当，而当点扩展的大小较大而非零元素个数较少时，本专利方法会比Xu Zhou的方法慢。

附图说明

图1示出的本发明的计算流程；

图2示出的是实验中所用的灰度图像Saturno图像，大小为512×512；

图3示出的是实验中所用的灰度图像Cameraman图像,大小为256×256；

图4示出是实验中图像Saturno所用的点扩展函数Saturnopsf,

图5、图6和图7示出的是实验中图像Cameraman所用的三个随机生成的点扩展函数PSFn1,PSFn2和PSFn3,大小分别为53×53,7×7，17×13和31×37。

具体实施方式

在以下描述中，出于对专利说明的目的，结合附图，对发明内容中所述实施例的具体内容进行阐述。为了表达清晰和方便于编程实现，一些计算公式中应用了Matlab的函数形式。第1实施例

如图1所示，本发明给出Replicated边界条件下乘积一g1＝Hf和乘积二g2＝H^Tf的一种替代计算方法，包括以下步骤：

(1)输入F和h，对F按行优先重排成图像矢量f；

(2)对g1和g2进行分解；

(3)对g1和g2分解公式中各分解部分涉及的模糊矩阵，构造其点扩展函数，据此，采用点扩展函数与图像或对应图像边界间卷积，计算各分解部分的矩阵形式；

(4)计算g1和g2的对应图像；

(5)对g1和g2的对应图像，按行优先重排成矢量，得到g1和g2。

第2实施例

本发明对g1的分解为g1＝Hf＝(HTT+HTL+HLT+HLL)f,对g2的分解为g2＝H^Tf＝(H^TTT+H^TTL+H^TLT+H^TLL)f，两个分解公式中：各矩阵中，上标T表示对矩阵的转置，第1个下标代表当前矩阵的块间结构，第2个下标代表当前矩阵的块内结构，T和L分别具体代表Toeplitz矩阵和L-type矩阵，这里不给出T与L矩阵的具体结构，这不会影响后续对计算方法的描述。

第3实施例

本发明给出乘积一分解公式g1＝Hf＝(HTT+HTL+HLT+HLL)f和乘积二分解公式g2＝H^Tf＝(H^TTT+H^TTL+H^TLT+H^TLL)f中，各矩阵的点扩展函数的构造方法。具体为：

记fliplr(.)表示对矩阵列方向即左右翻转，flipud(.)表示对矩阵行方向即上下翻转，conv2(.,.)表示二维卷积，i:j表示从i到j的所有下标，V(:,:)中的:标引矩阵V的所有行与列，V(end,end)中的end标引矩阵V的最后行与最后列；此外，记(注：我们将点扩展函数的行数p和列数q同假定为奇数，即p＝2p1+1，q＝2q1+1，即使真实点扩展函数的大小不同为奇数，也可以通过在点扩展函数的行或列尾，追加一个且仅一个零行或/和追加一个且仅一个零列，使之同成为奇数，相应修改后的点扩展函数，不会影响本专利的计算结果)，则：

(1)HTT的点扩展函数为h1＝h；

(2)构造HTL的点扩展函数，共2个核，分别对应每个L型子块首尾共2列中的非零元素，具体方法为和其中i＝1,...,q1 j＝1,...,p；

(3)构造HLT的点扩展函数，共2个核，分别对应HLT的首尾共两个非零列块，具体为：为和其中1≤i≤p1；

(4)构造HLL的点扩展函数，共4个核，分别对应HLL首尾共两个非零列块中的首尾非零列，(一共4个非零列)，相应计算式与HTL的核有关，具体为和

(5)H^TTT的点扩展函数为h2＝fliplr(flipud(h))；

(6)构造H^TTL的点扩展函数，共2个核，分别对应其的方法为h2_TL_1＝fliplr((h1_TL_1)^T)和h2_TL_2＝fliplr((h1_TL_2)^T)；

(7)构造H^TLT的点扩展函数，共2个核的方法为h2_LT_1＝fliplr(flipud(h1_LT_1))和h2_LT_2＝fliplr(flipud(h1_LT_2))；

(8)构造H^TLL的点扩展函数，共4个核的方法为h2_LL_1＝fliplr(flipud((h1_LL_1)^T))，h2_LL_2＝fliplr(flipud((h1_LL_2)^T))，h2_LL_3＝fliplr(flipud((h1_LL_3)^T))和h2_LL_4＝fliplr(flipud((h1_LL_4)^T))。

第4实施例

本发明基于卷积，以矩阵形式，给出乘积一分解公式g1＝Hf＝(HTT+HTL+HLT+HLL)f和乘积二g2＝H^Tf＝(H^TTT+H^TTL+H^TLT+H^TLL)f中，各矩阵与矢量乘积的计算方法，具体为：

(1)HTTf对应的矩阵形式G1_TT，大小为m×n，其计算方法具体为：G1_TT＝conv2(F,h,'same')；

(2)HTLf的矩阵形式G1_TL，大小为m×n，计算方法为：计算r1_TL_1＝conv2(flipud((h1_TL_1)^T),F(:,1))和r1_TL_2＝conv2(flipud((h1_TL_2)^T),F(:,n))，则G1_TL(:,1:q1)＝r1_TL_1(p1+1:end-p1,:),G1_TL(:,end-q1+1:end)＝r1_TL_2(p1+1:end-p1,:),G1_TL的其它区域中像素值为0；

(3)HLTf的矩阵形式G1_LT，大小为m×n，计算方法为：计算r1_LT_1＝conv2(h1_LT_1,F(1,:))和r1_LT_2＝conv2(h1_LT_2,F(m,:))，则G1_LT(1:p1,:)＝r1_LT_1(:,q1+1:end-q1)，G1_LT(end-p1+1:end,:)＝r1_LT_2(:,q1+1:end-q1)，G1_LT的其它区域中像素值为0；

(4)HLLf的矩阵形式G1_LL，大小为m×n，计算方法为：计算r1_LL_1＝conv2((h1_LL_1)^T,F(1,1))，r1_LL_2＝conv2((h1_LL_2)^T,F(1,n))，r1_LL_3＝conv2((h1_LL_3)^T,F(m,1))和r1_LL_4＝conv2((h1_LL_4)^T,F(m,n))，则G1_LL(1:p1,1:q1)＝r1_LL_1，G1_LL(1:p1,end-q1+1:end)＝r1_LL_2，G1_LL(end-p1+1:end,1:q1)＝r1_LL_3，G1_LL(end-p1+1:end,end-q1+1:end)＝r1_LL_4，G1_LL的其它区域中像素值为0；

(5)H^TTTf的矩阵形式G2_TT，大小为m×n，计算方法为,G2_TT＝conv2(F,fliplr(flipud(h)),'same')；

(6)H^TTLf的矩阵形式G2_TL，大小为m×n，计算方法为：计算r2_TL_1＝conv2(h2_TL_1,F(:,1:q1))，r2_TL_2＝conv2(h2_TL_2,F(:,end-q1+1:end))，则G2_TL(:,1)＝r2_TL_1(p1+1:end-p1,q1)，G2_TL(:,n)＝r2_TL_2(p1+1:end-p1,q1)，G2_TL的其他区域内像素值为0；

(7)H^TLTf的矩阵形式G2_LT，大小为m×n，计算方法为：计算r2_LT_1＝conv2(h2_LT_1,F(1:p1,:))，r2_LT_2＝conv2(h2_LT_2,F(end-p1+1:end,:))，则G2_LT(1,:)＝r2_LT_1(p1,q1+1:end-q1)，G2_LT(m,:)＝r2_LT_2(p1,q1+1:end-q1)，G2_LT的其他区域内像素值为0；

(8)H^TLLf的矩阵形式G2_LL，大小为m×n，计算方法为：计算r2_LL_1＝conv2(h2_LL_1,F(1:p1,1:q1))，r2_LL_2＝conv2(h2_LL_2,F(1:p1,end-q1+1:end))，r2_LL_3＝conv2(h2_LL_3,F(end-p1+1:end,1:q1))，r2_LL_4＝conv2(h2_LL_4,F(end-p1+1:end,end-q1+1:end))，则G2_LL(1,1)＝r2_LL_1(p1,q1)，G2_LL(1,n)＝r2_LL_2(p1,q1)，G2_LL(m,1)＝r2_LL_3(p1,q1)，G2_LL(m,n)＝r2_LL_4(p1,q1)，G2_LL的其它区域中像素值为0。

第5实施例

本发明在RPBC边界条件下，给出乘积一的矩阵结果为G1＝G1_TT+G1_TL+G1_LT+G1_LL，乘积二的矩阵结果为G2＝G2_TT+G2_TL+G2_LT+G2_LL，两式右边各个矩阵的计算方法，已由实施例4给出。

第6实施例

本发明以一个实验，从计算误差和计算时间两个方面，来验证上述实施例所给出方法的有效性；实验采用Matlab2014a，在CPU为2.0GHz,内存为4G的ThinkPad笔记本上完成。实验材料为：实验图像F包括大小为512×512的Saturno图像，和大小为256×256的Cammerman图像，分别见图2和图3；对实验图像1，点扩展函数h采用的是大小为53×53的Saturnopsf，其能量主要集中在的h(26:30,25:29)范围，见图4；对实验图像2，点扩展函数h采用的是3个随机生成的矩阵，大小分别为7×7，17×13和31×37，见图5，图6和图7。

下面表1和表2中，分别对比了在计算乘积一和乘积二的过程中，本发明方法和Xu zhou等的方法分别得到的均方根误差(Root mean square error,Rmse)和所用时间(单位为秒)。Rmse的计算方法为其中为采用本专利方法或Xu Zhou方法对乘积一或乘积二的替代计算结果，g为真实的乘积一或乘积二。

表1乘积一的实验结果

表2乘积二的实验结果

从上面表1和表2可以看出，在计算精度上，本发明方法计算得到的乘积一和乘积二的Rmse更小，精度更高。

上面表1和表2的前三行显示，当点扩展函数为前三个所代表的较常规的情形时，本方法与Xu Zhou的方法在时间耗费上相当；两个表格的最后一行显示，当点扩展函数变大，尤其是点扩展函数中的非零元素较少时，本方法的时间耗费比Xu Zhou等的方法要大，主要是耗费在了对G1_TT和G2_TT的卷积计算上。