CN102110289B - 基于变分框架的彩色图像对比度增强方法 - Google Patents

Abstract

一种基于变分框架的彩色图像对比度增强方法，其主要步骤：第一步，初始化读入图像，将输入的RGB图像转换到YCbCr颜色空间；第二步，对YCbCr空间中彩色图像的Y分量进行降采样；第三步，建立图像的成像模型，并将其变换到对数域；第四步，利用变分框架模型对亮度图像进行计算；第五步，用输入图像减去亮度图像得到图像的反射信息；第六步，将反射图像取指数；第七步，采用插值算法将反射图像还原到原始图像大小；第八步，将得到的反射图像与Cb和Cr分量结合，并反变换到RGB空间，作为输出结果。本发明采用变分框架模型对低能见度天气下的彩色图像进行处理，使得图像质量显著提高，相对于以往的图像增强方法，处理速度快、算法稳定性高、使用范围广、图像增强效果好。

Description

基于变分框架的彩色图像对比度增强方法

技术领域

本发明属于数字图像增强技术领域，具体涉及到一种基于变分框架的彩色图像增强方法。

背景技术

图像增强的主要目的在于使人眼获得更好的视觉效果，或更方便计算机系统识别以满足系统对图像信息的需求。图像特征增强主要包括对比度、亮度、颜色等多方面的处理，目前并没有统一的衡量图像的质量的客观标准，评价算法的有效性主要依靠人眼的客观感受。图像增强的方法基本上可以分为空域处理和频域处理两大类。空域处理是直接对原图像进行数据运算，频域处理是指先将原图像变换到特定的变换域，然后在该变换域上进行相关处理以增强有用的频率分量，最后再对处理后的结果作反变换，便得到增强后的图像。

传统的图像增强算法主要有：线性变换、直方图均衡化等。其中线性变换即比例变换，即因变量和自变量之间存在固定的比例系数，通常情况下线性变换是将图像中较小的灰度范围拉伸到较大的灰度范围，因此常称为灰度拉伸；直方图均衡化算法的基本思想是把原始图像的直方图变换为均匀分布的形式，这样能够增加像素灰度值的动态范围，从而达到增强图像整体对比度的效果。

近年来，基于人眼视觉系统模型的图像增强成为研究热点，在该理论的基础上提出了很多具有现实意义和应用价值的方法，而变分框架模型用于图像增强是一种新的尝试。

发明内容

本发明是一种基于变分框架的彩色图像增强方法，利用变分框架对低能见度环境中的彩色图像进行处理，提高图像质量使之能够用于其它用途，为了实现上述目标，本发明主要采取如下技术方案：

一种基于变分框架的彩色图像对比度增强方法，按照以下步骤进行：

步骤1：初始化，读入一帧大小为M₁×M₂×3的低能见度天气下彩色图像S_in，其中M₁和M₂为正整数，分别表示图像矩阵的行数和列数，然后将输入的彩色图像从RGB彩色空间转换到YCbCr彩色空间，转换后的图像记为S_T，大小为M₁×M₂×3，取S_T中Y分量图像，记为S_Y，大小为M₁×M₂，并将S_Y按照行排成大小为M₁M₂×1的列向量，重排后的图像记为

步骤2：对重排后的图像

先进行降采样得到大小为N₁N₂×1的列向量的降采样图像S_Down，然后再按照行排列成N₁×N₂的图像矩阵S，降采样图像S_Down与

的关系式为

*其中表示矩阵向量乘积，D表示大小为N₁N₂×M₁M₂的降采样矩阵，其具体表达式为

$D=\frac{1}{p^2}{\left(D_{\mathrm{ij}}\right)}_{N_1{N}_2\×M_1{M}_2}=\frac{1}{p^2}{\left[\begin{array}{cccc}{D}_{11}& D_{12}& L& D_{{1N}_2}\\ D_{21}& D_{22}& L& D_{2N_2}\\ M& M& M& M\\ D_{N_{1}1}& D_{N_{1}2}& L& D_{N_1{N}_2}\end{array}\right]}_{N_1{N}_2\×M_1{M}_2}---\left(1\right)$

对于i＝1，2，...，N₁N₂，

其中，p为降采样因子且为正整数，p＝1～4，i、j为降采样矩阵的坐标，D_i，j为降采样矩阵中坐标为i、j的点，N₁和N₂为正整数且分别为重采样的每帧低分辨率图像矩阵的行数和列数，

步骤3：建立图像成像模型，图像成像由亮度信息和反射信息两部分来表示，其成像过程表示为

S＝L·R            (2)

其中S是上步中降采样后得到的图像矩阵S，L表示图像中的亮度信息，R表示图像中的反射信息，将(2)式等号两边同时取自然对数，得到对数域中的图像成像模型：

s＝l+r                        (3)

其中，s＝logS，s是对数域中的图像矩阵，l＝logL，l是对数域中的亮度信息，r＝log R，r是对数域中的反射信息，

步骤4：利用变分方法对亮度信息l进行估计得到相应估计值

具体方法为：

根据对数域中的图像成像模型，将对数域成像模型中的亮度信息l进行变分优化：

$F\left(l\right)=\underset{\mathrm{\Ω}}{\∫}|\▿l|+\mathrm{\α}{(l-s)}^2+\mathrm{\β}{|\▿(l-s)|}^2\mathrm{dxdy}---\left(4\right)$

其中，α和β是任意非负的实数，Ω为图像区域，

为梯度算子，(4)式对应的欧拉-拉格朗日方程为

-Δl+α(l-s)-βΔ(l-s)＝0                        (5)

其中，Δ表示拉普拉斯算子，利用梯度下降流求解该欧拉-拉格朗日方程(5)，其对应的梯度下降流为

$\frac{\∂l(x,y,t)}{\∂t}=-(-\mathrm{\Δl}(x,y,t)+\mathrm{\α}(l(x,y,t)-s(x,y,t))-\mathrm{\β\Δ}(l(x,y,t)-s(x,y,t)))---\left(6\right)$

其中x、y表示图像的坐标，t表示时间，采用半点格式的中心差分来离散化上式，选择中心像素点(i，j)，为了利用像素周围更多的信息，选择中心像素点的八个相邻点(i-1，j-1)、(i-1，j)、(i-1，j+1)、(i，j-1)、(i，j+1)、(i+1，j-1)、(i+1，j)和(i+1，j+1)，设h为x，y方向上两像素间的网格步长，δt为时间步长，记亮度图像采样值

原始图像采样值

中心像素点的一阶导数对应的有限差分格式为：

$\frac{\∂l(x,y,t)}{\∂t}=\frac{l_{i,j}^{n+1}-l_{i,j}^n}{\mathrm{\δt}},$

${\left(l_x\right)}_{i,j}^n=\frac{l_{i+1,j}^n-l_{i-1,j}^n}{2h},$ ${\left(l_y\right)}_{i,j}^n=\frac{l_{i,j+1}^n-l_{i,j-1}^n}{2h},$

${\left(s_x\right)}_{i,j}^n=\frac{s_{i+1,j}^n-s_{i-1,j}^n}{2h},$ ${\left(s_y\right)}_{i,j}^n=\frac{s_{i,j+1}^n-s_{i,j-1}^n}{2h},$

二阶导数对应的有限差分格式：

${\left(l_{\mathrm{xx}}\right)}_{i,j}^n=\frac{l_{i+1,j}^n-2l_{i,j}^n+l_{i-1,j}^n}{h^2},$ ${\left(l_{\mathrm{yy}}\right)}_{i,j}^n=\frac{l_{i,j+1}^n-2l_{i,j}^n+l_{i,j-1}^n}{h^2},$

${\left(s_{\mathrm{xx}}\right)}_{i,j}^n=\frac{s_{i+1,j}^n-2s_{i,j}^n+s_{i-1,j}^n}{h^2}$ ${\left(s_y\right)}_{i,j}^n=\frac{s_{i,j+1}^n-2s_{i,j}^n+s_{i,j-1}^n}{h^2},$

${\left(\mathrm{\Δl}\right)}_{i,j}^n={\left(l_{\mathrm{xx}}\right)}_{i,j}^n+{\left(l_{\mathrm{yy}}\right)}_{i,j}^n,$ ${\left(\mathrm{\Δs}\right)}_{i,j}^n={\left(s_{\mathrm{xx}}\right)}_{i,j}^n+{\left(s_{\mathrm{yy}}\right)}_{i,j}^n;$

故上述梯度下降流(6)的离散迭代格式为

$\frac{l_{i,j}^{n+1}-l_{i,j}^n}{\mathrm{\δt}}=(-{\left(\mathrm{\Δl}\right)}_{i,j}^n+\mathrm{\α}(l_{i,j}^n-s_{i,j}^n)-\mathrm{\β}({\left(\mathrm{\Δl}\right)}_{i,j}^n-{\left(\mathrm{\Δs}\right)}_{i,j}^n))---\left(7\right)$

即

$l_{i,j}^{n+1}=l_{i,j}^n-\mathrm{\δt}(-{\left(\mathrm{\Δl}\right)}_{i,j}^n+\mathrm{\α}(l_{i,j}^n-s_{i,j}^n)-\mathrm{\β}({\left(\mathrm{\Δl}\right)}_{i,j}^n-{\left(\mathrm{\Δs}\right)}_{i,j}))---\left(8\right)$

迭代中止条件为：

$\frac{\left|\right|l^{n+1}-l^n\left|\right|}{\left|\right|l^n\left|\right|}\≤\mathrm{\ϵ}---\left(9\right)$

其中ε是预置的迭代终止系数，取0.0001，

步骤5：用对数域中的降采样图像s减去步骤4中得到的对数域亮度图像估计值

得到对数域反射图像r的估计值

$\hat{r}=s-\hat{l}---\left(10\right)$

步骤6：将上一步中得到的对数域反射图像估计值

逐点取指数exp，得到反射图像

$\hat{R}=\exp \left(\hat{r}\right)---\left(11\right)$

步骤7：对步骤6中的反射图像

采用公知的拉格朗日插值算法进行插值，插值后得到的图像为原始大小M₁×M₂，得到的结果记做

步骤8：将作为Y分量，并与S_T中Cb分量和Cr分量重新组合得到YCbCr空间的彩色图像，再将其转换到RGB空间。

与现有的技术相比，本发明的特点在于：

1、本发明对彩色图像进行变换，选取其中的一个分量进行处理，并且在处理过程中进行降采样和对数变换，使得算法的复杂度降低，提高了算法的运行速度，使之能够用于实时系统，

2、本发明对图像的对比度要求不高，能够有效地增强图像中暗区的信息并且抑制噪声，使得处理后的输出图像的动态范围大，图像细节部分对比度提高，图像质量得到明显改善，

3、本发明中提出的模型对于其各项参数具有良好的鲁棒性，当参数在较大范围内变化时，输出图像基本保持一致，并且对给定的参数能够快速收敛达到稳态，算法稳定性很好。

附图说明

图1是整个程序的流程图；

图2是处理前的沙尘天气下的图像；

图3是处理后的沙尘天气下的图像；

图4是处理前的夜晚环境下的图像；

图5是处理后的夜晚环境下的图像；

图6是处理前的雾天天气下的图像；

图7是处理后的雾天天气下的图像；

图8是处理前的逆光拍摄图像；

图9是处理后的逆光拍摄图像。

具体实施方式

在具体的实施方式中，将结合附图，清楚完整地描述基于变分框架模型的图像增强算法的详细过程，

一种基于变分框架的彩色图像对比度增强方法，按照以下步骤进行：

步骤1：初始化，读入一帧大小为M₁×M₂×3的低能见度天气下彩色图像S_in，其中M₁和M₂为正整数，分别表示图像矩阵的行数和列数，然后将输入的彩色图像从RGB彩色空间转换到YCbCr彩色空间，转换后的图像记为S_T，大小为M₁×M₂×3，取S_T中Y分量图像，记为S_Y，大小为M₁×M₂，并将S_Y按照行排成大小为M₁M₂×1的列向量，重排后的图像记为从RGB彩色空间转换到YCbCr彩色空间的具体过程为：

$\left[\begin{array}{c}Y\\ \mathrm{Cb}\\ \mathrm{Cr}\end{array}\right]=\left[\begin{array}{c}16\\ 128\\ 128\end{array}\right]+(1/256)\left[\begin{array}{ccc}65.738& 129.057& 25.06\\ -37.945& -74.494& 112.43\\ 112.439& -94.154& -18.28\end{array}\right]\left[\begin{array}{c}R\\ G\\ B\end{array}\right]---\left(1\right)$

步骤2：对重排后的图像

先进行降采样得到大小为N₁N₂×1的列向量的降采样图像S_Down，然后再按照行排列成N₁×N₂的图像矩阵S，降采样图像S_Down与

的关系式为

*其中表示矩阵向量乘积，D表示大小为N₁N₂×M₁M₂的降采样矩阵，其具体表达式为

$D=\frac{1}{p^2}{\left(D_{\mathrm{ij}}\right)}_{N_1{N}_2\×M_1{M}_2}=\frac{1}{p^2}{\left[\begin{array}{cccc}{D}_{11}& D_{12}& L& D_{{1N}_2}\\ D_{21}& D_{22}& L& D_{2N_2}\\ M& M& M& M\\ D_{N_{1}1}& D_{N_{1}2}& L& D_{N_1{N}_2}\end{array}\right]}_{N_1{N}_2\×M_1{M}_2}---\left(2\right)$

对于i＝1，2，...，N₁N₂，

其中，p为降采样因子且为正整数，在应用实例中选择p＝3，i、j为降采样矩阵的坐标，D_i，j为降采样矩阵中坐标为i、j的点，N₁和N₂为正整数且分别为重采样的每帧低分辨率图像矩阵的行数和列数，

步骤3：建立图像成像模型，图像成像由亮度信息和反射信息两部分来表示，其成像过程表示为

S＝L·R            (3)

其中S是上步中降采样后得到的图像矩阵S，L表示图像中的亮度信息，R表示图像中的反射信息，将(3)式等号两边同时取自然对数，得到对数域中的图像成像模型：

s＝l+r                (4)

其中，s＝log S，s是对数域中的图像矩阵，l＝log L，l是对数域中的亮度信息，r＝log R，r是对数域中的反射信息，

步骤4：利用变分方法对亮度信息l进行估计得到相应估计值

具体方法为：

根据对数域中的图像成像模型，将对数域成像模型中的亮度信息l进行变分优化：

$F\left(l\right)=\underset{\mathrm{\Ω}}{\∫}|\▿l|+\mathrm{\α}{(l-s)}^2+\mathrm{\β}{|\▿(l-s)|}^2\mathrm{dxdy}---\left(5\right)$

其中，α和β是任意非负的实数，在应用实例中选取α＝0.0001，β＝0.1，Ω为图像区域，

为梯度算子，(5)式对应的欧拉-拉格朗日方程为

-Δl+α(l-s)-βΔ(l-s)＝0                        (6)

其中，Δ表示拉普拉斯算子，利用梯度下降流求解该欧拉-拉格朗日方程(6)，其对应的梯度下降流为

$\frac{\∂l(x,y,t)}{\∂t}=-(-\mathrm{\Δl}(x,y,t)+\mathrm{\α}(l(x,y,t)-s(x,y,t))-\mathrm{\β\Δ}(l(x,y,t)-s(x,y,t)))---\left(7\right)$

其中x、y表示图像的坐标，t表示时间，采用半点格式的中心差分来离散化上式，选择中心像素点(i，j)，为了利用像素周围更多的信息，选择中心像素点的八个相邻点(i-1，j-1)、(i-1，j)、(i-1，j+1)、(i，j-1)、(i，j+1)、(i+1，j-1)、(i+1，j)和(i+1，j+1)，设h为x，y方向上两像素间的网格步长，δt为时间步长，记亮度图像采样值

原始图像采样值中心像素点的一阶导数对应的有限差分格式为：

$\frac{\∂l(x,y,t)}{\∂t}=\frac{l_{i,j}^{n+1}-l_{i,j}^n}{\mathrm{\δt}},$

${\left(l_x\right)}_{i,j}^n=\frac{l_{i+1,j}^n-l_{i-1,j}^n}{2h},$ ${\left(l_x\right)}_{i,j}^n=\frac{l_{i+1,j}^n-l_{i-1,j}^n}{2h},$

${\left(s_x\right)}_{i,j}^n=\frac{s_{i+1,j}^n-s_{i-1,j}^n}{2h},$ ${\left(s_y\right)}_{i,j}^n=\frac{s_{i,j+1}^n-s_{i,j-1}^n}{2h},$

中心像素点的二阶导数对应的有限差分格式：

${\left(l_{\mathrm{xx}}\right)}_{i,j}^n=\frac{l_{i+1,j}^n-2l_{i,j}^n+l_{i-1,j}^n}{h^2},$ ${\left(l_{\mathrm{yy}}\right)}_{i,j}^n=\frac{l_{i,j+1}^n-2l_{i,j}^n+l_{i,j-1}^n}{h^2},$

${\left(s_{\mathrm{xx}}\right)}_{i,j}^n=\frac{s_{i+1,j}^n-2s_{i,j}^n+s_{i-1,j}^n}{h^2},$ ${\left(s_y\right)}_{i,j}^n=\frac{s_{i,j+1}^n-2s_{i,j}^n+s_{i,j-1}^n}{h^2},$

${\left(\mathrm{\Δl}\right)}_{i,j}^n={\left(l_{\mathrm{xx}}\right)}_{i,j}^n+{\left(l_{\mathrm{yy}}\right)}_{i,j}^n,$ ${\left(\mathrm{\Δs}\right)}_{i,j}^n={\left(s_{\mathrm{xx}}\right)}_{i,j}^n+{\left(s_{\mathrm{yy}}\right)}_{i,j}^n;$

故上述梯度下降流(7)的离散迭代格式为

$\frac{l_{i,j}^{n+1}-l_{i,j}^n}{\mathrm{\δt}}=(-{\left(\mathrm{\Δl}\right)}_{i,j}^n+\mathrm{\α}(l_{i,j}^n-s_{i,j}^n)-\mathrm{\β}({\left(\mathrm{\Δl}\right)}_{i,j}^n-{\left(\mathrm{\Δs}\right)}_{i,j}^n))---\left(8\right)$

即

$l_{i,j}^{n+1}=l_{i,j}^n-\mathrm{\δt}(-{\left(\mathrm{\Δl}\right)}_{i,j}^n+\mathrm{\α}(l_{i,j}^n-s_{i,j}^n)-\mathrm{\β}({\left(\mathrm{\Δl}\right)}_{i,j}^n-{\left(\mathrm{\Δs}\right)}_{i,j}))---\left(9\right)$

迭代中止条件为：

$\frac{\left|\right|l^{n+1}-l^n\left|\right|}{\left|\right|l^n\left|\right|}\≤\mathrm{\ϵ}---\left(10\right)$

其中ε是预置的迭代终止系数，在应用实例中选取ε＝0.0001，

步骤5：用对数域中的降采样图像s减去步骤4中得到的对数域亮度图像估计值

得到对数域反射图像r的估计值

$\hat{r}=s-\hat{l}---\left(11\right)$

步骤6：将上一步中得到的对数域反射图像估计值

逐点取指数exp，得到反射图像

$\hat{R}=\exp \left(\hat{r}\right)---\left(12\right)$

步骤7：对步骤6中的反射图像

采用公知的拉格朗日插值算法进行插值，插值后得到的图像为原始大小M₁×M₂，得到的结果记做

步骤8：将

作为Y分量，并与ST中Cb分量和Cr分量重新组合得到YCbCr空间的彩色图像，再将其转换到RGB空间，所得结果作为输出图像R_out，从YCbCr彩色空间到RGB彩色空间的具体转换过程为：

$\left[\begin{array}{c}R\\ G\\ B\end{array}\right]=(1/125)*\left[\begin{array}{ccc}298.082& 0& 408.58\\ 298.082& -100.291& -208.12\\ 298.082& 516.411& 0\end{array}\right]\left[\begin{array}{c}Y\\ \mathrm{Cb}\\ \mathrm{Cr}\end{array}\right]-\left[\begin{array}{c}16\\ 128\\ 128\end{array}\right]---\left(13\right)$

变分框架模型的低能见度图像增强算法的应用试验

在应用试验中，采用了沙尘、夜晚、雾天和逆光拍摄等四种情况下的四幅图像进行处理，选取α＝0.0001，β＝0.1，ε＝0.0001，降采样因子p＝3，图2是处理前的沙尘天气下的图像，图3是处理后的沙尘天气下的图像，图4是处理前的夜晚环境下的图像，图5是处理后的夜晚环境下的图像，图6是处理前的雾天天气下的图像，图7是处理后的雾天天气下的图像，图8是处理前的逆光拍摄图像，图9是处理后的逆光拍摄图像，从前后对比中可以看出处理后的图像变得更加清晰、细节部分更加明显，经本发明处理后的图像取得了较好的直观效果。

Claims (1)

1.一种基于变分框架模型的低能见度天气下彩色图像对比度增强方法，其特征在于按照以下步骤进行：

步骤1：初始化，读入一帧大小为M₁×M₂×3的低能见度天气下彩色图像S_in，其中M₁和M₂为正整数，分别表示图像矩阵的行数和列数，然后将输入的彩色图像从RGB彩色空间转换到YCbCr彩色空间，转换后的图像记为S_T，大小为M₁×M₂×3，取S_T中Y分量图像，记为S_Y，大小为M₁×M₂，并将S_Y按照行排成大小为M₁M₂×1的列向量，重排后的图像记为

步骤2：对重排后的图像

先进行降采样得到大小为N₁N₂×1的列向量的降采样图像S_Down，然后再按照行排列成N₁×N₂的图像矩阵S，降采样图像S_Down与

的关系式为

*其中表示矩阵向量乘积，D表示大小为N₁N₂×M₁M₂的降采样矩阵，其具体表达式为

$D=\frac{1}{p^1}{\left(D_{\mathrm{ij}}\right)}_{N_1{N}_2\×M_1{M}_2}=\frac{1}{p^2}{\left[\begin{array}{cccc}{D}_{11}& D_{12}& L& D_{{1N}_2}\\ D_{21}& D_{22}& L& D_{2N_2}\\ M& M& M& M\\ D_{N_{1}1}& D_{N_{1}2}& L& D_{N_1{N}_2}\end{array}\right]}_{N_1{N}_2\×M_1{M}_2}---\left(1\right)$

对于i＝1，2，...，N₁N₂，

其中，p为降采样因子且为正整数，并选择p＝3，i、j为降采样矩阵的坐标，D_i，j为降采样矩阵中坐标为i、j的点，N₁和N₂为正整数且分别为重采样的每帧低分辨率图像矩阵的行数和列数，

步骤3：建立图像成像模型，图像成像由亮度信息和反射信息两部分来表示，其成像过程表示为

S＝L·R                       (2)

其中S是上步中降采样后得到的图像矩阵S，L表示图像中的亮度信息，R表示图像中的反射信息，将(2)式等号两边同时取自然对数，得到对数域中的图像成像模型：

s＝l+r                        (3)

其中，s＝logS，s是对数域中的图像矩阵，l＝logL，l是对数域中的亮度信息，r＝logR，r是对数域中的反射信息，

步骤4：利用变分方法对亮度信息l进行估计得到相应估计值

具体方法为：根据对数域中的图像成像模型，将对数域成像模型中的亮度信息l进行变分优化：

$F\left(l\right)=\underset{\mathrm{\Ω}}{\∫}|\▿l|\mathrm{\α}{(l-s)}^2+\mathrm{\β}{|\▿(l-s)|}^2\mathrm{dxdy}---\left(4\right)$

其中，α＝0.0001，β＝0.1，Ω为图像区域，为梯度算子，(4)式对应的欧拉拉格朗日方程为

-Δl+α(l-s)-βΔ(l-s)＝0            (5)

其中，Δ表示拉普拉斯算子，利用梯度下降流求解该欧拉拉格朗日方程(5)，其对应的梯度下降流为

$\frac{\∂l(x,y,t)}{\∂t}=-(-\mathrm{\Δl}(x,y,t)+\mathrm{\α}(l(x,y,t)-s(x,y,t))-\mathrm{\β\Δ}(l(x,y,t)-s(x,y,t)))---\left(6\right)$

其中x、y表示图像像素点的坐标，t表示时间，采用半点格式的中心差分来离散化上式，选择中心像素点(i，j)，为了利用像素周围更多的信息，选择中心像素点的八个相邻点(i-1，j-1)、(i-1，j)、(i-1，j+1)、(i，j-1)、(i，j+1)、(i+1，j-1)、(i+1，j)和(i+1，j+1)，设h为x，y方向上两像素间的网格步长，δt为时间步长，记亮度图像采样值

原始图像采样值

中心像素点的一阶导数对应的有限差分格式为：

$\frac{\∂l(x,y,t)}{\∂t}=\frac{l_{i,j}^{n+1}-l_{i,j}^n}{\mathrm{\δt}},$

${\left(l_x\right)}_{i,j}^n=\frac{l_{i+1,j}^n-l_{i-1,j}^n}{2h},$ ${\left(l_y\right)}_{i,j}^n=\frac{l_{i,j+1}^n-l_{i,j-1}^n}{2h},$

${\left(s_x\right)}_{i,j}^n=\frac{s_{i+1,j}^n-s_{i-1,j}^n}{2h},$ ${\left(s_y\right)}_{i,j}^n=\frac{s_{i,j+1}^n-s_{i,j-1}^n}{2h},$

二阶导数对应的有限差分格式：

${\left(l_{\mathrm{xx}}\right)}_{i,j}^n=\frac{l_{i+1,j}^n-2l_{i,j}^n+l_{i-1,j}^n}{h^2},$ ${\left(l_{\mathrm{yy}}\right)}_{i,j}^n=\frac{l_{i,j+1}^n-2l_{i,j}^n+l_{i,j-1}^n}{h^2},$

${\left(s_{\mathrm{xx}}\right)}_{i,j}^n=\frac{s_{i+1,j}^n-2s_{i,j}^n+s_{i-1,j}^n}{h^2},$ ${\left(s_{\mathrm{yy}}\right)}_{i,j}^n=\frac{s_{i,j+1}^n-2s_{i,j}^n+s_{i,j-1}^n}{h^2},$

${\left(\mathrm{\Δl}\right)}_{i,j}^n={\left(l_{\mathrm{xx}}\right)}_{i,j}^n+{\left(l_{\mathrm{yy}}\right)}_{i,j}^n,$ ${\left(\mathrm{\Δs}\right)}_{i,j}^n={\left(s_{\mathrm{xx}}\right)}_{i,j}^n+{\left(s_{\mathrm{yy}}\right)}_{i,j}^n;$

故上述梯度下降流(6)的离散迭代格式为

$\frac{l_{i,j}^{n+1}-l_{i,j}^n}{\mathrm{\δt}}=-(-{\left(\mathrm{\Δl}\right)}_{i,j}^n+\mathrm{\α}(l_{i,j}^n-s_{i,j}^n)-\mathrm{\β}({\left(\mathrm{\Δl}\right)}_{i,j}^n-{\left(\mathrm{\Δs}\right)}_{i,j}^n))---\left(7\right)$

即

$l_{i,j}^{n+1}=l_{i,j}^n-\mathrm{\δt}(-{\left(\mathrm{\Δl}\right)}_{i,j}^n+\mathrm{\α}(l_{i,j}^n-s_{i,j}^n)-\mathrm{\β}({\left(\mathrm{\Δl}\right)}_{i,j}^n-{\left(\mathrm{\Δs}\right)}_{i,j}))---\left(8\right)$

迭代中止条件为：

$\frac{\left|\right|l^{n+1}-l^n\left|\right|}{\left|\right|l^n\left|\right|}\≤\mathrm{\ϵ}---\left(9\right)$

其中ε是预置的迭代终止系数，ε取值越小效果越好，一般取0.0001，

步骤5：用对数域中的降采样图像s减去步骤4中得到的对数域亮度图像估计值得到对数域反射图像r的估计值

$\hat{r}=s-\hat{l}---\left(10\right)$

步骤6：将上一步中得到的对数域反射图像估计值

逐点取指数exp，得到反射图像

$\hat{R}=\exp \left(\hat{r}\right)---\left(11\right)$

步骤7：对步骤6中的反射图像

采用公知的拉格朗日插值算法进行插值，插值后得到的图像为原始大小M₁×M₂，得到的结果记做

步骤8：将

作为Y分量，并与S_T中Cb分量和Cr分量重新组合得到YCbCr空间的彩色图像，再将其转换到RGB空间。

Images