

%Read the mask file with class labels
maskfile = load('tm_train.mat');
maskim = maskfile.tm_train_mat;

%Read the 6 spectral bands
tm = load('tm.mat');
tm
b1 = tm.Converted_from_raw_data_1;
b2 = tm.Converted_from_raw_data_2;
b3 = tm.Converted_from_raw_data_3;
b4 = tm.Converted_from_raw_data_4;
b5 = tm.Converted_from_raw_data_5;
b6 = tm.Converted_from_raw_data_6;
nbands = 6;
[hgh,wid] = size(b1);







k1 = find(maskim==1);
[nofklass1,dummy] = size(k1);

k2 = find(maskim==2);
[nofklass2,dummy] = size(k2);

k3 = find(maskim==3);
[nofklass3,dummy] = size(k3);

k4 = find(maskim==4);
[nofklass4,dummy] = size(k4);


nofclasses = 4; 

klass1feat = zeros(nofklass1,nbands);
klass2feat = zeros(nofklass2,nbands);
klass3feat = zeros(nofklass3,nbands);
klass4feat = zeros(nofklass4,nbands);
currclass1 = 0;
currclass2 = 0;
currclass3 = 0;
currclass4 = 0;



%Put the training data into a data structure per class
for y=1:hgh
    for x=1:wid
            if maskim(y,x)==1
                currclass1 = currclass1 + 1;
                klass1feat(currclass1,1)=b1(y,x);
                klass1feat(currclass1,2)=b2(y,x);
                klass1feat(currclass1,3)=b3(y,x);
                klass1feat(currclass1,4)=b4(y,x);
                klass1feat(currclass1,5)=b5(y,x);
                klass1feat(currclass1,6)=b6(y,x);
            
            elseif maskim(y,x)==2
               currclass2 = currclass2 + 1;
                klass2feat(currclass2,1)=b1(y,x);
                klass2feat(currclass2,2)=b2(y,x);
                klass2feat(currclass2,3)=b3(y,x);
                klass2feat(currclass2,4)=b4(y,x);
                klass2feat(currclass2,5)=b5(y,x);
                klass2feat(currclass2,6)=b6(y,x);     
                
            elseif maskim(y,x)==3
               currclass3 = currclass3 + 1;
                klass3feat(currclass3,1)=b1(y,x);
                klass3feat(currclass3,2)=b2(y,x);
                klass3feat(currclass3,3)=b3(y,x);
                klass3feat(currclass3,4)=b4(y,x);
                klass3feat(currclass3,5)=b5(y,x);
                klass3feat(currclass3,6)=b6(y,x); 
            
                 
            elseif maskim(y,x)==4
                    currclass4 = currclass4 + 1;
                     klass4feat(currclass4,1)=b1(y,x);
                     klass4feat(currclass4,2)=b2(y,x);
                     klass4feat(currclass4,3)=b3(y,x);
                     klass4feat(currclass4,4)=b4(y,x);
                     klass4feat(currclass4,5)=b5(y,x);
                     klass4feat(currclass4,6)=b6(y,x);   
            end
    end
end


klassmean = zeros(nofclasses,nbands);
klasscov = zeros(nofclasses,nbands,nbands);


allklassfeat = [klass1feat; klass2feat; klass3feat; klass4feat];
% Compute mean and covariance matrixes

    for b=1:nbands
        klassmean(1,b)=mean(klass1feat(:,b));
        klassmean(2,b)=mean(klass2feat(:,b));
        klassmean(3,b)=mean(klass3feat(:,b));
        klassmean(4,b)=mean(klass4feat(:,b));
        
    end
    mu0 = mean(allklassfeat);
    
    mudiff1 = klassmean(1,:)-mu0;
    mudiff2 = klassmean(2,:)-mu0;
    mudiff3 = klassmean(3,:)-mu0;
    mudiff4 = klassmean(4,:)-mu0;
    
    nofpatt = nofklass1 + nofklass2 + nofklass3 + nofklass4;
    
    %Normalize the scatter by the number of samples in each class
    p1 = nofklass1/nofpatt;
    p2 = nofklass2/nofpatt;
    p3 = nofklass3/nofpatt;
    p4 = nofklass4/nofpatt;
    
    
    
    
    klasscov(1,:,:) = cov(klass1feat);
    klasscov(2,:,:) = cov(klass2feat);
    klasscov(3,:,:) = cov(klass3feat);
    klasscov(4,:,:) = cov(klass4feat);
    
    cov1 = cov(klass1feat);
    cov2 = cov(klass2feat);
    cov3 = cov(klass3feat);
    cov4 = cov(klass4feat);
    
    %Computation of the scatter matrices SB and SW is now very simple:
    SB = p1* mudiff1'*mudiff1 + p2* mudiff2'*mudiff2 + p3* mudiff3'*mudiff3 + p4* mudiff4'*mudiff4;
    
    SW = p1* cov1 + p2* cov2 + p3* cov3 + p4* cov4;
    
    %Verify that SB has rank nofclasses-1: 
    rank(SB)
    
    %Compute the critetion function and do eigenvector decomposition of it
    JMAT = inv(SW)*SB;
    [V,D] = eig(JMAT);
    
    %Form one feature at a time from the eigenvectors V
    FD1 = V(:,1);
    FD2 = [V(:,2)];
    FD3 = [V(:,3)];
    
    fd1im = zeros(size(b1));
    fd2im = zeros(size(b1));
    fd3im = zeros(size(b1));
    
    
    %Compute the 3 transformed Fisher features for our 4-class problem
    for y=1:hgh
        for x=1:wid
            currvec = [double(b1(y,x)); double(b2(y,x)); double(b3(y,x)); double(b4(y,x)); double(b5(y,x)); double(b6(y,x))];
            fd1im(y,x) = FD1'*currvec;
            fd2im(y,x) = FD2'*currvec;
            fd3im(y,x) = FD3'*currvec;
        end
    end
