BEEPS-仿美图秀秀磨皮算法，让美女的皮肤更光滑

  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);
 61 end
 62 for i=(len-1):-1:1
 63     Phi1(i)=(1-lambda*exp((-(X1(i)-Phi1(i+1))^2)/(2*sigma^2)))*X1(i)+lambda*exp((-(X1(i)-Phi1(i+1))^2)/(2*sigma^2))*Phi1(i+1);
 64     Phi2(i)=(1-lambda*exp((-(X2(i)-Phi2(i+1))^2)/(2*sigma^2)))*X2(i)+lambda*exp((-(X2(i)-Phi2(i+1))^2)/(2*sigma^2))*Phi2(i+1);
 65     Phi3(i)=(1-lambda*exp((-(X3(i)-Phi3(i+1))^2)/(2*sigma^2)))*X3(i)+lambda*exp((-(X3(i)-Phi3(i+1))^2)/(2*sigma^2))*Phi3(i+1);
 66 end
 67 for i=1:len
 68     Y1(i)=(Psi1(i)-(1-lambda)*X1(i)+Phi1(i))/(1+lambda);
 69     Y2(i)=(Psi2(i)-(1-lambda)*X2(i)+Phi2(i))/(1+lambda);
 70     Y3(i)=(Psi3(i)-(1-lambda)*X3(i)+Phi3(i))/(1+lambda);
 71 end
 72 for i=1:row
 73     for j=1:column
 74         tempimg1(i,j)=Y1((i-1)*column+j);
 75         tempimg3(i,j)=Y2((i-1)*column+j);
 76         tempimg5(i,j)=Y3((i-1)*column+j);
 77         
 78     end
 79 end
 80 
 81 % 2.vertical processing
 82 for j=1:column
 83     for i=1:row
 84         X1((j-1)*row+i)=tempimg1(i,j);
 85         X2((j-1)*row+i)=tempimg3(i,j);
 86         X3((j-1)*row+i)=tempimg5(i,j);
 87 
 88     end
 89 end
 90 Psi1(1)=X1(1);
 91 Phi1(len)=X1(len);
 92 Psi2(1)=X2(1);
 93 Phi2(len)=X2(len);
 94 Psi3(1)=X3(1);
 95 Phi3(len)=X3(len);
 96 for i=2:len
 97     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);
 98     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);
 99     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);
100 end
101 for i=(len-1):-1:1
102     Phi1(i)=(1-lambda*exp((-(X1(i)-Phi1(i+1))^2)/(2*sigma^2)))*X1(i)+lambda*exp((-(X1(i)-Phi1(i+1))^2)/(2*sigma^2))*Phi1(i+1);
103     Phi2(i)=(1-lambda*exp((-(X2(i)-Phi2(i+1))^2)/(2*sigma^2)))*X2(i)+lambda*exp((-(X2(i)-Phi2(i+1))^2)/(2*sigma^2))*Phi2(i+1);
104     Phi3(i)=(1-lambda*exp((-(X3(i)-Phi3(i+1))^2)/(2*sigma^2)))*X3(i)+lambda*exp((-(X3(i)-Phi3(i+1))^2)/(2*sigma^2))*Phi3(i+1);
105 end
106 for i=1:len
107     Y1(i)=(Psi1(i)-(1-lambda)*X1(i)+Phi1(i))/(1+lambda);
108     Y2(i)=(Psi2(i)-(1-lambda)*X2(i)+Phi2(i))/(1+lambda);
109     Y3(i)=(Psi3(i)-(1-lambda)*X3(i)+Phi3(i))/(1+lambda);
110 end
111 for j=1:column
112     for i=1:row
113         tempimg1(i,j)=Y1((j-1)*row+i);
114         tempimg3(i,j)=Y2((j-1)*row+i);
115         tempimg5(i,j)=Y3((j-1)*row+i);
116 
117     end
118 end
119 
120 % vertical-horizon
121 % 1.vertical processing
122 for j=1:column
123     for i=1:row
124         X1((j-1)*row+i)=Rimg(i,j);
125         X2((j-1)*row+i)=Gimg(i,j);
126         X3((j-1)*row+i)=Bimg(i,j);
127     end
128 end
129 Psi1(1)=X1(1);
130 Phi1(len)=X1(len);
131 Psi2(1)=X2(1);
132 Phi2(len)=X2(len);
133 Psi3(1)=X3(1);
134 Phi3(len)=X3(len);
135 for i=2:len
136     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);
137     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);
138     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);
139 end
140 for i=(len-1):-1:1
141     Phi1(i)=(1-lambda*exp((-(X1(i)-Phi1(i+1))^2)/(2*sigma^2)))*X1(i)+lambda*exp((-(X1(i)-Phi1(i+1))^2)/(2*sigma^2))*Phi1(i+1);
142     Phi2(i)=(1-lambda*exp((-(X2(i)-Phi2(i+1))^2)/(2*sigma^2)))*X2(i)+lambda*exp((-(X2(i)-Phi2(i+1))^2)/(2*sigma^2))*Phi2(i+1);
143     Phi3(i)=(1-lambda*exp((-(X3(i)-Phi3(i+1))^2)/(2*sigma^2)))*X3(i)+lambda*exp((-(X3(i)-Phi3(i+1))^2)/(2*sigma^2))*Phi3(i+1);
144 end
145 for i=1:len
146     Y1(i)=(Psi1(i)-(1-lambda)*X1(i)+Phi1(i))/(1+lambda);
147     Y2(i)=(Psi2(i)-(1-lambda)*X2(i)+Phi2(i))/(1+lambda);
148     Y3(i)=(Psi3(i)-(1-lambda)*X3(i)+Phi3(i))/(1+lambda);
149 end
150 for j=1:column
151     for i=1:row
152         tempimg2(i,j)=Y1((j-1)*row+i);
153         tempimg4(i,j)=Y2((j-1)*row+i);
154         tempimg6(i,j)=Y3((j-1)*row+i);
155 
156     end
157 end 
158 % 2.horizon processing
159 for i=1:row
160     for j=1:column
161         X1((i-1)*column+j)=tempimg2(i,j);
162         X2((i-1)*column+j)=tempimg4(i,j);
163         X3((i-1)*column+j)=tempimg6(i,j);
164 
165     end
166 end
167 Psi1(1)=X1(1);
168 Phi1(len)=X1(len);
169 Psi2(1)=X2(1);
170 Phi2(len)=X2(len);
171 Psi3(1)=X3(1);
172 Phi3(len)=X3(len);
173 for i=2:len
174     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);
175     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);
176     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);
177 end
178 for i=(len-1):-1:1
179     Phi1(i)=(1-lambda*exp((-(X1(i)-Phi1(i+1))^2)/(2*sigma^2)))*X1(i)+lambda*exp((-(X1(i)-Phi1(i+1))^2)/(2*sigma^2))*Phi1(i+1);
180     Phi2(i)=(1-lambda*exp((-(X2(i)-Phi2(i+1))^2)/(2*sigma^2)))*X2(i)+lambda*exp((-(X2(i)-Phi2(i+1))^2)/(2*sigma^2))*Phi2(i+1);
181     Phi3(i)=(1-lambda*exp((-(X3(i)-Phi3(i+1))^2)/(2*sigma^2)))*X3(i)+lambda*exp((-(X3(i)-Phi3(i+1))^2)/(2*sigma^2))*Phi3(i+1);
182 end
183 for i=1:len
184     Y1(i)=(Psi1(i)-(1-lambda)*X1(i)+Phi1(i))/(1+lambda);
185     Y2(i)=(Psi2(i)-(1-lambda)*X2(i)+Phi2(i))/(1+lambda);
186     Y3(i)=(Psi3(i)-(1-lambda)*X3(i)+Phi3(i))/(1+lambda);
187 end
188 for i=1:row
189     for j=1:column
190         tempimg2(i,j)=Y1((i-1)*column+j);
191         tempimg4(i,j)=Y2((i-1)*column+j);
192         tempimg6(i,j)=Y3((i-1)*column+j);
193         
194     end
195 end
196 
197 tempimg7=(tempimg1+tempimg2)/2;
198 tempimg8=(tempimg3+tempimg4)/2;
199 tempimg9=(tempimg5+tempimg6)/2;
200 
201 img(:,:,1)=tempimg7;
202 img(:,:,2)=tempimg8;
203 img(:,:,3)=tempimg9;
204 figure,imshow(uint8(img));