当前位置：首页 > news >正文

C++基于opencv实现的暗通道的先验图像去雾

news 2025/10/16 6:11:44

测试环境：
vs2019
opencv==4.9.0，支持opencv3.4以上版本

前言

目前的图像去雾算法很多，思路基本上两条：

1,基于非物理模型，本质上是增强图像的对比度与颜色，并没有对雾天图像的形成原因进行分析而补偿。代表方法是直方图均衡化。效果一般。
2，基于物理模型，现在效果好的去雾算法都是基于物理模型，利用大气物理散射规律来建立图像还原模型，而不同的论文算法，用的模型都不尽相同。基于物理模型的算法因为基于模型更好的分析了含雾图像。并且与现实贴近，效果都不错，只是算法复杂度不同，计算时间长短不同而已。代表方法是何恺明博士的Single Image Haze Removal Using Dark Channel Prior，即基于暗通道先验的去雾算法。
目前感觉效果最好的就是基于暗通道先验的去雾算法。

##Single Image Haze Removal Using Dark Channel Prior
在去雾算法中，利用的以下物理模型：![此处输入图片的描述][2]
其中I(x)是有雾后的图像，J(x)是无雾图像，A是全球大气光照值，t(x)是透射率图。

上式经过化简之后可得到：![][3]
其实就是已知I(x)，然后通过分析I(x)，算出J(x)。

##暗通道
首先看什么是暗通道，在何的论文中，对于暗通道的定义是：在无雾图像中，在大多数局部区域内，其中的一些像素会在某个通道内含有非常低的像素值(换句话说也就是，在某个区域内，所有像素的各个通道的最小值的像素值非常小（0~16）)。这些像素值的产生主要是由于阴影（shadow）, 彩色物体（colorful object）（某一个通道的值太大，导致其他通道的值小）, 黑色物体等。
而在含雾图像中，因为雾气的存在，暗通道并不是全黑，比如下图：
![][4]
可以看出来，上方含雾的部分偏白，而下方城市部分大部分为黑色。这样我们都能得用暗通道来得到含雾图像中雾的有无，雾的浓度。这个是作者分析了很多图像后得出来的结论，所以说讲先验。

暗通道的计算方法也很容易
![][5]

其中Ω为一个正方形区域，即对每一个像素计算每一个通道一个区域内最小值的那个像素值作为暗通道的值。

##全球大气光照值
因为在计算透射率图的时候要用到全球大气光照值，所以先计算A。
A的计算方法很简单，

the pixels with highest intensity in the input image I is selected as the atmopheric light.

这是文中给出的方法，在暗通道中取亮度最高的前0.1%像素，然后找出在原始图像中I亮度最高的像素作为A。

##透射率
对于透射率的计算有以下式子：
![][6]

其中，w是为了不使图像失真，而引人的控制保留雾的比重的参数。
如果只是这么做了话，因为透射率图的计算方法，我们是在一个个区域内算的，所以透射率图会是一块块的。所以我们还要对透射率图进行处理。计算出精细的透射率图。
在原论文中，采用的软扣图算法，十分慢。

##去雾
我们现在已经计算得出了A，t(x)，I(x)，可以进行去雾了。

如下式：
![][7]

其中t0是为了防止t(x)中的值为0以及防止t(x)过下而使J(x)过大，J(x)过大会使图像偏白。

#改进

原算法中，计算精细的透射率图是采用的软扣图，而我改用了同是何博士发的算法论文中的导向滤波。而导向滤波是一个快速的滤波算法，感兴趣的可以看看论文。而现在MATLAB与OpenCV中集成了此滤波算法。在我写的时候还没有，所以自己写了一个导向滤波。其实何博士的论文对于导向滤波讲的很清楚，推荐大家去读一读。
算法对于天空部分计算不是很好，特别是还有雾的天空，会使天空出现明显的不自然。
这里我采用了一种很简单的方法来改进：
引入一个参数K，如果全球大气光照值A与对应的暗通道值的绝对值小于K就对除以这个绝对值并乘以K。如果大于的话就不变。

可以在最后计算式中加上一个小参数来增加图像的亮度。因为去雾后的图像会速度偏暗。何的论文中的图像都经过了再次处理。

A的值可以取3个通道的平均值，也可以取前0.1%的平均值。也可以两者都取。

在计算导向滤波的时间可以用原图，也可以用原图的灰度图。不过用原图的话，因为是三个通道都要计算，所以会比单通道要慢一点。

其实还有一些改进的地方但是都是一些可说可不说的，感觉大家都讲的差不多了。

对了，何博士好像不久前又写了一篇论文，是对导向滤波的改进,叫快速导向滤波。感叹一下，又会有人水一篇论文了：用快速导向滤波来做去雾。

##感想
滤波算法还有不少。但是基本上都出不了这个框架了。并且我也看不了不少论文，自己也写了一些，感觉就这个算法比较稳定，优化优化，对于任何图像都有不错的效果，并且速度快，能用于视频处理。

其实现在国内很多文章都是在优化这个算法，点子基本上都是放在了透射率图的优化上，计算速度，去雾效果等等，比如用其它的滤波方法，对t(x)的值进行限定，等等。各种方法都有，效果都还不错，只是速度应该快不起来，这就是学术与工业的差别。

我在参考另外一个博客的时候，发现他写了一个时间复杂度为O(1)的去雾算法，我也看了看，感觉效果没有这个好，虽然速度快。具体看：[一种可实时处理 O(1)复杂度图像去雾算法的实现][16]。
另外一个效果比较好的，那位大神也写了：[优化的对比度增强算法用于有雾图像的清晰化处理][17]。

代码

#define _CRT_SECURE_NO_WARNINGS
#include "main.h"
#include "guidedfilter.h"
using namespace std;
using namespace cv;int _PriorSize = 15;		//窗口大小
double _topbright = 0.001;//亮度最高的像素比例
double _w = 0.95;			//w
float t0 = 0.1;			//T(x)的最小值   因为不能让tx小于0 等于0 效果不好
int SizeH = 0;			//图片高度
int SizeW = 0;			//图片宽度
int SizeH_W = 0;			//图片中的像素总 数 H*W
Vec<float, 3> a;//全球大气的光照值
Mat trans_refine;
Mat dark_out1;char img_name[100] = "ny_night.jpg";//文件名//读入图片
Mat ReadImage()
{Mat img = imread(img_name);SizeH = img.rows;SizeW = img.cols;SizeH_W = img.rows*img.cols;Mat real_img(img.rows, img.cols, CV_32FC3);img.convertTo(real_img, CV_32FC3);real_img = real_img / 255;return real_img;//读入图片 并其转换为3通道的矩阵后 //除以255 将其RBG确定在0-1之间
}//计算暗通道
//J^{dark}(x)=min( min( J^c(y) ) )
Mat DarkChannelPrior(Mat img)
{Mat dark = Mat::zeros(img.rows, img.cols, CV_32FC1);//新建一个所有元素为0的单通道的矩阵for (int i = 0; i<img.rows; i++){for (int j = 0; j<img.cols; j++){dark.at<float>(i, j) = min(min(img.at<Vec<float, 3>>(i, j)[0], img.at<Vec<float, 3>>(i, j)[1]),min(img.at<Vec<float, 3>>(i, j)[0], img.at<Vec<float, 3>>(i, j)[2]));//就是两个最小值的过程}}erode(dark, dark_out1, Mat::ones(_PriorSize, _PriorSize, CV_32FC1));//这个函数叫腐蚀 做的是窗口大小的模板运算 ,对应的是最小值滤波,即 黑色图像中的一块块的东西return dark_out1;//这里dark_out1用的是全局变量，因为在其它地方也要用到
}
Mat DarkChannelPrior_(Mat img)//这个函数在计算tx用到，因为与计算暗通道一样都用到了求最小值的过程，变化不多，所以改了改就用这里了
{double A = (a[0] + a[1] + a[2]) / 3.0;//全球大气光照值 此处是3通道的平均值Mat dark = Mat::zeros(img.rows, img.cols, CV_32FC1);Mat dark_out = Mat::zeros(img.rows, img.cols, CV_32FC1);for (int i = 0; i<img.rows; i++){for (int j = 0; j<img.cols; j++){dark.at<float>(i, j) = min(min(img.at<Vec<float, 3>>(i, j)[0] / A, img.at<Vec<float, 3>>(i, j)[1] / A),min(img.at<Vec<float, 3>>(i, j)[0] / A, img.at<Vec<float, 3>>(i, j)[2] / A));//同理}}erode(dark, dark_out, Mat::ones(_PriorSize, _PriorSize, CV_32FC1));//同上return dark_out;}//计算A的值
Vec<float, 3> Airlight(Mat img, Mat dark)//vec<float ,3>表示有3个大小的vector 类型为float
{int n_bright = _topbright*SizeH_W;Mat dark_1 = dark.reshape(1, SizeH_W);//这里dark_1是一个有图片像素那么多行的矩阵 方便下面循环计算vector<int> max_idx;float max_num = 0;Vec<float, 3> A(0, 0, 0);Mat RGBPixcels = Mat::ones(n_bright, 1, CV_32FC3);for (int i = 0; i<n_bright; i++){max_num = 0;max_idx.push_back(max_num);for (float * p = (float *)dark_1.datastart; p != (float *)dark_1.dataend; p++){if (*p>max_num){max_num = *p;//记录光照的最大值max_idx[i] = (p - (float *)dark_1.datastart);//位置RGBPixcels.at<Vec<float, 3>>(i, 0) = ((Vec<float, 3> *)img.data)[max_idx[i]];//对应 的三个通道的值给RGBPixcels}}((float *)dark_1.data)[max_idx[i]] = 0;//访问过的标记为0，这样就不会重复访问}for (int j = 0; j<n_bright; j++){A[0] += RGBPixcels.at<Vec<float, 3>>(j, 0)[0];A[1] += RGBPixcels.at<Vec<float, 3>>(j, 0)[1];A[2] += RGBPixcels.at<Vec<float, 3>>(j, 0)[2];}//将光照值累加A[0] /= n_bright;A[1] /= n_bright;A[2] /= n_bright;//除以总数   即取所有符合的点的平均值。return A;
}//Calculate Transmission Matrix
Mat TransmissionMat(Mat dark)
{double A = (a[0] + a[1] + a[2]) / 3.0;for (int i = 0; i < dark.rows; i++){for (int j = 0; j < dark.cols; j++){double temp = (dark_out1.at<float>(i, j));double B = fabs(A - temp);//	conut++;//cout << conut << endl;//if (B==)if (B - 0.3137254901960784 < 0.0000000000001)//K=80    80/255=0.31   这里浮点数要这样做减法才能正确的比较{dark.at<float>(i, j) = (1 - _w*dark.at<float>(i, j))*(0.3137254901960784 / (B));//此处为改过的式子部分}else{dark.at<float>(i, j) = 1 - _w*dark.at<float>(i, j);}if (dark.at<float>(i, j) <= 0.2)//保证Tx不失真，因为会以上除出的结果会有不对{dark.at<float>(i, j) = 0.5;}if (dark.at<float>(i, j) >= 1)//同上{dark.at<float>(i, j) = 1.0;}}}return dark;
}
Mat TransmissionMat1(Mat dark)
{double A = (a[0] + a[1] + a[2]) / 3.0;for (int i = 0; i < dark.rows; i++){for (int j = 0; j < dark.cols; j++){dark.at<float>(i, j) = (1 - _w*dark.at<float>(i, j));}}return dark;
}
//Calculate Haze Free Image
Mat hazefree(Mat img, Mat t, Vec<float, 3> a, float exposure = 0)//此处的exposure的值表示去雾后应该加亮的值。
{double AAA = a[0];if (a[1] > AAA)AAA = a[1];if (a[2] > AAA)AAA = a[2];//取a中的最大的值//新开一个矩阵Mat freeimg = Mat::zeros(SizeH, SizeW, CV_32FC3);img.copyTo(freeimg);//两个迭代器，这样的写法可以不用两层循环，比较快点Vec<float, 3> * p = (Vec<float, 3> *)freeimg.datastart;float * q = (float *)t.datastart;for (; p<(Vec<float, 3> *)freeimg.dataend && q<(float *)t.dataend; p++, q++){(*p)[0] = ((*p)[0] - AAA) / std::max(*q, t0) + AAA + exposure;(*p)[1] = ((*p)[1] - AAA) / std::max(*q, t0) + AAA + exposure;(*p)[2] = ((*p)[2] - AAA) / std::max(*q, t0) + AAA + exposure;}return freeimg;
}void printMatInfo(char * name, Mat m)
{cout << name << ":" << endl;cout << "\t" << "cols=" << m.cols << endl;cout << "\t" << "rows=" << m.rows << endl;cout << "\t" << "channels=" << m.channels() << endl;
}//Main Function
int main(int argc, char * argv[])
{Mat dark_channel;Mat trans;Mat img;Mat free_img;char filename[100];while (_access(img_name, 0) != 0)//检测图片是否存在{std::cout << "The image " << img_name << " don't exist." << endl << "Please enter another one:" << endl;cin >> filename;}clock_t start, finish;double duration1, duration3, duration4, duration7;//读入图片cout << "读入图片 ..." << endl;start = clock();img = ReadImage();cv::namedWindow("原图", cv::WINDOW_KEEPRATIO);imshow("原图", img);//printMatInfo("img", img);finish = clock();duration1 = (double)(finish - start) / CLOCKS_PER_SEC;cout << "Time Cost: " << duration1 << "s" << endl;//输出这一步的时间cout << endl;//计算暗通道cout << "计算暗通道 ..." << endl;start = clock();dark_channel = DarkChannelPrior(img);cv::namedWindow("Dark Channel Prior", cv::WINDOW_KEEPRATIO);imshow("Dark Channel Prior", dark_channel);printMatInfo("dark_channel", dark_channel);finish = clock();duration3 = (double)(finish - start) / CLOCKS_PER_SEC;cout << "Time Cost: " << duration3 << "s" << endl;cout << endl;//计算全球光照值cout << "计算A值 ..." << endl;start = clock();a = Airlight(img, dark_channel);cout << "Airlight:\t" << " B:" << a[0] << " G:" << a[1] << " R:" << a[2] << endl;finish = clock();duration4 = (double)(finish - start) / CLOCKS_PER_SEC;cout << "Time Cost: " << duration4 << "s" << endl;cout << endl;//计算txcout << "Reading Refine Transmission..." << endl;trans_refine = TransmissionMat(DarkChannelPrior_(ReadImage()));//printMatInfo("trans_refine", trans_refine);//imshow("Refined Transmission Mat",trans_refine);cout << endl;Mat tran = guidedFilter(img, trans_refine, 60, 0.0001);//导向滤波 得到精细的透射率图//imshow("fitler", tran);//去雾cout << "Calculating Haze Free Image ..." << endl;start = clock();free_img = hazefree(img, tran, a, 0);//此处 如果用tran的话就是导向滤波部分//如果是trans_refine 就没有用导向滤波 效果不是那么						的好/*上面第四个参数是用来增加亮度的，0.1比较好cv::namedWindow("去雾后",cv::WINDOW_KEEPRATIO);								 */imshow("去雾后", free_img);finish = clock();duration7 = (double)(finish - start) / CLOCKS_PER_SEC;cout << "Time Cost: " << duration7 << "s" << endl;cout << "Total Time Cost: " << duration1 + duration3 + duration4 + duration7 << "s" << endl;//保存图片的代码imwrite("output.jpg", free_img * 255);waitKey();cout << endl;return 0;
}