clear all;

% uses the errors_final_smaller.mat "training" matrix
% uses the "test1" SVM model

% 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];


scene = sceneCreate('frequencyOrientation');
scene = sceneSet(scene,'fov',4);

vcAddAndSelectObject(scene);
%sceneWindow

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

vcReplaceAndSelectObject(oi);
%oiWindow;

sensor = sensorCreate;
sensor = sensorSetSizeToFOV(sensor,5);
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');
whitePt1 = imageGet(vci,'displaywhitepoint');

srgb = xyz2srgb(xyz);
imshow(srgb./255);

[finalRow finalCol temp] = size(srgb);
sceneRow = sceneGet(scene, 'rows');    
sceneCol = sceneGet(scene, 'cols');
sFactor = [96/sceneRow 96/sceneCol];        %96x96 is the size of final image minus the border
refImg = sceneInterpolate(scene, sFactor);  %Resize reference image to match final image

fiveKxyz = [0.9803 1.0000 0.8603]; % 5000K white point
refImgBorder = zeros(100,120);              %Create border
refImgBorder(3:98,13:108)= refImg.data.luminance;

whitePt2 = 100*fiveKxyz;

refBlock = zeros(12,12,3);    %With border, final image contains 96/8 = 12 blocks for row and col
finalBlock = refBlock;
dEab = zeros(8);            %8x8 to store error values
s_dEab = zeros(8);

% set whitepoint
whitePt = {whitePt1; whitePt2};

i = 1;
j = 1;

params.sampPerDeg = 12/visDeg(3);
params.filterSize = 12/visDeg(3);
        
params.deltaEversion = '2000';  
params.filters = [];
params.imageFormat = 'xyz2';
params.dimension = 2;

for r = 3:12:finalRow-13,   %Row border consists of rows 1-2 and rows 99-100 
    for c = 13:12:finalCol-23,   %Col border consists of cols 1-12 and col 109-120
        blockLum = refImgBorder(r:r+11, c:c+11);    %Divide 64x64 into 8x8 blocks
        
        maxL = max(max(blockLum));
        refBlock(:,:,1) = fiveKxyz(1).*blockLum./maxL*100;
        refBlock(:,:,2) = fiveKxyz(2).*blockLum./maxL*100;
        refBlock(:,:,3) = fiveKxyz(3).*blockLum./maxL*100;
        
        finalBlock = xyz(r:r+11, c:c+11,:);
        tmp_1 = deltaEab(finalBlock, refBlock, whitePt);
        tmp_2 = scielab(finalBlock, refBlock, whitePt, params);
        
        h = fspecial('gaussian',5);
        tmp_1 = imfilter(tmp_1, h, 0, 'same');
        tmp_2 = imfilter(tmp_2, h, 0, 'same');
        % hasler
        tmp_srgb = xyz2srgb(finalBlock);
        rg = tmp_srgb(:,:,1)-tmp_srgb(:,:,2);
        yb = 0.5*(tmp_srgb(:,:,1)+tmp_srgb(:,:,2)) - tmp_srgb(:,:,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;
        
        hasler(i,j) = M;
        
        dEab(i,j)= mean2(tmp_1) + 3*std2(tmp_1);
        s_dEab(i,j)= mean2(tmp_2) + 3*std2(tmp_2);

        j = mod(j,8)+1;
    end
    i = mod(i,8)+1;
end

% sceneGet(referenceImg, 'rows')
% sceneGet(referenceImg, 'cols')
% vcAddAndSelectObject(referenceImg);
% sceneWindow
% At distance of 24 inches, these are the results (by Erin, Sandy, and 
% Eric), where 1=colored, 0=not colored (8x8)

truth1 = [ 0  0  0  2  2  2  2  2;
           0  0  1  2  2  2  2  2;
           0  0  1  1  2  2  1  2;
           0  0  1  1  2  2  2  2;
           0  0  1  1  2  2  2  2;
           0  0  1  1  2  2  2  2;
           0  0  1  2  2  2  2  2;
           0  0  1  2  2  2  2  1 ];
       
truth2 = [ 0  0  0  1  1  1  1  1;
           0  0  0  1  1  1  1  1;
           0  0  0  0  1  1  1  1;
           0  0  0  0  1  1  1  1;
           0  0  0  1  1  1  1  1;
           0  0  0  1  1  1  1  1;
           0  0  0  1  1  1  1  1;
           0  0  0  1  1  1  1  1 ];
       
%% Classify "test" patterns       
freq = [ 1 2 3 4 5 6 7 8 ] / .5;
ang = 0:100/8:90;
dist = 24*2.54;

dEab_test = reformatData(dist, freq, ang, dEab);
s_dEab_test = reformatData(dist, freq, ang, s_dEab);
hasler_test = reformatData(dist, freq, ang, hasler);

y = classifyLR('CIELAB',dEab_test);
y_dEab_LR = [y(1:8) y(9:16) y(17:24) y(25:32) y(33:40) y(41:48) y(49:56) y(57:64)];
y = classifySVM('CIELAB',dEab_test);
y(find(y==-1)) = 0;
y_dEab_SVM = [y(1:8) y(9:16) y(17:24) y(25:32) y(33:40) y(41:48) y(49:56) y(57:64)];
y = classifyNB('CIELAB',dEab_test);
y_dEab_NB = [y(1:8) y(9:16) y(17:24) y(25:32) y(33:40) y(41:48) y(49:56) y(57:64)];

y = classifyLR('sCIELAB',s_dEab_test);
y_s_dEab_LR = [y(1:8) y(9:16) y(17:24) y(25:32) y(33:40) y(41:48) y(49:56) y(57:64)];
y = classifySVM('sCIELAB',s_dEab_test);
y(find(y==-1)) = 0;
y_s_dEab_SVM = [y(1:8) y(9:16) y(17:24) y(25:32) y(33:40) y(41:48) y(49:56) y(57:64)];
y = classifyNB('sCIELAB',s_dEab_test);
y_s_dEab_NB = [y(1:8) y(9:16) y(17:24) y(25:32) y(33:40) y(41:48) y(49:56) y(57:64)];

y = classifyLR('hasler',hasler_test);
y_hasler_LR = [y(1:8) y(9:16) y(17:24) y(25:32) y(33:40) y(41:48) y(49:56) y(57:64)];
y = classifySVM('hasler',hasler_test);
y(find(y==-1)) = 0;
y_hasler_SVM = [y(1:8) y(9:16) y(17:24) y(25:32) y(33:40) y(41:48) y(49:56) y(57:64)];
y = classifyNB('hasler',hasler_test);
y_hasler_NB = [y(1:8) y(9:16) y(17:24) y(25:32) y(33:40) y(41:48) y(49:56) y(57:64)];


error_dEab_LR = sum(sum(abs(y_dEab_LR - truth2)));
error_dEab_SVM = sum(sum(abs(y_dEab_SVM - truth2)));
error_dEab_NB = sum(sum(abs(0.5*(y_dEab_NB - truth1))));

error_s_dEab_LR = sum(sum(abs(y_s_dEab_LR - truth2)));
error_s_dEab_SVM = sum(sum(abs(y_s_dEab_SVM - truth2)));
error_s_dEab_NB = sum(sum(abs(0.5*(y_s_dEab_NB - truth1))));

error_hasler_LR = sum(sum(abs(y_hasler_LR - truth2)));
error_hasler_SVM = sum(sum(abs(y_hasler_SVM - truth2)));
error_hasler_NB = sum(sum(abs(0.5*(y_hasler_NB - truth1))));

%%

[xe xs xh truth] = prepareData('LR');
y = classifyLR('CIELAB',xe);
train_dEab_LR = sum(abs(y-truth));
y = classifyLR('sCIELAB',xs);
train_s_dEab_LR = sum(abs(y-truth));
y = classifyLR('hasler',xh);
train_hasler_LR = sum(abs(y-truth));
    
%%

[xe xs xh truth] = prepareData('SVM');
y_CIE = truth;
y = classifySVM('CIELAB',xe);
if (~isempty(find(y_CIE==0)))
    fid = fopen('models.CIELAB.transduction');
    update = textscan(fid,'%n%n%n%n','delimiter',':');
    fclose(fid);
    y_new = update{2};
    y_CIE(find(y_CIE==0)) = y_new;
    y_CIE(find(y_CIE==-1)) = 0;
else
    y_CIE(find(y_CIE==-1)) = 0;
end
y(find(y==-1))=0;
train_dEab_SVM = sum(abs(y-y_CIE));

y_sCIE = truth;
y = classifySVM('sCIELAB',xs);
if (~isempty(find(y_sCIE==0)))
    fid = fopen('models.sCIELAB.transduction');
    update = textscan(fid,'%n%n%n%n','delimiter',':');
    fclose(fid);
    y_new = update{2};
    y_sCIE(find(y_sCIE==0)) = y_new;
    y_sCIE(find(y_sCIE==-1)) = 0;
else
    y_sCIE(find(y_sCIE==-1)) = 0;
end
y(find(y==-1))=0;
train_s_dEab_SVM = sum(abs(y-y_sCIE));

y_has = truth;
y = classifySVM('hasler',xh);
if (~isempty(find(y_has==0)))
    fid = fopen('models.hasler.transduction');
    update = textscan(fid,'%n%n%n%n','delimiter',':');
    fclose(fid);
    y_new = update{2};
    y_has(find(y_has==0)) = y_new;
    y_has(find(y_has==-1)) = 0;
else
    y_has(find(y_has==-1)) = 0;
end

y(find(y==-1))=0;
train_hasler_SVM = sum(abs(y-y_has));

%%

[xe xs xh truth] = prepareData('NB');
y = classifyNB('CIELAB',xe);
train_dEab_NB = 0.5*sum(abs(y-truth));
y = classifyNB('sCIELAB',xs);
train_s_dEab_NB = 0.5*sum(abs(y-truth));
y = classifyNB('hasler',xh);
train_hasler_NB = 0.5*sum(abs(y-truth));

error_matrix = [error_dEab_LR error_dEab_SVM error_dEab_NB;
                error_s_dEab_LR error_s_dEab_SVM error_s_dEab_NB;
                error_hasler_LR error_hasler_SVM error_hasler_NB;]'./64
            
train_err_matrix = [train_dEab_LR train_dEab_SVM train_dEab_NB;
                    train_s_dEab_LR train_s_dEab_SVM train_s_dEab_NB;
                    train_hasler_LR train_hasler_SVM train_hasler_NB;]'./400
