导向滤波算法的兑现

图片 1  (6)

图片 2
  图片 3

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

  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

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

先是种,高方差。若是图像P在窗口Wk中有诸多变迁,那么σk2>>ε,有ak≈1,bk≈0。

                         原图            
                                         
 半径5,ε=202

                              
原图                                                                                                
半径3,ε=52

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

图片 5图片 6

/* 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;
}

图片 7  (5)

图片 8  (5)

其中图片 9图片 10,/ai和/bi的结果要总括有所覆盖了像素i的窗口Wk的ak和bk的平均值。除了用平均值,在实际应用中,作者还看到过任何的总括/ai和/bi的方法。比如依据像素i在窗口Wk的职位,给予差其他权重。如若i距离窗口Wk的焦点岗位越远,则给予ak和bk越低的权重。假诺i位于窗口宗旨,则ak和bk有最高的权重。最常用的就是用高斯分布来予以差异的权重。那样考虑到了邻座像素距离远近对i影响的分寸,最终的结果会更准确一些。

把(1)代入第二个对象,使窗口内的像素同时知足下边五个标准化。

图片 11  
图片 12

图片 13  
图片 14

图片 15  (10)

图片 16  (7)

 总体来说,P=G时,导向滤波的保边平滑性情和带有保边效指点域平滑滤波有像样的功能。

图片 17

 那样导向滤波完成的首要难题都化解了,算法步骤如下:

 那样导向滤波完毕的要害难题都化解了,算法步骤如下:

图片 18  (12)

View
Code

                         原图            
                                         
 半径5,ε=202

本文从数学上演绎导向滤波的算法,其算法的求实实今后下一篇导向滤波算法的落成介绍。

图片 19
  图片 20

图片 21图片 22

图片 23  (8)

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

View
Code


图片 24  
图片 25

图片 26

其中图片 27图片 28

图片 29  (8)

图片 30  (7)

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

八个平方和,四个乘积和都得以用积分图来测算。只是要注意当图像充分大的时候,要用合适的数据类型。假如像素的数目范围是0~255的整型,假使平方积分图用3四人的整型数据,那么只好帮助最大256×256高低的图像。当先那些尺寸,就非得要用六十几位的整型了。上面给出使用模板函数的乘积积分图函数,可以按照需求选用差其余数据类型。p1和p2是图像数据指针,当它们对准相同的数码时,那几个函数就改为了平方积分图。注意积分图的长宽比原始数据的长宽都要大1。

图片 31

图片 32  (13)

图片 33

图片 34  (8)

中间ε是多少个惩治大的ak的正则化参数。使(2)最小,满意

 下边以饱含保边平滑天性的导向滤波为例,来看看效果。如上一篇所说,当输入图P和指引图G相同时,导向滤波显示保边平滑天性。

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

图片 35

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

    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;

积分拿到

/* 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;
}

 上边的参考代码中,pData存储输入和出口图像,pGuidedData指引图,radius领域半径

View
Code

 

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

 总体来说,P=G时,导向滤波的保边平滑个性和包罗保边成效领域平滑滤波有接近的作用。

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

求偏导

图片 36图片 37

  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

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

                              
原图                                                                                                
半径3,ε=52

View
Code

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

 

此处小编要么用最简便易行的平均值的艺术来总结/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。


此处本人依旧用最简便易行的平均值的法门来计算/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。

图片 38图片 39

图片 40  (1)

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

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

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

 导向滤波还有一种高效算法,基本思想是因此下采样输入图P和率领图G,得到较小的图像P’和G’。用它们来总括周到A’和B’。然后经过线性插值的不二法门復苏原来大小得到近似的A和B,用来和原始大小的指引图来总结输出Q。那样,ak和bk的盘算不是在全尺寸图像上,能节省恒河沙数运算量,而最终的结果不受相当的大影响。

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

图片 41  (14)

图片 42

图片 43  
图片 44

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

 上边的参照代码中,pData存储输入和输出图像,pGuidedData指引图,radius领域半径

图片 45  (6)

图片 46  (11)

    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;

图片 47  (5)

图片 48

其中图片 49图片 50,/ai和/bi的结果要计算有所覆盖了像素i的窗口Wk的ak和bk的平均值。除了用平均值,在实质上选用中,我还察看过其余的统计/ai和/bi的点子。比如按照像素i在窗口Wk的岗位,给予不一致的权重。若是i距离窗口Wk的为主岗位越远,则给予ak和bk越低的权重。如若i位于窗口中央,则ak和bk有最高的权重。最常用的就是用高斯分布来予以不同的权重。那样考虑到了紧邻像素距离远近对i影响的轻重缓急,最终的结果会更可信一些。

贰个平方和,贰个乘积和都得以用积分图来计量。只是要留心当图像充裕大的时候,要用合适的数据类型。假如像素的多寡范围是0~255的整型,若是平方积分图用31人的整型数据,那么只好辅助最大256×256大大小小的图像。超过这些尺寸,就不能够不要用6四个人的整型了。上面给出使用模板函数的乘积积分图函数,可以按照须要采取区其他数据类型。p1和p2是图像数据指针,当它们对准相同的数额时,那么些函数就改成了平方积分图。注意积分图的长宽比原始数据的长宽都要大1。

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

设想富含边缘阶跃信号的图像,尽管i,j在边缘的同一侧,(gik)和(gjk)符号相同,而当它们在边缘的两侧时则符号差异。所以Wij(G)当多少个像素点在边缘的两侧时比在同样边时要小得多。那就表达隔着边缘的时候,pj对结果的孝敬一点都不大,窗口像素不会平均到一起。而当j和i在边缘的同侧时,输出像素qi是同侧像素的加权平均值。加权周到由导向图G决定。因而可以看出,导向滤波确实可以起到保留边缘的意义。

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

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

图片 52

图片 53

计算出ak,bk后,就足以依据(1)来计量窗口Wk的出口像素。对于3个像素i,输出值qi和拥有覆盖像素i的窗口Wk有关。所以当Wk不同,qi的值也不等同。一个粗略的策略是平均全数大概的qi值。计算了富有覆盖i的窗口Wk的ak,bk,全部覆盖像素i的窗口Wk的个数为|W|,那么

图片 54  (2)

图片 55

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

特意的,当辅导图G和输入图像P相同的时候,导向滤波出现边缘保持平滑性情,分析如下。当G=P时,很强烈图片 56图片 57,由式(5)(6)得到图片 58图片 59。当ε=0时,ak=1,bk=0,即出口和输入图像相同。若是ε>0,考虑二种处境。

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

图片 60  (4)

图片 61  (9)

图片 62  (3)

图片 63  (6)

那般,当一个像素在高方差的窗口中时,它的输出值是不变的。在平坦区域中,它的输出值变成周围窗口像素的平均值。具体的,高方差和平坦的正规化是由参数ε控制的。假设窗口的方差比此参数小的多则被平整,那么方差大得多的则被保留。而窗口的大大小小决定了是参照周围多大范围的像一向计量方差和均值。

其中

图片 64  (7)

图片 65

相关文章