forked from yqueau/normal_integration
-
Notifications
You must be signed in to change notification settings - Fork 0
/
demo_2_quadratic.m
148 lines (125 loc) · 4.08 KB
/
demo_2_quadratic.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
140
141
142
143
144
145
146
147
148
clear
close all
addpath('Toolbox/');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load a dataset containing:
% -- p (nrows x ncols) : gradient in the u- (bottom) direction
% -- q (nrows x ncols) : gradient in the v- (right) direction
% -- u (nrows x ncols) : ground truth depth map
% -- mask (nrows x ncols) : mask of the pixels on the vase (binary)
load Datasets/vase
indices_mask = find(mask>0); % Indices of the pixel inside the mask
% Add zero-mean, Gaussian noise inside the mask
std_noise = 0.02*max(sqrt(p(indices_mask).^2+q(indices_mask).^2));
p = p+std_noise*randn(size(p));
q = q+std_noise*randn(size(q));
% Fill the gradient with 0 to test rectangular integration
p(mask==0) = 0;
q(mask==0) = 0;
% Remove the ground truth depth values outside the mask
u(mask==0) = NaN;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set optimization parameters
lambda = 1e-6*ones(size(p)); % Uniform field of weights (nrows x ncols)
z0 = zeros(size(p)); % Null depth prior (nrows x ncols)
solver = 'pcg'; % Solver ('pcg' means conjugate gradient, 'direct' means backslash i.e. sparse Cholesky)
precond = 'CMG'; % Preconditioner ('none' means no preconditioning, 'ichol' means incomplete Cholesky, 'CMG' means conjugate combinatorial multigrid -- the latter is fastest, but it need being installed, see README)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Integrate without mask (full rectangular domain)
Omega_1 = ones(size(p));
t_1 = tic;
z_1 = smooth_integration(p,q,Omega_1,lambda,z0,solver,precond);
t_1 = toc(t_1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Integrate with mask (reduced non-rectangular domain)
Omega_2 = mask;
t_2 = tic;
z_2 = smooth_integration(p,q,Omega_2,lambda,z0,solver,precond);
t_2 = toc(t_2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Evaluate both integration methods over the mask
% Find the constant of integration which minimizes RMSE wrt ground truth
lambda_1 = -mean(z_1(indices_mask)-u(indices_mask));
z_1 = z_1+lambda_1;
lambda_2 = -mean(z_2(indices_mask)-u(indices_mask));
z_2 = z_2+lambda_2;
% Calculate RMSEs
RMSE_1 = sqrt(mean((z_1(indices_mask)-u(indices_mask)).^2));
RMSE_2 = sqrt(mean((z_2(indices_mask)-u(indices_mask)).^2));
% Display evaluation results in terminal
disp('=============================');
disp('Integration without mask:');
disp(sprintf('CPU: %.4f',t_1));
disp(sprintf('RMSE over mask: %.2f',RMSE_1));
disp('');
disp('================================');
disp('Integration with mask:');
disp(sprintf('CPU: %.4f',t_2));
disp(sprintf('RMSE over mask: %.2f',RMSE_2));
disp('');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display a few things
figure('units','normalized','outerposition',[0 0 1 1])
% Input data: p, q and mask
subplot(3,3,1)
imagesc(p);
axis image
axis off
title('$$p$$','Interpreter','Latex','Fontsize',14)
subplot(3,3,2)
imagesc(q);
axis image
axis off
title('$$q$$','Interpreter','Latex','Fontsize',14)
subplot(3,3,3)
imagesc(mask);
axis image
axis off
colormap gray
title('$$\Omega$$','Interpreter','Latex','Fontsize',14)
subplot(3,3,4)
surfl(u,[-135 30]);
view(-35,20)
axis ij;
shading flat;
colormap gray;
axis equal;
grid off
axis off
title('Ground truth depth','Interpreter','Latex','Fontsize',14)
subplot(3,3,5)
surfl(z_1,[-135 30]);
view(-35,20)
axis ij;
shading flat;
colormap gray;
axis equal;
grid off
axis off
title('Reconstruction without mask','Interpreter','Latex','Fontsize',14)
subplot(3,3,6)
surfl(z_2,[-135 30]);
view(-35,20)
axis ij;
shading flat;
colormap gray;
axis equal;
grid off
axis off
title('Reconstruction with mask','Interpreter','Latex','Fontsize',14)
error_map_1 = abs(u-z_1);
error_map_1(mask==0) = NaN;
error_map_2 = abs(u-z_2);
error_map_2(mask==0) = NaN;
subplot(3,3,8)
imagesc(error_map_1,[0 5]);
axis image
axis off
colormap gray
title('Absolute error without mask','Interpreter','Latex','Fontsize',14)
subplot(3,3,9)
imagesc(error_map_2,[0 5]);
axis image
axis off
colormap gray
title('Absolute error with mask','Interpreter','Latex','Fontsize',14)