img = imread('marui2-warp.bmp');
high = imread('Z:\REAL_HOME\cnn\pred\marui2-warp-meannorm.jpg-hi28154.bmp');%Highlight-net output
shading = imread('Z:\REAL_HOME\cnn\pred\marui2-hf.bmp-sh-hf-2035.bmp');%albedo-shading net output: diffuse shading
albedo = imread('Z:\REAL_HOME\cnn\pred\marui2-hf.bmp-ab-hf-2035.bmp');%albedo-shading net output: albedo layer
warpmean = im2double(imread('marui2-normal.bmp'));%the surface normal, if you do not have other codes to estimate normals, you can use this one as approximation, since the landmarks are aligned
mask = warpmean(:,:,3)>0;
mask1 = reshape(mask,[90000 1]);
mask2 = rgb2gray(img)>0;
mask = mask.*mask2;
mask3 = repmat(mask,[1 1 3]);
albedo = im2double(albedo);
albedo = imresize(albedo,[300 300]);
albedo = reshape(albedo,[90000 3]);
albedo = albedo(find(mask==1),:);
medab = median(albedo);
load mean_diffuseab
load mean_rough
load mean_speab 
dis = repmat(medab,[6,1])-mean_diffuseab;
dis = sum(abs(dis),2);
skinty = find(dis==min(dis(:)))
channel = 3;
wm = warpmean.*mask3 - 0.5;
wm = wm*2;
[w1 w2 w3] = size(warpmean);
wm = reshape(wm,[w1*w2,3]);
d = wm.*wm;
d = sum(d,2);
d = sqrt(d);
d = repmat(d,[1 3]);
wm = wm./d;
wm = reshape(wm,[w1 w2 3]);
mask3 = repmat(mask,[1 1 3]);
wm = wm.*mask3;
highold=im2double(high).*mask3;
highg = rgb2gray(highold);
shading = imresize(im2double(shading),[300 300]);
highg = highg./max(highg(:));
higholdrescale=(im2double(shading).*repmat(highold(:,:,3),[1 1 3]));
higholdrescale =higholdrescale./max(higholdrescale(:))*max(max(highold(:,:,3)));
clear final
%%
% get masks for 10 face regions, if you are not using the mean normal file
% we provided, uncomment this section
%if you are using the normal file we provided, just load
% regionmask.mat and keep this section commented, but using mean normal may
% not get very accurate result
load regionmask
% regionmask = zeros(100,200,10);
% tic
% for iiiii=1:10
% lighting=zeros(100,200);
% name = ['region' num2str(iiiii)];
% highg = im2double(imread([name '-mask.bmp']));
% maskg = highg(:,:,1);
% highg = highg(:,:,1);
% figure,imshow(highg)
% highg = reshape(highg,[90000 1]);
% wmlarge = reshape(wm,[90000 3]);
% clear find
% [a b] = find(maskg==1);
% clear length
% maxtheta= 0;
% for indexi=1:100
%     for indexj = 1:200
% fi = (100-indexi)/100*pi;
% theta = indexj/200*2*pi;
% [x y z] = sph2xyz(theta,fi);
% ld= zeros(3,1);
% ld(1) = y;
% ld(2) = -z;
% ld(3)=-x;
% v = [0 0 1];
% nm = ld'+v;
% nm = nm./sqrt(sum(nm.*nm));
% lengnm = length(wmlarge);
% nm = repmat(nm,[lengnm 1]);
% dis = nm-wmlarge;
% dis = abs(dis);
% dis = sum(dis,2);
% tmp=sort(dis);
% for sorti=1:10
%     ind = find(dis==tmp(sorti));
%     lighting(indexi,indexj) = max(lighting(indexi,indexj),highg(ind(1)));
% end
%     end
% end
% regionmask(:,:,iiiii) =lighting;
% end
% toc
%%
% project to the environment map and deconvolution
finalr = zeros(100,200,10);
finalg = zeros(100,200,10);
finalb = zeros(100,200,10);
tic
for iiiii=1:3
highg= higholdrescale(:,:,iiiii);
lighting=zeros(100,200);
maskg = mask(:,:,1);
highg = reshape(highg,[90000 1]);
wmlarge = reshape(wm,[90000 3]);
clear length
maxtheta= 0;
for indexi=1:100
    for indexj = 1:200
fi = (100-indexi)/100*pi;
theta = indexj/200*2*pi;
[x y z] = sph2xyz(theta,fi);
ld= zeros(3,1);
ld(1) = y;
ld(2) = -z;
ld(3)=-x;
v = [0 0 1];
nm = ld'+v;
nm = nm./sqrt(sum(nm.*nm));
lengnm = length(wmlarge);
nm = repmat(nm,[lengnm 1]);
dis = nm-wmlarge;
dis = abs(dis);
dis = sum(dis,2);
tmp=sort(dis);
for sorti=1:10
    ind = find(dis==tmp(sorti));
    lighting(indexi,indexj) = max(lighting(indexi,indexj),highg(ind(1)));
end
    end
end
for region=1:10
    lighting1 =lighting/mean_speab(skinty,region);
    ksize = mean_rough(skinty,region);
    for i=1:100
        quickconv;
        
        if iiiii==1
            finalr(i,:,region)=J(i,:);
        end
        if iiiii==2
            finalg(i,:,region)=J(i,:);
        end
        if iiiii==3
            finalb(i,:,region)=J(i,:);
        end
    end
end
end
toc
% combine the results of 10 regions
i=1;
final(:,:,1) = finalr(:,:,1).*regionmask(:,:,i);
final(:,:,2) = finalg(:,:,1).*regionmask(:,:,i);
final(:,:,3) = finalb(:,:,1).*regionmask(:,:,i);
for i=2:10
    final(:,:,1) = max(final(:,:,1),finalr(:,:,i).*regionmask(:,:,i));
    final(:,:,2) = max(final(:,:,2),finalg(:,:,i).*regionmask(:,:,i));
    final(:,:,3) = max(final(:,:,3),finalb(:,:,i).*regionmask(:,:,i));
end
figure,imshow(final./max(final(:)))
final(find(final<0))=0;
hdrwrite(final,'result.hdr')