MSRCR(Multi-Scale Retinex with Color Restore)多尺度Retinex图像增强

时间:2023-02-04 23:42:12

引言

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

算法概述

Retinex 理论的基本内容是物体的颜色是由物体对长波(红)、中波(绿)和短波(蓝)光线的反射能力决定的,而不是由反射光强度的绝对值决定的;物体的色彩不受光照非均性的影响,具有一致性,即Retinex理论是以色感一致性(颜色恒常性)为基础的。如下图所示,观察者所看到的物体的图像S是由物体表面对入射光L反射得到的,反射率R由物体本身决定,不受入射光L变化。 
MSRCR(Multi-Scale Retinex with Color Restore)多尺度Retinex图像增强
对于观察图像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)=kexp(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))=Weight1log(Rσ1(x,y))+Weight2log(Rσ2(x,y))+Weight3log(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种不同尺度下获取的亮度图,等分权重。 
MSRCR(Multi-Scale Retinex with Color Restore)多尺度Retinex图像增强

色彩恢复

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实现了一遍,效果如下

宽动态图像

MSRCR(Multi-Scale Retinex with Color Restore)多尺度Retinex图像增强
MSRCR(Multi-Scale Retinex with Color Restore)多尺度Retinex图像增强

水下图像

MSRCR(Multi-Scale Retinex with Color Restore)多尺度Retinex图像增强
MSRCR(Multi-Scale Retinex with Color Restore)多尺度Retinex图像增强

低照度图像

MSRCR(Multi-Scale Retinex with Color Restore)多尺度Retinex图像增强

雾天图像

MSRCR(Multi-Scale Retinex with Color Restore)多尺度Retinex图像增强
MSRCR(Multi-Scale Retinex with Color Restore)多尺度Retinex图像增强

评论

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