s_evalGlobalBoothLEDCube.m

% This script runs the computational spectral flash algorithm on sets of
% images captured in the LED Cube booth. The booth LEDs are controlled to
% mimic the flash, and the ambient illuminant is generated by an external
% light bulb.
%
% Note that parts of the code in this script refer to image acquisition
% hardware that would capture images with esitmated flash weights. 
% By default these sections are not executed (capturePhotos =
% false), and instead data that had been captured is read in.
%
% Copyright, Henryk Blasinski 2017

close all;
clear all;
clc;

ieInit;
set(groot,'defaultAxesColorOrder',[1 0 0; 0 1 0; 0 0 1; 1 1 0; 1 0 1]);

%%

wave = 400:5:700;
nWaves = length(wave);



desiredIll = 'D65';

capturePhotos = false;

% Save plots and figures if the path is not empty.
destDir = fullfile(cmfRootPath,'..','Figures');
% resDir = [];


% Create base sensor model
fName = fullfile(cmfRootPath,'Parameters','XimeaSpectralResponsivities');
cameraResp = ieReadColorFilter(wave,fName);

sensorBase = sensorCreate('bayer (bggr)');
sensorBase = sensorSet(sensorBase,'wave',wave);
sensorBase = sensorSet(sensorBase,'size',[1024 1280]);
sensorBase = sensorSet(sensorBase,'noise flag',0);
sensorBase = sensorSet(sensorBase,'filter transmissivities',cameraResp);

% Load LEDCube led spectra. The LED cube contains a few wide-band LEDs, 
% these are removed from the flash design. We pick only the narrowband ones.
selLEDs = [2 4 6 7 8 9 10 11];
nLEDs = length(selLEDs);
fName = fullfile(cmfRootPath,'Parameters','LEDCubeSpectra.mat');
ledCubeLeds = ieReadSpectra(fName,wave);
ledCubeLeds = ledCubeLeds(:,selLEDs);
ledCubeLeds = Energy2Quanta(wave,ledCubeLeds);
ledCubeLeds = ledCubeLeds/max(ledCubeLeds(:));


ambientDeltaE = zeros(24,5);
matchingDeltaE = zeros(24,5);
complementDeltaE = zeros(24,5);

for jj=3
    
    switch jj
        case 1
            illCondition = 'Tungsten-LowPower';
            dataDir = fullfile(cmfRootPath,'..','Data','Booth',sprintf('2016-11-29 15-14-15-%s',illCondition));
            cp = [732 958;987 955;985 787;733 789];
            EV = 20;
        
        case 2
            illCondition = 'EcoBulb-LowPower';
            dataDir = fullfile(cmfRootPath,'..','Data','Booth',sprintf('2016-12-06 15-31-44-%s',illCondition));
            cp = [732 958;987 955;985 787;733 789];
            EV = 20;
        
        case 3
            illCondition = 'Green-LowPower';
            dataDir = fullfile(cmfRootPath,'..','Data','Booth',sprintf('2016-11-29 14-35-30-%s',illCondition));
            cp = [732 958;987 955;985 787;733 789];
            EV = 20;       
            
        case 4
            illCondition = 'Orange-LowPower';
            dataDir = fullfile(cmfRootPath,'..','Data','Booth',sprintf('2016-11-29 14-49-57-%s',illCondition));
            cp = [732 958;987 955;985 787;733 789];
            EV = 20;
            
        
    end
    
    % Load LED power setting
    load(fullfile(dataDir,'maxPower'));
    
    % Load ambient illuminant
    ambientSpectrum = ieReadSpectra(fullfile(dataDir,'ambient.mat'),wave);
    ambientSpectrum = Energy2Quanta(wave,ambientSpectrum);
    ambientSpectrum = ambientSpectrum/max(ambientSpectrum);
   
    % Load captured data    
    [ measurement, mask, cp ] = readXimeaImageStack(dataDir,EV,nLEDs,'cp',cp);
    
     
    %% Ambient estimate
    
    % Global: all image pixels
    [ ambientEst, ambientWghts, ambientPredictions ] = globalAmbientEst( measurement.downsampled.ambient,...
        measurement.downsampled.led, ...
        ledCubeLeds,...
        'alpha',1e2);
    
    % Patch: Macbeth patches only
    [ ambientEstPatch, ambientWghtsPatch, ambientPredictionsPatch ] = globalAmbientEst( measurement.patch.ambient,...
        measurement.patch.led, ...
        ledCubeLeds,...
        'alpha',1);
    
    % Plot the quality of the approximation in camera RGB space
    
    figure;
    hold on; grid on; box on;
    plot(squeeze(measurement.patch.ambient)',ambientPredictionsPatch','o');
    xlabel('Ambient appearance');
    ylabel('Approximated');
    title('RGB space');
    
    
    %% Render ambient approximation image
    
    [~, ip] = renderFlashImage(zeros(size(measurement.raw.ambient)),measurement.raw.led,ambientWghts,...
                                    sensorBase,'name','Ambient approximation');
    
    
    image = ipGet(ip,'data srgb');
    figure; imshow(image,'Border','tight');
    
    if isempty(resDir) == false
        fName = fullfile(resDir,sprintf('%s_ambientApprox.eps',illCondition));
        print('-depsc',fName);
    end
    
    
    %% Complement estimate
    
    ill = illuminantCreate(desiredIll,wave);
    desiredSpectrum = illuminantGet(ill,'photons');
    
    [ flashEst, flashWghts ] = globalComplementEst( desiredSpectrum, ambientEst, ledCubeLeds, cameraResp,...
        'flashMode',true);
    
    
    %% Best approximation to the ambient
    
    cvx_begin
        variable optWghts(nLEDs,1)
        minimize norm(ledCubeLeds*optWghts - ambientSpectrum*max(ambientEst),2)
        subject to
            ledCubeLeds*optWghts >= 0
    cvx_end
    
    
    % Plot the illuminant spectra
    figure;
    hold on; grid on; box on;
    plot(wave,ambientEst,'LineWidth',2);
    plot(wave,ambientSpectrum*max(ambientEst),'--','LineWidth',2);
    plot(wave,ledCubeLeds*optWghts,'-.','LineWidth',2);
    plot(wave,ledCubeLeds*flashWghts,'LineWidth',2);
    xlabel('Wavelength, nm','Interpreter','Latex','FontSize',6);
    set(gca,'TickLabelInterpreter','Latex');
    mx = max([ambientEst(:); ambientSpectrum*max(ambientEst); ledCubeLeds*flashWghts]);
    ylim([-0.05*mx 1.05*mx]);
    set(gcf,'Units','Centimeters');
    set(gca,'FontSize',6);
    set(gcf,'PaperPosition',[1 1 4 3.25]);
    
    if isempty(resDir) == false
        fName = fullfile(resDir,sprintf('%s_spectra.eps',illCondition));
        print('-depsc',fName);
    end
        
    
    %% Complementary flash
    % Render an image as if captured with the complementary flash.
    
    [sensor, ip] = renderFlashImage(measurement.raw.ambient,measurement.raw.led,flashWghts,...
                                    sensorBase,'name',sprintf('Computational-complement : %s',desiredIll));
    
    rendering = sensorGet(sensor,'volts');                           
   
    [complementLab, ~, complementDeltaE(:,jj)] = macbethColorError(ip,desiredIll,cp);
    
    image = ipGet(ip,'data srgb');
    figure; imshow(image,'Border','tight');
    
    if isempty(resDir) == false
        fName = fullfile(resDir,sprintf('%s_flashComplement.eps',illCondition));
        print('-depsc',fName);
    end
    
    
    if capturePhotos
        % Control the hardware (camera+LEDCube) to acquire the image with
        % the pre-computed weights. 
        pwr = zeros(11,1);
        pwr(selLEDs) = flashWghts;
        
        [Img, sh, g, bl] = getXimeaLEDCube(EV,pwr,'',maxPower);
        imwrite(Img,fullfile(dataDir,sprintf('Complement_%i.tiff',EV)));
        fid = fopen(fullfile(dataDir,sprintf('Complement_%i.txt',EV)),'w');
        fprintf(fid,'Shutter[us] Gain[dB] BlackLevel\n');
        fprintf(fid,'%i %f %i\n',sh,g,bl);
        fclose(fid);
    end
    
    % Read in the image and compare
    
    complementCapture = getImage(fullfile(dataDir,sprintf('Complement_%i.tiff',EV)));
    complementCapture = imageExpose(complementCapture,0.01);
    sensor = sensorSet(sensorBase,'volts',complementCapture);
    sensor = sensorSet(sensor,'name',sprintf('Captured-complement : %s',desiredIll));
    
    ieAddObject(sensor);
    sensorWindow;
    
    
    ip = ipCreate;
    ip = ipCompute(ip,sensor);
    ip = ipSet(ip,'name','Captured-complement');
    ieAddObject(ip);
    ipWindow;
    
    image = ipGet(ip,'data srgb');
    figure; imshow(image,'Border','tight');
    
    if isempty(resDir) == false
        fName = fullfile(resDir,sprintf('%s_capturedComplement.eps',illCondition));
        print('-depsc',fName);
    end
    
    
    % Create linearity plot
    ds = imresize(rendering,0.125,'nearest');
    ds2 = imresize(complementCapture,0.125,'nearest');
    
    figure;
    hold on; grid on; box on;
    plot(ds(:),ds2(:),'k.');
    xlabel('Computational','Interpreter','Latex','FontSize',6);
    ylabel('Captured','Interpreter','Latex','FontSize',6);
    xlim([0 1]);
    ylim([0 1]);
    set(gca,'TickLabelInterpreter','Latex');
    set(gcf,'Units','Centimeters');
    set(gca,'FontSize',6);
    set(gcf,'PaperPosition',[1 1 4 4]);
    
    if isempty(resDir) == false
        fName = fullfile(resDir,sprintf('%s_linearity.eps',illCondition));
        print('-depsc',fName);
    end
    
    
    
    %% Matching flash
    
    matchingAmbWghts = ambientWghts./max(ambientWghts);
    
    [sensor, ip] = renderFlashImage(measurement.raw.ambient,measurement.raw.led,matchingAmbWghts,...
                                    sensorBase,'name',sprintf('Computational-matching : %s',desiredIll));
    
    [matchingLab, ~, matchingDeltaE(:,jj)] = macbethColorError(ip,desiredIll,cp);
    
    
    image = ipGet(ip,'data srgb');
    figure; imshow(image,'Border','tight');
    
    if isempty(resDir) == false
        fName = fullfile(resDir,sprintf('%s_flashMatch.eps',illCondition));
        print('-depsc',fName);
    end
    
    if capturePhotos
        % Control the hardware (camera+LEDCube) to acquire the image with
        % the pre-computed weights. 
        cvx_begin
            variables mW(nLEDs,1) w
            minimize norm(ambientEst - ledCubeLeds*mW)
            subject to
            mW >= 0
        cvx_end
        
        pwr = zeros(11,1);
        pwr(selLEDs) = mW;
        
        [Img, sh, g, bl] = getXimeaLEDCube(EV,pwr,'rescale',maxPower);
        imwrite(Img,fullfile(dataDir,sprintf('Matching_%i.tiff',EV)));
        fid = fopen(fullfile(dataDir,sprintf('Matching_%i.txt',EV)),'w');
        fprintf(fid,'Shutter[us] Gain[dB] BlackLevel\n');
        fprintf(fid,'%i %f %i\n',sh,g,bl);
        fclose(fid);
    end
    
    % Display the matching image
    
    complementCapture = getImage(fullfile(dataDir,sprintf('Matching_%i.tiff',EV)));
    complementCapture = imageExpose(complementCapture,0.01);
    sensor = sensorSet(sensor,'volts',complementCapture);
    sensor = sensorSet(sensor,'name',sprintf('Captured-matching : %s',desiredIll));
    
    ieAddObject(sensor);
    sensorWindow;
    
    
    ip = ipCreate;
    ip = ipCompute(ip,sensor);
    ip = ipSet(ip,'name','Captured-matching');
    ieAddObject(ip);
    
    ipWindow;
    
    
    
    %% Ambient baseline
    
    rendering = measurement.raw.ambient;
    rendering = imageExpose(rendering,0.01);
    
    sensor = sensorSet(sensorBase,'volts',rendering);
    sensor = sensorSet(sensor,'name',sprintf('Ambient'));
    
    ieAddObject(sensor);
    sensorWindow;
    
    
    
    ip = ipCreate;
    % ip = ipSet(ip,'sensor conversion method','none');
    ip = ipSet(ip,'name','Ambient');
    ip = ipCompute(ip,sensor);
    ieAddObject(ip);
    
    ipWindow;
    
    [ambientLab, ~, ambientDeltaE(:,jj)] = macbethColorError(ip,desiredIll,cp);
    
    
    image = ipGet(ip,'data srgb');
    figure; imshow(image,'Border','tight');
    
    if isempty(resDir) == false
        fName = fullfile(resDir,sprintf('%s_ambient.eps',illCondition));
        print('-depsc',fName);
    end
    
    
end

%% D65 baseline


if capturePhotos
    % Control the hardware (camera+LEDCube) to acquire the image with
    % the pre-computed weights. 
    
    d65 = illuminantCreate('d65',wave);
    d65 = illuminantGet(d65,'photons');
    d65 = d65/max(d65);
    
    cvx_begin
        variables mW(nLEDs,1) w
        minimize norm(d65 - ledCubeLeds*mW)
        subject to
            mW >= 0
    cvx_end
    
    
    
    pwr = zeros(11,1);
    pwr(selLEDs) = mW;
    
    [Img, sh, g, bl] = getXimeaLEDCube(EV,pwr,'rescale',maxPower);
    imwrite(Img,fullfile(dataDir,sprintf('D65_%i.tiff',EV)));
    fid = fopen(fullfile(dataDir,sprintf('D65_%i.txt',EV)),'w');
    fprintf(fid,'Shutter[us] Gain[dB] BlackLevel\n');
    fprintf(fid,'%i %f %i\n',sh,g,bl);
    fclose(fid);
    
    
    TL_Cube_awakeMode('com3',[0;0;0;0],1,2);
    TL_Cube_lightRecipe('com3',[0, 0, 0, 0]',zeros(11,1),2);
    
end

%% Read D65 data
dataDir = fullfile(cmfRootPath,'..','Data','Booth','D65');
complementCapture = getImage(fullfile(dataDir,sprintf('D65_%i.tiff',EV)));
complementCapture = imageExpose(complementCapture,0.01);
sensor = sensorSet(sensor,'volts',complementCapture);
sensor = sensorSet(sensor,'name',sprintf('D65 : %s',desiredIll));

ieAddObject(sensor);
sensorWindow;

ip = ipCreate;
ip = ipCompute(ip,sensor);
ip = ipSet(ip,'name','D65');
ieAddObject(ip);

ipWindow;

[d65Lab, ~, d65DeltaE] = macbethColorError(ip,desiredIll,cp);

    

%% Genrate plots

data = [mean(ambientDeltaE); mean(matchingDeltaE); mean(complementDeltaE)]';

figure;
hold on; grid on; box on;
bar(data)
plot(0:5,mean(d65DeltaE)*ones(6,1),'r','LineWidth',2);
ylabel('Average CIE $\Delta E$','Interpreter','Latex','FontSize',8);
xlabel('Scene','Interpreter','latex','FontSize',8);
xlim([0.5 4.5]);
set(gca,'XTick',1:4);
set(gca,'XTickLabel',{'Tungsten','Fluorescent','Orange','Green'})
set(gca,'TickLabelInterpreter','Latex');
set(gcf,'Units','Centimeters');
set(gca,'FontSize',6);
set(gcf,'PaperPosition',[1 1 14 5]);
legend({'Ambient','Match','Complement'},'Location','NorthWest','Interpreter','latex')

if isempty(resDir) == false
    fName = fullfile(resDir,'CIEdeltaE.eps');
    print('-depsc',fName);
end