compMultispFlash/s_evalNatural100NumPatches.m at master · hblasins/compMultispFlash · 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
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
% Evaluate the illuminant estimation algorithm using pixel data from
% different number of distinct patches of uniform reflectance. To simplify
% the analysis we pick the illuminant spectrum in terms of CIE xy
% chromaticity, a camera model, and regularization.
%
% Copyright, Henryk Blaisnski 2017


close all;
clear all;
clc;

%%
fName = fullfile(cmfRootPath,'Results','simNatural100V2.mat');
load(fName);

% Pick number of LEDs
ledSets = {[1,2,3,4,5,6,7]};
nLedSets = length(ledSets);

% Illuminant chromaticity
xChromaticities = [3];
yChromaticities = [3];

% Camera models
cameraIDs = [3];

% Number of bootstrap samples
nDraws = 10;

% Number of patches in an evaluation experiment
numPatches = [1:2:10, 20, 30];
nPatchSets = length(numPatches);

alphaVec = [0.06];
nAlpha = length(alphaVec);


%% Estimate the illuminant
%  From measurements

channelEstLin =  cell(nX,nY,nCameras,nLedSets,nAlpha,nDraws,nPatchSets);
wghtEst = cell(nX,nY,nCameras,nLedSets,nAlpha,nDraws,nPatchSets);
selectedPatches = cell(nDraws,nPatchSets);

% Random draw
for d=1:nDraws
    % Patch set
    for p=1:nPatchSets

        nPatches = numPatches(p);
        patchIndices = randi(110,[nPatches, 1]);

        selectedPatches{d,p} = patchIndices;

        %LED set
        for s=1:nLedSets

            ledIndices = ledSets{s};
            nLEDs = length(ledIndices);

            currentFlash = flashNorm(:,ledIndices);

            % Chromaticity (we only analyze a subset)
            for xx=xChromaticities
                for yy=yChromaticities

                    [ spd ] = xy2Spectrum( xVec(xx), yVec(yy), wave );
                    if sum(isnan(spd)) > 0
                        continue;
                    end

                    % Regularization
                    for a=1:nAlpha
                        % Camera model
                        for cam=cameraIDs


                             measurementSet = measurement{xx,yy,cam};
                             measurementSet = [measurementSet.patch];
                             ambient = measurementSet.ambient(:,:,patchIndices);

                             measurementSet = [measurementSet.led];
                             measurementSet = measurementSet(:,ledIndices,patchIndices);

                            [ ambientEst, ambientWghts, ambientPredictions ] = globalAmbientEst( ambient, measurementSet, currentFlash, 'alpha', alphaVec(a) );


                            channelEstLin{xx,yy,cam,s,a,d,p} = ambientPredictions;
                            wghtEst{xx,yy,cam,s,a,d,p} = ambientWghts;

                        end
                    end

                end
            end
        end

    end
end

% Save data
fName = fullfile(cmfRootPath,'Results','numPatches.mat');
save(fName);


%% Compute estimation errors


estimationError = zeros(nX,nY,cam,s,a,d,p);
referenceError = zeros(nX,nY,nLedSets,1);

% Random draw
for d=1:nDraws
    % Patch set
    for p=1:nPatchSets

        patchIndices = selectedPatches{d,p};

        %LED set
        for s=1:nLedSets

            ledIndices = ledSets{s};
            nLEDs = length(ledIndices);

            % Chromaticity (we only analyze a subset)
            for xx=xChromaticities
                for yy=yChromaticities

                    [ spd ] = xy2Spectrum( xVec(xx), yVec(yy), wave );
                    if sum(isnan(spd)) > 0
                        continue;
                    end

                    cvx_begin
                        variables vv(nLEDs,1)
                        minimize norm(spd - flashNorm(:,ledIndices)*vv)
                        subject to
                            flashNorm(:,ledIndices)*vv >= 0
                    cvx_end

                    fNorm = flashNorm(:,ledIndices)*vv;
                    fNorm = fNorm/max(fNorm);

                    cvx_begin
                        variable g(1,1)
                        minimize norm(fNorm - g*spd)
                    cvx_end

                    referenceError(xx,yy,s) = rms(g*spd - fNorm);

                    % Regularization
                    for a=1:nAlpha
                        % Camera model
                        for cam=cameraIDs

                            estimatedSpectrum = flashNorm(:,ledIndices)*wghtEst{xx,yy,cam,s,a,d,p};
                            estimatedSpectrum = estimatedSpectrum/max(estimatedSpectrum);

                            cvx_begin
                                variable g(1,1)
                                minimize norm(estimatedSpectrum - g*spd)
                            cvx_end

                            err = rms(estimatedSpectrum - g*spd);

                            estimationError(xx,yy,cam,s,a,d,p) = err;

                        end
                    end

                end
            end
        end

    end
end

% Save data
fName = fullfile(cmfRootPath,'Results','numPatches.mat');
save(fName);