导向图滤波与积分图算法

导向图滤波中需要用到一个滤波窗口W大小内图像的方差和均值的大量计算，而积分图可以快速的对图像内任何矩形所覆盖的像素值的和进行计算，从而使计算时间达到常数级别。

积分图

      对于积分图的概念，该篇文章解释的相当不错，我在这里只是做了一些总结和添加一些没有提到的东西。

在积分图像上任意位置(x, y)处的ii(x, y)表示该点左上角所有像素之和，表示如下：

从给定图像I从上到下、从左到右计算得到和的积分图像公式如下：

其中(x<0 || y<0) 时ii(x,y)=0, i(x,y)=0。其实，可以很明显的看出来，ii(x-1,y-1),是ii(x-1,y)和ii(x,y-1)重叠的部分。并且，当(x<0 || y<0) 时ii(x,y)=0, i(x,y)=0。

另外需要提一下的是，由以上定义可知，当所有的ii(x,y)都计算好以后组成的一张图是要比原图多一行和一列。在程序进行编程的时候还需要注意的是，一般的图像之是0-255，用int8类型即可保存，但是计算积分图以后，越往右下角其值越大。最好用int64的整型才能不会超出其数值范围。

定义已经解释清楚，那么任意矩形区域的像素值和的计算公式可以表达成：

导向图滤波：

        普通滤波，如均值滤波，或者高斯滤波等平滑滤波中，确定滤波后第i个像素点的值是根据以第i个像素点为中心的窗口w所覆盖的原图上的像素加权求均值而来，即： 

     其中，W,ij，是权重值。如果是均值滤波，则其值相同，都为1/n，n为窗口覆盖的像素个数。如果是高斯滤波器，则权重值W,ij服从二维高斯分布，即越往中心值越大。

    类似以上两种滤波方法，其实有很大的缺点就是都是各向同性滤波。即对图像中，无论是在纹理边缘区域还是在平缓区域，他们的权重值既然确定了就不在改变。从而导致对于图像边缘和平缓地方的处理方式都是一样的。但是，对于图像去噪声的滤波中，我们希望边缘（梯度变化较大）得到保持，平缓区域的噪声（梯度变化也很大）被滤除。仔细研究噪声和边缘的梯度就会发现，边缘的梯度变化在切线方向变化较小而在法线方向变化较大，噪声的梯度各个方向变化都很大。因此，形成了为了保持边缘梯度而去除噪声的梯度而形成的保边滤波算法，如双边滤波算法，维纳滤波算法等。导向图滤波算法也是包边滤波算法的一种。但是导向图滤波还可以用在其他图像处理方法中，如图像的去雾等

导向图滤波从名字上就能猜出，对一张图进行滤波时，还需要另一张图做引导。而且，导向图可以就是原图本身，这样导向滤波就变成一种保边滤波器，从而做到对图像的平滑降噪。

导向图滤波可以由上图所示，并且其依赖一个重要的假设：导向滤波器在导向图像I和滤波输出q之间在一个二维窗口内是一个局部线性模型，a和b是当窗口中心位于i时该线性函数的系数。由图中q = aI + b可知，导向图I是提供指导信息的，指导滤波器哪里是边缘（保留）哪里是平缓区域（磨平），指导方法就是动态改变a和b的值，如何改变或确定a和b的值就是其核心所在。

设引导图I，假设输入图像P，输出图像Q。导向滤波的目标是使得输入P和输出Q尽可能相同，同时纹理部分和引导图I相似。单个像素点用小写字母表示，所以需满足：

并且：  

通过积分可得：

考虑一个小窗口Wk，在Wk内认为a，b保持不变，设为ak，bk。Wk内的像素满足：

由上面这个公式可知，一个导向图上的像素点对应被滤波后的图像上的像素点。走到这一步，先暂停一下下面的推导。解释一下滤波后的一个像素到底是怎么来的。在一个窗口内，认定a，b不变，假设窗口的大小为3*3。当窗口对原图进行滑动滤波时，对于普通的滤波器，如高斯或者均值滤波器（全面已经介绍过），一个像素会被滤波器多次滤过，因此为了更加精确a和b的值，在得出A和B后，再次利用一次均值滤波确定ak和bk的值。滤过后的像素点的位置是对应窗口W中心的点对应的原图像所在的像素点。

继续推导，将上面三个公式综合可得：

 根据无约束图像复原方法，可对下面这个函数求最小值：

分别为a或b进行求导，并领求导后的函数为零，最后得到重要公式（一）：

其中，μk是导向图I中窗口Wk中的平均值，σ2k是导向图I中窗口Wk中的方差，|W|是窗口Wk中像素的数量，pk是待滤波图像p在窗口Wk中的均值，即。并且， 为导向图I和输入图P在Wk窗口下的协方差。从这个公式中可以看到，几乎所有的操作都与乘积和有关，并且乘积和的范围都是窗口Wk的大小。因此，文章刚开始介绍的积分图算法就派上大用场了。窗口Wk在确定的位置就会确定所覆盖像素点的线性函数，即：qi=a*Ii + b，每个中心像素可以看成是由多个线性函数共同描述，那么确定某一像素点输出值时，只需要将这些线性函数值求平均即可，重要公式（二）：

  

其中，

每一个像素点都是由a和b确定，而a和b的值在每一个像素都不相同，由此可以由a和b组成长宽都和原图相同的矩阵。

保边滤波性：

        由重要公式（二）发现，导向图I对滤波后的像素值起着很大的作用。如果导向图就是原图的话，怎么就会起到保边滤波性呢？。如果导向图是原图，那么重要公式（一）中，求解ak时就变成在Wk窗口范围下像素点的方差。那么就可以表示成：，在窗口Wk范围下，几个像素点的分布会出现两种情况：

一：高方差，在这种情况下如果，那么，那么输出值qi就是导向图的值，也就是原图像素点值，相当于没有起到滤波作用。

二：平坦，在这种情况下，，那么，那么qi等于Uk，即Wk窗口下这一块像素的均值，就变成了均值滤波。

这样是不是起到了有纹理的地方就保持原图，相对平坦区域就是滤波，起到了保边的作用。这一块到底是高方差还是平坦快，这是由ε来控制。

将该方法封装成一个类的话，Python的实现代码可参考如下：

    def _initFilter(self):
        I = self._I
        r = self._radius
        self._I_mean = cv2.blur(I, (r, r))
        I_mean_sq = cv2.blur(I ** 2, (r, r))
        self._I_var = I_mean_sq - self._I_mean ** 2

    def _computeCoefficients(self, p):
        r = self._radius
        p_mean = cv2.blur(p, (r, r))
        IP_cor = cv2.blur(p*self._I, (r, r))
        p_cov = IP_cor - self._I_mean * p_mean
        a = p_cov / (self._I_var + self._epsilon)
        b = p_mean - a * self._I_mean
        a_mean = cv2.blur(a, (r, r))
        b_mean = cv2.blur(b, (r, r))
        return a_mean, b_mean

    def _computeOutput(self, ab, I):
        a_mean, b_mean = ab
        return a_mean * I + b_mean

   def _initFilter(self):
        I = self._I
        r = self._radius
        eps = self._epsilon

        Ir, Ig, Ib = I[:, :, 0], I[:, :, 1], I[:, :, 2]

        self._Ir_mean = cv2.blur(Ir, (r, r))
        self._Ig_mean = cv2.blur(Ig, (r, r))
        self._Ib_mean = cv2.blur(Ib, (r, r))

        Irr_var = cv2.blur(Ir ** 2, (r, r)) - self._Ir_mean ** 2 + eps
        Irg_var = cv2.blur(Ir * Ig, (r, r)) - self._Ir_mean * self._Ig_mean
        Irb_var = cv2.blur(Ir * Ib, (r, r)) - self._Ir_mean * self._Ib_mean
        Igg_var = cv2.blur(Ig * Ig, (r, r)) - self._Ig_mean * self._Ig_mean + eps
        Igb_var = cv2.blur(Ig * Ib, (r, r)) - self._Ig_mean * self._Ib_mean
        Ibb_var = cv2.blur(Ib * Ib, (r, r)) - self._Ib_mean * self._Ib_mean + eps

        Irr_inv = Igg_var * Ibb_var - Igb_var * Igb_var
        Irg_inv = Igb_var * Irb_var - Irg_var * Ibb_var
        Irb_inv = Irg_var * Igb_var - Igg_var * Irb_var
        Igg_inv = Irr_var * Ibb_var - Irb_var * Irb_var
        Igb_inv = Irb_var * Irg_var - Irr_var * Igb_var
        Ibb_inv = Irr_var * Igg_var - Irg_var * Irg_var

        I_cov = Irr_inv * Irr_var + Irg_inv * Irg_var + Irb_inv * Irb_var
        Irr_inv /= I_cov
        Irg_inv /= I_cov
        Irb_inv /= I_cov
        Igg_inv /= I_cov
        Igb_inv /= I_cov
        Ibb_inv /= I_cov

        self._Irr_inv = Irr_inv
        self._Irg_inv = Irg_inv
        self._Irb_inv = Irb_inv
        self._Igg_inv = Igg_inv
        self._Igb_inv = Igb_inv
        self._Ibb_inv = Ibb_inv

    def _computeCoefficients(self, p):
        r = self._radius
        I = self._I
        Ir, Ig, Ib = I[:, :, 0], I[:, :, 1], I[:, :, 2]

        p_mean = cv2.blur(p, (r, r))

        Ipr_mean = cv2.blur(Ir * p, (r, r))
        Ipg_mean = cv2.blur(Ig * p, (r, r))
        Ipb_mean = cv2.blur(Ib * p, (r, r))

        Ipr_cov = Ipr_mean - self._Ir_mean * p_mean
        Ipg_cov = Ipg_mean - self._Ig_mean * p_mean
        Ipb_cov = Ipb_mean - self._Ib_mean * p_mean

        ar = self._Irr_inv * Ipr_cov + self._Irg_inv * Ipg_cov + self._Irb_inv * Ipb_cov
        ag = self._Irg_inv * Ipr_cov + self._Igg_inv * Ipg_cov + self._Igb_inv * Ipb_cov
        ab = self._Irb_inv * Ipr_cov + self._Igb_inv * Ipg_cov + self._Ibb_inv * Ipb_cov
        b = p_mean - ar * self._Ir_mean - ag * self._Ig_mean - ab * self._Ib_mean

        ar_mean = cv2.blur(ar, (r, r))
        ag_mean = cv2.blur(ag, (r, r))
        ab_mean = cv2.blur(ab, (r, r))
        b_mean = cv2.blur(b, (r, r))

        return ar_mean, ag_mean, ab_mean, b_mean

    def _computeOutput(self, ab, I):
        ar_mean, ag_mean, ab_mean, b_mean = ab

        Ir, Ig, Ib = I[:, :, 0], I[:, :, 1], I[:, :, 2]

        q = (ar_mean * Ir +
             ag_mean * Ig +
             ab_mean * Ib +
             b_mean)

        return q