小时候就觉得人的大脑是个极其之神奇的东西，然后高中学了生物化学那些东西后，觉得脑子居然可以产生思想，感情，宇宙中居然有这样的东西的存在，然后就感慨不已，上了大学，搞了一些算法学术后，意识到一些以前并没有怎么意识到的一些大脑的功能【虽说不管怎么说其神奇程度都不及会产生情感，思想。。】，比如说人大脑的模糊算法，我们看一局棋，可以快速的判断出谁处于优势谁处于劣势，有经验的人瞄一眼天空，就知道是不是要马上下雨了。“大概怎么样”这个概念，对应于属于“模糊”，人的大脑先天就有这种神奇的功能，对于计算机这种连真正的随机数都产生不了的“渣渣”而言，人脑实在过于强大【我一直觉得，计算机之所以会有存在的价值，就在于他可以高速的做一些简单的运算，如果人脑可以快速的在1纳秒内算出根号π精确到小数点后100位的话，那计算机对我们现实的价值可能会少很多。。】。再比如说学习功能，现在多少所谓的机器学习算法其实被发明的起源都是来源于对大脑的研究。

好吧，其实我扯了好远。。真正跟本文有关的事，大脑的知觉恒常性问题【这个问题应该不能归属到模糊逻辑里面去】。什么是知觉恒常性，比如说大小恒常性，一个物体，放在离你很近的地方，或者离你很远的地方，你都可以大概知道这个物体有多大。这个就是我们大脑的大小恒常机制带来的。计算机就完全办不到了，应该。。。办不到吧。。。远处的东西虽然看上去变小了，但是你大脑还是会对得到的图像进行缩放，因为你可以从眼睛得到的图像里面得出物体的深度信息。而计算机图片就是一张二维的照片，没有深度的，即便如此，我们看照片也可以得到深度信息，而计算机却不行。。。吧。。。

知觉恒常性还包括很多，方向恒常，形状恒常等。。

好吧，其实我还是扯远了，本文要讲的是另一个知觉恒常问题——颜色恒常性！！一个绿色的叶子，你在室外，你看的出来它是绿色的，然后你走进一间屋子，打开一盏红光灯，你依旧知道它是绿色的，即便肉眼看上去不是那么的“纯绿”，但是你的大脑知道，那是一片绿色的叶子。计算机呢？用摄像头拍下红光下的绿叶，你看着照片，你也可以判断出它是绿叶。计算机自己呢，他会计算出叶子对应位置的R,G,B三种颜色的像素，然后通过一个表，告诉你这个RGB值对应什么颜色，显然不会再是绿色。

那么计算机颜色恒常性问题就是要帮助计算机解决这么一个问题！

好吧，这个问题现在解决方法多到要死，下图是以前做的一个给老师的报告里面的一个图的一个小部分。。。

我这里想讲的是。。。Retinex算法，算法介绍什么的自己谷歌，其实就是很多很多年前的几个人根据一些仿生的东西退出来一些理论，然后现在发现可以运用于颜色恒常性。Retinex是retina与cortex的合成缩写。

虽说要搞颜色恒常性这个问题，嗯，水很深，但是要弄懂最基本的Retinex算法，却是一件很简单的事情。。。

先讲一下它的基本思想，Retinex是基于这样一个世界观（果然我动画看多了么。。），他认为，我们获得的图像可以表示成下面的形式：

S(x,y) = R(x,y)×L(x,y)

其中R代表入射光，就是上面例子中的红光，L代表反射性质，相当于背景原本应该具有的颜色，即例子中的叶子，S代表反射光，也就是我们获得的图像，我们知道S，但是不知道R和L，我们只要知道R或者L中的一个，就可以知道另一个了~

乘法计算起来不方便，于是我们变到对数域上去：

log(S(x,y)) = log(R(x,y))+log(L(x,y))

Retinex的思路就是：获得S，估算出R，得出物体本身的颜色L~~~~~

怎么估算出L，我们要先研究一下R和L的数据的一些特性，好吧，我直接告诉你好了，R是低频信号，L是高频信号。

其实这很好理解，我们红光照在一间屋子里，得到屋子的图像，可以分解成一副纯红光的图像和一副屋子在白光下的图片，纯红光图片对应频谱是啥，肯定是低频啦，与此相比，背景就是高频信号。

不一定是有色光下的图片，还可以是夜间图片，相当于“黑光”照耀下的图片。。

这样我们就知道要怎么算啦，对L进行低通滤波，这样就得到R，然后对数域下S减去R，就得到L了，嗯，没错，这就是Retinex！！

一般我们看到的论文文献里面采用的滤波器都是高斯滤波器。【高斯是Retinex算法的官配吧。。】

高斯滤波器模板就是：

G(x,y) = λ exp[-(x²+y²)/σ²]  

其中σ是尺度参数，越大锐化越厉害；λ归一化常数，存在的价值就是使得下式成立:

∫∫G(x,y)dxdy = 1

然后Retinex里面常用的三个模型，SSR,MSR,MSRCR；SSR是单尺度Retinex(Single Scale Rtinex),MSR是多尺度的(Multi Scale Retinex)，MSRCR是带彩色恢复的多尺度Retinex(MSR with color restore)。

上面讲的最基本的算法呢，就是SSR了，算法就是：

1. 用高斯模板和图像卷积，给图像滤波，得到滤波后的图像
2. 将图像和滤波后的图像都Log一下，转到对数域
3. 相减
4. 做Exp指数运算，恢复到正常的图像来。

简单吧，简单到你不敢相信吧~

好吧，算法呢，就是这样，但是你自己编程的时候才会发现种种问题，比如Exp回来之后数据不在0~255之间要怎么办？？高斯模板的参数要怎么选取云云。自己查看相关论文吧【或者参考本文最后的code】。

上面就是SSR了，MSR，多尺度是怎么回事呢？其实也很简单，就是弄多几个不同参数的高斯模板，不同参数导致锐化程度不一样，然后每个模板相当于做SSR，得到的几个结果加权相加，一般都去权重相同。

算法不用我重复吧，还是不懂的，读下文的code去吧。。

再说MSRCR了，我其实不知道这个算法的发明人是怎么想出来的，或者我对其理解不完全，而且我以前实验的结果也不知道和MSR比优势在哪。。。好吧，我没做太深入的研究，这方面。。

算法呢，是这样的，好吧，还是他妹的要敲公式。。。

先回到MSR，假设我们得到了MSR处理的结果了，设为I_MSR，这里面有R,G,B三层数据，假设每一层分别是I_{MSR_i}，i=1,2,3分别表示R,G,B三层的数据，每个点就可以表示成I_{MSR_i}(x,y)了，然后我们的MSRCR的结果是怎么算的呢？

I_{MSRCR_i}(x,y) = C_i(x,y)×I_{MSR_i}(x,y)       i=1,2,3

什么意思呢，就是说每一层的每个点都要系数的修正，系数的修正是怎么算的呢，这个有很多算法了，但是基本思想都是考虑在(x,y)这个点上R,G,B像素的比值。好吧，说不清，看公式。。

 C_i(x,y) = log[α I_i(x,y)] – log[ ∑I_i(x,y) ]

其中I_i(x,y)是原始图像中(x,y)这个点的R,G,B像素值，这里涉及到一个参数α，自己看论文选吧。

为什么要这么做？为了满足灰色世界理论，就是说令修复后的图像R,G,B三个的像素值趋于1:1:1。【灰色世界理论的适用范围本文不讨论。】

 C_i(x,y) 还有一些别的计算公式，比如：

 C_i(x,y) =   {  log[α I_i(x,y)] – log[ ∑I_i(x,y) ]  + b} × G

其实就是第一个 C_i(x,y) 的结果加b之后乘G，很无聊。。是吧。。

嘛，算法就讲到这里，下面先给副结果图，一个比较暗的环境下的照片，经过三种算法处理的结果。

结果表示，三种算法，差别不大。。【不要吐槽图片的画质。。。因为滤波过程计算复杂度较大，耗时较长，所以选了一副低像素的图片】

这里只想再补充一句，Retinex算法的一个很麻烦的问题，就是处理后会有白膜，上面的图还看不出来。学术界有很多对付白膜的研究，想深入研究的人再去慢慢深入研究吧~

下面是代码。。

Code:

close all;
clc;
I = imread('1.jpg');
Ir=I(:,:,1);
Ig=I(:,:,2);
Ib=I(:,:,3); 

%%%%%%%%%%设定所需参数%%%%%% 
G = 192;
b = -30;
alpha = 125;
beta = 46;
Ir_double=double(Ir);
Ig_double=double(Ig);
Ib_double=double(Ib); 

%%%%%%%%%%设定高斯参数%%%%%% 
sigma_1=15; %三个高斯核 
sigma_2=80;
sigma_3=250;
[x y]=meshgrid((-(size(Ir,2)-1)/2):(size(Ir,2)/2),(-(size(Ir,1)-1)/2):(size(Ir,1)/2));
gauss_1=exp(-(x.^2+y.^2)/(2*sigma_1*sigma_1)); %计算高斯函数 
Gauss_1=gauss_1/sum(gauss_1(:)); %归一化处理 
gauss_2=exp(-(x.^2+y.^2)/(2*sigma_2*sigma_2));
Gauss_2=gauss_2/sum(gauss_2(:));
gauss_3=exp(-(x.^2+y.^2)/(2*sigma_3*sigma_3));
Gauss_3=gauss_3/sum(gauss_3(:)); 

%%%%%%%%%%对R分量操作%%%%%%% 
% MSR部分 
Ir_log=log(Ir_double+1); %将图像转换到对数域 
f_Ir=fft2(Ir_double); %对图像进行傅立叶变换,转换到频域中 

%sigam=15的处理结果 
fgauss=fft2(Gauss_1,size(Ir,1),size(Ir,2));
fgauss=fftshift(fgauss); %将频域中心移到零点 
Rr=ifft2(fgauss.*f_Ir); %做卷积后变换回空域中 
min1=min(min(Rr));
Rr_log= log(Rr - min1+1);
Rr1=Ir_log-Rr_log; 

%sigam=80 
fgauss=fft2(Gauss_2,size(Ir,1),size(Ir,2));
fgauss=fftshift(fgauss);
Rr= ifft2(fgauss.*f_Ir);
min1=min(min(Rr));
Rr_log= log(Rr - min1+1);
Rr2=Ir_log-Rr_log; 

 %sigam=250 
fgauss=fft2(Gauss_3,size(Ir,1),size(Ir,2));
fgauss=fftshift(fgauss);
Rr= ifft2(fgauss.*f_Ir);
min1=min(min(Rr));
Rr_log= log(Rr - min1+1);
Rr3=Ir_log-Rr_log; 

Rr=0.33*Rr1+0.34*Rr2+0.33*Rr3; %加权求和 
MSR1 = Rr;
SSR1 = Rr2;
%计算CR 
CRr = beta*(log(alpha*Ir_double+1)-log(Ir_double+Ig_double+Ib_double+1)); 

%SSR:
min1 = min(min(SSR1));
max1 = max(max(SSR1));
SSR1 = uint8(255*(SSR1-min1)/(max1-min1)); 

%MSR
min1 = min(min(MSR1));
max1 = max(max(MSR1));
MSR1 = uint8(255*(MSR1-min1)/(max1-min1)); 

%MSRCR 
Rr = G*(CRr.*Rr+b);
min1 = min(min(Rr));
max1 = max(max(Rr));
Rr_final = uint8(255*(Rr-min1)/(max1-min1)); 

%%%%%%%%%%对g分量操作%%%%%%% 
Ig_double=double(Ig);
Ig_log=log(Ig_double+1); %将图像转换到对数域 
f_Ig=fft2(Ig_double); %对图像进行傅立叶变换,转换到频域中 

fgauss=fft2(Gauss_1,size(Ig,1),size(Ig,2));
fgauss=fftshift(fgauss); %将频域中心移到零点 
Rg= ifft2(fgauss.*f_Ig); %做卷积后变换回空域中 
min2=min(min(Rg));
Rg_log= log(Rg-min2+1);
Rg1=Ig_log-Rg_log; %sigam=15的处理结果 

fgauss=fft2(Gauss_2,size(Ig,1),size(Ig,2));
fgauss=fftshift(fgauss);
Rg= ifft2(fgauss.*f_Ig);
min2=min(min(Rg));
Rg_log= log(Rg-min2+1);
Rg2=Ig_log-Rg_log; %sigam=80 

fgauss=fft2(Gauss_3,size(Ig,1),size(Ig,2));
fgauss=fftshift(fgauss);
Rg= ifft2(fgauss.*f_Ig);
min2=min(min(Rg));
Rg_log= log(Rg-min2+1);
Rg3=Ig_log-Rg_log; %sigam=250 

Rg=0.33*Rg1+0.34*Rg2+0.33*Rg3; %加权求和 
SSR2 = Rg2;
MSR2 = Rg;
%计算CR 
CRg = beta*(log(alpha*Ig_double+1)-log(Ir_double+Ig_double+Ib_double+1)); 

%SSR:
min2 = min(min(SSR2));
max2 = max(max(SSR2));
SSR2 = uint8(255*(SSR2-min2)/(max2-min2)); 

%MSR
min2 = min(min(MSR2));
max2 = max(max(MSR2));
MSR2 = uint8(255*(MSR2-min2)/(max2-min2)); 

%MSRCR 
Rg = G*(CRg.*Rg+b);
min2 = min(min(Rg));
max2 = max(max(Rg));
Rg_final = uint8(255*(Rg-min2)/(max2-min2)); 

%%%%%%%%%%对B分量操作同R分量%%%%%%% 
Ib_double=double(Ib);
Ib_log=log(Ib_double+1);
f_Ib=fft2(Ib_double); 

fgauss=fft2(Gauss_1,size(Ib,1),size(Ib,2));
fgauss=fftshift(fgauss);
Rb= ifft2(fgauss.*f_Ib);
min3=min(min(Rb));
Rb_log= log(Rb-min3+1);
Rb1=Ib_log-Rb_log; 

fgauss=fft2(Gauss_2,size(Ib,1),size(Ib,2));
fgauss=fftshift(fgauss);
Rb= ifft2(fgauss.*f_Ib);
min3=min(min(Rb));
Rb_log= log(Rb-min3+1);
Rb2=Ib_log-Rb_log; 

fgauss=fft2(Gauss_3,size(Ib,1),size(Ib,2));
fgauss=fftshift(fgauss);
Rb= ifft2(fgauss.*f_Ib);
min3=min(min(Rb));
Rb_log= log(Rb-min3+1);
Rb3=Ib_log-Rb_log; 

Rb=0.33*Rb1+0.34*Rb2+0.33*Rb3; 

%计算CR 
CRb = beta*(log(alpha*Ib_double+1)-log(Ir_double+Ig_double+Ib_double+1));
SSR3 = Rb2;
MSR3 = Rb;
%SSR:
min3 = min(min(SSR3));
max3 = max(max(SSR3));
SSR3 = uint8(255*(SSR3-min3)/(max3-min3));

%MSR
min3 = min(min(MSR3));
max3 = max(max(MSR3));
MSR3 = uint8(255*(MSR3-min3)/(max3-min3));

%MSRCR
Rb = G*(CRb.*Rb+b);
min3 = min(min(Rb));
max3 = max(max(Rb));
Rb_final = uint8(255*(Rb-min3)/(max3-min3)); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ssr = cat(3,SSR1,SSR2,SSR3);
msr = cat(3,MSR1,MSR2,MSR3);
msrcr=cat(3,Rr_final,Rg_final,Rb_final); %将三通道图像合并 

subplot(2,2,1);imshow(I);title('原图') %显示原始图像 
subplot(2,2,2);imshow(ssr);title('SSR')
subplot(2,2,3);imshow(msr);title('MSR')
subplot(2,2,4);imshow(msrcr);title('MSRCR') %显示处理后的图像

【完】

本文内容遵从CC版权协议，转载请注明出自http://www.kylen314.com