一、单尺度SSR(Single Scale Retinex)理论

Retinex理论始于Land和McCann于20世纪60年代作出的一系列贡献,其基本思想是人感知到某点的颜色和亮度并不仅仅取决于该点进入人眼的绝对光线，还和其周围的颜色和亮度有关。Retinex这个词是由视网膜(Retina)和大脑皮层(Cortex)两个词组合构成的.Land之所以设计这个词，是为了表明他不清楚视觉系统的特性究竟取决于此两个生理结构中的哪一个，抑或是与两者都有关系。

Land的Retinex模型是建立在以下的基础之上的：

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

Retinex 理论的基本内容

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

如下图所示，观察者所看到的物体的图像S是由物体表面对入射光L反射得到的，反射率R由物体本身决定，不受入射光L变化。

基于Retinex的图像增强的目的就是从原始图像S中估计出光照L，从而分解出R，消除光照不均的影响，以改善图像的视觉效果,正如人类视觉系统那样。
Retinex方法的核心就是估测照度L，从图像S中估测L分量，并去除L分量，得到原始反射分量R，即：

二、Retinex理论的理解

下面介绍两个经典的Retinex算法：
基于路径的Retinex以及基于中心/环绕Retinex。

Retinex理论，与降噪类似，该理论的关键就是合理地假设了图像的构成。如果将观察者看到的图像看成是一幅带有乘性噪声的图像，那么入射光的分量就是一种乘性的，相对均匀，且变换缓慢的噪声。Retinex算法所做的就是合理地估计图像中各个位置的噪声，并除去它。

在极端情况下，我们大可以认为整幅图像中的分量都是均匀的，那么最简单的估计照度L的方式就是在将图像变换到对数域后对整幅图像求均值。因此，我设计了以下算法来验证自己的猜想，流程如下：
（1） 将图像变换到对数域 $s=\log S$；

（2） 归一化去除加性分量 $r^{\prime }=\frac{r}{\max (r)}$；

（3） 对步骤3得到的结果求指数，反变换到实数域 $R=\exp (r\cdot \log 255)$。

• 计算流程如下：
  $\mathrm{L}\mathrm{o}\mathrm{g}(\mathrm{R}(\mathrm{x},\mathrm{y}))=(\mathrm{L}\mathrm{o}\mathrm{g}(\mathrm{I}(\mathrm{x},\mathrm{y}))-\mathrm{L}\mathrm{o}\mathrm{g}(\mathrm{I}(\mathrm{x},\mathrm{y})\ast \mathrm{G}(\mathrm{x},\mathrm{y})))（\ast 乘）(1-1)$
  $\mathrm{R}(\mathrm{x},\mathrm{y})=(\mathrm{V}\mathrm{a}\mathrm{l}\mathrm{u}\mathrm{e}-\mathrm{M}\mathrm{i}\mathrm{n})/(\mathrm{M}\mathrm{a}\mathrm{x}-\mathrm{M}\mathrm{i}\mathrm{n})\ast (255-0)(1-2)$

其中：

$I(x,y)$代表被观察或照相机接收到的图像信号； $L(x,y)$代表环境光的照射分量 ； $R(x,y)$表示携带图像细节信息的目标物体的反射分量 。
$L(x,y)$是不能够求得的，只能近似求出。我们用 $I(x,y)$和一个高斯核的卷积来近似表示 $L(x,y)$。

import numpy as np
import cv2def replaceZeroes(data):min_nonzero = min(data[np.nonzero(data)])data[data == 0] = min_nonzeroreturn datadef SSR(src_img, size):L_blur = cv2.GaussianBlur(src_img, (size, size), 0)img = replaceZeroes(src_img)L_blur = replaceZeroes(L_blur)dst_Img = cv2.log(img/255.0)dst_Lblur = cv2.log(L_blur/255.0)dst_IxL = cv2.multiply(dst_Img,dst_Lblur)log_R = cv2.subtract(dst_Img, dst_IxL)dst_R = cv2.normalize(log_R,None,0,255,cv2.NORM_MINMAX)log_uint8 = cv2.convertScaleAbs(dst_R)return log_uint8if __name__ == '__main__':img = './gggg/20190701133802.png'size = 3src_img = cv2.imread(img)b_gray, g_gray, r_gray = cv2.split(src_img)b_gray = SSR(b_gray, size)g_gray = SSR(g_gray, size)r_gray = SSR(r_gray, size)result = cv2.merge([b_gray, g_gray, r_gray])cv2.imshow('img',src_img)cv2.imshow('result',result)cv2.waitKey(0)cv2.destroyAllWindows()

三、多尺度MSR(Multi-Scale Retinex)

MSR是在SSR基础上发展来的，优点是可以同时保持图像高保真度与对图像的动态范围进行压缩的同时，MSR也可实现色彩增强、颜色恒常性、局部动态范围压缩、全局动态范围压缩，也可以用于X光图像增强。最为经典的就是3尺度的，大、中、小，既能实现图像动态范围的压缩，又能保持色感的一致性较好。

计算流程

(1) 需要对原始图像进行每个尺度的高斯模糊，得到模糊后的图像 $L_i(x,y)$,其中小标 $i$表示尺度数。
(2) 对每个尺度下进行累加计算
$\log [R(x,y)]=\log [R(x,y)]+Weight(i)\ast (\log [I_i(x,y)]-\log [L_i(x,y)]);$
其中Weight(i)表示每个尺度对应的权重，要求各尺度权重之和必须为1，经典的取值为等权重。

import numpy as np
import cv2def replaceZeroes(data):min_nonzero = min(data[np.nonzero(data)])data[data == 0] = min_nonzeroreturn datadef MSR(img, scales):weight = 1 / 3.0scales_size = len(scales)h, w = img.shape[:2]log_R = np.zeros((h, w), dtype=np.float32)for i in range(scales_size):img = replaceZeroes(img)L_blur = cv2.GaussianBlur(img, (scales[i], scales[i]), 0)L_blur = replaceZeroes(L_blur)dst_Img = cv2.log(img/255.0)dst_Lblur = cv2.log(L_blur/255.0)dst_Ixl = cv2.multiply(dst_Img, dst_Lblur)log_R += weight * cv2.subtract(dst_Img, dst_Ixl)dst_R = cv2.normalize(log_R,None, 0, 255, cv2.NORM_MINMAX)log_uint8 = cv2.convertScaleAbs(dst_R)return log_uint8if __name__ == '__main__':img = './gggg/20190701133802.png'scales = [15,101,301]  # [3,5,9]  #看不出效果有什么差别src_img = cv2.imread(img)b_gray, g_gray, r_gray = cv2.split(src_img)b_gray = MSR(b_gray, scales)g_gray = MSR(g_gray, scales)r_gray = MSR(r_gray, scales)result = cv2.merge([b_gray, g_gray, r_gray])cv2.imshow('img',src_img)cv2.imshow('MSR_result',result)cv2.waitKey(0)cv2.destroyAllWindows()

四、MSRCR & MSRCP

由于R是对数域的输出，要转换为数字图像，必须将他们量化为[0,255]的数字图像范畴，关于这个量化的算法，有这极为重要的意义，他的好坏直接决定了最终输出的图像的品质。

目前，有4种方式进行处理：

• 第一种，直接线性量化，在单尺度SSR中采用的就是这种方法。
  这种方式，由于Retinex数据处理后的高动态特性，数据分布很广，会出现严重的两极化现象，一般难以获得满意的结果。

• 第二种，就是在经典的MSRCR(Multi-Scale Retinex with Color Restoration)文章《A Multiscale Retinex for Bridging the Gap Between Color Images and the Human Observation of Scenes》中提出的Canonical Gain/set 算法。计算公式如：
  $R_{MSRC{R}_i}(x,y)=G[R_{MSRC{R}_i}(x,y)-b]$
  其中G和b为经验参数。

• 第三种，实在上述文章中提到的Simplest Color Balance（我简写为SCR）方式，这种方式的处理类似于Photoshop中的自动色阶，他把数据按照一定的百分比去除最小和最大的部分，然后中间的部分重新线性量化到0和255之间。

• 第四种，就是GIMP的Retinex算法，这个可详见 带色彩恢复的多尺度视网膜增强算法（MSRCR）的原理、实现及应用 一文的描述。

还有一种方式，就是大家知道HDR的过程吧，他也是将高动态的数据量化到图像的可视范围，因此可以直接将这类算法应用与这个问题上。我也做了实验，效果似乎一般。

在用第二种或第三种方式处理时，最好还需要有个Color Restoration的过程，因为如果直接对MSR处理的结果进行量化，得到的图像往往整体偏灰度，这是由于原始的彩色值经过log处理后的数据范围就比较小了，这样各通道之间的差异也很小，而之后的线性量化比log曲线要平滑很多，因此整体就丧失了彩色。

论文中提出了修正方式如下：　
$$I_i^{\prime }(x,y)=\frac{I_i(x,y)}{{\sum }_{j=1}^S{I}_j(x,y)}$$

$$C_i(x,y)=\beta \log [\alpha I_i^{\prime }(x,y)]$$

$$R_{MSRC{R}_i}(x,y)=C_i(x,y)\ast R_{MS{R}_i}(x,y)$$

其中 $\beta =46，\alpha =125$为经验参数，测试取1试验表明效果还不错。
对于一些原始图像HUE较为合理的图，如果用经典的MSRCR算法，会导致处理后的图容易偏色，上述论文提出了对图像的Intensity数据进行Retinex处理，然后再把数据根据原始的RGB的比例映射到每个通道，这样就能在保留原始颜色分布的基础上增强图像，文章中称其为MSRCP。

MSRCR 算法步骤

该代码有误待修改

import numpy as np
import cv2def replaceZeroes(data):min_nonzero = min(data[np.nonzero(data)])data[data == 0] = min_nonzeroreturn datadef simple_color_balance(input_img, s1, s2):h, w = input_img.shape[:2]temp_img = input_img.copy()one_dim_array = temp_img.flatten()sort_array = sorted(one_dim_array)per1 = int((h * w) * s1 / 100)minvalue = sort_array[per1]per2 = int((h * w) * s2 / 100)maxvalue = sort_array[(h * w) - 1 - per2]# 实施简单白平衡算法if (maxvalue <= minvalue):out_img = np.full(input_img.shape, maxvalue)else:scale = 255.0 / (maxvalue - minvalue)out_img = np.where(temp_img < minvalue, 0)    # 防止像素溢出out_img = np.where(out_img > maxvalue, 255)   # 防止像素溢出out_img = scale * (out_img - minvalue)        # 映射中间段的图像像素out_img = cv2.convertScaleAbs(out_img)return out_imgdef MSRCR(img, scales, s1, s2):h, w = img.shape[:2]scles_size = len(scales)log_R = np.zeros((h, w), dtype=np.float32)img_sum = np.add(img[:,:,0],img[:,:,1],img[:,:,2])img_sum = replaceZeroes(img_sum)gray_img = []for j in range(3):img[:, :, j] = replaceZeroes(img[:, :, j])for i in range(0, scles_size):L_blur = cv2.GaussianBlur(img[:, :, j], (scales[i], scales[i]), 0)L_blur = replaceZeroes(L_blur)dst_img = cv2.log(img[:, :, j]/255.0)dst_Lblur = cv2.log(L_blur/255.0)dst_ixl = cv2.multiply(dst_img, dst_Lblur)log_R += cv2.subtract(dst_img, dst_ixl)MSR = log_R / 3.0MSRCR = MSR * (cv2.log(125.0 * img[:, :, j]) - cv2.log(img_sum))gray = simple_color_balance(MSRCR, s1, s2)gray_img.append(gray)return gray_imgif __name__ == '__main__':img = './gggg/20190701133802.png'scales = [15,101,301]s1, s2 = 2,3src_img = cv2.imread(img)MSRCR_Out = MSRCR(src_img, scales, s1, s2)result = cv2.merge(MSRCR_Out[0],MSRCR_Out[1],MSRCR_Out[2])cv2.imshow('img',src_img)cv2.imshow('MSR_result',result)cv2.waitKey(0)cv2.destroyAllWindows()

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;
• 就是经过这么简单的处理，实践证明可以取得非常好的效果，下面贴出一些处理后的效果。

MSRCP算法实现：（代码有问题，改进待续）

import numpy as np
import cv2def replaceZeroes(data):min_nonzero = min(data[np.nonzero(data)])data[data == 0] = min_nonzeroreturn datadef simple_color_balance(input_img, s1, s2):h, w = input_img.shape[:2]temp_img = input_img.copy()one_dim_array = temp_img.flatten()sort_array = sorted(one_dim_array)per1 = int((h * w) * s1 / 100)minvalue = sort_array[per1]per2 = int((h * w) * s2 / 100)maxvalue = sort_array[(h * w) - 1 - per2]# 实施简单白平衡算法if (maxvalue <= minvalue):out_img = np.full(input_img.shape, maxvalue)else:scale = 255.0 / (maxvalue - minvalue)out_img = np.where(temp_img < minvalue, 0,temp_img)   # 防止像素溢出out_img = np.where(out_img > maxvalue, 255,out_img)   # 防止像素溢出out_img = scale * (out_img - minvalue)                # 映射中间段的图像像素out_img = cv2.convertScaleAbs(out_img)return out_imgdef MSRCP(img, scales, s1, s2):h, w = img.shape[:2]scales_size = len(scales)B_chan = img[:, :, 0]G_chan = img[:, :, 1]R_chan = img[:, :, 2]log_R = np.zeros((h, w), dtype=np.float32)array_255 = np.full((h, w),255.0,dtype=np.float32)I_array = (B_chan + G_chan + R_chan) / 3.0I_array = replaceZeroes(I_array)for i in range(0, scales_size):L_blur = cv2.GaussianBlur(I_array, (scales[i], scales[i]), 0)L_blur = replaceZeroes(L_blur)dst_I = cv2.log(I_array/255.0)dst_Lblur = cv2.log(L_blur/255.0)dst_ixl = cv2.multiply(dst_I, dst_Lblur)log_R += cv2.subtract(dst_I, dst_ixl)MSR = log_R / 3.0Int1 = simple_color_balance(MSR, s1, s2)B_array = np.maximum(B_chan,G_chan,R_chan)A = np.minimum(array_255 / B_array, Int1/I_array)R_channel_out = A * R_chanG_channel_out = A * G_chanB_channel_out = A * B_chanMSRCP_Out_img = cv2.merge([B_channel_out, G_channel_out, R_channel_out])MSRCP_Out = cv2.convertScaleAbs(MSRCP_Out_img)return MSRCP_Outif __name__ == '__main__':img = './gggg/20190701133802.png'scales = [15,101,301]s1, s2 = 2,3src_img = cv2.imread(img)result = MSRCP(src_img, scales, s1, s2)cv2.imshow('img',src_img)cv2.imshow('MSR_result',result)cv2.waitKey(0)cv2.destroyAllWindows()

特别鸣谢
https://blog.csdn.net/lz0499/article/details/81840201
http://www.cnblogs.com/sleepwalker/p/3676600.html?utm_source=tuicool&utm_medium=referral
https://blog.csdn.net/weixin_38285131/article/details/88097771
https://blog.csdn.net/ajianyingxiaoqinghan/article/details/71435098
https://blog.csdn.net/lz0499/article/details/81154937