1、何恺明的暗通道去雾算法
论文原文：Single Image Haze Removal Using Dark Channel Prior | IEEE Journals & Magazine | IEEE Xplore
参考博客：[论文阅读] (11)ACE算法和暗通道先验图像去雾算法（Rizzi | 何恺明老师）_暗通道去雾算法_Eastmount的博客-CSDN博客
参考公众号：GiantPandaCV
将各位大佬的实现代码粘出来，供大家学习！
#include <opencv2/opencv.hpp>
#include <iostream>
#include <Windows.h>
#include <vector>
using namespace std;
using namespace cv;
int getMax(Mat src)
{
int row = src.rows;
int col = src.cols;
int temp = 0;
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++)
{
temp = max((int)src.at<uchar>(i, j), temp);
}
if (temp == 255) return temp;
}
return temp;
}
Mat dehaze(Mat src)
{
double eps;
int row = src.rows;
int col = src.cols;
Mat M = Mat::zeros(row, col, CV_8UC1);
Mat M_max = Mat::zeros(row, col, CV_8UC1);
Mat M_ave = Mat::zeros(row, col, CV_8UC1);
Mat L = Mat::zeros(row, col, CV_8UC1);
Mat dst = Mat::zeros(row, col, CV_8UC3);
double m_av, A;
double sum = 0;
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++)
{
uchar r, g, b, temp1, temp2;
b = src.at<Vec3b>(i, j)[0];
g = src.at<Vec3b>(i, j)[1];
r = src.at<Vec3b>(i, j)[2];
temp1 = min(min(b, g), r);
temp2 = max(max(b, g), r);
M.at<uchar>(i, j) = temp1;
M_max.at<uchar>(i, j) = temp2;
sum += temp1;
}
}
m_av = sum / (row * col * 255);
eps = 0.75 / m_av;
boxFilter(M, M_ave, CV_8UC1, Size(15, 15)); // 中值滤波
double delta = min(0.9, eps * m_av);
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++)
{
L.at<uchar>(i, j) = min((int)(delta * M_ave.at<uchar>(i, j)), (int)M.at<uchar>(i, j));
}
}
A = (getMax(M_max) + getMax(M_ave)) * 0.5;
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++)
{
int temp = L.at<uchar>(i, j);
for (int k = 0; k < 3; k++)
{
int val = A * (src.at<Vec3b>(i, j)[k] - temp) / (A - temp);
if (val > 255) val = 255;
if (val < 0) val = 0;
dst.at<Vec3b>(i, j)[k] = val;
}
}
}
return dst;
}
// 获取最小值矩阵
int row, col;
int** getMinChannel(Mat src)
{
if (src.channels() != 3)
{
cout << "input error!" << endl;
exit(-1);
}
int** imgGray;
imgGray = new int* [row];
for (int i = 0; i < row; i++)
{
imgGray[i] = new int[col];
}
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++)
{
int localmin = 255;
for (int k = 0; k < 3; k++)
{
if (src.at<Vec3b>(i, j)[k] < localmin)
{
localmin = src.at<Vec3b>(i, j)[k];
}
}
imgGray[i][j] = localmin;
}
}
return imgGray;
}
// 求暗通道 通过最小矩阵，及poolsize，获取最小矩阵在poolsize范围内的最小值
int** getDarkChannel(int** imgGray, int blocksize = 3)
{
if (blocksize % 2 == 0 || blocksize < 3)
{
fprintf(stderr, "blocksize is not odd or too small!");
exit(-1);
}
int poolsize = (blocksize - 1) / 2;
int newheight = row + blocksize - 1;
int newwidth = col + blocksize - 1;
int** imgMiddle;
imgMiddle = new int* [newheight];
for (int i = 0; i < newheight; i++)
{
imgMiddle[i] = new int[newwidth];
}
for (int i = 0; i < newheight; i++)
{
for (int j = 0; j < newwidth; j++)
{
if (i < row && j < col)
{
imgMiddle[i][j] = imgGray[i][j];
}
else
{
imgMiddle[i][j] = 255;
}
}
}
int** imgDark;
imgDark = new int* [row];
for (int i = 0; i < row; i++)
{
imgDark[i] = new int[col];
}
int localmin = 255;
for (int i = poolsize; i < newheight - poolsize; i++)
{
for (int j = poolsize; j < newwidth - poolsize; j++)
{
localmin = 255;
for (int k = i - poolsize; k < i + poolsize + 1; k++)
{
for (int l = j - poolsize; l < j + poolsize + 1; l++)
{
if (imgMiddle[k][l] < localmin)
{
localmin = imgMiddle[k][l];
}
}
}
imgDark[i - poolsize][j - poolsize] = localmin;
}
}
return imgDark;
}
struct node {
int x, y, val;
node() {}
node(int _x,int _y,int _val):x(_x),y(_y),val(_val){}
bool operator<(const node& rhs) {
return val > rhs.val;
}
};
// 估算全局大气光值
int getGlobalAtmosphericLightValue(int** imgDark, Mat src, bool meanmodel = false, float percent = 0.001)
{
int size = row * col;
vector<node> nodes;
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++)
{
node temp;
temp.x = i, temp.y = j, temp.val = imgDark[i][j];
nodes.push_back(temp);
}
}
sort(nodes.begin(), nodes.end()); // 从小到大排序
int atmosphericLight = 0;
if (int(percent * size) == 0)
{
for (int i = 0; i < 3; i++)
{
if (src.at<Vec3b>(nodes[0].x, nodes[0].y)[i] > atmosphericLight) // imgDark最小处的原图中像素点中的最大值
{
atmosphericLight = src.at<Vec3b>(nodes[0].x, nodes[0].y)[i];
}
}
}
// 开启均值模式
//if (meanmodel == true)
//{
//int sum = 0;
//for (int i = 0; i < int(percent * size); i++)
//{
//for (int j = 0; j < 3; j++)
//{
//sum = sum + src.at<Vec3b>(nodes[i].x, nodes[i].y)[j];
//}
//}
//}
// 获取暗通道在前0.1%的位置的像素点在原图像中的最高亮度值
for (int i = 0; i < int(percent * size); i++)
{
for (int j = 0; j < 3; j++)
{
if (src.at<Vec3b>(nodes[i].x, nodes[i].y)[j] > atmosphericLight)
{
atmosphericLight = src.at<Vec3b>(nodes[i].x, nodes[i].y)[j];
}
}
}
return atmosphericLight;
}
// 恢复原图 omega 去雾比例 t0 最小透射率值
Mat getRecoverScene(Mat src,float omega=0.85,float t0 =0.5,int blocksize=15,bool meanmodel=false,float percent=0.001)
{
int** imgGray = getMinChannel(src); // 获取每个像素点b,g,r中的最小值
int** imgDark = getDarkChannel(imgGray, blocksize = blocksize); // 通过最小矩阵，及poolsize，获取imgGray在poolsize范围内的最小值
int atmosphericLight = getGlobalAtmosphericLightValue(imgDark, src, meanmodel = meanmodel, percent = percent); // 估算全局大气光值
float** imgDark2, ** transmission;
/*imgDark2 = new float* [row];
for (int i = 0; i < row; i++)
{
imgDark2[i] = new float[col];
}*/
transmission = new float* [row]; // t(x)
for (int i = 0; i < row; i++)
{
transmission[i] = new float[col];
}
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++)
{
//imgDark2[i][j] = float(imgDark[i][j]);
transmission[i][j] = 1 - omega * imgDark[i][j] / atmosphericLight;
if (transmission[i][j] < 0.1)
{
transmission[i][j] = 0.1;
}
}
}
Mat dst(src.rows, src.cols, CV_8UC3);
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++)
{
for (int channel = 0; channel < 3; channel++)
{
int temp = (src.at<Vec3b>(i, j)[channel] - atmosphericLight) / transmission[i][j] + atmosphericLight; // J(x) = (I(x) - A)/t(x) + A
if (temp > 255) temp = 255;
if (temp < 0) temp = 0;
dst.at<Vec3b>(i, j)[channel] = temp;
}
}
}
return dst;
}
Mat MedianFilterFogRemove(Mat src, float p = 0.85, int kernelsize = 15, int blocksize = 3, bool meanmodel = false, float percent = 0.001)
{
int** imgGray = getMinChannel(src); // 获取每个像素点b,g,r中的最小值
int** imgDark = getDarkChannel(imgGray, blocksize = blocksize); // 通过最小矩阵，及poolsize，获取imgGray在poolsize范围内的最小值
//int atmosphericLight = getGlobalAtmosphericLightValue(imgDark, src, meanmodel = meanmodel, percent = percent); // 估算全局大气光值
int Histgram[256] = {0};
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++)
{
Histgram[imgDark[i][j]]++;
}
}
int sum = 0, atmosphericLight = 0;
for (int i = 255; i >= 0; i--)
{
sum += Histgram[i]; // // sum的最大值为row*col
if (sum > row * col * 0.01)
{
atmosphericLight = i; // 从大到小的前0.1%所在的值作为大气光值
break;
}
}
int sumb = 0, sumg = 0, sumr = 0, amount = 0;
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++)
{
if (imgDark[i][j] >= atmosphericLight)
{
sumb += src.at<Vec3b>(i, j)[0];
sumg += src.at<Vec3b>(i, j)[1];
sumr += src.at<Vec3b>(i, j)[2];
amount++;
}
}
}
sumb /= amount;
sumg /= amount;
sumr /= amount; // 从所有最小值矩阵范围内的最小值矩阵中大于大气光值的，获取原图中的所有b,g,r值，并取均值
Mat filter(row, col, CV_8UC1);
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++)
{
filter.at<uchar>(i, j) = imgDark[i][j]; // 将imgDark赋给filter
}
}
Mat A(row, col, CV_8UC1);
medianBlur(filter, A, kernelsize); // 中值滤波，得到A
Mat temp(row, col, CV_8UC1);
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++)
{
int diff = filter.at<uchar>(i, j) - A.at<uchar>(i, j);
if (diff < 0) diff = -diff;
temp.at<uchar>(i, j) = diff;
}
}
medianBlur(temp, temp, kernelsize); // 中值滤波
Mat B(row, col, CV_8UC1);
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++)
{
int diff = A.at<uchar>(i, j) - temp.at<uchar>(i, j);
if (diff < 0) diff = 0;
B.at<uchar>(i, j) = diff;
}
}
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++)
{
int minn = B.at<uchar>(i, j) * p;
if (imgDark[i][j] > minn)
{
B.at<uchar>(i,j) = minn;
}
else
{
B.at<uchar>(i, j) = imgDark[i][j];
}
}
}
Mat dst(row, col, CV_8UC3);
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++)
{
int F = B.at<uchar>(i, j);
int value;
if (sumb != F)
{
value = sumb * (src.at<Vec3b>(i, j)[0] - F) / (sumb - F);
}
else
{
value = src.at<Vec3b>(i, j)[0];
}
if (value < 0) value = 0;
else if (value > 255) value = 255;
dst.at<Vec3b>(i, j)[0] = value;
if (sumg != F)
{
value = sumg * (src.at<Vec3b>(i, j)[1] - F) / (sumg - F);
}
else
{
value = src.at<Vec3b>(i, j)[1];
}
if (value < 0) value = 0;
else if (value > 255) value = 255;
dst.at<Vec3b>(i, j)[1] = value;
if (sumr != F)
{
value = sumr * (src.at<Vec3b>(i, j)[2] - F) / (sumr - F);
}
else
{
value = src.at<Vec3b>(i, j)[2];
}
if (value < 0) value = 0;
else if (value > 255) value = 255;
dst.at<Vec3b>(i, j)[2] = value;
}
}
return dst;
}
int main()
{
Mat image = imread("E:/dataset/python_dataset/VOC07+12+test/VOCdevkit/VOC2007/JPEGImages/2008_000084.jpg");
clock_t t0,t1,t2;
t0 = clock();
// 方式1 快速去雾
Mat dst = dehaze(image);
cout << "the method1 run time is :" << clock() - t0 << "ms" << endl;
// 方式2 暗通道去雾法
row = image.rows;
col = image.cols;
t1 = clock();
Mat dst1 = getRecoverScene(image);
cout << "the method2 run time is :" << clock() - t1 << "ms" << endl;
// 方式3 中值滤波进行去雾
t2 = clock();
Mat dst2 = MedianFilterFogRemove(image);
cout << "the method3 run time is :" << clock() - t2 << "ms" << endl;
imshow("origin", image);
imshow("result1", dst);
imshow("result2", dst1);
imshow("result3", dst2);
char key = waitKey();
if (key == 'q')
{
return 0;
}
system("pause");
return 0;
}
输入输出结果：
