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算法之前,我先来讲一点个人的理解,以便第一次接触该理论的朋友能够更快速地理解。当然,如果我的理解有问题,也请大家帮忙指出。
在极端情况下,我们大可以认为整幅图像中的分量都是均匀的,那么最简单的估计照度L的方式就是在将图像变换到对数域后对整幅图像求均值。因此,我设计了以下算法来验证自己的猜想,流程如下:
这里为了简化描述,省略了对图像本身格式的变换,算法用Matlab实现:
![](/assets/blank.gif)
% ImOriginal:原始图像 % type:'add'表示分量是加性的,如雾天图像;'mult'表示分量是乘性的,如对照度的估计 [m,n,z] = size(ImOriginal); ImOut = uint8(zeros(m,n,z)); for i = 1:zif strcmp(type,'add')ImChannel = double(ImOriginal(:,:,i))+eps; elseif strcmp(type,'mult')ImChannel = log(double(ImOriginal(:,:,i))+eps); elseerror('type must be ''add'' or ''mult'''); endImOut(:,:,i) = EnhanceOneChannel(ImChannel); end ImOut = max(min(ImOut,255), 0); endfunction ImOut = EnhanceOneChannel(ImChannel) % 计算计算单个通道的反射分量 % 1.对全图进行照射分量估计 % 2.减去照射分量 % 3.灰度拉伸 ImChannel = ImChannel./max(ImChannel(:)); ImRetinex = round(exp(ImChannel.*5.54)); ImOut = uint8(ImRetinex); end
![](/assets/blank.gif)
为了验证算法的有效性,这里使用经典的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算法
1. 将原图像变换到对数域,彩色图像对各通道都进行对数变换;
3. 计算路径,如上图8所示,这里令为路径上的点,从远到近排列;
公式所表示的大致意思为:从远到近,中心点像素值减去路径上的像素值得到的差值的一半与前一时刻的估计值之间的和。最终,中心像素点的像素大致的形式为
图9.算法的测试图(来自Retinex in Matlab)
![](/assets/blank.gif)
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:nIterationsCompareWith(0, shift); % horizontal stepCompareWith(shift, 0); % vertical step endshift = -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)); elseIP(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
![](/assets/blank.gif)
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. 当最底层计算完毕得到的即最终增强后的图像。
为了方便各位自己研究,下面我给出该算法源码供大家参考:
![](/assets/blank.gif)
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
![](/assets/blank.gif)
总结
在上述两个经典的估计算法之后经过反对数变换就可以恢复成原格式的图像。通过试验可以发现,这些算法都还是有一定缺陷的,对于照度的估计后续还有很多算法如:单尺度Retinex (SSR)、多尺度Retinex (MSR)、可变框架的Retinex等等;还有针对增强后的光晕现象先用进行照度分割,然后在再计算等等方法,本质上都大同小异。这些改进算法这里不再进行介绍,有兴趣的朋友可以去下载些论文看看。有不明白的或者本文有错误的地方也希望各位能够指出,以免误导后面的读者。
以上为参考的文章,代码如下:
主程序
clear all;
rgb = imread('2.jpg');%需要处理的图片
m=size(rgb,1);
n=size(rgb,2);
% 初始化矩阵
rr = zeros(m,n);
gg = zeros(m,n);
bb = zeros(m,n);
% eps的目的是防止分母出现0的时候产生正负无穷
% 如果反射率分量是乘性的,那需要变换到对数域
rr = log(double(rgb(:,:,1))+eps);
gg = log(double(rgb(:,:,2))+eps);
bb = log(double(rgb(:,:,3))+eps);
% 如果反射率分量是加性的,那只需要归一化就可以
% rr = double(rgb(:,:,1))+eps;
% gg = double(rgb(:,:,2))+eps;
% bb = double(rgb(:,:,3))+eps;
% 归一化
rr=rr/max(rr(:));
gg=gg/max(gg(:));
bb=bb/max(bb(:));
% Retinex算法
rrr = retinex_frankle_mccann(rr, 8);
ggg = retinex_frankle_mccann(gg, 8);
bbb = retinex_frankle_mccann(bb, 8);
% e^5.54约等于255,这里Retinex算法输出的是指数分量,每个像素亮度与该点计算结果x的关系为:I=255^x
rrr = round(exp(rrr.*5.54));
ggg = round(exp(ggg.*5.54));
bbb = round(exp(bbb.*5.54));
RGB = cat(3,uint8(rrr),uint8(ggg),uint8(bbb));
RGB = max(min(RGB,255),0);
subplot(2,2,1);imshow(rgb);title('原图像');
subplot(2,2,2),imhist(rgb2gray(rgb));title('原图像的直方图');
subplot(2,2,3);imshow(RGB);title('新图像');
subplot(2,2,4),imhist(rgb2gray(RGB));title('新图像的直方图');
调用函数:
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
算法的主体思想是对照度的估计,mccan的算法思想是从外围的像素点作为估计初始值,不断的向内收敛,对当前估计值和选取点的像素值做平均后作为下一点的估计值。最终得到的估计值便是照度的估计值。
shift为距离中心的偏移,逐渐的向内收缩。
Retinex去照度算法相关推荐
- [论文阅读] (11)ACE算法和暗通道先验图像去雾算法(Rizzi | 何恺明老师)
<娜璋带你读论文>系列主要是督促自己阅读优秀论文及听取学术讲座,并分享给大家,希望您喜欢.由于作者的英文水平和学术能力不高,需要不断提升,所以还请大家批评指正,非常欢迎大家给我留言评论,学 ...
- 单幅图像去雾算法研究综述
来源 <计算机工程与应用>北大核心期刊,CSCD数据库. 影响因子:2.348 简介 图像去雾算法是以满足特定场景需求,突出图片细节并增强图片质量为目的的图像分析与处理方法.在雾霾天气下, ...
- 图像和视频的快速去雾算法研究
王昕, 孙莹莹, 李影昉. 图像和视频的快速去雾算法研究[J]. 影像科学与光化学, 2016, 34(1): 82-87. WANG Xin, SUN Yingying, LI Yingfang ...
- c++ opencv编程实现暗通道图像去雾算法_OpenCV图像处理专栏十五 |一种基于亮度均衡的图像阈值分割技术...
前言 对于光照不均匀的图像,用通常的图像分割方法不能取得满意的效果.为了解决这个问题,论文<一种基于亮度均衡的图像阈值分割技术>提出了一种实用而简便的图像分割方法.该方法针对图像中不同亮度 ...
- c++ opencv编程实现暗通道图像去雾算法_OpenCV图像处理专栏十三 | 利用多尺度融合提升图像细节...
前言 今天为大家介绍一个利用多尺度来提升图像细节的算法.这个算法来自于论文<DARK IMAGE ENHANCEMENT BASED ON PAIRWISE TARGET CONTRAST AN ...
- 图像去雾算法(二)基于暗通道先验算法学习笔记
在http://write.blog.csdn.net/postedit/78301999中介绍了图像去雾的相关研究方法,发现目前为止在图像去雾方面,何凯明博士基于暗通道先验的算法具有很好的效果,关于 ...
- 综述:视频和图像去雾算法以及相关的图像恢复和增强研究
综述:视频和图像去雾算法以及相关的图像恢复和增强研究 翻译自IEEE的一篇文章<Review of Video and Image Defogging Algorithms and Relate ...
- Python基于OpenCV的图像去雾算法[完整源码&部署教程]
1.图片识别 2.视频展示 [项目分享]Python基于OpenCV的图像去雾算法[完整源码&部署教程]_哔哩哔哩_bilibili 3.算法原理 图像增强算法常见于对图像的亮度.对比度.饱和 ...
- 第一篇 Frankle-Mccan去雾算法
因为对图像处理感兴趣,就去本部蹭了丁老师的<数字图像处理>课程,感觉比自学效率高多了: 业余时间打算把课上部分代码实现,于是就有了这个博客:这里会不定期更新一些基础算法,大都是老师课上讲的 ...
最新文章
- c# 一些控件常用屬性
- 10个机器学习的JavaScript示例
- strcmp可以比较数组么_005 继承、封装、多态及数组初识
- python 打包文件
- Spring-Cloud 学习笔记-(4)负载均衡器Ribbon
- python基础教程(十一)
- 分治应用--万里挑一 找假硬币
- python程序文件是什么_.py文件是什么?
- 【TensorFlow】TensorFlow函数精讲之tf.constant()
- 配置tomcat用户
- 提高MySQL性能的方法
- JSPs only permit GET POST or HEAD的解决方案(REST风格)
- Springboot初始化过程(1.5.9.RELEASE)(一)
- 有一种VR电影比爱情动作片更“爽”
- Arcgis字段计算器实现自动顺序编号
- GEE拼接字符串出错,原因是忘了加getInfo()
- 百度HI QQ和MSN 阿里旺旺贸易通MSN在线客服聊天代码
- 在Windows Server 2008上用Windows Media Service打造流媒体直播系统
- Python效率之王之多进程和多线程详解
- 博士申请 | 香港城市大学计算机学院徐伟涛老师组招收人工智能全奖博士生