Retinex系列之McCann99 Retinex 分类：图像处理 ...-白红宇

Retinex系列之McCann99 Retinex 分类：图像处理 ...

阅读量：5945 次

发布时间：2019-06-19

本文共 5665 字，大约阅读时间需要 18 分钟。

一、McCann99 Retinex

McCann99利用金字塔模型建立对图像的多分辨率描述,自顶向下逐层迭代，提高增强效率。对输入图像的长宽有

严格的限制，要求可表示成 $nrows\cdot 2^{n}\times ncols\cdot2^{n}$ ，且 $cols\geq rows$ ， $1\leq cols,rows\leq5$ 。

上述限制来源于金字塔模型的结构要求，由于要对输入图像进行下采样，金字塔中上层低分辨率图像的宽分别为下

层高分辨率图像的1/2，顶层(第n层)大小为 $rows\times cols$ ，底层(第0层)为原图像。金字塔结构如下图所示。

McCann99算法对输入图像的尺寸要求过于严格，以至于大部分图像不能直接用此算法进行增强，后续有很多改进

措施，此处暂不考虑。

算法从顶层开始，将每个像素与其8领域像素比较，估计每个像素点的亮度值(lightness)，比较的迭代次数nIterator

由用户决定，每次迭代有四步操作：比例(ratio)、乘积(product)、重置(reset)、平均(average)。随着迭代次数的增

大，像素受到邻域的影响范围就会扩大。每层计算结束后对该层所得图像进行插值运算，使其尺寸与下一层相同，并

将插值结果作为下一层处理的初始值。对下一层进行相同操作，这样自顶向下直至金字塔底层得到最终增强结果。

对每金字塔每一层：设OP为上一步迭代的乘积；NP为当前迭代的乘积；IP为中间乘积结果；R为该层输入图像；符号

*表示重置操作。其中对于每个像素执行的四步操作使用和Frankle-McCann相同的公式：

$NP(x,y)=\sqrt{OP(x,y)IP^{*}(x,y)}=\sqrt{OP(x,y)[OP(xs,ys)\frac{R(x,y)}{R(xs,ys)}]^{*}}$

变量初始化：

金字塔层数由原始图像尺寸决定，大小为

$nrows\cdot 2^{n}\times ncols\cdot2^{n}$ ，层数为n+1；

OP初始尺寸为 $rows\times cols$ ，各像素值一般设为原始图像中的最大值；

R：第k层输入图为像原始图像的下采样，大小为 $nrows\cdot 2^{n-k}\times ncols\cdot2^{n-k}$ ；

迭代次数nIterator一般设为4；

图像金字塔模型

二、Matlab实现

function Test()ImOriginal=imread('fig5.tif');[m,n,z] = size(ImOriginal);ImOut = zeros(m,n,z);for i = 1:z    ImChannel = log(double(ImOriginal(:,:,i))+eps);    ImOut(:,:,i)=retinex_mccann99(ImChannel,4);     ImOut(:,:,i)=exp(ImOut(:,:,i));    a=min(min(ImOut(:,:,i)));    b=max(max(ImOut(:,:,i)));    ImOut(:,:,i)=((ImOut(:,:,i)-a)/(b-a))*255;     endImOut=uint8(ImOut);figure(1);imshow(ImOriginal);figure(2);imshow(ImOut);function Retinex = retinex_mccann99(L, nIterations)% INPUT:  L           - logarithmic single-channel intensity image to be processed%         nIterations - number of Retinex iterations%% OUTPUT: Retinex     - raw Retinex outputglobal OPE RRE Maximum[nrows ncols] = size(L);                             % get size of the input imagenLayers = ComputeLayers(nrows, ncols);               % compute the number of pyramid layersnrows = nrows/(2^nLayers);                           % size of image to process for layer 0ncols = ncols/(2^nLayers);if (nrows*ncols > 25)                                % not processing images of area > 25  error('invalid image size.')                       % at first layerendMaximum = max(L(:));                                 % maximum color value in the imageOP = Maximum*ones([nrows ncols]);                    % initialize Old Productfor 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 MATLABendRetinex = NP;%将当前像素与八邻域比较function CompareWithNeighbor(dif_row, dif_col)global OPE RRE Maximum% Ratio-Product operationIP = 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%因OPE边界处填充了0，故IP对应的边界处结果无意义，直接置成原值if (dif_col == -1) IP(:,1) = OPE(2:(end-1),2); endif (dif_col == +1) IP(:,end) = OPE(2:(end-1),end-1); endif (dif_row == -1) IP(1,:) = OPE(2, 2:(end-1)); endif (dif_row == +1) IP(end,:) = OPE(end-1, 2:(end-1)); endNP = (OPE(2:(end-1),2:(end-1)) + IP)/2;              % The Averaging operationOPE(2:(end-1), 2:(end-1)) = NP;%power:nrows,ncols的最大公约数且是2的整数次方function Layers = ComputeLayers(nrows, ncols)power = 2^fix(log2(gcd(nrows, ncols)));              % start from the Greatest Common Divisorwhile(power > 1 && ((rem(nrows, power) ~= 0) || (rem(ncols, power) ~= 0)))   power = power/2;                                  % and find the greatest common divisorend                                                  % that is a power of 2Layers = log2(power);%下采样，将blocksize*blocksize区域映射为一个像素点function Result = ImageDownResolution(A, blocksize)[rows, cols] = size(A);                              % the input matrix A is viewed asresult_rows = rows/blocksize;                        % a series of square blocksresult_cols = cols/blocksize;                        % of size = blocksizeResult = 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));   endend

测试结果：

输入

输出

注：输出时只是简单的进行线性拉伸，使得灰度值落在[0-255],没有使用更好调整方法。

参考：

[1] J.J. McCann, “LessonsLearned from Mondrians Applied to Real Images and Color Gamuts”, Proc.IS&T/SID Seventh Color Imaging Conference, pp. 1-8, 1999.

[2] Brian Funt, FlorianCiurea, and John McCann "Retinex in Matlab," Proceedings of the IS&T/SIDEighth Color Imaging Conference: Color Science, Systems and Applications, 2000.

转载于:https://www.cnblogs.com/luo-peng/p/4646232.html

你可能感兴趣的文章