msr颜色复原matlab程序,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的一些列算法实现。