蛋疼之余，仿真了一下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));
