几种去雾算法介绍

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

         下面给出本文方法的多组运行结果：