OpenCV—Python Retinex图像增强算法_opencv 增加反光

Retinex图像增强算法

一、单尺度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的方式就是在将图像变换到对数域后对整幅图像求均值。因此，我设计了以下算法来验证自己的猜想，流程如下：

import numpy as np
import cv2
 
 
def replaceZeroes(data):
    min_nonzero = min(data[np.nonzero(data)])
    data[data == 0] = min_nonzero
    return data
 
def 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_uint8
 
 
if __name__ == '__main__':
    img = './gggg/20190701133802.png'
    size = 3
    src_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尺度的，大、中、小，既能实现图像动态范围的压缩，又能保持色感的一致性较好。

import numpy as np
import cv2
 
 
def replaceZeroes(data):
    min_nonzero = min(data[np.nonzero(data)])
    data[data == 0] = min_nonzero
    return data
 
 
def MSR(img, scales):
    weight = 1 / 3.0
    scales_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_uint8
 
 
if __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]的数字图像范畴，关于这个量化的算法，有这极为重要的意义，他的好坏直接决定了最终输出的图像的品质。

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

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

论文中提出了修正方式如下

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

MSRCR 算法步骤

该代码有误待修改

import numpy as np
import cv2
 
 
def replaceZeroes(data):
    min_nonzero = min(data[np.nonzero(data)])
    data[data == 0] = min_nonzero
    return data
 
 
def 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_img
 
 
def 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.0
        MSRCR = MSR * (cv2.log(125.0 * img[:, :, j]) - cv2.log(img_sum))
        gray = simple_color_balance(MSRCR, s1, s2)
        gray_img.append(gray)
    return gray_img
 
 
if __name__ == '__main__':
    img = './gggg/20190701133802.png'
    scales = [15,101,301]
    s1, s2 = 2,3
    src_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 cv2
 
 
def replaceZeroes(data):
    min_nonzero = min(data[np.nonzero(data)])
    data[data == 0] = min_nonzero
    return data
 
def 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_img
 
 
def 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.0
    I_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.0
    Int1 = 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_chan
    G_channel_out = A * G_chan
    B_channel_out = A * B_chan
 
    MSRCP_Out_img = cv2.merge([B_channel_out, G_channel_out, R_channel_out])
    MSRCP_Out = cv2.convertScaleAbs(MSRCP_Out_img)
 
    return MSRCP_Out
 
if __name__ == '__main__':
    img = './gggg/20190701133802.png'
    scales = [15,101,301]
    s1, s2 = 2,3
    src_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
————————————————
版权声明：本文为CSDN博主「SongpingWang」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/wsp_1138886114/article/details/83096109