使输入P和输出Q尽只怕相似,导向滤波算法的兑现

本文从数学上演绎导向滤波的算法,其算法的现实性实未来下一篇导向滤波算法的达成介绍。

正文从数学上演绎导向滤波的算法,其算法的求实实以后下一篇导向滤波算法的兑现介绍。

由上篇导向滤波算法分析,根据(5)~(8)式就可以计算输出图像Q

设指引图G,输入图像P,输出图像Q。导向滤波的对象是驱动输入P和出口Q尽恐怕相同,同时纹理部分和指点图G相似。

设指点图G,输入图像P,输出图像Q。导向滤波的对象是驱动输入P和输出Q尽或者相同,同时纹理部分和指导图G相似。

图片 1  (5)

为了满足第1个对象,使输入P和出口Q尽或者相似,大家渴求最小化平方差

为了满意首个对象,使输入P和出口Q尽恐怕相似,大家渴求最小化平方差

图片 2  (6)

图片 3

图片 4

图片 5  (7)

为了满意第3个对象,使出口图像Q的纹理和引导图G的纹理相似,大家要

为了满意第四个对象,使出口图像Q的纹理和指点图G的纹路相似,大家要

图片 6  (8)

图片 7

图片 8

其中图片 9图片 10,/ai和/bi的结果要总结有所覆盖了像素i的窗口Wk的ak和bk的平均值。除了用平均值,在其实使用中,小编还看到过其余的总计/ai和/bi的章程。比如依据像素i在窗口Wk的职责,给予差别的权重。如若i距离窗口Wk的骨干地方越远,则给予ak和bk越低的权重。要是i位于窗口中央,则ak和bk有参天的权重。最常用的就是用高斯分布来予以不一致的权重。那样考虑到了隔壁像素距离远近对i影响的深浅,最终的结果会更纯粹一些。

积分得到

积分得到

那边自身要么用最简便的平均值的主意来总结/ai和/bi。大家日前早已假定了在窗口Wk内,ak和bk是常数,因此ak和bk只和Wk的地方有关。取Wk为半径为r的方形窗口,使窗口的主导像素地点遍历整个图像,那么就使Wk取到了差距的享有职位,在各样岗位统计出相应的ak和bk。所有的ak和bk结缘了和输入图像P相同长宽维度的多寡集合,记为A和B。对于自由像素i,/ai和/bi即分别为以i为主题半径r的窗口Wk内A和B的数目均值,这不正是大家耳熟能详的图像均值模糊的算法吗?而图像均值模糊有两种相比较早熟的火速算法,比如积分图算法和依照直方图的高速算法。只要有了A和B,就足以一本万利的使用均值模糊得出/ai和/bi,从而采用(8)计算出输出图像Q。

图片 11

图片 12

为了总括A和B,从(6)式看到,须要各自总结输入图像P和导向图G的均值模糊结果。而(5)式必要统计导向图G的方差,还有P和G的协方差。方差和协方差都事关到乘积的求和测算,可以由下边的公式,通过积分图来总结。那七个公式很不难推导出来,就不赘述了。

考虑3个小窗口Wk,在Wk内觉得a,b保持不变,设为ak,bk。Wk内的像素满意

设想贰个小窗口Wk,在Wk内觉得a,b保持不变,设为ak,bk。Wk内的像素满足

图片 13

图片 14  (1)

图片 15  (1)

图片 16

把(1)代入第三个目的,使窗口内的像素同时知足上边三个尺码。

把(1)代入第多少个目的,使窗口内的像素同时满意上边三个规范。

3个平方和,三个乘积和都得以用积分图来计量。只是要专注当图像丰盛大的时候,要用合适的数据类型。假若像素的数量范围是0~255的整型,假使平方积分图用三十五位的整型数据,那么只可以援救最大256×256尺寸的图像。超越这么些尺寸,就亟须要用6三个人的整型了。下边给出使用模板函数的乘积积分图函数,可以依照须求采用不一致的数据类型。p1和p2是图像数据指针,当它们对准相同的数据时,那么些函数就改为了平方积分图。注意积分图的长宽比原始数据的长宽都要大1。

图片 17  (2)

图片 18  (2)

图片 19图片 20

里面ε是贰个收拾大的ak的正则化参数。使(2)最小,满意

其间ε是三个收拾大的ak的正则化参数。使(2)最小,满意

/* Cumulative image of the multiplication of p1.*p2.
 p1 and p2 can be the same pointer and it becomes square cumulative image.
 The returned cumulative image MUST be freed by the caller! */
template <class T1, class T2, class T3>
BOOL MultiplyCumImage(T1 *p1, T2 *p2, T3 **ppCum)
{
    long i, j;
    long Width = GetWidth();
    long Height = GetHeight();
    long integral_len = (Width + 1)*(Height + 1);

    // Only allocate cumulative image memory when *ppCum is NULL
    if (*ppCum == NULL)
    {
        try { *ppCum = new T3[integral_len]; }
        catch (CException *pe)
        {
            pe->Delete();
            return FALSE;
        }
    }
    memset(*ppCum, 0, sizeof(T3)*integral_len);
    // The cumulative values of the leftmost and the topmost pixels are always 0.
    for (i = 1; i <= Height; i++)
    {
        T3 *prow, *puprow;
        prow = *ppCum + i*(Width + 1);
        puprow = *ppCum + (i - 1)*(Width + 1);
        T3 sum = 0;
        long up_row_idx = (i - 1)*Width;
        for (j = 1; j <= Width; j++)
        {
            long idx = up_row_idx + j - 1;
            sum += p1[idx] * p2[idx];
            prow[j] = puprow[j] + sum;
        }
    }
    return TRUE;
}

图片 21

图片 22

View
Code

其中|W|是窗口Wk的像素总数。解得

其中|W|是窗口Wk的像素总数。解得

 这样导向滤波完毕的重大难题都消除了,算法步骤如下:

图片 23  (3)

图片 24  (3)

  1. 总括指引图G的积分图和平方积分图
  2. 算算输入图像P的积分图, P和G的乘积积分图
  3. 用上两步得出的积分图计算P和G的均值,G的方差,P和G的协方差,窗口半径为r
  4. 接下来用(5)(6)式总括周到图A和B
  5. 总括A和B的积分图
  6. 计算A和B的窗口半径r的均值,并用(8)式总括输出图Q

图片 25  (4)

图片 26  (4)

 下边的参阅代码中,pData存储输入和输出图像,pGuidedData指点图,radius领域半径

如果设/pk是输入图P在窗口Wk的平均值,μk和σk2是指点图G在窗口Wk的平均值和方差。大家发现

如果设/pk是输入图P在窗口Wk的平均值,μk和σk2是指导图G在窗口Wk的平均值和方差。大家发现

图片 27图片 28

图片 29  (5)

图片 30  (5)

    long len = Width*Height;

    // Cululative image and square cululative for guided image G
    UINT64 *pCum = NULL, *pSquareCum = NULL;
    CumImage(pGuidedData, &pCum);
    MultiplyCumImage(pGuidedData, pGuidedData, &pSquareCum);

    // Allocate memory for a and b
    float *pa, *pb;
    pa = new float[len];
    pb = new float[len];
    memset(pa, 0, sizeof(float)*len);
    memset(pb, 0, sizeof(float)*len);

    UINT64 *pInputCum = NULL, *pGPCum = NULL;
    CumImage(pData, &pInputCum);
    MultiplyCumImage(pGuidedData, pData, &pGPCum);

    int field_size;
    UINT64 cum, square_cum;
    long uprow, downrow, upidx, downidx;        // In cumulative image
    long leftcol, rightcol;
    float g_mean, g_var;        // mean and variance of guided image
    long row_idx = 0;
    UINT64 p_cum, gp_cum;
    float p_mean;
    // Calculate a and b
    // Since we're going to calculate cumulative image of a and b, we have to calculate the whole image of a and b.
    for (i = 0; i < Height; i++)
    {
        // Check the boundary for radius
        if (i < radius) uprow = 0;
        else uprow = i - radius;
        upidx = uprow*(Width + 1);
        if (i + radius >= Height) downrow = Height;
        else downrow = i + radius + 1;
        downidx = downrow*(Width + 1);
        for (j = 0; j < Width; j++)
        {
            // Check the boundary for radius
            if (j < radius) leftcol = 0;
            else leftcol = j - radius;
            if (j + radius >= Width) rightcol = Width;
            else rightcol = j + radius + 1;
            field_size = (downrow - uprow)*(rightcol - leftcol);
            long p1, p2, p3, p4;
            p1 = downidx + rightcol;
            p2 = downidx + leftcol;
            p3 = upidx + rightcol;
            p4 = upidx + leftcol;
            // Guided image summary in the field
            cum = pCum[p1] - pCum[p2] - pCum[p3] + pCum[p4];
            // Guided image square summary in the field
            square_cum = pSquareCum[p1] - pSquareCum[p2] - pSquareCum[p3] + pSquareCum[p4];
            // Field mean
            g_mean = (float)(cum) / field_size;
            // Field variance
            g_var = float(square_cum) / field_size - g_mean * g_mean;
            // Summary of input image in the field
            p_cum = pInputCum[p1] - pInputCum[p2] - pInputCum[p3] + pInputCum[p4];
            // Input image field mean
            p_mean = float(p_cum) / field_size;
            // Multiply summary in the field
            gp_cum = pGPCum[p1] - pGPCum[p2] - pGPCum[p3] + pGPCum[p4];
            long idx = row_idx + j;
            pa[idx] = (float(gp_cum) / field_size - g_mean*p_mean) / (g_var + epsilon);
            pb[idx] = p_mean - g_mean*pa[idx];
        }
        row_idx += Width;
    }
    // not needed after this
    delete[] pCum;
    delete[] pSquareCum;
    delete[] pInputCum;
    delete[] pGPCum;

    // Cumulative image of a and b
    float *pCuma = NULL, *pCumb = NULL;
    CumImage(pa, &pCuma);
    CumImage(pb, &pCumb);

    // Finally calculate the output image q=ag+b
    float mean_a, mean_b;
    row_idx = Hstart*Width;
    for (i = Hstart; i < Hend; i++)
    {
        // Check the boundary for radius
        if (i < radius) uprow = 0;
        else uprow = i - radius;
        upidx = uprow*(Width + 1);
        if (i + radius >= Height) downrow = Height;
        else downrow = i + radius + 1;
        downidx = downrow*(Width + 1);
        for (j = Wstart; j < Wend; j++)
        {
            // Check the boundary for radius
            if (j < radius) leftcol = 0;
            else leftcol = j - radius;
            if (j + radius >= Width) rightcol = Width;
            else rightcol = j + radius + 1;
            field_size = (downrow - uprow)*(rightcol - leftcol);
            long p1, p2, p3, p4;
            p1 = downidx + rightcol;
            p2 = downidx + leftcol;
            p3 = upidx + rightcol;
            p4 = upidx + leftcol;
            // Field mean
            mean_a = (pCuma[p1] - pCuma[p2] - pCuma[p3] + pCuma[p4]) / field_size;
            // Field mean
            mean_b = (pCumb[p1] - pCumb[p2] - pCumb[p3] + pCumb[p4]) / field_size;
            // New pixel value
            long idx = row_idx + j;
            int value = int(mean_a*pGuidedData[idx] + mean_b);
            CLAMP0255(value);
            pData[idx] = value;
        }
        row_idx += Width;
    }

    delete[] pa;
    delete[] pb;
    delete[] pCuma;
    delete[] pCumb;

图片 31  (6)

图片 32  (6)

View
Code

其中图片 33是率领图G和输入图P在Wk的协方差。

其中图片 34是指引图G和输入图P在Wk的协方差。

 导向滤波还有一种高效算法,基本思想是经过下采样输入图P和指引图G,拿到较小的图像P’和G’。用它们来计算周到A’和B’。然后通过线性插值的法门復苏原来大小得到近似的A和B,用来和原来大小的教导图来计量输出Q。那样,ak和bk的一个钱打二拾肆个结不是在全尺寸图像上,能节省多如牛毛运算量,而结尾的结果不受很大影响。

计算出ak,bk后,就可以依据(1)来计量窗口Wk的输出像素。对于两个像素i,输出值qi和拥有覆盖像素i的窗口Wk有关。所以当Wk不同,qi的值也不均等。贰个简短的政策是平均全数或许的qi值。统计了拥有覆盖i的窗口Wk的ak,bk,全体覆盖像素i的窗口Wk的个数为|W|,那么

计算出ak,bk后,就足以依据(1)来计量窗口Wk的出口像素。对于一个像素i,输出值qi和持有覆盖像素i的窗口Wk有关。所以当Wk不同,qi的值也不雷同。二个简单的国策是平均全体或然的qi值。总结了拥有覆盖i的窗口Wk的ak,bk,全部覆盖像素i的窗口Wk的个数为|W|,那么

 上边以带有保边平滑性格的导向滤波为例,来探望效果。如上一篇所说,当输入图P和指导图G相同时,导向滤波突显保边平滑本性。

图片 35  (7)

图片 36  (7)

图片 37  
图片 38

图片 39  (8)

图片 40  (8)

                              
原图                                                                                                
半径3,ε=52

其中图片 41图片 42

其中图片 43图片 44

图片 45  
图片 46



                   
半径3,ε=102                                                
                                          
半径3,ε=202

尤其的,当率领图G和输入图像P相同的时候,导向滤波出现边缘保持平滑性格,分析如下。当G=P时,很明显图片 47图片 48,由式(5)(6)得到图片 49图片 50。当ε=0时,ak=1,bk=0,即出口和输入图像相同。即使ε>0,考虑二种情形。

专门的,当辅导图G和输入图像P相同的时候,导向滤波出现边缘保持平滑本性,分析如下。当G=P时,很明朗图片 51图片 52,由式(5)(6)得到图片 53图片 54。当ε=0时,ak=1,bk=0,即出口和输入图像相同。假诺ε>0,考虑二种处境。

图片 55
  图片 56

第叁种,高方差。如若图像P在窗口Wk中有许多转变,那么σk2>>ε,有ak≈1,bk≈0。

先是种,高方差。倘诺图像P在窗口Wk中有这一个变迁,那么σk2>>ε,有ak≈1,bk≈0。

                         原图            
                                         
 半径5,ε=202

第二种,平坦块。那么σk2<<ε,有ak≈0,bk≈μk。要是一切输入图像都如窗口Wk一般很平整的,当ak,bk被平均后收获/ak≈0,/bk≈μk,qi≈μk

第二种,平坦块。那么σk2<<ε,有ak≈0,bk≈μk。假使整个输入图像都如窗口Wk诚如很平整的,当ak,bk被平均后收获/ak≈0,/bk≈μk,qi≈μk

 总体来说,P=G时,导向滤波的保边平滑天性和带有保边功效领域平滑滤波有相近的效能。

这么,当三个像素在高方差的窗口中时,它的输出值是不变的。在平坦区域中,它的输出值变成周围窗口像素的平均值。具体的,高方差和平坦的正规化是由参数ε控制的。即使窗口的方差比此参数小的多则被平整,那么方差大得多的则被封存。而窗口的轻重缓急决定了是参照周围多大范围的像向来总结方差和均值。

诸如此类,当1个像素在高方差的窗口中时,它的输出值是不变的。在平坦区域中,它的输出值变成周围窗口像素的平均值。具体的,高方差和平坦的正式是由参数ε控制的。尽管窗口的方差比此参数小的多则被平整,那么方差大得多的则被封存。而窗口的高低决定了是参照周围多大范围的像一向测算方差和均值。

 



时至前天,已经可以依照(5)~(8)式总结导向滤波的参数,从而总结输出图像Q。导向滤波的的算法实以往下一篇作品介绍。下边把导向滤波用通用滤波器核来表明举办特别分析。

距今,已经足以依照(5)~(8)式总括导向滤波的参数,从而总括输出图像Q。导向滤波的的算法实今后下一篇小说介绍。上边把导向滤波用通用滤波器核来表明举行更为分析。

导向滤波在像素点i的滤波结果可以发布为贰个加权平均

导向滤波在像素点i的滤波结果可以发挥为3个加权平均

图片 57  (9)

图片 58  (9)

内部i,j都以像素下标。滤波器核Wij是指点图G的函数并且与P独立。下边总括滤波器核。把(6)带入(8)消去b得到:

其中i,j都是像素下标。滤波器核Wij是指导图G的函数并且与P独立。上边计算滤波器核。把(6)带入(8)消去b拿到:

图片 59  (10)

图片 60  (10)

求偏导

求偏导

图片 61  (11)

图片 62  (11)

其中

其中

图片 63  (12)

图片 64  (12)

图片 65当j处于窗口Wk时,否则为0。

图片 66当j处于窗口Wk时,否则为0。

图片 67

图片 68

图片 69  (13)

图片 70  (13)

把(12)(13)带入(11)得到

把(12)(13)带入(11)得到

图片 71

图片 72

图片 73  (14)

图片 74  (14)

图片 75

图片 76

考虑富含边缘阶跃信号的图像,假如i,j在边缘的同一侧,(gik)和(gjk)符号相同,而当它们在边缘的两侧时则符号差异。所以Wij(G)当三个像素点在边缘的两侧时比在同等边时要小得多。那就印证隔着边缘的时候,pj对结果的孝敬很小,窗口像素不会平均到联合。而当j和i在边缘的同侧时,输出像素qi是同侧像素的加权平均值。加权周到由导向图G决定。由此得以见见,导向滤波确实可以起到保留边缘的成效。

设想富含边缘阶跃信号的图像,假若i,j在边缘的同一侧,(gik)和(gjk)符号相同,而当它们在边缘的两侧时则符号不一样。所以Wij(G)当八个像素点在边缘的两侧时比在相同边时要小得多。那就认证隔着边缘的时候,pj对结果的贡献很小,窗口像素不会平均到三头。而当j和i在边缘的同侧时,输出像素qi是同侧像素的加权平均值。加权全面由导向图G决定。因而可以看看,导向滤波确实可以起到保留边缘的意义。