几种去雾算法介绍
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; }
对于上述代码，由于输入雾天图像不同，调整对应雾天图像的无穷远处大气光值即可。如下两图分别为雾天图像与去雾结果图像：
Single image dehazing via reliability guided fusion
Irfan Riaz, Teng Yu, Yawar Rehman, Hyunchul Shin
Hanyang University
本文的去雾方法本质上是暗通道去雾方法的一种改进，效果很不错，如果你追求去雾效果的话，本文的算法是一种很好的选择。本文主要介绍的是对传输函数的优化，作者构造一种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
下面给出本文方法的多组运行结果：
Efficient Image Dehazing with Boundary Constraint and Contextual Regularization
Gaofeng Meng, Ying Wang, Jiangyong Duan, Shiming Xiang, Chunhong Pan
National Laboratory of Pattern Recognition, Institute of Automation
Chinese Academy of Science, Beijing, P.R. China
本文依然依据大气散射模型进行图像去雾，则大气散射模型在此处不再叙述，直接介绍去雾算法。雾天图像中某一点像素处的雾气比例由传输函数决定，根据大气散射模型变形可得传输函数的倒数的表达式为：
考虑到图像中像素的灰度值总是有一定的界限，设其上界和下界分别为
和
，则无雾场景用边界条件界定为：
图像中的各个元素分量及其示意图如下图所示：
对于每一个
，均对应一个
，由上图中的几何关系可得：
其中
为
的下界，由下式给出：
再对作最大值滤波和最小值滤波，可得到传输函数的粗略估计值：
下面通过加权LI范数正则化对初始估计的传输函数进行优化：通常，图像中矩形窗口中邻近像素具有相似的景深，可用下式描述该规律：
其中
为权重函数，有很多种权重函数可以选择，本文选取关于雾天图像像素灰度值的高斯分布函数：
其中
为预设系数。则关于传输函数
的正则化表达式为：
其中，
表示整幅图像区域，
为像素
的领域。上式表示连续积分的过程，由于图像为离散化的数据，则上式的离散化形式为：
其中，
为图像像素的索引，
为的离散化形式。
若引入微分算子，则正则化的离散形式公式为：
该式进一步精简为：
其中，
表示矩形领域窗口中所有像素的集合，运算符
表示像素灰度值的点乘，
为卷积算符，
为一阶微分算子，
表示权重矩阵。
将前面所述的高斯分布函数的权重函数用微分算子的形式重新表述为（即
）：
文中选择了多个微分算子，如下：
作者在估计大气光向量
时，选取了何恺明的方法：首先计算雾天图像的暗通道图像，然后取对应于暗通道图像的最亮的0.1%像素处的原雾天图像像素，计算这些像素处三个通道灰度值的平均值，即可得到大气光向量
。
接下来着重介绍传输函数的优化过程。首先定义关于传输函数
的目标函数：
其中第一项为数据项，衡量与边界条件得到的传输函数
之间的相似程度，
为平衡数据项与正则项的系数。
为了优化该目标函数，引入辅助项
以构造一个新的代价函数（cost function）：
其中，
为权重系数。
对于固定系数
，最小化损失函数的过程即为分别关于
和
求解最优化的过程，求解过程为：首先固定
求解的最优解，然后固定求解的最优解。
具体步骤如下：1，优化
：若固定，则正则化函数变为：
该优化过程可归结为如下形式的最优化问题：
其中，系数
，
和
已经给出，则该最优化问题可直接通过下式解出：
2，优化：若固定损失函数中的
分量，则损失函数可进一步简化为：
对上式关于求导并变形，且令求解结果为0，可得：
其中
通过算子
关于中心像素作镜像得到，对上式两边作2D傅里叶变换和傅里叶逆变换，可解得的最优解
为：
其中
表示傅里叶变换，
表示傅里叶逆变换，
表示复数共轭。
需要不断迭代参数
以求解上述最优化过程，将
从
增加到
，且增加的步长为
。
在实验过程中，参数选择分别为：
，
和
。
使用C++实现本文的代码如下：
#include <iostream> #include <stdlib.h> #include <time.h> #include <cmath> #include <algorithm> #include <opencv2/opencv.hpp> #include <opencv2/highgui/highgui.hpp> #include <opencv2/imgproc/imgproc.hpp> using namespace cv; using namespace std; Mat minFilter(Mat& input, int kernelSize); Mat boundCon(Mat& srcimg, float *airlight, int C0, int C1); Mat calTrans(Mat& srcimg, Mat& t, float lambda, float param); Mat calWeight(Mat& srcimg, Mat& kernel, float param); Mat fft2(Mat I, Size size); void ifft2(const Mat &src, Mat &Fourier); void psf2otf(Mat& src, int rows, int cols, Mat& dst); void circshift(Mat& img, int dw, int dh, Mat& dst); Mat sign(Mat& input); Mat fliplr(Mat& input); Mat flipud(Mat& input); Mat recover(Mat& srcimg, Mat& t, float *airlight, float delta); int main() { clock_t start, finish; double duration; cout << "A defog program" << endl << "----------------" << endl; Mat image = imread("forest.png"); imshow("hazyiamge", image); cout << "input hazy image" << endl; Mat resizedImage; int originRows = image.rows; int originCols = image.cols; double scale = 1.0; 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); vector<Mat> channels(3); split(convertImage, channels); Mat R = channels[2]; Mat G = channels[1]; Mat B = channels[0]; int kernelSize = 15; Mat minImgR = minFilter(channels[2], kernelSize); Mat minImgG = minFilter(channels[1], kernelSize); Mat minImgB = minFilter(channels[0], kernelSize); double RminValue = 0, RmaxValue = 0; Point RminPt, RmaxPt; minMaxLoc(minImgR, &RminValue, &RmaxValue, &RminPt, &RmaxPt); cout << RmaxValue << endl; double GminValue = 0, GmaxValue = 0; Point GminPt, GmaxPt; minMaxLoc(minImgG, &GminValue, &GmaxValue, &GminPt, &GmaxPt); cout << GmaxValue << endl; double BminValue = 0, BmaxValue = 0; Point BminPt, BmaxPt; minMaxLoc(minImgB, &BminValue, &BmaxValue, &BminPt, &BmaxPt); cout << BmaxValue << endl; float A[3]; A[0] = BmaxValue; A[1] = GmaxValue; A[2] = RmaxValue; Mat trans = boundCon(convertImage, A, 30, 300); imshow("t", trans); int s = 3; Mat element; Mat transrefine; element = getStructuringElement(MORPH_ELLIPSE, Size(s, s)); morphologyEx(trans, transrefine, CV_MOP_CLOSE, element); imshow("transrefine", transrefine); float lambda = 2.0; float param = 0.5; Mat finaltrans = calTrans(convertImage, transrefine, lambda, param); imshow("finaltrans", finaltrans); float delta = 0.85; Mat recoverimg = recover(convertImage, finaltrans, A, delta); imshow("recover", recoverimg); recoverimg.convertTo(recoverimg, CV_8UC3, 255.0, 0); imwrite("recover.png", recoverimg); finish = clock(); duration = (double)(finish - start); cout << "defog used " << duration << "ms time;" << endl; waitKey(0); return 0; } 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 boundCon(Mat& srcimg, float *airlight, int C0, int C1) { int nr = srcimg.rows; int nl = srcimg.cols; float c0 = C0 / 255.0; float c1 = C1 / 255.0; Mat trans = Mat::zeros(nr, nl, CV_32FC1); float airlight0 = airlight[0]; float airlight1 = airlight[1]; float airlight2 = airlight[2]; for (unsigned int r = 0; r < nr; r++) { float* srcPtr = srcimg.ptr<float>(r); float* tPtr = trans.ptr<float>(r); float B, G, R; float t_b, t_g, t_r; float tb; float tmax = 1.0; for (unsigned int c = 0; c < nl; c++) { B = *srcPtr++; G = *srcPtr++; R = *srcPtr++; t_b = std::max((airlight0 - B) / (airlight0 - c0), (B - airlight0) / (c1 - airlight0)); t_g = std::max((airlight1 - G) / (airlight1 - c0), (G - airlight1) / (c1 - airlight1)); t_r = std::max((airlight2 - R) / (airlight2 - c0), (R - airlight2) / (c1 - airlight2)); tb = t_b > t_g ? t_b : t_g; tb = tb > t_r ? tb : t_r; tb = std::min(tb, tmax); *tPtr++ = tb; } } return trans; } Mat calTrans(Mat& srcimg, Mat& t, float lambda, float param) { int nsz = 3; int NUM = nsz * nsz; int nr = t.rows; int nl = t.cols; Size size = t.size(); Mat kernel1 = (Mat_<float>(3, 3) << 5, 5, 5, -3, 0, -3, -3, -3, -3); Mat kernel2 = (Mat_<float>(3, 3) << -3, 5, 5, -3, 0, 5, -3, -3, -3); Mat kernel3 = (Mat_<float>(3, 3) << -3, -3, 5, -3, 0, 5, -3, -3, 5); Mat kernel4 = (Mat_<float>(3, 3) << -3, -3, -3, -3, 0, 5, -3, 5, 5); Mat kernel5 = (Mat_<float>(3, 3) << 5, 5, 5, -3, 0, -3, -3, -3, -3); Mat kernel6 = (Mat_<float>(3, 3) << -3, -3, -3, 5, 0, -3, 5, 5, -3); Mat kernel7 = (Mat_<float>(3, 3) << 5, -3, -3, 5, 0, -3, 5, -3, -3); Mat kernel8 = (Mat_<float>(3, 3) << 5, 5, -3, 5, 0, -3, -3, -3, -3); normalize(kernel1, kernel1, 1.0, 0.0, NORM_L2); normalize(kernel2, kernel2, 1.0, 0.0, NORM_L2); normalize(kernel3, kernel3, 1.0, 0.0, NORM_L2); normalize(kernel4, kernel4, 1.0, 0.0, NORM_L2); normalize(kernel5, kernel5, 1.0, 0.0, NORM_L2); normalize(kernel6, kernel6, 1.0, 0.0, NORM_L2); normalize(kernel7, kernel7, 1.0, 0.0, NORM_L2); normalize(kernel8, kernel8, 1.0, 0.0, NORM_L2); Mat d = Mat::zeros(3, 3, CV_32FC(8)); vector<Mat> dchannels(8); split(d, dchannels); dchannels[0] = kernel1; dchannels[1] = kernel2; dchannels[2] = kernel3; dchannels[3] = kernel4; dchannels[4] = kernel5; dchannels[5] = kernel6; dchannels[6] = kernel7; dchannels[7] = kernel8; Mat wfun = Mat::zeros(nr, nl, CV_32FC(8)); vector<Mat> wfunchannels(8); for (int k = 0; k < 8; k++) { wfunchannels[k] = calWeight(srcimg, dchannels[k], param); } Mat Tf; Tf=fft2(t,size); Mat D = Mat::zeros(nr, nl, CV_32FC(8)); Mat DS = Mat::zeros(nr, nl, CV_32FC1); vector<Mat> Dchannels(8); split(D, Dchannels); for (int k = 0; k < 8; k++) { psf2otf(dchannels[k], nr, nl, Dchannels[k]); Mat Dchannels_temp = Mat::zeros(nr, nl, CV_32FC1); vector<Mat> Dchannels_temp1(2); split(Dchannels[k], Dchannels_temp1); magnitude(Dchannels_temp1[0], Dchannels_temp1[1], Dchannels_temp); pow(Dchannels_temp, 2, Dchannels_temp); DS += Dchannels_temp; } float beta = 1.0; float beta_rate = 2 * sqrt(2); float beta_max = 256; while (beta < beta_max) { float gamma = lambda / beta; Mat dt = Mat::zeros(nr, nl, CV_32FC(8)); Mat u = Mat::zeros(nr, nl, CV_32FC(8)); Mat DU = Mat::zeros(nr, nl, CV_32FC2); vector<Mat> dtchannels(8); split(dt, dtchannels); vector<Mat> uchannels(8); split(u, uchannels); for (int k = 0; k < 8; k++) { float zero = 0.0; filter2D(t, dtchannels[k], t.depth(), dchannels[k]); uchannels[k] = max(abs(dtchannels[k]) - (1 / (8 * beta)) * wfunchannels[k], zero).mul(sign(dtchannels[k])); Mat dfliplr = fliplr(dchannels[k]); Mat dflipud = flipud(dfliplr); filter2D(uchannels[k], uchannels[k], uchannels[k].depth(), dflipud); Mat DUtemp; DUtemp=fft2(uchannels[k],size); DU += DUtemp; } Mat ifft; Mat tfftup = gamma * Tf + DU; Mat tfftdown = gamma * Mat::ones(nr, nl, CV_32FC1) + DS; vector<Mat> tfft_temp(2); split(tfftup, tfft_temp); tfft_temp[0] = tfft_temp[0] / tfftdown; tfft_temp[1] = tfft_temp[1] / tfftdown; Mat final_tfft; merge(tfft_temp, final_tfft); ifft2(final_tfft, ifft); t = abs(ifft); beta = beta * beta_rate; } Mat trans(nr, nl, CV_32FC1); trans = t; return trans; } Mat calWeight(Mat& srcimg, Mat& kernel, float param) { int nr = srcimg.rows; int nl = srcimg.cols; Mat wfun; float weight = -1 / (param * 2); vector<Mat> channels(3); split(srcimg, channels); Mat d_b; Mat d_g; Mat d_r; filter2D(channels[0], d_b, channels[0].depth(), kernel); filter2D(channels[1], d_g, channels[1].depth(), kernel); filter2D(channels[2], d_r, channels[2].depth(), kernel); exp(weight * (d_b.mul(d_b) + d_g.mul(d_g) + d_r.mul(d_r)), wfun); return wfun; } Mat fft2(Mat I, Size size) { Mat If = Mat::zeros(I.size(), I.type()); Size dftSize; // compute the size of DFT transform dftSize.width = getOptimalDFTSize(size.width); dftSize.height = getOptimalDFTSize(size.height); // allocate temporary buffers and initialize them with 0's Mat tempI(dftSize, I.type(), Scalar::all(0)); //copy I to the top-left corners of temp Mat roiI(tempI, Rect(0, 0, I.cols, I.rows)); I.copyTo(roiI); if (I.channels() == 1) { dft(tempI, If, DFT_COMPLEX_OUTPUT); } else { vector<Mat> channels; split(tempI, channels); for (int n = 0; n<I.channels(); n++) { dft(channels[n], channels[n], DFT_COMPLEX_OUTPUT); } cv::merge(channels, If); } return If(Range(0, size.height), Range(0, size.width)); } void ifft2(const Mat &src, Mat &Fourier) { int mat_type = src.type(); assert(mat_type < 15); if (mat_type < 7) { Mat planes[] = { Mat_<double>(src), Mat::zeros(src.size(), CV_64F) }; merge(planes, 2, Fourier); dft(Fourier, Fourier, DFT_INVERSE + DFT_SCALE, 0); } else { Mat tmp; dft(src, tmp, DFT_INVERSE + DFT_SCALE, 0); vector<Mat> planes; split(tmp, planes); magnitude(planes[0], planes[1], planes[0]); Fourier = planes[0]; } } void psf2otf(Mat& src, int rows, int cols, Mat& dst) { Mat src_fill = Mat::zeros(rows, cols, CV_32FC1); Mat src_fill_out = Mat::zeros(rows, cols, CV_32FC1); for (int i = 0; i < src.rows; i++) { float* data = src_fill.ptr<float>(i); float* data2 = src.ptr<float>(i); for (int j = 0; j < src.cols; j++) { data[j] = data2[j]; } } Size size; size.height = rows; size.width = cols; circshift(src_fill, -int(src.rows / 2), -int(src.cols / 2), src_fill_out); dst = fft2(src_fill_out, size); return; } void circshift(Mat& img, int dw, int dh, Mat& dst) { int rows = img.rows; int cols = img.cols; dst = img.clone(); if (dw < 0 && dh < 0) { for (int i = 0; i < rows; i++) { int t = i + dw; if (t >= 0) { float* data = img.ptr<float>(i); float* data2 = dst.ptr<float>(t); for (int j = 0; j < cols; j++) { data2[j] = data[j]; } } else { float* data = img.ptr<float>(i); float* data2 = dst.ptr<float>(dst.rows + t); for (int j = 0; j < cols; j++) { data2[j] = data[j]; } } } for (int j = 0; j < cols; j++) { int t = j + dh; if (t >= 0) { for (int i = 0; i < rows; i++) { float* data = img.ptr<float>(i); float* data2 = dst.ptr<float>(i); data2[t] = data[j]; } } else { for (int i = 0; i < rows; i++) { float* data = img.ptr<float>(i); float* data2 = dst.ptr<float>(i); data2[dst.cols + t] = data[j]; } } } } return; } Mat sign(Mat& input) { int nr = input.rows; int nl = input.cols; Mat output(nr, nl, CV_32FC1); for (int i = 0; i < nr; i++) { const float* inData = input.ptr<float>(i); float* outData = output.ptr<float>(i); for (int j = 0; j < nl; j++) { if (*inData > 0) { *outData = 1; } else if (*inData < 0) { *outData = -1; } else { *outData = 0; } outData++; } } return output; } Mat fliplr(Mat& input) { int nr = input.rows; int nl = input.cols; Mat output(nr, nl, CV_32FC1); float temp; for (int i = 0; i < nr; i++) { for (int j = 0; j < (nl - 1) / 2 + 1; j++) { temp = input.at<float>(i, j); output.at<float>(i, j) = input.at<float>(i, nl - j - 1); output.at<float>(i, nl - j - 1) = temp; } } return output; } Mat flipud(Mat& input) { int nr = input.rows; int nl = input.cols; Mat output(nr, nl, CV_32FC1); float temp; for (int i = 0; i < (nr - 1) / 2 + 1; i++) { for (int j = 0; j < nl; j++) { temp = input.at<float>(i, j); output.at<float>(i, j) = input.at<float>(nr - 1 - i, j); output.at<float>(nr - 1 - i,j) = temp; } } return output; } Mat recover(Mat& srcimg, Mat& t, float *airlight, float delta) { int nr = srcimg.rows, nl = srcimg.cols; float tnow = t.at<float>(0, 0); float t0 = 0.0001; 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); pow(tnow, delta); for (int i = 0; i < 3; i++) { val = (*srcPtr++ - airlight[i]) / tnow + airlight[i]; *outPtr++ = val; } } } return finalimg; }
为了节省篇幅，在此仅举一例实验结果：
雾天图像 原始传输函数 正则化后的传输函数 复原结果
