-
Notifications
You must be signed in to change notification settings - Fork 21
/
demo_1_calibrated_gray_LS_ps.m
128 lines (115 loc) · 3.96 KB
/
demo_1_calibrated_gray_LS_ps.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
clear
close all
addpath('Toolbox/');
dataset = 'Statuette';
do_write_obj = 1; % Set to 0 to de-active OBJ writing, to 1 to activate it
%%% Set dataset
% File containing the parameters of the LEDs
light_file = 'Datasets/light.mat';
% File containing the camera's intrinsics
camera_file = 'Datasets/camera.mat';
% Folder containing the photometric stereo images
images_folder = sprintf('Datasets/%s',dataset);
%%% Read calibration files
disp('Reading calibration files');
load(light_file);
load(camera_file);
if(K(1,3)==0) K = K';end % Matlab uses a transposed definition of intrinsics matrix
% Store in a compact structure
calib.S = S; clear S; % Locations (nimgs x 3)
calib.Dir = Dir; clear Dir; % Orientations (nimgs x 3)
calib.mu = mu; clear mu; % Anisotropy factors (nimgs x 1)
calib.Phi = Phi; clear Phi; % Intensities (nimgs x nchannels)
calib.K = K; clear K; % Intrinsics (3 x 3)
%%% Read dataset
disp('Reading photometric stereo images');
% Get number of images from light calibration
nimgs = size(calib.S,1);
% Read ambient light image
Iamb = double(imread(sprintf('%s/photometric_sample_raw_ambient.png',images_folder)));
[nrows,ncols,nchannels] = size(Iamb);
% Read mask
mask = double(imread(sprintf('%s/photometric_sample_mask_raw.png',images_folder)));
mask = mask(:,:,1)>0;
% Create the cos^4 alpha map to compensate for darkening at the borders
[xx,yy] = meshgrid(1:ncols,1:nrows); % pixel coordinates
xx = xx-calib.K(1,3); % x coordinates w.r.t. principal point
yy = yy-calib.K(2,3); % y coordinates w.r.t. principal point
ff = mean([calib.K(1,1),calib.K(2,2)]); % focal length
cos4a = (ff./sqrt(xx.^2+yy.^2+ff^2)).^4; % cosine-forth law of attenuation
clear xx yy ff
% Read each image, substract ambient light, and compensate for cos^4 alpha
I = zeros(nrows,ncols,nchannels,nimgs);
for i = 1:nimgs
Ii = double(imread(sprintf('%s/photometric_sample_raw_%04d.png',images_folder,i)));
Ii = bsxfun(@rdivide,Ii-Iamb,cos4a);
I(:,:,:,i) = max(0,Ii);
end
clear i Ii cos4a Iamb light_file camera_file images_folder nrows ncols nchannels
% Store in a compact structure
data.I = I; clear I; % Images (nrows x ncols x nchannels x nimgs)
data.mask = mask; clear mask; % Images (nrows x ncols x nchannels x nimgs)
%%% Convert images to graylevel for this demo
data.I = squeeze(mean(data.I,3));
calib.Phi = mean(calib.Phi,2);
%%% Show the input images
disp('Displaying data');
figure(1001)
subplot(3,1,1)
imshow(uint8(255*data.I(:,:,1)./max(data.I(:))))
title('$$I^1$$','Interpreter','Latex','Fontsize',18);
subplot(3,1,2)
imshow(uint8(255*data.I(:,:,2)./max(data.I(:))))
title('$$I^2$$','Interpreter','Latex','Fontsize',18);
subplot(3,1,3)
imshow(uint8(255*data.mask))
title('$$\Omega$$','Interpreter','Latex','Fontsize',18);
drawnow
%%% Set parameters
disp('Setting parameters');
params.precond = 'ichol'; % Preconditioner
params.z0 = 500*ones(size(data.mask)); % Initial depth map: a plane at 500mm from camera
params.estimator = 'LS'; % Least-squares estimation
params.scales = 8; % Number of pyramid levels
params.ratio = 1; % Downsample by ratio>1 for faster results
params.self_shadows = 1; % Do or do not explicitly take into account self-shadows
params.display = 1; % Do not display result at each iteration
params.zeta = 0.0; % Fidelity to z0 ==> CHRISTOPH: THIS PARAMETER CONTROLS WHAT YOU WANT
%%% Solve photometric stereo
disp('Solving photometric stereo');
[XYZ,N,rho,Phi,mask] = near_ps(data,calib,params);
% Scale albedo for visualization
rho = rho./max(rho(:));
rho = uint8(255*rho);
%%% Show the result
disp('Displaying results');
figure(1002);
surfl(-XYZ(:,:,1),-XYZ(:,:,2),-XYZ(:,:,3),[0 90]);
shading flat;
colormap gray
view(-220,30);
axis ij
axis image
title('Shape')
figure(1003)
imagesc(rho)
if(size(rho,3)==1)
colormap gray
end
axis image
axis off
title('Albedo')
figure(1004)
Ndisp = N;
Ndisp(:,:,3) = -Ndisp(:,:,3);
Ndisp = 0.5*(Ndisp+1);
imagesc(Ndisp)
axis image
axis off
title('Normals')
drawnow
%%% Save to an obj file
if(do_write_obj)
disp('Saving results');
export_obj(XYZ,N,rho,mask,dataset);
end