拜读 何恺明 博士的论文，醍醐灌顶。
    Single_Image_Haze_Removal_Using_Dark_Channel_Prior.pdf
    暗通道先验滤波方法
    水下去雾.pdf
    实现代码

    1. #include <opencv2/opencv.hpp>
    2. #include<opencv2/imgproc/imgproc.hpp>
    3. #include<opencv2/highgui/highgui.hpp>
    4. #include<opencv2/core/core.hpp>
    5. #include<iostream>
    7. void makeDepth32f(Mat& source, Mat& output)
    8. {
    9.     if ((source.depth() != CV_32F) > FLT_EPSILON)
    10.     source.convertTo(output, CV_32F);
    11.     else
    12.         output = source;
    13. }
    16. void guidedFilter(Mat& source, Mat& guided_image, Mat& output, int radius, float epsilon)
    17. {
    18.     //CV_Assert(radius >= 2 && epsilon > 0);
    19.     CV_Assert(source.data != NULL && source.channels() == 1);//可改变输入图像类型（通道）
    20.     CV_Assert(guided_image.channels() == 1);                 //导向图一般为单通道
    21.     CV_Assert(source.rows == guided_image.rows && source.cols == guided_image.cols);
    23.     Mat guided;
    24.     if (guided_image.data == source.data) 
    25.         guided_image.copyTo(guided);
    26.     else
    27.         guided = guided_image;
    29.     Mat source_32f, guided_32f;
    30.     makeDepth32f(source, source_32f);//将输入扩展为32位浮点型，以便以后做乘法
    31.     makeDepth32f(guided, guided_32f);
    33.     Mat mat_Ip, mat_I2;   //计算I*p和I*I
    34.     multiply(guided_32f, source_32f, mat_Ip);
    35.     multiply(guided_32f, guided_32f, mat_I2);
    37.     Mat mean_p, mean_I, mean_Ip, mean_I2;   //计算各种均值
    38.     Size win_size(2 * radius + 1, 2 * radius + 1);
    39.     boxFilter(source_32f, mean_p, CV_32F, win_size);
    40.     boxFilter(guided_32f, mean_I, CV_32F, win_size);
    41.     boxFilter(mat_Ip, mean_Ip, CV_32F, win_size);
    42.     boxFilter(mat_I2, mean_I2, CV_32F, win_size);
    44.     Mat cov_Ip = mean_Ip - mean_I.mul(mean_p);//计算Ip的协方差和I的方差
    45.     Mat var_I = mean_I2 - mean_I.mul(mean_I);
    46.     var_I += epsilon;
    48.     Mat a, b;   //求a和b
    49.     divide(cov_Ip, var_I, a);
    50.     b = mean_p - a.mul(mean_I);
    52.     Mat mean_a, mean_b;  //对包含像素i的所有a、b做平均
    53.     boxFilter(a, mean_a, CV_32F, win_size);
    54.     boxFilter(b, mean_b, CV_32F, win_size);
    56.     output = mean_a.mul(guided_32f) + mean_b;//计算输出 (depth == CV_32F)
    57. }
    60. void HazeRemoval(Mat& image, Mat& imageRGB)
    61. {
    62.     CV_Assert(!image.empty() && image.channels() == 3);
    63.     Mat fImage;
    64.     image.convertTo(fImage, CV_32FC3, 1.0 / 255, 0);//图片归一化
    66.     Mat fImageBorder;
    67.     int hPatch = 15,vPatch = 15;//设定最小滤波patch的大小,且均为奇数
    68.     copyMakeBorder(fImage, fImageBorder, vPatch / 2, vPatch / 2, hPatch / 2, hPatch / 2, BORDER_REPLICATE);//给归一化的图片添加边界
    69.     vector<Mat> fImageBorderVector(3);
    70.     split(fImageBorder, fImageBorderVector);//分离通道
    72.     Mat darkChannel(image.rows, image.cols, CV_32FC1);//创建darkChannel
    73.     double minTemp, minPixel;
    74.     for (unsigned int r = 0; r < darkChannel.rows; r++)
    75.     {
    76.         for (unsigned int c = 0; c < darkChannel.cols; c++)
    77.         {
    78.             minPixel = 1.0;
    79.             for (vector<Mat>::iterator it = fImageBorderVector.begin(); it != fImageBorderVector.end(); it++)
    80.             {
    81.                 Mat roi(*it, Rect(c, r, hPatch, vPatch));
    82.                 minMaxLoc(roi, &minTemp);
    83.                 minPixel = min(minPixel, minTemp);
    84.             }
    85.             darkChannel.at<float>(r, c) = float(minPixel);
    86.         }
    87.     }
    88.     //darkChannel.convertTo(darkChannel8U, CV_8UC1, 255, 0);
    90.     //求出A(global atmospheric light),计算出darkChannel中前top个亮的值,论文中取值为0.1%
    91.     float top = 0.001;
    92.     float numberTop = top * darkChannel.rows*darkChannel.cols;
    93.     Mat darkChannelVector;
    94.     darkChannelVector = darkChannel.reshape(1, 1);//reshape的一个参数表示通道数，第二个表示矩阵行数
    95.     Mat_<int> darkChannelVectorIndex;
    96.     sortIdx(darkChannelVector, darkChannelVectorIndex, CV_SORT_EVERY_ROW + CV_SORT_DESCENDING);//降序，返回像素索引
    98.     int count = 0, temp = 0;
    99.     unsigned int x, y; //映射回暗通道图的像素位置
    100.     Mat mask(darkChannel.rows, darkChannel.cols, CV_8UC1);//制作掩码,注意mask的类型必须是CV_8UC1
    101.     for (unsigned int r = 0; r < darkChannelVectorIndex.rows; r++)
    102.     {
    103.         for (unsigned int c = 0; c < darkChannelVectorIndex.cols; c++)
    104.         {
    105.             temp = darkChannelVectorIndex.at<int>(r, c);
    106.             x = temp / darkChannel.cols;
    107.             y = temp % darkChannel.cols;
    109.             if (count < numberTop) {
    110.                 mask.at<uchar>(x, y) = 1;
    111.                 count++;
    112.             }
    113.             else
    114.                 mask.at<uchar>(x, y) = 0;
    116.         }
    117.     }
    119.     vector<double> A(3);                //分别存取B,G,R通道的最大A值
    120.     vector<Mat> fImageBorderVectorA(3);//在求第三步的t(x)时，会用到以下的矩阵，这里可以提前求出
    121.     vector<double>::iterator itA = A.begin();
    122.     vector<Mat>::iterator it = fImageBorderVector.begin();
    123.     vector<Mat>::iterator itAA = fImageBorderVectorA.begin();
    124.     for (; it != fImageBorderVector.end() && itA != A.end() && itAA != fImageBorderVectorA.end(); it++, itA++, itAA++)
    125.     {
    126.         Mat roi(*it, Rect(0, 0, darkChannel.cols, darkChannel.rows));
    127.         minMaxLoc(roi, 0, &(*itA), 0, 0, mask);
    128.         (*itAA) = (*it) / (*itA); //注意：这个地方有除号，但是没有判断是否等于0,*itA等于零的可能性很小
    129.     }
    131.     //求出t(x)
    132.     Mat darkChannelA(darkChannel.rows, darkChannel.cols, CV_32FC1);
    133.     float omega = 0.95;      //论文中取值为0.95
    134.     for (unsigned int r = 0; r < darkChannel.rows; r++)
    135.     {
    136.         for (unsigned int c = 0; c < darkChannel.cols; c++)
    137.         {
    138.             minPixel = 1.0;
    139.             for (itAA = fImageBorderVectorA.begin(); itAA != fImageBorderVectorA.end(); itAA++)
    140.             {
    141.                 Mat roi(*itAA, Rect(c, r, hPatch, vPatch));
    142.                 minMaxLoc(roi, &minTemp);
    143.                 minPixel = min(minPixel, minTemp);
    144.             }
    145.             darkChannelA.at<float>(r, c) = float(minPixel);
    146.         }
    147.     }
    148.     Mat tx1 = 1.0 - omega * darkChannelA;
    149.     Mat tx(darkChannel.rows, darkChannel.cols, CV_32FC1);
    151.     guidedFilter(tx1, tx1,tx, 8, 500);
    153.     namedWindow("tx", CV_WINDOW_AUTOSIZE);
    154.     imshow("tx", tx);
    157.     //求出J(x)
    158.     float t0 = 0.1;//论文中取0.1
    159.     //Mat jx(image.rows, image.cols, CV_32FC3);
    160.     for (size_t r = 0; r < imageRGB.rows; r++)
    161.     {
    162.         for (size_t c = 0; c < imageRGB.cols; c++)
    163.         {
    164.             imageRGB.at<Vec3f>(r, c) = Vec3f((fImage.at<Vec3f>(r, c)[0] - A[0]) / max(tx.at<float>(r, c), t0) + A[0], (fImage.at<Vec3f>(r, c)[1] - A[1]) / max(tx.at<float>(r, c), t0) + A[1], (fImage.at<Vec3f>(r, c)[2] - A[2]) / max(tx.at<float>(r, c), t0) + A[2]);
    165.         }
    166.     }
    168. }
    171. int main()
    172. {
    173.     Mat image = imread("C:\\AllProgram\\testimage\\haze1.bmp");//设置0写入灰度图像
    174.     if (image.empty())
    175.     {
    176.         std::cout << "打开图片失败,请检查" << std::endl;
    177.         system("pause");
    178.         return -1;
    179.     }
    180.     Mat testimage(image.size(), CV_32FC3);
    181.     HazeRemoval(image,testimage);
    183.     namedWindow("原图", CV_WINDOW_AUTOSIZE);
    184.     namedWindow("变换后的图", CV_WINDOW_AUTOSIZE);
    185.     imshow("原图", image);
    186.     imshow("变换后的图", testimage);
    187.         waitKey(0);
    188.     destroyAllWindows();
    190.     system("pause");
    191.     return 0;
    192. }