基于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> _centroid;
std::vector<double> _radii;
std::vector<std::vector<double>> _pose;
long int _iters;

void decompose(Eigen::MatrixXd& A,
Eigen::VectorXd& s,
Eigen::MatrixXd& U,
Eigen::MatrixXd& V);

int findIdx(Eigen::VectorXd& vec, double val);

public:
Mvee();

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

std::vector<double> centroid() const;
std::vector<double> radii() const;
std::vector<std::vector<double>> pose() 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& A,
        Eigen::VectorXd& s,
        Eigen::MatrixXd& U,
        Eigen::MatrixXd& V) {
        Eigen::JacobiSVD<Eigen::MatrixXd> SVD(A,
            Eigen::ComputeThinU | Eigen::ComputeThinV);
        s = SVD.singularValues();
        V = SVD.matrixV();
        U = SVD.matrixU();
    }

    int Mvee::findIdx(Eigen::VectorXd& vec, double val) {
        for (int i = 0; i < vec.rows(); i++) {
            if (vec(i) == val) {
                return i;
            }
        }
        return -1;
    }

    void Mvee::compute(Eigen::MatrixXd& data,
        double eps,
        double lmcoeff) {
        _iters = 0;
        double err = 1 + eps;
        double alpha;
        int n = data.rows();
        int d = data.cols();
        Eigen::MatrixXd X = data.transpose();
        Eigen::MatrixXd Q(X.rows() + 1, X.cols());
        Q.row(0) = X.row(0);
        Q.row(1) = X.row(1);
        Q.row(2).setOnes();

        Eigen::VectorXd u = (1.0 / n) * Eigen::VectorXd::Ones(n);
        Eigen::VectorXd uhat = u;
        Eigen::MatrixXd G(d + 1, d + 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() * Q).diagonal();
            m = g.maxCoeff(&idx);
            alpha = (m - d - 1) / ((d + 1) * (m - 1));
            uhat = (1 - alpha) * u;
            uhat(idx) += alpha;
            err = (uhat - u).norm();
            u = uhat;
            _iters++;
        }

        Eigen::MatrixXd Einv = (1.0 / d) * (
            X * u.asDiagonal() * X.transpose() -
            (X * u) * (X * u).transpose()
            );

        Eigen::VectorXd x = X * u;
        Eigen::VectorXd s;
        Eigen::MatrixXd V_mat;
        Eigen::MatrixXd U_mat;
        decompose(Einv, s, U_mat, V_mat);

        // 存储结果
        _centroid.resize(x.size());
        for (int i = 0; i < x.size(); i++) {
            _centroid[i] = x(i);
        }

        _radii.resize(s.size());
        for (int i = 0; i < s.size(); i++) {
            _radii[i] = 1.0 / std::sqrt(s(i));
        }

        _pose.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++) {
                _pose[i][j] = V_mat(i, j);
            }
        }
    }

    std::vector<double> Mvee::centroid() const {
        return _centroid;
    }

    std::vector<double> Mvee::radii() const {
        return _radii;
    }

    std::vector<std::vector<double>> Mvee::pose() const {
        return _pose;
    }

    long int Mvee::iters() const {
        return _iters;
    }

    Mvee::~Mvee() {}

} // namespace mvee

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

#include "mvee.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-18

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);
cv::Mat flatten_transform(cv::Mat& input);

int main() {
    // 读取图像
    std::string pic_name = "pumpkins.png";

    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." << 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);
    for (int r = 0; r < rows; r++) {
        double* tPtr = trans.ptr<double>(r);
        for (int c = 0; c < cols; c++) {
            // ROI（Region of Interest）
            cv::Rect roi(c, r, patchSize, patchSize);
            cv::Mat image_roi = parseImage(roi);
            cv::Mat image_roi_flatten = flatten_transform(image_roi);

            // 使用Mvee类计算最小包围椭球
            Eigen::MatrixXd region_roi;
            cv::cv2eigen(image_roi_flatten, region_roi);

            mvee::Mvee Khachiyan_solver;
            Khachiyan_solver.compute(region_roi, EPS, LMC);

            std::vector<double> center = Khachiyan_solver.centroid();
            std::vector<std::vector<double>> E = Khachiyan_solver.pose();
            std::vector<double> radii = Khachiyan_solver.radii();

            double thetaVal;
            // 检查计算结果是否有效
            bool valid = true;
            for (double val : center) {
                if (std::isnan(val)) {
                    valid = false;
                    break;
                }
            }
            if (valid) {
                for (const auto& row : E) {
                    for (double val : row) {
                        if (std::isnan(val)) {
                            valid = false;
                            break;
                        }
                    }
                    if (!valid) break;
                }
            }

            if (!valid) {
                thetaVal = 0.1;
            }
            else {
                // 修复：使用原始代码中的计算逻辑
                double d1 = radii[0]; // 原始代码使用第一个特征值
                std::vector<double> v1 = { E[0][0], E[1][0], E[2][0] }; // 第一列特征向量

                // 计算系数
                double coeffi = 1.0 / d1; // 因为radii[i] = 1/sqrt(s(i)), 所以sqrt(s(i)) = 1/radii[i]

                // 计算g向量
                std::vector<double> g(3);
                g[0] = center[0] - coeffi * v1[0];
                g[1] = center[1] - coeffi * v1[1];
                g[2] = center[2] - coeffi * v1[2];

                // 计算thetaVal
                double norm = std::sqrt(g[0] * g[0] + g[1] * g[1] + g[2] * g[2]);
                thetaVal = norm / std::sqrt(3.0);
            }

            double omega = 0.95;
            *tPtr = 1 - omega * thetaVal;
            tPtr++;
        }
    }

    // 使用导向滤波对初始估计得到的传输函数进行细化
    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_show.png", trans_show);

    cv::Mat refine_trans_show;
    refine_trans.convertTo(refine_trans_show, CV_8UC1, 255.0, 0);
    cv::imwrite("refine_trans_show.png", refine_trans_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();
    /************************************************** 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;
}

cv::Mat flatten_transform(cv::Mat& input) {
    int rows = input.rows; int cols = input.cols; int length = rows * cols;
    cv::Mat output = cv::Mat::zeros(length, 3, CV_64FC1);

    for (int ker_row = 0; ker_row < rows; ker_row++) {
        for (int ker_col = 0; ker_col < cols; ker_col++) {
            cv::Vec3d pixel = input.at<cv::Vec3d>(ker_row, ker_col);
            output.at<double>(ker_row * cols + ker_col, 0) = pixel[0];
            output.at<double>(ker_row * cols + ker_col, 1) = pixel[1];
            output.at<double>(ker_row * cols + ker_col, 2) = pixel[2];
        }
    }

    return output;
}

          下面给出该代码的多组运行结果：

   

    

            附录（Khachiyan算法）