带色彩恢复的视网膜增强算法实现 （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的噪声有时候难以忍受%}clcclear 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=20n=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); %这个系数不知道有没有影响了    endend % 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           % 不知道什么地方出现了错误    endendalf1=5000; %定义标准差alf=a^2/2  a=100n=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); %这个系数不知道有没有影响了    endend 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       endendalf1=45000; %定义标准差alf=a^2/2  a=300n=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); %    endend % 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           % 不知道什么地方出现了错误    endendimgout_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);% 实践证明，上面这个拉伸和直方图拉效果很像。。figuresubplot(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');