Guided Image Filtering --引导图像滤波 解析+图像去雾应用
论文提出了一种新型图像滤波算法，基于局部线性模型，灵感来源于边缘保留滤波技术。该方法通过定义权重矩阵，结合导向图像和输入图像，实现边缘保留和去噪效果。论文详细解释了权重的定义、优化过程以及在图像去雾中的应用，并提供了基于OpenCV的代码实现。该方法在保持边缘细节的同时，显著提升了去噪性能，具有较高的实用价值。
论文地址以及代码：http://kaiminghe.com/eccv10/
Abstract: 本文中，作者阐述了提出了一种创新的图像滤波算法。该方法的灵感来源于局部线性模型。本研究引入了一个引导滤波器，该引导滤波器可以是输入图像本身，也可以是另一幅图像。算法通过利用引导滤波器的纹理特征对目标图像进行滤波处理。与传统的双边滤波器类似，该方法也是一种边缘保持平滑运算器，但其在边缘保留性能上优于双边滤波器。此外，该方法能够保留引导滤波器的结构特征，从而可用于去雾、羽化等图像处理任务。与现有技术相比，该算法在处理效率方面表现出色，堪称当前最有效的边缘保留滤波技术之一。
一、图像滤波的简要介绍
在前面部分，我们探讨了权重与边缘保留之间的关系，该文的阐述非常清晰，值得参考。
首先，我们定义一个线性转换的滤波器：
这个公式的意思是说，输出像素
的值是围绕像素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)仅作为粗略估计值，而所采用的软边缘方法显得过于复杂，且其计算复杂度同样较高。在本文中提出的方法导向图滤波，能够替代软边缘方法，以更精确地优化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);
}五、效果图
