磨皮

% 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));