iterSpEst/s_analyzeTime.m at master · hblasins/iterSpEst · 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
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
149
150
151
152
153
% This is a demo script that runs the spectral estimatino algorithm on a
% subsamples image of a Macbeth chart and compares the execution time of a
% conventional approach (CVX) and the proposed iterative method (ADMM).
%
% Copyright, Henryk Blasinski, 2015


close all;
clear all;
clc;

% Choose the number of basis functions
nBasis = 6;
averageError = 0.001;
wave = 380:5:950; wave = wave(:);

deltaL = wave(2) - wave(1);
nWaves = length(wave);


% The camera is the product of the quantum efficiency and the illuminant
fName = fullfile(iterSpEstRoot,'Parameters','sensorQe');
sensorQe = readSpectralQuantity(fName,wave);

fName = fullfile(iterSpEstRoot,'Parameters','illuminant');
illuminantEnergy = readSpectralQuantity(fName,wave);
illuminantQuanta = energy2Quanta(illuminantEnergy,wave);

camera = illuminantQuanta'*diag(sensorQe)*deltaL;

% Reference reflectance
fName = fullfile(iterSpEstRoot,'Parameters','macbethRefl');
reflRef = readSpectralQuantity(fName,wave);

% Reflectance basis functions
fName = fullfile(iterSpEstRoot,'Parameters','reflBasis');
reflBasis = readSpectralQuantity(fName,wave);
reflBasis = reflBasis(:,1:nBasis);

% Load the sample macbeth Image
fName = fullfile(iterSpEstRoot,'Data','macbethImage');
load(fName);
h = size(Img,1);
w = size(Img,2);
nChannels = size(Img,3);

Img = max(Img - repmat(shiftdim(cameraOffset',-2),[h w 1]),0);
cameraMat = diag(cameraGain)*camera;

reflRefImg = imresize(reshape(reflRef',[4 6 nWaves]),[h w],'nearest');


figure;
imshow(Img(:,:,5),[]);
title('Sample image channel');


%% Solve spectral smoothness constraint only
%  This is a least-squares problem that can be solved quite efficiently
%  with Matlab functions.
alpha = 0.1;


tic
XlsCvx = spectralEstCVX(Img, cameraMat, reflBasis, alpha, 0, 0);
tlsCvx = toc;


tic;
Rlambda = [eye(nWaves-1) zeros(nWaves-1,1)] - [zeros(nWaves-1,1) eye(nWaves-1)];
measVals = reshape(Img,h*w,nChannels)';

mhat = [measVals; zeros(nWaves-1,size(measVals,2))];
Aalpha = [cameraMat*reflBasis; sqrt(alpha)*Rlambda*reflBasis];
Xls = Aalpha\mhat;
Xls = reshape(Xls',[h w nBasis]);
tls = toc;


%% Spectral smoothness + box constraint
alpha = 0.1;
beta = 0;
gamma = 1;

tic
XBoxCvx = spectralEstCVX(Img, cameraMat, reflBasis, alpha, beta, gamma);
tBoxCvx = toc;

tstart = tic;
[ XBoxAdmm, histBox ] = spectralEstADMM( Img, cameraMat, reflBasis, alpha, beta, gamma,...
    'maxIter',1000,...
    'verbose',true,...
    'reference',XBoxCvx,...
    'tolConv',averageError*(h*w));
tBoxAdmm = toc(tstart);

showResults( XBoxCvx, reflBasis, cameraMat, Img, wave, reflRefImg );

%% Spectral and spatial smoothness
alpha = 0.1;
beta = 0.1;
gamma = 0;

tic
XSpatialCvx = spectralEstCVX(Img, cameraMat, reflBasis, alpha, beta, gamma);
tSpatialCvx = toc;


tstart = tic;
[ XSpatialAdmm, histSpatial ] = spectralEstADMM( Img, cameraMat, reflBasis, alpha, beta, gamma,...
    'maxIter',1000,...
    'verbose',true,...
    'reference',XSpatialCvx,...
    'tolConv',averageError*(h*w));
tSpatialAdmm = toc(tstart);

showResults( XSpatialAdmm, reflBasis, cameraMat, Img, wave, reflRefImg );


%% Spectral and spatial smoothness + box
alpha = 0.1;
beta = 0.1;
gamma = 1;

tic
XBoxSpatialCvx = spectralEstCVX(Img, cameraMat, reflBasis, alpha, beta, gamma);
tBoxSpatialCvx = toc;

tstart = tic;
[ XBoxSpatialAdmm, histBoxSpatial ] = spectralEstADMM( Img, cameraMat, reflBasis, alpha, beta, gamma,...
    'maxIter',1000,...
    'verbose',true,...
    'reference',XSpatialCvx,...
    'tolConv',averageError*(h*w));
tBoxSpatialAdmm = toc(tstart);

% Timing results
fprintf('Spectrally smooth (least-squares):            %f sec.\n',tls);
fprintf('Spectrally smooth + box (ADMM):               %f sec\n',tBoxAdmm);
fprintf('Spectrally and spatially smooth (ADMM):       %f sec\n',tSpatialAdmm);
fprintf('Spectrally and spatially smooth + Box (ADMM): %f sec\n',tBoxSpatialAdmm);

dName = fullfile(iterSpEstRoot,'Results');
if ~exist(dName,'dir');
    mkdir(dName);
end

fName = fullfile(iterSpEstRoot,'Results','macbethResults.mat');
save(fName,'tBoxAdmm','histBox','tSpatialAdmm','histSpatial','tBoxSpatialAdmm','histBoxSpatial',...
                'tBoxCvx','tSpatialCvx','tBoxSpatialCvx','tlsCvx','tls');