暗通道先驗去霧算法介紹及其幾何意義的解釋過程

         何愷明的暗通道先驗（dark channel prior）去霧算法是CV界去霧領域很有名的算法，關於該算法的論文"Single Image Haze Removal Using Dark Channel Prior"一舉獲得2009年CVPR最佳論文。作者統計了大量的無霧圖像，發現一條規律：每一幅圖像的每一個像素的RGB三個顏色通道中，總有一個通道的灰度值很低。基於這個幾乎可以視作是定理的先驗知識，作者提出暗通道先驗的去霧算法。

       作者首先介紹去霧模型如下：

如果你不是CV界新手的話，應該對上式倒背如流，在此不再對上式中的各個參量作詳細介紹。對於任意一幅輸入圖像，定義其暗通道的數學表達式為：

其中c表示rgb三通道中的某一通道。上式表示在一幅輸入圖像中，先取圖像中每一個像素的三通道中的灰度值的最小值，得到一幅灰度圖像，再在這幅灰度圖像中，以每一個像素為中心取一定大小的矩形窗口，取矩形窗口中灰度值最小值代替中心像素灰度值（最小值濾波），從而得到該霧天圖像的暗通道圖像。暗通道圖像為灰度圖像，通過大量統計並觀察發現，暗通道圖像的灰度值是很低的，所以將整幅暗通道圖像中所有像素的灰度值近似為0，即：

        作者在文中假設大氣光A為已知量，以便節省篇幅介紹傳輸函數的求解方法。在此介紹一種簡單普遍的求解大氣光的方法：對於任意一幅霧天圖像，取其暗通道圖像灰度值最大的0.1%的像素點對應於原霧天圖像的像素位置的每個通道的灰度值的平均值，從而計算出每個通道的大氣光值，即大氣光值A是一個三元素向量，每一個元素對應於每一個顏色通道。

       對於成像模型，將其歸一化，即兩邊同時除以每個通道的大氣光值：

假設在圖像中一定大小的矩形窗口內，傳輸函數的值為定值，對上式兩邊用最小化算子（minimum operators）作最小化運算：

 由於在矩形區域內為定值，故將其拿出運算符外部。由於場景輻射（scene radiance）是無霧圖像，將暗通道先驗應用於J，則有：

由於總是正值，則有：

 將上式代入到最小化運算的式子中，即可得到傳輸函數的估計值為：

為了防止去霧太過徹底，恢復出的景物不自然，應引入參數，重新定義傳輸函數為：

        對於求解得到的傳輸函數，應用濾波器對其進行細化，文章中介紹的方法是軟摳圖的方法，此方法過程復雜，速度緩慢，因此采用導向濾波對傳輸函數進行濾波。導向濾波的原理此處不再贅述，其偽代碼為:

 上述偽代碼中，I表示導向圖像（guided image）,p為輸入圖像（input image）,q為輸出圖像（output image），表示均值濾波，r為窗口半徑。

        將求解得到的去窮遠處大氣光值和傳輸函數代入到大氣散射模型之中，反解即可得到恢復的目標圖像，即場景輻射。本例使用C++實現上述去霧過程，代碼如下：

#include <iostream>
#include <stdlib.h>
#include <time.h>
#include <cmath>
#include <algorithm>
#include <opencv2\opencv.hpp>

using namespace cv;
using namespace std;

typedef struct Pixel 
{
    int x, y;
    int data;
}Pixel;

bool structCmp(const Pixel &a, const Pixel &b) 
{
    return a.data > b.data;//descending降序
}

Mat minFilter(Mat srcImage, int kernelSize);
void makeDepth32f(Mat& source, Mat& output);
void guidedFilter(Mat& source, Mat& guided_image, Mat& output, int radius, float epsilon);
Mat getTransmission_dark(Mat& srcimg, Mat& darkimg, int *array, int windowsize);
Mat recover(Mat& srcimg, Mat& t, float *array, int windowsize);

int main() 
{
    string loc = "E:\\darkchannel\\firstopencv\\sky.jpg";
    double scale = 1.0;
    string name = "forest";
    clock_t start, finish;
    double duration;

    cout << "A defog program" << endl
        << "----------------" << endl;

    Mat image = imread(loc);
    Mat resizedImage;
    int originRows = image.rows;
    int originCols = image.cols;

    if (scale < 1.0) 
    {
        resize(image, resizedImage, Size(originCols * scale, originRows * scale));
    }
    else 
    {

        scale = 1.0;
        resizedImage = image;
    }

    int rows = resizedImage.rows;
    int cols = resizedImage.cols;
    Mat convertImage;
    resizedImage.convertTo(convertImage, CV_32FC3, 1 / 255.0, 0);
    int kernelSize = 15 ? max((rows * 0.01), (cols * 0.01)) : 15 < max((rows * 0.01), (cols * 0.01));
    //int kernelSize = 15;
    int parse = kernelSize / 2;
    Mat darkChannel(rows, cols, CV_8UC1);
    Mat normalDark(rows, cols, CV_32FC1);
    int nr = rows;
    int nl = cols;
    float b, g, r;
    start = clock();
    cout << "generating dark channel image." << endl;
    if (resizedImage.isContinuous()) 
    {
        nl = nr * nl;
        nr = 1;
    }
    for (int i = 0; i < nr; i++) 
    {
        float min;
        const uchar* inData = resizedImage.ptr<uchar>(i);
        uchar* outData = darkChannel.ptr<uchar>(i);
        for (int j = 0; j < nl; j++) 
        {
            b = *inData++;
            g = *inData++;
            r = *inData++;
            min = b > g ? g : b;
            min = min > r ? r : min;
            *outData++ = min;
        }
    }
    darkChannel = minFilter(darkChannel, kernelSize);

    imshow("darkChannel", darkChannel);
    cout << "dark channel generated." << endl;

    //estimate Airlight
    //開一個結構體數組存暗通道，再sort，取最大0.1%，利用結構體內存儲的原始坐標在原圖中取點
    cout << "estimating airlight." << endl;
    rows = darkChannel.rows, cols = darkChannel.cols;
    int pixelTot = rows * cols * 0.001;
    int *A = new int[3];
    Pixel *toppixels, *allpixels;
    toppixels = new Pixel[pixelTot];
    allpixels = new Pixel[rows * cols];


    for (unsigned int r = 0; r < rows; r++) 
    {
        const uchar *data = darkChannel.ptr<uchar>(r);
        for (unsigned int c = 0; c < cols; c++) 
        {
            allpixels[r*cols + c].data = *data;
            allpixels[r*cols + c].x = r;
            allpixels[r*cols + c].y = c;
        }
    }
    std::sort(allpixels, allpixels + rows * cols, structCmp);

    memcpy(toppixels, allpixels, pixelTot * sizeof(Pixel));

    float A_r, A_g, A_b, avg, maximum = 0;
    int idx, idy, max_x, max_y;
    for (int i = 0; i < pixelTot; i++) 
    {
        idx = allpixels[i].x; idy = allpixels[i].y;
        const uchar *data = resizedImage.ptr<uchar>(idx);
        data += 3 * idy;
        A_b = *data++;
        A_g = *data++;
        A_r = *data++;
        //cout << A_r << " " << A_g << " " << A_b << endl;
        avg = (A_r + A_g + A_b) / 3.0;
        if (maximum < avg) 
        {
            maximum = avg;
            max_x = idx;
            max_y = idy;
        }
    }

    delete[] toppixels;
    delete[] allpixels;

    for (int i = 0; i < 3; i++) 
    {
        A[i] = resizedImage.at<Vec3b>(max_x, max_y)[i];
    }
    cout << "airlight estimated as: " << A[0] << ", " << A[1] << ", " << A[2] << endl;
    //cout << endl;

    //暗通道歸一化操作（除A）
    //(I / A)
    cout << "start normalization of input image I." << endl;
    float tmp_A[3];
    tmp_A[0] = A[0] / 255.0;
    tmp_A[1] = A[1] / 255.0;
    tmp_A[2] = A[2] / 255.0;
    for (int i = 0; i < nr; i++) 
    {
        float min = 1.0;
        const float* inData = convertImage.ptr<float>(i);
        float* outData = normalDark.ptr<float>(i);
        for (int j = 0; j < nl; j++) 
        {
            b = *inData++ / tmp_A[0];
            g = *inData++ / tmp_A[1];
            r = *inData++ / tmp_A[2];
            min = b > g ? g : b;
            min = min > r ? r : min;
            *outData++ = min;
        }
    }
    cout << "normalization finished." << endl << "generating relative dark channel image." << endl;
    //暗通道最小濾波
    normalDark = minFilter(normalDark, kernelSize);
    cout << "dark channel image generated." << "start estimating transmission and guided image filtering." << endl;
    imshow("normal", normalDark);
    int kernelSizeTrans = std::max(3, kernelSize);
    //求t與將t進行導向濾波

    Mat trans = getTransmission_dark(convertImage, normalDark, A, kernelSizeTrans);
    cout << "tansmission estimated and guided filtered." << endl;
    imshow("filtered t", trans);
    cout << "start recovering." << endl;
    Mat finalImage = recover(convertImage, trans, tmp_A, kernelSize);
    cout << "recovering finished." << endl;
    Mat resizedFinal;
    if (scale < 1.0) 
    {
        resize(finalImage, resizedFinal, Size(originCols, originRows));
        imshow("final", resizedFinal);
    }
    else 
    {
        imshow("final", finalImage);
    }
    finish = clock();
    duration = (double)(finish - start);
    cout << "defog used " << duration << "ms time;" << endl;
    waitKey(0);

    finalImage.convertTo(finalImage, CV_8UC3, 255);
    imwrite(name + "_refined.png", finalImage);
    destroyAllWindows();
    image.release();
    resizedImage.release();
    convertImage.release();
    darkChannel.release();
    trans.release();
    finalImage.release();
    return 0;
}

Mat minFilter(Mat srcImage, int kernelSize) 
{
    int radius = kernelSize / 2;

    int srcType = srcImage.type();
    int targetType = 0;
    if (srcType % 8 == 0) 
    {
        targetType = 0;
    }
    else 
    {
        targetType = 5;
    }
    Mat ret(srcImage.rows, srcImage.cols, targetType);
    Mat parseImage;
    copyMakeBorder(srcImage, parseImage, radius, radius, radius, radius, BORDER_REPLICATE);
    for (unsigned int r = 0; r < srcImage.rows; r++) 
    {
        float *fOutData = ret.ptr<float>(r);
        uchar *uOutData = ret.ptr<uchar>(r);
        for (unsigned int c = 0; c < srcImage.cols; c++) 
        {
            Rect ROI(c, r, kernelSize, kernelSize);
            Mat imageROI = parseImage(ROI);
            double minValue = 0, maxValue = 0;
            Point minPt, maxPt;
            minMaxLoc(imageROI, &minValue, &maxValue, &minPt, &maxPt);
            if (!targetType) 
            {
                *uOutData++ = (uchar)minValue;
                continue;
            }
            *fOutData++ = minValue;
        }
    }
    return ret;
}

void makeDepth32f(Mat& source, Mat& output)
{
    if ((source.depth() != CV_32F) > FLT_EPSILON)
        source.convertTo(output, CV_32F);
    else
        output = source;
}

void guidedFilter(Mat& source, Mat& guided_image, Mat& output, int radius, float epsilon)
{
    CV_Assert(radius >= 2 && epsilon > 0);
    CV_Assert(source.data != NULL && source.channels() == 1);
    CV_Assert(guided_image.channels() == 1);
    CV_Assert(source.rows == guided_image.rows && source.cols == guided_image.cols);

    Mat guided;
    if (guided_image.data == source.data)
    {
        //make a copy
        guided_image.copyTo(guided);
    }
    else
    {
        guided = guided_image;
    }

    //將輸入擴展為32位浮點型，以便以后做乘法
    Mat source_32f, guided_32f;
    makeDepth32f(source, source_32f);
    makeDepth32f(guided, guided_32f);

    //計算I*p和I*I
    Mat mat_Ip, mat_I2;
    multiply(guided_32f, source_32f, mat_Ip);
    multiply(guided_32f, guided_32f, mat_I2);

    //計算各種均值
    Mat mean_p, mean_I, mean_Ip, mean_I2;
    Size win_size(2 * radius + 1, 2 * radius + 1);
    boxFilter(source_32f, mean_p, CV_32F, win_size);
    boxFilter(guided_32f, mean_I, CV_32F, win_size);
    boxFilter(mat_Ip, mean_Ip, CV_32F, win_size);
    boxFilter(mat_I2, mean_I2, CV_32F, win_size);

    //計算Ip的協方差和I的方差
    Mat cov_Ip = mean_Ip - mean_I.mul(mean_p);
    Mat var_I = mean_I2 - mean_I.mul(mean_I);
    var_I += epsilon;

    //求a和b
    Mat a, b;
    divide(cov_Ip, var_I, a);
    b = mean_p - a.mul(mean_I);

    //對包含像素i的所有a、b做平均
    Mat mean_a, mean_b;
    boxFilter(a, mean_a, CV_32F, win_size);
    boxFilter(b, mean_b, CV_32F, win_size);

    //計算輸出 (depth == CV_32F)
    output = mean_a.mul(guided_32f) + mean_b;
}

Mat getTransmission_dark(Mat& srcimg, Mat& darkimg, int *array, int windowsize)
{
    //t = 1 - omega * minfilt(I / A);
    float avg_A;
    //convertImage是一個CV_32FC3的圖
    if (srcimg.type() % 8 == 0) {
        avg_A = (array[0] + array[1] + array[2]) / 3.0;
    }
    else {
        avg_A = (array[0] + array[1] + array[2]) / (3.0 * 255.0);
    }
    float w = 0.95;
    int radius = windowsize / 2;
    int nr = srcimg.rows, nl = srcimg.cols;
    Mat transmission(nr, nl, CV_32FC1);

    for (int k = 0; k<nr; k++) 
    {
        const float* inData = darkimg.ptr<float>(k);
        float* outData = transmission.ptr<float>(k);
        float pix[3] = { 0 };
        for (int l = 0; l<nl; l++)
        {
            *outData++ = 1.0 - w * *inData++;
        }
    }
    imshow("t", transmission);

    Mat trans(nr, nl, CV_32FC1);
    Mat graymat(nr, nl, CV_8UC1);
    Mat graymat_32F(nr, nl, CV_32FC1);

    if (srcimg.type() % 8 != 0) 
    {
        cvtColor(srcimg, graymat_32F, CV_BGR2GRAY);
        guidedFilter(transmission, graymat_32F, trans, 6 * windowsize, 0.001);
    }
    else 
    {
        cvtColor(srcimg, graymat, CV_BGR2GRAY);

        for (int i = 0; i < nr; i++) 
        {
            const uchar* inData = graymat.ptr<uchar>(i);
            float* outData = graymat_32F.ptr<float>(i);
            for (int j = 0; j < nl; j++)
                *outData++ = *inData++ / 255.0;
        }
        guidedFilter(transmission, graymat_32F, trans, 6 * windowsize, 0.001);
    }
    return trans;
}

Mat recover(Mat& srcimg, Mat& t, float *array, int windowsize)
{
    //J(x) = (I(x) - A) / max(t(x), t0) + A;
    int radius = windowsize / 2;
    int nr = srcimg.rows, nl = srcimg.cols;
    float tnow = t.at<float>(0, 0);
    float t0 = 0.1;
    Mat finalimg = Mat::zeros(nr, nl, CV_32FC3);
    float val = 0;

    //Be aware that transmission is a grey image
    //srcImg is a color image
    //finalImg is a color image
    //Mat store color image a pixel per 3 position
    //store grey image a pixel per 1 position
    for (unsigned int r = 0; r < nr; r++) {
        const float* transPtr = t.ptr<float>(r);
        const float* srcPtr = srcimg.ptr<float>(r);
        float* outPtr = finalimg.ptr<float>(r);
        for (unsigned int c = 0; c < nl; c++) {
            //transmission image is grey, so only need 
            //to move once per calculation, using index 
            //c(a.k.a. columns) to move is enough 
            tnow = *transPtr++;
            tnow = std::max(tnow, t0);
            for (int i = 0; i < 3; i++) {
                //so to calculate every color channel per pixel
                //move the ptr once after one calculation.
                //after 3 times move, calculation for a pixel is done
                val = (*srcPtr++ - array[i]) / tnow + array[i];
                *outPtr++ = val + 10 / 255.0;
            }
        }
    }
    cout << finalimg.size() << endl;
    return finalimg;
}

        運行代碼的測試結果為：

        可以很明顯地觀察到，所選素材中含有大量的天空區域。利用暗通道先驗規律對含有大量天空區域的圖像作暗通道圖像運算時，會發現天空區域的暗通道的灰度值並非趨向於零，因此造成在天空區域的傳輸函數估計值很低，造成誤差。為了彌補這個誤差，作者對傳輸函數的值設定了一個下界，這個下界通常取0.1，然而設定這個下界后，會發現天空區域如上述去霧后的結果那樣出現光圈（halo）效應，天空部分的去霧效果很不自然。這也是暗通道先驗去霧算法的劣處所在。

          

 

         基於暗通道先驗的去霧算法所得到的去霧結果令人驚嘆，其理論過程簡單易行，至今仍是去霧領域的經典算法。近年來，許多學者對該算法不斷從技術上改進，也有學者試圖從別的角度解釋暗通道先驗規律的成立原因，下面從幾何的角度解釋暗通道先驗去霧算法之所以成立的原因。為了與先前介紹的暗通道先驗算法中的符號區分，接下來的文中出現的特定的符號的意義如下：表示圖像中特定大小的矩陣窗口內的像素集合；表示將對應於中的像素的灰度值在RGB空間中表示的坐標的集合。

       本例介紹的論文為“Single Image Dehazing Using Color Ellipsoid Prior”，作者在這篇文章中介紹了暗通道先驗去霧算法的幾何數學解釋過程。

   無論是何種去霧方法，將去霧前后圖像的同一片區域投射在RGB坐標空間中，形成如上圖所示的點坐標的分布。將去霧前的坐標和去霧后的坐標相連，形成RGB空間的向量，將不同的去霧方法對應得到的向量定義為先驗向量(prior vector)。上圖中，通過先驗向量P去霧擴大的無霧坐標分布要比Q擴大的要小，這意味着通過先驗向量P去霧得到的結果保持了飽和度，卻喪失了對比度。反之，通過先驗向量Q去霧得到的結果保持了對比度，卻產生了過度飽和的現象,所以通過先驗向量Q的去霧結果中的樹葉部分出現黑點(dark spots)。所有的去霧算法歸結為設定合適的先驗向量，求解合適的傳輸函數，滿足：盡可能地增大圖像對比度，同時抑制過度飽和。

        上圖表示霧天圖像兩個相近區域的像素灰度值分布情況，其中一區域不含有不規則像素值(irregular pixel values)，另一區域則含有不規則像素值。不規則像素值的分布是零星的，兩個區域的大部分像素值的分布區域幾乎一致，所以兩個區域的去霧結果大體是相同的。然而在選取先驗向量時，若從不規則像素值區域選取的話，則會產生很大的誤差，因此設計的先驗向量應該具有數據魯棒性(statistically robust)以克服不規則像素值的干擾。

 

        如上圖所示，霧天圖像的局部區域中含有樹葉部分和背景部分，樹葉部分和背景部分的像素灰度值的RGB空間中如圖所示標出。在兩個部分的邊緣地帶，有大量的像素灰度值是由兩個區域混合的，這就造成選取混合區域的先驗向量進行去霧時，會在邊緣地帶造成光圈效應(halo effect)。后續算法會解釋如何去除這種光圈效應。

          下面介紹顏色橢球先驗(color ellipsoid prior)：在RGB三維空間中，含霧圖像的向量往往聚集在一起並有序的排列着。在三維空間中，向量可以用空間中的坐標表示，所以向量和顏色空間中的點是等價的。向量聚集區域可以統計地看做是一個橢球，橢球的形狀由向量的聚集狀態決定。這種方法的缺點是無法統計隨意分布的向量區域。

 

   如圖（a）和圖（b）所示，其為圖像區域對應的向量分布區域的示意圖。前者為沒有隨機分布點的規則橢球向量分布，后者是含有不規則向量的橢球向量分布，但是兩者可以視為同等橢球向量分布。

       現在根據橢球模型構造數學表達式，設為RGB空間的向量變量，對應於圖像中像素區域（i為序數標記）的向量橢球區域可以表示為：

 

 其中，表示霧天圖像在區域內的灰度值的均值對應的向量，其表達式為：

其中為霧天圖像的歸一化像素值，為三維向量。為無窮遠處的大氣光值，此處假設其已知：

方差定義為：

在RGB空間中，均值決定了球的中心位置，方差決定了球的方向和大小。

        上述各式為橢球模型的基本定義的參量，現在分析CEP值和傳輸函數的計算過程。文中有一段文字“The proposed method calculates the vector having the smallest color component among the vectors in the ellipsoid  to maximally stretch the hazy vector cluster while preventing the color components of the majority of streched vectors from being negative”,這段話的意思是，通過先驗向量將原來的圖像霧天的像素灰度值分布的橢球拉伸映射到灰度值很低的橢球中（無霧圖像的像素灰度值相對於原霧天圖像是比較低的）。故在設計合適的先驗向量和先驗值（本例中為CEP值）去霧時，應滿足如前所述的條件。通過對大氣散射模型反解，可以得到目標輻射景象為：

為了使得去霧圖像，即目標輻射場景的灰度值盡可能低，則傳輸函數的值應該盡可能地大，又因為傳輸函數與先驗值的關系為，所以先驗值應該盡可能地小。所以本例中定義的先驗值，即CEP值的數學表達式為：

該式表示，在橢球中的每一個向量取其三顏色通道分量的最小值，再取這些最小值當中的最小值作為CEP值（這與暗通道先驗規律具有一致性）。

       下面從幾何的角度求解CEP值：首先定義三個坐標軸的單位向量，，分別為：，，；將定義為與單位向量垂直且通過坐標原點的平面；如此，則為向量的R，G，B某一分量的值，並且該數值表示向量對應於空間中的坐標（即在橢球面上或者橢球內的點）到平面的距離；設為使得距離最小的向量，即：

之所以有標記c，是因為表示橢球面上或橢球內部中到平面距離最小的點對應的向量，中的c的取值有三個：，所以對應的中c的取值也有三個且和中c的取值一致，分別為，，，三者分別表示橢球中到，，平面距離最小的點坐標對應的向量。

         如上所示的（a）,（b），（c）分別為，，對應在空間中的坐標點和其到平面，，的距離的示意圖。容易觀察到，到三個平面距離最近的點均在橢球面上，且在點處橢球的切面與對應的平面平行。由於點在橢球面上，則滿足表達式：

上式表示點在以為中心，為協方差矩陣（該協方差矩陣規定了橢球的形狀和方向）逆矩陣的橢球面上。則關於CEP值，有如下推導：

 由於點處橢球的切面與對應的平面平行，所以在點出的外單位法向量（outward normal vector）為上圖所示的。根據文中所述定理：The gradient direction of a vector on the ellipse surface is perpendicular to the tangent plane,and so the outward normal vector must have the same direction as the gradient of the ellipse surface，即在球面處的點的梯度方向和在該點處的切平面垂直，在該點處的外單位法向量和梯度方向是一致的。根據該定理，可得下式：

其中，梯度部分的計算過程為：

則有：

所以有：

其中，。由此重新表達為：

所以中的向量到平面距離最小值為：

所以，CEP值的表達式為：

故對應的傳輸函數的表達式為：

        下面根據上述推導過程介紹如何設計算法，設是使得CEP值在三通道中最小的那一通道，即：

則傳輸函數可以重新表示為：

若直接計算每一點鄰域矩形的均值與方差的差值，再取其三通道中的最小值較為繁瑣，作者提出一種近似方法，利用這種近似方法與原原本本的計算結果相差無幾，但是大大提高了速度和效率。近似方法為：

其中，表示輸入霧天圖像在像素點j處c通道的灰度值。

       設大氣光已知，按照上述過程計算每一像素點處的均值和方差，並將其代入傳輸函數公式中以計算傳輸函數，由此即可恢復目標圖像，處理結果為：

上述結果在邊緣部分出現了很模糊的現象，這是由於在選擇橢球區域時，所選區域的像素梯度過大，某一點的像素灰度值與其周圍像素的灰度值差別很大。因此有必要對求取的傳輸函數進行細化，作者分別對均值和方差進行細化（即濾波）：對於均值部分，作者采用了與導向濾波十分類似的濾波器，但不完全是導向濾波；對於方差部分，作者采用了均值濾波。濾波過程的偽代碼如下：

 

 

       用MATLAB對上述偽代碼的仿真為：

clear all;
I=imread('E:\picture\t.bmp');
I=double(I);
I=I./255;
dark=darkfunction(I);
[m,n,o]=size(I);
imsize=m*n;
numpx=floor(imsize/1000);
J_dark_vec=reshape(dark,imsize,1);
I_vec=reshape(I,imsize,3);
[J_dark_vec,indices]=sort(J_dark_vec);
indices=indices(imsize-numpx+1:end);
atmSum=zeros(1,3);
for ind=1:numpx
    atmSum=atmSum+I_vec(indices(ind),:);
end;
atmospheric=atmSum/numpx;
A=atmospheric;

omega=0.95;
r=10;epsilon=0.002;
epsilon=(1/mean(A))^2*epsilon;
R=I(:,:,1)/A(1);G=I(:,:,2)/A(2);B=I(:,:,3)/A(3);
Imin=min(min(R,G),B);

N=boxfilter(ones(m,n),r);
u_m=boxfilter(Imin,r)./N;
%對均值濾波，與導向濾波相似，但不是導向濾波
u_mm=boxfilter(Imin.*Imin,r)./N;
sigma_m=u_mm-u_m.*u_m;
a=sigma_m./(sigma_m+epsilon);
b=(1-a).*u_m;
u_a=boxfilter(a,r)./N;
u_b=boxfilter(b,r)./N;
u_m=u_a.*Imin+u_b;

sigma_mm=boxfilter((Imin-u_m).*(Imin-u_m),r)./N;
%對方差濾波，此處為均值濾波，也可仿上述濾波過程進行濾波
sigma_mm=boxfilter(sigma_mm,r)./N;

t=1-omega*(u_m-sqrt(abs(sigma_mm)));

for i=1:m
    for j=1:n
        for k=1:3
            I(i,j,k)=(I(i,j,k)-A(k))/max(t(i,j),0.1)+A(k);
        end;
    end;
end;
imshow(I);

function [dark] =darkfunction(I)
R=I(:,:,1);
G=I(:,:,2);
B=I(:,:,3);
[m,n]=size(R);
a=zeros(m,n);
for i=1:m
    for j=1:n
        a(i,j)=min(R(i,j),G(i,j));
        a(i,j)=min(a(i,j),B(i,j));
    end;
end;
d=ones(15,15);
fun=@(block_struct)min(min(block_struct.data))*d;
dark=blockproc(a,[15,15],fun);
dark=dark(1:m,1:n);
end

function [J] = boxfilter(I,r)
[m,n]=size(I);
J=zeros(size(I));
Icum=cumsum(I,1);
J(1:r+1,:)=Icum(r+1:2*r+1,:);
J(r+2:m-r,:)=Icum(2*r+2:m,:)-Icum(1:m-2*r-1,:);
J(m-r+1:m,:)=repmat(Icum(m,:),[r,1])-Icum(m-2*r:m-r-1,:);
Icum=cumsum(J,2);
J(:,1:r+1)=Icum(:,r+1:2*r+1);
J(:,r+2:n-r)=Icum(:,2*r+2:n)-Icum(:,1:n-2*r-1);
J(:,n-r+1:n)=repmat(Icum(:,n),[1,r])-Icum(:,n-2*r:n-r-1);
end

        用C++實現算法為：

#include <iostream>
#include <stdlib.h>
#include <time.h>
#include <cmath>
#include <algorithm>

#include <opencv2\opencv.hpp>
#include <opencv2\imgproc\imgproc.hpp>
#include <opencv2\highgui\highgui.hpp>

using namespace std;
using namespace cv;

int main()
{
    Mat GF_smooth(Mat& src, int s, double epsilon);
    Mat staticMin(Mat& I, int s, double eeps, double alpha);
    int est_air(Mat& src, int s, double *A_r, double *A_g, double *A_b);
    Mat est_trans_fast(Mat& src, int s, double eeps, double k, double A_r, double A_g, double A_b);
    Mat rmv_haze(Mat& src, Mat& t, double A_r, double A_g, double A_b);

    string loc = "E:\\sphere\\firstopencv\\t.bmp";
    string name = "forest";
    clock_t start, finish;
    double duration;

    cout << "A defog program" << endl
         << "---------------" << endl;

    start = clock();
        Mat image = imread(loc);
    imshow("hazyimage", image);
    cout << "input hazy image" << endl;
    
    int s = 15;
    double eeps = 0.002, omega = 0.95;
    double air_r, air_g, air_b;
    est_air(image, s, &air_r, &air_g, &air_b);
    cout << "airlight estimation as:" << air_r << ", " << air_g << ", " << air_b << endl;

    Mat t;
    t=est_trans_fast(image, s, eeps, omega, air_r, air_g, air_b);

    Mat dehazedimage;
    dehazedimage = rmv_haze(image, t, air_r, air_g, air_b);
    imshow("dehaze", dehazedimage);
    finish = clock();
    duration = (double)(finish - start);
    cout << "defog used" << duration << "ms time;" << endl;
    waitKey(0);

    imwrite(name + "_refined.png", dehazedimage);
    destroyAllWindows();
    image.release();
    dehazedimage.release();
    return 0;
}

//---------------------- GUIDED FILTER -------------------//
Mat GF_smooth(Mat& src, int s, double epsilon)
{
    src.convertTo(src, CV_64FC1);

    Mat mean_I;
    blur(src, mean_I, Size(s, s), Point(-1, -1));

    Mat II = src.mul(src);
    Mat var_I;
    blur(II, var_I, Size(s, s), Point(-1, -1));
    var_I = var_I - mean_I.mul(mean_I);

    Mat a = var_I / ((var_I + epsilon));
    Mat b = mean_I - a.mul(mean_I);

    Mat mean_a;
    blur(a, mean_a, Size(s, s), Point(-1, -1));
    Mat mean_b;
    blur(b, mean_b, Size(s, s), Point(-1, -1));

    return mean_a.mul(src) + mean_b;
}


Mat staticMin(Mat& I, int s, double eps, double alpha)
{
    Mat mean_I = GF_smooth(I, s, eps);

    Mat var_I;
    blur((I - mean_I).mul(I - mean_I), var_I, Size(s, s), Point(-1, -1));

    Mat mean_var_I;
    blur(var_I, mean_var_I, Size(s, s), Point(-1, -1));

    Mat z_I;
    sqrt(mean_var_I, z_I);

    return mean_I - alpha*z_I;
}


//---------------------- DEHAZING FUNCTIONS -------------------//
int est_air(Mat& src, int s, double *A_r, double *A_g, double *A_b)
    {
        /// Split RGB to channels
        src.convertTo(src, CV_64FC3);
        vector<Mat> channels(3);
        split(src, channels);        // separate color channels

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

        Mat Im = min(min(R, G), B);

        /// Estimate airlight
        Mat blur_Im;
        blur(Im, blur_Im, Size(s, s), Point(-1, -1));

        int maxIdx[2] = { 0, 0 };
        minMaxIdx(blur_Im, NULL, NULL, NULL, maxIdx);

        int width = R.cols;
        *A_r = ((double*)R.data)[maxIdx[0] * R.cols + maxIdx[1]];
        *A_g = ((double*)G.data)[maxIdx[0] * R.cols + maxIdx[1]];
        *A_b = ((double*)B.data)[maxIdx[0] * R.cols + maxIdx[1]];

        return 0;
    }

    Mat est_trans_fast(Mat& src, int s, double eeps, double k, double A_r, double A_g, double A_b)
    {
        /// Split RGB to channels
        src.convertTo(src, CV_64FC3);
        vector<Mat> channels(3);
        split(src, channels);        // separate color channels

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

        /// Estimate transmission
        Mat R_n = R / A_r;
        Mat G_n = G / A_g;
        Mat B_n = B / A_b;

        Mat Im = min(min(R_n, G_n), B_n);

        eeps = (3 * 255 / (A_r + A_g + A_b))*(3 * 255 / (A_r + A_g + A_b))*eeps;
        double alpha = 2;
        Mat z_Im = staticMin(Im, s, eeps, alpha);

        return min(max(0.001, 1 - k*z_Im), 1);
    }


    Mat rmv_haze(Mat& src, Mat& t, double A_r, double A_g, double A_b)
    {
        /// Split RGB to channels
        src.convertTo(src, CV_64FC3);
        vector<Mat> channels(3);
        split(src, channels);        // separate color channels

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

        /// Remove haze
        channels[2] = (R - A_r) / t + A_r;
        channels[1] = (G - A_g) / t + A_g;
        channels[0] = (B - A_b) / t + A_b;

        Mat dst;
        merge(channels, dst);

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

        用MATLAB或者C++運行的實驗結果為：

        結果分析：通過該實驗結果發現，顏色橢球先驗（CEP）的去霧算法的運行結果和暗通道先驗（DCP）算法的運行結果基本上是一致的，這也驗證了CEP是DCP在幾何上的解釋，兩者具有內在一致性，在本質上是相同的。但是，兩者也有細微的差別。由於兩種算法的實現結果一致，不好比較，現在將兩者在傳輸函數優化之前的去霧效果進行比較，下圖為DCP所求的傳輸函數在導向濾波之前的去霧效果：

下圖為CEP所求傳輸函數在濾波優化前的去霧效果：

很明顯地觀察到：DCP的邊緣效應比CEP的邊緣效應更明顯一些，說明基於CEP計算的傳輸函數要比基於DCP計算的傳輸函數更精確一些。這也體現了文章中所說的，與之前的DCP相比較，CEP去霧效果更具有數據魯棒性（statistically robust）。