fiveKxyz = [0.9803 1.0000 0.8603]; % 5000K white point
truth = zeros(24,30,3);
%mask = zeros(24,30);
%mask(5:21, 8:24) = 1;

% degrees per box at 6 in = 1.4917
% deg per box at 12 in = 0.7460
% deg per box at 24 in = 0.3730
% deg per box at 36 in = 0.2487

visDeg = [1.4917 0.7460 0.3730 0.2487];

%err = zeros(17,17);
% maybe store "different" errors for different distances (they should all
% be roughly the same, but at least there's some "variation" there....)

%error_sCIE = zeros(10,10);
error_CIE = zeros(10,10);
error_hasler = zeros(10,10);
whitePt1 = [];
for f=1:10,
    for a=0:10:90,
        f
        a
        p.freq = f; p.contrast = 1; p.ph = 0;
         p.ang= a*pi/180; p.row = 64; p.col = 64; p.GaborFlag=0;
         [scene,p] = sceneCreate('harmonic',p);

        scene = sceneSet(scene,'fov',0.5); % could vary FOV here

        vcAddAndSelectObject(scene);
        %sceneWindow

        oi = oiCreate('default');
        oi = oiCompute(scene,oi);

        vcReplaceAndSelectObject(oi);
        %oiWindow;

        sensor = sensorCreate;
        sensor = sensorSetSizeToFOV(sensor, 1.2);
        %sensorGet(sensor, 'size')
        sensor = sensorCompute(sensor, oi);

        vcReplaceAndSelectObject(sensor);
        %sensorImageWindow

        vci = vcImageCreate;
        vci = imageSet(vci, 'colorbalancemethod','gray world');
        vci = imageSet(vci,'colorconversionmethod','MCC optimized');
        vci = imageSet(vci,'demosaicmethod', 'adaptive laplacian');
        vci = imageSet(vci,'internalcolorspace', 'XYZ');
        vci = vcimageCompute(vci, sensor);
        vcReplaceAndSelectObject(vci);
        %vcimageWindow
        
        xyz = imageGet(vci,'dataXYZ');
        if (isempty(whitePt1)),
            whitePt1 = imageGet(vci,'displaywhitepoint');
        end
        srgb = xyz2srgb(xyz);
        
        %figure(1)
        %imshow(srgb./255);
        

        t.freq = f; t.contrast = 1; t.ph = 0;
         t.ang= a*pi/180; t.row = 12; t.col = 12; t.GaborFlag=0;
         [scene,t] = sceneCreate('harmonic',t);
        
        % create "reference" image
        maxL = max(max(scene.data.luminance));
        truth(7:18, 10:21,1) = fiveKxyz(1).*scene.data.luminance./maxL*100;         
        truth(7:18, 10:21,2) = fiveKxyz(2).*scene.data.luminance./maxL*100;
        truth(7:18, 10:21,3) = fiveKxyz(3).*scene.data.luminance./maxL*100;
        
        whitePt2 = fiveKxyz.*100;
        srgb_true = xyz2srgb(truth);
        
        %figure(2)
        %imshow(srgb_true./255);
        
        % set whitepoint
        whitePt = {whitePt1; whitePt2};
        
        % calculate CIELAB "error" image
        tmp = deltaEab(xyz(7:18, 10:21,:), truth(7:18,10:21,:), whitePt);
        % remove "noise" in "black" area
        % dEab = tmp(5:21, 8:24);
        dEab = tmp;

        
        %figure(3)
        %colormap('gray');
        %imagesc(dEab);
        %title('CIELAB error');
        
        % calculate sCIELAB "error" image
        % set params for sCIELAB
   %     params.sampPerDeg = 12/visDeg(3);
    %    params.filterSize = 12/visDeg(3);
        
        params.deltaEversion = '2000';  
        params.filters = [];
        params.imageFormat = 'xyz2';
        params.dimension = 2;       

        %tmp = scielab(xyz, truth, whitePt, params);
        %s_dEab = tmp(5:21,8:24);
        
        
        %figure(4)
        %colormap('gray');
        %imagesc(s_dEab);
        %title('sCIELAB error');
        
        % filter and "collect" data
        h = fspecial('gaussian',5);
        dEab = imfilter(dEab, h, 0, 'same');
        %s_dEab = imfilter(s_dEab, h, 0, 'same');
        
        for ind = 1:4,
            params.sampPerDeg = 12/visDeg(ind); % it's really the "degree" we see at the back of our eye
            params.filterSize = 12/visDeg(ind);
            tmp = scielab(xyz(7:18, 10:21,:), truth(7:18,10:21,:), whitePt, params);
            %tmp = scielab(xyz, truth, whitePt, params);
            % remove "noise" in "black" area
            %eval(['s_dEab_' num2str(ind) '= tmp(7:18, 10:21);']);
            eval(['s_dEab_' num2str(ind) '= tmp;']);
            eval(['s_dEab_' num2str(ind) '=imfilter(s_dEab_' num2str(ind) ',h,0,''same'');']);
            eval(['error_sCIE_' num2str(ind) '(a/10+1, f) = mean2(s_dEab_' num2str(ind) ') + 3*std2(s_dEab_' num2str(ind) ');'])
        end
        

        % Colorfulness metric (Hasler, Susstrunk)
        srgb_tmp = xyz2srgb(xyz(7:18, 10:21,:));
        rg = srgb_tmp(:,:,1)-srgb_tmp(:,:,2);
        yb = 0.5*(srgb_tmp(:,:,1)+srgb_tmp(:,:,2)) - srgb_tmp(:,:,3);
    
        stddev_rgyb = sqrt( std2(rg)^2+ std2(yb)^2 );
        mean_rgyb = sqrt( mean2(rg)^2 + mean2(yb)^2 );
        M = stddev_rgyb + 0.3*mean_rgyb;
        
        error_hasler(a/10+1,f) = M;
        
        %error_sCIE(a/10+1, f) = mean2(s_dEab) +3*std2(s_dEab);%mean2(s_dEab) + 3*std2(s_dEab);
        %error_sCIE_std(a/10+1, f) = std2(s_dEab);
        error_CIE(a/10+1, f) = mean2(dEab) + 3*std2(dEab);%mean2(dEab) + 3*std2(dEab);
        %error_CIE_std(a/10+1, f) = std2(dEab);
    end
end

filename = 'errors_final.mat';
disp(sprintf('Attempting to write data to %s...\n', filename));
if(~exist(which(filename)))
    eval(['save ' filename ' error_sCIE_1 error_sCIE_2 error_sCIE_3 error_sCIE_4 error_hasler error_CIE']);
    disp('Data written.');
else
    disp('Data not written. File already exists.');
end