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等等；還有針對增強后的光暈現象先用進行照度分割，然后在再計算等等方法，本質上都大同小異。這些改進算法這里不再進行介紹，有興趣的朋友可以去下載些論文看看。有不明白的或者本文有錯誤的地方也希望各位能夠指出，以免誤導后面的讀者。