本文主要是介绍Retinex系列之McCann99 Retinex,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!
一、McCann99 Retinex
McCann99利用金字塔模型建立对图像的多分辨率描述,自顶向下逐层迭代,提高增强效率。对输入图像的长宽有
严格的限制,要求可表示成 ,且 ,。
上述限制来源于金字塔模型的结构要求,由于要对输入图像进行下采样,金字塔中上层低分辨率图像的宽分别为下
层高分辨率图像的1/2,顶层(第n层)大小为,底层(第0层)为原图像。金字塔结构如下图所示。
McCann99算法对输入图像的尺寸要求过于严格,以至于大部分图像不能直接用此算法进行增强,后续有很多改进
措施,此处暂不考虑。
算法从顶层开始,将每个像素与其8领域像素比较,估计每个像素点的亮度值(lightness),比较的迭代次数nIterator
由用户决定,每次迭代有四步操作:比例(ratio)、乘积(product)、重置(reset)、平均(average)。随着迭代次数的增
大,像素受到邻域的影响范围就会扩大。每层计算结束后对该层所得图像进行插值运算,使其尺寸与下一层相同,并
将插值结果作为下一层处理的初始值。对下一层进行相同操作,这样自顶向下直至金字塔底层得到最终增强结果。
对每金字塔每一层:设OP为上一步迭代的乘积;NP为当前迭代的乘积;IP为中间乘积结果;R为该层输入图像;符号
*表示重置操作。其中对于每个像素执行的四步操作使用和Frankle-McCann相同的公式:
变量初始化:
金字塔层数由原始图像尺寸决定,大小为 ,层数为n+1;OP初始尺寸为,各像素值一般设为原始图像中的最大值;
R:第k层输入图为像原始图像的下采样,大小为;
迭代次数nIterator一般设为4;
图像金字塔模型
二、Matlab实现
function Test()
ImOriginal=imread('fig5.tif');
[m,n,z] = size(ImOriginal);
ImOut = zeros(m,n,z);
for i = 1:zImChannel = 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;
end
ImOut=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 output
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:nLayersRR = ImageDownResolution(L, 2^(nLayers-layer)); % reduce input to required layer sizeOPE = [zeros(nrows,1) OP zeros(nrows,1)]; % pad OP with additional columnsOPE = [zeros(1,ncols+2); OPE; zeros(1,ncols+2)]; % and rowsRRE = [RR(:,1) RR RR(:,end)]; % pad RR with additional columnsRRE = [RRE(1,:); RRE; RRE(end,:)]; % and rowsfor iter = 1:nIterationsCompareWithNeighbor(-1, 0); % NorthCompareWithNeighbor(-1, 1); % North-EastCompareWithNeighbor(0, 1); % EastCompareWithNeighbor(1, 1); % South-EastCompareWithNeighbor(1, 0); % SouthCompareWithNeighbor(1, -1); % South-WestCompareWithNeighbor(0, -1); % WestCompareWithNeighbor(-1, -1); % North-WestendNP = 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 Imagenrows = 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
%因OPE边界处填充了0,故IP对应的边界处结果无意义,直接置成原值
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)); endNP = (OPE(2:(end-1),2:(end-1)) + IP)/2; % The Averaging operation
OPE(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 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);%下采样,将blocksize*blocksize区域映射为一个像素点
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 asfor crt_col = 1:result_cols % the average of each such blockResult(crt_row, crt_col) = mean2(A(1+(crt_row-1)*blocksize:crt_row*blocksize, ...1+(crt_col-1)*blocksize:crt_col*blocksize));end
end
输入
输出
注:输出时只是简单的进行线性拉伸,使得灰度值落在[0-255],没有使用更好调整方法。
参考:
http://www.cnblogs.com/sleepwalker/p/3676600.html
[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.
这篇关于Retinex系列之McCann99 Retinex的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!