Retinex 图像去模糊（含MATLAB代码）

Retinex 图像去模糊

Retinex 理论：
Retinex是由 Retina + Cortex 组成，被称为视网膜皮层理论。 该理论指出 物体能够被观察到的信息是由两个因素来决定的：物体本身的反射性质和物体周围的光照强度。其中，光照强度决定了原始图像中所有像素点的动态范围大小，而原始图像的固有属性（颜色）是由物体本身的反射系数决定的。 因此该方法是将一幅图像分成两个不同的图像：反射图像和亮度图像，然后去除光照的影响，保留物体的固有属性。

该方法目前主要有三个分支：单尺度Retinex（SSR），多尺度Retinex（MSR）和带颜色恢复的多尺度Retinex（MSRCR）。
Retinex理论的基本思路：在原始图像中，通过某种方法去除或者降低入射图像的影响，从而保留物体本质的反射属性图像。

单尺度Retinex（SSR）：

算法思路：
1 读取原始图像，分别取出RGB三个通道的数值矩阵，转化成 double 型。
2 利用高斯模板对各个通道的进行卷积操作处理。
3 将RGB图像以及卷积后的图像转化到对数域，并相减，然后利用反对数转化到 实数域。
4 对获得的各个通道的图像进行线性拉伸，合并RGB就可的处理后的图像。
算法实现（MATLAB）：

close all; clear all; clc
I = imread('E:\水下目标识别\chinamm2019uw_train\jpegimage\000000.jpg');
I_r = double(I(:,:,1));
I_g = double(I(:,:,2));
I_b = double(I(:,:,3));

I_r_log = log(I_r+1);
I_g_log = log(I_g+1);
I_b_log = log(I_b+1);

Rfft1 = fft2(I_r);
Gfft1 = fft2(I_g);
Bfft1 = fft2(I_b);

%  SSR算法
[m,n] = size(I_r);
sigma = 200;
f = fspecial('gaussian', [m, n], sigma);
efft1 = fft2(double(f));

D_r = ifft2(Rfft1.*efft1);
D_g = ifft2(Gfft1.*efft1);
D_b = ifft2(Bfft1.*efft1);

D_r_log = log(D_r + 1);
D_g_log = log(D_g + 1);
D_b_log = log(D_b + 1);

R = I_r_log - D_r_log;
G = I_g_log - D_g_log;
B = I_b_log - D_b_log;

R = exp(R);
MIN = min(min(R)); 
MAX = max(max(R));
R = (R - MIN)/(MAX - MIN);
R = adapthisteq(R);
G = exp(G);
MIN = min(min(G)); 
MAX = max(max(G));
G = (G - MIN)/(MAX - MIN);
G = adapthisteq(G);
B = exp(B);
MIN = min(min(B)); 
MAX = max(max(B));
B = (B - MIN)/(MAX - MIN);
B = adapthisteq(B);

J = cat(3, R, G, B);

figure;
subplot(121);imshow(I);
subplot(122);imshow(J);
figure;imshow(J)

结果：

多尺度Retinex（MSR）：

由于SSR需要在颜色保真度和细节保持上追求一个完美的平衡，而这个平衡不宜实现。MSR的提出就是解决该问题，它是针对一幅图像在不同尺度上利用高斯滤波处理，然后将多个尺度下的图像进行加权叠加。
它其实就是 SSR 的一中叠加操作。
代码：

close all; clear all; clc
I = imread('E:\水下目标识别\chinamm2019uw_train\jpegimage\000000.jpg');
I_r = double(I(:,:,1));
I_g = double(I(:,:,2));
I_b = double(I(:,:,3));

I_r_log = log(I_r+1);
I_g_log = log(I_g+1);
I_b_log = log(I_b+1);

Rfft1 = fft2(I_r);
Gfft1 = fft2(I_g);
Bfft1 = fft2(I_b);[m,n] = size(I_r);
sigma1 = 15;
sigma2 = 80;
sigma3 = 200;
f1 = fspecial('gaussian', [m, n], sigma1);
f2 = fspecial('gaussian', [m, n], sigma2);
f3 = fspecial('gaussian', [m, n], sigma3);
efft1 = fft2(double(f1));
efft2 = fft2(double(f2));
efft3 = fft2(double(f3));

D_r1 = ifft2(Rfft1.*efft1);
D_g1 = ifft2(Gfft1.*efft1);
D_b1 = ifft2(Bfft1.*efft1);
D_r_log1 = log(D_r1 + 1);
D_g_log1 = log(D_g1 + 1);
D_b_log1 = log(D_b1 + 1);
R1 = I_r_log - D_r_log1;
G1 = I_g_log - D_g_log1;
B1 = I_b_log - D_b_log1;

D_r2 = ifft2(Rfft1.*efft2);
D_g2 = ifft2(Gfft1.*efft2);
D_b2 = ifft2(Bfft1.*efft2);
D_r_log2 = log(D_r2 + 1);
D_g_log2 = log(D_g2 + 1);
D_b_log2 = log(D_b2 + 1);
R2 = I_r_log - D_r_log2;
G2 = I_g_log - D_g_log2;
B2 = I_b_log - D_b_log2;

D_r3 = ifft2(Rfft1.*efft3);
D_g3 = ifft2(Gfft1.*efft3);
D_b3 = ifft2(Bfft1.*efft3);
D_r_log3 = log(D_r3 + 1);
D_g_log3 = log(D_g3 + 1);
D_b_log3 = log(D_b3 + 1);
R3 = I_r_log - D_r_log3;
G3 = I_g_log - D_g_log3;
B3 = I_b_log - D_b_log3;

R = 0.1*R1 + 0.4*R2 + 0.5*R3;
G = 0.1*G1 + 0.4*G2 + 0.5*G3;
B = 0.1*B1 + 0.4*B2 + 0.5*B3;

R = exp(R);
MIN = min(min(R)); 
MAX = max(max(R));
R = (R - MIN)/(MAX - MIN);
R = adapthisteq(R);
G = exp(G);
MIN = min(min(G)); 
MAX = max(max(G));
G = (G - MIN)/(MAX - MIN);
G = adapthisteq(G);
B = exp(B);
MIN = min(min(B)); 
MAX = max(max(B));
B = (B - MIN)/(MAX - MIN);
B = adapthisteq(B);

J = cat(3, R, G, B);

figure;
subplot(121);imshow(I);
subplot(122);imshow(J,[]);

figure;imshow(J)

结果：

带颜色恢复的多尺度Retinex（MSRCR）：

由于在采用SSR和MSR的过程会产生色差，MSRCR将MSR得到的结果按照一定的比例进行调整以求恢复原来的比例数值，具体是通过引入颜色恢复因子C，来弥补由于图像局部区域对比度增强而导致的图像颜色失真的缺陷。
代码：

% close all; clear all; clc
% I = imread('E:\水下目标识别\chinamm2019uw_train\jpegimage\000000.jpg');
% figure;imshow(I)
% I_r = I(:,:,1);
% I_g = I(:,:,2);
% I_b = I(:,:,3);
% figure;
% subplot(131);imshow(I_r);
% subplot(132);imshow(I_g);
% subplot(133);imshow(I_b);
% 
% I_r = double(I_r);
% I_g = double(I_g);
% I_b = double(I_b);
% 
% [m, n] = size(I_r);
% Rlog = log(I_r+1);
% Rfft2 = fft2(I_r);
% 
% Glog = log(I_g+1);
% Gfft2 = fft2(I_g);
% 
% Blog = log(I_b+1);
% Bfft2 = fft2(I_b);
% 
% sigma1 = 128;
% f = fspecial('gaussian', [m,n], sigma1);
% Efft1 = fft2(double(f));
% 
% DI_r = Rfft2.* Efft1;
% DR = ifft2(DI_r);
% DRlog = log(DR+1);
% Rr1 = Rlog - DRlog;
% figure;imshow(Rr1);
% 
% DI_g = Gfft2.* Efft1;
% DG = ifft2(DI_g);
% DGlog = log(DG+1);
% Gg1 = Glog - DGlog;
% figure;imshow(Gg1);
% 
% DI_b = Bfft2.* Efft1;
% DB = ifft2(DI_b);
% DBlog = log(DB+1);
% Bb1 = Blog - DBlog;
% figure;imshow(Bb1);
% 
% Rr1 = Rr1 .* 255/max(max(Rr1));
% Gg1 = Gg1 .* 255/max(max(Gg1));
% Bb1 = Bb1 .* 255/max(max(Bb1));
% 
% II(:,:,1) = Rr1;
% II(:,:,2) = Gg1;
% II(:,:,3) = Bb1;
% figure;
% subplot(121);imshow(I);
% subplot(122);imshow(II)

I = imread('E:\水下目标识别\chinamm2019uw_train\jpegimage\000000.jpg');
figure;imshow(I)

% 提取RGB
R = I(:, :, 1);
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 = 125;
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);

Blog = log(B0+1); 
Bfft2 = fft2(B0); 

DB0 = Bfft2.* Efft1; 
DB = ifft2(DB0); 

DBlog = log(DB +1); 
Bb1 = Blog - DBlog; 

DB0 = Gfft2.* Efft2; 
DB = ifft2(DB0); 
DBlog = log(DB +1); 
Bb2 = Blog - DBlog;
DB0 = Gfft2.* Efft3; 
DB = ifft2(DB0);
DBlog = log(DB +1); 
Bb3 = Blog - DBlog;
Bb = (Bb1 + Bb2 + Bb3)/3; 

Ib = immultiply(B0, a);
C = imdivide(Ib, II); 
C = log(C+1);
Bb = immultiply(C, Bb); 
EXPBb= exp(Bb);
MIN = min(min(EXPBb)); 
MAX = max(max(EXPBb));
EXPBb = (EXPBb - MIN)/(MAX - MIN);
EXPBb = adapthisteq(EXPBb);

result = cat(3, EXPRr, EXPGg, EXPBb);  
figure;
subplot(121), imshow(I);
subplot(122), imshow(result);
figure;imshow(result)

结果：

上述是Retinex的一些列算法实现。