水下图像增强与复原技术：暗通道先验算法详解及代码实现

水下图像增强与复原技术：暗通道先验算法详解及代码实现

引用

CSDN

1.

https://blog.csdn.net/2301_79776310/article/details/144070138

水下图像由于光的吸收和散射作用，往往存在色彩失真和模糊等问题。本文将介绍一种常用的水下图像增强算法——暗通道先验算法，并提供Matlab和Python两种语言的实现代码。

水下光图像成像模型

在水下环境中，光的能量会因为反射和散射而大量损失，这导致水下图像质量下降。Jaffe-McGlamery提出的水下光图像成像模型，将到达摄像机的光分为三部分：直接光、前向散射光和后向散射光。

直接光

直接光是指直接从物体反射进入摄像机的光线，其表达式为：

其中，表示水中物体所在处位置光强，表示摄像机直接接受的光强，表示总的衰减系数，包含由吸收和散射造成的衰减，为图像RGB三个颜色通道中的某一通道，考虑到水对不同光谱的吸收程度不一样，因此总衰减系数根据颜色通道的变化而变化，为物体与摄像机之间的距离。

前向散射光

前向散射光是物体反射光经小角度散射后形成的光，会导致图像模糊，可用卷积表示为：

其中，“*”表示卷积运算，是点扩散函数，定义如下：

后向散射光

后向散射光是由周围环境光照被水中物体如浮游生物、悬浮颗粒等散射后进入摄像机的光，可用如下公式表示：

其中，称之为背景光。则总光强的详细表达式为：

当物体与摄像机相距很小时，可以忽略前向散射光引起的模糊，因此上式可简化为：

而在水下图像中，该模型的常见形式为：

其中 I(x) 代表退化的水下图像，J(x) 代表清晰的图像，B 是水下环境中的环境光， t(x) 是水下场景的光传输函数。

暗通道先验算法

建立完整模型

为了从退化的水下图像I(x)中恢复出清晰的图像J(x)，我们需要求解环境光B和光传输函数t(x)。

求解环境光B

通过分析大量清晰的水下图像，可以发现它们的暗通道像素值区域接近零。暗通道是指在清晰无雾的图片中，除天空区域外的任一局部区域像素至少有一个通道值（红或绿或蓝）很低，几乎趋近于零。

因此，我们可以通过计算暗场图像中亮度最高的几个像素值来估计环境光B。具体步骤如下：

1. 计算输入图像I(x)的暗场图像
2. 在暗场图像中寻找最亮的几个像素值
3. 这些像素值在原图像J(x)中对应的像素值应该接近零
4. 将这些像素值代入成像模型公式，可以得到环境光B的估计值

求解光传输函数t(x)

求解t(x)的过程相对简单，主要涉及公式推导：

1. 将成像模型公式两边同时除以环境光B
2. 计算上式的暗通道，即求最小通道和作最小值滤波
3. 由于J(x)在暗通道处接近零，且B已经求得，可以推导出t(x)的表达式

代码实现

以下是暗通道先验算法的Matlab实现代码：

function dark = DarkChannel(im, sz)
    b = im(:,:,1);
    g = im(:,:,2);
    r = im(:,:,3);
    dc = min([min(r, g), b]);
    kernel = strel('rectangle', [sz sz]);
    dark = imerode(dc, kernel);
end

function A = AtmLight(im, dark)
    [h, w] = size(im);
    imsz = h * w;
    numpx = max(floor(imsz / 1000), 1);
    darkvec = reshape(dark, imsz, 1);
    imvec = reshape(im, imsz, 3);
    [~, indices] = sort(darkvec);
    indices = indices(end-numpx+1:end);
    atmsum = zeros(1, 3);
    for ind = 1:numpx
        atmsum = atmsum + imvec(indices(ind), :);
    end
    A = atmsum / numpx;
end

function transmission = TransmissionEstimate(im, A, sz)
    omega = 0.95;
    im3 = zeros(size(im));
    for ind = 1:3
        im3(:,:,ind) = im(:,:,ind) / A(ind);
    end
    dark = DarkChannel(im3, sz);
    transmission = 1 - omega * dark;
end

function q = Guidedfilter(im, p, r, eps)
    mean_I = imfilter(double(im), fspecial('average', [r r]));
    mean_p = imfilter(double(p), fspecial('average', [r r]));
    mean_Ip = imfilter(double(im) .* double(p), fspecial('average', [r r]));
    cov_Ip = mean_Ip - mean_I .* mean_p;
    mean_II = imfilter(double(im) .* double(im), fspecial('average', [r r]));
    var_I = mean_II - mean_I .* mean_I;
    a = cov_Ip ./ (var_I + eps);
    b = mean_p - a .* mean_I;
    mean_a = imfilter(a, fspecial('average', [r r]));
    mean_b = imfilter(b, fspecial('average', [r r]));
    q = mean_a .* double(im) + mean_b;
end

function t = TransmissionRefine(im, et)
    gray = rgb2gray(im);
    gray = im2double(gray);
    r = 60;
    eps = 0.0001;
    t = Guidedfilter(gray, et, r, eps);
end

function res = Recover(im, t, A, tx)
    if nargin < 4
        tx = 0.1;
    end
    t = max(t, tx);
    res = zeros(size(im));
    for ind = 1:3
        res(:,:,ind) = (im(:,:,ind) - A(ind)) ./ t + A(ind);
    end
end

以下是暗通道先验算法的Python实现代码：

import math
import numpy as np
import sys
import cv2
import os

def DarkChannel(im, sz):
    b, g, r = cv2.split(im)
    dc = cv2.min(cv2.min(r, g), b);
    kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (sz, sz))
    dark = cv2.erode(dc, kernel)
    return dark

def AtmLight(im, dark):
    [h, w] = im.shape[:2]
    imsz = h * w
    numpx = int(max(math.floor(imsz / 1000), 1))
    darkvec = dark.reshape(imsz);
    imvec = im.reshape(imsz, 3);
    indices = darkvec.argsort();
    indices = indices[imsz - numpx::]
    atmsum = np.zeros([1, 3])
    for ind in range(1, numpx):
        atmsum = atmsum + imvec[indices[ind]]
    A = atmsum / numpx;
    return A

def TransmissionEstimate(im, A, sz):
    omega = 0.95;
    im3 = np.empty(im.shape, im.dtype);
    for ind in range(0, 3):
        im3[:, :, ind] = im[:, :, ind] / A[0, ind]
    transmission = 1 - omega * DarkChannel(im3, sz);
    return transmission

def Guidedfilter(im, p, r, eps):
    mean_I = cv2.boxFilter(im, cv2.CV_64F, (r, r));
    mean_p = cv2.boxFilter(p, cv2.CV_64F, (r, r));
    mean_Ip = cv2.boxFilter(im * p, cv2.CV_64F, (r, r));
    cov_Ip = mean_Ip - mean_I * mean_p;
    mean_II = cv2.boxFilter(im * im, cv2.CV_64F, (r, r));
    var_I = mean_II - mean_I * mean_I;
    a = cov_Ip / (var_I + eps);
    b = mean_p - a * mean_I;
    mean_a = cv2.boxFilter(a, cv2.CV_64F, (r, r));
    mean_b = cv2.boxFilter(b, cv2.CV_64F, (r, r));
    q = mean_a * im + mean_b;
    return q;

def TransmissionRefine(im, et):
    gray = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY);
    gray = np.float64(gray) / 255;
    r = 60;
    eps = 0.0001;
    t = Guidedfilter(gray, et, r, eps);
    return t;

def Recover(im, t, A, tx=0.1):
    res = np.empty(im.shape, im.dtype);
    t = cv2.max(t, tx);
    for ind in range(0, 3):
        res[:, :, ind] = (im[:, :, ind] - A[0, ind]) / t + A[0, ind]
    return res

if __name__ == '__main__':
    try:
        fn = sys.argv[1]
    except:
        fn = r'D:\math\picture'  # 修改为文件夹路径
    def nothing(*argv):
        pass
    output_dir = 'D:\\math\\youhua'  # 输出文件夹路径，可以根据需要自行修改
    if not os.path.exists(output_dir):  # 如果输出文件夹不存在，创建它
        os.makedirs(output_dir)
    img_dir = os.listdir(fn)  # 获取文件夹下的所有文件名
    for img_file in img_dir:
        if img_file.endswith('.png') or img_file.endswith('.jpg'):  # 仅处理png和jpg格式的图片
            src = cv2.imread(os.path.join(fn, img_file))  # 使用os.path.join连接路径和文件名
            ...  # 后续处理代码不变
            I = src.astype('float64') / 255;
            dark = DarkChannel(I, 15);
            A = AtmLight(I, dark);
            te = TransmissionEstimate(I, A, 15);
            t = TransmissionRefine(src, te);
            J = Recover(I, t, A, 0.1)  # 恢复图片并乘以255
            output_path = os.path.join(output_dir,
                                       os.path.splitext(img_file)[0] + '.png')  # 构建输出文件路径，以原文件名（去掉扩展名）加上输出路径和新的扩展名.png
            cv2.imwrite(output_path, J * 255)  # 写入输出文件夹中

实验结果

以下是使用暗通道先验算法处理后的水下图像示例：

需要注意的是，暗通道先验算法在水下图像处理中存在一些局限性，例如：

1. 水下环境中红光最先衰减，导致暗通道特征与空气中的暗原色先验理论不完全适用
2. 去雾后的图像可能存在过饱和问题

因此，在实际应用中可能需要进一步优化算法，例如采用红色逆通道先验算法或中值暗通道先验算法等。