/
demo.m
139 lines (110 loc) · 4.64 KB
/
demo.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
clear all
%%
addpath('ZhuRamananDetector','optimisations','utils','comparison');
% YOU MUST set this to the base directory of the Basel Face Model
BFMbasedir = '';
% Load morphable model
load(strcat(BFMbasedir,'01_MorphableModel.mat'));
% Important to use double precision for use in optimisers later
shapeEV = double(shapeEV);
shapePC = double(shapePC);
shapeMU = double(shapeMU);
% We need edge-vertex and edge-face lists to speed up finding occluding boundaries
% Either: 1. Load precomputed edge structure for Basel Face Model
load('BFMedgestruct.mat');
% Or: 2. Compute lists for new model:
%TR = triangulation(tl,ones(k,1),ones(k,1),ones(k,1));
%Ev = TR.edges;
%clear TR;
%Ef = meshFaceEdges( tl,Ev );
%save('edgestruct.mat','Ev','Ef');
%% ADJUSTABLE PARAMETERS
% Number of model dimensions to use
ndims = 60;
% Prior weight for initial landmark fitting
w_initialprior=0.7;
% Number of iterations for iterative closest edge fitting
icefniter=7;
options.Ef = Ef;
options.Ev = Ev;
% w1 = weight for edges
% w2 = weight for landmarks
% w3 = 1-w1-w2 = prior weight
options.w1 = 0.45;
options.w2 = 0.15;
%% Setup basic parameters
testdir='testImages/';
im = imread(strcat(testdir,'image_0018.png'));
edgeim = edge(rgb2gray(im),'canny',0.15);
scale=1;
testdir='C:\Users\qcri\Documents\Ajay\Face\SOA\pami09\';
testfname = strcat(testdir,'ss3.PNG');
im = imread(testfname);
im_size = size(im(:,:,1));
im = imresize(im,scale);
edgeim = edge(rgb2gray(im),'canny',0.15);
ZRtimestart = tic;
bs = LandmarkDetector(im);
ZRtime = toc(ZRtimestart);
disp(['Time for landmark detection: ' num2str(ZRtime) ' seconds']);
[x_landmarks,landmarks]=ZR2BFM( bs,im );
%% Initialise using only landmarks
disp('Fitting to landmarks only...');
[b,R,t,s] = FitSingleSOP( x_landmarks,shapePC,shapeMU,shapeEV,ndims,landmarks,w_initialprior );
FV.vertices=reshape(shapePC(:,1:ndims)*b+shapeMU,3,size(shapePC,1)/3)';
FV.faces = tl;
%% Initialise using iterative closest edge fitting (ICEF)
disp('Fitting to edges with iterative closest edge fitting...');
[b,R,t,s] = FitEdges(im,x_landmarks,landmarks,shapePC,shapeMU,shapeEV,options.Ef,options.Ev,tl,ndims, w_initialprior, options.w1, options.w2,icefniter);
FV.vertices = reshape(shapePC(:,1:ndims)*b+shapeMU,3,size(shapePC,1)/3)';
%% Run final optimisation of hard edge cost
disp('Optimising non-convex edge cost...');
maxiter = 5;
iter = 0;
diff = 1;
eps = 1e-9;
[r,c]=find(edgeim);
r = size(edgeim,1)+1-r;
while (iter<maxiter) && (diff>eps)
FV.vertices=reshape(shapePC(:,1:ndims)*b+shapeMU,3,size(shapePC,1)/3)';
[ options.occludingVertices ] = occludingBoundaryVertices( FV,options.Ef,options.Ev,R );
X = reshape(shapePC(:,1:ndims)*b+shapeMU,3,size(shapePC(:,1:ndims),1)/3);
% Compute position of projected occluding boundary vertices
x_edge = R*X(:,options.occludingVertices);
x_edge = x_edge(1:2,:);
x_edge(1,:)=s.*(x_edge(1,:)+t(1));
x_edge(2,:)=s.*(x_edge(2,:)+t(2));
% Find edge correspondences
[idx,d] = knnsearch([c r],x_edge');
% Filter edge matches - ignore the worse 5%
sortedd=sort(d);
threshold = sortedd(round(0.95*length(sortedd)));
idx = idx(d<threshold);
options.occludingVertices = options.occludingVertices(d<threshold);
b0 = b;
[ R,t,s,b ] = optimiseHardEdgeCost( b0,x_landmarks,shapeEV,shapeMU,shapePC,R,t,s,r,c,landmarks,options,tl,false );
diff = norm(b0-b);
disp(num2str(diff));
iter = iter+1;
end
% Run optimisation for a final time but without limit on number of
% iterations
[ R,t,s,b ] = optimiseHardEdgeCost( b,x_landmarks,shapeEV,shapeMU,shapePC,R,t,s,r,c,landmarks,options,tl,true );
disp('Rendering final results...');
FV.vertices=reshape(shapePC(:,1:ndims)*b+shapeMU,3,size(shapePC,1)/3)';
figure; subplot(1,3,1); patch(FV, 'FaceColor', [1 1 1], 'EdgeColor', 'none', 'FaceLighting', 'phong'); light; axis equal; axis off;
subplot(1,3,2); imshow(renderFace(FV,im,R,t,s,false));
subplot(1,3,3); imshow(renderFace(FV,im,R,t,s,true));
%% save data 3dmodel
save([testfname(1:end-3) 'mat'],'R','t','s','b');
% %% save data for sig17 MTP
% [normals,mask,depthMap]=render_face_normals(FV,im,R,t,s);
% depthMap(:,:,1)=(depthMap(:,:,1)./max(max(depthMap(:,:,1))));
% depthMap(:,:,2)=(depthMap(:,:,2)./max(max(depthMap(:,:,2))));
% depthMap(:,:,3)=(depthMap(:,:,3)/scale/scale);
% depthMap = imresize(depthMap,im_size);
% imwrite(depthMap,'C:\Users\qcri\Documents\Ajay\Face\SOA\sig_2017_mtp\img\depth_010.png');
% nornmals = imresize(normals,im_size);
% imwrite((normals+1)/2,'C:\Users\qcri\Documents\Ajay\Face\SOA\sig_2017_mtp\img\normals_010.png');
% mask = imresize(mask,im_size);
% imwrite(mask(:,:,1),'C:\Users\qcri\Documents\Ajay\Face\SOA\sig_2017_mtp\img\masks_010.png');