MATLAB代码Retinex算法解决光照不均问题

Retinex算法解决光照不均问题的MATLAB代码实现，可直接运行出结果，有注释。 （另外，代码最后的注释部分为利用顶帽变换方法处理光照不均问题。在本人处理的实际问题中，顶帽变换算法不理想，但Retinex算法的处理结果很好。）

%
% Retinex 补偿光照不均匀
% 由于卷积时需要保证精度，应先把[0255]的unit8图像归一化成[0.01.0]的浮点型。
% 那么利用matlab实现步骤为：读入图像→归一化→设置高斯函数参数及矩阵→高斯函数和输入图像矩阵卷积→取对数→和输入图像矩阵的对数相差→取指数→输出图像
name = [‘g-4-1.png‘];
f = imread（name）;  

% f=imread（‘H25.jpg‘）;
fr=f（: : 1）; fg=f（: : 2）; fb=f（: : 3）;%RGB通道
mr=mat2gray（im2double（fr））; mg=mat2gray（im2double（fg））; mb=mat2gray（im2double（fb））;%数据类型归一化
alf1=1458; %定义标准差alf=a^2/2  a=54
n=161;%定义模板大小 
n1=floor（（n+1）/2）;%计算中心 
for i=1:n 
    for j=1:n 
      b（ij） =exp（-（（i-n1）^2+（j-n1）^2）/（4*alf1））/（pi*alf1）; %高斯函数
    end 
end 
nr1 = imfilter（mrb‘conv‘ ‘replicate‘）;ng1 = imfilter（mgb‘conv‘ ‘replicate‘）;nb1 = imfilter（mbb‘conv‘ ‘replicate‘）;%卷积滤波
ur1=log（nr1）; ug1=log（ng1）; ub1=log（nb1）;
tr1=log（mr）;tg1=log（mg）;tb1=log（mb）;
yr1=（tr1-ur1）/3;yg1=（tg1-ug1）/3;yb1=（tb1-ub1）/3;
alf2=53.38; %定义标准差alf=a^2/2    a=10.3325
x=31;%定义模板大小 
x1=floor（（n+1）/2）;%计算中心 
for i=1:n 
    for j=1:n 
      a（ij） =exp（-（（i-n1）^2+（j-n1）^2）/（4*alf2））/（6*pi*alf2）; %高斯函数
    end 
end 
nr2 = imfilter（mra‘conv‘ ‘replicate‘）;ng2 = imfilter（mga‘conv‘ ‘replicate‘）;nb2 = imfilter（mba‘conv‘ ‘replicate‘）;%卷积滤波
ur2=log（nr2）; ug2=log（ng2）; ub2=log（nb2）;
tr2=log（mr）;tg2=log（mg）;tb2=log（mb）;
yr2=（tr2-ur2）/3;yg2=（tg2-ug2）/3;yb2=（tb2-ub2）/3;
alf3=13944.5; %定义标准差alf=a^2/2  a=167
l=501;%定义模板大小 
l1=floor（（n+1）/2）;%计算中心 
for i=1:n 
    for j=1:n 
      e（ij） =exp（-（（i-n1）^2+（j-n1）^2）/（4*alf3））/（4*pi*alf3）; %高斯函数
    end 
end 
nr3 = imfilter（mre‘conv‘ ‘replicate‘）;ng3 = imfilter（mge‘conv‘ ‘replicate‘）;nb3 = imfilter（mbe‘conv‘ ‘replicate‘）;%卷积滤波
ur3=log（nr3）; ug3=log（ng3）; ub3=log（nb3）;
tr3=log（mr）;tg3=log（mg）;tb3=log（mb）;
yr3=（tr3-ur3）/3;yg3=（tg3-ug3）/3;yb3=（tb3-ub3）/3;
dr=yr1+yr2+yr3;dg=yg1+yg2+yg3;db=yb1+yb2+yb3;
    
cr=im2uint8（dr）; cg=im2uint8（dg）; cb=im2uint8（db）;
z=cat（3 cr cg cb）;          
figure（1） imshow（z）
figure（2） imshow（f）
ibw = im2bw（z graythresh（z））;
figure（3）; imshow（ibw）;

I2 = imclose（ibwones（55））;
figure（4）; imshow（I2）;
I2 = imopen（I2ones（55））;
figure（5）; imshow（I2）;

bw = im2bw（f graythresh（f））;
subplot（221）;imshow（f）;
% title（‘原图像‘）;
subplot（223）;imshow（z）;
% title（‘增强图像‘）;
subplot（222）;imshow（bw）;
% title（‘原图像二值化后的图像‘）;
subplot（224）;imshow（I2）;
% title（‘增强图像二值化后的图像‘）;
imwrite（z‘增强后.png‘）;
%}

%{
% 测试顶帽变换
for i = 1:1
%    II = imread（[‘C:\Documents and Settings\Administrator\桌面\tophat\‘ int2str（i） ‘.png‘]）; 
   name = [‘C:\Documents and Settings\Administrator\桌面\tophat\g-4-1.png‘];
   II = imread（name）;  
   I = rgb2gray（II）;
%    name = [‘C:\Documents and Settings\Administrator\桌面\tophat\g-‘ int2str（i） ‘.png‘];
%    imwrite（II name）;
%    subplot（141）;
   figure（1）;
   imshow（I）;
   
   thresh = graythresh（I）;
   ibw = im2bw（I thresh）;
%    subplot（142）; 
   figure（2）;
   imshow（ibw）;

   bg = imopen（I strel（‘disk‘ 50））;
   figure（6）;imshow（bg）;
   Itophat = imsubtract（I bg）;
%    subplot（143）; 
   figure（3）;
   im

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件     165374  2016-08-31 13:50  g-4-1.PNG

     文件       3479  2016-10-11 09:23  Retinex.m

----------- ---------  ---------- -----  ----

               168853                    2