MSRCR（Multi-Scale Retinex with Color Restore）

引言

始于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)]}。

流程

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

入射分量获取

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

色彩恢复

GIMP中色彩恢复代码如下：

<pre name="code" class="cpp">/*  Final calculation with original value and cumulated filter values.  The parameters gain, alpha and offset are constants.  */   /* Ci(x,y)=log[a Ii(x,y)]-log[ Ei=1-s Ii(x,y)] */   alpha  = 128.0f;   gain   = 1.0f;   offset = 0.0f;   for ( i = 0; i < size; i += bytes )   {   float logl;   psrc = src+i;   pdst = dst+i;   logl = (float)log( (float)psrc[0] + (float)psrc[1] + (float)psrc[2] + 3.0f );   pdst[0] = gain * ((float)(log(alpha * (psrc[0]+1.0f)) - logl) * pdst[0]) + offset;   pdst[1] = gain * ((float)(log(alpha * (psrc[1]+1.0f)) - logl) * pdst[1]) + offset;   pdst[2] = gain * ((float)(log(alpha * (psrc[2]+1.0f)) - logl) * pdst[2]) + offset;   }

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

 <pre name="code" class="cpp">  /*  Adapt the dynamics of the colors according to the statistics of the first and second order.  The use of the variance makes it possible to control the degree of saturation of the colors.  */   pdst = dst;compute_mean_var( pdst, &mean, &var, size, bytes );   mini = mean - rvals.cvar*var;   maxi = mean + rvals.cvar*var;   range = maxi - mini;   if ( !range ) range = 1.0;   for ( i = 0; i < size; i+= bytes )   {   psrc = src + i;   pdst = dst + i;   for (j = 0 ; j < 3 ; j++)   {   float c = 255 * ( pdst[j] - mini ) / range;   psrc[j] = (unsigned char)clip( c, 0, 255 );   }   }   free (dst);

效果

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

宽动态图像

水下图像

低照度图像

雾天图像

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

转载请保留以下信息

作者	日期	联系方式
风吹夏天	2015年5月13日	wincoder@qq.com