CN103839220A - 基于矩阵分解的图像无缝克隆方法 - Google Patents

Abstract

本发明公开了一种基于矩阵分解的图像无缝克隆方法，该方法使用梯度域方法构建原始线性系统方程，对原始线性系统的系数矩阵进行矩阵分解分解成带状对角矩阵，并将带状对角矩阵分解为若干个小带状对角矩阵，根据各个小带状对角矩阵构建小型系统方程，并通过GPU的不同线程中求解该小型系统方程得到各个像素点的初始像素补偿值，将所有像素点的初始像素补偿值作为原始系统的初始值代入原始线性系统方程直接求解得到最终像素补偿值，进一步得到融合区域中各个像素点的像素值，完成无缝克隆。本发明的图像无缝克隆方法充分利用了GPU的数据并行处理能力，大大提高了运算速率，使实时图像处理成为可能，且内存消耗低，可扩展到广泛的运用程序中。

Description

基于矩阵分解的图像无缝克隆方法

技术领域

本发明涉及计算机图像领域，尤其涉及一种基于矩阵分解的图像无缝克隆方法。

背景技术

近几年，梯度域方法被广泛运用于图像处理领域中，包括本征图像复原、无缝克隆和图像拼接领域。这些算法的求解通常是解一个大型稀疏线性系统：泊松方程。但是，解泊松方程是一个计算密集型和内存密集型任务，不适合实时图像编辑。同时，在处理大规模图像时，如百万像素图，甚至是，亿像素图像，这个问题就显得异常突出。为了解决这个问题，研究人员已经提出了各种各样的方法。这些方法大致可以分为两类：基于退化空间的求解器和基于数值分析的求解器。

对第一类，基于退化空间的求解器，一些算法将全分辨率图像转换到退化空间，从而大量减少内存开销和计算开销。第二类基于数值分析的求解器可以细分为两个子类：直接解决器和迭代解决器。Pardiso，MUMPS，和SuperLU属于直接解决器。而共轭梯度，多重网格法和串流多重网格法是比较常用的迭代解决器。直接法的思想是，经过有限步算术计算，可求得方程的精确解的方法，但由于实践中有四舍五入的存在，这类方法也只能得到近似解，直接法最基本的是高斯消去法。而直接解决器的可扩展性差，不适合求解大型线性系统。迭代法的思想是，采用某种极限过程逐步逼近线性方程的精确解，迭代法需要的计算机存储单元较少、程序设计简单，但是迭代解决器不是很健壮，收敛速度慢。

发明内容

针对现有图像无缝克隆方法技术的不足，本发明提供了一种基于矩阵分解的图像无缝克隆方法。

一种基于矩阵分解的图像无缝克隆方法，包括：

（1）根据源图像中克隆区域的梯度域和目标图像中融合区域的梯度域构建泊松方程，并根据泊松方程得到对应的拉普拉斯方程；

（2）根据所述的拉普拉斯方程构建原始线性系统和原始线性系统方程；

（3）对所述原始线性系统方程的系数矩阵进行分解得到一个带状对角矩阵和残余矩阵，并根据克隆区域的形状将所述的带状对角矩阵划分为若干个小带状对角矩阵；

（4）将各个小带状对角矩阵分配到GPU的各个线程，结合残余矩阵求解得到融合区域中各个像素点相对于克隆区域中相应像素点的初始像素补偿值；

（5）将融合区域中所有像素点相对于克隆区域中相应像素点的初始像素补偿值作为原始线性系统的初始输入，通过直接法迭代得到融合区域中所有像素点相对于克隆区域中相应像素点的最终像素补偿值；

（6）以融合区域中各个像素点相对于克隆区域中相应像素点的最终像素补偿值与源图像中该像素点的像素值的和作为融合区域中相应像素点的像素值，完成图像的无缝克隆。

所述步骤（1）中源图像中的克隆区域与目标图像中的融合区域的大小和形状相同。

无缝克隆法的关键是，寻找一个函数满足目标图像中克隆区域（融合区域）的边界条件，且目标图像中克隆区域的梯度域尽可能与源图像的梯度域相同。

所述步骤（5）中通过直接法迭代得到融合区域中所有像素点相对于克隆区域中相应像素点的最终像素补偿值，即将内部迭代输出的克融合区域中所有像素点相对于克隆区域中相应像素点的最终像素补偿值作为原始线性系统的初始输入。

本发明的图像无缝克隆方法，使用梯度域方法构建泊松方程，并得到原始线性系统方程，对原始线性系统的系数矩阵进行分解，产生大小大致相同的带状对角矩阵和残余矩阵。根据克隆区域的形状将带状对角矩阵划分为小带状对角矩阵并分配给GPU的不同线程中进行计算，充分利用了GPU的数据并行处理能力，大大提高了运算速率，使实时图像处理成为可能，且内存消耗低，可扩展到广泛的运用程序中。且通过将泊松方程产生的线性系统方程的系数矩阵进行分解，对每一个分解得到的小型矩阵采用并行处理迭代操作，并在得到融合区域中各个像素点的初始像素补偿值，将所有像素点的初始像素补偿值作为原始系统的初始值（即输入）代入原始线性系统方程直接求解，在外部迭代器中执行一个直接法计算全解，当外部迭代收敛时，完成求解，可以克服直接法和迭代法的缺点：比直接解决器的可扩展性更好，比传统预处理迭代解决器更健壮。

所述步骤（3）中根据公式：

A=D+R

将原始线性系统的系数矩阵A划分为带状对角矩阵D和残余矩阵R，具体如下：将系数矩阵A的每一行的第一个和最后一个不为零的元素设为零，得到带状对角矩阵D，保留系数矩阵A中每一行的第一个和最后一个不为零的元素，其余值设为零，得到残余矩阵R。

所述步骤（3）通过以下步骤将带状对角矩阵划分为若干个小带状对角矩阵：

（3-1）判断克隆区域是否为规则区域：

若克隆区域为矩形，则克隆区域为规则区域，

否则，克隆区域为不规则区域；

（3-2）根据判断结果进行如下操作：

（a）对规则区域，将带状对角矩阵划分为若干个小带状对角矩阵；

（b）对不规则区域，确定所述带状对角矩阵中每一行的带宽，并确定最大带宽W_max，从第一行开始逐行进行宽度累加直至：

${\mathrm{\Σ}}_{i=s}^t{\mathrm{\ω}}_i\≥W_{\max },$

则进一步判断：

若

则将带状对角矩阵中的第s行至第t行划分为一个小带状对角矩阵，

否则，将带状对角矩阵中的第s行至第t-1行划分为一个小带状对角矩阵，其中，带状对角矩阵中每一行的带宽为第一个不为零的元素与最后一个不为零的元素的列号的差值，ω_i为带状对角矩阵中第i行的带宽，s为每次宽度累加时起始行的行号，每次进行宽度累加均以尚未划分的首行作为起始行，

若分块对角矩阵中尚未划分的所有行的宽度累加和小于W_max时，直接将剩余的所有尚未划分的行作为一个小带状对角矩阵。

为了平衡GPU线程任务，保证各个线程上的计算量大致相同，因此根据克隆区域的形状将带状对角矩阵划分为小带状对角矩阵，对于矩形的克隆区域，各个小分块对角矩阵的大小相同，对于非矩形的克隆区域，采用归并方法将对角块矩阵划分为若干个的大小大致相同的小带状对角矩阵，使不同线程上处理计算的数据量大致相同，从而平衡计算速度。

对于不规则非矩形的克隆区域（即不规则克隆区域），首先从带状对角矩阵的第一行开始划分，即第一次划分时，令起始行号s=1，此后每次划分时，以上一次划分结束后，小带状对角矩阵中尚未参与划分（即没有划分到上一个小带状对角矩阵中的）行号为起始行号。

将不规则区域的带状对角矩阵划分为小带状对角矩阵时，若分块对角矩阵中尚未划分的所有行的宽度累加和小于W_max时，直接将剩余的所有尚未划分的行作为一个小带状对角矩阵。

所述步骤（3-2）中对于规则区域，所述的小带状对角矩阵的个数取决于GPU的内核数。具体如下：

如果GPU是双核则h=2，GPU是4核则h=4，……，GPU是c核则h=c），

各个小带状对角矩阵D_i（i=1,2,...,c）的具体取值方式为：

确定小带状对角矩阵D_i（i=1,2，……，c-1）中的行数为C=[m/c]，其中[m/c]表示对m/c向下取整运算。根据各个小带状对角矩阵中的行数将带状对角矩阵D划分为c个小带状对角矩阵，其中小带状对角矩阵D_i为由带状对角矩阵D中的第(i-1)×CC1到第i×C行中每一行第一个不为零的元素到最后一个不为零的元素构成的矩阵，小带状对角矩阵D_i（i=c）为由带状对角矩阵D中的第i×C+1到第m行中每一行第一个不为零的元素到最后一个不为零的元素构成的矩阵。

所述步骤（4）包括以下步骤：

（4-1）采用丢弃奇数行策略将残余矩阵近似为近似残余矩阵；

（4-2）并将各个小带状对角矩阵分配到GPU的不同线程进行内部迭代得到融合区域中各个像素点相对于克隆区域相应像素点的初始像素补偿值。

所述步骤（4-2）中的任意一个小带状对角矩阵内部迭代如下：

（4-21）将当前小带状对角矩阵扩展为与原始系数矩阵大小相同的矩阵作为当前小带状对角矩阵的近似矩阵；

（4-22）将所述的近似矩阵和近似残余矩阵代入原始线性系统方程得到当前小带状对角矩阵对应的近似线性系统方程：

$(\hat{D}+\hat{R})x=b$

其中，

为任意一个小带状对角矩阵的近似矩阵，

为近似残余矩阵，

x是融合区域中所有像素点相对于克隆区域中相应像素点的像素补偿值组成的向量，

b为原始线性系统方程的常数项向量；

（4-23）根据公式：

$G={\hat{D}}^{-1}\·\hat{R}$

计算得到矩阵G，将矩阵G中所有非零元素对应的像素点作为非零像素点集合，其中，

是

的逆矩阵；

（4-24）根据所述的非零像素点集合构建小型系统方程：

$\hat{G}{\hat{x}}_0=\hat{g}$

其中， $\hat{G}=I_Z+G_Z,$

I_Z是非零像素点集合中所有像素点在克隆区域中对应的像素点的像素值构成的矩阵；

G_Z是非零像素点集合中所有像素点在融合区域中对应的像素点的像素值构成的矩阵；

是融合区域中与非零像素点集合中所有像素点对应的像素点相对于克隆区域中对应的像素点的像素补偿值组成的向量；

$\hat{g}=(g_1,g_2\·\·\·\·\·\·,g_s),$ 其中，g₁，g₂,…g_s分别是 $g={\hat{D}}_1^{-1}b$ 的第1,2，……s个分量的值，s为非零像素点集合中的像素点的个数；

（4-25）使用预处理共轭梯度迭代方法求解小型系统方程得到融合区域中与非零像素点集合中所有像素点对应的像素点相对于克隆区域中对应的像素点的初始像素补偿值。

本发明中根据矩阵G的计算方法，可以看出矩阵G中每一行中至多只有一个非零元素。

所有非零元素对应的像素点作为非零像素点集合，其中，

是

的逆矩阵

所述步骤（4-21）中扩展时令比当前小带状对角矩阵中多出的元素为零。

进行这样的扩展，不会大量增加计算操作，同时通过进行零与非零的比较可以记录当前元素与扩充元素的分割，在进行代码编写时可以采用索引的方式进行记录，并且如果进行随机填充会导致不同小带状对角矩阵的数据冲突导致运算结果不准确。

所述步骤（4-22）中原始线性系统方程的常数项向量b中的各个分量的取值由对应像素点的位置确定，其中第i个分量b_i的取值如下：

若对应的像素点P_i不在源图像克隆区域的边界上，则b_i=0；

否则：

b_i=G(P_i)-I(P_i)，

其中，G(P_i)为像素点P_i在目标图像中的像素值，I(P_i)为像素点P_i在源图像中的像素值。

满足该条件能够保证源图像中的克隆区域在边界上与目标图像的克隆区域无缝的融合。

本发明的图像无缝克隆方法使用梯度域方法构建泊松方程，并得到原始线性系统方程，对原始线性系统的系数矩阵进行分解，产生大小大致相同的分块对角矩阵和残余矩阵，根据克隆区域的形状将分块对角矩阵划分为小分块对角矩阵并分配给GPU的不同线程中进行计算，结合直接法和迭代法，充分利用了GPU的数据并行处理能力，大大提高了运算速率，使实时图像处理成为可能，且内存消耗低，可扩展到广泛的运用程序中。

具体实施方式

下面将结合具体实施方式对本发明进行详细说明。

本实施例中用S代表源图像，T代表目标图像。用P_s表示源图像中的克隆区域，P_t表示目标图像中的融合区域。克隆区域P_s和融合区域P_t的形状和大小相同。

本实施例的基于矩阵分解的图像无缝克隆方法，包括：

（1）根据源图像中克隆区域的梯度域和目标图像中融合区域的梯度域构建泊松方程，并根据泊松方程得到对应的拉普拉斯方程。

构建得到的泊松方程如下：

$f(x,y)=g(x,y),\∀(x,y)\∈{\∂P}_t$

其中，f(x,y)是目标图像中融合区域中坐标为(x,y)的像素点的像素值；

g(x,y)是目标图像中像素点（x，y）（即坐标为(x,y)的像素点）的像素值；

是P_t的边界，

表示的像素点(x,y)为待克隆区域的边界上的像素点；

Δf(x,y)为源图像中像素点（x，y）的像素值的拉普拉斯算子；

是源图像中像素点（x，y）的像素值的梯度，

$\▿I(x,y)=\left[I\right(x+1,y)-I(x,y),I(x,y+1)-I(x,y)],$

其中，I(x+1,y)是像素点（x+1，y）（即坐标为(x+1,y)的像素点）的像素值，

I(x,y)是坐标为(x,y)的像素点的强度，I(x,y+1)是像素点（x+1，y）的像素值；

表示对

求散度；

根据上述泊松方程产生等价Laplace方程：

$f\′(x,y)=g(x,y)-I(x,y),\∀(x,y)$

Δf′(x,y)=0,其它

其中，f′(x,y)是融合区域（目标图像）中像素点(x,y)相对于克隆区域（源图像）中的像素点(x,y)的像素补偿值；

g(x,y)是融合区域中像素点(x,y)的像素值；

I(x，y)是克隆区域中像素点(x,y)的像素值；

??是源图像中的克隆区域；

表示源图像中像素点(x,y)在克隆区域的边界上；

Δf′(x,y)是克隆区域中像素点(x,y)的像素值的拉普拉斯算子。

在离散情况下，源图像中像素点(x,y)的像素值的拉普拉斯算子为：

$\mathrm{\Δf}\′(x,y)=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}f\′\left(n_i\right)-4\×\′(x,y)$

其中，n_i表示像素点(x,y)的第i个四邻域像素点；

f′(n_i)是四邻域像素点n_i在源图像与目标图像中的像素值的差值。

（2）根据拉普拉斯方程构建原始线性系统和原始线性系统方程，原始线性系统方程如下：

Ax=b

其中，A为原始线性系统方程的系数矩阵，

x是源图像中的待克隆区域中每一个像素点对应的强度补偿值所组成的向量，

b为原始线性系统方程的常数项向量，常数项向量b中的各个分量的取值由对应像素点的位置确定，其中第i个分量b_i的取值如下：

若对应的像素点P_i不在源图像的边界上，则b_i=0；

否则：

b_i=G(P_i)-I(P_i)，

其中，G(P_i)为像素点P_i在目标图像中的像素值，I(P_i)为像素点P_i在源图像中的像素值。

本实施例中克隆区域的像素点个数为m，则系数矩阵A是一个m×m的方阵，系数矩阵A中第i（i=1,2,...,m）行与克隆区域中第i个像素点P_i相互对应：

令第i行第i列元素为-4（在离散情况下，拉普拉斯算子计算时f′(x,y)的系数）；

令第i行中像素点P_i的四个四邻域像素点所在的列的元素为1；

令第i行中其余列的取值为0。

（3）对原始线性系统方程的系数矩阵进行分解得到一个带状对角矩阵和残余矩阵，并根据克隆区域的形状将带状对角矩阵划分为若干个小带状对角矩阵。

将系数矩阵A分解带状对角矩阵和残余矩阵的方法如下：将系数矩阵A中的每一行的第一个和最后一个不为零的元素设为零，得到带状对角矩阵D，保留系数矩阵A中每一行的第一个和最后一个不为零的元素，其余值设为零，得到残余矩阵R。带状对角矩阵和残余矩阵满足条件：

A=D+R，

D为带状对角矩阵，且带状对角矩阵D尽量逼近原始线性系统的系数矩阵A。进一步根据待克隆区域的形状将紧密带状对角矩阵D划分为若干个小带状对角矩阵，具体如下：

（3-1）判断待克隆区域是否为规则区域：

若待克隆区域为矩形，则待克隆区域为规则区域，

否则，待克隆区域为不规则区域；

（3-2）根据判断结果进行如下操作：

（a）对规则区域，根据泊松方程的特性，通过上述操作，带状对角矩阵D可表示为：

$D=\left(\begin{array}{cccc}{D}_1& 0& 0& 0\\ 0& D_2& 0& 0\\ 0& 0& \·\·\·& 0\\ 0& 0& 0& D_h\end{array}\right)$

其中，D₁、D₂...D_h则为带状对角矩阵D划分得到的小带状对角矩阵，h为小带状矩阵的个数，小带状对角矩阵的个数h的取值取决于GPU的内核数。具体如下：

如果GPU是c核则h=c，各个小带状对角矩阵D_i（i=1,2,...,h）的具体取值方式为：

确定小带状对角矩阵D_i（i=1,2，……，h-1）中的行数为C=[m/h]，其中[m/h]表示对m/h向下取整运算。根据各个小带状对角矩阵中的行数将带状对角矩阵D划分为c个小带状对角矩阵，其中小带状对角矩阵D_i为由带状对角矩阵D中的第(i-1)×C+1到第i×C行中每一行第一个不为零的元素到最后一个不为零的元素构成的矩阵，小带状对角矩阵D_i（i=h）为由带状对角矩阵D中的第i×C+1到第m行中每一行第一个不为零的元素到最后一个不为零的元素构成的矩阵。

（b）对不规则区域，确定带状对角矩阵中每一行的带宽，并确定最大带宽为W_max，从第一行开始逐行进行宽度累加直至：

${\mathrm{\Σ}}_{i=s}^t{\mathrm{\ω}}_i\≥W_{\max },$

则进一步判断：

若

则将带状对角矩阵中的第s行至第t行划分为一个小带状对角块矩阵，

否则，将带状对角矩阵中的第s行至第t-1行划分为一个小带状对角块矩阵，

其中，带状对角矩阵中每一行的带宽为第一个不为零的元素与最后一个不为零的元素的列号的差值，

ω_i为带状对角矩阵中第i行的带宽，

s为每次宽度累加时起始行的行号，每次进行宽度累加均以尚未划分的首行作为起始行，

若分块对角矩阵中尚未划分的所有行的宽度累加和小于W_max时，直接将剩余的所有尚未划分的行作为一个小带状对角矩阵。

即说明克隆区域为不规则区域时，带状对角矩阵D的每个小带状对角矩阵的会完全不同，因此需要根据（b）中的方法将带状对角矩阵D划分为若干个小带状对角矩阵。

本实施例中克隆区域为不规则区域，因此需要根据（b）将带状对角矩阵D划分为h个小带状对角矩阵。

（4）将各个小带状对角矩阵分配到GPU的各个线程，结合残余矩阵求解得到融合区域中各个像素点相对于克隆区域中相应像素点的初始像素补偿值，包括以下步骤：

（4-1）采用丢弃奇数行策略（即令残余矩阵R中奇数行的元素的取值为零）将残余矩阵R近似为近似残余矩阵

（4-2）并将各个小分块对角矩阵分配到GPU的不同线程进行内部迭代，得到克隆区域中各个像素点相对于融合区域（目标图像）中相应像素点的初始像素补偿值，具体如下：

（4-21）将当前小带状对角矩阵扩展为与原始系数矩阵大小相同的矩阵作为当前小带状对角矩阵的近似矩阵；

（4-22）将所述的近似矩阵和近似残余矩阵代入原始线性系统方程得到当前小带状对角矩阵对应的近似线性系统方程：

$(\hat{D}+\hat{R})x=b$

其中，

为任意一个小带状对角矩阵的近似矩阵，

为近似残余矩阵，

x是融合区域中所有像素点相对于克隆区域中相应像素点的像素补偿值组成的向量，

b为原始线性系统方程的常数项向量；

（4-23）根据公式：

$G={\hat{D}}^{-1}\·\hat{R}$

计算得到矩阵G，将矩阵G中所有非零元素对应的像素点作为非零像素点集合，其中，是

的逆矩阵；

（4-24）根据所述的非零像素点集合构建小型系统方程：

$\hat{G}{\hat{x}}_0=\hat{g}$

其中， $\hat{G}=I_Z+G_Z,$

I_Z是非零像素点集合中所有像素点在克隆区域中对应的像素点的像素值构成的矩阵；

G_Z是非零像素点集合中所有像素点在融合区域中对应的像素点的像素值构成的矩阵；

是融合区域中与非零像素点集合中所有像素点对应的像素点相对于克隆区域中对应的像素点的像素补偿值组成的向量；

$\hat{g}=(g_1,g_2\·\·\·\·\·\·,g_s),$ 其中，g₁,g₂,…g_s分别是 $g={\hat{D}}_1^{-1}b$ 的第1,2，……s个分量的值，s为非零像素点集合中的像素点的个数；

（4-25）使用预处理共轭梯度迭代方法求解小型系统方程

得到融合区域中与非零像素点集合中所有像素点对应的像素点相对于克隆区域中对应的像素点的初始像素补偿值。

对于小带状对角矩阵D₁对应的所有像素点中不属于非零像素点集合中的像素点（即在

中不存在对应元素的像素点），以非零像素点集合中与该像素点最邻近（距离最小的像素点）的像素点的初始像素补偿值作为该像素点的初始像素补偿值，即可得到小带状对角矩阵D₁中包含的所有像素点的初始猜测补偿值并形成相应的初始像素补偿值向量x₁。

GPU的线程分别对小带状对角矩阵D₁、小带状对角矩阵D₂...小带状对角矩阵D_h进行（4-22）到（4-25）的处理，得到各个小带状对角矩阵对应的像素点的初始像素补偿值向量x₁，x₂，…x_h。

（5）将待克隆区域中所有像素点的初始强度补偿值作为原始线性系统的初始输入进行外部迭代，当算法收敛时即完成泊松方程的求解得到各个像素点的最终强度补偿值。

外部迭代采用直接求解Ax=b的思想，x为向量，其初始值为步骤（4）中内部迭代时GPU的c个线程分别处理的小带状对角矩阵D₁、D₂...D_c，内部迭代得到的初始像素补偿值向量x₁，x₂，…，x_h组合得到（初始像素补偿值向量x₁，x₂，…，x_h组合覆盖带状对角矩阵D中所有行，因此初始像素补偿值向量x₁，x₂，…，x_h组合可以得到融合区域中所有像素点的初始像素补偿值）。具体方法如下：将x₁，x₂，…，x_h的每个分量对应的像素点的初始强度补偿值作为向量x中该像素点对应的分量的取值，如：x₁中的某个元素对应融合区域中的第i个像素点，则说明将元素为融合区域中的第i个像素点的初始强度补偿值，因此，将该分量作为向量x中第i行的分量。

（6）以各个像素点的最终像素补偿值与源图像中该像素点的像素值的和作为目标图像中相应像素点的像素值，即完成图像无缝克隆。

根据f（x，y）-I（x，y）+f′（x，y），计算目标图像融合区域的强度，其中，f（x，y）表示目标图像中融合区域像素点（x，y）最终的像素值；

I（x，y）表示目标图像融合区域中像素点（x，y)在源图像中对应的像素点的像素值；

f′（x，y）表示目标图像融合区域中像素点（x，y）的最终像素补偿值，通过外部求解得到（为外部迭代解得的向量x中与像素点（x，y)对应的元素）。

以上公开的仅为本发明的具体实施例，但是本发明的保护范围并不局限于此，任何熟悉本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围，都应涵盖在本实用发明的保护范围之内。

Claims (7)

1.一种基于矩阵分解的图像无缝克隆方法，其特征在于，包括：

（1）根据源图像中克隆区域的梯度域和目标图像中融合区域的梯度域构建泊松方程，并根据泊松方程得到对应的拉普拉斯方程；

（2）根据所述的拉普拉斯方程构建原始线性系统和原始线性系统方程；

（3）对所述原始线性系统方程的系数矩阵进行分解得到一个带状对角矩阵和残余矩阵，并根据克隆区域的形状将所述的带状对角矩阵划分为若干个小带状对角矩阵；

（4）将各个小带状对角矩阵分配到GPU的各个线程，结合残余矩阵求解得到融合区域中各个像素点相对于克隆区域中相应像素点的初始像素补偿值；

（5）将融合区域中所有像素点相对于克隆区域中相应像素点的初始像素补偿值作为原始线性系统的初始输入，通过直接法迭代得到融合区域中所有像素点相对于克隆区域中相应像素点的最终像素补偿值；

（6）以融合区域中各个像素点相对于克隆区域中相应像素点的最终像素补偿值与源图像中该像素点的像素值的和作为融合区域中相应像素点的像素值，完成图像的无缝克隆。

2.如权利要求1所述的基于矩阵分解的图像无缝克隆方法，其特征在于，所述步骤（3）通过以下步骤将带状对角矩阵划分为若干个小带状对角矩阵：

（3-1）判断克隆区域是否为规则区域：

若克隆区域为矩形，则克隆区域为规则区域，

否则，克隆区域为不规则区域；

（3-2）根据判断结果进行如下操作：

（a）对规则区域，将带状对角矩阵划分为若干个小带状对角矩阵；

（b）对不规则区域，确定所述带状对角矩阵中每一行的带宽，并确定最大带宽W_max，从第一行开始逐行进行宽度累加直至：

${\mathrm{\Σ}}_{i=s}^t{\mathrm{\ω}}_i\≥W_{\max },$

则进一步判断：

若

则将带状对角矩阵中的第s行至第t行划分为一个小带状对角矩阵，

否则，将带状对角矩阵中的第s行至第t-1行划分为一个小带状对角矩阵，其中，带状对角矩阵中每一行的带宽为第一个不为零的元素与最后一个不为零的元素的列号的差值，ω_i为带状对角矩阵中第i行的带宽，s为每次宽度累加时起始行的行号，每次进行宽度累加均以尚未划分的首行作为起始行，

若分块对角矩阵中尚未划分的所有行的宽度累加和小于W_max时，直接将剩余的所有尚未划分的行作为一个小带状对角矩阵。

3.如权利要求2所述的基于矩阵分解的图像无缝克隆方法，其特征在于，所述步骤（3-2）中对于规则区域，所述的小带状对角矩阵的个数等于GPU的内核数。

4.如权利要求3所述的基于矩阵分解的图像无缝克隆方法，其特征在于，所述步骤（4）包括以下步骤：

（4-1）采用丢弃奇数行策略将残余矩阵近似为近似残余矩阵；

（4-2）并将各个小带状对角矩阵分配到GPU的不同线程进行内部迭代得到融合区域中各个像素点相对于克隆区域相应像素点的初始像素补偿值。

5.如权利要求4所述的基于矩阵分解的图像无缝克隆方法，其特征在于，所述步骤（4-2）中的任意一个小带状对角矩阵内部迭代如下：

（4-21）将当前小带状对角矩阵扩展为与原始系数矩阵大小相同的矩阵作为当前小带状对角矩阵的近似矩阵；

（4-22）将所述的近似矩阵和近似残余矩阵代入原始线性系统方程得到当前小带状对角矩阵对应的近似线性系统方程：

$(\hat{D}+\hat{R})x=b$

其中，

为任意一个小带状对角矩阵的近似矩阵，

为近似残余矩阵，

x是融合区域中所有像素点相对于克隆区域中相应像素点的像素补偿值组成的向量，

b为原始线性系统方程的常数项向量；

（4-23）根据公式：

$G={\hat{D}}^{-1}\·\hat{R}$

计算得到矩阵G，将矩阵G中所有非零元素对应的像素点作为非零像素点集合，其中，是

的逆矩阵；

（4-24）根据所述的非零像素点集合构建小型系统方程：

$\hat{G}{\hat{x}}_0=\hat{g}$

其中， $\hat{G}=I_Z+G_Z,$

I_Z是非零像素点集合中所有像素点在克隆区域中对应的像素点的像素值构成的矩阵；

G_Z是非零像素点集合中所有像素点在融合区域中对应的像素点的像素值构成的矩阵；

是融合区域中与非零像素点集合中所有像素点对应的像素点相对于克隆区域中对应的像素点的像素补偿值组成的向量；

$\hat{g}=(g_1,g_2\·\·\·\·\·\·,g_s),$ 其中，g1,g2,…gs分别是 $g={\hat{D}}_1^{-1}b$ 的第1,2，……s个分量的值，s为非零像素点集合中的像素点的个数；

（4-25）使用预处理共轭梯度迭代方法求解小型系统方程

得到融合区域中与非零像素点集合中所有像素点对应的像素点相对于克隆区域中对应的像素点的初始像素补偿值。

6.如权利要求5所述的基于矩阵分解的图像无缝克隆方法，其特征在于，所述步骤（4-21）中扩展时令比当前小带状对角矩阵中多出的元素为零。

7.如权利要求6所述的基于矩阵分解的图像无缝克隆方法，其特征在于，所述步骤（4-22）中原始线性系统方程的常数项向量b中的各个分量的取值由对应像素点的位置确定，其中第i个分量b_i的取值如下：

若对应的像素点P_i不在源图像克隆区域的边界上，则b_i=0；

否则：

b_i=G(P₁)-I(P_i)，

其中，G(P_i)为像素点P_i在目标图像中的像素值，I(P_i)为像素点P_i在源图像中的像素值。