Retinex算法

概述

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

（1）真实世界是无颜色的，我们所感知的颜色是光与物质的相互作用的结果。我们见到的水是无色的，但是水膜—肥皂膜却是显现五彩缤纷，那是薄膜表面光干涉的结果。


（2）每一颜色区域由给定波长的红、绿、蓝三原色构成的；

（3）三原色决定了每个单位区域的颜色。

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

不同于传统的线性、非线性的只能增强图像某一类特征的方法，Retinex可以在动态范围压缩、边缘增强和颜色恒常三个方面达到平衡，因此可以对各种不同类型的图像进行自适应的增强。

原理

一幅给定的图像S(x,y)可以分解为两个不同的图像：反射图像R(x,y)和亮度图像（也有人称之为入射图像）L(x,y)，其原理如下图所示：

L又称为光照分量，R又称为反射分量。Retinex的基本思想是，在S中，降低L的影响，从而尽量保留物体本质的反射属性，即R。

求解R是一个奇异问题，只能通过数学方法近似估算。Retinex是一个框架，采用用不同的亮度图像估计算法，其具体形式也不同，但基本流程如下：

一般直接通过空间平滑得到照射图像估计，常用的是高斯核的低通滤波。

retinex在对数域里进行：

以上是单尺度retinex（SSR），多尺度（MSR）的表达式如下：

K是尺度数，通常取K为3，并且权值w1=w2=w3=1/3，F是中心环绕函数，即高斯核函数、卷积核：

其中λ是归一化常数，确保卷积核内的积分为1；c为卷积核尺寸，retinex中常用的卷积核半径为15,80,120（小、中、大3个尺度）

中心环绕函数可以估计出光照分量L所对应图像的低频成分，并从S中去除低频成分，剩下来的就是高频分量，即反射分量R，保留了图像的边缘细节。

以上是经典Retinex算法，其缺点：都假设光照是低频的，即光照是缓慢变化的，但实际情况并非如此。在亮度相差很大的边缘处，光照变化是非平滑的。此时，在亮度差异较大的区域会产生光晕现象。此外，边缘锐化不足、阴影边界突兀、部分颜色扭曲、纹理不清晰、高光区域（过曝）的细节没有明显改善也是其缺点。

光晕现象是经典retinex算法的最大缺点，它使图像发生了失真，而在图像增强中，失真是无论如何不允许出现的。光晕可以理解为光照变化剧烈的地方，被当作细节增强了。

对光晕现象更深层次的解释为：实际上，由一幅图像分解成两幅图像在数学上是一个“病态”问题。因为单一的像素亮度表达了照度信息和图像内容双重含义，在这种情况下要想准确的估计山照度情况只能依靠计算像素本身和图像中其它像素的关系。空域滤波在计算照度时主要根据剧围像素的位置远近赋予不同的权重来估算当前像素的照度（空域=二维特征空间），没有充分考虑到像素亮度本身的意义（空域+亮度=三维特征空间）。这直接导致了在估算图像高对比度边缘区域照度的时候，边缘两边的高低值像素会相互影响：高值像素的照度受相邻低值像素的影响导致所估计的照度较低；低值像素的照度受相邻高值像素的影响导致所估计的照度较高。这会造成这一区域照度估计失真，出现结果图像中的“光晕伪影”。

人们为了解决这一问题，相继提出了改进的措施。主要思路是改进算法框图中的“照度图像估计”这一块。经典算法使用的是高斯低通滤波，对图像的每一处都使用同一个滤波模板；而改进的方法，将考虑邻近像素间的灰度值差异。常见的方法有，双边滤波、各向异性扩散。双边滤波有快速算法，能够满足实时性要求。

经过以上步骤获得了反射分量，一般还会加入增强环节，例如对数变换、指数变换、伽马校正等一些非线性的灰度变换，进一步提高细节的呈现。

matlab程序

产生中心环绕函数

sigma=1; %sigma
delta=1; %取值步长为1
g_width_1=11; %卷积核尺寸
gauss_1=zeros(g_width_1,g_width_1);
for i=-1*floor(g_width_1/2):floor(g_width_1/2)
    for j=-1*floor(g_width_1/2):floor(g_width_1/2)
        gauss_1(i+1+floor(g_width_1/2),j+1+floor(g_width_1/2))=exp(-1*((i*delta)^2+(j*delta)^2)/(2*sigma^2))/(2*sigma^2);
    end
end
gauss_1=gauss_1/sum(sum(gauss_1));