MSRCR（Multi-Scale Retinex with Color Restore）多尺度Retinex图像增强

引言

始于Edwin Herbert Land（埃德温·赫伯特·兰德）于1971年提出的一种被称为色彩恒常的理论，并基于此理论的图像增强方法。Retinex这个词由视网膜(Retina)和大脑皮层(Cortex)合成而来.之所以这样设计，表明Land他也不清楚视觉系统的特性究竟取决于此两个生理结构中的哪一个，抑或两者都有关系。不同于传统的图像增强算法，如线性、非线性变换、图像锐化等只能增强图像的某一类特征，如压缩图像的动态范围，或增强图像的边缘等，Retinex可以在动态范围压缩、边缘增强和颜色恒常三方面达到平衡，可以对各种不同类型的图像进行自适应性地增强，在很多方面得到了广泛的应用。

算法概述

Retinex 理论的基本内容是物体的颜色是由物体对长波（红）、中波（绿）和短波（蓝）光线的反射能力决定的，而不是由反射光强度的绝对值决定的；物体的色彩不受光照非均性的影响，具有一致性，即Retinex理论是以色感一致性（颜色恒常性）为基础的。如下图所示，观察者所看到的物体的图像S是由物体表面对入射光L反射得到的，反射率R由物体本身决定，不受入射光L变化。 

对于观察图像S，有公式表示为： 

S(x,y)=R(x,y)L(x,y)(1)
其中， L(x,y) 表示亮度分量， R(x,y) 表示物体反射分量， S(x,y)表示观测到的图像 。 
Retinex理论增强算法的思想就是利用式(1)，去除亮度分量 L 求得反射分量 R ，从而达到图像增强效果。 
对(1)式两边取对数 
log(S(x,y))=log(L(x,y))+log(R(x,y))(2)
由(2)得 
log(R(x,y))=log(S(x,y))−log(L(x,y))(3)
由式(3)可知，只需要估计亮度分量 L 就能求得反射分量，因此L的估计直接决定图像去雾恢复效果.Jobson等论证了高斯卷积函数可以从已知图像 S 中更好地估计出亮度分量，即 
L(x,y)=S(x,y)∗G(x,y)(4)
其中，’*’代表卷积操作，高斯函数 G(x,y)=k⋅exp(−x2+y2σ2) ， σ 是高斯函数尺度参数， k 为归一化因子，使 ∬G(x,y)dxdy=1 。对于RGB彩色图像的某一个通道， 
log(Rc(x,y))=log(Sc(x,y))−log(Sc(x,y))∗G(x,y)(5)
其中， c 表示RGB3个通道的某一个通道； σ  为尺度函数，取值大时，颜色失真小但细节恢复差，取值小时，细节恢复好，但颜色失真大。 
为弥补SSR算法不足，MSR算法对每一个通道进行3次不同尺度滤波，加权求和，处理时间也会变长，且使用不同的尺度，导致恢复的RGB比值与原图比值不太一样，颜色失真。 
log(R(x,y))=Weight1⋅log(Rσ1(x,y))+Weight2⋅log(Rσ2(x,y))+Weight3⋅log(Rσ3(x,y))(6)
Rehman等提出MSRCR算法，引入分量比值调整因子从而降低色彩失真的影响。 
RMSRCR(x,y)=C(x,y)RMSR(x,y)(7)
其中 C(x,y)=G{log[αIc(x,y)]−log[∑c=13Ic(x,y)]} 。

流程

1. 按照（4）式，计算入射分量 L ；
2. 按照（3）式，计算单一尺度下的 log(Rσ(x,y)) ；
3. 按照（6）式，对不同尺度进行加权求和；
4. 按照（7）式，进行色彩恢复；
5. 将 RMSRCR 量化为0到255范围，输出。

入射分量获取

Retinex算法核心在于入射分量的获取，简单来说，这个可以对原图像进行高斯卷积获取，卷积核的大小与原图像大小相等。当图像比较大的时候，其计算过程是非常复杂的， O(M2N2) ， M ， N 是原图像的高度和宽度。实现的时候可以使用递归高斯滤波，复杂度 O(MN) （参见paper：Recursive implementation of the Gaussian filter，源码在GIMP中有）。下面我给出在3种不同尺度下获取的亮度图，等分权重。 

色彩恢复

GIMP中色彩恢复代码如下：

<code class="hljs cpp has-numbering" style="display: block; padding: 0px; background-color: transparent; color: inherit; box-sizing: border-box; font-family: 'Source Code Pro', monospace;font-size:undefined; white-space: pre; border-top-left-radius: 0px; border-top-right-radius: 0px; border-bottom-right-radius: 0px; border-bottom-left-radius: 0px; word-wrap: normal; background-position: initial initial; background-repeat: initial initial;">    <span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">/*  
    Final calculation with original value and cumulated filter values.  
    The parameters gain, alpha and offset are constants.  
    */</span>   
    <span class="hljs-comment" style="color: rgb(136, 0, 0); box-sizing: border-box;">/* Ci(x,y)=log[a Ii(x,y)]-log[ Ei=1-s Ii(x,y)] */</span>   

    alpha  = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">128.0f</span>;   
    gain   = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1.0f</span>;   
    offset = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0.0f</span>;   

    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">for</span> ( i = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0</span>; i < size; i += bytes )   
    {   
        <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">float</span> logl;   

        psrc = src+i;   
        pdst = dst+i;   

        logl = (<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">float</span>)<span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">log</span>( (<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">float</span>)psrc[<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0</span>] + (<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">float</span>)psrc[<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>] + (<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">float</span>)psrc[<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">2</span>] + <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">3.0f</span> );   

        pdst[<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0</span>] = gain * ((<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">float</span>)(<span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">log</span>(alpha * (psrc[<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0</span>]+<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1.0f</span>)) - logl) * pdst[<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0</span>]) + offset;   
        pdst[<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>] = gain * ((<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">float</span>)(<span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">log</span>(alpha * (psrc[<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>]+<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1.0f</span>)) - logl) * pdst[<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1</span>]) + offset;   
        pdst[<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">2</span>] = gain * ((<span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">float</span>)(<span class="hljs-built_in" style="color: rgb(102, 0, 102); box-sizing: border-box;">log</span>(alpha * (psrc[<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">2</span>]+<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1.0f</span>)) - logl) * pdst[<span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">2</span>]) + offset;   
    }</code><ul class="pre-numbering" style="box-sizing: border-box; position: absolute; width: 50px; background-color: rgb(238, 238, 238); top: 0px; left: 0px; margin: 0px; padding: 6px 0px 40px; border-right-width: 1px; border-right-style: solid; border-right-color: rgb(221, 221, 221); list-style: none; text-align: right;"><li style="box-sizing: border-box; padding: 0px 5px;">1</li><li style="box-sizing: border-box; padding: 0px 5px;">2</li><li style="box-sizing: border-box; padding: 0px 5px;">3</li><li style="box-sizing: border-box; padding: 0px 5px;">4</li><li style="box-sizing: border-box; padding: 0px 5px;">5</li><li style="box-sizing: border-box; padding: 0px 5px;">6</li><li style="box-sizing: border-box; padding: 0px 5px;">7</li><li style="box-sizing: border-box; padding: 0px 5px;">8</li><li style="box-sizing: border-box; padding: 0px 5px;">9</li><li style="box-sizing: border-box; padding: 0px 5px;">10</li><li style="box-sizing: border-box; padding: 0px 5px;">11</li><li style="box-sizing: border-box; padding: 0px 5px;">12</li><li style="box-sizing: border-box; padding: 0px 5px;">13</li><li style="box-sizing: border-box; padding: 0px 5px;">14</li><li style="box-sizing: border-box; padding: 0px 5px;">15</li><li style="box-sizing: border-box; padding: 0px 5px;">16</li><li style="box-sizing: border-box; padding: 0px 5px;">17</li><li style="box-sizing: border-box; padding: 0px 5px;">18</li><li style="box-sizing: border-box; padding: 0px 5px;">19</li><li style="box-sizing: border-box; padding: 0px 5px;">20</li><li style="box-sizing: border-box; padding: 0px 5px;">21</li><li style="box-sizing: border-box; padding: 0px 5px;">22</li><li style="box-sizing: border-box; padding: 0px 5px;">23</li></ul>

另外，最后一步代码又增加了调节项，用于控制饱和度和颜色，最后将其量化到[0~255]，详细代码如下。

<code class="hljs applescript has-numbering" style="display: block; padding: 0px; background-color: transparent; color: inherit; box-sizing: border-box; font-family: 'Source Code Pro', monospace;font-size:undefined; white-space: pre; border-top-left-radius: 0px; border-top-right-radius: 0px; border-bottom-right-radius: 0px; border-bottom-left-radius: 0px; word-wrap: normal; background-position: initial initial; background-repeat: initial initial;">    /*  
    Adapt <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">the</span> dynamics <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">of</span> <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">the</span> colors according <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">to</span> <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">the</span> statistics <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">of</span> <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">the</span> <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">first</span> <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">and</span> <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">second</span> order.  
    The use <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">of</span> <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">the</span> variance makes <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">it</span> possible <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">to</span> control <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">the</span> degree <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">of</span> saturation <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">of</span> <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">the</span> colors.  
    */   
    pdst = dst;

    compute_mean_var( pdst, &mean, &var, size, bytes );   
    mini = mean - rvals.cvar*var;   
    maxi = mean + rvals.cvar*var;   
    range = maxi - mini;   

    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">if</span> ( !range ) range = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">1.0</span>;   

    <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">for</span> ( i = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0</span>; i < size; i+= bytes )   
    {   
        psrc = src + i;   
        pdst = dst + i;   

        <span class="hljs-keyword" style="color: rgb(0, 0, 136); box-sizing: border-box;">for</span> (j = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0</span> ; j < <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">3</span> ; j++)   
        {   
            float c = <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">255</span> * ( pdst[j] - mini ) / range;   

            psrc[j] = (unsigned char)clip( c, <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">0</span>, <span class="hljs-number" style="color: rgb(0, 102, 102); box-sizing: border-box;">255</span> );   
        }   
    }   

    free (dst);</code><ul class="pre-numbering" style="box-sizing: border-box; position: absolute; width: 50px; background-color: rgb(238, 238, 238); top: 0px; left: 0px; margin: 0px; padding: 6px 0px 40px; border-right-width: 1px; border-right-style: solid; border-right-color: rgb(221, 221, 221); list-style: none; text-align: right;"><li style="box-sizing: border-box; padding: 0px 5px;">1</li><li style="box-sizing: border-box; padding: 0px 5px;">2</li><li style="box-sizing: border-box; padding: 0px 5px;">3</li><li style="box-sizing: border-box; padding: 0px 5px;">4</li><li style="box-sizing: border-box; padding: 0px 5px;">5</li><li style="box-sizing: border-box; padding: 0px 5px;">6</li><li style="box-sizing: border-box; padding: 0px 5px;">7</li><li style="box-sizing: border-box; padding: 0px 5px;">8</li><li style="box-sizing: border-box; padding: 0px 5px;">9</li><li style="box-sizing: border-box; padding: 0px 5px;">10</li><li style="box-sizing: border-box; padding: 0px 5px;">11</li><li style="box-sizing: border-box; padding: 0px 5px;">12</li><li style="box-sizing: border-box; padding: 0px 5px;">13</li><li style="box-sizing: border-box; padding: 0px 5px;">14</li><li style="box-sizing: border-box; padding: 0px 5px;">15</li><li style="box-sizing: border-box; padding: 0px 5px;">16</li><li style="box-sizing: border-box; padding: 0px 5px;">17</li><li style="box-sizing: border-box; padding: 0px 5px;">18</li><li style="box-sizing: border-box; padding: 0px 5px;">19</li><li style="box-sizing: border-box; padding: 0px 5px;">20</li><li style="box-sizing: border-box; padding: 0px 5px;">21</li><li style="box-sizing: border-box; padding: 0px 5px;">22</li><li style="box-sizing: border-box; padding: 0px 5px;">23</li><li style="box-sizing: border-box; padding: 0px 5px;">24</li><li style="box-sizing: border-box; padding: 0px 5px;">25</li><li style="box-sizing: border-box; padding: 0px 5px;">26</li><li style="box-sizing: border-box; padding: 0px 5px;">27</li></ul>

效果

我依据GIMP中插件contrast-retinex.c用Matlab实现了一遍，效果如下

宽动态图像

水下图像

低照度图像

雾天图像

评论

Retinex算法在图像增强诸多领域应用广泛，如上述宽动态图像增强，水下图像增强，雾天图像增强，低照度图像增强等等，更多可以参考NASA关于Retinex的介绍。Retinex的实现，C语言用了递归高斯滤波器，我不知道Matlab是否有对应的函数，我在这里是将图像转换到频率域高斯低通滤波处理的，因此算法时间耗费也并不大。OpenCV版本的Retinex算法已传到CSDN上共享，下载请点击这里，记得给好评哦。

更多阅读

http://dragon.larc.nasa.gov/ 
http://www.cnblogs.com/Imageshop/archive/2013/04/17/3026881.html 
A multiscale retinex for bridging the gap between color images and the human observation of scenes