clear

% actual bleb data created from manuallyFindBlebs.m 
filenameTrain = 'BLEBS_agar_2per_4.mat'
filenameTest = 'BLEBS_agar_2per_5.mat'

load( filenameTrain )  
% contains the following variables
%       actualblebdata (this contains the following:
%          bleb:           4D array
%                          bleb(:,:,j,i) is the jth bleb at time frame i
%          blebsize:       blebsize(j,i) is the size of bleb j at time fram i
%          alpha:          alpha(j,i) is the "roundness" of bleb j at time frame i
%          showAllBlebs:   showAllBlebs(:,:,i) is the union of all bleb(:,:,j,i)
%          numberofblebs:  numberofblebs(i) is the number of blebs at time frame i
%          numberofimages: total number of images
%          images:         structure containing the image matices) 
%       notblebdata 
%          blebsize, alpha:    (j,i) is the size and roundness of bleb j at time frame i                  

totalNumberOfBlebs = sum( actualblebdata.numberofblebs )
totalNumberOfNotBlebs = sum( notblebdata.alpha(notblebdata.alpha ~= 0) )

[BLEB NOTBLEB] = calculations( actualblebdata , notblebdata );
   
% display to console
BLEB
NOTBLEB

% % plot Gaussian contours for no blebs
X = [BLEB.sizes BLEB.alpha] ;
Y = [NOTBLEB.sizes NOTBLEB.alpha] ;
% x1 = -10:.1:500; x2 = 0:.01:1;
% [X1,X2] = meshgrid(x1,x2);
% F = mvnpdf([X1(:) X2(:)],mean(Y),cov(Y));
% F = reshape(F,length(x2),length(x1));
% contour(x1,x2,F,'LineColor','r');

% % same for blebs
% x1 = -200:.5:5000; x2 = 0:.01:1;
% [X1,X2] = meshgrid(x1,x2);
% G = mvnpdf([X1(:) X2(:)],mean(X),cov(X));
% G = reshape(G,length(x2),length(x1));
% contour(x1,x2,G,'LineColor','b')

% create Gaussian Bayes classifier using the training data
BAYES_CLASSIFIER.numberofblebs = size( BLEB.sizes );
BAYES_CLASSIFIER.numberofnotblebs = size( NOTBLEB.sizes );
BAYES_CLASSIFIER.meanblebs = mean(X);
BAYES_CLASSIFIER.covblebs = cov(X);
BAYES_CLASSIFIER.meannotblebs = mean(Y);
BAYES_CLASSIFIER.covnotblebs = cov(Y);
BAYES_CLASSIFIER.type = 0;

% training error for the Gaussian Bayes classifier
BAYESTRAININGERROR = calculateError( BLEB , NOTBLEB , BAYES_CLASSIFIER )

% define the Simple Model
SIMPLEMODEL.uppersizethreshold = max( NOTBLEB.sizes ) + 1;
SIMPLEMODEL.lowersizethreshold = min( BLEB.sizes ) - 1;
% calculate the alpha threshold
n = 1; count = 0;
for n = 1 : size( BLEB.sizes )
   if BLEB.sizes(n) < SIMPLEMODEL.uppersizethreshold && BLEB.sizes(n) >= SIMPLEMODEL.lowersizethreshold
    count = count + 1;
    A( count ) = BLEB.alpha( n );
   end
end
n = 1; count = 0;
for n = 1 : size( NOTBLEB.sizes )
   if NOTBLEB.sizes(n) < SIMPLEMODEL.uppersizethreshold && NOTBLEB.sizes(n) >= SIMPLEMODEL.lowersizethreshold
    count = count + 1;
    B( count ) = NOTBLEB.alpha( n );
   end
end
SIMPLEMODEL.alphathreshold = (1/2) * (mean(A) + mean(B));
SIMPLEMODEL.type = 1;

% plot blebs in blue in blebsize x alpha space
b1 = plot( BLEB.sizes , BLEB.alpha , 'b.', 'MarkerSize', 30) ;
hold on 
% plot not blebs in red
n1 = plot(NOTBLEB.sizes , NOTBLEB.alpha, 'r.', 'MarkerSize', 30);
% plot average points
b2 = plot( BLEB.avSize, BLEB.avAlpha, 'b.','MarkerSize',60);
n2 = plot( NOTBLEB.avSize , NOTBLEB.avAlpha , 'r.', 'MarkerSize',60);
title('Bleb Size verses Alpha')
xlabel('Bleb Size (/pixels)')
ylabel('Alpha')

b3 = plot( (max( NOTBLEB.sizes ) + min( BLEB.sizes ))/2, mean(A), 'bo','MarkerSize',30);
n3 = plot( (max( NOTBLEB.sizes ) + min( BLEB.sizes ))/2, mean(B), 'ro', 'MarkerSize',30);

% draw classifying regions on the plot
legend([b1 n1],{'Blebs','Not blebs'})
 rectangle('Position',[min( BLEB.sizes )-1 (1/2)*(mean(A) + mean(B)) (max(NOTBLEB.sizes)-min( BLEB.sizes )+1) (0.35-(1/2)*(mean(A) + mean(B)))], 'FaceColor',[0 0 1 0.1],'EdgeColor',[0 0 1 0.1])
 rectangle('Position',[min( BLEB.sizes )-1 0 (max( NOTBLEB.sizes )-min( BLEB.sizes )+1) (1/2)*(mean(A) + mean(B))], 'FaceColor',[1 0 0 0.1],'EdgeColor',[1 0 0 0.1])
 rectangle('Position',[max( NOTBLEB.sizes )+1 0 1800 0.35], 'FaceColor',[0 0 1 0.1],'EdgeColor',[0 0 1 0.1])
 rectangle('Position',[0 0 min( BLEB.sizes )-1 0.35], 'FaceColor',[1 0 0 0.1],'EdgeColor',[0 0 1 0.1])
 alpha(.5)

% display to console
SIMPLEMODEL

TRAININGERROR = calculateError( BLEB, NOTBLEB , SIMPLEMODEL )

% calculate the generalisation error
clearvars -except BAYES_CLASSIFIER SIMPLEMODEL filenameTest
load( filenameTest );

totalNumberOfBlebs = sum( actualblebdata.numberofblebs );

[BLEB NOTBLEB] = calculations(actualblebdata, notblebdata);

GENERALISATIONERROR = calculateError( BLEB, NOTBLEB, SIMPLEMODEL )
BAYESGENERALISATIONERROR = calculateError( BLEB , NOTBLEB, BAYES_CLASSIFIER )

function ERROR = calculateError( BLEB , NOTBLEB , MODEL)  

n = 1;
count1 = 0; % correctly identified blebs
for n = 1 : size( BLEB.sizes )
   
    if MODEL.type == 0
        if probBleb( BLEB.sizes(n), BLEB.alpha(n), MODEL) > 0.5
            count1 = count1 + 1;
        end
    elseif MODEL.type == 1
        if  BLEB.sizes(n) > MODEL.uppersizethreshold
            count1 = count1 + 1;
        elseif BLEB.sizes(n) > MODEL.lowersizethreshold && BLEB.alpha(n) > MODEL.alphathreshold
            count1 = count1 + 1;
        end
    end
end

n = 1;
count2 = 0; % correctly identified as not blebs
for n = 1 : size( NOTBLEB.sizes )
   
     if MODEL.type == 0
        if probBleb( NOTBLEB.sizes(n), NOTBLEB.alpha(n),MODEL) < 0.5
        count2 = count2 + 1;
    end
     elseif MODEL.type == 1
        if NOTBLEB.sizes(n) < MODEL.lowersizethreshold
            count2 = count2 + 1;
        elseif NOTBLEB.sizes(n) < MODEL.uppersizethreshold && NOTBLEB.alpha(n) < MODEL.alphathreshold
            count2 = count2 + 1;
        end
    end
    
end

ERROR.positivepositives = strcat( int2str(count1), {' out of '}, int2str(size(BLEB.sizes,1)));
ERROR.negativenegatives = strcat( int2str(count2), {' out of '}, int2str(size(NOTBLEB.sizes,1)));
ERROR.blebdetectpercentage = count1 / size(BLEB.sizes,1) ;
ERROR.notblebdetectpercentage = count2 / size(NOTBLEB.sizes,1) ;
ERROR.accuracy = (count1 + count2) / (size(BLEB.sizes,1) + size(NOTBLEB.sizes,1)) ;
ERROR.weightedaccuracy = ((count1 / size(BLEB.sizes,1)) + (count2 / size(NOTBLEB.sizes,1) )) / 2;

end

function ret =  probBleb ( xsize , alpha , model)

    p = model.numberofblebs / ( model.numberofblebs + model.numberofnotblebs) ;
    ret = p * mvnpdf([xsize alpha],model.meanblebs,model.covblebs) / ...
        (p * mvnpdf([xsize alpha],model.meanblebs,model.covblebs) + ...
        (1 - p) * mvnpdf([xsize alpha],model.meannotblebs,model.covnotblebs))  ;

end

function [BLEB NOTBLEB] = calculations( actualblebdata , notblebdata )

% converts the 2d blebsize matrix into a 1d array
BLEB.sizes = actualblebdata.blebsize(:);
BLEB.sizes = BLEB.sizes(BLEB.sizes ~= 0);
BLEB.alpha = actualblebdata.alpha(:);
BLEB.alpha = BLEB.alpha(BLEB.alpha ~= 0);

NOTBLEB.sizes = notblebdata.blebsize(:);
NOTBLEB.sizes = NOTBLEB.sizes(NOTBLEB.sizes ~= 0);
NOTBLEB.alpha = notblebdata.alpha(:);
NOTBLEB.alpha = NOTBLEB.alpha(NOTBLEB.alpha ~= 0);

% some calculations
BLEB.avSize = mean (BLEB.sizes);
BLEB.stdSize = std ( BLEB.sizes );
BLEB.minSize = min( BLEB.sizes );
BLEB.maxSize = max(  BLEB.sizes );

BLEB.avAlpha = mean (BLEB.alpha );
BLEB.stdAlpha = std ( BLEB.alpha  );
BLEB.minAlpha = min(BLEB.alpha );
BLEB.maxAlpha = max(BLEB.alpha );

NOTBLEB.avSize = mean (NOTBLEB.sizes);
NOTBLEB.stdSize = std ( NOTBLEB.sizes );
NOTBLEB.minSize = min( NOTBLEB.sizes );
NOTBLEB.maxSize = max(  NOTBLEB.sizes );

NOTBLEB.avAlpha = mean (NOTBLEB.alpha );
NOTBLEB.stdAlpha = std ( NOTBLEB.alpha  );
NOTBLEB.minAlpha = min(NOTBLEB.alpha );
NOTBLEB.maxAlpha = max(NOTBLEB.alpha );

end