（二）Retinex原理以及python实现

一. Retinex的世界观

Retinex的世界观

上图告诉我们，我是看到的图像 $S(x, y)$ ，实际上是由入射光 $L(x, y)$ 与反射图像 $R(x, y)$ 结合得到。复原的目标就是从 $S(x, y)$ 中分离出 $R(x, y)$ 。

这是一个奇异问题，正常情况下是无法解的，不过好在我们有这样一个认识，入射光 $L(x, y)$ 通常在照射面上的强度变化比较缓慢（具体来说，图像上各个像素位置的入射光强度差异不会很大），而反射图像 $R(x, y)$ 由于跟物体性质关系密切，它的变化会比较陡峭。这样，我们就知道，入射光 $L(x, y)$ 是一种 低频量，而反射图像 $R(x, y)$ 是 高频量 。这样，我们就可以通过滤波的方法来恢复出 $R(x, y)$ 。

最基础的解决方法：
$R=\log S - \log (F*S)$
其中， $F*S$ 是用于获取低频量（也就是入射光），通常 $F$ 是高斯函数，上式表示对看到的图像 $S(x, y)$ 进行高斯模糊，这样就得到了低频量。通过 $\log S$ 减去 $\log (F*S)$ 就能滤除低频量。

二. 代码的实现

代码转载自Retinex
Retinex.py

import numpy as np
import cv2

def singleScaleRetinex(img, sigma):

    retinex = np.log10(img) - np.log10(cv2.GaussianBlur(img, (0, 0), sigma))

    return retinex

def multiScaleRetinex(img, sigma_list):

    retinex = np.zeros_like(img)
    for sigma in sigma_list:
        retinex += singleScaleRetinex(img, sigma)

    retinex = retinex / len(sigma_list)

    return retinex

def colorRestoration(img, alpha, beta):

    img_sum = np.sum(img, axis=2, keepdims=True)

    color_restoration = beta * (np.log10(alpha * img) - np.log10(img_sum))

    return color_restoration

def simplestColorBalance(img, low_clip, high_clip):    

    total = img.shape[0] * img.shape[1]
    for i in range(img.shape[2]):
        unique, counts = np.unique(img[:, :, i], return_counts=True)
        current = 0
        for u, c in zip(unique, counts):            
            if float(current) / total < low_clip:
                low_val = u
            if float(current) / total < high_clip:
                high_val = u
            current += c
                
        img[:, :, i] = np.maximum(np.minimum(img[:, :, i], high_val), low_val)

    return img    

def MSRCR(img, sigma_list, G, b, alpha, beta, low_clip, high_clip):

    img = np.float64(img) + 1.0

    img_retinex = multiScaleRetinex(img, sigma_list)    
    img_color = colorRestoration(img, alpha, beta)    
    img_msrcr = G * (img_retinex * img_color - b)

    for i in range(img_msrcr.shape[2]):
        img_msrcr[:, :, i] = (img_msrcr[:, :, i] - np.min(img_msrcr[:, :, i])) / \
                             (np.max(img_msrcr[:, :, i]) - np.min(img_msrcr[:, :, i])) * \
                             255
    
    img_msrcr = np.uint8(np.minimum(np.maximum(img_msrcr, 0), 255))
    img_msrcr = simplestColorBalance(img_msrcr, low_clip, high_clip)       

    return img_msrcr

def automatedMSRCR(img, sigma_list):

    img = np.float64(img) + 1.0

    img_retinex = multiScaleRetinex(img, sigma_list)

    for i in range(img_retinex.shape[2]):
        unique, count = np.unique(np.int32(img_retinex[:, :, i] * 100), return_counts=True)
        for u, c in zip(unique, count):
            if u == 0:
                zero_count = c
                break
            
        low_val = unique[0] / 100.0
        high_val = unique[-1] / 100.0
        for u, c in zip(unique, count):
            if u < 0 and c < zero_count * 0.1:
                low_val = u / 100.0
            if u > 0 and c < zero_count * 0.1:
                high_val = u / 100.0
                break
            
        img_retinex[:, :, i] = np.maximum(np.minimum(img_retinex[:, :, i], high_val), low_val)
        
        img_retinex[:, :, i] = (img_retinex[:, :, i] - np.min(img_retinex[:, :, i])) / \
                               (np.max(img_retinex[:, :, i]) - np.min(img_retinex[:, :, i])) \
                               * 255

    img_retinex = np.uint8(img_retinex)
        
    return img_retinex

def MSRCP(img, sigma_list, low_clip, high_clip):

    img = np.float64(img) + 1.0

    intensity = np.sum(img, axis=2) / img.shape[2]

    retinex = multiScaleRetinex(intensity, sigma_list)

    intensity = np.expand_dims(intensity, 2)
    retinex = np.expand_dims(retinex, 2)

    intensity1 = simplestColorBalance(retinex, low_clip, high_clip)  # Limit the value in a range

    intensity1 = (intensity1 - np.min(intensity1)) / \
                 (np.max(intensity1) - np.min(intensity1)) * \
                 255.0 + 1.0

    img_msrcp = np.zeros_like(img)
    
    for y in range(img_msrcp.shape[0]):
        for x in range(img_msrcp.shape[1]):
            B = np.max(img[y, x])
            A = np.minimum(256.0 / B, intensity1[y, x, 0] / intensity[y, x, 0])
            img_msrcp[y, x, 0] = A * img[y, x, 0]
            img_msrcp[y, x, 1] = A * img[y, x, 1]
            img_msrcp[y, x, 2] = A * img[y, x, 2]

    img_msrcp = np.uint8(img_msrcp - 1.0)

    return img_msrcp

run.py

import sys
import os

import cv2
import json

import retinex

import matplotlib.pyplot as plt
import numpy as np
from PIL import Image

data_path = 'data'
img_list = os.listdir(data_path)
if len(img_list) == 0:
    print('Data directory is empty.')
    exit()

with open('config.json', 'r') as f: config = json.load(f)

for img_name in img_list:
    if img_name == '.gitkeep':
        continue

    img = Image.open(os.path.join(data_path, img_name))
    img = np.array(img)

    img_msrcr = retinex.MSRCR(
        img,
        config['sigma_list'],
        config['G'],
        config['b'],
        config['alpha'],
        config['beta'],
        config['low_clip'],
        config['high_clip']  )
   
    img_amsrcr = retinex.automatedMSRCR(
        img,
        config['sigma_list'] )

    img_msrcp = retinex.MSRCP(
        img,
        config['sigma_list'],
        config['low_clip'],
        config['high_clip']   )

    # visualize the result
    plt.figure()
    plt.subplots_adjust(wspace=0.3, hspace=0)  # 调整子图间距
    plt.subplot(1, 4, 1); plt.title('img'); plt.imshow(img)
    plt.subplot(1, 4, 2); plt.title('img_msrcr'); plt.imshow(img_msrcr)
    plt.subplot(1, 4, 3); plt.title('img_amsrcr'); plt.imshow(img_amsrcr)
    plt.subplot(1, 4, 4); plt.title('img_msrcp'); plt.imshow(img_msrcp)
    plt.show()