Retinex图像增强算法

前一段时间研究了一下图像增强算法，发现Retinex理论在彩色图像增强、图像去雾、彩色图像恢复方面拥有很好的效果，下面介绍一下我对该算法的理解。

Retinex理论

Retinex理论始于Land和McCann于20世纪60年代作出的一系列贡献,其基本思想是人感知到某点的颜色和亮度并不仅仅取决于该点进入人眼的绝对光线，还和其周围的颜色和亮度有关。Retinex这个词是由视网膜(Retina)和大脑皮层(Cortex)两个词组合构成的.Land之所以设计这个词，是为了表明他不清楚视觉系统的特性究竟取决于此两个生理结构中的哪一个，抑或是与两者都有关系。

Land的Retinex模型是建立在以下的基础之上的：

一、真实世界是无颜色的，我们所感知的颜色是光与物质的相互作用的结果。我们见到的水是无色的，但是水膜—肥皂膜却是显现五彩缤纷，那是薄膜表面光干涉的结果；

二、每一颜色区域由给定波长的红、绿、蓝三原色构成的；

三、三原色决定了每个单位区域的颜色。

Retinex 理论的基本内容是物体的颜色是由物体对长波（红）、中波（绿）和短波（蓝）光线的反射能力决定的，而不是由反射光强度的绝对值决定的；物体的色彩不受光照非均性的影响，具有一致性，即Retinex理论是以色感一致性（颜色恒常性）为基础的。如下图所示，观察者所看到的物体的图像S是由物体表面对入射光L反射得到的，反射率R由物体本身决定，不受入射光L变化。

图1.Retinex理论中图像的构成

Retinex理论的基本假设是原始图像S是光照图像L和反射率图像R的乘积,即可表示为下式的形式：

基于Retinex的图像增强的目的就是从原始图像S中估计出光照L，从而分解出R，消除光照不均的影响，以改善图像的视觉效果,正如人类视觉系统那样。在处理中，通常将图像转至对数域，即，从而将乘积关系转换为和的关系：

Retinex方法的核心就是估测照度L，从图像S中估测L分量，并去除L分量，得到原始反射分量R，即：

函数实现对照度L的估计（可以去这么理解，实际很多都是直接估计r分量）。

Retinex理论的理解

如果大家看论文，那么在接下去的篇幅当中，肯定会介绍两个经典的Retinex算法：基于路径的Retinex以及基于中心/环绕Retinex。在介绍两个经典的Retinex算法之前，我先来讲一点个人的理解，以便第一次接触该理论的朋友能够更快速地理解。当然，如果我的理解有问题，也请大家帮忙指出。

Retinex理论就我理解，与降噪类似，该理论的关键就是合理地假设了图像的构成。如果将观察者看到的图像看成是一幅带有乘性噪声的图像，那么入射光的分量就是一种乘性的，相对均匀，且变换缓慢的噪声。Retinex算法所做的就是合理地估计图像中各个位置的噪声，并除去它。

在极端情况下，我们大可以认为整幅图像中的分量都是均匀的，那么最简单的估计照度L的方式就是在将图像变换到对数域后对整幅图像求均值。因此，我设计了以下算法来验证自己的猜想，流程如下：

（1） 将图像变换到对数域 ；

（2） 归一化去除加性分量 ；

（3） 对步骤3得到的结果求指数，反变换到实数域。

这里为了简化描述，省略了对图像本身格式的变换，算法用Matlab实现：

% ImOriginal:原始图像
% type:'add'表示分量是加性的，如雾天图像；'mult'表示分量是乘性的，如对照度的估计
[m,n,z] = size(ImOriginal);
ImOut = uint8(zeros(m,n,z));
for i = 1:z
    if strcmp(type,'add')
        ImChannel = double(ImOriginal(:,:,i))+eps;
    elseif strcmp(type,'mult')
        ImChannel = log(double(ImOriginal(:,:,i))+eps);
    else
        error('type must be ''add'' or ''mult''');
    end
    ImOut(:,:,i) = EnhanceOneChannel(ImChannel);
end
ImOut = max(min(ImOut,255), 0);
end

function ImOut = EnhanceOneChannel(ImChannel)
% 计算计算单个通道的反射分量
% 1.对全图进行照射分量估计
% 2.减去照射分量
% 3.灰度拉伸
ImChannel = ImChannel./max(ImChannel(:));
ImRetinex = round(exp(ImChannel.*5.54));
ImOut = uint8(ImRetinex);
end

为了验证算法的有效性，这里使用经典的Retinex算法与我所用的算法进行对比试验，效果如下：

图2.测试原图

图3.经典Retinex算法结果

图4.上述方法结果

从对比中可以看到，对于去除照度，还原图像本身来讲，效果还可以，并且不会在边缘位置产生光晕现象。缺点就是在去除照度分量L过程中，保留的反射分量R我在上述算法中使用归一化后直接进行反变换。这一步的作用可以近似看成去除一个均匀的直流分量，即均匀的照度分量。由于操作都是全局的，这里默认假设了所有位置的照射分量都是相同的，因此在灰度拉伸的时候没有照顾到局部的特性，图像整体亮度偏暗。当然，全局的照度估计对于图像的增强肯定有相当的局限性，其增强效果在色彩的还原和亮度处理等方面还是有一定缺陷的。

个人认为，Retinex算法的关键还是正确的分析了噪声的性质。相信很多人都看到利用基于Retinex的图像水下增强、基于Retinex的图像去雾等等，我也好奇，那就试试吧。大雾图片嘛谁没有，前一阵子大雾天，没事拍了几张照片，终于用上了，请看对比图：

图5.有雾原图

图6.经典Retinex去雾效果

图7.上述方法去雾效果

还是老规矩，这时候对比试验还是最能说明效果的，为此选了一幅干扰很大的图像。基本上各位要是显示器比较差一点，从原图当中是很难看出大雾后面的东西。从去雾效果来看，上述方法的效果并不比经典算法要差，至少在去雾的效果上，本实验结果从主观上给人的感觉还是不错的。

在上述案例中，最重要的就是正确分析有雾图像的结构，与Retinex理论一开始的核心思想有区别的是，在针对这种加性的干扰时，经典的Retinex算法在处理过程中，其实仅仅是利用其估计加性的干扰分量；当然，抛开Retinex理论对照度、反射率对最终图像形成的核心思想（如图1），后续最重要的就是对这个加性的干扰的估计了。

对于有雾的图像，我们大可以看作透过一块磨砂玻璃去看一幅清晰的图像，这样大家就能很好理解为什么认为在这个案例中，将雾的干扰认为是一个加性的了。诸如后面两个经典的算法，所有的这类算法归根结底就是更好地利用原图像中的像素点去估计原始照度。从上面例程上可以看出，使用一个全局估计对局部的增强是比较差的，如果存在照度不均匀（雾的浓度不均匀），或者背景颜色亮度很高等情况时，处理结果会趋向恶劣，效果比较差。

当然，经典也不是完美，从图3中可以看到，经典的算法容易出现光晕效果（蓝色书本文字周围一圈白色），因此后续对于照度估计和去除光晕等问题又有很多基于Retinex理论的变体算法，这里暂不进行介绍。下面开始介绍两个经典的算法，查看Matlab代码下载点击：

http://www.cs.sfu.ca/~colour/publications/IST-2000/

McCann Retinex算法

McCann Retinex 算法是McCann和Frankle一起提出的一种Retinex算法，该算法是一种基于多重迭代策略的Retinex算法，单个点的像素值取决于一条特定路径的环绕的结果，经过多次迭代逼近理想值。通过比较螺旋式路径上的各像素点的灰度值来估计和去除图像的照度分量。

图8.McCann Retinex算法路径选择示意图

这个图是参照人家论文的，不过我准备像论文那样讲，因为太复杂了，不容易懂。从图中我们可以看到，算法沿着这个螺旋式的路径选取用于估计的点，并且越靠近预测的中心点选取的点数越多。这个从实际物理上也解释的通，靠的近的像素点与中心像素点的相关性肯定要要比远处的点要高。

该算法估测图像经过以下几个步骤：

1. 将原图像变换到对数域，彩色图像对各通道都进行对数变换；

2. 初始化常数图像矩阵，该矩阵作为进行迭代运算的初始值；

3. 计算路径，如上图8所示，这里令为路径上的点，从远到近排列；

4. 对路径上的像素点按照如下公式运算：

公式所表示的大致意思为：从远到近，中心点像素值减去路径上的像素值得到的差值的一半与前一时刻的估计值之间的和。最终，中心像素点的像素大致的形式为

其中表示中心位置最终的反射率估计，为常数值为转换后的图像中的最大值，在步骤2中被确定。从这里将按照Retinex理论进行分解，最终公式中可以看到，最终照度分量被去除了，而中心位置的反射率由路径上各点的反射率之间的差值进行估计，并且从轨迹上可以看到，靠的越近的在最终估计的时候所占比重越大。

就我个人理解，这类估计算法，实质是将中心位置和周围亮度之间的差异作为最终中心位置的反射率的估计。如果中心位置亮度本身高，那么最终的结果还是高亮度的；如果中心位置亮度低，那么中心与其它点的差肯定是负值，最终的结果就比较小，亮度就比较低。这就要求路径上的点需要能够很好地代表整幅图像的特性，如果图像本身很规则，那么可能中央与周围计算的结果无法估测该点像素本身的灰度特性，结果就与预期的可能不一样了。如下图9可以看到，算法实质是估计中心和周边的差值，因此图中原本黑色区域的图像由于与周边差别很小，因此呈现高亮度；而在小型的矩形周围随着距离的增大，差别渐渐变小，因此亮度逐渐升高，无法体现原本黑色像素点处原本的亮度特性。

图9.算法的测试图（来自Retinex in Matlab）

从原始算法中可以看到，还有一个重要的步骤，就是迭代。迭代的作用是尽可能保留中央周边差的分量，原本每次只保留中央周边差的一半（见步骤4中最后的除2的处理），迭代次数越多，保留的分量就越多，迭代n次保留的分量就是，这样局部的信息就更多，相当于降低动态范围的压缩。这样的操作可以使图像更加自然，但是会增加运算量，下图可以更加明显地看出迭代次数的影响：

图10.迭代次数的影响

为了方便各位自己研究，下面我给出该算法源码供大家参考：

function Retinex = retinex_frankle_mccann(L, nIterations)

% RETINEX_FRANKLE_McCANN: 
%         Computes the raw Retinex output from an intensity image, based on the
%         original model described in:
%         Frankle, J. and McCann, J., "Method and Apparatus for Lightness Imaging"
%         US Patent #4,384,336, May 17, 1983
%
% INPUT:  L           - logarithmic single-channel intensity image to be processed
%         nIterations - number of Retinex iterations
%
% OUTPUT: Retinex     - raw Retinex output
%
% NOTES:  - The input image is assumed to be logarithmic and in the range [0..1]
%         - To obtain the retinex "sensation" prediction, a look-up-table needs to
%         be applied to the raw retinex output
%         - For colour images, apply the algorithm individually for each channel
%
% AUTHORS: Florian Ciurea, Brian Funt and John McCann. 
%          Code developed at Simon Fraser University.
%
% For information about the code see: Brian Funt, Florian Ciurea, and John McCann
% "Retinex in Matlab," by Proceedings of the IS&T/SID Eighth Color Imaging 
% Conference: Color Science, Systems and Applications, 2000, pp 112-121.
%
% paper available online at http://www.cs.sfu.ca/~colour/publications/IST-2000/
%
% Copyright 2000. Permission granted to use and copy the code for research and 
% educational purposes only.  Sale of the code is not permitted. The code may be 
% redistributed so long as the original source and authors are cited.

global RR IP OP NP Maximum
RR = L;
Maximum = max(L(:));                                 % maximum color value in the image
[nrows, ncols] = size(L);

shift = 2^(fix(log2(min(nrows, ncols)))-1);          % initial shift
OP = Maximum*ones(nrows, ncols);                     % initialize Old Product

while (abs(shift) >= 1)
   for i = 1:nIterations
      CompareWith(0, shift);                         % horizontal step
      CompareWith(shift, 0);                         % vertical step
   end
   shift = -shift/2;                                 % update the shift
end
Retinex = NP;

function CompareWith(s_row, s_col)
global RR IP OP NP Maximum
IP = OP;
if (s_row + s_col > 0)
   IP((s_row+1):end, (s_col+1):end) = OP(1:(end-s_row), 1:(end-s_col)) + ...
   RR((s_row+1):end, (s_col+1):end) - RR(1:(end-s_row), 1:(end-s_col));
else
   IP(1:(end+s_row), 1:(end+s_col)) = OP((1-s_row):end, (1-s_col):end) + ...
   RR(1:(end+s_row),1:(end+s_col)) - RR((1-s_row):end, (1-s_col):end);
end
IP(IP > Maximum) = Maximum;                          % The Reset operation
NP = (IP + OP)/2;                                    % average with the previous Old Product
OP = NP;                                             % get ready for the next comparison

McCann99 Retinex算法

McCann99 Retinex算法本质上与McCann Retinex算法没有区别，两者不同的就是在McCann99算法中，不再是使用螺旋式的路径来选取估计像素点，而是使用图像金字塔的方式逐层选取像素。算法同样经取点、比较、平均这几步进行迭代运算。

图像金字塔的概念很简单，就是对图像进行下采样，以多分辨率的形式表示图像。最顶层的图像分辨率最低，底层的最高（一般为原图）。

图11.图像金字塔示意图（来自网络）

如上图所示，McCann99算法就是从顶层开始让每一个像素点与其8个相邻的像素进行比较，估计反射率分量；前一层计算结束之后对估计的反射率分类进行插值运算，使上一层估计的结果的图像与金字塔下一层图像尺寸相同，并再一次进行相同的比较运算；最终，对原始图像进行8邻域的比较结束后就能得到最终的结果，即增强后的图像了。其中，比较操作与McCann算法相同，都是相减后求平均，见上一节步骤4中描述。

因此，McCann99算法，此处简化描述为以下几步：

1. 将原图像变换到对数域，彩色图像对各通道都进行对数变换；

2. 初始化（计算图像金字塔层数；初始化常数图像矩阵，该矩阵作为进行迭代运算的初始值）；

3. 从顶层开始，到最后一层进行8邻域比较运算，运算规则与MccCann Retinex算法相同，见上一节步骤4；

4. 第n层运算结束后对第n层的运算结果进行插值，变成原来的两倍，与n+1层大小相同（此处默认n越大越靠近底层）；

5. 当最底层计算完毕得到的即最终增强后的图像。

为了方便各位自己研究，下面我给出该算法源码供大家参考：

function Retinex = retinex_mccann99(L, nIterations)

global OPE RRE Maximum
[nrows ncols] = size(L);                             % get size of the input image
nLayers = ComputeLayers(nrows, ncols);               % compute the number of pyramid layers
nrows = nrows/(2^nLayers);                           % size of image to process for layer 0
ncols = ncols/(2^nLayers);
if (nrows*ncols > 25)                                % not processing images of area > 25
  error('invalid image size.')                       % at first layer
end
Maximum = max(L(:));                                 % maximum color value in the image
OP = Maximum*ones([nrows ncols]);                    % initialize Old Product
for layer = 0:nLayers
   RR = ImageDownResolution(L, 2^(nLayers-layer));   % reduce input to required layer size
  
   OPE = [zeros(nrows,1) OP zeros(nrows,1)];         % pad OP with additional columns
   OPE = [zeros(1,ncols+2); OPE; zeros(1,ncols+2)];  % and rows
   RRE = [RR(:,1) RR RR(:,end)];                     % pad RR with additional columns
   RRE = [RRE(1,:); RRE; RRE(end,:)];                % and rows
  
   for iter = 1:nIterations
     CompareWithNeighbor(-1, 0);                     % North
     CompareWithNeighbor(-1, 1);                     % North-East
     CompareWithNeighbor(0, 1);                      % East
     CompareWithNeighbor(1, 1);                      % South-East
     CompareWithNeighbor(1, 0);                      % South
     CompareWithNeighbor(1, -1);                     % South-West
     CompareWithNeighbor(0, -1);                     % West
     CompareWithNeighbor(-1, -1);                    % North-West
   end
   NP = OPE(2:(end-1), 2:(end-1));
   OP = NP(:, [fix(1:0.5:ncols) ncols]);             %%% these two lines are equivalent with
   OP = OP([fix(1:0.5:nrows) nrows], :);             %%% OP = imresize(NP, 2) if using Image
   nrows = 2*nrows; ncols = 2*ncols;                 % Processing Toolbox in MATLAB
end
Retinex = NP;

function CompareWithNeighbor(dif_row, dif_col)
global OPE RRE Maximum

% Ratio-Product operation
IP = OPE(2+dif_row:(end-1+dif_row), 2+dif_col:(end-1+dif_col)) + ...
     RRE(2:(end-1),2:(end-1)) - RRE(2+dif_row:(end-1+dif_row), 2+dif_col:(end-1+dif_col));
    
IP(IP > Maximum) = Maximum;                          % The Reset step

% ignore the results obtained in the rows or columns for which the neighbors are undefined
if (dif_col == -1) IP(:,1) = OPE(2:(end-1),2); end
if (dif_col == +1) IP(:,end) = OPE(2:(end-1),end-1); end
if (dif_row == -1) IP(1,:) = OPE(2, 2:(end-1)); end
if (dif_row == +1) IP(end,:) = OPE(end-1, 2:(end-1)); end
NP = (OPE(2:(end-1),2:(end-1)) + IP)/2;              % The Averaging operation
OPE(2:(end-1), 2:(end-1)) = NP;

function Layers = ComputeLayers(nrows, ncols)
power = 2^fix(log2(gcd(nrows, ncols)));              % start from the Greatest Common Divisor
while(power > 1 & ((rem(nrows, power) ~= 0) | (rem(ncols, power) ~= 0)))
   power = power/2;                                  % and find the greatest common divisor
end                                                  % that is a power of 2
Layers = log2(power);

function Result = ImageDownResolution(A, blocksize)
[rows, cols] = size(A);                              % the input matrix A is viewed as
result_rows = rows/blocksize;                        % a series of square blocks
result_cols = cols/blocksize;                        % of size = blocksize
Result = zeros([result_rows result_cols]);
for crt_row = 1:result_rows                          % then each pixel is computed as
   for crt_col = 1:result_cols                       % the average of each such block
      Result(crt_row, crt_col) = mean2(A(1+(crt_row-1)*blocksize:crt_row*blocksize, ...
                                       1+(crt_col-1)*blocksize:crt_col*blocksize));
   end
en

总结

      在上述两个经典的估计算法之后经过反对数变换就可以恢复成原格式的图像。通过试验可以发现，这些算法都还是有一定缺陷的，对于照度的估计后续还有很多算法如：单尺度Retinex (SSR)、多尺度Retinex (MSR)、可变框架的Retinex等等；还有针对增强后的光晕现象先用进行照度分割，然后在再计算等等方法，本质上都大同小异。这些改进算法这里不再进行介绍，有兴趣的朋友可以去下载些论文看看。有不明白的或者本文有错误的地方也希望各位能够指出，以免误导后面的读者。