磨皮

% Reference: % 该算法可以用来去噪，同时还可以用来磨皮，具有较好的保护边缘细节的能力 % Philippe Thevenaz et al,"Bi-Exponential Edge-Preserving Smoother". clear all; close all; [FileName,PathName]=uigetfile({'*.jpg';'*.png';'*.bmp';'*.jpeg'},'Open an Image File'); img = imread([PathName FileName]); figure,imshow(img); img=double(img); Rimg=img(:,:,1); Gimg=img(:,:,2); Bimg=img(:,:,3); % 初始化参数 [row,column]=size(Rimg); tempimg1=zeros(row,column); tempimg2=zeros(row,column); tempimg3=zeros(row,column); tempimg4=zeros(row,column); tempimg5=zeros(row,column); tempimg6=zeros(row,column); len=numel(Rimg(:)); Psi1=zeros(1,len); Phi1=zeros(1,len); Psi2=zeros(1,len); Phi2=zeros(1,len); Psi3=zeros(1,len); Phi3=zeros(1,len); X1=zeros(1,len); X2=zeros(1,len); X3=zeros(1,len); Y1=zeros(1,len); Y2=zeros(1,len); Y3=zeros(1,len); lambda=1; % 注：参数sigma的值越大，图像越模糊，如果需要较大程度地磨皮，应该在保持sigma值较小的前提下，逐渐增大lambda的值 % 这样才能使图像不会变得太模糊 sigma=12; % 对RGB三个通道分别处理 % horizon-vertical processing % 1.horizon processing %行处理 图像像素行存放 for i=1:row for j=1:column X1((i-1)*column+j)=Rimg(i,j); X2((i-1)*column+j)=Gimg(i,j); X3((i-1)*column+j)=Bimg(i,j); end end Psi1(1)=X1(1); Phi1(len)=X1(len); Psi2(1)=X2(1); Phi2(len)=X2(len); Psi3(1)=X3(1); Phi3(len)=X3(len); for i=2:len 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); 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); 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); end for i=(len-1):-1:1 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); 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); 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); end for i=1:len Y1(i)=(Psi1(i)-(1-lambda)*X1(i)+Phi1(i))/(1+lambda); Y2(i)=(Psi2(i)-(1-lambda)*X2(i)+Phi2(i))/(1+lambda); Y3(i)=(Psi3(i)-(1-lambda)*X3(i)+Phi3(i))/(1+lambda); end %转化成矩阵 for i=1:row for j=1:column tempimg1(i,j)=Y1((i-1)*column+j); tempimg3(i,j)=Y2((i-1)*column+j); tempimg5(i,j)=Y3((i-1)*column+j); end end %列处理 像素按列存放 % 2.vertical processing for j=1:column for i=1:row X1((j-1)*row+i)=tempimg1(i,j); X2((j-1)*row+i)=tempimg3(i,j); X3((j-1)*row+i)=tempimg5(i,j); end end Psi1(1)=X1(1); Phi1(len)=X1(len); Psi2(1)=X2(1); Phi2(len)=X2(len); Psi3(1)=X3(1); Phi3(len)=X3(len); for i=2:len 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); 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); 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); end for i=(len-1):-1:1 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); 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); 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); end for i=1:len Y1(i)=(Psi1(i)-(1-lambda)*X1(i)+Phi1(i))/(1+lambda); Y2(i)=(Psi2(i)-(1-lambda)*X2(i)+Phi2(i))/(1+lambda); Y3(i)=(Psi3(i)-(1-lambda)*X3(i)+Phi3(i))/(1+lambda); end for j=1:column for i=1:row tempimg1(i,j)=Y1((j-1)*row+i); tempimg3(i,j)=Y2((j-1)*row+i); tempimg5(i,j)=Y3((j-1)*row+i); end end %先列后行 % vertical-horizon % 1.vertical processing for j=1:column for i=1:row X1((j-1)*row+i)=Rimg(i,j); X2((j-1)*row+i)=Gimg(i,j); X3((j-1)*row+i)=Bimg(i,j); end end Psi1(1)=X1(1); Phi1(len)=X1(len); Psi2(1)=X2(1); Phi2(len)=X2(len); Psi3(1)=X3(1); Phi3(len)=X3(len); for i=2:len 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); 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); 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); end for i=(len-1):-1:1 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); 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); 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); end for i=1:len Y1(i)=(Psi1(i)-(1-lambda)*X1(i)+Phi1(i))/(1+lambda); Y2(i)=(Psi2(i)-(1-lambda)*X2(i)+Phi2(i))/(1+lambda); Y3(i)=(Psi3(i)-(1-lambda)*X3(i)+Phi3(i))/(1+lambda); end for j=1:column for i=1:row tempimg2(i,j)=Y1((j-1)*row+i); tempimg4(i,j)=Y2((j-1)*row+i); tempimg6(i,j)=Y3((j-1)*row+i); end end % 2.horizon processing for i=1:row for j=1:column X1((i-1)*column+j)=tempimg2(i,j); X2((i-1)*column+j)=tempimg4(i,j); X3((i-1)*column+j)=tempimg6(i,j); end end Psi1(1)=X1(1); Phi1(len)=X1(len); Psi2(1)=X2(1); Phi2(len)=X2(len); Psi3(1)=X3(1); Phi3(len)=X3(len); for i=2:len 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); 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); 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); end for i=(len-1):-1:1 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); 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); 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); end for i=1:len Y1(i)=(Psi1(i)-(1-lambda)*X1(i)+Phi1(i))/(1+lambda); Y2(i)=(Psi2(i)-(1-lambda)*X2(i)+Phi2(i))/(1+lambda); Y3(i)=(Psi3(i)-(1-lambda)*X3(i)+Phi3(i))/(1+lambda); end for i=1:row for j=1:column tempimg2(i,j)=Y1((i-1)*column+j); tempimg4(i,j)=Y2((i-1)*column+j); tempimg6(i,j)=Y3((i-1)*column+j); end end tempimg7=(tempimg1+tempimg2)/2; tempimg8=(tempimg3+tempimg4)/2; tempimg9=(tempimg5+tempimg6)/2; img(:,:,1)=tempimg7; img(:,:,2)=tempimg8; img(:,:,3)=tempimg9; figure,imshow(uint8(img));