你好,游客 登录
背景:
阅读新闻

图像保边滤波器---双边滤波算法

[日期:2018-05-25] 来源:  作者: [字体: ]

双边滤波算法

  1. #include "string.h" 
  2. #include "stdio.h" 
  3. #include "stdlib.h" 
  4. #include "math.h" 
  5. #include"SoftSkin.h" 
  6.  
  7.  
  8.  
  9. //垂直方向递归 
  10. void runVerticalHorizontal(double *data,int width,int height,double spatialDecay,double *exp_table,double *g_table) 
  11.     int length0=height*width; 
  12.     double* g= new double[length0]; 
  13.     int m = 0; 
  14.     for (int k2 = 0;k2<height;++k2) 
  15.     { 
  16.         int n = k2; 
  17.         for (int k1 = 0;k1<width;++k1) 
  18.         { 
  19.             g[n]=data[m++]; 
  20.             n += height; 
  21.         } 
  22.     } 
  23.     double*p = new double[length0]; 
  24.     double*r = new double[length0]; 
  25.     memcpy(p, g, sizeof(double) * length0); 
  26.     memcpy(r, g, sizeof(double) * length0); 
  27.     for (int k1 = 0;k1<width; ++k1) 
  28.     { 
  29.         int startIndex=k1 * height; 
  30.         double mu = 0.0; 
  31.         for (int k=startIndex+1,K =startIndex+height;k<K;++k) 
  32.         { 
  33.             int div0=fabs(p[k]-p[k-1]); 
  34.             mu =exp_table[div0]; 
  35.             p[k] = p[k - 1] * mu + p[k] * (1.0 - mu);//文献中的公式1,这里做了一下修改,效果影响不大 
  36.         } 
  37.         for (int k =startIndex+height-2;startIndex <= k;--k) 
  38.         { 
  39.             int div0=fabs(r[k]-r[k+1]); 
  40.             mu =exp_table[div0]; 
  41.             r[k] = r[k+1] * mu + r[k] * (1.0-mu) ;//文献公式3 
  42.         } 
  43.     } 
  44.     double rho0=1.0/(2-spatialDecay); 
  45.     for (int k = 0;k <length0;++k) 
  46.     { 
  47.         r[k]= (r[k]+p[k])*rho0-g_table[(int)g[k]]; 
  48.     } 
  49.     m = 0; 
  50.     for (int k1=0;k1<width;++k1) 
  51.     { 
  52.         int n = k1; 
  53.         for (int k2 =0;k2<height;++k2) 
  54.         { 
  55.             data[n] = r[m++]; 
  56.             n += width; 
  57.         } 
  58.     } 
  59.     memcpy(p,data, sizeof(double) * length0); 
  60.     memcpy(r,data, sizeof(double) * length0); 
  61.     for (int k2 = 0; k2<height;++k2) 
  62.     { 
  63.  
  64.         int startIndex=k2 * width; 
  65.         double mu = 0.0; 
  66.         for (int k=startIndex+1, K=startIndex+width;k<K;++k) 
  67.         { 
  68.             int div0=fabs(p[k]-p[k-1]); 
  69.             mu =exp_table[div0]; 
  70.             p[k] = p[k - 1] * mu + p[k] * (1.0 - mu); 
  71.         } 
  72.         for (int k=startIndex+width-2;startIndex<=k;--k) 
  73.         { 
  74.             int div0=fabs(r[k]-r[k+1]); 
  75.             mu =exp_table[div0]; 
  76.             r[k] = r[k + 1] * mu + r[k] * (1.0 - mu) ; 
  77.         } 
  78.     } 
  79.  
  80.     double init_gain_mu=spatialDecay/(2-spatialDecay); 
  81.     for (int k = 0; k <length0; k++) 
  82.     { 
  83.         data[k]=(p[k]+r[k])*rho0-data[k]*init_gain_mu;//文献中的公式5 
  84.     } 
  85.     delete p; 
  86.     delete r; 
  87.     delete g; 
  88. //水平方向递归 
  89. void runHorizontalVertical(double *data,int width,int height,double spatialDecay,double *exptable,double *g_table) 
  90.     int length=width*height; 
  91.     double* g = new double[length]; 
  92.     double* p = new double[length]; 
  93.     double* r = new double[length]; 
  94.     memcpy(p,data, sizeof(double) * length); 
  95.     memcpy(r,data, sizeof(double) * length); 
  96.     double rho0=1.0/(2-spatialDecay); 
  97.     for (int k2 = 0;k2 < height;++k2) 
  98.     { 
  99.         int startIndex=k2 * width; 
  100.         for (int k=startIndex+1,K=startIndex+width;k<K;++k) 
  101.         { 
  102.             int div0=fabs(p[k]-p[k-1]); 
  103.             double mu =exptable[div0]; 
  104.             p[k] = p[k - 1] * mu + p[k] * (1.0 - mu);//文献公式1 
  105.  
  106.         } 
  107.         for (int k =startIndex + width - 2;startIndex <= k;--k) 
  108.         { 
  109.             int div0=fabs(r[k]-r[k+1]); 
  110.             double mu =exptable[div0]; 
  111.             r[k] = r[k + 1] * mu + r[k] * (1.0 - mu);//文献公式3 
  112.         } 
  113.         for (int k =startIndex,K=startIndex+width;k<K;k++) 
  114.         { 
  115.             r[k]=(r[k]+p[k])*rho0- g_table[(int)data[k]]; 
  116.         } 
  117.     } 
  118.      
  119.     int m = 0; 
  120.     for (int k2=0;k2<height;k2++) 
  121.     { 
  122.         int n = k2; 
  123.         for (int k1=0;k1<width;k1++) 
  124.         { 
  125.             g[n] = r[m++]; 
  126.             n += height; 
  127.         } 
  128.     } 
  129.     memcpy(p, g, sizeof(double) * height * width); 
  130.     memcpy(r, g, sizeof(double) * height * width); 
  131.     for (int k1=0;k1<width;++k1) 
  132.     { 
  133.         int startIndex=k1 * height; 
  134.         double mu = 0.0; 
  135.         for (int k =startIndex+1,K =startIndex+height;k<K;++k) 
  136.         { 
  137.             int div0=fabs(p[k]-p[k-1]); 
  138.             mu =exptable[div0]; 
  139.             p[k] = p[k - 1] * mu + p[k] * (1.0 - mu); 
  140.         } 
  141.         for (int k=startIndex+height-2;startIndex<=k;--k) 
  142.         { 
  143.             int div0=fabs(r[k]-r[k+1]); 
  144.             mu =exptable[div0]; 
  145.             r[k] = r[k + 1] * mu + r[k] * (1.0 - mu); 
  146.         } 
  147.     } 
  148.  
  149.     double init_gain_mu=spatialDecay/(2-spatialDecay); 
  150.     for (int k = 0;k <length;++k) 
  151.     { 
  152.         r[k]= (r[k]+p[k])*rho0- g[k]*init_gain_mu; 
  153.     } 
  154.     m = 0; 
  155.     for (int k1=0;k1<width;++k1) 
  156.     { 
  157.         int n = k1; 
  158.         for (int k2=0;k2<height;++k2) 
  159.         { 
  160.             data[n]=r[m++]; 
  161.             n += width; 
  162.         } 
  163.     } 
  164.     delete p; 
  165.     delete r; 
  166.     delete g; 
  167.  
  168. void ApplyBiExponentialEdgePreservingSmoother(double photometricStandardDeviation, double spatialDecay, unsigned char* m_pImage, int  m_nWidth, int m_nHeight,int m_nStride) 
  169.     double m_exp_table[256]; 
  170.     double m_g_table[256]; 
  171.     memset(m_exp_table, 0, 256); 
  172.     memset(m_g_table, 0, 256); 
  173.     double c=-0.5/(photometricStandardDeviation * photometricStandardDeviation); 
  174.     double mu=spatialDecay/(2-spatialDecay); 
  175.     for (int i=0;i<=255;i++) 
  176.     { 
  177.         float a=exp(c*i*i); 
  178.         m_exp_table[i]=(1-spatialDecay)* exp(c*i*i); 
  179.         m_g_table[i]=mu*i; 
  180.     } 
  181.     unsigned char* p0 =m_pImage; 
  182.     const int nChannel = 4; 
  183.     int m_length =m_nWidth*m_nHeight; 
  184.     float maxerror=0; 
  185.     float sum=0; 
  186.     // 对每个channel进行处理 
  187.  
  188.     for (int idxChannel=0;idxChannel <nChannel; idxChannel++) 
  189.     { 
  190.         double *data1 = new double[m_length]; 
  191.         double* data2 = new double[m_length]; 
  192.  
  193.         unsigned char *p1=p0+idxChannel; 
  194.         for (int i = 0; i < m_length;++i) 
  195.         { 
  196.             data1[i] = p1[i * nChannel]; 
  197.         } 
  198.         memcpy(data2,data1, sizeof(double) * m_length); 
  199.  
  200.         runHorizontalVertical(data1, m_nWidth, m_nHeight,spatialDecay,m_exp_table,m_g_table); 
  201.         runVerticalHorizontal(data2, m_nWidth, m_nHeight,spatialDecay,m_exp_table,m_g_table); 
  202.         sum=0; 
  203.         for (int i =0;i<m_length;++i) 
  204.         { 
  205.             double val=(data1[i] + data2[i]) * 0.5; 
  206.  
  207.             if(255.0<val)val=255.0; 
  208.             p1[i * nChannel]=(unsigned char)val; 
  209.          } 
  210.         delete data1; 
  211.         delete data2; 
  212.     } 
  213.      
  214.  
  215. void f_Bilateralfilter(unsigned char* pImage, int nWidth, int nHeight, int nStride, double std) 
  216.     if (pImage == NULL || std == 0) 
  217.     { 
  218.         return
  219.     } 
  220.     ApplyBiExponentialEdgePreservingSmoother(std,std * 0.001,pImage, nWidth, nHeight, nStride); 

效果如下:

原图

双边滤波效果图(Delta=20)

递归双边滤波是一种快速双边滤波算法,这里本人是对四通道处理的,大家可以在YUV颜色空间对Y通道处理,UV不变,这样速度会更快一些。

推荐 打印 | 录入:Cstor | 阅读:
相关新闻      
本文评论   
评论声明
  • 尊重网上道德,遵守中华人民共和国的各项有关法律法规
  • 承担一切因您的行为而直接或间接导致的民事或刑事法律责任
  • 本站管理人员有权保留或删除其管辖留言中的任意内容
  • 本站有权在网站内转载或引用您的评论
  • 参与本评论即表明您已经阅读并接受上述条款