前言

先讓我們來膜拜何凱明大佬這篇獲得CVPR 2009的Best Paper Award論文。這篇論文的靈感來自自於作者2個觀察，第一個是在3D遊戲中的霧使得作者堅信人眼有特殊的東西去感知霧，而不僅僅是靠對比度。第二個是作者觀察了之前的一篇去霧論文:《Single Image Dehazing》，發現這篇paper中的Dark Object Subtraction可以處理均勻的霧，但是非均勻的就處理不好了，所以作者在區域性使用了Dark Object Subtraction，結果得到了驚人的效果。。。觀察能力無敵

原理

• 暗通道先驗：首先說在絕大多數非天空的區域性區域裡，某一些畫素總會有至少一個顏色通道具有很低的值，也就是說該區域光強是一個很小的值。所以給暗通道下了個數學定義，對於任何輸入的影象J,其暗通道可以用下面的公式來表示：
  ，其中$J^C$表示彩色影象每個通道，$Ω(x)$表示以畫素X為中心的一個視窗。要求暗通道的影象是比較容易的，先求出每個畫素在3個通道的最小值，存到一個二維Mat中(灰度圖)，然後做一個最小值濾波，濾波的半徑由視窗大小決定，這裡視窗大小為WindowSize，有公式表示WindowsSize=2*Radius+1。
• 暗通道先驗指出的結論，這個我不知道如何證明，不過作者給出了幾個原因。a)汽車、建築物和城市中玻璃窗戶的陰影，或者是樹葉、樹與岩石等自然景觀的投影；b)色彩鮮豔的物體或表面，在RGB的三個通道中有些通道的值很低（比如綠色的草地／樹／植物，紅色或黃色的花朵／葉子，或者藍色的水面）；c)顏色較暗的物體或者表面，例如灰暗色的樹幹和石頭。總之，自然景物中到處都是陰影或者彩色，這些景物的影象的暗原色總是很灰暗的。作者在論文中，統計了5000多副影象的特徵，也都基本符合這個先驗。因此，我們可以認為其實一條定理。
• 基於這個先驗，接下來就是該論文中最核心的部分了，首先，在計算機視覺和影象處理中，公式(1)這個霧形成模型被廣泛的應用：$I(x)=J(x)t(x)+A(1-t(x))$，其中I(x)是我們待處理的影象，J(x)是我們要恢復的沒有霧的影象，A是全球大氣光成分，t(x)為透射率。現在已知了I(X)，我們需要求取J(X)，顯然這個不定方程有無數解，所以需要定義一些先驗。將上式處理變形得到:$\frac{I^c(x)}{A^c}=t(x)\frac{J^c(x)}{A^c}+1-t(x)$
  ，其中上標C代表R,G,B3個通道。然後假設在每一個視窗中透射率t(x)是一個常數，定義為$t‘(x)$並且A已經值已經給定，然後對這個式子左右同時取2次最小值，得到下面的式子：那個t’(x)就是公式(8)中那個t^{}~(x)。上式中，J是待求的無霧的影象，根據前述的暗原色先驗理論有： 因此，可推匯出：

把式(10)帶入式(8)中，得到：這就是透射率$t'(x)$的預估值。 在現實生活中，即使是晴天白雲，空氣中也存在著一些顆粒，因此，看遠處的物體還是能感覺到霧的影響，另外，霧的存在讓人類感到景深的存在，因此，有必要在去霧的時候保留一定程度的霧，這可以通過在式（11）中引入一個在[0,1] 之間的因子，則式（11）修正為： 本文中所有的測試結果依賴於： ω=0.95。

• 上述的推導是基於A已知的情況下，然而事實是A還不知道呢？A怎麼計算呢？在實際中，我們可以藉助於暗通道圖來從有霧影象中獲取該值。具體步驟如下：(1)從暗通道圖中按照亮度的大小取前0.1%的畫素。(2)在這些位置中，在原始有霧影象I中尋找對應的具有最高亮度的點的值，作為A值。 到這一步，我們就可以進行無霧影象的恢復了。由$I(x)=J(x)t(x)+A(1-t(x))$，推出 $J(x)=(I(x)-A)/t(x)+A$，現在I,A,t都已經求得了，因此，完全可以進行J的計算。當投射圖t 的值很小時，會導致J的值偏大，從而使得影象整體向白場過度，因此一般可設定一閾值T0，當t值小於T0時，令t=T0，本文中所有效果圖均以T0=0.1為標準計算。
• 最終的結論為：
• 現在準備用C++拿這個公式去復現一下論文，之後會放上原始碼。

按照上面的公式復現了論文，給幾張圖片測試的結果

原圖2
程式碼如下：

#include <opencv2/opencv.hpp>
#include <iostream>
#include <algorithm>
#include <vector>
using namespace cv;
using namespace std;

int rows, cols;
//獲取最小值矩陣
int **getMinChannel(cv::Mat img){
    rows = img.rows;
    cols = img.cols;
    if(img.channels() != 3){
        fprintf(stderr, "Input Error!");
        exit(-1);
    }
    int **imgGray;
    imgGray = new int *[rows];
    for(int i = 0; i < rows; i++){
        imgGray[i] = new int [cols];
    }
    for(int i = 0; i < rows; i++){
        for(int j = 0; j < cols; j++){
            int loacalMin = 255;
            for(int k = 0; k < 3; k++){
                if(img.at<Vec3b>(i, j)[k] < loacalMin){
                    loacalMin = img.at<Vec3b>(i, j)[k];
                }
            }
            imgGray[i][j] = loacalMin;
        }
    }
    return imgGray;
}

//求暗通道
int **getDarkChannel(int **img, int blockSize = 3){
    if(blockSize%2 == 0 || blockSize < 3){
        fprintf(stderr, "blockSize is not odd or too small!");
        exit(-1);
    }
    //計算pool Size
    int poolSize = (blockSize - 1) / 2;
    int newHeight = rows + blockSize - 1;
    int newWidth = cols + 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 < rows && j < cols){
                imgMiddle[i][j] = img[i][j];
            }else{
                imgMiddle[i][j] = 255;
            }
        }
    }
    int **imgDark;
    imgDark = new int *[rows];
    for(int i = 0; i < rows; i++){
        imgDark[i] = new int [cols];
    }
    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 **darkChannel, cv::Mat img, bool meanMode = false, float percent = 0.001){
    int size = rows * cols;
    std::vector <node> nodes;
    for(int i = 0; i < rows; i++){
        for(int j = 0; j < cols; j++){
            node tmp;
            tmp.x = i, tmp.y = j, tmp.val = darkChannel[i][j];
            nodes.push_back(tmp);
        }
    }
    sort(nodes.begin(), nodes.end());
    int atmosphericLight = 0;
    if(int(percent*size) == 0){
        for(int i = 0; i < 3; i++){
            if(img.at<Vec3b>(nodes[0].x, nodes[0].y)[i] > atmosphericLight){
                atmosphericLight = img.at<Vec3b>(nodes[0].x, nodes[0].y)[i];
            }
        }
    }
    //開啟均值模式
    if(meanMode == true){
        int sum = 0;
        for(int i = 0; i < int(percent*size); i++){
            for(int j = 0; j < 3; j++){
                sum = sum + img.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(img.at<Vec3b>(nodes[i].x, nodes[i].y)[j] > atmosphericLight){
                atmosphericLight = img.at<Vec3b>(nodes[i].x, nodes[i].y)[j];
            }
        }
    }
    return atmosphericLight;
}

//恢復原影象
// Omega 去霧比例 引數
//t0 最小透射率值
cv::Mat getRecoverScene(cv::Mat img, float omega=0.95, float t0=0.1, int blockSize=15, bool meanModel=false, float percent=0.001){
    int** imgGray = getMinChannel(img);
    int **imgDark = getDarkChannel(imgGray, blockSize=blockSize);
    int atmosphericLight = getGlobalAtmosphericLightValue(imgDark, img, meanModel=meanModel, percent=percent);
    float **imgDark2, **transmission;
    imgDark2 = new float *[rows];
    for(int i = 0; i < rows; i++){
        imgDark2[i] = new float [cols];
    }
    transmission = new float *[rows];
    for(int i = 0; i < rows; i++){
        transmission[i] = new float [cols];
    }
    for(int i = 0; i < rows; i++){
        for(int j = 0; j < cols; 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;
            }
        }
    }
    cv::Mat dst(img.rows, img.cols, CV_8UC3);
    for(int channel = 0; channel < 3; channel++){
        for(int i = 0; i < rows; i++){
            for(int j = 0; j < cols; j++){
                int temp = (img.at<Vec3b>(i, j)[channel] - atmosphericLight) / transmission[i][j] + atmosphericLight;
                if(temp > 255){
                    temp = 255;
                }
                if(temp < 0){
                    temp = 0;
                }
                dst.at<Vec3b>(i, j)[channel] = temp;
            }
        }
    }
    return dst;
}

int main(){
    cv::Mat src = cv::imread("/home/zxy/CLionProjects/Acmtest/4.jpg");
    rows = src.rows;
    cols = src.cols;
    cv::Mat dst = getRecoverScene(src);
    cv::imshow("origin", src);
    cv::imshow("result", dst);
    cv::imwrite("../zxy.jpg", dst);
    waitKey(0);
}

效果確實和論文的效果幾乎一致，在640*480的影象上的執行速度是250ms，這個時間如果想要實時的話仍然需要走到優化的道路上。已經有一篇可以實時的部落格了，大家可以去看https://www.cnblogs.com/Imageshop/p/3281703.html