基于压缩感知的高动态范围图像去伪影融合方法与流程

本发明涉及一种将低动态范围图像合成高动态范围图像的方法，具体来说，它涉及一种利用k-svd字典学习、压缩感知和去伪影融合算法完成对多曝光的低动态范围图像序列的压缩感知和融合，最后生成没有伪影的高动态范围图像的方法。

背景技术：

高动态范围成像已经开始成为一种商业产品，比如hdr相机、hdr电视机等。在相同真实的场景中，大多数成像传感器的有限动态范围，往往无法捕捉到场景完整动态范围的亮度，然而使用一种相对简单廉价的方式可以解决这个限制，就是捕获若干对同一场景不同曝光时间的图像然后把它们融合成一张记录场景亮度的高动态范围图像，因此有效的扩展图像动态范围。然而传统的图像获取方法需先采集含冗余的大量数据，再进行压缩处理以提取有用信息，效率低并且需要较大存储空间。而压缩感知可避免对冗余数据的采样，在远低于奈奎斯特采样率的条件下，随机采样获取离散信号，并通过重构方法高概率的重建原信号。此外由于相机的抖动或者场景中物体的移动，往往会导致融合后的图像出现伪影或者模糊的现象。

技术实现要素：

针对以上的不足，本发明提供了一种基于压缩感知的高动态范围图像去伪影融合方法。本发明有效降低图像获取端的采样率和计算复杂度，节约存储空间，并解决现有融合后的图像出现伪影或者模糊的现象。

为了解决上述技术问题，本发明的技术方案为：

本发明基于压缩感知的高动态范围图像去伪影融合方法，包括有如下步骤：

1)对输入多曝光图像序列进行分块压缩采样；

2)对压缩采样后的图像块进行ldr图像序列的重构；

3)对经过压缩感知后的多曝光图像序列进行高动态范围图像的去伪影融合。

所述步骤1)的具体实现过程为：

11)以当前输入的ldr图像为参考，进行k-svd字典学习，得到用于稀疏表示的超完备字典ψ；

12)对分块后的图像信号矩阵f进行稀疏表示，得到f＝ψθ，设计观测矩阵φ，则压缩采样的过程重写为：

y＝φf＝φψθ＝θθ\*mergeformat(1)。

所述步骤11)还包括如下内容：

111)判断当前的ldr图像是否为灰度图像，若不是，则转换为灰度图像，并进一步转换为double类型的图像；

112)选择所需的图像分块大小，将111)中获取的灰度图像进行分块，并重排图像块为矩阵列，形成字典学习的数据集矩阵，也即待稀疏表示的信号矩阵y，进一步将矩阵y表示为其中n为矩阵y的列数；

113)k-svd算法参数初始化，包括：

numatom：需要训练的字典元素个数；

numiteration：迭代次数

errorflag：等于0则表示每个信号的稀疏系数个数固定，需要配置参数l；否则表示信号的稀疏系数不固定，需要配置参数errorgoal；

errorgoal：最大允许表示误差

preservedcatom：等于1则表示字典的第一个原子为常量；

initializationmethod：若为“dataelements”则用信号本身作为初始字典，若为“givenmatrix”则用给定的矩阵作为初始字典，并对初始字典进行归一化处理，初始字典表示为ψ⁽⁰⁾∈r^n×k；

初始化算法迭代次数j＝1；

114)稀疏编码阶段：使用正交匹配追踪算法omp来计算y中每一个列向量在当前字典上的稀疏表示所对应的系数向量xi，通过求解下式中的目标函数：

其中为由系数向量所构成的矩阵，i表示列向量的索引，||||f为frobenius范数，||||0为l0范数，t0为稀疏表示稀疏中非零分量的数目的上限，即系数向量中的最大差异度；

115)字典更新阶段：初始字典往往不是最优的，满足稀疏性的稀疏矩阵表示的数据和原数据可能会有比较大的误差，因此需在满足稀疏的条件下逐行逐列进行字典的更新优化，减小整体误差，逼近最优字典，对于ψ^(j-1)中的每一列k＝1,2,...k，定义使用到字典中第k个原子ψk的所有信号集合{yi}的索引所构成的集合为：

其中为系数矩阵x中原子ψk所对应的第k行，不同于x的第k列xk的转置，则目标函数式(2)重写为：

上式中，字典与系数矩阵的乘积ψx被分解为k个秩为1的矩阵的和的形式，当前更新的矩阵列为第k列，则其余k-1个项固定，剥离第k个原子的贡献，得到的矩阵ek为去除原子ψk的成分其余原子所造成的误差；

116)由于中0的存在，用svd得到的更新向量中的非零元素的数量和位置和原中的都不相同，即产生了“发散”，因此需要去除中所有0元素，定义矩阵ωk为n×ωk的矩阵，它在(ωk(i),i)处的值为1，其他位置上的值都为0，则我们定义y、ek去除零项后的收缩结果为的长度为|ωk|，是用到原子ψk的信号的集合，为剥离原子ψk的误差，两者都为n×|ωk|大小的矩阵，则式(4)进一步转换为：

117)对进行svd分解，即矩阵u的第一列表示为即为列ψk的更新项，更新的为矩阵v的第一列乘上δ(1,1)，j＝j+1进入下一次迭代直至满足迭代条件为止，迭代条件根据设定，选择为选定的迭代次数或者是误差值，迭代完成后，即得到最优的超完备字典ψ。

所述步骤12)还包括：

121)观测矩阵φ的行数为选择的图像块大小的平方乘上采样率，用m表示，列数定义为图像块的大小的平方，用n表示，此外传感矩阵θ＝φψ应满足rip性质(有限等距性质)：保证观测矩阵不会把两个不同的稀疏信号映射到同一个集合中，保证原空间到稀疏空间的一一映射关系，定义观测矩阵φ的rip常数δf为满足下式的最小值：

其中f即为分块后的图像信号矩阵，也即稀疏表示后的信号，假设为f稀疏的，即含有f个非零元素，若δf≤1，测成矩阵满则f阶rip条件。

所述步骤2)的具体实现过程为：

21)初始化：令标签索引集残差向量r0＝y，y即为步骤12)中所获得的观测向量，迭代次数t＝1；

22)辨识：在传感矩阵θ中求出与残差向量rt-1最相关的所对应的列向量λt，即内积值最大时所对应的索引(列序号)：

其中n为传感矩阵的列数，θj表示传感矩阵的第j列；

23)更新：更新标签索引集λt＝λt-1∪{λt}，将找到的原子所对应的列向量添加到集合其中，t＞0，λt表示t次迭代的索引，θt表示按索引λt选出的矩阵θ的列集合(大小为m×t的矩阵)；

24)估计：求y＝θθ的最小二乘解：

25)更新残差向量

26)令t＝t+1；并不断重复步骤22)至步骤26)，若满足某个迭代停止条件，则停止迭代并进入第b7)步；

27)求得输出系数向量

所述步骤26)还包括：

在上述算法的步骤26)中，其迭代停止条件有以下常见的3种情况：

261)当运行到t＞s时，迭代停止，其中，s表示固定的迭代步数，如m/4；

262)残差向量的能量值小于某个预先给定的常数ε，例如取1e-6；

||rt||2≤ε\*mergeformat(11)

263)当传感矩阵θ的任何一列都没有残差向量rt的明显能量时：

||θ^trt||∞≤ε\*mergeformat(12)。

所述步骤3)的具体实现过程为：

31)定义输入的经过压缩感知后的图像中任意一点处的像素值为b，对灰度值b归一化处理之后该点的灰度值为i：

i＝b/255\*mergeformat(13)

32)现假设源图像为s，参考图像r，源图像和参考图像合成的图像为l，patchmatch算法就是一个以参考图像为模板，配准源图像生成图像l的过程；由于patchmatch算法是处理一对图像，现假设输入图像为i1...in，以n＝5为例，首先令i3为参考图像r，则i3和i4作为其源图像s，然后令i2和i4作为参考图像r，i1和i5分别作为i2和i4对应的源图像s；

33)现定义patchmatch算法合成图像l的二次函数：

其中τ为灰度映射函数，ω为图像r和图像s的图像域，i为图像域上的任意一个像素点，n(i)为以i中心p×p的邻域，其中p为邻域的大小，故j是邻域n(i)上的像素点，r(i)是图像r上的第i个像素点，s(i+u(j))是图像s上的第(i+u(j))个像素点，其中u(j)表示从图像l上的像素点j映射到图像s的偏移量，α为一个归一化的因子其中wτ和wu为一对权重函数，wτ(i)表示图像r中像素点i映射到图像l的比重，wu(j)表示图像s中像素点j加上偏移量u(j)映射到图像l的比重；

34)灰度映射函数τ定义如下：

其中灰度映射函数的导数τ′≥0，τ()∈[0,1]，i为图像域ω上的像素点，故l(i)表示为图像l上的第i个像素点，使用迭代重加权最小二乘法算法求解灰度映射函数，则将灰度映射函数τ的目标函数重写为：

其中求解目标函数τ过程中，τ和权重因子ω更新为：

其中n表示迭代次数，δ为一个很小的正常数，令δ＝10^-10；

35)权重函数wτ定义如下：

当图像r中图像域上亮度太暗或太亮(即图像域的像素点的灰度值小于3/255，或者大于252/255)时，像素点将会被clipped，否则，就不clipped；

36)权重函数wμ定义如下：

其中d(·,·)表示输入参数之间的空间距离；υ1,υ2为两个归一化参数，分别取对应空间距离的75百分位数；

37)对于参数x和y，d(x,y)＝||x-y||²，而对于表示为：在图像r和图像s的图像域上的任意一个像素点i，取以i为中心，大小为p×p的邻域，得到图像块和然后图像块经过灰度映射函数得到而图像块相对于i平移了u(i)得到最后求两者的空间距离；同理得其中τ-¹()灰度映射函数的逆函数；

38)通过上面定义的函数知，patchmatch算法实际上就是求解二次函数的过程，输入图像r和图像s，并分别对两幅图像进行向下采样，分别得到图像r和图像s的金字塔图像集，从金字塔顶端向下迭代，求得在对应每层金字塔图像下合成的图像l和灰度映射函数τ，将此结果作为下一次的迭代的初始值，当迭代完成后即得到最终的配准图像l，依照此方法，即得到输入图像i1in配准后的图像l1ln；

39)利用伽马曲线对已经配准好的图像集进行辐射校准，消除潜在的相机的移动而产生的噪音.定义伽马(gamma)函数为：

gamma＝cr^γ\*mergeformat(21)

其中c和γ为常数，c取为1，γ取为2.2；这里将经patchmatch配准过的图像集用gamma曲线进行辐射校准；

310)对配准后的图像集使用秩最小化算法得到批量的对齐图像，首先列向量化所有输入的矩阵得到矩阵d，并初始化低秩矩阵a和噪声矩阵e，根据增广拉格朗日乘子法经过内迭代和外迭代得到最佳低秩矩阵a，则噪声矩阵e＝d-a，最后对得到的低秩矩阵和噪声矩阵进行调整m×n大小的图像，即得到输入图像l对应的低秩图像和噪声图像；

311)输入对齐后的图像集a，将图像a合成目标hdr图像：

其中nimg表示为输入图像的数量，x∈{r,g,b}，r，g，b为彩色图像的三个通道；a(x)和h(x)分别输入图像和hdr图像的x通道图像，最后融合h(x)即得到hdr图像h。

本发明的有益成果

本发明提出基于k-svd字典学习和压缩感知的高动态范围图像去伪影融合的方法。k-svd字典学习生成压缩感知中第一步信号稀疏表示所需要的字典，通常情况下采用的固定字典，如dct字典、haar字典、小波字典等等，虽然使用简单计算量低，但是不能保证表达的稀疏程度，只适用于部分类型的图像。而通过k-svd字典学习生成的字典是全局最优的，因此压缩感知中的信号稀疏表示也是最优的，为后续重构提高了精确度。压缩感知突破了奈奎斯特采样率的限制，边采样变压缩，实现的不再是模拟数字转换(adc)，而是模拟信息转换(aic)，对信号进行高度不完备线性测量高精确的重建，极大的减少了采样率，降低了数据存储和传输的代价。基于patchmatch和秩最小算法对压缩感知后的多曝光图像序列进行去伪影融合，能够很好的去除因相机的抖动和图像中物体的运动造成融合后图像伪影和模糊的现象。本发明能在极大降低采样率，同时，能很好去除融合动态场景时伪影和模糊的问题。

附图说明

图1为本发明基于k-svd字典学习和压缩感知的高动态范围图像去伪影融合方法的原理框图；

图2为k-svd字典学习流程图；

图3为压缩感知重构方法流程图；

图4为hdr图像去伪影融合流程图；

图5为本发明待处理的动态场景的多曝光图像序列；

图6为经过压缩感知重构后的多曝光图像序列

图7为融合后去伪影的hdr图像。

具体实施方式

下面结合附图对本发明进行进一步阐述。

如图1所示，本发明基于k-svd字典学习和压缩感知的高动态范围图像去伪影融合方法：1)对输入多曝光图像序列进行分块压缩采样；2)对压缩采样后的图像块使用正交匹配追踪方法(orthogonalmatchingpursuit,omp)实现ldr图像序列的重构；3)对经过压缩感知后的多曝光图像序列使用基于patchmatch和秩最小化算法进行高动态范围图像的去伪影融合。

下面对各步骤进行详细阐述：

1)如图5所示为动态场景的多曝光图像序列，对输入多曝光图像序列进行压缩采样和重构。

11)以当前输入的ldr图像为参考，进行k-svd字典学习，得到用于稀疏表示的超完备字典ψ，如图2所示；

12)对分块后的图像信号矩阵f进行稀疏表示，得到f＝ψθ，设计观测矩阵φ，则压缩采样的过程可重写为：

y＝φf＝φψθ＝θθ\*mergeformat(1)

所述步骤11)还包括：

111)判断当前的ldr图像是否为灰度图像，若不是，则转换为灰度图像，并进一步转换为double类型的图像；

112)选择所需的图像分块大小，图像块大小通常可选为8×8，16×16，32×32等，将111)中获取的灰度图像进行分块，并重排图像块为矩阵列，形成字典学习的数据集矩阵，也即待稀疏表示的信号矩阵y，可进一步将矩阵y表示为其中n为矩阵y的列数；

113)k-svd算法参数初始化，包括：

numatom：需要训练的字典元素个数；

numiteration：迭代次数

errorflag：等于0则表示每个信号的稀疏系数个数固定，需要配置参数l；

否则表示信号的稀疏系数不固定，需要配置参数errorgoal；

errorgoal：最大允许表示误差

preservedcatom：等于1则表示字典的第一个原子为常量；

initializationmethod：若为“dataelements”则用信号本身作为初始字典，若为“givenmatrix”则用给定的矩阵作为初始字典，并对初始字典进行归一化处理，初始字典表示为ψ⁽⁰⁾∈r^n×k；

初始化算法迭代次数j＝1；

114)稀疏编码阶段：使用正交匹配追踪算法omp来计算y中每一个列向量在当前字典上的稀疏表示所对应的系数向量xi，通过求解下式中的目标函数：

其中为由系数向量所构成的矩阵，i表示列向量的索引，||||f为frobenius范数，||||0为l0范数，t0为稀疏表示稀疏中非零分量的数目的上限，即系数向量中的最大差异度；

115)字典更新阶段：初始字典往往不是最优的，满足稀疏性的稀疏矩阵表示的数据和原数据可能会有比较大的误差，因此需在满足稀疏的条件下逐行逐列进行字典的更新优化，减小整体误差，逼近最优字典。对于ψ^(j-1)中的每一列k＝1,2,...k，定义使用到字典中第k个原子ψk的所有信号集合{yi}的索引所构成的集合为：

其中为系数矩阵x中原子ψk所对应的第k行，不同于x的第k列xk的转置，则目标函数式(2)可重写为：

上式中，字典与系数矩阵的乘积ψx被分解为k个秩为1的矩阵的和的形式。当前更新的矩阵列为第k列，则其余k-1个项固定，剥离第k个原子的贡献，得到的矩阵ek为去除原子ψk的成分其余原子所造成的误差；

116)由于中0的存在，用svd得到的更新向量中的非零元素的数量和位置和原中的都不相同，即产生了“发散”，因此需要去除中所有0元素，定义矩阵ωk为n×ωk的矩阵，它在(ωk(i),i)处的值为1，其他位置上的值都为0，则我们定义y、ek去除零项后的收缩结果为的长度为|ωk|，是用到原子ψk的信号的集合，为剥离原子ψk的误差，两者都为n×|ωk|大小的矩阵。则式(4)可进一步转换为：

117)对进行svd分解，即矩阵u的第一列表示为即为列ψk的更新项，更新的为矩阵v的第一列乘上δ(1,1)。j＝j+1进入下一次迭代直至满足迭代条件为止，迭代条件根据设定，可选择为选定的迭代次数或者是误差值，迭代完成后，即得到最优的超完备字典ψ。

所述步骤12)还包括：

121)观测矩阵φ的行数为选择的图像块大小的平方乘上采样率，用m表示，列数定义为图像块的大小的平方，用n表示。此外传感矩阵θ＝φψ应满足rip性质(有限等距性质)：保证观测矩阵不会把两个不同的稀疏信号映射到同一个集合中(保证原空间到稀疏空间的一一映射关系)。定义观测矩阵φ的rip常数δf为满足下式的最小值：

其中f即为分块后的图像信号矩阵，也即稀疏表示后的信号，假设为f稀疏的，即含有f个非零元素，若δf≤1，测成矩阵满则f阶rip条件。

2)对压缩采样后的图像块进行ldr图像序列的重构，重构流程如图3所示，压缩感知重构后的多曝光图像序列如图6所示。

21)初始化：令标签索引集残差向量r0＝y，y即为步骤a2)中所获得的观测向量，迭代次数t＝1；

22)辨识：在传感矩阵θ中求出与残差向量rt-1最相关的所对应的列向量λt，即内积值最大时所对应的索引(列序号)：

其中n为传感矩阵的列数，θj表示传感矩阵的第j列；

23)更新：更新标签索引集λt＝λt-1∪{λt}，将找到的原子所对应的列向量添加到集合其中，t＞0，λt表示t次迭代的索引，θt表示按索引λt选出的矩阵θ的列集合(大小为m×t的矩阵)；

24)估计：求y＝θθ的最小二乘解：

25)更新残差向量

26)令t＝t+1；并不断重复步骤22)至步骤26)。若满足某个迭代停止条件，则停止迭代并进入第27)步；

27)求得输出系数向量

所述步骤26)还包括：

在上述算法的步骤26)中，其迭代停止条件有以下常见的3种情况：

261)当运行到t＞s时，迭代停止，其中，s表示固定的迭代步数，如m/4；

262)残差向量的能量值小于某个预先给定的常数ε，例如取1e-6；

||rt||2≤ε\*mergeformat(11)

263)当传感矩阵θ的任何一列都没有残差向量rt的明显能量时：

||θ^trt||∞≤ε\*mergeformat(12)

3)对经过压缩感知后的多曝光图像序列进行高动态范围图像的去伪影融合，如图4所示。

31)定义输入的经过压缩感知后的图像中任意一点处的像素值为b，对灰度值b归一化处理之后该点的灰度值为i：

i＝b/255\*mergeformat(13)

32)现假设源图像为s，参考图像r，源图像和参考图像合成的图像为l，patchmatch算法就是一个以参考图像为模板，配准源图像生成图像l的过程；由于patchmatch算法是处理一对图像，现假设输入图像为i1...in，以n＝5为例，首先令i3为参考图像r，则i3和i4作为其源图像s，然后令i2和i4作为参考图像r，i1和i5分别作为i2和i4对应的源图像s；

33)现定义patchmatch算法合成图像l的二次函数：

其中τ为灰度映射函数。ω为图像r和图像s的图像域，i为图像域上的任意一个像素点，n(i)为以i中心p×p的邻域，其中p为邻域的大小，取p＝7，故j是邻域n(i)上的像素点，r(i)是图像r上的第i个像素点，s(i+u(j))是图像s上的第(i+u(j))个像素点，其中u(j)表示从图像l上的像素点j映射到图像s的偏移量。α为一个归一化的因子其中wτ和wu为一对权重函数，wτ(i)表示图像r中像素点i映射到图像l的比重，wu(j)表示图像s中像素点j加上偏移量u(j)映射到图像l的比重；

34)灰度映射函数τ定义如下：

其中灰度映射函数的导数τ′≥0，τ()∈[0,1]，i为图像域ω上的像素点，故l(i)表示为图像l上的第i个像素点。使用迭代重加权最小二乘法算法求解灰度映射函数，则可将灰度映射函数τ的目标函数重写为：

其中求解目标函数τ过程中，τ和权重因子ω更新为：

其中n表示迭代次数，δ为一个很小的正常数，令δ＝10^-10；

35)权重函数wτ定义如下：

当图像r中图像域上亮度太暗或太亮(即图像域的像素点的灰度值小于3/255，或者大于252/255)时，像素点将会被clipped，否则，就不clipped；

36)权重函数wμ定义如下：

其中d(·,·)表示输入参数之间的空间距离；υ1,υ2为两个归一化参数，分别取对应空间距离的75百分位数；

37)对于参数x和y，d(x,y)＝||x-y||²，而对于表示为：在图像r和图像s的图像域上的任意一个像素点i，取以i为中心，大小为p×p的邻域，得到图像块和然后图像块经过灰度映射函数得到而图像块相对于i平移了u(i)得到最后求两者的空间距离；同理可得其中τ^-1()灰度映射函数的逆函数；

38)通过上面定义的函数可知，patchmatch算法实际上就是求解二次函数的过程。输入图像r和图像s，并分别对两幅图像进行向下采样，分别得到图像r和图像s的金字塔图像集，从金字塔顶端向下迭代，求得在对应每层金字塔图像下合成的图像l和灰度映射函数τ，将此结果作为下一次的迭代的初始值，当迭代完成后即可得到最终的配准图像l，依照此方法，即可以得到输入图像i1in配准后的图像l1ln；

39)利用伽马曲线对已经配准好的图像集进行辐射校准，消除潜在的相机的移动而产生的噪音.定义伽马(gamma)函数为：

gamma＝cr^γ\*mergeformat(21)

其中c和γ为常数，c取为1，γ取为2.2；这里将经patchmatch配准过的图像集用gamma曲线进行辐射校准；

310)对配准后的图像集使用秩最小化算法得到批量的对齐图像。首先列向量化所有输入的矩阵得到矩阵d，并初始化低秩矩阵a和噪声矩阵e，根据增广拉格朗日乘子法经过内迭代和外迭代得到最佳低秩矩阵a，则噪声矩阵e＝d-a。最后对得到的低秩矩阵和噪声矩阵进行调整m×n大小的图像，即可得到输入图像l对应的低秩图像和噪声图像；

311)输入对齐后的图像集a，将图像a合成目标hdr图像：

其中nimg表示为输入图像的数量，x∈{r,g,b}，r，g，b为彩色图像的三个通道；a(x)和h(x)分别输入图像和hdr图像的x通道图像，最后融合h(x)即可得到hdr图像h，如图7所示。

以上所述仅为本发明的较佳实施方式，本发明并不局限于上述实施方式，在实施过程中可能存在局部微小的结构改动，如果对本发明的各种改动或变型不脱离本发明的精神和范围，且属于本发明的权利要求和等同技术范围之内，则本发明也意图包含这些改动和变型。