pseudo-array/main.m at master · photon-efficient-imaging/pseudo-array · GitHub

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
%%
% Sample code and data for
%   Photon-efficient computational 3D and reflectivity
%   imaging with single-photon detectors
% by D.Shin, A.Kirmani, V.Goyal, and J.H.Shapiro, IEEE TCI 2015

% Uses variants of SPIRAL-TAP software developed in the paper
%   This is SPIRAL-TAP: Sparse Poisson intensity reconstruction algorithms?theory and practice
% by Z.Harmany, R.Marcia, and R.Willett, IEEE TIP 2012

% When using this software in your publications, please cite both papers..

%%

clc; clear; close all;


load('data_face');
[nr,nc] = size(Ts);
C0 = zeros(nr,nc);
D0 = zeros(nr,nc);
for i=1:nr
    for j=1:nc
        ts = Ts{i,j};
        C0(i,j) = length(ts);
        D0(i,j) = mean(ts);
    end
end

addpath(genpath([pwd '/fcns']));
AT = @(x) x; A = @(x) x;
C1 = SPIRALTAP_BIN(C0,A,0.5,100, ...
    'noisetype','binomial', ... % approx for binomial
    'penalty','tv', ...
    'maxiter',60,... % `120
    'Initialization',C0,...
    'AT',AT,...
    'monotone',1,...
    'miniter',1,...
    'stopcriterion',3,...
    'tolerance',1e-12,...
    'alphainit',0.1,...
    'alphaaccept',1e80,...
    'logepsilon',1e-10,...
    'saveobjective',1,...
    'savereconerror',1,...
    'savecputime',1,...
    'savesolutionpath',0,...
    'truth',C0,...
    'verbose',0);
max_count = 20;
scales = C1/max_count;
scales(scales>1) = 1;
img_filt = get_rom(Ts);
D_mean = get_thres(Ts,img_filt,scales);
D_mean(isnan(D_mean)) = 0;
AT = @(x) x; A = @(x) x;
D1 = SPIRALTAP_blank(D_mean,A,10, ...
    'noisetype','gaussian', ...
    'penalty','tv', ...
    'maxiter',40,...
    'Initialization',img_filt,...
    'AT',AT,...
    'monotone',1,...
    'miniter',1,...
    'stopcriterion',3,...
    'tolerance',1e-12,...
    'alphainit',0.01,...
    'alphaaccept',1e80,...
    'logepsilon',1e-10,...
    'saveobjective',1,...
    'savereconerror',1,...
    'savecputime',1,...
    'savesolutionpath',0,...
    'truth',zeros(nr,nc),...
    'verbose',0);
range_C = [0,max_count];
range_D = [3550,3700];

figure;
subplot(221); imagesc(C0,range_C); axis image;
title('reflectivity - ML');
subplot(222); imagesc(D0,range_D); axis image;
title('depth - ML');
subplot(223); imagesc(C1,range_C); axis image;
title('reflectivity - proposed');
subplot(224); imagesc(D1,range_D); axis image;
title('depth - proposed');
colormap('gray');