基于RGB颜色空间椭球模型的去雾算法

On the effectiveness of the dark channel prior for single image dehazing by approximating with minimum volume ellipsoid

Kristofer B.Gibson, Truong Q.Nguyen

University of California San Diego

对于雾天成像模型：，利用主成分分析法分析其像素灰度值的椭球聚集状态。

对于雾天图像的某一窗口，且中总共有个像素，存在一系列位于窗口中的无雾图像的向量（均为列向量）。

设为其中的一个3x1的颜色向量，则其3x3的自相关矩阵为：

则有K-L变换为：

对角矩阵中的元素值为的特征值的递减排列，且有，其中，与为自相关矩阵的特征值，且3x3矩阵为标准正交矩阵。

设在矩形窗口内的传输函数为一定值，则有雾天像素值的自相关矩阵为：

所以雾天图像的窗口的像素灰度值的自相关矩阵与无雾图像的窗口的像素灰度值的自相关矩阵之间是比例系数关系。同时两者对应的对角矩阵也满足如上比例系数关系，推导过程为：

从而有，。

对于雾天成像模型：，其传输函数，当时，；当时，。所以传输函数的取值范围为：。

对于大气散射模型，两边取期望算符以求取窗口内部的向量的数学期望，则有：

则对于聚集的向量的期望，当时，；当时，。所以在与之间变化。由此可得结论：当景深时，，且聚集的向量的平均向量越接近于大气光向量。

又有：

由此得出结论：当景深时，，则窗口中任一向量的自相关矩阵对应的特征值矩阵中的特征值随着减小而减小，从而向量聚集的椭球的尺寸也随之减小，即相较于图像中近处区域的像素灰度值聚集的椭球，图像中远处区域的像素灰度值的聚集椭球在空间中的位置更加接近于大气光向量，且其体积更小。

上面是从理论的角度分析图像中不同区域的像素灰度值的聚集状态，下面通过实验验证理论的正确性。图（a）中，选取两个大小相同的矩形区域，红色矩形框内为远处景物的像素，蓝色矩形框内为近处景物的像素。将两个矩形框中的像素的灰度值投射到坐标空间中，观察近处景物与远处景物对应的像素灰度值的分布情况，可以得出结论：远处景物的像素的聚集向量在位置上更加接近大气光向量，且聚集的椭球的尺寸比近处景物的聚集椭球尺寸要小。这与上述用主成分分析法分析像素灰度值的聚集状态的规律是一致的。

将图像的某一矩形窗口中像素的灰度值投影到RGB颜色空间中，观察其分布状态的matlab代码如下：

clear all;
I = imread('C:\Users\Desktop\qiu\sky.jpg');
[m, n, o] = size(I);
I = double(I);
I = I ./ 255;
[region, ~] = imcrop(I);
[h, w, ch] = size(region);
col_img = zeros(h * w, 1, ch);
for i = 1 : 1 : w
    col_img(1 + (i - 1) * h : i * h, 1, :) = region(1 : h, i, :);
end

x = col_img(:, :, 1); y = col_img(:, :, 2); z = col_img(:, :, 3); 
scatter3(x, y, z, 'p');
xlabel('R'); ylabel('G'); zlabel('B');
title('RGB Color Space');
axis([0 1 0 1 0 1]);

将探测系统获取图像某一矩形窗口的像素灰度值投射到空间中，其向量分布情况如上述模型所述，设坐标空间中的任一向量为，图像的一个矩形窗口区域为，其对应的在空间中的向量聚集区域为，可由下式表示：

其中表示椭球的中心向量，其表示椭球中心在空间中的位置；为椭球的协方差矩阵，其规定了椭球的形状和方向，可用下式表示协方差矩阵：

设椭球中的向量的集合为（为椭球中向量点坐标的个数），则上述构造的椭球模型的协方差矩阵和中心向量应满足：包含集合中的所有向量，且使得椭球的体积最小。因此，可将椭球模型的求解过程转化为最小椭球体积（Minimum Volume Enclosing Ellipse, MVEE）的求解问题，即求解如下最小化问题：

通过Khachiyan算法对该最小化问题进行求解，可解得该椭球的协方差矩阵的逆矩阵与中心向量。

由于的性质，可通过奇异值分解获得其逆矩阵的特征值和特征向量：

其中为特征值递减排列构成的正交矩阵，特征值对应的特征向量构成的矩阵为。

得到上述的各个分量，就可以估计每一个像素处的传输函数值。作者提出椭球先验（Ellipsoid Prior）的传输函数估计方法：

其中为修正系数，通常取0.95，向量为(为作者笔误)。对求解得到的传输函数进一步利用导向滤波等方法对其进行优化，即可得到最终估计的传输函数。

下面给出上述求解过程的C++代码，注意要配置opencv库和eigen库，另外需在release下运行：

首先是：mvee.h

#ifndef MVEE_H
#define MVEE_H

#include <cmath>
#include <vector>
#include <Eigen/SVD>
#include <Eigen/Dense>

namespace mvee {

class Mvee {
private:
    std::vector<double> _center;        // 几何中心
    std::vector<double> _semi_axes;     // 半轴长度
    std::vector<std::vector<double>> _eigen_vecs; // 特征向量
    long int _iters;

    void decompose(Eigen::MatrixXd& mat, Eigen::VectorXd& eigen_vals, Eigen::MatrixXd& U, Eigen::MatrixXd& V);

public:
    Mvee();

    void compute(Eigen::MatrixXd& data, double eps = 0.1, double lmcoeff = 1e-8);

    std::vector<double> center() const;
    std::vector<double> semi_axes() const;
    std::vector<std::vector<double>> eigen_vecs() const;

    long int iters() const;

    ~Mvee();
};

} // namespace mvee

#endif // MVEE_H

然后是：mvee.cpp

#include "mvee.h"

namespace mvee {

Mvee::Mvee() : _iters(0) {}

void Mvee::decompose(Eigen::MatrixXd& mat, Eigen::VectorXd& eigen_vals, Eigen::MatrixXd& U, Eigen::MatrixXd& V) 
{
    Eigen::JacobiSVD<Eigen::MatrixXd> SVD(mat, Eigen::ComputeThinU | Eigen::ComputeThinV);
    eigen_vals = SVD.singularValues();
    V = SVD.matrixV();
    U = SVD.matrixU();
}

void Mvee::compute(Eigen::MatrixXd& data, double eps, double lmcoeff) 
{
    _iters = 0;
    double err = 1 + eps;
    double alpha;
    int n = data.rows();
    int p = data.cols(); // 维度p (RGB为3)
    
    Eigen::MatrixXd X = data.transpose();
    Eigen::MatrixXd Q(X.rows() + 1, X.cols());
    
    for(int i = 0; i < X.rows(); ++i) {
        Q.row(i) = X.row(i);
    }
    Q.row(X.rows()).setOnes();

    Eigen::VectorXd u = (1.0 / n) * Eigen::VectorXd::Ones(n);
    Eigen::VectorXd uhat = u;
    Eigen::MatrixXd G(p + 1, p + 1);
    Eigen::MatrixXd noiseye = Eigen::MatrixXd::Identity(G.rows(), G.cols()) * lmcoeff;
    Eigen::VectorXd g(n);
    double m;
    int idx;

    while (err > eps) {
        G = Q * u.asDiagonal() * Q.transpose() + noiseye;
        g = (Q.transpose() * G.inverse()).cwiseProduct(Q.transpose()).rowwise().sum();
        m = g.maxCoeff(&idx);

        if (m <= p + 1.0 + 1e-12) break;

        alpha = (m - p - 1) / ((p + 1) * (m - 1));
        uhat = (1 - alpha) * u;
        uhat(idx) += alpha;
        err = (uhat - u).norm();
        u = uhat;
        _iters++;
    }

    // cov_mat: 协方差矩阵 (Covariance Matrix)
    Eigen::MatrixXd cov_mat = (1.0 / p) * (
        X * u.asDiagonal() * X.transpose() -
        (X * u) * (X * u).transpose()
    );

    // center_vec: 几何中心向量
    Eigen::VectorXd center_vec = X * u;
    
    Eigen::VectorXd eigen_vals; 
    Eigen::MatrixXd V_mat;  
    Eigen::MatrixXd U_mat;
    
    decompose(cov_mat, eigen_vals, U_mat, V_mat);

    _center.resize(center_vec.size());
    for (int i = 0; i < center_vec.size(); i++) {
        _center[i] = center_vec(i);
    }

    _semi_axes.resize(eigen_vals.size());
    for (int i = 0; i < eigen_vals.size(); i++) {
        _semi_axes[i] = std::sqrt(eigen_vals(i));
    }

    _eigen_vecs.resize(V_mat.rows(), std::vector<double>(V_mat.cols()));
    for (int i = 0; i < V_mat.rows(); i++) {
        for (int j = 0; j < V_mat.cols(); j++) {
            _eigen_vecs[i][j] = V_mat(i, j);
        }
    }
}

std::vector<double> Mvee::center() const { return _center; }
std::vector<double> Mvee::semi_axes() const { return _semi_axes; }
std::vector<std::vector<double>> Mvee::eigen_vecs() const { return _eigen_vecs; }
long int Mvee::iters() const { return _iters; }
Mvee::~Mvee() {}

} // namespace mvee

最后在 main.cpp 中对其进行调用，并进行图像去雾复原：

#include "mvee.h"

#include <atomic>
#include <omp.h>

#include <iostream>
#include <algorithm>

#include <chrono>
#include <iomanip>

#include <Eigen/Dense>
#include <opencv2/opencv.hpp>
#include <opencv2/core/eigen.hpp>

#define EPS 0.1
#define LMC 1e-8

using namespace cv;
using namespace std;
using namespace Eigen;
using namespace mvee;

// 函数声明
double* estimate_airlight(const cv::Mat& image, int kernelSize);
cv::Mat guided_filter(cv::Mat I, cv::Mat p, int radius, double eps);
cv::Mat dehaze(cv::Mat& src, cv::Mat& t, double A_r, double A_g, double A_b);

int main(int argc, char** argv) 
{
    if (argc < 2) {
        std::cerr << "Usage: " << argv[0] << " <image_path>" << std::endl;
        std::cerr << "Example: ./dehaze input.jpg" << std::endl;
        return -1;
    }

    // 读取图像
    std::string pic_name = argv[1];

    std::chrono::high_resolution_clock::time_point start, end;
    // 记录开始时间
    start = std::chrono::high_resolution_clock::now();

    cv::Mat image = cv::imread(pic_name, cv::IMREAD_UNCHANGED);
    if (image.empty()) {
        std::cerr << "Error: Could not read image from path: \"" << pic_name << "\"." << std::endl;
        return -1;
    }

    int rows = image.rows;
    int cols = image.cols;

    // 获取雾天图像的灰度图 并归一化
    cv::Mat gray_img;
    cv::cvtColor(image, gray_img, cv::COLOR_BGR2GRAY);
    gray_img.convertTo(gray_img, CV_64FC1, 1.0 / 255.0);

    // 转换为浮点型并归一化
    cv::Mat convertImage;
    image.convertTo(convertImage, CV_64FC3, 1.0 / 255.0);

    // 估计大气光
    double* airlight = estimate_airlight(image, 15);

    int patchSize = 15; int radius = patchSize / 2;
    cv::Mat parseImage;
    cv::copyMakeBorder(convertImage, parseImage, radius, radius, radius, radius, cv::BORDER_REPLICATE);

    cv::Mat trans(rows, cols, CV_64FC1);
    cv::Mat theta_E_mat(rows, cols, CV_64FC1);

    // 用于进度条的原子计数器，保证多线程下计数安全
    std::atomic<int> progress_counter(0);
    std::cout << "开始估算传输函数，正在多线程并行计算..." << std::endl;

    int patch_area = patchSize * patchSize;
    // 启用 OpenMP 并行化按行处理 采用动态调度
    #pragma omp parallel for schedule(dynamic)
    for (int r = 0; r < rows; r++) {
        double* tPtr = trans.ptr<double>(r);
        // 将矩阵的实例化放在外层循环
        // 这样每个线程只会 new 一次内存 然后在列循环中被反复复用
        Eigen::MatrixXd region_roi(patch_area, 3);
        mvee::Mvee Khachiyan_solver;
        
        for (int c = 0; c < cols; c++) {
            // 避免 cv::Mat 的拷贝和 flatten_transform
            // 直接通过底层指针运算把 RGB 塞进 Eigen 矩阵 0 次内存分配
            int idx = 0;
            for (int kr = 0; kr < patchSize; kr++) {
                // 获取当前行向下的第 kr 行的指针起始位置
                const cv::Vec3d* ptr = parseImage.ptr<cv::Vec3d>(r + kr);
                for (int kc = 0; kc < patchSize; kc++) {
                    const cv::Vec3d& pixel = ptr[c + kc];
                    region_roi(idx, 0) = pixel[0];
                    region_roi(idx, 1) = pixel[1];
                    region_roi(idx, 2) = pixel[2];
                    idx++;
                }
            }

            //mvee::Mvee Khachiyan_solver;
            Khachiyan_solver.compute(region_roi, EPS, LMC);
            std::vector<double> b = Khachiyan_solver.center(); // 中心点 b
            std::vector<std::vector<double>> V = Khachiyan_solver.eigen_vecs(); // 特征矩阵 V
            std::vector<double> semi_axes = Khachiyan_solver.semi_axes();

            bool valid = true;
            for (double val : b) { 
                if (std::isnan(val)) { 
                    valid = false; 
                    break; 
                } 
            }
            if (valid) {
                for (const auto& row : V) {
                    for (double val : row) { 
                        if (std::isnan(val)) { 
                            valid = false; 
                            break; 
                        }
                    }
                    if (!valid) break;
                }
            }

            double theta_E; // 论文公式(11): Ellipsoid Prior
            if (!valid) theta_E = 0.1;
            else {
                // principal_radius: 最大半轴长
                // principal_vector: 第一主方向特征向量
                // 舍弃了原论文公式 (11) 中错误写法
                double principal_radius = semi_axes[0];
                std::vector<double> principal_vector = { V[0][0], V[1][0], V[2][0] };

                // g_pos: 主轴正向端点
                double g_pos_0 = b[0] + principal_radius * principal_vector[0];
                double g_pos_1 = b[1] + principal_radius * principal_vector[1];
                double g_pos_2 = b[2] + principal_radius * principal_vector[2];
                double norm_g_pos = std::sqrt(g_pos_0 * g_pos_0 + g_pos_1 * g_pos_1 + g_pos_2 * g_pos_2);

                // g_neg: 主轴反向端点
                double g_neg_0 = b[0] - principal_radius * principal_vector[0];
                double g_neg_1 = b[1] - principal_radius * principal_vector[1];
                double g_neg_2 = b[2] - principal_radius * principal_vector[2];
                double norm_g_neg = std::sqrt(g_neg_0 * g_neg_0 + g_neg_1 * g_neg_1 + g_neg_2 * g_neg_2);
            
                // norm_g: 论文公式(11)中的 ||g||
                double norm_g = std::min(norm_g_pos, norm_g_neg);
            
                // norm_f: 论文公式(11)中的 ||f||
                double norm_f = std::sqrt(3.0); 

                // 论文公式(11): theta_E = ||g|| / ||f||
                theta_E = norm_g / norm_f;
                theta_E_mat.at<double>(r, c) = theta_E;
            }
            
            // t_E = 1 - w * theta_E
            double w = 0.95; 
            *tPtr = 1.0 - w * theta_E;
            tPtr++;
        }

        // 进度条渲染逻辑 (每个线程完成一行后更新)
        int current_row = ++progress_counter;
        if (current_row % 5 == 0 || current_row == rows) {
            #pragma omp critical 
            {
                float progress = (float)current_row / rows;
                int barWidth = 50;
                std::cout << "[";
                int pos = barWidth * progress;
                for (int i = 0; i < barWidth; ++i) {
                    if (i < pos) std::cout << "=";
                    else if (i == pos) std::cout << ">";
                    else std::cout << " ";
                }
                std::cout << "] " << int(progress * 100.0) << " %\r";
                std::cout.flush();
            }
        }
    }
    std::cout << std::endl; // 进度条结束后换行

    // 使用导向滤波对初始估计得到的传输函数进行细化
    cv::Mat refine_trans = guided_filter(gray_img, trans, 15, 0.000001);
    // 将细化后的传输函数 以及大气光向量代入复原公式 复原雾霾图像
    double A_b = airlight[0] * 255.0; double A_g = airlight[1] * 255.0; double A_r = airlight[2] * 255.0;
    cv::Mat dehaze_result = dehaze(image, refine_trans, A_r, A_g, A_b);

    // 记录结束时间
    end = std::chrono::high_resolution_clock::now();
    // 计算时间间隔 毫秒和秒
    auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
    // 转换为秒 精确到小数点后三位
    double duration_sec = duration.count() / 1000.0;

    // 打印耗时统计
    std::cout << "\n============ 耗时统计 ============\n";
    std::cout << "传输函数估计耗时: "
        << duration.count() << " ms ("
        << std::fixed << std::setprecision(3) << duration_sec << " s)\n";

    cv::imwrite("dehaze_result.png", dehaze_result);

    /************************************************** FOR SHOW **************************************************/
    cv::Mat trans_show;
    trans.convertTo(trans_show, CV_8UC1, 255.0, 0);
    cv::imwrite("trans.png", trans_show);

    cv::Mat refine_trans_show;
    refine_trans.convertTo(refine_trans_show, CV_8UC1, 255.0, 0);
    cv::imwrite("refine_trans.png", refine_trans_show);
    
    cv::Mat theta_E_show;
    theta_E_mat.convertTo(theta_E_show, CV_8UC1, 255.0, 0);
    cv::imwrite("theta_E.png", theta_E_show);
    /************************************************** FOR SHOW **************************************************/
    
    image.release();
    gray_img.release();
    convertImage.release(); 
    parseImage.release(); 
    trans.release();  
    refine_trans.release();
    dehaze_result.release();

    /************************************************** FOR SHOW **************************************************/
    trans_show.release();
    refine_trans_show.release();
    theta_E_show.release();
    /************************************************** FOR SHOW **************************************************/

    return 0;
}

// 大气光估计函数
double* estimate_airlight(const cv::Mat& image, int kernelSize) 
{
    int rows = image.rows, cols = image.cols;
    cv::Mat darkChannel(rows, cols, CV_8UC1);

    // 计算每个像素的最小通道值
    for (int i = 0; i < rows; i++) {
        const uchar* inData = image.ptr<uchar>(i);
        uchar* outData = darkChannel.ptr<uchar>(i);
        for (int j = 0; j < cols; j++) {
            uchar b = *inData++;
            uchar g = *inData++;
            uchar r = *inData++;
            *outData++ = std::min(r, std::min(g, b));
        }
    }

    // 使用形态学操作 最小值滤波
    cv::Mat kernel = cv::getStructuringElement(cv::MORPH_RECT, cv::Size(kernelSize, kernelSize));
    cv::erode(darkChannel, darkChannel, kernel);

    int pixelTot = static_cast<int>(rows * cols * 0.001);
    std::vector<cv::Point> allPixels(rows * cols);

    // 收集所有像素位置
    for (int r = 0; r < rows; r++) {
        const uchar* data = darkChannel.ptr<uchar>(r);
        for (int c = 0; c < cols; c++) {
            allPixels[r * cols + c] = cv::Point(c, r);
        }
    }

    // 按暗通道值降序排序
    std::sort(allPixels.begin(), allPixels.end(), [&](const cv::Point& a, const cv::Point& b) {
        return darkChannel.at<uchar>(a) > darkChannel.at<uchar>(b);
        });

    // 计算前0.1%像素的BGR平均值
    double A_b_sum = 0.0, A_g_sum = 0.0, A_r_sum = 0.0;
    for (int i = 0; i < pixelTot; i++) {
        cv::Vec3b pixel = image.at<cv::Vec3b>(allPixels[i].y, allPixels[i].x);
        A_b_sum += pixel[0];
        A_g_sum += pixel[1];
        A_r_sum += pixel[2];
    }

    static double A[3];
    A[0] = A_b_sum / (pixelTot * 255.0);
    A[1] = A_g_sum / (pixelTot * 255.0);
    A[2] = A_r_sum / (pixelTot * 255.0);

    std::cout << "Estimated Airlight: "
        << A[0] << ", " << A[1] << ", " << A[2] << std::endl;

    return A;
}

cv::Mat guided_filter(cv::Mat I, cv::Mat p, int radius, double eps) 
{
    int win_size = 2 * radius + 1;

    cv::Mat mean_I;
    cv::boxFilter(I, mean_I, CV_64FC1, cv::Size(win_size, win_size));

    cv::Mat mean_p;
    cv::boxFilter(p, mean_p, CV_64FC1, cv::Size(win_size, win_size));

    cv::Mat mean_Ip;
    cv::boxFilter(I.mul(p), mean_Ip, CV_64FC1, cv::Size(win_size, win_size));

    cv::Mat cov_Ip = mean_Ip - mean_I.mul(mean_p);

    cv::Mat mean_II;
    cv::boxFilter(I.mul(I), mean_II, CV_64FC1, cv::Size(win_size, win_size));

    cv::Mat var_I = mean_II - mean_I.mul(mean_I);
    cv::Mat a = cov_Ip / (var_I + eps);
    cv::Mat b = mean_p - a.mul(mean_I);

    cv::Mat mean_a;
    cv::boxFilter(a, mean_a, CV_64FC1, cv::Size(win_size, win_size));

    cv::Mat mean_b;
    cv::boxFilter(b, mean_b, CV_64FC1, cv::Size(win_size, win_size));

    cv::Mat q = mean_a.mul(I) + mean_b;

    return q;
}

cv::Mat dehaze(cv::Mat& src, cv::Mat& t, double A_r, double A_g, double A_b)
{
    // 确保 t 的最小值为 0.1
    cv::Mat t_clamped = cv::max(t, 0.1);

    src.convertTo(src, CV_64FC3);
    std::vector<cv::Mat> channels(3);
    cv::split(src, channels);

    cv::Mat R = channels[2];
    cv::Mat G = channels[1];
    cv::Mat B = channels[0];

    // 使用下界约束后的 t_clamped 进行计算
    channels[2] = (R - A_r) / t_clamped + A_r;
    channels[1] = (G - A_g) / t_clamped + A_g;
    channels[0] = (B - A_b) / t_clamped + A_b;

    cv::Mat dst;
    cv::merge(channels, dst);

    dst.convertTo(dst, CV_8UC3);
    return dst;
}