代碼開源在 https://github.com/BBuf/Image-processing-algorithm，感興趣給我來個星星唄。

1. 前言

這是OpenCV圖像處理算法樸素實現專欄的第17篇文章。今天為大家帶來一篇之前看到的用于單幅圖像去霧的算法，作者來自清華大學，論文原文見附錄。

2. 霧天退化模型

之前在介紹何凱明博士的暗通道去霧論文OpenCV圖像處理專欄六 | 來自何凱明博士的暗通道去霧算法(CVPR 2009最佳論文)的時候已經講到了這個霧天退化模型，我們這里再來回顧一下。在計算機視覺領域，通常使用霧天圖像退化模型來描述霧霾等惡劣天氣條件對圖像造成的影響，該模型是McCartney首先提出。該模型包括衰減模型和環境光模型兩部分。模型表達式為：

$I(x)=J(x)e^{?rd(x)}+A(1?e^{?rd(x)})...........................(1)$

其中，$x$是圖像像素的空間坐標，$H$是觀察到的有霧圖像，$F$是待恢復的無霧圖像，$r$表示大氣散射系數，$d$代表景物深度，$A$是全局大氣光，通常情況下假設為全局常量，與空間坐標$x$無關。

公式(1)中的$e^{?r(dx)}$表示坐標空間$x$處的透射率，我們使用$t(x)$來表示透射率，于是得到公式（2）：

$I(x)=J(x)t(x)+A(1?t(x)).............................(2)$

由此可見，圖像去霧過程就是根據$I(x)$求解$J(x)$的過程。要求解出$J(x)$，還需要根據$I(x)$求解出透射率$t(x)$和全局大氣光$A$。實際上，所有基于霧天退化模型的去霧算法就是是根據已知的有霧圖像$I(x)$求解出透射率$t(x)$和全局大氣光$A$。

對于暗通道去霧算法來說，先從暗原色通道中選取最亮的0.1%比例的像素點，然后選取原輸入圖像中這些像素具有的最大灰度值作為全局大氣光值$A$。RGB三通道中每一個通道都有一個大氣光值。

然后根據公式(2)可以得出：

$t(x)=\frac{A?I(x)}{A?J(x)}...........................(3)$

首先可以確定的是$t(x)$的范圍是$[0,1]$，$I(x)$的范圍是$[0,255]$，$J(x)$的范圍是$[0,255]$。$A$和$I(x)$是已知的，可以根據$J(x)$的范圍從而確定$t(x)$的范圍。已知的條件有：

$0<=J(x)<=255,0<=I(x)<=A,0<=J(x)<=A,0<=t(x)<=1...............(4)$

$=>t(x)>=\frac{A?I(x)}{A?0}=\frac{A?I(x)}{A}=1?\frac{I(x)}{A}..................(5)$

根據(4)和(5)推出：
$1?\frac{I(x)}{A}<=t(x)<=1..................................(6)$

因此初略估計透射率的計算公式：

$t(x)=1?\frac{I(x)}{A}...........................................(7)$

最后為了保證圖片的自然性，增加一個參數$w$來調整透射率 ：

$t(x)=1?w\frac{I(x)}{A}..............................(8)$

好了，上面復習完了何凱明博士的暗通道去霧，我們一起來看看清華大學這篇論文。

3. 算法流程

實際上有了這個算法流程就可以寫出代碼了，不過為了加深理解可以看下面的一些推導。

4. 一些推導

我們知道去霧的步驟主要就是估計全局大氣光值$A$和透射率$t(x)$，因此，本文就是根據輸入圖像估計$A$和$L(x)$（這篇論文使用了$L(x)$來代替$A(1?t(x))$），然后根據霧天退化模型求取去霧后的圖像。

4.1 估計透射率t(x)

從第二節的介紹我們知道

$t(x)>=1?\frac{H(x)}{A}...............................(2)$

然后這篇論文使用了$L(x)$來代替$A(1?t(x))$，即：

$L(x)=A(1?t(x)).............................(1)$。

我們取$H(x)$三個通道的最小值并記為：

$M(x)=mi{n}_{c\in r,g,b}(H^c(x)).......................(3)$

所以公式2變換為

$t(x)>=1?\frac{M(x)}{A}...................................(4)$

對公式(4)右邊進行均值濾波：

$meadia{n}_{s_a}(1?\frac{M(x)}{A})=1?\frac{media{n}_{s_a}(M(x))}{A}=1?\frac{{\sum }_{y\in \Omega (x)M(y)}}{A{s}_a^2}...........(5)$

其中$s_a$代表均值濾波的窗口大小，$\Omega (x)$表示像素$x$的$s_a\times s_a$的鄰域。

均值濾波后的結果可以反映$t(x)$的大致趨勢，但與真實的$t(x)$還差一定的絕對值，因此，我們先得出透射率的粗略估計值：

$t(x)=1?\frac{M_{ave}(x)}{A}+\phi \frac{M_{ave}(x)}{A}=1?\delta \frac{M_{ave}(x)}{A}................(6)$

其中$M_{ave}(x)=media{n}_{s_a}(M(x)),\delta =1?\phi ,\phi \in [0,1]$，因此$\delta \in [0,1]$。

為了防止去霧后圖像出現整體畫面偏暗，這里根據圖像的均值來調整$\delta$，即：

$\delta =\rho m_{ave}$

其中$m_{ave}$是$M(x)$中所有元素的均值，$\rho$是調節因子。

因此可以得到透射率的計算公式:

$t(x)=max(1?min(\rho m_{av},0.9)\frac{M_{ave}(x)}{A},1?\frac{M(x)}{A})...................(7)$

結合公式(1)推出：

$L(x)=min(min(\rho m_{av},0.9)M_{ave}(x),M(x))..........(8)$。

4.2 估計全球大氣光值

公式(5)中第一個等式左側的表達式取值范圍為$[0,1]$，由此得出

$A>=max(M_{ave}(x))$

一般情況下又存在

$A<=max(ma{x}_{c\in r,g,b(H^c(x))})$

(KaiMing He的暗通道先驗理論)。這樣就初步確定了全局大氣光的范圍，為了能快速獲取全局大氣光，文章直接取兩者的平均值作為全局大氣光值，即：

$A=1/2(max(H(x))+max(M_{ave}(x)))$…(9)。

然后大氣光值$A$和$L(x)$都搞定了，那么帶入算法流程中的最后一個公式就可以獲取最后的圖像了。

5. 代碼實現

下面是代碼實現。

#include <opencv2/opencv.hpp> #include <iostream> #include <algorithm> #include <vector> using namespace cv; using namespace std;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;//get Mdouble 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(r, g), b);temp2 = max(max(r, g), b);M.at<uchar>(i, j) = temp1;M_max.at<uchar>(i, j) = temp2;sum += temp1;}}m_av = sum / (row * col * 255);eps = 0.85 / m_av;boxFilter(M, M_ave, CV_8UC1, Size(51, 51));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 main() {Mat src = imread("F:\\fog\\1.jpg");Mat dst = dehaze(src);cv::imshow("origin", src);cv::imshow("result", dst);cv::imwrite("F:\\fog\\res.jpg", dst);waitKey(0);return 0; }

6. 結果

7. 結論

算法里面有2個參數可以自己調節，濾波的半徑和$\rho$。具體如何調節？我就不放在這里說了，這個算法后面會在我的新專題里面進行一遍優化，到時候再來回答這個問題。如果你迫切需要這個算法的實現或者對它感興趣，可以自己嘗試調整這兩個參數獲得想要的效果。這里的均值濾波也可以換成我們之前講的Side Window Filter說不定可以獲得更好的效果。

8. 參考

• https://blog.csdn.net/u013684730/article/details/76640321
• https://www.cnblogs.com/Imageshop/p/3410279.html

---

歡迎關注GiantPandaCV, 在這里你將看到獨家的深度學習分享，堅持原創，每天分享我們學習到的新鮮知識。( ? ?ω?? )?

有對文章相關的問題，或者想要加入交流群，歡迎添加BBuf微信：

總結

以上是生活随笔為你收集整理的OpenCV图像处理专栏十七 | 清华大学《基于单幅图像的快速去雾》C++复现（有一定工程意义）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。