单尺度Retinex多尺度Retinex（matlab）

说明：此代码是针对灰度图像的
单尺度

clc;clear;
I=imread('DR_ck_0001.tiff');
 
% Ir=I(:,:,1);%提取红色分量
% Ig=I(:,:,2);%提取绿色分量
% Ib=I(:,:,3);%提取蓝色分量
% I1=SSR(Ir); %自定义单尺度Retinex处理函数
% I2=SSR(Ig);
% I3=SSR(Ib);
% In=cat(3,I1,I2,I3);          %cat用于构造多维数组
In = SSR(I);
t = exp(-10);
enhanced =  log(double(I)+t)-log(In+t);
enhanced = normal(enhanced);
subplot(1,3,1);imshow(I);
title('A). 原始图像');  
subplot(1,3,2);imshow(In,[]);
title('B). 单尺度Retinex');
subplot(1,3,3);imshow(enhanced,[]);
title('C). 单尺度Retinex');
 
%——————————————划重点！！！！！！！！！！！！！！！！！！！！————————————————
 
function G = SSR(I)
%G=I;
[m,n]=size(I);
I=double(I);
x0=floor((m+1)/2);y0=floor((n+1)/2);
c = 70;
K=1/(sqrt(2*pi)*c);
for x=1:m
    for y=1:n
            G(x,y)=K*(exp(-((x)^2+(y)^2)/(2*c^2))); %高斯函数
    end
end
 
Gfft=fft2(G);
Ifft=fft2(I);
Lfft=Gfft.*Ifft;
L=ifft2(Lfft);
 
r=log(I+1)-log(L+1);                   %加1是为了防止对数为0时，log0没有定义
 
%线性拉伸，将灰度范围转换到0-255，（直接求反对数的效果不好）
MIN=min(min(r));
MAX=max(max(r));
new=((r-MIN)*255/(MAX-MIN));
A=new;
end
function out =normal(img)
maxx = max(max(img));
minn = min(min(img));
out=(img-minn)./(maxx-minn);
end

多尺度

clc;clear;
I = imread('DR_ck_0001.tiff');
 
R = I(:, :);
% G = I(:, :, 2);
% B = I(:, :, 3);
R0 = double(R);
% G0 = double(G);
% B0 = double(B);
 
[N1, M1] = size(R);
 
Rlog = log(R0+1);
Rfft2 = fft2(R0);
 
sigma1 = 128;
F1 = fspecial('gaussian', [N1,M1], sigma1);
Efft1 = fft2(double(F1));
 
DR0 = Rfft2.* Efft1;
DR = ifft2(DR0);
 
DRlog = log(DR +1);
Rr1 = Rlog - DRlog;
 
sigma2 = 256;
F2 = fspecial('gaussian', [N1,M1], sigma2);
Efft2 = fft2(double(F2));
 
DR0 = Rfft2.* Efft2;
DR = ifft2(DR0);
 
DRlog = log(DR +1);
Rr2 = Rlog - DRlog;
 
sigma3 = 512;
F3 = fspecial('gaussian', [N1,M1], sigma3);
Efft3 = fft2(double(F3));
 
DR0 = Rfft2.* Efft3;
DR = ifft2(DR0);
 
DRlog = log(DR +1);
Rr3 = Rlog - DRlog;
 
Rr = (Rr1 + Rr2 +Rr3)/3;
 
% a = 1;
% II = R0;
% %II = imadd(R0, G0);
% %II = imadd(II, B0);
% Ir = immultiply(R0, a);
% C = imdivide(Ir, II);
% C = log(C+1);
%  
% Rr = immultiply(C, Rr);
EXPRr = exp(Rr);
MIN = min(min(EXPRr));
MAX = max(max(EXPRr));
EXPRr = (EXPRr - MIN)/(MAX - MIN);
%EXPRr = adapthisteq(EXPRr);
 
% Glog = log(G0+1);
% Gfft2 = fft2(G0);
%  
% DG0 = Gfft2.* Efft1;
% DG = ifft2(DG0);
%  
% DGlog = log(DG +1);
% Gg1 = Glog - DGlog;
%  
%  
% DG0 = Gfft2.* Efft2;
% DG = ifft2(DG0);
%  
% DGlog = log(DG +1);
% Gg2 = Glog - DGlog;
%  
%  
% DG0 = Gfft2.* Efft3;
% DG = ifft2(DG0);
%  
% DGlog = log(DG +1);
% Gg3 = Glog - DGlog;
%  
% Gg = (Gg1 + Gg2 +Gg3)/3;
%  
% Ig = immultiply(G0, a);
% C = imdivide(Ig, II);
% C = log(C+1);
%  
% Gg = immultiply(C, Gg);
% EXPGg = exp(Gg);
% MIN = min(min(EXPGg));
% MAX = max(max(EXPGg));
% EXPGg = (EXPGg - MIN)/(MAX - MIN);
% EXPGg = adapthisteq(EXPGg);
 
% B通道的处理方法与R和G类似，这里省略
 
% result = cat(3, EXPRr, EXPGg, EXPBb);  
result =  uint8(255.*EXPRr);
DRlog = normal(DRlog);
imwrite(result,'result2.png');
imwrite(uint8(255.*DRlog),'result3.png');
subplot(131), imshow(I);
subplot(132), imshow(result);
subplot(133), imshow(DRlog);
function out =normal(img)
maxx = max(max(img));
minn = min(min(img));
out=(img-minn)./(maxx-minn);
end