%% Gaussain Fitting for single LED Point Spectral Function
% The output x of this script is used in brsGetBacklightImage2 as PSF
% curve fitting

nLED = 1380;
gaussrad = [4 6 8 12 16 32 64 128];
gausswidth = min(89*3,10*gaussrad); % Filter constants
led_coefficients = [0;103.195;6.464;44.012;78.929;71.404;0;64.792;]; 
% relative weights of the test gaussians derived using brsPSF fscript
% and measured PSF data

%% Center LED response
backlight2 = zeros(31*5,89*3);
backlight2(16*5,44*3) = 1;
backlight2_sum = zeros(31*5,89*3);
for idx = 1:length(gaussrad)
    backlight2_gauss(:,:,idx) = led_coefficients(idx)*imfilter(backlight2,fspecial('gaussian',gausswidth(idx),gaussrad(idx)),'replicate');
    backlight2_sum = backlight2_sum + backlight2_gauss(:,:,idx);
end
% Resize backlight to full resolution
backlight2 = imresize(backlight2_sum,[31*2 89*2],'bicubic');
backlight2 = backlight2./max(backlight2(:));
%figure;imagesc(backlight2);colormap('gray');
%% Get Each LED responses
LED_response = zeros(31*2*89*2,nLED);
idx=1;
for m = 1:31
    for n = 1:44+2*(m/2-floor(m/2))
        col = 2*n-2*(m/2-floor(m/2));
        Test_LED = zeros(31*2,89*2);
        Test_LED(m*2,col*2) = 1;
        backlight = imfilter(backlight2,Test_LED,'conv');
        %figure;imagesc(backlight);colormap('gray');
        LED_response(:,idx) = backlight(:);
        idx = idx+1;
    end
end
save LED_response_matrix_2_2 LED_response

%% Display LED response in intermediate space, and Evaluate
load hdrCalibrationData\LED_response_matrix_2_2;
figure(); imagesc(reshape(LED_response(:,1380),31*2,89*2)); colormap('gray');
TEST_LED_value = x;
est_backlight = reshape(LED_response*TEST_LED_value,31*2,89*2);
figure(); imagesc(ideal_backlight); colormap('gray');
figure(); imagesc(est_backlight); colormap('gray');
figure(); plot(ideal_backlight(:));
figure(); plot(est_backlight(:));
figure(); plot(abs(ideal_backlight(:)-est_backlight(:)));

idx=1;
for m = 1:31
    for n =1:89
        if (mod(m+n,2)==0)
            x(idx) = hexagonal_backlight(m,n);
            idx = idx+1;
        end
    end
end
