%% read the data

%load training data
% x_train = h5read('camelyonpatch_level_2_split_train_x.h5', '/x');
% y_train = h5read('camelyonpatch_level_2_split_train_y.h5', '/y');

%load testing data
x_test = h5read('camelyonpatch_level_2_split_test_x.h5', '/x');
y_test = double(squeeze(h5read('camelyonpatch_level_2_split_test_y.h5', '/y')));

% %example to read an image
% id = 100;%id number of the data
% im = permute(x_test(:,:,:,id),[2 3 1]); %read the image
% label = y_test(id); %label for the image
% imshow(im)
% title(mat2str(label))

xx = (zeros(size(x_test,2)*size(x_test,3)*size(x_test,4),3));

temp = 1;

for k = 1:size(x_test,4)
   V = permute(x_test(:,:,:,k),[2,3,1]);
   temp2 =V(:,:,1); 
   xx(temp:size(temp2,1)*size(temp2,1)+temp-1,1)= temp2(:);
   temp2 =V(:,:,2); 
   xx(temp:size(temp2,1)*size(temp2,1)+temp-1,2)= temp2(:);
   temp2 =V(:,:,3); 
   xx(temp:size(temp2,1)*size(temp2,1)+temp-1,3)= temp2(:);
   
   temp = temp+size(temp2,1)*size(temp2,1);
end
 
ScaleF = .01;
xx = xx(randperm(size(xx,1),round(size(xx,1)*ScaleF)),:);
[~,C]= kmeans(xx,3,'Display','iter');
clear xx
%% example training
si = size(x_test);
trainD = zeros(si(4),25); %preallocate space
avgHSV = zeros(si(4),3);

%% populate data sets
for k = 1 : 100%si(4)
    im = permute(x_test(:,:,:,k),[2 3 1]);          %read the image
    imRedChannel = double(im(:,:,1));
    imGreenChannel = double(im(:,:,2));
    imBlueChannel = double(im(:,:,3));
    
    %% get mean for each color channel
    trainD(k,1) = mean(imRedChannel(:));            %mean of R channel
    trainD(k,2) = mean(imGreenChannel(:));          %mean of G channel
    trainD(k,3) = mean(imBlueChannel(:));           %mean of B channel
    
    %% get kurtosis for each color channel
    trainD(k,4) = kurtosis(imRedChannel(:));        %kurtosis of R channel
    trainD(k,5) = kurtosis(imGreenChannel(:));      %kurtosis of G channel
    trainD(k,6) = kurtosis(imBlueChannel(:));       %kurtosis of B channel
    
    %% get skewness for each color channel
    trainD(k,7) = skewness(imRedChannel(:));        %skewness of R channel
    trainD(k,8) = skewness(imGreenChannel(:));      %skewness of G channel
    trainD(k,9) = skewness(imBlueChannel(:));       %skewness of B channel
    
    %% get variance for each color channel
    trainD(k,7) = var(imRedChannel(:));             %variance of R channel
    trainD(k,8) = var(imGreenChannel(:));           %variance of G channel
    trainD(k,9) = var(imBlueChannel(:));            %variance of B channel
    
    %% get entropy for each color channel
     trainD(k,10) = entropy(imRedChannel(:));             %entropy of R channel
     trainD(k,11) = entropy(imGreenChannel(:));           %entropy of G channel
     trainD(k,12) = entropy(imBlueChannel(:));            %entropy of B channel
    
    %% get mode for each color channel
     trainD(k,13) = mode(imRedChannel(:));             %entropy of R channel
     trainD(k,14) = mode(imGreenChannel(:));           %entropy of G channel
     trainD(k,15) = mode(imBlueChannel(:));            %entropy of B channel
    
    %% get median for each color channel
    trainD(k,13) = median(imRedChannel(:));             %entropy of R channel
    trainD(k,14) = median(imGreenChannel(:));           %entropy of G channel
    trainD(k,15) = median(imBlueChannel(:));            %entropy of B channel
    
    %% get standard deviation for each color channel
    trainD(k,16) = std(imRedChannel(:));             %entropy of R channel
    trainD(k,17) = std(imGreenChannel(:));           %entropy of G channel
    trainD(k,18) = std(imBlueChannel(:));            %entropy of B channel
    
    % convert image to double precision
    tempim = rgb2gray(im);
    A = im2double(tempim(:));
    
    % define custom functions (add functions here)
    entFunc = @(x)entropy(x(:));
    stdFunc = @(x)std(x(:));
    meanFunc = @(x)mean(x(:));
    kurtosisFunc = @(x)kurtosis(x(:));
    
    % apply nlfilter and add to trainD 
%      B = nlfilter(A, [3, 3], entFunc);
%      trainD(k, 19) = entropy(B);
%      
%      B = nlfilter(A, [5, 5], stdFunc);
%      trainD(k, 20) = std(B);
%      
%      B = nlfilter(A, [7, 7], meanFunc);
%      trainD(k, 21) = mean(B);
%      
%      B = nlfilter(A, [7, 7], kurtosisFunc);
%      trainD(k, 22) = kurtosis(B);

     imFlat = zeros(size(im,1)*size(im,2),3);
     imFlat(:,1) = reshape(im(:,:,1),size(im,1)*size(im,2),1);
     imFlat(:,2) = reshape(im(:,:,2),size(im,1)*size(im,2),1);
     imFlat(:,3) = reshape(im(:,:,3),size(im,1)*size(im,2),1);
     
     idx = predictKmeans(imFlat,C);
     a = sum(idx==1);
     b = sum(idx==2);
     c = sum(idx==3);
      
     trainD(k, 19) = a/b;
     trainD(k, 20) = a/c;
     trainD(k, 21) = b/c;
      
    
    %% kmeans segmentation
%     cform = makecform('srgb2lab');
%     lab_he = applycform(im,cform);
%     ab = double(lab_he(:,:,2:3));
%     nrows = size(ab,1);
%     ncols = size(ab,2);
%     ab = reshape(ab,nrows*ncols,2);
% 
%     nColors = 5;
%     % repeat the clustering 3 times to avoid local minima
%     [cluster_idx, cluster_center] = kmeans(ab,nColors,'distance','sqEuclidean','Replicates',3);
%     
%     % show images after kmeans applied
%     pixel_labels = reshape(cluster_idx,nrows,ncols);
%     subplot(1,2,1)
%     imshow(pixel_labels,[]), title('image labeled by cluster index');
%     
%     subplot(1,2,2)
%     imshow(im), title('image labeled by cluster index');
    
    %% get region properties
%     level = graythresh(im);
%     BW = imbinarize(im, level);
%     imshowpair(im, BW, 'montage')
end

% YIQ = rgb2ntsc(trainD);                                         %average YIQ values from average rgb
% LAB = rgb2lab(trainD);                                          %average Lab values from average rgb
% YCBCR = rgb2ycbcr(trainD);                                      %average YCBCR values from average rgb

%% clear original data before training saves memory
clear x_test
i = 1;
xlearn = 0;
yACC = 0;
numLearners = 100;
xFPR = 0;
yTPR = 0;

while numLearners <= 1000
    %train model
    Mdl = fitensemble(trainD,y_test,'AdaBoostM1',numLearners,'Tree','NPrint',5);

    %evaluate model
    Yfit = predict(Mdl,trainD);

    % confusion matrix (Mean)
    C = confusionmat(y_test, Yfit)
    TP = C(1);
    FN = C(2);
    FP = C(3);
    TN = C(4);
    
    % TPR
    TPR = TP/(TP+FN);
    FPR = FP/(FP+TP);
    ACC = (TP+TN)/(TP+TN+FP+FN);
    
    xlearn(i) = numLearners;
    yACC(i) = ACC;
    yTPR(i) = TPR;
    xFPR(i) = FPR;
    
    fprintf('TPR = %0.3f%%\n', TPR*100);
    fprintf('FPR = %0.3f%%\n', FPR*100);
    fprintf('ACC = %0.3f%%\n', ACC*100);
    
    numLearners = numLearners + 50;
    i = i + 1;
end

%% plot graphs

% Line graph (Mean)
subplot(1, 2, 1);
plot(xlearn, yACC);
xlabel('# of learners');
ylabel('Accuracy Rate');
title('# Learners vs. ACC (Mean)');

% Scatter plot (Mean)
subplot(1,2,2);
scatter(xFPR, yTPR);
xlabel('FALSE positive');
ylabel('TRUE positive');
title('TPR vs. FPR (Mean)');


% % HSV Histogram
% subplot(3,3,3);
% hist(HSV, 0:.1:255);
% title('HSV');
% 
% % YIQ Histogram
% subplot(3,3,4);
% hist(YIQ);
% title('YIQ');
% 
% % Lab Histogram
% subplot(3,3,5);
% hist(LAB);
% title('LAB');
% 
% % YCBCR Histogram
% subplot(3,3,6);
% hist(YCBCR);
% title('YCBCR');



%% TEST
% disp(HSV)

% plot hsv
% subplot(3,3,7);
% scatter(HSV(1), HSV(2));
% title('HSV');

% L = kmeans(LAB, 5, 'Replicates', 5);
% disp(L)
