这里主要是main代码。
实现的cpp在我资源里可以下载。
// day11.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include <iostream>;
#include "opencv2/core/core.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "guidedfilter.h"
#define ATV at<Vec3d>
#define ATD at<double>
#define ATU at<int>
#define win_size 5
#define omegar 0.95
#define gfR 15
#define gfE 0.001
using namespace cv;
using namespace std;
Mat pick_patch(Mat& mat, int r, int c, int patch_size)
{
 return mat(Range(r, r + patch_size), Range(c, c + patch_size));
}
double min3(double a, double b, double c)
{
 if (a > b)
 {
  if (b > c)
  {
   return c;
  }
  else
  {
   return b;
  }
 }
 else
 {
  if (a < c)
  {
   return a;
  }
  else
  {
   return c;
  }
 }
}
double calMin(Mat &img)
{
 double result;
 vector<Mat> temp;
 split(img, temp);
 Point minPoint;
 double minVal1 = 0;
 minMaxLoc(temp[0], &minVal1, NULL, &minPoint, NULL);
 double minVal2 = 0;
 minMaxLoc(temp[1], &minVal2, NULL, &minPoint, NULL);
 double minVal3 = 0;
 minMaxLoc(temp[2], &minVal3, NULL, &minPoint, NULL);
 result = min3(minVal1, minVal2, minVal3);
 return result;
}
Mat getDarkChannel(Mat &image, int winSize)
{
 int nRows = image.rows;
 int nCols = image.cols;
 Mat result = Mat::zeros(nRows, nCols, CV_64FC1);
 /*
 Rect roi(winSize / 2, winSize / 2, nCols, nRows);
 Mat tempI = Mat::ones(nRows + (winSize / 2) * 2, nCols + (winSize / 2) * 2, CV_64FC3);
 image.copyTo(tempI(roi));
 */
 Mat tempI;
 copyMakeBorder(image, tempI, winSize / 2, winSize / 2, winSize / 2, winSize / 2, BORDER_REPLICATE);
 double temp = 0.0;
 Mat patch;
 for (int i = 0; i < nRows; ++i)
 {
  for (int j = 0; j < nCols; ++j)
  {
   patch = pick_patch(tempI, i, j, winSize);
   temp = calMin(patch);
   result.ATD(i, j) = temp;
  }
 }
 return result;
}
vector<double> getAtmosphere(Mat &img, Mat &darkChannel)
{
 int m = img.rows;
 int n = img.cols;
 vector<double> result(3);
 int nPixels = m*n;
 int nSearchPixels = nPixels*0.01;
 Mat darkVect = darkChannel.reshape(0, nPixels);
 Mat imgVect = img.reshape(0, nPixels);
 Mat indices;
 sortIdx(darkVect, indices, SORT_EVERY_COLUMN + SORT_DESCENDING);
 vector<double> accumulator(3);
 accumulator[0] = 0.0;
 accumulator[1] = 0.0;
 accumulator[2] = 0.0;
 for (int i = 0; i < nSearchPixels; ++i)
 {
  accumulator[0] += imgVect.ATD(indices.ATU(i, 0), 0);
  accumulator[1] += imgVect.ATD(indices.ATU(i, 0), 1);
  accumulator[2] += imgVect.ATD(indices.ATU(i, 0), 2);
 }
 result[0] = accumulator[0] / nSearchPixels;
 result[1] = accumulator[1] / nSearchPixels;
 result[2] = accumulator[2] / nSearchPixels;
 return result;
}
Mat getTransmissionEstimate(Mat &img, vector<double> &atmosphere, double omega, int winSize)
{
 Mat result;
 int m = img.rows;
 int n = img.cols;
 Mat tempAtm(1, 1, CV_64FC3);
 tempAtm.ATV(0, 0)[0] = atmosphere[0];
 tempAtm.ATV(0, 0)[1] = atmosphere[1];
 tempAtm.ATV(0, 0)[2] = atmosphere[2];
 Mat repAtmosphere = repeat(tempAtm, m, n);
 Mat tempDiv;
 divide(img, repAtmosphere, tempDiv);
 result = 1.0 - omega * getDarkChannel(tempDiv, win_size);
 return result;
}
Mat guidedfilter(Mat& source, Mat& guidedImage, int radius, double epsilon)
{
 Mat result;
 Mat IP, I2;
 multiply(guidedImage, source, IP);
 multiply(guidedImage, guidedImage, I2);
 Mat meanI, meanP, meanIP, meanI2;
 Size wSize(2 * radius + 1, 2 * radius + 1);
 boxFilter(source, meanP, CV_64FC1, wSize);
 boxFilter(guidedImage, meanI, CV_64FC1, wSize);
 boxFilter(IP, meanIP, CV_64FC1, wSize);
 boxFilter(I2, meanI2, CV_64FC1, wSize);
 Mat varI = meanI2 - meanI.mul(meanI);
 Mat covIP = meanIP - meanI.mul(meanP);
 varI += epsilon;
 Mat a, b;
 divide(covIP, varI, a);
 b = meanP - a.mul(meanI);
 Mat meanA, meanB;
 boxFilter(a, meanA, CV_64FC1, wSize);
 boxFilter(b, meanB, CV_64FC1, wSize);
 result = meanA.mul(guidedImage) + meanB;
 return result;
}
Mat getRadiance(Mat &img, Mat &transmission, vector<double> &atmosphere)
{
 Mat result;
 int m = img.rows;
 int n = img.cols;
 Mat tempAtm(1, 1, CV_64FC3);
 tempAtm.ATV(0, 0)[0] = atmosphere[0];
 tempAtm.ATV(0, 0)[1] = atmosphere[1];
 tempAtm.ATV(0, 0)[2] = atmosphere[2];
 Mat repAtmosphere = repeat(tempAtm, m, n);
 Mat maxTransmission(m, n, CV_64FC3);
 for (int i = 0; i < m; ++i)
 {
  for (int j = 0; j < n; ++j)
  {
   if (transmission.ATD(i, j) <= 0.1)
   {
    transmission.ATD(i, j) == 0.1;
   }
   maxTransmission.ATV(i, j)[0] = transmission.ATD(i, j);
   maxTransmission.ATV(i, j)[1] = transmission.ATD(i, j);
   maxTransmission.ATV(i, j)[2] = transmission.ATD(i, j);
  }
 }
 Mat tempDiv = img - repAtmosphere;
 divide(tempDiv, maxTransmission, result);
 result = result + repAtmosphere;
 for (int i = 0; i < m; ++i)
 {
  for (int j = 0; j < n; ++j)
  {
   if (result.ATV(i, j)[0] < 0.0)
   {
    result.ATV(i, j)[0] = 0.0;
   }
   if (result.ATV(i, j)[0] > 1.0)
   {
    result.ATV(i, j)[0] = 1.0;
   }
   if (result.ATV(i, j)[1] < 0.0)
   {
    result.ATV(i, j)[1] = 0.0;
   }
   if (result.ATV(i, j)[1] > 1.0)
   {
    result.ATV(i, j)[1] = 1.0;
   }
   if (result.ATV(i, j)[2] < 0.0)
   {
    result.ATV(i, j)[2] = 0.0;
   }
   if (result.ATV(i, j)[2] > 1.0)
   {
    result.ATV(i, j)[2] = 1.0;
   }
  }
 }
 return result;
}
int main()
{
 Mat src = imread("009.bmp");// IMG_8766.jpg train.bmp tiananmen.png pumpkins.jpg IMG_8766.jpg
 int ratio = 1.0;
 resize(src, src, cv::Size(src.cols / ratio, src.rows / ratio));
 Mat tempImg;
 cvtColor(src, tempImg, CV_RGB2GRAY);
 int m = src.rows;
 int n = src.cols;
 src.convertTo(src, CV_64FC3, 1.0 / 255.0, 0);
 tempImg.convertTo(tempImg, CV_64FC1, 1.0 / 255.0, 0);
 Mat darkChannel;
 darkChannel = getDarkChannel(src, win_size);
 Mat temp;
 darkChannel.convertTo(temp, CV_8UC3, 255.0, 0);
// imshow("temp", temp);
 vector<double> atmosphere(3);
 atmosphere = getAtmosphere(src, darkChannel);
 cout << atmosphere[0] << " " << atmosphere[1] << " " << atmosphere[2] << endl;
 Mat transEst = getTransmissionEstimate(src, atmosphere, omegar, win_size);
 Mat x = guidedFilter(tempImg, transEst, gfR, gfE);
// imshow("x ", x);
 Mat result;
 result = getRadiance(src, x, atmosphere);
 result.convertTo(result, CV_8UC3, 255.0, 0);
 imwrite("myresult.jpg", result);
 imshow("result", result);
 waitKey();
 return 0;
}

导向滤波是处理速度较快的一种保边滤波器，适用范围也很广泛，去噪，细节增强以及HDR等地方都可以得到较好的结果，先从原理公式上进行推到，详细如下：
我在学习导向滤波物理意义的时候参考了下面这个网页：

引导滤波/导向滤波（Guided Filter）

我在学习导向滤波公式推导的时候参考了下面这个网页，讲得挺详细的，大家可以参考学习：

导向滤波(Guided Filter)公式详解
为了给自己留个纪念，我推导的过程如下：

前言
在图像处理上，导向滤波器（Guided Image Filter）是一种能使图像平滑化的非线性滤波器。与双边滤波器（Bilateral Filter）相同，这个滤波器同样能够在清楚保持图像边界的情况下，达到让图像平滑的效果。


但不同于双边滤波器，导向滤波器有两个优点：

首先，双边滤波器有非常大的计算复杂度（O(N^2)），但导向滤波器因为并未用到过于复杂的数学计算，有线性的计算复杂度。
双边滤波器因为数学模型的缘故，在某些时候会发生梯度反转（gradient reverse）的状况，出现图像有损；而导向滤波器因为在数学上以线性组合为基础出发，输出图片（Output Image）与引导图片（Guidance Image）的梯度方向一致，不会出现梯度反转的问题（大概率不出现，某些条件下必定不出现）。

可以说，导向滤波相比双边滤波的两大优势就是速度快和不会有梯度反转。
实际的应用场景除了去噪平滑外，还可以用于细节加强（detail smoothing/enhancement，如“羽化”）、HDR compression、image matting/feathering、haze removal（去雾）、joint upsampling、深度图修整等功能。


原理
为了达到图像平滑去噪效果，首先定义输出的结果图是输入图减去噪声后的结果。同时，为了让输出图保持引导图的边界，将输出图定为引导图的线性组合。
可以说，导向滤波核心原理是假设导向图I与滤波结果输出图q符合局部（以像素$k$为中心的$w_k$窗口内）线性模型：



局部线性模型（local linear model）保证了结果图与导向图的edge一致（$∇q = a∇I$）。
为了得到线性系数，需要构建方程求解。论文采用的是最小化输出q与输入图p之间的差异，即最小化窗口内的代价函数：



其中 $\epsilon$ 是防止 $a_k$ 过大的正则化参数。
方程的解可以根据 linear regression 求得，细节见参考资料[7]或[18]中推导：



其中，$\mu_k$和$\sigma^2_k$是导向图I在窗口$w_k$内的均值和方差，$|w|$是窗口$w_k$内的像素数目，$\overline{p}_k=\frac{1}{|w|}\sum_{i\in w_k}{p_i}$是窗口$w_k$内的均值。
基本上，根据得到的$a_k$和$b_k$就可以计算得出窗口$w_k$内的每一个$q_i$。但是进一步考虑，由于每一个像素不一定只被一个窗口$w_k$所包含，例如九宫格情况下中心像素点就被9个3x3的$w_k$窗口包含。



所以最简单的方式则是对这9个$w_k$窗口得到的$q_i$做一个加权平均，得到的最终$q_i$才是真正的结果值。



经过对所有$q_i$的加权平均（实际上用的是均值滤波），$∇q$不再是 $∇I$线性关系。但是由于$(\overline a_i, \overline b_i)$是经过均值滤波得到，在导向图的强边界处，输出图的梯度会比导向图小。这种情况下可以认为$∇q \approx \overline a∇I$，表示导向图I边界的强变化还能被输出图q维持。



算法伪码如下：



其中，$f_{mean}(·, r)$ 是半径为r的均值滤波器。
而方差和协方差定义如下：





对式子 (5) 进行变换，


则可以得到算法伪码中的：

特别说明：

通过参数$\epsilon$ 定义什么是“平坦区块(patch)”或“高变化区块”。若一个区块的方差远低于参数$\epsilon$ ，其通过滤波器后将被平滑；反之，方差远高于$\epsilon$的区块将被视为边界而被保留。
双边滤波中的范围方差(range variance)参数$\sigma _{r}^{2}$的功能和导向滤波的$\epsilon$相似。它们都定义了什么样的区块应该被平滑，而什么样的区块应该被保留。
实现
OpenCV中对导向滤波有CPU实现。


核心代码如下：
void GuidedFilterImpl::filter(InputArray src, OutputArray dst, int dDepth /*= -1*/)
{
    CV_Assert( !src.empty() && (src.depth() == CV_32F || src.depth() == CV_8U) );
    if (src.rows() != h || src.cols() != w)
    {
        CV_Error(Error::StsBadSize, "Size of filtering image must be equal to size of guide image");
        return;
    }

    if (dDepth == -1) dDepth = src.depth();
    int srcCnNum = src.channels();

    vector<Mat> srcCn(srcCnNum);
    vector<Mat>& srcCnMean = srcCn;
    split(src, srcCn);

    if (src.depth() != CV_32F)
    {
        parConvertToWorkType(srcCn, srcCn);
    }

    vector<vector<Mat> > covSrcGuide(srcCnNum);
    computeCovGuideAndSrc(srcCn, srcCnMean, covSrcGuide);

    vector<vector<Mat> > alpha(srcCnNum);
    for (int si = 0; si < srcCnNum; si++)
    {
        alpha[si].resize(gCnNum);
        for (int gi = 0; gi < gCnNum; gi++)
            alpha[si][gi].create(h, w, CV_32FC1);
    }
    runParBody(ComputeAlpha_ParBody(*this, alpha, covSrcGuide));
    covSrcGuide.clear();

    vector<Mat>& beta = srcCnMean;
    runParBody(ComputeBeta_ParBody(*this, alpha, srcCnMean, beta));

    parMeanFilter(beta, beta);
    parMeanFilter(alpha, alpha);

    runParBody(ApplyTransform_ParBody(*this, alpha, beta));
    if (dDepth != CV_32F)
    {
        for (int i = 0; i < srcCnNum; i++)
            beta[i].convertTo(beta[i], dDepth);
    }
    merge(beta, dst);
}


具体文件参考GitHub的OpenCV Contrib包实现。
GPU版导向滤波实现参考GitHub - TracelessLe/pybind11_guidedfilter_cuda。
cv::cuda::GpuMat GuidedFilterMono::filterSingleChannel(const cv::cuda::GpuMat &p, cv::cuda::Stream &stream) const {
  cv::cuda::GpuMat mean_p, mean_Ip, cov_Ip;
  box_filter->apply(p, mean_p, stream);
  cv::cuda::multiply(I, p, mean_Ip, 1, -1, stream);
  box_filter->apply(mean_Ip, mean_Ip, stream);
  cv::cuda::multiply(mean_I, mean_p, cov_Ip, 1, -1, stream);
  cv::cuda::subtract(mean_Ip,
                     cov_Ip,
                     cov_Ip,
                     cv::noArray(),
                     -1,
                     stream); // this is the covariance of (I, p) in each local patch.

  cv::cuda::GpuMat a, b;
  cv::cuda::add(var_I, cv::Scalar(eps), a, cv::noArray(), -1, stream);
  cv::cuda::divide(cov_Ip, a, a, 1, -1, stream); // Eqn. (5) in the paper;

  cv::cuda::multiply(a, mean_I, b, 1, -1, stream);
  cv::cuda::subtract(mean_p, b, b, cv::noArray(), -1, stream); // Eqn. (6) in the paper;

  box_filter->apply(a, a, stream);
  box_filter->apply(b, b, stream);

  cv::cuda::multiply(a, I, a, 1, -1, stream);
  cv::cuda::add(a, b, a, cv::noArray(), -1, stream);

  return a;
}

扩展讨论
（1）相比双边滤波，导向滤波有速度快和避免梯度反转等优势。



（2）基于原始的导向滤波算法引入resize得到的Fast导向滤波能够将时间复杂度从O(N)降到O(N / r^2)，同时保证滤波结果图像质量损失不大。其中r是resize（或称之为scale）的倍数。



参考资料
[1] wikipedia - Edge-preserving smoothing

[2] 维基百科 - 引导影像滤波器

[3] Guided Image Filtering - Kaiming He

[4] GitHub - opencv_contrib/modules/ximgproc/src/guided_filter.cpp

[5] OpenCV Docs - GuidedFilter

[6] 知乎 - 导向滤波原理（Guided Filter）

[7] 知乎 - 引导滤波guideFilter原理推导与实验

[8] 维基百科 - 方差

[9] 维基百科 - 协方差

[10]  知网 - 引导滤波算法的CUDA加速实现

[11] 豆丁网 - 引导滤波算法的CUDA加速实现

[12] GitHub - acstacey/GLFCV/src/guidedfilter.cpp

[13] GitHub - xxxzhou/oeip/oeip-win-cuda/GuidedFilterLayer.cpp

[14] cnblogs - CUDA加opencv复现导向滤波算法

[15] GitHub - foowaa/3DVisionUnit/GuidedFitlerOptimzation_CUDA/GuidedFilter.cu

[16] GitHub - TracelessLe/pybind11_guidedfilter_cuda

[17] csdn - 双边滤波原理浅析

[18] 导向滤波 Guided Image Filtering

        最近两天看了何凯明大神的导向滤波器（Guided Image Filtering），刚刚才顿悟整篇文章的核心。本文只针对于$I=p$的情况，即对图像$I$进行滤波处理。个人认为导向滤波器其实就是对均值滤波器的改进。由于均值滤波器会使图像模糊，尤其是对图像的边缘影响较大。然而，对图像进行处理时，图像的主要信息都隐藏在边缘处。所以为了保护图像的边缘信息，又能具有良好的滤波效果，是不是有一种滤波器在图像平滑的地方进行均值滤波，而在图像边缘的地方不进行滤波，或者进行轻微的滤波，那么图像的边缘信息不就保留下来了么，而且平滑的地方还有很好的均值滤波效果。

        那么问题来了，如何才能知道到底是不是图像的边缘呢？这里何凯明大神选取的指标是方差。方差是表示被统计信息变化特征的量。方差越大，表示图片在这里变化越大，也就表示此处是图像里某一物体的边缘，此处的信息需要进行保护。方差越小，表示图片在这里变化较小，该处所含信息量较少，可以进行滤波平滑处理。

        边缘参数找到了，如何将边缘参数的信息反映到对图像的滤波处理上呢？何凯明大神此处选取的是对图片每一像素值的线性处理，方程为：

$\begin{aligned}
q_i=a_iI_i+b_i && (1)
\end{aligned}$

        这里$q_i$是处理后的像素值，$I_i$是待处理的像素值。 $a_i$和$b_i$是对像素值 进行线性处理的参数。且

$\begin{aligned}
a_i=\frac{\sigma_i^2}{\sigma_i^2+\epsilon} && (2)
\end{aligned}$

$\begin{aligned}
b_i=\mu_i-a_i\mu_i && (3)
\end{aligned}$

        将（3）代入（1）得

$\begin{aligned}
q_i=a_i(I_i-\mu_i)+\mu_i && (4)
\end{aligned}$
        这里$\mu_i$为该点及附近点的均值（即：该点经过均值滤波器处理后的值）， $\sigma_i^2$为该点及附近点的方差。 $\epsilon$为一设置参数，用来表示图片此处平坦还是边缘的方差阈值。由式（2）可知，$a_i$的范围为（0,1）。当$\sigma \ll \epsilon$时，表示图片此处平坦，则$a_i \approx 0$，$q_i \approx \mu_i$ 。此处的滤波效果为经过均值滤波器的效果。当$\sigma_i \gg \epsilon$时，表示图片此处变化较为剧烈，为边缘，则$a_i \approx 1$, $q_i \approx I_i$结果为图片未经过任何处理。

细读论文会发现，滤波的最终公式为:

$\begin{aligned}
q_i = \bar a_iI_i+\bar b_i && (5)
\end{aligned}$

        分别对参数$a_i$和$b_i$进行了一次均值滤波处理。这里我没太想明白。可能是作者嫌处理后的图像不够平滑，又做了一次平滑处理吧。近似相当于对式（1）得到的结果进行了一次均值滤波处理，但是作者可能又觉得对式（1）得到的结果进行均值滤波处理边缘信息丢失的比较严重，所以作者只对参数$a_i$和$b_i$分别进行均值滤波处理。这样得到的图片会更加光滑一些，且边缘信息丢失的也不严重吧。