基于Retinex的图像去雾算法（MATLAB实现）

在【图像处理中的数学原理】专栏（该专栏中的文章已经结集出版，书名为《图像处理中的数学修炼》）之前的一些文章中，我们已经讨论了诸多非常有用的图像增强算法，例如直方图均衡算法以及更加强大的CLAHE。通常图像增强算法或多或少都有一定的去雾效果，只是这个效果有强有若罢了。本文将讨论另外一类十分重要的图像增强算法——Retinex算法。并在MATLAB中实验一下这类方法的去雾效果。

提及图像去雾算法，【图像处理中的数学原理】专栏中也已经讨论过两个非常有名的算法：基于暗通道先验的图像去雾算法和基于优化对比度增强的图像去雾算法，对图像去雾这一话题仍然感觉比较陌生的读者不妨先仔细读读上面两篇文章。

欢迎关注白马负金羁的博客 http://blog.csdn.net/baimafujinji，为保证公式、图表得以正确显示，强烈建议你从该地址上查看原版博文。本博客主要关注方向包括：数字图像处理、算法设计与分析、数据结构、机器学习、数据挖掘、统计分析方法、自然语言处理。

最基本的Retinex理论与算法

Retinex是一种常用的建立在科学实验和科学分析基础上的图像增强方法，它是以人类视觉系统为出发点发展而来的一套理论方法，最早由埃德温•兰德（Edwin. H. Land）于1963年提出。Retinex是由两个单词合成的一个词语，他们分别是retina 和cortex，即视网膜和皮层。

Retinex理论的基本内容是物体的颜色是由物体对长波（红色）、中波（绿色）、短波（蓝色）光线的反射能力来决定的，而不是由反射光强度的绝对值来决定的，物体的色彩不受光照非均匀性的影响，具有一致性，即Retinex是以色感一致性（颜色恒常性）为基础的。

根据兰德提出的理论，一幅给定的图像S(x,y)可以分解为两个不同的图像：反射图像R(x,y)和亮度图像（或称之为入射图像）L(x,y)，其原理如下图所示。


对于给定图像S中的每个点（x，y），用公式可以表示为：S(x,y)=R(x,y)⋅L(x,y)实际上，Retinex理论就是通过图像S来得到物体的反射性质R，也就是设法去除（或者降低）入射光L的影响从而得到物体原本该有的样子。但是具体该如何来估计R并没有一个明确的答案，因此根据不同的估计方法，也就产生了各种各样的Retinex算法。

单尺度的Retinex算法（SSR, Single Scale Retinex）是最基础、最简单的一种Retinex算法，而且这个算法也给出了广义上Retinex算法的大致框架。

Step1：利用取对数的方法将照射光分量和反射光分量分离，即logS(x,y)=logR(x,y)+logL(x,y)

Step2：一般我们会把最终的反射图像假设地估计为空间平滑图像（其物理解释就是通过计算图像中像素点与周围区域中像素的加权平均来对图像中照度变化做估计，并将其去除，最后只保留图像中物体的反射属性），所以可以用高斯模板对原图像作卷积，即相当于对原图像作低通滤波，得到低通滤波后的图像D(x,y)，F(x,y)表示高斯滤波函数：D(x,y)=S(x,y)∗F(x,y)

Step3：在对数域中，用原图像减去低通滤波后的图像，得到高频增强的图像G(x,y)：G(x,y)=logS(x,y)−logD(x,y)

Step4：对G(x,y)取反对数，得到增强后的图像R(x,y)：R(x,y)=expG(x,y)

基于SSR算法便可以实现一个基本的图像去雾程序。下面的MATLAB代码是完全按照上面的思路来实现的。只是在最后，对R(x,y)作对比度增强，以得到最终的去雾图像。此外，因为这里处理的是彩色图像，所以我们对R、G、B三个通道分别进行了处理。

I = imread('canon.jpg');

R = I(:, :, 1);
[N1, M1] = size(R);
R0 = double(R);
Rlog = log(R0+1);
Rfft2 = fft2(R0);

sigma = 250;
F = fspecial('gaussian', [N1,M1], sigma);
Efft = fft2(double(F));

DR0 = Rfft2.* Efft;
DR = ifft2(DR0);

DRlog = log(DR +1);
Rr = Rlog - DRlog;
EXPRr = exp(Rr);
MIN = min(min(EXPRr));
MAX = max(max(EXPRr));
EXPRr = (EXPRr - MIN)/(MAX - MIN);
EXPRr = adapthisteq(EXPRr);

G = I(:, :, 2);

G0 = double(G);
Glog = log(G0+1);
Gfft2 = fft2(G0);

DG0 = Gfft2.* Efft;
DG = ifft2(DG0);

DGlog = log(DG +1);
Gg = Glog - DGlog;
EXPGg = exp(Gg);
MIN = min(min(EXPGg));
MAX = max(max(EXPGg));
EXPGg = (EXPGg - MIN)/(MAX - MIN);
EXPGg = adapthisteq(EXPGg);

B = I(:, :, 3);

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

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

DBlog = log(DB+1);
Bb = Blog - DBlog;
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);
subplot(121), imshow(I);
subplot(122), imshow(result);

这里使用了两幅常用的有雾图像来做实验，其中左侧的图像为有雾图像，右侧图像是基于SSR实现的去雾后的效果图，可见SSR除了对图像进行了一定的增强之外，也有一定的去雾效果。





更多图像处理中的数学问题请参考《图像处理中的数学修炼》，也欢迎广大读者到图像处理书籍读者群中参与讨论学习，并交流关于图像增强算法和图像复原算法的研究心得。

多尺度Retinex算法与MSRCR

多尺度Retinex算法（MSR）是从SSR发展而来的一种Retinex算法，它的基本公式如下：R(x,y)=∑kKwk{logS(x,y)−logFk(x,y)∗S(x,y)}其中，R(x,y)是Retinex的输出（这与之前介绍SSR时的情况一致），i∈R,G,B表示三个颜色通道，F(x,y)是高斯滤波函数。wk表示尺度的权重因子，K表示使用的尺度数目。K=3 即表示彩色图像（此时通道i有红、绿、蓝三个）；K=1 表示灰度图像，即只有一个颜色通道。上式公式也揭示了MSR算法的特点，也就是在输出图像时，能够兼顾到色调再现和动态范围压缩两个特性。

在MSR算法的增强过程中，图像可能会因为增加了噪声而导致图像中局部区域的色彩失真，使得物体真正的颜色效果不能很好地显现出来，从而影响了整体的视觉观感。为了改进这方面的不足，一般情况下会使用带色彩恢复因子C的多尺度算法（MSRCR）来解决（Rehman等提出）。带色彩恢复因子C的多尺度算法是在多个固定尺度的基础上考虑色彩不失真恢复的结果，在多尺度Retinex算法过程中，通过引入一个色彩因子C来弥补由于图像局部区域对比度增强而导致的图像颜色失真的缺陷，通常情况下所引入的色彩恢复因子C的表达式为：RMSRCRi(x,y)=Ci(x,y)RMSRi(x,y)Ci(x,y)=f[Ii(x,y)∑Nj=1Ij(x,y)]其中，Ci表示第i个通道的色彩恢复系数，它的作用是调节3个通道颜色的比例，f(⋅)表示颜色空间的映射函数，通常可以采用下面的形式：Ci(x,y)=βlogαIi(x,y)∑Nj=1Ij(x,y)=β{log[αIi]−log[∑j=1NIj(x,y)]}其中，β is gain constant，α controls the strength of non-linearity。带色彩恢复的多尺度Retinex算法通过色彩恢复因子C这个系数来调整原始图像中3个颜色通道之间的比例关系，从而把相对有点暗的区域的信息凸显出来，以达到消除图像色彩失真缺陷的目的。处理后的图像局域对比度得以提高，而且其亮度与真实的场景很相似，图像在人们的视觉感知下显得更为逼真。因此，MSRCR算法会具有比较好的颜色再现性、亮度恒常性与动态范围压缩等特性。

下面就在MATLAB中来实现一下MSRCR，并验证一下它的去雾效果。

I = imread('toys.jpg');

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);

% B通道的处理方法与R和G类似，这里省略

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

从下面的结果来看，MSRCR比SSR的图像增强效果更佳，色彩也更逼真。当然，上面实现的MSRCR并没有为去雾做过多的特别设计，所以这也不完全算是一种很完善的去雾算法。如果实验更多的有雾图片，就会发现它的一些弱点，例如对天空部分的处理效果还不尽如人意，还有某些色彩比较深的地方会变得更暗而导致辨识度下降。有兴趣的读者也可以在这些地方加以改进，以期获得实用性更强的、普适性更强的图像去雾算法。

参考文献及其他推荐阅读材料

[1] 一篇介绍Retinex算法的博客文章，个人感觉应该算是比较详细的。

[2] 本文代码的最初实现来自《数字图像处理高级应用：基于MATLAB与CUDA的实现》，笔者有修改。