beeps磨皮算法的Matlab实现：双指数平滑处理

 1 % Reference:
  2 % 该算法可以用来去噪，同时还可以用来磨皮，具有较好的保护边缘细节的能力
  3 % Philippe Thevenaz et al,"Bi-Exponential Edge-Preserving Smoother".
  4 clear all;
  5 close all;
  6 [FileName,PathName]=uigetfile({'*.jpg';'*.png';'*.bmp';'*.jpeg'},'Open an Image File');
  7 img = imread([PathName FileName]);
  8 figure,imshow(img);
  9 img=double(img);
 10 Rimg=img(:,:,1);
 11 Gimg=img(:,:,2);
 12 Bimg=img(:,:,3);
 13 % 初始化参数
 14 [row,column]=size(Rimg);
 15 tempimg1=zeros(row,column);
 16 tempimg2=zeros(row,column);
 17 tempimg3=zeros(row,column);
 18 tempimg4=zeros(row,column);
 19 tempimg5=zeros(row,column);
 20 tempimg6=zeros(row,column);
 21 len=numel(Rimg(:));
 22 Psi1=zeros(1,len);
 23 Phi1=zeros(1,len);
 24 Psi2=zeros(1,len);
 25 Phi2=zeros(1,len);
 26 Psi3=zeros(1,len);
 27 Phi3=zeros(1,len);
 28 X1=zeros(1,len);
 29 X2=zeros(1,len);
 30 X3=zeros(1,len);

 31 Y1=zeros(1,len);
 32 Y2=zeros(1,len);
 33 Y3=zeros(1,len);
 34 
 35 lambda=1.05;
 36 % 注：参数sigma的值越大，图像越模糊，如果需要较大程度地磨皮，应该在保持sigma值较小的前提下，逐渐增大lambda的值
 37 % 这样才能使图像不会变得太模糊
 38 sigma=14;
 39 % 对RGB三个通道分别处理
 40 % horizon-vertical processing
 41 % 1.horizon processing
 42 for i=1:row
 43     for j=1:column
 44         
 45         X1((i-1)*column+j)=Rimg(i,j);
 46         X2((i-1)*column+j)=Gimg(i,j);
 47         X3((i-1)*column+j)=Bimg(i,j);
 48 
 49     end
 50 end
 51 Psi1(1)=X1(1);
 52 Phi1(len)=X1(len);
 53 Psi2(1)=X2(1);
 54 Phi2(len)=X2(len);
 55 Psi3(1)=X3(1);
 56 Phi3(len)=X3(len);
 57 for i=2:len
 58     Psi1(i)=(1-lambda*exp((-(X1(i)-Psi1(i-1))^2)/(2*sigma^2)))*X1(i)+lambda*exp((-(X1(i)-Psi1(i-1))^2)/(2*sigma^2))*Psi1(i-1);
 59     Psi2(i)=(1-lambda*exp((-(X2(i)-Psi2(i-1))^2)/(2*sigma^2)))*X2(i)+lambda*exp((-(X2(i)-Psi2(i-1))^2)/(2*sigma^2))*Psi2(i-1);
 60     Psi3(i)=(1-lambda*exp((-(X3(i)-Psi3(i-1))^2)/(2*sigma^2)))*X3(i)+lambda*exp((-(X3(i)-Psi3(i-1))^2)/(2*sigma^2))*Psi3(i-1);