Guided Image Filtering --引导图像滤波 解析+图像去雾应用

论文地址以及代码：http://kaiminghe.com/eccv10/

Abstract: 这篇文章中，作者提出了一种新颖的图像滤波算法。这种方法的灵感来自于局部线性模型。本方法引入了一个导向滤波图，这个导向滤波图可以是图片本身也可以是另一个图片，然后算法利用导向滤波图的纹路来对目标图像进行滤波。这个方法类似于双边滤波器是一个边缘保留平滑算子，但是在边缘保留的表现上优于双边滤波器。此外，他还可以将导向图的结构保留到滤波的结果中，可以用来去雾、羽化等。该方法处理速度快，应是现在最快的边缘保留滤波器。
一、图像滤波的简要介绍
说在前面：关于权重与边缘保留的关系，这篇文章解释的挺清楚的：http://blog.csdn.net/baimafujinji/article/details/74750283

首先，我们定义一个线性转换的滤波器：



这个公式的意思是说，输出像素

的值是输入图像p以像素j为中心的一个小窗口和其对应像素的权重相乘后的和。这是最简单的滤波器，根据不同的权重有不同的作用。例如高斯滤波器，它的权重服从二维的高斯分布。也就是说，离中心像素越近的点就拥有更高的权重。这样的滤波器拥有很好的平滑作用，但也因此会使边缘变得模糊。
再例如双边滤波器的权重，如公式（2）所示，相比起高斯滤波就更加的科学，也拥有更好的边缘保留能力。



其中k表示了窗口中的像素数目，

表示了两个像素之间的空间距离，

表示了两个像素之间的颜色差距，



这两个是常数系数。这个权重的意思是说，空间距离越近，颜色越相近，这个像素所占的权重就越大，反之则越小。由于边缘两边的像素通常很接近，但颜色相差却比较大，这样的权重就会相对小一些，也就达到了边缘保留的目的。
二、导向图滤波
符号定义：

I：引导图像
p：输入图像
q：输出图像
    公式定义：

，其中，



都是常数系数，要通过计算获得。
    这个公式非常简单却具有边缘保留的能力以及不逊色的平滑去噪能力。

        1. 首先来看一下为什么这个公式具有边缘保留能力：
        对公式（4）求导可以得到：

，也就是说，导向图I有梯度的地方，输出图像q就有梯度，即，他们拥有相同的边缘。
        2. 如何平滑去噪
        在（4）式的线性模型下，我们在保留边缘的时候还要实现去噪，用公式表达的话就是：
     

其中的

就表示噪声。
        由（4）（5）可得我们需要最小化的损失函数为：
        


        其中的

为惩罚项。
        对（6）式求解的结果如下：
        


        也就是说，在系数



为（7）（8）式的时候，整个线性模型不仅能实现平滑去噪，还能实现边缘保留。
       

        然而，上面的公式都是针对一个像素点而言的。当我们真正处理起来的时候一般是使用滑动窗口的形式，于是，线性模型就变成了如下的形式：


        把前面的系数放进括号里面，就可以简单表达为：


        这里需要注意的是，经过（9）（10）式的变换以后，

不在是

的简单缩放了。因为

空间可变，现在相当于被一个均值滤波器的滤过，是一个平均值，所以输出图像的梯度将会比输入图像的梯度小一些。在这种情况下，

，也就意味着I中突然的梯度变化大部分能被q保留。
三、运用到图像去雾中
    导向图滤波还能应用到图像去雾中去。上一篇文章《Single Image Haze Removal Using Dark Channel Prior》中，求出的t(x)只是粗略估计，而文章中的softmatting的方法过于复杂，时间复杂度也很高。在这篇文章中提出的导向图滤波就能代替softmatting的方法精细化t(x）。
    在图像去雾中，导向图就是有雾图像的暗通道图，输入图像就是有雾图像，算法如下图所示：
    


    四、代码
opencv版#include<opencv2/opencv.hpp>
#include<opencv2/highgui/highgui.hpp>
#include<opencv2/imgproc/imgproc.hpp>
//#include<iostream>
using namespace cv;

typedef struct Pixel
{
int x;
int y;
int value;
}Pixel;

inline bool cmp(Pixel a, Pixel b)
{
return a.value>b.value;
}

//最小值滤波
//dark_img = minfliter(dark_img, windowsize);
Mat minfilter(Mat& D, int windowsize) {
int dr = D.rows;
int dc = D.cols;
int r = (windowsize - 1) / 2; //windowsize是奇数
Mat dark_img(dr, dc, CV_8UC1); //暗通道图
dark_img = D.clone();
for (int i = r;i <= dr - r - 1;i++) {
for (int j = r;j <= dc - r - 1;j++) {
int min = 255;
for (int m = i - r;m <= i + r;m++) {
for (int n = j - r;n <= j + r;n++) {
if (D.at<uchar>(m, n) < min)
min = D.at<uchar>(m, n);
}
}
dark_img.at<uchar>(i, j) = min;
}
}
//std::cout <<"暗通道图像"<< dark_img;
imshow("暗通道图", dark_img);
return dark_img;
}

//暗通道的计算
Mat Producedarkimg(Mat& I, int windowsize)
{
int min = 255;
Mat dark_img(I.rows, I.cols, CV_8UC1); //单通道
int radius = (windowsize - 1) / 2;
int nr = I.rows; // number of rows
int nc = I.cols;
int b, g, r;
if (I.isContinuous()) { //当图像连通时，我们就可以把图像完全展开，看成是一行
nc = nr * nc;
nr = 1;
}
for (int i = 0;i<nr;i++)
{
const uchar* inData = I.ptr<uchar>(i);
uchar* outData = dark_img.ptr<uchar>(i);
for (int j = 0;j<nc;j++)
{
b = *inData++;
g = *inData++;
r = *inData++;
min = min>b ? b : min; //寻找三通道中的最小值
min = min>g ? g : min;
min = min>r ? r : min;
*outData++ = min;
min = 255;
}
}
dark_img = minfilter(dark_img, windowsize);
return dark_img;
}

//计算大气光值A
int* getatmospheric_light(Mat& darkimg, Mat& srcimg, int windowsize)
{
int radius = (windowsize - 1) / 2;
int nr = darkimg.rows, nl = darkimg.cols;
int darksize = nr*nl;
int topsize = darksize / 1000; //暗通道图中0.1%最亮像素的个数
int *A = new int[3];
int sum[3] = { 0,0,0 };
Pixel *toppixels, *allpixels;
toppixels = new Pixel[topsize];
allpixels = new Pixel[darksize];

for (int i = 0;i<nr;i++) {
const uchar* outData = darkimg.ptr<uchar>(i);
for (int j = 0;j<nl;j++)
{
allpixels[i*nl + j].value = *outData++;
allpixels[i*nl + j].x = i;
allpixels[i*nl + j].y = j;
}
}
//std::qsort(allpixels,darksize,sizeof(Pixel),qcmp);
std::sort(allpixels, allpixels + darksize, cmp);

memcpy(toppixels, allpixels, (topsize) * sizeof(Pixel)); //找到了在darkimg中最亮的0.1%个

int val0, val1, val2, avg, max = 0, maxi, maxj, x, y;
for (int i = 0;i<topsize;i++)
{
x = allpixels[i].x;y = allpixels[i].y; //暗通道中最亮像素的坐标
const uchar* outData = srcimg.ptr<uchar>(x);
outData += 3 * y;
val0 = *outData++;
val1 = *outData++;
val2 = *outData++;
avg = (val0 + val1 + val2) / 3;
if (max<avg) { max = avg;maxi = x;maxj = y; }
}
for (int i = 0;i<3;i++)
{
A[i] = srcimg.at<Vec3b>(maxi, maxj)[i];
//A[i]=srcimg.at<Vec4b>(maxi,maxj)[i];
//A[i]=A[i]>220?220:A[i];
}
return A; //返回A的3个b,g,r分量
}

//导向图滤波
Mat guidedFilter(Mat& transmission, Mat& graymat, Mat& trans, int windowsize, double ap) {
windowsize = 6 * windowsize;
Mat mean_P(transmission.rows, transmission.cols, CV_32FC1);
Mat mean_I(graymat.rows, graymat.cols, CV_32FC1);
Mat corr_I(graymat.rows, graymat.cols, CV_32FC1);
Mat corr_IP(graymat.rows, graymat.cols, CV_32FC1);
Mat var_I(graymat.rows, graymat.cols, CV_32FC1);
Mat cov_IP(graymat.rows, graymat.cols, CV_32FC1);
Mat a(graymat.rows, graymat.cols, CV_32FC1);
Mat b(graymat.rows, graymat.cols, CV_32FC1);
Mat mean_a(graymat.rows, graymat.cols, CV_32FC1);
Mat mean_b(graymat.rows, graymat.cols, CV_32FC1);
blur(transmission, mean_P, Size(windowsize,windowsize));
blur(graymat, mean_I, Size(windowsize, windowsize));
corr_I = graymat.mul(graymat);
corr_IP = graymat.mul(transmission);
var_I = corr_I - mean_I.mul(mean_I);
cov_IP = corr_IP - mean_I.mul(mean_P);
a = cov_IP / (var_I + ap);
b = mean_P - a.mul(mean_I);
blur(a, mean_a, Size(windowsize, windowsize));
blur(b, mean_b, Size(windowsize, windowsize));
trans = mean_a.mul(graymat) + mean_b;
return trans;
}

//计算透射图t并精细化
Mat getTransmission_dark(Mat& srcimg, Mat& darkimg, int *array, int windowsize)
{
float test;
float avg_A = (array[0] + array[1] + array[2]) / 3.0;
float w = 0.95;
int radius = (windowsize - 1) / 2;
int nr = srcimg.rows, nl = srcimg.cols;
Mat transmission(nr, nl, CV_32FC1);

for (int k = 0;k<nr;k++) {
const uchar* inData = darkimg.ptr<uchar>(k);
for (int l = 0;l<nl;l++)
{
transmission.at<float>(k, l) = 1 - w*(*inData++ / avg_A);
}
}
Mat trans(nr, nl, CV_32FC1);
Mat graymat(nr, nl, CV_8UC1);
Mat graymat_32F(nr, nl, CV_32FC1);
cvtColor(srcimg, graymat, CV_BGR2GRAY);
for (int i = 0;i<nr;i++) {
const uchar* inData = graymat.ptr<uchar>(i);
for (int j = 0;j<nl;j++)
graymat_32F.at<float>(i, j) = *inData++ / 255.0;
}
guidedFilter(transmission, graymat_32F, trans, 6 * windowsize, 0.001);

return trans;
}

//计算J（X）
Mat recover(Mat& srcimg, Mat& t, int *array, int windowsize)
{
int test;
int radius = (windowsize - 1) / 2;
int nr = srcimg.rows, nl = srcimg.cols;
float tnow = t.at<float>(radius, radius);
float t0 = 0.1;
Mat finalimg = Mat::zeros(nr, nl, CV_8UC3);
int val = 0;
for (int i = 0;i<3;i++) {
for (int k = radius;k<nr - radius;k++) {
const float* inData = t.ptr<float>(k); inData += radius;
const uchar* srcData = srcimg.ptr<uchar>(k); srcData += radius * 3 + i;
uchar* outData = finalimg.ptr<uchar>(k); outData += radius * 3 + i;
for (int l = radius;l<nl - radius;l++)
{
tnow = *inData++;
tnow = tnow>t0 ? tnow : t0;
val = (int)((*srcData - array[i]) / tnow + array[i]);
srcData += 3;
val = val<0 ? 0 : val;
*outData = val>255 ? 255 : val;
outData += 3;
}
}

}
return finalimg;
}

int main() {
int windowsize = 15;
Mat srcimg=imread("1.jpg");
imshow("原图", srcimg);
Mat darkimg(srcimg.rows,srcimg.cols, CV_8UC1);
Mat trans(srcimg.rows, srcimg.cols, CV_32FC1);
Mat pic(srcimg.rows, srcimg.cols, CV_8UC3);
int *A=new int[3];

//暗通道的计算
darkimg=Producedarkimg(srcimg, windowsize);
//计算大气光值A
A=getatmospheric_light(darkimg, srcimg,windowsize);
//计算透射图t并精细化
trans = getTransmission_dark(srcimg, darkimg, A, windowsize);
//计算J（X）
pic = recover(srcimg, trans, A, windowsize);
//显示去雾后的图像
imshow("去雾后",pic);
waitKey(0);
}五、效果图