%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% EE 684 - Spring 2016
%%%% Simple Code for CAD of mammography data sets 
%%%%
%%%% Your project should try to classify different ROIs (N/B/M) 
%%%% Target is to design a system with the best diagnostic performance
%%%%
%%%% We start from given ROIs - extracted from MiniMIAS database
%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

% Load all your learn/test data files
load ('EE684_CADe.mat');                % roi data for CADe problem
%load ('EE684_CADx.mat');                % roi data for CADx problem

% check Workspace on Matlab to see that the following new variables:
% EE684_CADe.mat:  roi_n (normal) , roi_a (Abnormal)  (size: 72x32x32)
% EE684_CADx.mat:  roi_n (normal) , roi_b (benign), roi_m (Malignant) (size: 50x32x32)
%

% Uncomment to use Cadx database
% roi_n= roi_b;
% roi_a= roi_m;


% Uncomment this to change ROI size - here is example for 8x8
% roi_n= roi_n(:,1:8,1:8);
% roi_a= roi_a(:,1:8,1:8);

% To simplify our project (without loss of generality), 
%  consider two class problem (Normal vs. Abnormal)
%
% Now, prepare the data files that represent your 2 classes

size1= size(roi_n);
Nim= size1(1);                      % number of images in data set
Roi_dim= size1(2);                  % size of ROI - Assume square ROI
Nlrn=   Nim/2;                      % set number of learning samples here
Ntst=  Nim-Nlrn;                    % rest is for testing
roi_n = double(roi_n);
roi_a = double(roi_a);

% Class 1 must be Normal - divide data between learning and testing sets
Nsig_learn= roi_n(1:Nlrn,:,:);       % select images are for learning
Nsig_test= roi_n(Nlrn+1:end,:,:);    % select images are for testing

% Class 2 must be abnormality
Asig_learn= roi_a(1:Nlrn,:,:);       % select images are for learning
Asig_test= roi_a(Nlrn+1:end,:,:);    % select images are for testing



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Now we start the feature extraction process
%%% Feature vector has dimension of Number of ROIs x Number of features
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% init feature vectors
Nsig_lrn_feat= zeros(Nlrn,1);
Nsig_tst_feat= zeros(Ntst,1);
Asig_lrn_feat= zeros(Nlrn,1);
Asig_tst_feat= zeros(Ntst,1);

% True class membership  - 1=Normal , 2=Abnormal
Nsig_lrn_class= zeros(Nlrn,1)+1;      
Nsig_tst_class= zeros(Ntst,1)+1;
Asig_lrn_class= zeros(Nlrn,1)+2;
Asig_tst_class= zeros(Ntst,1)+2;

%=====================================================
% Spatial Domain Features - Features directly from ROI 
%=====================================================

%---------------------------------------------------
%%% basic first order Statistical features
%---------------------------------------------------
% mean
Nsig_lrn_feat= mean(reshape(Nsig_learn,Nlrn,Roi_dim*Roi_dim),2);
Nsig_tst_feat= mean(reshape(Nsig_test,Ntst,Roi_dim*Roi_dim),2);
Asig_lrn_feat= mean(reshape(Asig_learn,Nlrn,Roi_dim*Roi_dim),2);
Asig_tst_feat= mean(reshape(Asig_test,Ntst,Roi_dim*Roi_dim),2);

% std dev
Nsig_lrn_feat= [Nsig_lrn_feat , std(reshape(Nsig_learn,Nlrn,Roi_dim*Roi_dim)')'];
Nsig_tst_feat= [Nsig_tst_feat , std(reshape(Nsig_test,Ntst,Roi_dim*Roi_dim)')'];
Asig_lrn_feat= [Asig_lrn_feat , std(reshape(Asig_learn,Nlrn,Roi_dim*Roi_dim)')'];
Asig_tst_feat= [Asig_tst_feat , std(reshape(Asig_test,Ntst,Roi_dim*Roi_dim)')'];

% Percentiles (Quantiles): 10, 30, 50, 70, 90
Nsig_lrn_feat= [Nsig_lrn_feat ,prctile(reshape(Nsig_learn,Nlrn,Roi_dim*Roi_dim),[10 30 50 70 90],2)];
Nsig_tst_feat= [Nsig_tst_feat ,prctile(reshape(Nsig_test,Ntst,Roi_dim*Roi_dim),[10 30 50 70 90],2)];
Asig_lrn_feat= [Asig_lrn_feat ,prctile(reshape(Asig_learn,Nlrn,Roi_dim*Roi_dim),[10 30 50 70 90],2)];
Asig_tst_feat= [Asig_tst_feat ,prctile(reshape(Asig_test,Ntst,Roi_dim*Roi_dim),[10 30 50 70 90],2)];


% compute Histogram
hist_levels= 32;
hist_n= zeros(Nim, hist_levels); 
hist_a= zeros(Nim, hist_levels); 

for k1= 1:Nim
    hist_n(k1,:) = hist(reshape(roi_n(k1,:,:),Roi_dim*Roi_dim,1),hist_levels);
    hist_a(k1,:) = hist(reshape(roi_a(k1,:,:),Roi_dim*Roi_dim,1),hist_levels);
    
    % normalize
    hist_n(k1,:)= hist_n(k1,:)/sum(hist_n(k1,:),2);
    hist_a(k1,:)= hist_a(k1,:)/sum(hist_a(k1,:),2);
end

Nhist_learn= double(hist_n(1:Nlrn,:,:));       % select GLCMs are for learning
Nhist_test= double(hist_n(Nlrn+1:end,:,:));    % select GLCMs are for testing
Ahist_learn= double(hist_a(1:Nlrn,:,:));       % select GLCMs are for learning
Ahist_test= double(hist_a(Nlrn+1:end,:,:));    % select GLCMs are for testing

% First Order Uniformity 
Nsig_lrn_feat= [Nsig_lrn_feat , sum(Nhist_learn.^2,2)];
Nsig_tst_feat= [Nsig_tst_feat , sum(Nhist_test.^2,2)];
Asig_lrn_feat= [Asig_lrn_feat , sum(Ahist_learn.^2,2)];
Asig_tst_feat= [Asig_tst_feat , sum(Ahist_test.^2,2)];

% First Order Entropy
Nsig_lrn_feat= [Nsig_lrn_feat , sum(Nhist_learn.*log(Nhist_learn+0.001),2)];
Nsig_tst_feat= [Nsig_tst_feat , sum(Nhist_test.*log(Nhist_test+0.001),2)];
Asig_lrn_feat= [Asig_lrn_feat , sum(Ahist_learn.*log(Ahist_learn+0.001),2)];
Asig_tst_feat= [Asig_tst_feat , sum(Ahist_test.*log(Ahist_test+0.001),2)];




%---------------------------------------------------
% Higher order Statistical features
%---------------------------------------------------
% compute GLCM 
GLCM_levels= 8;
GLCM_n= zeros(Nim, GLCM_levels, GLCM_levels); 
GLCM_a= zeros(Nim, GLCM_levels, GLCM_levels); 

for k1= 1:Nim
    GLCM_n(k1,:,:) = graycomatrix(squeeze(roi_n(k1,:,:)),'NumLevels', GLCM_levels, 'GrayLimits', [], 'Offset', [0 2]);
    GLCM_a(k1,:,:) = graycomatrix(squeeze(roi_a(k1,:,:)),'NumLevels', GLCM_levels, 'GrayLimits', [], 'Offset', [0 2]);
    
    % normalize
    GLCM_n(k1,:,:)= GLCM_n(k1,:,:)/sum(sum(GLCM_n(k1,:,:),3),2);
    GLCM_a(k1,:,:)= GLCM_a(k1,:,:)/sum(sum(GLCM_a(k1,:,:),3),2);
end

NGLCM_learn= double(GLCM_n(1:Nlrn,:,:));       % select GLCMs are for learning
NGLCM_test= double(GLCM_n(Nlrn+1:end,:,:));    % select GLCMs are for testing
AGLCM_learn= double(GLCM_a(1:Nlrn,:,:));       % select GLCMs are for learning
AGLCM_test= double(GLCM_a(Nlrn+1:end,:,:));    % select GLCMs are for testing

% compute GLCM features
% First Order Uniformity 
Nsig_lrn_feat= [Nsig_lrn_feat , sum(sum(NGLCM_learn.^2,3),2)];
Nsig_tst_feat= [Nsig_tst_feat , sum(sum(NGLCM_test.^2,3),2)];
Asig_lrn_feat= [Asig_lrn_feat , sum(sum(AGLCM_learn.^2,3),2)];
Asig_tst_feat= [Asig_tst_feat , sum(sum(AGLCM_test.^2,3),2)];

% First Order Entropy
Nsig_lrn_feat= [Nsig_lrn_feat , sum(sum(NGLCM_learn.*log(NGLCM_learn+0.001),3),2)];
Nsig_tst_feat= [Nsig_tst_feat , sum(sum(NGLCM_test.*log(NGLCM_test+0.001),3),2)];
Asig_lrn_feat= [Asig_lrn_feat , sum(sum(AGLCM_learn.*log(AGLCM_learn+0.001),3),2)];
Asig_tst_feat= [Asig_tst_feat , sum(sum(AGLCM_test.*log(AGLCM_test+0.001),3),2)];



%=====================================================
% Transform Domain Features - Transform then extract 
%=====================================================

% compute Fourier and wavelet transforms

% compute any of the spatial domain features in the transform domain



%##################################################################
% Basic Feature Selection
%##################################################################

%%%% combine all training samples/labels each in one array 
x_train= [Nsig_lrn_feat ; Asig_lrn_feat];
y_train= [Nsig_lrn_class ; Asig_lrn_class];

%%%% combine all training samples/labels each in one array 
x_test= [Nsig_tst_feat ; Asig_tst_feat];
y_test= [Nsig_tst_class ; Asig_tst_class];

%%%%%% t-Test
size1= size(Nsig_lrn_feat);
Nfeat= size1(2);
Pval= zeros(Nfeat,1);
for k1=1:Nfeat
    [H,P] = ttest2(Nsig_lrn_feat(:,k1),Asig_lrn_feat(:,k1));
    Pval(k1)= P;
end
        
figure(2)
plot(Pval);
textout= sprintf('T-Test: Number of useful features (P-Value < 0.05) = %d out of %d',sum(Pval<0.05), Nfeat);
disp(textout);

%%%%%% Correlation coefficient
% Corval= zeros(Nfeat,1);
% for k1=1:Nfeat
%     R= corrcoef(x_train(:,k1),y_train);
%     Corval(k1)= R(2,1);
% end
%         
% figure(3)
% plot(Corval);
% textout= sprintf('Correlation: Number of useful features (Correlation Coef > 0.5) = %d out of %d',sum(abs(Corval)>0.5), Nfeat);
% disp(textout);

%##################################################################
% Classification using several classical methods
%##################################################################

 [c1,err,P1] = classify(x_test,x_train,y_train,'linear');     % LDA - Linear
 c2 = classify(x_test,x_train,y_train,'quadratic');  % QDA - Quadratic
% c3 = knnclassify(x_test,x_train,y_train',1);        % KNN K=1
% c4 = knnclassify(x_test,x_train,y_train',3);        % KNN K=3
% c5 = knnclassify(x_test,x_train,y_train',5);        % KNN K=3

% -------------------------- New KNN ----------------------------------
mdl = fitcknn(x_train,y_train,'NumNeighbors',1);
c3= predict(mdl,x_test);
mdl.NumNeighbors = 3;
c4= predict(mdl,x_test);
mdl.NumNeighbors = 5;
c5= predict(mdl,x_test);
%----------------------------------------------------------------------

svmStruct = svmtrain(x_train,y_train);              % SVM - linear (default)
c6 = svmclassify(svmStruct,x_test);                 % SVM
svmStruct = svmtrain(x_train,y_train,'kernel_function','quadratic'); % SVM - quadratic
c7 = svmclassify(svmStruct,x_test);                 % SVM
svmStruct = svmtrain(x_train,y_train,'kernel_function','rbf'); % SVM - RBF
c8 = svmclassify(svmStruct,x_test);                 % SVM
svmStruct = svmtrain(x_train,y_train,'kernel_function','mlp'); % SVM - mlp
c9 = svmclassify(svmStruct,x_test);                 % SVM
%c10 = classify(x_test,x_train,y_train,'diaglinear');     % Naive Bayes

err1 = sum(y_test~= c1)/Ntst;         % mis-classification rate
err2 = sum(y_test~= c2)/Ntst;         % mis-classification rate
err3 = sum(y_test~= c3)/Ntst;         % mis-classification rate
err4 = sum(y_test~= c4)/Ntst;         % mis-classification rate
err5 = sum(y_test~= c5)/Ntst;         % mis-classification rate
err6 = sum(y_test~= c6)/Ntst;         % mis-classification rate
err7 = sum(y_test~= c7)/Ntst;         % mis-classification rate
err8 = sum(y_test~= c8)/Ntst;         % mis-classification rate
err9 = sum(y_test~= c9)/Ntst;         % mis-classification rate
err10 = sum(y_test~= c9)/Ntst;         % mis-classification rate

Error_Rate= [err1 , err2 , err3 , err4 , err5 , err6 , err7 , err8 , err9, err10]

%##################################################################
% Performance Evauation Example
%##################################################################

%%%%%% Sample confusion matrix computation 
%%% check that the total of the off-diagonal elements = error rate
conMat1 = confusionmat(y_test,c1)/Ntst      % confusion matrix for c1
TP= conMat1(1,1);
TN= conMat1(2,2);
FP= conMat1(1,2);
FN= conMat1(2,1);

Sens= TP/(TP+FN)
Spec= TN/(TN+FP)
PPV= TP/(TP+FP)
NPV= TN/(TN+FN)
Acc= (TP+TN)/(TP+TN+FP+FN)
FNrate= (1- Sens)
FPrate= (1- Spec)
PLR= Sens/(1-Spec)
NLR= (1-Sens)/Spec
DPrev= (TP+FN)/(TP+FN+TN+FP) 


%figure(1);
% [X,Y,T,AUC] = perfcurve(y_test,P1(:,2),2);
% plot(X,Y)
% xlabel('False positive rate'); ylabel('True positive rate')
% title('ROC for Classification')
% AUC



