图像处理（六）递归双边滤波磨皮

递归双边滤波是双边滤波的一种加速算法，加速比非常大，就像高斯模糊快速算法一样，加速起来，速度跟原算法相比，提高了十几倍。双边滤波的加速算法很多，文献都比较老旧，我这边主要讲一下比较新的算法：递归双边滤波，对应的paper为：《Recursive Bilateral Filtering》

这个算法比我另外一篇博文将的双指数滤波快一点，两篇文献的基本思想一样：



递归公式跟双指数的递归公式基本上一样。于是写代码就轻轻松松了，之前写过双指数的算法，接着只要在双指数的基础上，把代码改一改就OK了，这里就不详细讲解算法流程了，有兴趣的同学可以到paper的作者主页上，里面有代码，不过他的代码好像是要实现立体匹配的，把里面的相关的代码扣出来就可以了，下面是部分重要函数的代码。

[cpp] view
plain copy

void CBRF::recursive_bilateral_filter(double *data,unsigned char**disparity,double***in,unsigned char***texture,double sigma_spatial,double sigma_range,int height,int width,int nr_channel,double***temp,double***temp_2w)  

{  

    double yp[100];  

    double alpha=exp(-sqrt(2.0)/(sigma_spatial*width));//filter kernel size  

    double inv_alpha=(1-alpha);  

    double range_table[256];  

    for(int i=0;i<=255;i++)   

        range_table[i]=exp(-double(i)/(sigma_range*255));  

  

  

    BYTE* p0 =m_pImage;  

    const int nChannel = 4;  

    int m_length =m_nWidth*m_nHeight;  

//水平迭代滤波  

    //double***in_=in;/*horizontal filtering*/  

    //double***out_=temp;  

    double ***out=new double**[m_nHeight];  

    for (int i=0;i<m_nHeight;i++)  

    {  

        out[i]=new double*[m_nWidth];  

    }  

    for (int i=0;i<m_nHeight;i++)  

    {  

        for (int j=0;j<m_nWidth;j++)  

        {  

            out[i][j]=new double [3];  

        }  

    }  

    double ***out=new double**[m_nHeight];  

  

    for(int i=0;i<m_nHeight;i++)  

    {  

        BYTE* p0 =m_pImage+i*m_nWidth;  

        memcpy(out[i][0],p0,sizeof(double)*m_nWidth*nChannel);  

        for(int j=1;j<width;j++)   

        {  

            double weight=range_table[euro_dist_rgb_max(out[i][j],out[i][j-1])];  

            double alpha_=weight*alpha;  

            for(int k=0;k<nChannel;k++)   

            {  

                double ycc=inv_alpha*m_pImage[i*m_nWidth*nChannel+j*nChannel+k]+alpha_*yp[k];  

                out[i][j][k]=ycc  

                yp[k]=ycc;  

            }  

        }  

  

  

        int w1=width-1;  

        for(int c=0;c<nr_channel;c++)   

            out[y][w1][c]=0.5*(out[y][w1][c]+in[y][w1][c]);  

        memcpy(yp,out_[y][w1],sizeof(double)*nr_channel);  

        for(int x=width-2;x>=0;x--)   

        {  

            double weight=range_table[euro_dist_rgb_max(texture[y][x],texture[y][x+1])];  

            double alpha_=weight*alpha;  

            for(int c=0;c<nr_channel;c++)   

            {  

                double ycc=inv_alpha*in_[y][x][c]+alpha_*yp[c];  

                out_[y][x][c]=0.5*(out_[y][x][c]+ycc);  

                yp[c]=ycc;  

            }  

        }  

    }  

  

  

  

  

//垂直迭代滤波  

    /*in_=temp; 

    alpha=exp(-sqrt(2.0)/(sigma_spatial*height));//filter kernel size 

    inv_alpha=(1-alpha); 

    double**ycy,**ypy,**xcy,**xpy; 

    unsigned char**tcy,**tpy; 

    memcpy(out[0][0],in_[0][0],sizeof(double)*width*nr_channel); 

    for(int y=1;y<height;y++) 

    { 

        tpy=texture[y-1]; 

        tcy=texture[y]; 

        xcy=in_[y]; 

        ypy=out[y-1]; 

        ycy=out[y]; 

        for(int x=0;x<width;x++) 

        { 

            double weight=range_table[euro_dist_rgb_max(tcy[x],tpy[x])]; 

            double alpha_=weight*alpha; 

            for(int c=0;c<nr_channel;c++)  

                ycy[x][c]=inv_alpha*xcy[x][c]+alpha_*ypy[x][c]; 

        } 

    } 

    int h1=height-1; 

    ycy=temp_2w[0]; 

    ypy=temp_2w[1]; 

    memcpy(ypy[0],in_[h1][0],sizeof(double)*width*nr_channel); 

    for(int x=0;x<width;x++)  

    { 

        unsigned char disp=0;  

        double min_cost=0.5*(out[h1][x][0]+ypy[x][0]); 

        for(int c=1;c<nr_channel;c++) 

        { 

            double cost=0.5*(out[h1][x][c]+ypy[x][c]); 

            if(cost<min_cost) 

            { 

                min_cost=cost; 

                disp=c; 

            } 

        } 

        disparity[h1][x]=disp; 

    } 

    for(int y=h1-1;y>=0;y--) 

    { 

        tpy=texture[y+1]; 

        tcy=texture[y]; 

        xcy=in_[y]; 

        for(int x=0;x<width;x++) 

        { 

            double weight=range_table[euro_dist_rgb_max(tcy[x],tpy[x])]; 

            double alpha_=weight*alpha; 

            unsigned char disp=0; 

            double min_cost=255; 

            for(int c=0;c<nr_channel;c++)  

            { 

                ycy[x][c]=inv_alpha*xcy[x][c]+alpha_*ypy[x][c]; 

                double cost=0.5*(out[y][x][c]+ycy[x][c]); 

                if(cost<min_cost) 

                { 

                    min_cost=cost; 

                    disp=c; 

                } 

            } 

            disparity[y][x]=disp; 

        } 

        memcpy(ypy[0],ycy[0],sizeof(double)*width*nr_channel); 

    }*/  

}  

unsigned char CBRF::euro_dist_rgb_max(unsigned char *a,unsigned char *b)  

{  

    unsigned char x,y,z;  

    x=abs(a[0]-b[0]);  

    y=abs(a[1]-b[1]);  

    z=abs(a[2]-b[2]);  

    return(max(max(x,y),z));  

}  

看一些我用这个算法写的demo的测试结果：



原图



美图秀秀



递归双边滤波

最后人一定要靠自己，想要深入理解一个算法，让记忆更深，一定要自己把代码敲过一遍，既然选择了IT，就要坚持写代码。本文地址：http://blog.csdn.net/hjimce/article/details/45421207  
 作者：hjimce     联系qq：1393852684    更多资源请关注我的博客：http://blog.csdn.net/hjimce  原创文章，转载请保留这几行信息