磨皮对于现在的图像处理软件，可以说是一项重要的功能，在天天P图，可牛，ps，美图秀秀等软件中随处可见，有可能即使你非常熟悉图像处理的算法，然而却不懂磨皮怎么实现。其实磨皮就是所谓的保边缘滤波，也就是说在图像处理领域只要是滤波算法都可以实现磨皮，只是效果好坏的区别，然而现在对于大部分，都要求具有保细节的功能，这边先给大家介绍一种算法：双指数保边缘滤波，对应的 外围文献为：《Bi-Exponential Edge-Preserving Smoother》
下面是这篇文献最重要的部分，算法整个过程分为三个步骤：
（1）水平方向递归：包括对原始图像的水平方向进行向前递归（公式1）、对原始图片的水平方向进行向后递归（公式3），然后把向前递归结果、向后递归结果、原图像数据做一个加权组合(公式5)，得到新的像素值。
（2）垂直方向递归：同样的对原始图像数据进行与水平方向的递归方式类似的方法，计算得到新的像素值
（3）把垂直递归与水平递归的结果相加，然后除以2，得到最后的结果
这个算法具有高度并行的效果，因此启用多线程毫无压力，而且算法是图像中常用的递归加速，因此速度挺快的，如果还觉得速度不够快，可以采用查询表的方式
接着贴一下部分重要函数的代码，供参考，下面代码中我用了查询表的方式进行加速，速度可以比之前提高一倍，效果和原算法肉眼看不出来有什么区别，差别非常小，因此建议用查询表进行加速，还有其它的一些小细节我也没有根据原文进行写代码。
- void CBeeps::Run(BYTE* pImage, int nWidth, int nHeight, int nStride)
- {
- if (pImage == NULL)
- {
- return;
- }
- m_pImage=pImage;
- m_nWidth=nWidth;
- m_nHeight=nHeight;
- m_nStride=nStride;
- ApplyBiExponentialEdgePreservingSmoother(20,0.02);
- }
- void CBeeps::ApplyBiExponentialEdgePreservingSmoother(double photometricStandardDeviation, double spatialDecay)
- {
- m_exp_table=new double[256];
- m_g_table=new double[256];
- double c=-0.5/(photometricStandardDeviation * photometricStandardDeviation);
- double mu=spatialDecay/(2-spatialDecay);
- for (int i=0;i<=255;i++)
- {
- float a=exp(c*i*i);
- m_exp_table[i]=(1-spatialDecay)* exp(c*i*i);
- m_g_table[i]=mu*i;
- }
- BYTE* p0 =m_pImage;
- const int nChannel = 4;
- int m_length =m_nWidth*m_nHeight;
- float maxerror=0;
- float sum=0;
- // 对每个channel进行处理
- #pragma omp parallel for
- for (int idxChannel=0;idxChannel <nChannel; idxChannel++)
- {
- double *data1 = new double[m_length];
- double* data2 = new double[m_length];
- //double* data3 = new double[m_length];
- //double* data4 = new double[m_length];
- BYTE *p1=p0+idxChannel;
- for (int i = 0; i < m_length;++i)
- {
- data1[i] = p1[i * nChannel];
- }
- memcpy(data2,data1, sizeof(double) * m_length);
- //memcpy(data3,data1, sizeof(double) * m_length);
- //memcpy(data4,data1, sizeof(double) * m_length);
- //CBEEPSHorizontalVertical hv(data3, m_nWidth, m_nHeight, photometricStandardDeviation, spatialDecay);
- //hv.run();
- //CBEEPSVerticalHorizontal vh(data4, m_nWidth, m_nHeight, photometricStandardDeviation, spatialDecay);
- //vh.run();
- runHorizontalVertical(data1, m_nWidth, m_nHeight,spatialDecay,m_exp_table,m_g_table);
- runVerticalHorizontal(data2, m_nWidth, m_nHeight,spatialDecay,m_exp_table,m_g_table);
- sum=0;
- for (int i =0;i<m_length;++i)
- {
- //double val = (data3[i] + data4[i]) * 0.5;
- double val=(data1[i] + data2[i]) * 0.5;
- //double error0=abs((int)val2-(int)val);
- //sum=sum+error0;
- //if (error0>maxerror)
- //{
- //maxerror=error0;
- //}
- if(255.0<val)val=255.0;
- p1[i * nChannel]=(BYTE)val;
- }
- //sum=sum/m_length;
- delete data1;
- delete data2;
- }//for
- delete m_exp_table;m_exp_table=NULL;
- delete m_g_table;m_g_table=NULL;
- }
- //垂直方向递归
- void CBeeps::runVerticalHorizontal(double *data,int width,int height,double spatialDecay,double *exp_table,double *g_table)
- {
- int length0=height*width;
- double* g= new double[length0];
- int m = 0;
- for (int k2 = 0;k2<height;++k2)
- {
- int n = k2;
- for (int k1 = 0;k1<width;++k1)
- {
- g[n]=data[m++];
- n += height;
- }
- }
- double*p = new double[length0];
- double*r = new double[length0];
- memcpy(p, g, sizeof(double) * length0);
- memcpy(r, g, sizeof(double) * length0);
- for (int k1 = 0;k1<width; ++k1)
- {
- int startIndex=k1 * height;
- double mu = 0.0;
- for (int k=startIndex+1,K =startIndex+height;k<K;++k)
- {
- int div0=fabs(p[k]-p[k-1]);
- mu =exp_table[div0];
- p[k] = p[k - 1] * mu + p[k] * (1.0 - mu);//文献中的公式1，这里做了一下修改，效果影响不大
- }
- for (int k =startIndex+height-2;startIndex <= k;--k)
- {
- int div0=fabs(r[k]-r[k+1]);
- mu =exp_table[div0];
- r[k] = r[k+1] * mu + r[k] * (1.0-mu) ;//文献公式3
- }
- }
- double rho0=1.0/(2-spatialDecay);
- for (int k = 0;k <length0;++k)
- {
- r[k]= (r[k]+p[k])*rho0-g_table[(int)g[k]];
- }
- m = 0;
- for (int k1=0;k1<width;++k1)
- {
- int n = k1;
- for (int k2 =0;k2<height;++k2)
- {
- data[n] = r[m++];
- n += width;
- }
- }
- memcpy(p,data, sizeof(double) * length0);
- memcpy(r,data, sizeof(double) * length0);
- for (int k2 = 0; k2<height;++k2)
- {
- int startIndex=k2 * width;
- double mu = 0.0;
- for (int k=startIndex+1, K=startIndex+width;k<K;++k)
- {
- int div0=fabs(p[k]-p[k-1]);
- mu =exp_table[div0];
- p[k] = p[k - 1] * mu + p[k] * (1.0 - mu);
- }
- for (int k=startIndex+width-2;startIndex<=k;--k)
- {
- int div0=fabs(r[k]-r[k+1]);
- mu =exp_table[div0];
- r[k] = r[k + 1] * mu + r[k] * (1.0 - mu) ;
- }
- }
- double init_gain_mu=spatialDecay/(2-spatialDecay);
- for (int k = 0; k <length0; k++)
- {
- data[k]=(p[k]+r[k])*rho0-data[k]*init_gain_mu;//文献中的公式5
- }
- delete p;
- delete r;
- delete g;
- }
- //水平方向递归
- void CBeeps::runHorizontalVertical(double *data,int width,int height,double spatialDecay,double *exptable,double *g_table)
- {
- int length=width*height;
- double* g = new double[length];
- double* p = new double[length];
- double* r = new double[length];
- memcpy(p,data, sizeof(double) * length);
- memcpy(r,data, sizeof(double) * length);
- double rho0=1.0/(2-spatialDecay);
- for (int k2 = 0;k2 < height;++k2)
- {
- int startIndex=k2 * width;
- for (int k=startIndex+1,K=startIndex+width;k<K;++k)
- {
- int div0=fabs(p[k]-p[k-1]);
- double mu =exptable[div0];
- p[k] = p[k - 1] * mu + p[k] * (1.0 - mu);//文献公式1
- }
- for (int k =startIndex + width - 2;startIndex <= k;--k)
- {
- int div0=fabs(r[k]-r[k+1]);
- double mu =exptable[div0];
- r[k] = r[k + 1] * mu + r[k] * (1.0 - mu);//文献公式3
- }
- for (int k =startIndex,K=startIndex+width;k<K;k++)
- {
- r[k]=(r[k]+p[k])*rho0- g_table[(int)data[k]];
- }
- }
- int m = 0;
- for (int k2=0;k2<height;k2++)
- {
- int n = k2;
- for (int k1=0;k1<width;k1++)
- {
- g[n] = r[m++];
- n += height;
- }
- }
- memcpy(p, g, sizeof(double) * height * width);
- memcpy(r, g, sizeof(double) * height * width);
- for (int k1=0;k1<width;++k1)
- {
- int startIndex=k1 * height;
- double mu = 0.0;
- for (int k =startIndex+1,K =startIndex+height;k<K;++k)
- {
- int div0=fabs(p[k]-p[k-1]);
- mu =exptable[div0];
- p[k] = p[k - 1] * mu + p[k] * (1.0 - mu);
- }
- for (int k=startIndex+height-2;startIndex<=k;--k)
- {
- int div0=fabs(r[k]-r[k+1]);
- mu =exptable[div0];
- r[k] = r[k + 1] * mu + r[k] * (1.0 - mu);
- }
- }
- double init_gain_mu=spatialDecay/(2-spatialDecay);
- for (int k = 0;k <length;++k)
- {
- r[k]= (r[k]+p[k])*rho0- g[k]*init_gain_mu;
- }
- m = 0;
- for (int k1=0;k1<width;++k1)
- {
- int n = k1;
- for (int k2=0;k2<height;++k2)
- {
- data[n]=r[m++];
- n += width;
- }
- }
- delete p;
- delete r;
- delete g;
- }
最后看一下我用这个算法，实现磨皮的结果：
原图
美图秀秀磨皮结果
双指数递归算法磨皮结果
扫一扫
1689
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
