1 背景介绍

学习参考来自：OpenCV基础（24）改善夜间图像的照明

源码：

链接：https://pan.baidu.com/s/1cMDB-fDu9CXc0zJ97RT8eA 
提取码：123a

论文：
Shi Z, Zhu M M, Guo B, et al. Nighttime low illumination image enhancement with single image using bright/dark channel prior[J]. EURASIP Journal on Image and Video Processing, 2018, 2018: 1-15.

方法描述：
基于双通道（bright chennel 和 dark channel）先验的方法对单个图像进行低照度图像增强。

流程简介：

2 代码实现

python"># -*- coding: utf-8 -*-
"""
Created on Tue Sep  8 01:00:10 2020

@author: hp
"""

import cv2 as cv
import numpy as np
from guidedfilter import guided_filter


def get_illumination_channel(I, w):
    M, N, _ = I.shape  # 900, 1440
    padded = np.pad(I, ((int(w / 2), int(w / 2)), (int(w / 2), int(w / 2)), (0, 0)), 'edge')  # (914, 1454, 3)
    darkch = np.zeros((M, N))
    brightch = np.zeros((M, N))

    for i, j in np.ndindex(darkch.shape):
        darkch[i, j] = np.min(padded[i:i + w, j:j + w, :])
        brightch[i, j] = np.max(padded[i:i + w, j:j + w, :])

    return darkch, brightch  # (900, 1440)


def get_atmosphere(I, brightch, p=0.1):
    M, N = brightch.shape  # 900, 1440
    flatI = I.reshape(M * N, 3)  # (1296000, 3)
    flatbright = brightch.ravel()  # (1296000,)

    searchidx = (-flatbright).argsort()[:int(M * N * p)]  # - min to max = max to min, get top 10% max, 129600
    A = np.mean(flatI.take(searchidx, axis=0), dtype=np.float64, axis=0)  # array([0.32973554, 0.33248184, 0.29444357])
    return A


def get_initial_transmission(A, brightch):
    A_c = np.max(A)  # 0.33248184458985186
    init_t = (brightch - A_c) / (1. - A_c)  # (900, 1440)
    return (init_t - np.min(init_t)) / (np.max(init_t) - np.min(init_t))


def get_corrected_transmission(I, A, darkch, brightch, init_t, alpha, omega, w):
    im3 = np.empty(I.shape, I.dtype)  # (900, 1440, 3)
    for ind in range(0, 3):
        im3[:, :, ind] = I[:, :, ind] / A[ind]  # array([0.32973554, 0.33248184, 0.29444357])
    dark_c, _ = get_illumination_channel(im3, w)
    dark_t = 1 - omega * dark_c
    corrected_t = init_t
    diffch = brightch - darkch

    for i in range(diffch.shape[0]):
        for j in range(diffch.shape[1]):
            if (diffch[i, j] < alpha):
                corrected_t[i, j] = dark_t[i, j] * init_t[i, j]

    return np.abs(corrected_t)


def get_final_image(I, A, refined_t, tmin):
    refined_t_broadcasted = np.broadcast_to(refined_t[:, :, None], (refined_t.shape[0], refined_t.shape[1], 3))  # (900, 1440)->(900, 1440, 3)
    J = (I - A) / (np.where(refined_t_broadcasted < tmin, tmin, refined_t_broadcasted)) + A

    return (J - np.min(J)) / (np.max(J) - np.min(J))


def dehaze(I, tmin, w, alpha, omega, p, eps, reduce=False):
    m, n, _ = I.shape  # 900, 1440
    Idark, Ibright = get_illumination_channel(I, w)  # cv.imwrite("xxx.jpg", xxx*255)
    A = get_atmosphere(I, Ibright, p)  # array([0.32973554, 0.33248184, 0.29444357])

    init_t = get_initial_transmission(A, Ibright)  # (900, 1440) cv.imwrite("init_t.jpg", init_t*255)
    if reduce:
        init_t = reduce_init_t(init_t)
    corrected_t = get_corrected_transmission(I, A, Idark, Ibright, init_t, alpha, omega, w)  # cv.imwrite("corrected_t.jpg", corrected_t*255)

    normI = (I - I.min()) / (I.max() - I.min())
    refined_t = guided_filter(normI, corrected_t, w, eps) # cv.imwrite("refined_t.jpg", refined_t*255)
    J_refined = get_final_image(I, A, refined_t, tmin)  # cv.imwrite("J_refined.jpg", J_refined*255)

    enhanced = (J_refined * 255).astype(np.uint8)  # cv.imwrite("enhanced.jpg", enhanced)
    f_enhanced = cv.detailEnhance(enhanced, sigma_s=10, sigma_r=0.15)  # cv.imwrite("f_enhanced.jpg", f_enhanced)
    f_enhanced = cv.edgePreservingFilter(f_enhanced, flags=1, sigma_s=64, sigma_r=0.2) # cv.imwrite("f_f_enhanced.jpg", f_enhanced)
    return f_enhanced


def reduce_init_t(init_t):
    init_t = (init_t * 255).astype(np.uint8)
    xp = [0, 32, 255]
    fp = [0, 32, 48]
    x = np.arange(256)
    table = np.interp(x, xp, fp).astype('uint8')
    init_t = cv.LUT(init_t, table)
    init_t = init_t.astype(np.float64) / 255
    return init_t


if __name__ == "__main__":
    im = cv.imread('./shop.png')  # (900, 1440, 3)
    # im = cv.imread('./vast_of_night.jpg')  # (900, 1440, 3)
    orig = im.copy()

    tmin = 0.1  # minimum value for t to make J image
    w = 15  # window size, which determine the corseness of prior images
    alpha = 0.4  # threshold for transmission correction
    omega = 0.75  # this is for dark channel prior
    p = 0.1  # percentage to consider for atmosphere
    eps = 1e-3  # for J image

    I = np.asarray(im, dtype=np.float64)  # Convert the input to an array.
    I = I[:, :, :3] / 255

    f_enhanced = dehaze(I, tmin, w, alpha, omega, p, eps)
    f_enhanced2 = dehaze(I, tmin, w, alpha, omega, p, eps, reduce=True)
    cv.imshow('original', orig)
    cv.imshow('F_enhanced', f_enhanced)
    cv.imshow('F_enhanced2', f_enhanced2)

    cv.imwrite("original.jpg", orig)
    cv.imwrite('F_enhanced.jpg', f_enhanced)
    cv.imwrite('F_enhanced2.jpg', f_enhanced2)

    cv.waitKey(0)
    cv.destroyAllWindows()

其中 guidedfilter 实现如下

python"># -*- coding: utf-8 -*-
"""
Created on Tue Sep  8 00:59:12 2020

@author: hp
"""

from itertools import combinations_with_replacement
from collections import defaultdict

import numpy as np
from numpy.linalg import inv

B, G, R = 0, 1, 2  # index for convenience

def boxfilter(I, r):
    """Fast box filter implementation.
    Parameters
    ----------
    I:  a single channel/gray image data normalized to [0.0, 1.0]
    r:  window radius
    Return
    -----------
    The filtered image data.
    """
    M, N = I.shape
    dest = np.zeros((M, N))
    #print(I)
    
    # cumulative sum over Y axis (tate-houkou no wa)
    sumY = np.cumsum(I, axis=0)
    #print('sumY:{}'.format(sumY))
    # difference over Y axis
    dest[:r + 1] = sumY[r:2*r + 1] # top r+1 lines
    dest[r + 1:M - r] = sumY[2*r + 1:] - sumY[:M - 2*r - 1]
    #print(sumY[2*r + 1:]) # from 2*r+1 to end lines
    #print(sumY[:M - 2*r - 1]) # same lines of above, from start
    #tile replicate sumY[-1] and line them up to match the shape of (r, 1)
    dest[-r:] = np.tile(sumY[-1], (r, 1)) - sumY[M - 2*r - 1:M - r - 1] # bottom r lines

    # cumulative sum over X axis
    sumX = np.cumsum(dest, axis=1)
    #print('sumX:{}'.format(sumX))
    # difference over X axis
    dest[:, :r + 1] = sumX[:, r:2*r + 1] # left r+1 columns
    dest[:, r + 1:N - r] = sumX[:, 2*r + 1:] - sumX[:, :N - 2*r - 1]
    dest[:, -r:] = np.tile(sumX[:, -1][:, None], (1, r)) - sumX[:, N - 2*r - 1:N - r - 1] # right r columns

    #print(dest)

    return dest


def guided_filter(I, p, r=15, eps=1e-3):
    """Refine a filter under the guidance of another (RGB) image.
    Parameters
    -----------
    I:   an M * N * 3 RGB image for guidance.
    p:   the M * N filter to be guided. transmission is used for this case.
    r:   the radius of the guidance
    eps: epsilon for the guided filter
    Return
    -----------
    The guided filter.
    """
    M, N = p.shape
    base = boxfilter(np.ones((M, N)), r) # this is needed for regularization
    
    # each channel of I filtered with the mean filter. this is myu.
    means = [boxfilter(I[:, :, i], r) / base for i in range(3)]
    
    # p filtered with the mean filter
    mean_p = boxfilter(p, r) / base

    # filter I with p then filter it with the mean filter
    means_IP = [boxfilter(I[:, :, i]*p, r) / base for i in range(3)]

    # covariance of (I, p) in each local patch
    covIP = [means_IP[i] - means[i]*mean_p for i in range(3)]

    # variance of I in each local patch: the matrix Sigma in ECCV10 eq.14
    var = defaultdict(dict)
    for i, j in combinations_with_replacement(range(3), 2):
        var[i][j] = boxfilter(I[:, :, i]*I[:, :, j], r) / base - means[i]*means[j]

    a = np.zeros((M, N, 3))
    for y, x in np.ndindex(M, N):
        #         rr, rg, rb
        # Sigma = rg, gg, gb
        #         rb, gb, bb
        Sigma = np.array([[var[B][B][y, x], var[B][G][y, x], var[B][R][y, x]],
                          [var[B][G][y, x], var[G][G][y, x], var[G][R][y, x]],
                          [var[B][R][y, x], var[G][R][y, x], var[R][R][y, x]]])
        cov = np.array([c[y, x] for c in covIP])
        a[y, x] = np.dot(cov, inv(Sigma + eps*np.eye(3)))  # eq 14

    # ECCV10 eq.15
    b = mean_p - a[:, :, R]*means[R] - a[:, :, G]*means[G] - a[:, :, B]*means[B]

    # ECCV10 eq.16
    q = (boxfilter(a[:, :, R], r)*I[:, :, R] + boxfilter(a[:, :, G], r)*I[:, :, G] + boxfilter(a[:, :, B], r)*I[:, :, B] + boxfilter(b, r)) / base

    return q

3 原理分析

输入图片

一、获取明暗通道先验

get_illumination_channel

先padding，然后滑动窗口 $\Omega$ 扫描每个像素点，size 为 w=15，窗口内像素值最小的为 dark，最大的为 bright

每个空间位置，3通道对比，实现速度还是超级慢的

Ibright.jpg

Idark.jpg

二、计算全局大气光照

明亮通道通过取前 10% 强度的平均值来计算的

get_atmosphere

r g b 每个通道中，原图（3通道）取 bright channel（单通道） 值最大的前 10% 对应空间位置的值，求平均值

python">array([0.32973554, 0.33248184, 0.29444357])

三、查找初始透射图

get_initial_transmission

处理后并归一化，和 bright 通道的差异并不是很大

init_t.jpg

和 Ibright.jpg 的差异

四、使用暗通道估计校正后的透射图

get_corrected_transmission

进行此计算是为了纠正从明亮通道先验获得的潜在错误透射估计

原图先除以全局大气光照得到新图 im3，再获取 im3 的 dark 通道 dark_c，dark_t = 1-0.75*dark_c，获取原图双通道的差距，如果差异较小，进行修复，修复的方式为 dark_t * init_t，最后结果取绝对值

corrected_t.jpg

corrected-t和init_t的比较

五、使用Guided Filter平滑透射图

refined_t.jpg

六、计算结果图像

get_final_image

J_enhanced.jpg

enhance.jpg

再进行 cv.detailEnhance 和 cv.edgePreservingFilter 加强

f_enhanced.jpg

f_f_enhanced.jpg

七、局限性

python">reduce = True

减少这些导致这些区域过度曝光的强烈白色斑点。这可以通过将值从 255 限制为某个最小值来完成

实现原理可以参考 【python】OpenCV—Brightness and Contrast adjustments

F_enhanced2.jpg

4 效果展示

输入

F_enhanced.jpg

F_enhanced2.jpg

---

输入

F_enhanced.jpg

F_enhanced2.jpg

---

输入

F_enhanced.jpg

F_enhanced2.jpg

---

输入

F_enhanced.jpg

F_enhanced2.jpg

---

输入

F_enhanced.jpg

F_enhanced2.jpg

---

输入

F_enhanced.jpg

F_enhanced2.jpg

---

输入

F_enhanced.jpg

F_enhanced2.jpg

5 附录

np.ndindex

np.ndindex 是 NumPy 库中的一个函数，它的主要功能是生成一个迭代器，该迭代器遍历给定维度的所有索引。这些索引以元组的形式返回，元组中的每个元素对应一个维度的索引。

一、功能描述

np.ndindex 函数接收一个或多个整数作为输入，这些整数代表数组的各个维度的大小。然后，它生成一个迭代器，该迭代器遍历这些维度上的所有可能索引组合。每个迭代返回的索引元组可以用于访问多维数组中的元素。

二、使用方法

示例 1：基础用法
假设我们有一个二维数组，其形状为 (3, 2)，即 3 行 2 列。我们可以使用 np.ndindex(3, 2) 来生成这个数组所有元素的索引：

python">import numpy as np  
  
for index in np.ndindex(3, 2):  
    print(index)

输出将是：

python">(0, 0)  
(0, 1)  
(1, 0)  
(1, 1)  
(2, 0)  
(2, 1)

numpy.ravel

numpy.ravel() 是 NumPy 库中的一个函数，它的主要作用是将多维数组（ndarray）展平（flatten）成一维数组，但会保持原数组数据的“C风格”内存布局。这意味着 ravel() 函数返回的数组是原数组在内存中的连续元素的视图（view），如果可能的话，不会复制数据。然而，如果需要对原数组进行任何修改而不影响 ravel() 返回的数组，或者需要确保返回的数组是原数组的副本，那么可能需要使用 .copy() 方法或者 numpy.flatten() 函数（尽管 flatten() 总是返回数组的副本）。

语法

numpy.ravel(a, order='C')

a：输入的多维数组。

order：{‘C’, ‘F’, ‘A’, ‘K’}，可选。指定了多维数组在内存中的存储顺序，以及如何将数组展平。

• ‘C’：意味着按行优先顺序进行展平（C风格的连续），即先遍历最外层的索引，然后是次外层，依此类推。
• ‘F’：意味着按列优先顺序进行展平（Fortran风格的连续），即先遍历最内层的索引，然后是次内层，依此类推。
• ‘A’：如果数组在内存中是连续的，则按‘F’顺序展平，否则按‘C’顺序展平。
• ‘K’：保留元素在内存中的顺序（这主要是为了与旧的NumPy代码兼容）。

numpy.argsort

numpy.argsort() 是 NumPy 库中的一个函数，它的作用是返回数组元素值从小到大的索引数组。换句话说，argsort() 并不是直接对数组进行排序，而是返回一个索引数组，这个索引数组指明了原数组中元素从小到大排列后的位置。

语法

numpy.argsort(a, axis=-1, kind='quicksort', order=None)

a：要排序的数组。

axis：沿着它进行排序的轴。如果没有指定或轴为-1，则数组将被展平，然后排序。

kind：{‘quicksort’, ‘mergesort’, ‘heapsort’, ‘stable’}，排序算法的选择。注意，‘stable’ 排序算法在 NumPy 中可能不是所有平台都可用，并且它的性能可能不如其他排序算法。

order：如果数组包含字段，则是要排序的字段名。否则，该参数将被忽略。

返回一个索引数组，表示原数组 a 中元素沿着指定轴排序后的索引。

cv2.detailEnhance

cv2.detailEnhance 是 OpenCV 库中的一个函数，用于增强图像的细节，使图像更加清晰和鲜明。这个函数基于多尺度Retinex算法（Multi-Scale Retinex）来实现图像的细节增强。

一、原理

多尺度Retinex算法采用了多个尺度的高斯滤波器对图像进行滤波，从而在不同尺度下获取图像的亮度信息。然后，通过将不同尺度下的图像进行加权平均来获取图像的整体亮度信息。通过将原始图像与多尺度Retinex算法得到的亮度信息相结合，可以得到增强了细节的图像，使图像的细节更加清晰和鲜明。

二、使用方法

在 OpenCV 中使用 cv2.detailEnhance 函数时，通常需要传入待处理的图像作为输入，并可以选择性地设置一些参数，如 sigma_s 和 sigma_r，这些参数分别用于控制图像的平滑程度和在平滑图像时保留边缘的重要性。

三、示例代码

下面是一个使用 cv2.detailEnhance 函数的示例代码片段（注意：cv2.detailEnhance 可能不是 OpenCV 的标准函数，这里假设其用法类似于 OpenCV 中的其他图像处理函数）：

python">import cv2  
  
# 读取图像  
img = cv2.imread('path_to_image.jpg')  
  
# 使用 detailEnhance 函数增强图像细节  
# 注意：这里的 sigma_s 和 sigma_r 是假设的参数，具体值需要根据实际情况调整  
enhanced_img = cv2.detailEnhance(img, sigma_s=10, sigma_r=0.15)  
  
# 显示原始图像和增强后的图像  
cv2.imshow('Original Image', img)  
cv2.imshow('Enhanced Image', enhanced_img)  
cv2.waitKey(0)  
cv2.destroyAllWindows()

四、注意事项

cv2.detailEnhance 可能不是 OpenCV 标准库中的直接可用函数，但在某些 OpenCV 的扩展库或特定版本中可能存在。如果标准 OpenCV 库中没有此函数，可能需要查找相关的扩展或第三方库。

在使用该函数时，需要根据具体的应用场景调整参数以获得最佳效果。
由于多尺度Retinex算法对计算资源的需求较高，处理大尺寸图像时可能会消耗较多的时间和内存。

五、展望

随着计算机视觉和图像处理技术的不断发展，cv2.detailEnhance 或类似的图像细节增强算法将会得到进一步的改进和优化。未来，可以针对其局限性进行改进，如加强对过度曝光或过度阴影图像的处理能力，优化算法性能，提高处理效率等。同时，结合深度学习技术，可以进一步提高图像细节增强的效果，使其在更多领域得到应用。

cv2.edgePreservingFilter

cv2.edgePreservingFilter 是 OpenCV 库中的一个函数，用于实现边缘保留滤波（Edge Preserving Filter，EPF）。这种滤波技术在平滑图像噪声的同时，能够尽可能地保留图像的边缘信息，从而避免图像边缘的模糊。

一、函数原型
在 Python 中使用 OpenCV 时，cv2.edgePreservingFilter 函数的原型大致如下（注意：

python

python">cv2.edgePreservingFilter(src, dst=None, flags=1, sigma_s=60, sigma_r=0.4)

二、参数解释

• src：输入图像，应为8位3通道图像。
• dst：输出图像，与输入图像具有相同的尺寸和类型。这是一个可选参数，如果不提供，函数将返回一个新的图像。
• flags：滤波器类型，是一个整数。常用的值有 cv2.RECURS_FILTER（递归滤波）和 cv2.NORMCONV_FILTER（归一化卷积滤波）。递归滤波通常比归一化卷积滤波更快，但归一化卷积滤波可能会产生更锐化的边缘。
• sigma_s：控制空间邻域大小的参数，其值越大，平滑效果越明显，边缘保持效果可能减弱。通常的取值范围在0到200之间。
• sigma_r：控制颜色相似性的参数，其值越大，颜色变化较大的区域越容易被平滑，边缘保持效果也可能减弱。通常的取值范围在0到1之间。

三、工作原理

• cv2.edgePreservingFilter 通过结合空间域和颜色域的信息来实现边缘保留滤波。在平滑图像时，它会考虑像素之间的空间距离和颜色差异。对于空间距离较近且颜色相似的像素，滤波器会倾向于将它们平滑处理；而对于颜色差异较大的像素（即边缘），滤波器则会尽量保留其原有的颜色信息，从而避免边缘的模糊。

四、应用场景

边缘保留滤波在图像处理中有广泛的应用，包括但不限于：

• 图像去噪：在去除图像噪声的同时，保持图像的边缘信息，使图像看起来更加清晰。
• 图像增强：通过平滑图像的非边缘区域，增强图像的对比度，使图像更加鲜明。
• 计算机视觉：在目标检测、图像分割等任务中，边缘信息的保留对于提高算法的准确性至关重要。

五、注意事项

• 在使用 cv2.edgePreservingFilter 函数时，需要根据具体的应用场景和图像特性调整 sigma_s 和 sigma_r 参数，以达到最佳的滤波效果。
• 由于边缘保留滤波涉及到复杂的计算，因此处理大尺寸图像时可能会消耗较多的计算资源和时间。
• OpenCV 的不同版本之间可能存在函数参数和性能的差异，因此在使用时请参考您所使用的 OpenCV 版本的官方文档。