带色彩恢复的视网膜增强算法实现 （MATLAB版本）

带色彩恢复的视网膜增强算法实现 （MATLAB版本）

%{
工程名称：带色彩恢复的视网膜增强算法实现
修改日期：2015年11月13日10:40:24
修改思路：http://www.cnblogs.com/Imageshop/archive/2013/04/17/3026881.html
测试结果：
下一步方向：
参考论文：
归纳总结：
  （1）S=L*R
  这个地方我认为S是原始图像，范围0-255.L是入射光，范围也是0-255
  而反射光是一个系数，有S和L的值决定。
  （2）这个算法那在水下图像处理中，有时候会产生颜色失真。
  （3）尺度系数选择20 100 300
  （4）这个算法产生的噪声比CLAHE小一点。。CLAHE的噪声有时候难以忍受
%}
clc
clear all;
close all
%第一部分： 程序部分 
img=imread('f:\\water1\\f6.jpg');% 修改成你的图像所在的地址
imgr=double(img(:,:,1));
imgg=double(img(:,:,2));
imgb=double(img(:,:,3));
mr=mat2gray(im2double(imgr)); 
mg=mat2gray(im2double(imgg)); 
mb=mat2gray(im2double(imgb));%数据类型归一化 
alf1=200; %定义标准差alf=a^2/2  a=20
n=351;%模板越大越好，161的时候好像效果不是很好
n1=floor((n+1)/2);%计算中心
for i=1:n
    for j=1:n
      b(i,j) =exp(-((i-n1)^2+(j-n1)^2)/(4*alf1))/(4*pi*alf1); %这个系数不知道有没有影响了
    end
end % 
sum1=sum(sum(b));%
b=b/sum1;% 归一化处理
nr1 = double(imfilter(imgr,b,'conv', 'replicate'));% 滤波器？？
ng1 = double(imfilter(imgg,b,'conv', 'replicate'));
nb1 = double(imfilter(imgb,b,'conv', 'replicate'));%卷积滤波
fil=cat(3,nr1,ng1,nb1);% 模糊结果
[h0,w0]=size(imgr); 
for i=1:h0
    for j=1:w0      
        % 通过循环进行控制
        if(imgr(i,j)==0)||(nr1(i,j)==0)
           % 这个地方透射率必然是0
           yr1(i,j)=0;
        else 
           yr1(i,j)=(log(imgr(i,j))-log(nr1(i,j)));% 这个地方计算的不准确啊。。。
        end        
        if(imgg(i,j)==0)||(ng1(i,j)==0)
           % 这个地方透射率必然是0
           yg1(i,j)=0;
        else 
           yg1(i,j)=(log(imgg(i,j))-log(ng1(i,j)));
        end        
        if(imgb(i,j)==0)||(nb1(i,j)==0)
           % 这个地方透射率必然是0
           yb1(i,j)=0;
        else 
           yb1(i,j)=(log(imgb(i,j))-log(nb1(i,j)));
        end   
        % 不知道什么地方出现了错误
    end
end
alf1=5000; %定义标准差alf=a^2/2  a=100
n=351;%
n1=floor((n+1)/2);%计算中心
for i=1:n
    for j=1:n
      b(i,j) =exp(-((i-n1)^2+(j-n1)^2)/(4*alf1))/(4*pi*alf1); %这个系数不知道有没有影响了
    end
end 
sum1=sum(sum(b));
b=b/sum1;% 
nr1 = double(imfilter(imgr,b,'conv', 'replicate'));% 滤波器？？
ng1 = double(imfilter(imgg,b,'conv', 'replicate'));
nb1 = double(imfilter(imgb,b,'conv', 'replicate'));%卷积滤波
fil=cat(3,nr1,ng1,nb1);% 这个地方进行模糊是没有问题的
[h0,w0]=size(imgr); 
for i=1:h0
    for j=1:w0      
        % 通过循环进行控制
        if(imgr(i,j)==0)||(nr1(i,j)==0)
           % 这个地方透射率必然是0
           yr2(i,j)=0;
        else 
           yr2(i,j)=(log(imgr(i,j))-log(nr1(i,j)));% 
        end        
        if(imgg(i,j)==0)||(ng1(i,j)==0)
           % 这个地方透射率必然是0
           yg2(i,j)=0;
        else 
           yg2(i,j)=(log(imgg(i,j))-log(ng1(i,j)));
        end        
        if(imgb(i,j)==0)||(nb1(i,j)==0)
           % 这个地方透射率必然是0
           yb2(i,j)=0;
        else 
           yb2(i,j)=(log(imgb(i,j))-log(nb1(i,j)));
        end   
    end
end
alf1=45000; %定义标准差alf=a^2/2  a=300
n=351;%模板的大小好像没什么大的影响
n1=floor((n+1)/2);%计算中心
for i=1:n
    for j=1:n
      b(i,j) =exp(-((i-n1)^2+(j-n1)^2)/(4*alf1))/(4*pi*alf1); %
    end
end % 
sum1=sum(sum(b));% 
b=b/sum1;% 归一化
nr1 = double(imfilter(imgr,b,'conv', 'replicate'));% 滤波器？？
ng1 = double(imfilter(imgg,b,'conv', 'replicate'));
nb1 = double(imfilter(imgb,b,'conv', 'replicate'));%卷积滤波
fil=cat(3,nr1,ng1,nb1);%
[h0,w0]=size(imgr); 
for i=1:h0
    for j=1:w0      
        % 通过循环进行控制
        if(imgr(i,j)==0)||(nr1(i,j)==0)
           % 这个地方透射率必然是0
           yr3(i,j)=0;
        else 
           yr3(i,j)=(log(imgr(i,j))-log(nr1(i,j)));% 这个地方计算的不准确啊。。。
           % 看来还是这个地方计算有问题
        end        
        if(imgg(i,j)==0)||(ng1(i,j)==0)
           % 这个地方透射率必然是0
           yg3(i,j)=0;
        else 
           yg3(i,j)=(log(imgg(i,j))-log(ng1(i,j)));
        end        
        if(imgb(i,j)==0)||(nb1(i,j)==0)
           % 这个地方透射率必然是0
           yb3(i,j)=0;
        else 
           yb3(i,j)=(log(imgb(i,j))-log(nb1(i,j)));
        end   
        % 不知道什么地方出现了错误
    end
end
imgout_r=(yr1+yr2+yr3)/3;% 如果不去拉伸，亮度会很低的
imgout_g=(yg1+yg2+yg3)/3;
imgout_b=(yb1+yb2+yb3)/3;

mean_r=mean2(imgout_r);% 对视网膜增强之后的图像进行拉伸处理
mean_g=mean2(imgout_g);
mean_b=mean2(imgout_b);
var_r=std2(imgout_r);
var_g=std2(imgout_g);
var_b=std2(imgout_b);
min_r=mean_r-2*var_r;  
max_r=mean_r+2*var_r;  
min_g=mean_g-2*var_g;  
max_g=mean_g+2*var_g; 
min_b=mean_b-2*var_b;  
max_b=mean_b+2*var_b;  
imgoutr=255*(imgout_r-min_r)/(max_r-min_r);
imgoutg=255*(imgout_g-min_g)/(max_g-min_g);
imgoutb=255*(imgout_b-min_b)/(max_b-min_b);
res_out=cat(3,imgoutr,imgoutg,imgoutb);% 实践证明，上面这个拉伸和直方图拉效果很像。。
figure
subplot(121),imshow(uint8(img)),title('原始输入图像');
subplot(122),imshow(uint8(res_out)),title('增强之后的图像');%
figure,
subplot(321),imhist(uint8(img(:,:,1))),title('原始图像直方图R');
subplot(322),imhist(uint8(res_out(:,:,1))),title('处理图像直方图R');
subplot(323),imhist(uint8(img(:,:,2))),title('原始图像直方图G');
subplot(324),imhist(uint8(res_out(:,:,2))),title('处理图像直方图G');
subplot(325),imhist(uint8(img(:,:,3))),title('原始图像直方图B');
subplot(326),imhist(uint8(res_out(:,:,3))),title('处理图像直方图B');