Retinex图像增强算法(SSR, MSR, MSRCR)详解及其OpenCV源码

Retinex是一种常用的建立在科学实验和科学分析基础上的图像增强方法，它是Edwin.H.Land于1963年提出的。就跟Matlab是由Matrix和Laboratory合成的一样，Retinex也是由两个单词合成的一个词语，他们分别是retina 和cortex，即：视网膜和皮层。Land的retinex模式是建立在以下三个假设之上的：

真实世界是无颜色的，我们所感知的颜色是光与物质的相互作用的结果。我们见到的水是无色的，但是水膜—肥皂膜却是显现五彩缤纷，那是薄膜表面光干涉的结果。
每一颜色区域由给定波长的红、绿、蓝三原色构成的；
三原色决定了每个单位区域的颜色。

Retinex理论的基础理论是物体的颜色是由物体对长波（红色）、中波（绿色）、短波（蓝色）光线的反射能力来决定的，而不是由反射光强度的绝对值来决定的，物体的色彩不受光照非均匀性的影响，具有一致性，即retinex是以色感一致性（颜色恒常性）为基础的。不同于传统的线性、非线性的只能增强图像某一类特征的方法，Retinex可以在动态范围压缩、边缘增强和颜色恒常三个方面达到平衡，因此可以对各种不同类型的图像进行自适应的增强。

40多年来，研究人员模仿人类视觉系统发展了Retinex算法，从单尺度Retinex算法，改进成多尺度加权平均的MSR算法，再发展成彩色恢复多尺度MSRCR算法。笔者在本文中，从原理阐述算法的流程，并提供源码地址。

参考链接：
http://blog.csdn.net/carson2005/article/details/9502053
http://www.cnblogs.com/Imageshop/archive/2013/04/17/3026881.html
参考论文：
《一种结合直方图均衡化和MSRCR的图像增强新算法》——李锦等人

一. 单尺度SSR(Single Scale Retinex)

1. 原理

一幅给定的图像S(x,y)可以分解为两个不同的图像：反射图像R(x,y)和入射图像（也有人称之为亮度图像）L(x,y)，其原理图如下所示：
这里写图片描述
如上图所示，图像可以看做是入射图像和反射图像构成，入射光照射在反射物体上，通过反射物体的反射，形成反射光进入人眼。最后形成的图像可以如下公式表示：
$r(x, y) = logR(x, y) = log \frac{S(x, y)}{L(x, y)}$

其中，R(x, y)表示了物体的反射性质，即图像内在属性，我们应该最大程度的保留；而L(x, y)表示入射光图像，决定了图像像素能达到的动态范围，我们应该尽量去除。
一般，我们把照射图像假设估计为空间平滑图像，原始图像为S(x, y)，反射图像为R(x, y)，亮度图像为L(x, y)，可以得出上面的公式(1)，以及下面的公式：

$r(x, y) = logR(x, y) = log \frac{S(x, y)}{L(x, y)}$ …………(2)

$r(x, y) = \log S(x, y) - \log [F(x, y) \bigotimes S(x, y) ]$ ………….(3)

这里，r(x, y)是输出图像，式(3)中后面中括号里的运算是卷积运算。F(x, y)是中心环绕函数，表示为：

$F(x, y) = \lambda e^{\frac{-(x^2 + y^2)}{c^{2}}}$ ………….(4)

式(4)中的C是高斯环绕尺度，λ是一个尺度，它的取值必须满足下式：

$\int \int F(x, y)dxdy = 1$ ………….(5)

上面的式中可以看出，SSR算法中的卷积是对入射图像的计算，其物理意义是通过计算像素点与周围区域在加权平均的作用下，估计图像中照度的变化，并将 $L(x, y)$ 去除，只保留 $S(x, y)$ 属性。

2. 算法的实现流程

单尺度Retinex算法SSR的实现流程可以概括如下：

读原图S(x, y)：
- 若原图为灰度图：将图像各像素的灰度值由整数型(int)转换为浮点数(float)，并转换到对数域；
- 若原图为彩色图：将颜色分通道处理，每个分量像素值由整数型(int)转换为浮点数(float)，并转换到对数域；
输入高斯环绕尺度C，把积分运算离散化，转为求和运算，通过上式(4)(5)确定λ的值；
由式(3)得r(x, y)；
- 若原图是灰度图，则只有一个r(x, y)；
- 若原图为彩色图，则每个通道都有一个对应的r(x, y)；
将r(x, y)从对数域转换到实数域，得到输出图像R(x, y)；
此时的R(x, y)值的范围并不是0–255，所以还需要进行线性拉伸并转换成相应的格式输出显示。

前面的公式中，中心环绕函数F(x, y)用的是低通函数，这样能够在算法中估计出入射图像对应原始图像的低频部分。从原始图像中除去低频照射部分，就会留下原始图像所对应的高频分量。高频分量很有价值，因为在人类的视觉系统中，人眼对边缘部分的高频信息相当敏感，所以SSR算法可以较好的增强图像中的边缘信息。
由于SSR算法中所选用的高斯函数特点，对于动态范围大幅度压缩和对比度增强两个指标，增强后的图像不能同时保证。但是为了平衡两种增强效果，就必须选择一个较为恰当的高斯尺度常量C。C值一般取值在80–100之间。

3. OpenCV源码

见源码附录的Retenix()函数。

二. 多尺度MSR(Multi-Scale Retinex)

1. 原理

MSR是在SSR基础上发展来的，优点是可以同时保持图像高保真度与对图像的动态范围进行压缩的同时，MSR也可实现色彩增强、颜色恒常性、局部动态范围压缩、全局动态范围压缩，也可以用于X光图像增强。
MSR计算公式如下：

$r(x, y) = \sum_{k}^{K}w_{k} \{ \log S(x, y) - \log [F_{k}(x, y) \cdot S(x, y) ] \}$ …………(6)

式中，K是高斯中心环绕函数的个数。当K=1时，MSR退化为SSR。
通常来讲，为了保证兼有SSR高、中、低三个尺度的优点来考虑，K取值通常为3，且有：

$w_{1} = w_{2} = w_{3} = \frac{1}{3}$ …………(7)

此外，实验表明，ci分别取15, 80, 200可以得到较好效果。

一般的Retinex算法对光照图像估计时，都会假设初始光照图像是缓慢变化的，即光照图像是平滑的。但实际并非如此，亮度相差很大区域的边缘处，图像光照变化并不平滑。所以在这种情况下，Retinuex增强算法在亮度差异大区域的增强图像会产生光晕。
另外MSR常见的缺点还有边缘锐化不足，阴影边界突兀，部分颜色发生扭曲，纹理不清晰，高光区域细节没有得到明显改善，对高光区域敏感度小等。

2. 效果对比

对于SSR与MSR的效果如下所示：

原图

SSR（尺度为300）增强

MSR（尺度为300，尺度数为3）增强

原图

SSR（尺度为300）增强

MSR（最大尺度为300，尺度数为3）增强

3. OpenCV源码

见源码附录的MultiScaleRetinex()函数。

三. 带颜色恢复的MSR方法MSRCR(Multi-Scale Retinex with Color Restoration)

在以上的两幅测试图像中，特别是第二幅，我们看到明显的偏色效果，这就是SSR和MSR普遍都存在的问题。

为此，研究者又开发出一种称之为带色彩恢复的多尺度视网膜增强算法（MSRCR,Multi-Scale Retinex with Color Restoration) ，具体讨论的过程详见《A Multiscale Retinex for Bridging the Gap Between Color Images and the Human Observation of Scenes》论文。

1. 原理

在前面的增强过程中，图像可能会因为增加了噪声，而使得图像的局部细节色彩失真，不能显现出物体的真正颜色，整体视觉效果变差。针对这一点不足，MSRCR在MSR的基础上，加入了色彩恢复因子C来调节由于图像局部区域对比度增强而导致颜色失真的缺陷。
改进算法如下所示：

$R_{MSRCR_{i}}(x, y) = C_{i}(x, y) R_{MSR_{i}}(x, y)$ …………(８)
$C_{i}(x, y) = f[I_{i}^{'}(x, y)] = f[\frac{I_{i}(x, y)}{ \sum_{j=1}^{N} I_{j}(x, y)}]$ …………(９)
$f[I_{i}^{'}(x, y)] = \beta \log [\alpha I_{i}^{'}(x, y)] = \beta \{ \log [\alpha I_{i}^{'}(x, y)] - \log [\sum_{j=1}^{N}I_{j}(x, y)] \}$ …………(10)
其中参数说明如下：

Ii(x, y)表示第i个通道的图像
Ci表示第i个通道的彩色回复因子，用来调节3个通道颜色的比例；
f(·)表示颜色空间的映射函数；
β是增益常数；
α是受控制的非线性强度；

MSRCR算法利用彩色恢复因子C，调节原始图像中3个颜色通道之间的比例关系，从而把相对较暗区域的信息凸显出来，达到了消除图像色彩失真的缺陷。
处理后的图像局部对比度提高，亮度与真实场景相似，在人们视觉感知下，图像显得更加逼真。

但是MSRCR算法处理图像后，像素值一般会出现负值。所以从对数域r(x, y)转换为实数域R(x, y)后，需要通过改变增益Gain，偏差Offset对图像进行修正。使用公式可以表示为：
$R_{MSRCR_{i}}(x, y)' = G\cdot R_{MSRCR_{i}}(x, y) + O$ …………(11)
式(11)中，G表示增益Gain，O表示偏差Offset。它们的值取决于软件中的算法实现。

2. MSRCR其他实现方法

MSRCR方法如上，且获得了不错的效果。但在博文《MSRCR》中，作者认为论文里的方法不起任何作用，并且论文里为了这个又引入了太多的可调参数，增加了算法的复杂性，不利于自动化实现。
博文作者认为，GIMP的contrast-retinex.c文件里使用的算法很好，效果也很好。他直接从量化的方式上入手，引入了均值和均方差的概念，再加上一个控制图像动态的参数来实现无色偏的调节过程，简要描述如下：

计算出 log[R(x,y)]中R/G/B各通道数据的均值Mean和均方差Var（注意是均方差）。
类似下述公式计算各通道的Min和Max值。
- Min = Mean - Dynamic * Var;
- Max = Mean + Dynamic * Var;
对Log[R(x,y)]的每一个值Value，进行线性映射：
R(x,y) = ( Value - Min ) / (Max - Min) * (255 - 0), 同时要注意增加一个溢出判断，即：
- if (R(x, y) > 255) R(x,y) = 255;
- else if (R(x,y) < 0) R(x,y) = 0;

就是经过这么简单的处理，实践证明可以取得非常好的效果，下面贴出一些处理后的效果。

MSR（尺度为300，尺度数为3）增强

MSRCR（尺度为300,尺度数3,Dynamic = 2）增强

MSRCR（尺度为300,尺度数6,Dynamic = 2）增强

MSR（尺度300，尺度数为3）增强
MSRCR（尺度300,尺度为3,Dynamic = 2）增强图像
MSRCR（尺度300, 尺度为3, Dynamic = 2）增强图像
MSRCR（尺度300,尺度为6,Dynamic = 2）增强图像
MSRCR（尺度300, 尺度为6, Dynamic = 2）增强图像

由以上三幅图的效果得出的结论：

MSRCR效果要比MSR好很多，基本消除了色偏。
对于MSRCR，尺度数对结果的影像不是特别大，但是随着尺度数的增加，算法耗时会线性增加，因此，一般尺度数取3就较为合适了。

继续做比较：

MSRCR（Dynamic = 1）增强图像

MSRCR（Dynamic = 2）增强图像

MSRCR（Dynamic = 5）增强图像

MSRCR（Dynamic = 1）增强图像

MSRCR（Dynamic = 2）增强图像

MSRCR（Dynamic = 5）增强图像

由以上三幅图的效果得出的结论：

Dynamic取值越小，图像的对比度越强。
一般来说Dynamic取值2-3之间能取得较为明显的增强效果，即能取得很自然过渡效果，又能保持图像的清晰度适度增强。

关于尺度，个人建议取值以大于100为佳。

retinex算法的效果对于一些正常的图像处理后的效果并不佳，我们可以认为它就是为那些在外界环境不理想的状态下拍摄的图像增强而设计的，如对于航拍的雾天图片，医学上的成像图片等成像条件恶劣的图有很明显的效果。笔者的毕业设计涉及了水下图像的修复，水下图像相对于空气中成像性低，应该会有很大的效果。

NASA在处理航天相关照片时也使用了Retinex技术。关于NASA对Retinex技术的应用，可以参考：http://dragon.larc.nasa.gov/retinex/

关于GIMP的实现代码，笔者提供从网上找到的链接。如果要想完全看懂代码的意思很难，但算法部分还是不难提取的（不过笔者也懒得把它移植到OpenCV上……）下载链接如下；
http://files.cnblogs.com/Imageshop/contrast-retinex.rar

四. Retinex图像增强算法的其他算法

先贴出前面MSRCR的基本公式如下：

R M S R C R i (x, y) = 255 R M S R C R i ( x , y ) - m i n i ( m i n ( x , y ) R M S R C R i ( x , y ) ) m a x i ( m a x ( x , y ) R M S R C R i ( x , y ) - m i n i ( m i n ( x , y ) R M S R C R ( x , y ) ) )

$R_{MSRCR_{i}}(x, y) = 255 \frac{R_{MSRCR_{i}}(x, y) - min_{i}(min_{(x, y)} R_{MSRCR_{i}}(x, y))}{max_{i}(max_{(x, y)} R_{MSRCR_{i}}(x, y) - min_{i}(min{(x, y)} R_{MSRCR_{(x, y)}}))}$
…………(12)

由于Retinex数据处理后的高动态特性，所以这种方式处理后，数据分布很广，会出现严重两极化现象，一般情况难以得到满意的结果。

第二种方法：经典Canonical Gain/Offset算法：
即在经典MSRCR文章《A Multiscale Retinex for Bridging the Gap Between Color Images and the Human Observation of Scenes》论文中提到的算法。公式如下：
$R_{MSRCR_{i}}(x, y) = G [R_{MSRCR_{i}}(x, y) - b]$ …………(13)
式(12)中的G与b均为经验参数。