close all
clear all

PCALHClassification = 1;
PCACDFClassification = 0;
LHNearestNeighborClassification = 0;
CDFNearestNeighborClassification = 0;


%Load images for testing and training
% load('C:\MelodysFiles\Histology\4X\4X\TERRY4-4X-HE-VASC.mat')
% load('C:\MelodysFiles\Histology\4X\4X\TERRY6-4X-HE.mat')

%Training Image
load('TERRY4-4X-HE-VASC.mat')
TrMask=imageData.roiMask;
TrImage1=imageData.originalImage;

%Testing Image
load('TERRY4-4X-HE-VASC.mat')
TeMask=imageData.roiMask;
TeImage1=imageData.originalImage;

%Specify what tissues
TissueTypes=[3 16 12];
%Specify number of lh's or cdf's to compute in PCA
K=2^6;
%Specify number of color bins
Bins=2^5;
%Support of window used to compute average lh's (or cdf's) and covariance matrix
Support = 2^5;
%Specify the scale of the Gaussian window
scale = 2^5;
%Specifiy number of eigenvalues to compute
NumberOfEigs = 5;


NumberOfTissues=length(TissueTypes);
%load training/testing images and masks
TrMask=TrMask(:,160:1359);
%Use this Mask to gaurantee non circular convolution in PCA.  Gaurantees
%pixels randomly chosen are not on the edge of image.
TrMask2=20*ones(size(TrMask));
TrMask2(Support+1:end-Support,Support+1:end-Support)=TrMask(Support+1:end-Support,Support+1:end-Support);
TrRGBImage=double(TrImage1(:,160:1359,:));

TeMask=TeMask(:,160:1359);
TeRGBImage=double(TeImage1(:,160:1359,:));

% imwrite(uint8(TrRGBImage),'NonBinnedTrainImage.png','png')
% imwrite(uint8(TeRGBImage),'NonBinnedTestImage.png','png')


clear('TrImage1','TeImage1')

[M,N]=size(TrRGBImage(:,:,1));

%Quantize images
DivideImageBy=(2^8)/Bins;

TrRGBImage=floor(TrRGBImage/DivideImageBy);
TeRGBImage=floor(TeRGBImage/DivideImageBy);

TrRImage = TrRGBImage(:,:,1);
TrBImage = TrRGBImage(:,:,3);

TeRImage = TeRGBImage(:,:,1);
TeBImage = TeRGBImage(:,:,3);

imwrite(uint8(DivideImageBy*TrRGBImage),'TrainImage.png','png')
imwrite(uint8(DivideImageBy*TeRGBImage),'TestImage.png','png')

clear('TrRGBImage','TeRGBImage')

%%%The Gaussian window
[Q,R]=ndgrid([0:M/2, -M/2+1:-1],[0:N/2,-N/2+1:-1]);
g=exp(-(Q.^2+R.^2)/scale^2);
g(M/2+1,:)=zeros(size(g(M/2+1,:)));
g(:,N/2+1)=zeros(size(g(:,N/2+1)));
clear('Q','R')

%Window for training
TrWindowSum=sum(sum(g(1:Support+1,1:Support+1)))+sum(sum(g(end-Support+1:end,1:Support+1)))+sum(sum(g(1:Support+1,end-Support+1:end)))+sum(sum(g(end-Support+1:end,end-Support+1:end)));
Trg=g/TrWindowSum;

%Window for testing
Support2=M/2;
g=fftshift(g);
g=padarray(g,[Support2 Support2]);
g=fftshift(g);
g=g/sum(sum(g));

FFTg=fft2(g);

%Terms for normalizing the star convolution
PaddedOnes=ones(M,N);
PaddedOnes=padarray(PaddedOnes,[Support2 Support2]);
ConvNormalizationTerm=ifft2(fft2(PaddedOnes).*FFTg);
ConvNormalizationTerm=ConvNormalizationTerm(Support2+1:end-Support2,Support2+1:end-Support2);
clear('PaddedOnes')

load('C:\Users\MMASSAR\Documents\Melodys Files\Histogram\MatlabPCA\PCA\PCAjournal1\PCAjournal.mat')

if PCALHClassification == 1
    V2=zeros(Bins^2,NumberOfEigs,NumberOfTissues);
    V=zeros(Bins,Bins,NumberOfEigs,NumberOfTissues);
    for a =1:NumberOfTissues
        [V2(:,:,a),D]=eigs(LHZZtranspose(:,:,a),NumberOfEigs);
        V(:,:,:,a)=reshape(V2(:,:,a),[Bins,Bins,NumberOfEigs]);
    end
    clear('V2')
    Temp=zeros(M,N);
    TissueProjections = zeros(M,N,NumberOfEigs,NumberOfTissues);
    TissueNormsSquared = zeros(M,N,NumberOfTissues);
     
    for RValue=0:Bins-1
        RValue
        for BValue=0:Bins-1
            TempR=(TeRImage==RValue);
            TempB=(TeBImage==BValue);
            Temp=TempR&TempB;
            Temp=padarray(Temp,[Support2 Support2]);
            if sum(sum(Temp))~=0
                gConvWithTemp=ifft2(FFTg.*fft2(Temp));
                gConvWithTemp=gConvWithTemp(Support2+1:end-Support2,Support2+1:end-Support2);
                gConvWithTemp=gConvWithTemp./ConvNormalizationTerm;
                for a=1:NumberOfTissues
                    TissueTemp = gConvWithTemp-LHxbar(RValue+1,BValue+1,a);
                    TissueNormsSquared(:,:,a)=TissueTemp.*TissueTemp+TissueNormsSquared(:,:,a);
                    for b=1:NumberOfEigs
                        TissueProjections(:,:,b,a)=TissueTemp*V(RValue+1,BValue+1,b,a)+TissueProjections(:,:,b,a);
                    end
                    if LHNearestNeighborClassification==1 
                        DistanceSquared(:,:,a) = DistanceSquared(:,:,a)+(gConvWithTemp-LHxbar(RValue+1,BValue+1,a)).^2;
                    end
                end
            end
        end
    end
end
    clear('LHZZtranspose','LHxbar','V','TissueTemp','TempR','TempB','Temp','gConvWithTemp')
   PCALHClassifiedImage=zeros(M,N,NumberOfEigs);
   PCALHClassificationConfusionMatrix = zeros(NumberOfTissues,NumberOfTissues,NumberOfEigs);
  for c=1:NumberOfEigs  
    TissueProjectionsNormSquared = zeros(M,N,NumberOfTissues);
    Data = zeros(M,N,NumberOfTissues);
    
    for a = 1:NumberOfTissues
        TissueProjectionsNormSquared(:,:,a) = sum((abs(TissueProjections(:,:,1:c,a))).^2,3);
        Data(:,:,a)=TissueNormsSquared(:,:,a)-TissueProjectionsNormSquared(:,:,a);
    end
    
    %%%Check each pixel to see what was best match and then assign label
    ClassifiedImage=zeros(M,N);
    [C,I]=min(Data,[],3);
    for d=1:NumberOfTissues
        ClassifiedImage(I==d)=TissueTypes(d);
    end
    
    %%%%%%%%%Confusion Matrix
    InMask=zeros(M,N,NumberOfTissues);
    ClassifiedAs=zeros(M,N,NumberOfTissues);
    ConfusionMatrix=zeros(NumberOfTissues,NumberOfTissues);
    for a=1:NumberOfTissues
        label=TissueTypes(a);
        InMask(:,:,a)=(TeMask==label);
        ClassifiedAs(:,:,a)=(ClassifiedImage==label);
    end
    
    for a=1:NumberOfTissues
        SumOfTissue=sum(sum(InMask(:,:,a)));
        for b=1:NumberOfTissues
            ConfusionMatrix(a,b)=sum(sum(InMask(:,:,a).*ClassifiedAs(:,:,b)))/SumOfTissue;
        end
    end
    %Relabel classified image and mask
%     ClassifiedImage((TeMask~=3)&(TeMask~=12)&(TeMask~=16))=5;
    ClassifiedImage(ClassifiedImage==3)=0;
    ClassifiedImage(ClassifiedImage==12)=1;
    ClassifiedImage(ClassifiedImage==16)=2;
    
    PCALHClassifiedImage(:,:,c)=ClassifiedImage;
    PCALHClassificationConfusionMatrix(:,:,c) = ConfusionMatrix;
    
end
toc    

imwrite(.95*(1/5)*PCALHClassifiedImage(:,:,1),'LHEig1ClassifiedImage.png','png')
imwrite(.95*(1/5)*PCALHClassifiedImage(:,:,2),'LHEig2ClassifiedImage.png','png')
imwrite(.95*(1/5)*PCALHClassifiedImage(:,:,3),'LHEig3ClassifiedImage.png','png')
imwrite(.95*(1/5)*PCALHClassifiedImage(:,:,4),'LHEig4ClassifiedImage.png','png')
imwrite(.95*(1/5)*PCALHClassifiedImage(:,:,5),'LHEig5ClassifiedImage.png','png')
clear('PCALHClassifiedImage')


% TeMask((TeMask~=3)&(TeMask~=12)&(TeMask~=16))=5;
TeMask(TeMask==3)=0;
TeMask(TeMask==12)=1;
TeMask(TeMask==16)=2;

imwrite(.95*(1/5)*TeMask,'TruthMask.png','png')