OpenCV图像处理专栏十七 | 清华大学《基于单幅图像的快速去雾》C++复现(有一定工程意义)
代碼開源在 https://github.com/BBuf/Image-processing-algorithm,感興趣給我來(lái)個(gè)星星唄。
1. 前言
這是OpenCV圖像處理算法樸素實(shí)現(xiàn)專欄的第17篇文章。今天為大家?guī)?lái)一篇之前看到的用于單幅圖像去霧的算法,作者來(lái)自清華大學(xué),論文原文見附錄。
2. 霧天退化模型
之前在介紹何凱明博士的暗通道去霧論文OpenCV圖像處理專欄六 | 來(lái)自何凱明博士的暗通道去霧算法(CVPR 2009最佳論文)的時(shí)候已經(jīng)講到了這個(gè)霧天退化模型,我們這里再來(lái)回顧一下。在計(jì)算機(jī)視覺(jué)領(lǐng)域,通常使用霧天圖像退化模型來(lái)描述霧霾等惡劣天氣條件對(duì)圖像造成的影響,該模型是McCartney首先提出。該模型包括衰減模型和環(huán)境光模型兩部分。模型表達(dá)式為:
I(x)=J(x)e?rd(x)+A(1?e?rd(x))...........................(1)I(x)=J(x)e^{-rd(x)}+A(1-e^{-rd(x)})...........................(1)I(x)=J(x)e?rd(x)+A(1?e?rd(x))...........................(1)
其中,xxx是圖像像素的空間坐標(biāo),HHH是觀察到的有霧圖像,FFF是待恢復(fù)的無(wú)霧圖像,rrr表示大氣散射系數(shù),ddd代表景物深度,AAA是全局大氣光,通常情況下假設(shè)為全局常量,與空間坐標(biāo)xxx無(wú)關(guān)。
公式(1)中的e?r(dx)e^{-r(dx)}e?r(dx)表示坐標(biāo)空間xxx處的透射率,我們使用t(x)t(x)t(x)來(lái)表示透射率,于是得到公式(2):
I(x)=J(x)t(x)+A(1?t(x)).............................(2)I(x)=J(x)t(x)+A(1-t(x)).............................(2)I(x)=J(x)t(x)+A(1?t(x)).............................(2)
由此可見,圖像去霧過(guò)程就是根據(jù)I(x)I(x)I(x)求解J(x)J(x)J(x)的過(guò)程。要求解出J(x)J(x)J(x),還需要根據(jù)I(x)I(x)I(x)求解出透射率t(x)t(x)t(x)和全局大氣光AAA。實(shí)際上,所有基于霧天退化模型的去霧算法就是是根據(jù)已知的有霧圖像I(x)I(x)I(x)求解出透射率t(x)t(x)t(x)和全局大氣光AAA。
對(duì)于暗通道去霧算法來(lái)說(shuō),先從暗原色通道中選取最亮的0.1%比例的像素點(diǎn),然后選取原輸入圖像中這些像素具有的最大灰度值作為全局大氣光值AAA。RGB三通道中每一個(gè)通道都有一個(gè)大氣光值。
然后根據(jù)公式(2)可以得出:
t(x)=A?I(x)A?J(x)...........................(3)t(x)=\frac{A-I(x)}{A-J(x)}...........................(3)t(x)=A?J(x)A?I(x)?...........................(3)
首先可以確定的是t(x)t(x)t(x)的范圍是[0,1][0, 1][0,1],I(x)I(x)I(x)的范圍是[0,255][0,255][0,255],J(x)J(x)J(x)的范圍是[0,255][0, 255][0,255]。AAA和I(x)I(x)I(x)是已知的,可以根據(jù)J(x)J(x)J(x)的范圍從而確定t(x)t(x)t(x)的范圍。已知的條件有:
0<=J(x)<=255,0<=I(x)<=A,0<=J(x)<=A,0<=t(x)<=1...............(4)0<=J(x)<=255, 0<=I(x)<=A,0<=J(x)<=A,0<=t(x)<=1...............(4)0<=J(x)<=255,0<=I(x)<=A,0<=J(x)<=A,0<=t(x)<=1...............(4)
=>t(x)>=A?I(x)A?0=A?I(x)A=1?I(x)A..................(5)=>t(x)>=\frac{A-I(x)}{A-0}=\frac{A-I(x)}{A}=1-\frac{I(x)}{A}..................(5)=>t(x)>=A?0A?I(x)?=AA?I(x)?=1?AI(x)?..................(5)
根據(jù)(4)和(5)推出:
 1?I(x)A<=t(x)<=1..................................(6)1-\frac{I(x)}{A}<=t(x)<=1..................................(6)1?AI(x)?<=t(x)<=1..................................(6)
因此初略估計(jì)透射率的計(jì)算公式:
t(x)=1?I(x)A...........................................(7)t(x)=1-\frac{I(x)}{A}...........................................(7)t(x)=1?AI(x)?...........................................(7)
最后為了保證圖片的自然性,增加一個(gè)參數(shù)www來(lái)調(diào)整透射率 :
t(x)=1?wI(x)A..............................(8)t(x)=1-w\frac{I(x)}{A}..............................(8)t(x)=1?wAI(x)?..............................(8)
好了,上面復(fù)習(xí)完了何凱明博士的暗通道去霧,我們一起來(lái)看看清華大學(xué)這篇論文。
3. 算法流程
 實(shí)際上有了這個(gè)算法流程就可以寫出代碼了,不過(guò)為了加深理解可以看下面的一些推導(dǎo)。
4. 一些推導(dǎo)
我們知道去霧的步驟主要就是估計(jì)全局大氣光值AAA和透射率t(x)t(x)t(x),因此,本文就是根據(jù)輸入圖像估計(jì)AAA和L(x)L(x)L(x)(這篇論文使用了L(x)L(x)L(x)來(lái)代替A(1?t(x))A(1-t(x))A(1?t(x))),然后根據(jù)霧天退化模型求取去霧后的圖像。
4.1 估計(jì)透射率t(x)
從第二節(jié)的介紹我們知道
t(x)>=1?H(x)A...............................(2)t(x)>=1-\frac{H(x)}{A}...............................(2)t(x)>=1?AH(x)?...............................(2)
然后這篇論文使用了L(x)L(x)L(x)來(lái)代替A(1?t(x))A(1-t(x))A(1?t(x)),即:
L(x)=A(1?t(x)).............................(1)L(x)=A(1-t(x)).............................(1)L(x)=A(1?t(x)).............................(1)。
我們?nèi)?span id="vt6mr5x" class="katex--inline">H(x)H(x)H(x)三個(gè)通道的最小值并記為:
M(x)=minc∈r,g,b(Hc(x)).......................(3)M(x)=min_{c\in r,g,b}(H^c(x)).......................(3)M(x)=minc∈r,g,b?(Hc(x)).......................(3)
所以公式2變換為
t(x)>=1?M(x)A...................................(4)t(x)>=1-\frac{M(x)}{A}...................................(4)t(x)>=1?AM(x)?...................................(4)
對(duì)公式(4)右邊進(jìn)行均值濾波:
meadiansa(1?M(x)A)=1?mediansa(M(x))A=1?∑y∈Ω(x)M(y)Asa2...........(5)meadian_{s_a}(1-\frac{M(x)}{A})=1-\frac{median_{s_a}(M(x))}{A}=1-\frac{\sum_{y\in\Omega(x)M(y)}}{As_a^2}...........(5)meadiansa??(1?AM(x)?)=1?Amediansa??(M(x))?=1?Asa2?∑y∈Ω(x)M(y)??...........(5)
其中sas_asa?代表均值濾波的窗口大小,Ω(x)\Omega(x)Ω(x)表示像素xxx的sa×sas_a\times s_asa?×sa?的鄰域。
均值濾波后的結(jié)果可以反映t(x)t(x)t(x)的大致趨勢(shì),但與真實(shí)的t(x)t(x)t(x)還差一定的絕對(duì)值,因此,我們先得出透射率的粗略估計(jì)值:
t(x)=1?Mave(x)A+φMave(x)A=1?δMave(x)A................(6)t(x)=1-\frac{M_{ave}(x)}{A}+\varphi\frac{M_{ave}(x)}{A}=1-\delta\frac{M_{ave}(x)}{A}................(6)t(x)=1?AMave?(x)?+φAMave?(x)?=1?δAMave?(x)?................(6)
其中Mave(x)=mediansa(M(x)),δ=1?φ,φ∈[0,1]M_{ave}(x)=median_{s_a}(M(x)),\delta=1-\varphi,\varphi\in[0,1]Mave?(x)=mediansa??(M(x)),δ=1?φ,φ∈[0,1],因此δ∈[0,1]\delta \in[0,1]δ∈[0,1]。
為了防止去霧后圖像出現(xiàn)整體畫面偏暗,這里根據(jù)圖像的均值來(lái)調(diào)整δ\deltaδ,即:
δ=ρmave\delta=\rho m_{ave}δ=ρmave?
其中mavem_{ave}mave?是M(x)M(x)M(x)中所有元素的均值,ρ\rhoρ是調(diào)節(jié)因子。
因此可以得到透射率的計(jì)算公式:
t(x)=max(1?min(ρmav,0.9)Mave(x)A,1?M(x)A)...................(7)t(x)=max(1-min(\rho m_{av}, 0.9)\frac{M_{ave}(x)}{A}, 1-\frac{M(x)}{A})...................(7)t(x)=max(1?min(ρmav?,0.9)AMave?(x)?,1?AM(x)?)...................(7)
結(jié)合公式(1)推出:
L(x)=min(min(ρmav,0.9)Mave(x),M(x))..........(8)L(x)=min(min(\rho m_{av}, 0.9)M_{ave}(x), M(x))..........(8)L(x)=min(min(ρmav?,0.9)Mave?(x),M(x))..........(8)。
4.2 估計(jì)全球大氣光值
公式(5)中第一個(gè)等式左側(cè)的表達(dá)式取值范圍為[0,1][0, 1][0,1],由此得出
A>=max(Mave(x))A>=max(M_{ave}(x))A>=max(Mave?(x))
一般情況下又存在
A<=max(maxc∈r,g,b(Hc(x)))A<=max(max_{c\in r,g,b(H^c(x))})A<=max(maxc∈r,g,b(Hc(x))?)
(KaiMing He的暗通道先驗(yàn)理論)。這樣就初步確定了全局大氣光的范圍,為了能快速獲取全局大氣光,文章直接取兩者的平均值作為全局大氣光值,即:
A=1/2(max(H(x))+max(Mave(x)))A = 1/2(max(H(x))+max(M_{ave}(x)))A=1/2(max(H(x))+max(Mave?(x)))…(9)。
然后大氣光值AAA和L(x)L(x)L(x)都搞定了,那么帶入算法流程中的最后一個(gè)公式就可以獲取最后的圖像了。
5. 代碼實(shí)現(xiàn)
下面是代碼實(shí)現(xiàn)。
#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. 結(jié)果
7. 結(jié)論
算法里面有2個(gè)參數(shù)可以自己調(diào)節(jié),濾波的半徑和ρ\rhoρ。具體如何調(diào)節(jié)?我就不放在這里說(shuō)了,這個(gè)算法后面會(huì)在我的新專題里面進(jìn)行一遍優(yōu)化,到時(shí)候再來(lái)回答這個(gè)問(wèn)題。如果你迫切需要這個(gè)算法的實(shí)現(xiàn)或者對(duì)它感興趣,可以自己嘗試調(diào)整這兩個(gè)參數(shù)獲得想要的效果。這里的均值濾波也可以換成我們之前講的Side Window Filter說(shuō)不定可以獲得更好的效果。
8. 參考
- https://blog.csdn.net/u013684730/article/details/76640321
 - https://www.cnblogs.com/Imageshop/p/3410279.html
 
歡迎關(guān)注GiantPandaCV, 在這里你將看到獨(dú)家的深度學(xué)習(xí)分享,堅(jiān)持原創(chuàng),每天分享我們學(xué)習(xí)到的新鮮知識(shí)。( ? ?ω?? )?
有對(duì)文章相關(guān)的問(wèn)題,或者想要加入交流群,歡迎添加BBuf微信:
總結(jié)
以上是生活随笔為你收集整理的OpenCV图像处理专栏十七 | 清华大学《基于单幅图像的快速去雾》C++复现(有一定工程意义)的全部?jī)?nèi)容,希望文章能夠幫你解決所遇到的問(wèn)題。
                            
                        - 上一篇: 基于C#.NET的高端智能化网络爬虫(二
 - 下一篇: 小黄鸭调试法-程序猿修炼之道