前言
先让我们来膜拜何凯明大佬这篇获得CVPR 2009的Best Paper Award论文。这篇论文的灵感来自自于作者2个观察，第一个是在3D游戏中的雾使得作者坚信人眼有特殊的东西去感知雾，而不仅仅是靠对比度。第二个是作者观察了之前的一篇去雾论文:《Single Image Dehazing》，发现这篇paper中的Dark Object Subtraction可以处理均匀的雾，但是非均匀的就处理不好了，所以作者在局部使用了Dark Object Subtraction，结果得到了惊人的效果。。。观察能力无敌
原理
- 暗通道先验：首先说在绝大多数非天空的局部区域里，某一些像素总会有至少一个颜色通道具有很低的值，也就是说该区域光强是一个很小的值。所以给暗通道下了个数学定义，对于任何输入的图像J,其暗通道可以用下面的公式来表示：
，其中 JC表示彩色图像每个通道， Ω(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)，显然这个不定方程有无数解，所以需要定义一些先验。将上式处理变形得到: AcIc(x)=t(x)AcJc(x)+1−t(x)，其中上标C代表R,G,B3个通道。然后假设在每一个窗口中透射率t(x)是一个常数，定义为 t‘(x)并且A已经值已经给定，然后对这个式子左右同时取2次最小值，得到下面的式子：
那个t’(x)就是公式(8)中那个t^{}~(x)。上式中，J是待求的无雾的图像，根据前述的暗原色先验理论有：
因此，可推导出：
- 上述的推导是基于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++拿这个公式去复现一下论文，之后会放上源码。
参考文章：https://www.cnblogs.com/Imageshop/p/3281703.html
按照上面的公式复现了论文，给几张图片测试的结果
代码如下：
#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
京公网安备 11010502036488号