幾種去霧算法介紹

Single Image Dehazing

Raanan Fattal

Hebrew University of Jerusalem,Israel

      這篇文章提出一種新的從單幅輸入圖像中估計傳輸函數的方法。新方法中，重新定義了大氣傳輸模型，大氣散射模型中除了傳輸函數（transmission function）這個變量外，還增加了表面陰影（surface shading）這個變量。作者假設一個前提，表面陰影和傳輸函數是統計無關的，根據這一前提對大氣散射模型進行運算分析，即可求得傳輸函數並對圖像去霧。

      作者首先介紹了大氣散射模型：

      該式定義域RGB三顏色通道空間，表示探測系統獲取的圖像，是無窮遠處的大氣光，表示目標輻射光，即需要回復的目標圖像，表示傳輸函數，即光在散射介質中傳輸經過衰減等作用能夠到達探測系統的那一部分的光的比例。坐標向量表示探測系統獲取的圖像中每一個像素的坐標位置。

      對大氣散射模型進行變形，將需要恢復的目標圖像視作表面反射系數（surface albedo coefficients）和陰影系數（shading factor）的按坐標的點乘，即，其中為三通道向量，是描述在表面反射的光的標量。即的尺度與相同，為彩色圖像，為灰度圖像。為了簡化，假設在某區域內為常數，即在像素區域內，為常數。則大氣散射模型變為：

       將向量分解成兩個部分，一部分為與大氣光平行的向量，另一部分為與大氣光垂直的殘留向量（residual vector），記作，且，表示與大氣光向量垂直的所有向量構成的向量空間。如圖所示，向量與向量之間的夾角為，從的端點引垂線到，垂足與的端點的連線即為向量。

       對於重新定義的大氣散射模型中的，將其寫成平行於的向量於平行於的向量之和：

         其中，記作，為表面反射和大氣光的相關量或相關系數，表示在RGB空間中的兩個三維向量的點積。

       為了獲得獨立的方程，求取輸入圖像沿着大氣光向量的那一分量（標量）為：

 

則輸入圖像沿着方向的那一分量（標量）為：

（因為向量和向量垂直，所以） 。則有：

 

 由上式解得傳輸函數為：

         若已知無窮遠出的大氣光，則與均可求，唯一未知量為，所以求解的問題就歸結為求解內的問題。

       由於傳輸函數，所以傳輸函數與散射系數和景深有關，而表面陰影與場景的光照（illuminance in the scene）、表面反射屬性（surface reflectance properties）和場景幾何結構（scene geometry）有關。所以假設，陰影函數和傳輸函數不具有局部相關性，用協方差表示這種無關性為：

其中表示為在區域內兩變量的協方差（covariance），協方差的計算公式為：

均值的計算公式為：

 

為了使計算簡便，考慮和的協方差無關性，即通過解出的表達式。重新表達和為：

上述兩式中，除了參數和為常量外，其余參數均為變量，式中定義為：

根據協方差公式的性質，和（a為常量），可以得到：

所以有，該式中由於和均為常量，所以可得，即，從而得到：

將解得的代入到傳輸函數的表達式中，即可解析去霧模型中的參量。

        本例中為了方便計算，所選塊狀區域的大小為整幅輸入圖像的尺寸；本文注重介紹傳輸函數的求解方法，所以無窮遠處大氣光的求解可以參考暗通道先驗模型進行求解；最后回復出的無霧圖像需要進行一次線性拉伸，才能顯示出去霧結果。本實驗的C++代碼如下：

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

using namespace cv;
using namespace std;

float sqr(float x);
float norm(float *array);
float avg(float *vals, int n);
float conv(float *xs, float *ys, int n);
Mat stress(Mat& input);
Mat getDehaze(Mat& scrimg, Mat& transmission, float *array);
Mat getTransmission(Mat& input, float *airlight);

int main()
{
    string loc = "E:\\fattal\\project\\project\\house.bmp";
    double scale = 1.0;

    clock_t start, finish;
    double duration;

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

    Mat image = imread(loc);
    imshow("hazyiamge", image);
    cout << "input hazy image" << endl;
    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;
    }

    start = clock();
    int rows = resizedImage.rows;
    int cols = resizedImage.cols;
    int nr = rows; int nl = cols;
    Mat convertImage(nr, nl, CV_32FC3);
    resizedImage.convertTo(convertImage, CV_32FC3, 1 / 255.0, 0);
    int kernelSize = 15;

    float tmp_A[3];
    tmp_A[0] = 0.84;
    tmp_A[1] = 0.83;
    tmp_A[2] = 0.80;

    Mat trans = getTransmission(convertImage, tmp_A);
    cout << "tansmission estimated." << endl;
    imshow("t", trans);

    cout << "start recovering." << endl;
    Mat finalImage = getDehaze(convertImage, trans, tmp_A);
    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);
    const char* path;
    path = "C:\\Users\\Administrator\\Desktop\\recover.jpg";
    imwrite(path, finalImage);
    destroyAllWindows();
    image.release();
    resizedImage.release();
    convertImage.release();
    trans.release();
    finalImage.release();
    return 0;
}

float sqr(float x)
{
    return x * x;
}

float norm(float *array)
{
    return sqrt(sqr(array[0]) + sqr(array[1]) + sqr(array[2]));
}

float avg(float *vals, int n)
{
    float sum = 0;
    for (int i = 0; i < n; i++)
    {
        sum += vals[i];
    }
    return sum / n;
}

float conv(float *xs, float *ys, int n)
{
    float ex = avg(xs, n);
    float ey = avg(ys, n);
    float sum = 0;
    for (int i = 0; i < n; i++)
    {
        sum += (xs[i] - ex)*(ys[i] - ey);
    }
    return sum / n;
}

Mat getDehaze(Mat& scrimg, Mat& transmission, float *array)
{
    int nr = transmission.rows; int nl = transmission.cols;
    Mat result = Mat::zeros(nr, nl, CV_32FC3);
    Mat one = Mat::ones(nr, nl, CV_32FC1);
    vector<Mat> channels(3);
    split(scrimg, channels);

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

    channels[2] = (R - (one - transmission)*array[2]) / transmission;
    channels[1] = (G - (one - transmission)*array[1]) / transmission;
    channels[0] = (B - (one - transmission)*array[0]) / transmission;

    merge(channels, result);
    return result;
}

Mat getTransmission(Mat& input, float *airlight)
{
    float normA = norm(airlight);
    //Calculate Ia
    int nr = input.rows, nl = input.cols;
    Mat Ia(nr, nl, CV_32FC1);
    for (int i = 0; i < nr; i++)
    {
        const float* inPtr = input.ptr<float>(i);
        float* outPtr = Ia.ptr<float>(i);
        for (int j = 0; j < nl; j++)
        {
            float dotresult = 0;
            for (int k = 0; k < 3; k++)
            {
                dotresult += (*inPtr++)*airlight[k];
            }
            *outPtr++ = dotresult / normA;
        }
    }
    imshow("Ia", Ia);

    //Calculate Ir
    Mat Ir(nr, nl, CV_32FC1);
    for (int i = 0; i < nr; i++)
    {
        Vec3f* ptr = input.ptr<Vec3f>(i);
        float* irPtr = Ir.ptr<float>(i);
        float* iaPtr = Ia.ptr<float>(i);
        for (int j = 0; j < nl; j++)
        {
            float inNorm = norm(ptr[j]);
            *irPtr = sqrt(sqr(inNorm) - sqr(*iaPtr));
            iaPtr++; irPtr++;
        }
    }
    imshow("Ir", Ir);

    //Calculate h
    Mat h(nr, nl, CV_32FC1);
    for (int i = 0; i < nr; i++)
    {
        float* iaPtr = Ia.ptr<float>(i);
        float* irPtr = Ir.ptr<float>(i);
        float* hPtr = h.ptr<float>(i);
        for (int j = 0; j < nl; j++)
        {
            *hPtr = (normA - *iaPtr) / *irPtr;
            hPtr++; iaPtr++; irPtr++;
        }
    }
    imshow("h", h);

    //Estimate the eta
    int length = nr * nl;
    float* Iapix = new float[length];
    float* Irpix = new float[length];
    float* hpix = new float[length];
    for (int i = 0; i < nr; i++)
    {
        const float *IaData = Ia.ptr<float>(i);
        const float *IrData = Ir.ptr<float>(i);
        const float *hData = h.ptr<float>(i);
        for (int j = 0; j < nl; j++)
        {
            Iapix[i*nl + j] = *IaData++;
            Irpix[i*nl + j] = *IrData++;
            hpix[i*nl + j] = *hData++;
        }
    }
    float eta = conv(Iapix, hpix, length) / conv(Irpix, hpix, length);
    cout << "the value of eta is:"  << eta << endl;

    //Calculate the transmission
    Mat t(nr, nl, CV_32FC1);
    for (int i = 0; i < nr; i++)
    {
        float* iaPtr = Ia.ptr<float>(i);
        float* irPtr = Ir.ptr<float>(i);
        float* tPtr = t.ptr<float>(i);
        for (int j = 0; j < nl; j++)
        {
            *tPtr = 1 - (*iaPtr - eta * (*irPtr)) / normA;
            tPtr++; iaPtr++; irPtr++;
        }
    }
    imshow("t1", t);
    Mat trefined;
    trefined = stress(t);
    return trefined;
}

Mat stress(Mat& input)
{
    float data_max = 0.0, data_min = 5.0;
    int nr = input.rows; int nl = input.cols;
    Mat output(nr, nl, CV_32FC1);
    for (int i = 0; i < nr; i++)
    {
        float* data = input.ptr<float>(i);
        for (int j = 0; j < nl; j++)
        {
            if (*data > data_max) data_max = *data;
            if (*data < data_min) data_min = *data;
            data++;
        }
    }
    float temp = data_max - data_min;
    for (int i = 0; i < nr; i++)
    {
        float* indata = input.ptr<float>(i);
        float* outdata = output.ptr<float>(i);
        for (int j = 0; j < nl; j++)
        {
            *outdata = (*indata - data_min) / temp;
            if (*outdata < 0.1) *outdata = 0.1;
            indata++; outdata++;
        }
    }
    return output;
}

        對於上述代碼，由於輸入霧天圖像不同，調整對應霧天圖像的無窮遠處大氣光值即可。如下兩圖分別為霧天圖像與去霧結果圖像：

                                                                        

 

 A fast single image haze removal algorithm using color attenuation prior

Qingsong Zhu, Jiaming Mai, Ling Shao

       這篇文章從統計的角度出發，提出顏色衰減先驗的規律（Color Attenuation Prior），將傳輸函數的景深以解析解的形式給出，而非將傳輸函數看作是一個整體，具有很高的創意。

       文章首先介紹霧天圖像的成像模型，即大氣散射模型：

式中，表示霧天圖像，表示無霧圖像，表示大氣光向量，表示在像素處的傳輸函數值，為散射系數，為景深，表示場景距離成像系統的距離。

       對於大氣散射模型中的大氣光向量，作者的估計方法為：

 本文主要介紹求解景深，當景深求解得到后，即可依據景深求解霧天圖像的大氣光向量：設定一個關於景深的閾值，即某個像素的景深距離達到一定程度時，該處的像素灰度值即為大氣光向量。

         下面介紹顏色衰減先驗模型：作者做了一個實驗，在同一幅霧天圖像選取三個區域，即高濃度霧區域（景深較大區域），中等濃度霧區域（中等景深區域）和低濃度霧區域（近景區域），分析景深和亮度與飽和度之間的關系。

         如上圖所示，c區域為近景無霧部分，b區域遠近適中，有一定霧，a區域為遠景霧氣濃重。由c右邊圖表可知，c圖的亮度46.18%，飽和度很高，兩者幾乎沒差。通過圖表的比較可以得出以下幾個結論：1.飽和度很容易受霧氣的影響，一旦有一點霧，下降的很快。2.亮度在有霧的情況下（有散熱光)反而會更亮。3.無霧情況下，亮度和飽和度幾乎沒差，受霧的影響下，亮度和飽和度之差懸殊。並且霧越濃重，兩者相差越懸殊，也就是說亮度和飽和度之差和霧濃度正相關。

      基於此，可以得到下面的數學關系：

其中，表示在像素處的霧氣濃度，和分別表示場景的亮度和飽和度。

         根據該數學關系，將霧天圖像從RGB顏色空間轉換到HSV顏色空間，將景深表示為關於HSV空間中的S分量和V分量：

其中，，與為待求系數，為隨機噪聲分布，本例中其為高斯噪聲，即其均值為0，方差為。

       根據高斯分布的性質，於是有：

        為了求解未知參數，作者通過合成霧天圖像數據集的方法求解： 清晰無霧圖+隨機深度圖（服從（0,1）均勻分布）+樣本大氣光=樣本圖像。首先，對每一張清晰圖像，生成一張相同大小的隨機深度圖，合成的深度圖服從在開區間（0,1）上的標准均勻分布；然后，生成（0.85,1）之間的隨機大氣光。將隨機生成的深度度和大氣光結合，即可得到合成的霧天圖像。本次實驗選取500幅圖像作為實驗樣本。

      對於聯合概率密度函數，我們假設每個像素處的概率是相互獨立的，構造最大似然函數如下：

 由於每個像素處均服從正太分布，則有：

其中為在第個場景點處的真實景深，所以最大化以求解參數，，與。為了計算簡便，對取自然對數：

       對上式關於求導並令其為0，可得：

       為了求解參數，，，采用梯度下降法，對分別關於，，求導：

則迭代計算參數，，如下：

       上述算法的算法流程如下：

        該實驗用了500個訓練樣本和1.2億像素點訓練線性模型，經過517次迭代得出系數：，，和。確定了相關系數，就可以恢復霧天圖像的深度。

        在得到景深圖后，需要進一步優化，作者首先對景深圖像作最小值濾波，然后再對其進行導向濾波，從而得到優化的景深圖像。

        將估計得到的景深代入到傳輸函數及成像模型之中，即可復原出無霧圖像：

為了計算簡便，散射系數設置為定值1.0。

         下面給出上述算法的C++代碼：

#include <iostream>
#include <stdlib.h>
#include <time.h>
#include <cmath>
#include <algorithm>
#include <random>
#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 getDepthmap(Mat& input, float theta0, float theta1, float theta2);
Mat normrnd(float aver, float sigma, int row, int col);
Mat minFilter(Mat& input, int kernelSize);
Mat guidedfilter(Mat& srcImage, Mat& srcClone, int r, double eps);
Mat recover(Mat& srcimg, Mat& t, float *array);

int main() 
{
    string loc = "E:\\code\\CAP\\project\\project\\cones.bmp";
    double scale = 1.0;
    clock_t start, finish;
    double duration;

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

    Mat image = imread(loc);
    imshow("hazyiamge", image);
    cout << "input hazy image" << endl;
    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;
    }

    start = clock();
    int rows = resizedImage.rows;
    int cols = resizedImage.cols;
    Mat convertImage; 
    resizedImage.convertTo(convertImage, CV_32FC3, 1 / 255.0, 0);
    int kernelSize = 7; float eps = 0.0001;
    int radius = kernelSize / 2;

    Mat depmap(rows, cols, CV_32FC1);
    float the0 = 0.121779, the1 = 0.959710, the2 = -0.780245;
    float aveg = 0.0; float sigma = 0.041337;
    depmap = getDepthmap(convertImage, the0, the1, the2);
    Mat guassian = normrnd(aveg, sigma, rows, cols);
    depmap += guassian;

    imshow("depmap", depmap);

    Mat refdep = minFilter(depmap, kernelSize);
    Mat finaldep(rows, cols, CV_32FC1);
    Mat graymat(rows, cols, CV_8UC1);
    Mat graymat_32F(rows, cols, CV_32FC1);

    cvtColor(image, graymat, CV_BGR2GRAY);

    for (int i = 0; i < rows; i++)
    {
        const uchar* inData = graymat.ptr<uchar>(i);
        float* outData = graymat_32F.ptr<float>(i);
        for (int j = 0; j < cols; j++)
        {
            *outData++ = *inData++ / 255.0;
        }
    }

    float epsilon = 0.0001;
    finaldep = guidedfilter(image, refdep, 6 * kernelSize, epsilon);
    
    //estimate Airlight
    cout << "estimating airlight." << endl;
    rows = depmap.rows, cols = depmap.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 = depmap.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++;
        
        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;

    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;

    float beta = 1.0;
    Mat trans;
    cv::exp(-beta * finaldep, trans);
    cout << "tansmission estimated." << endl;
    imshow("t", trans);

    cout << "start recovering." << endl;
    Mat finalImage = recover(convertImage, trans, tmp_A);
    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("refined.bmp", finalImage);
    destroyAllWindows();
    image.release();
    resizedImage.release();
    convertImage.release();
    trans.release();
    finalImage.release();
    return 0;
}

Mat getDepthmap(Mat& input, float theta0, float theta1,float theta2)
{
    Mat Ihsv; Mat output; Mat depth;

    cv::cvtColor(input, Ihsv, cv::COLOR_BGR2HSV);
    std::vector<cv::Mat>hsv_vec;
    cv::split(Ihsv, hsv_vec);
    cv::addWeighted(hsv_vec[1], theta2, hsv_vec[2], theta1, theta0, output);
    depth = output;

    return depth;
}

Mat normrnd(float aver, float sigma, int row, int col)
{
    Mat p(row, col, CV_32FC1);
    
    random_device rd;
    mt19937 gen(rd());
    
    for (int i = 0; i < row; i++)
    {
        float *pData = p.ptr<float>(i);
        for (int j = 0; j < col; j++)
        {
            normal_distribution<float> normal(aver, sigma);
            *pData = normal(gen);
            pData++;
        }
    }
    return p;
}
 
Mat minFilter(Mat& input, int kernelSize)
{
    int row = input.rows; int col = input.cols;
    int radius = kernelSize / 2;
    Mat parseImage;
    copyMakeBorder(input, parseImage, radius, radius, radius, radius, BORDER_REPLICATE);
    Mat output = Mat::zeros(input.rows, input.cols, CV_32FC1);

    for (unsigned int r = 0; r < row; r++)
    {
        float *fOutData = output.ptr<float>(r);

        for (unsigned int c = 0; c < col; 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);

            *fOutData++ = minValue;
        }
    }

    return output;
}

Mat guidedfilter(Mat& srcImage, Mat& srcClone, int r, double eps)
{
    Mat graymat;
    cvtColor(srcImage, graymat, CV_BGR2GRAY);
    graymat.convertTo(srcImage, CV_32FC1, 1 / 255.0);
    //srcClone.convertTo(srcClone, CV_64FC1);
    int nRows = srcImage.rows;
    int nCols = srcImage.cols;
    Mat boxResult;

    boxFilter(Mat::ones(nRows, nCols, srcImage.type()), boxResult, CV_32FC1, Size(r, r));

    Mat mean_I;
    boxFilter(srcImage, mean_I, CV_32FC1, Size(r, r));

    Mat mean_p;
    boxFilter(srcClone, mean_p, CV_32FC1, Size(r, r));

    Mat mean_Ip;
    boxFilter(srcImage.mul(srcClone), mean_Ip, CV_32FC1, Size(r, r));
    Mat cov_Ip = mean_Ip - mean_I.mul(mean_p);

    Mat mean_II;
    boxFilter(srcImage.mul(srcImage), mean_II, CV_32FC1, Size(r, r));

    Mat var_I = mean_II - mean_I.mul(mean_I);
    Mat var_Ip = mean_Ip - mean_I.mul(mean_p);

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

    Mat mean_a;
    boxFilter(a, mean_a, CV_32FC1, Size(r, r));
    mean_a = mean_a / boxResult;
    Mat mean_b;
    boxFilter(b, mean_b, CV_32FC1, Size(r, r));
    mean_b = mean_b / boxResult;

    Mat resultMat = mean_a.mul(srcImage) + mean_b;
    return resultMat;
}

Mat recover(Mat& srcimg, Mat& t, float *array)
{
    int nr = srcimg.rows, nl = srcimg.cols;
    float tnow = t.at<float>(0, 0);
    Mat finalimg = Mat::zeros(nr, nl, CV_32FC3);
    float val = 0;

    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++) 
        {
            tnow = *transPtr++;

            if (tnow < 0.1)
            {
                tnow = 0.1;
            }
            else if (tnow > 0.9)
            {
                tnow = 0.9;
            }

            for (int i = 0; i < 3; i++) 
            {
                val = (*srcPtr++ - array[i]) / tnow + array[i];
                *outPtr++ = val + 10 / 255.0;
            }

        }
    }
    return finalimg;
}

         

    

                                (a)                                                                            (b)                                                                      (c)                                                                    (d)

       上述一組實驗中，（a）圖為霧天圖像，(b)圖為本次實驗估計得到的景深圖像，(c)圖為估計得到的傳輸函數圖像，(d)圖為復原結果。從(b)圖中觀察可以看到，依據各個參數得到的景深圖像由於服從高斯噪聲分布，所以圖像中看一觀察到噪聲現象，所以有必要對該景深圖像進一步優化。作者先對景深圖作最小值濾波，然后再對其進行導向濾波，由此獲得優化的傳輸函數圖像(c)。

        如下給出幾組本次實驗的運行結果以供參考：       

   

   

 

Single image dehazing via reliability guided fusion

Irfan Riaz, Teng Yu, Yawar Rehman, Hyunchul Shin

       本文的去霧方法本質上是暗通道去霧方法的一種改進，效果很不錯，如果你追求去霧效果的話，本文的算法是一種很好的選擇。本文主要介紹的是對傳輸函數的優化，作者構造一種reliability map，將該圖像與暗通道得到的傳輸函數融合，從而得到優化的傳輸函數，以提升暗通道去霧方法的效果。

       文中將本文的暗通道優化方法以流程圖的形式給出，並和暗通道方法作了比較：

     對於圖像的固定窗口的尺寸，首先計算其暗通道圖像（c是RGB三個通道的某一個通道），可得：

對暗通道圖像作窗口為r x r (r=3)的均值濾波得到，然后對作窗口尺寸為的最大值濾波，可得block dark channel為：

 然后計算pixel dark channel為：

其中，與分別為對和作均值濾波，其窗口大小為r x r（r=3），為很小的數以防止除數為零。

       文中給出上述多圖對應的計算結果以及其復原效果：

        下面介紹計算reliability map以及其與暗通道圖像融合的方法。首先計算reliability map如下：

其中系數取0.0025，則將與暗通道圖像融合，並估計傳輸函數如下：

       文中給出上述各個參量的估計結果特例，以及復原結果：

          至此，優化的傳輸函數已經得到，但是將上述估計的傳輸函數代入到某些含有大量天空區域的霧天圖像進行圖像復原時，會發現天空區域的處理結果並不理想，因此有必要對天空區域作進一步處理。構造權重函數並修改如下：

將修改后的代替其原始的計算公式，即可解決圖像的天空區域復原的問題。式中參數與分別為10和0.2。

       下面給出本文的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);
Mat maxFilter(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(Mat& srcimg, Mat& transmission, int windowsize);
Mat recover(Mat& srcimg, Mat& t, float *array, int windowsize);

int main() 
{
    string loc = "E:\\code\\reliability\\project\\project\\cones.bmp";
    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;

    imshow("hazyimg", image);

    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);
    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);
    Mat normal(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);
    darkChannel.convertTo(normal, CV_32FC1, 1 / 255.0);

    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;

    //暗通道歸一化操作（除A）
    //(I / A)
    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;

    int radius = 3; int kernel = 2 * radius+1;
    Size win_size(kernel, kernel);

    Mat S(rows, cols, CV_32FC1);
    float w1 = 10.0; float w2 = 0.2; float min = 1.0;
    float b_A, g_A, r_A; float pixsum;

    for (int i = 0; i < nr; i++) 
    {
        const float* inData = convertImage.ptr<float>(i);
        float* outData = normalDark.ptr<float>(i);
        float* sData = S.ptr<float>(i);
        for (int j = 0; j < nl; j++) 
        {
            b = *inData++; g = *inData++; r = *inData++;

            b_A = b / tmp_A[0];
            g_A = g / tmp_A[1];
            r_A = r / tmp_A[2];

            min = b_A > g_A ? g_A : b_A;
            min = min > r_A ? r_A : min;
            *outData++ = min;

            pixsum = (b - tmp_A[0]) * (b - tmp_A[0]) + (g - tmp_A[1]) * (g - tmp_A[1]) + (r - tmp_A[2]) * (b - tmp_A[2]);
            *sData++ = exp((-1 * w1) * pixsum);
        }
    }

    imshow("S", S);

    //calculate the Iroi map
    Mat Ic = normalDark; Mat Icest;
    Mat Imin; Mat umin; Mat Ibeta;

    Ic = Ic.mul(Mat::ones(rows, cols, CV_32FC1) - w2 * S);
    imshow("Ic", Ic);
    Imin = minFilter(Ic, kernel);
    imshow("Imin", Imin);
    
    boxFilter(Imin, umin, CV_32F, win_size);
    Ibeta = maxFilter(umin, kernel);
    imshow("Ibeta", Ibeta);

    Mat ubeta; Mat uc;
    boxFilter(Ibeta, ubeta, CV_32F, win_size);
    boxFilter(Ic, uc, CV_32F, win_size);

    float fai = 0.0001; Mat Iroi;
    Mat weight = (Mat::ones(rows, cols, CV_32FC1)) * fai;
    divide((Ic.mul(ubeta)), (uc + weight), Iroi);
    imshow("Iroi", Iroi);

    //calculate the reliability map alpha
    Mat uepsilon; Mat alpha;
    Mat m = Ibeta - umin; Mat n = Ibeta - Iroi;
    boxFilter(m.mul(m) + n.mul(n), uepsilon, CV_32F, win_size);

    float zeta = 0.0025;
    uepsilon / (uepsilon + Mat::ones(rows, cols, CV_32FC1) * zeta);
    alpha = Mat::ones(rows, cols, CV_32FC1) - uepsilon / (uepsilon + Mat::ones(rows, cols, CV_32FC1) * zeta);
    imshow("alpha", alpha);

    //calculate the Idark map
    Mat Ialbe; Mat ualpha; Mat ualbe; Mat Idark;
    Ialbe = alpha.mul(Ibeta);

    boxFilter(alpha, ualpha, CV_32F, win_size);
    boxFilter(Ialbe, ualbe, CV_32F, win_size);

    Idark = Iroi.mul(Mat::ones(rows, cols, CV_32FC1) - ualpha) + ualbe;
    imshow("Idark", Idark);

    float w = 0.95;
    Mat t; t = Mat::ones(rows, cols, CV_32FC1) - w*Idark;
    int kernelSizeTrans = std::max(3, kernelSize);
    Mat trans = getTransmission(convertImage, t, kernelSizeTrans);
    imshow("t",trans);

    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("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;
}

Mat maxFilter(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)maxValue;
                continue;
            }
            *fOutData++ = maxValue;
        }
    }
    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(Mat& srcimg, Mat& transmission,  int windowsize)
{
    int nr = srcimg.rows, nl = srcimg.cols;
    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;

    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++) 
        {
            tnow = *transPtr++;
            tnow = std::max(tnow, t0);
            for (int i = 0; i < 3; i++)
            {
                val = (*srcPtr++ - array[i]) / tnow + array[i];
                *outPtr++ = val + 10 / 255.0;
            }
        }
    }
    return finalimg;
}

          以一組實驗結果為例，驗證實驗是否正確：

   

                                      hazy image                                                                                           Ibeta                                                                                         Iroi

  

                                       alpha                                                                                             transmission                                                                                     output

         下面給出本文方法的多組運行結果：