retinex 的水下圖像增強算法_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代碼下載下傳點選：

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 describedin:

% 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 andJohn McCann.

% Code developed at Simon Fraser University.

%

% For information about the codesee: Brian Funt, Florian Ciurea, andJohn 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 athttp://www.cs.sfu.ca/~colour/publications/IST-2000/

%

% Copyright2000. Permission granted to use and copy the code for research and% educational purposes only. Sale of the code isnotpermitted. The code may be

% redistributed so long as the original sourceandauthors 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 > 25error('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 obtainedin the rows orcolumns 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 2Layers = 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等等；還有針對增強後的光暈現象先用進行照度分割，然後在再計算等等方法，本質上都大同小異。這些改進算法這裡不再進行介紹，有興趣的朋友可以去下載下傳些論文看看。有不明白的或者本文有錯誤的地方也希望各位能夠指出，以免誤導後面的讀者。