% read images


function [Y_encoded, PSNR, rate] = NoMC_3D_PCA_DCT(all_stepsizes,numberofPCAcoefs,PCAcoefsstepsize,filename,endframe,movietype,Block_Size_x,Block_Size_y,Block_Size_t,Support_XY_begin,Support_XY_end, Support_t_begin,stepsize,outputfile,StartSupport);

Y = ReadVideo(filename, endframe, movietype);

Support_t_end = 0;

sizex = size(Y{1,1},1);
sizey = size(Y{1,1},2);

Y_encoded = cell(1,endframe);
Y_encoded_PCA = cell(1,endframe);

PSNR = [];
rate=[];

% First find the Principal Components

All_Blocks =[]; 
for k=1:Block_Size_t:endframe-Block_Size_t+1    
    k
    for i=1:Block_Size_t
      IC(:,:,i) = Y{1,k+i-1};
    end
    kk=1+(k-1)/Block_Size_t; 
    for i=Support_XY_end/Block_Size_x+1:1:(sizex/Block_Size_x-Support_XY_end/Block_Size_x)
        for j=Support_XY_end/Block_Size_y+1:1:(sizey/Block_Size_y-Support_XY_end/Block_Size_y)
            All_Blocks =[ All_Blocks reshape(IC((i-1)*Block_Size_x+1:i*Block_Size_x,(j-1)*Block_Size_y+1:j*Block_Size_y,:),Block_Size_x*Block_Size_y*Block_Size_t,1) ];
        end
    end        
end

end_ = size(All_Blocks,2);
if (end_>5000)
    end_=5000;
end
All_Blocks = All_Blocks(:,1:end_);

average = mean(All_Blocks')';
Subtracted_All_Blocks = All_Blocks - average*ones(1,end_);

[ U,S,V ] =svd(Subtracted_All_Blocks);


% Get PCA prediction

for k=1:Block_Size_t:endframe-Block_Size_t+1
    k
    kk=1+(k-1)/Block_Size_t;
    
    % initialization of the current and previous matrices
    for i=1:Block_Size_t
      IC(:,:,i) = Y{1,k+i-1};
    end
    
    % DCT
    for i=1:Block_Size_t
        Y_encoded_PCA{1,k+i-1}=IC(:,:,i);
    end   
    
    for i=Support_XY_end/Block_Size_x+1:(sizex/Block_Size_x-Support_XY_end/Block_Size_x)
        for j=Support_XY_end/Block_Size_y+1:(sizey/Block_Size_y-Support_XY_end/Block_Size_y)
            PCACOEFs(:,i,j,kk) = U(:,1:numberofPCAcoefs)'*(reshape(IC((i-1)*Block_Size_x+1:i*Block_Size_x,(j-1)*Block_Size_y+1:j*Block_Size_y,:),Block_Size_x*Block_Size_y*Block_Size_t,1)-average);
            % Quantize PCA COEFS
            [PCACOEFsCompressed(:,i,j,kk)] = round(PCACOEFs(:,i,j,kk)/PCAcoefsstepsize);
            Inverse_PCACOEFs(:,i,j,kk) = (U(:,1:numberofPCAcoefs)*(PCACOEFsCompressed(:,i,j,kk)*PCAcoefsstepsize))+average;
            Inverse = reshape(Inverse_PCACOEFs(:,i,j,kk), Block_Size_x, Block_Size_y, Block_Size_t);
            for ii=1:Block_Size_t
                Y_encoded_PCA{1,k+ii-1}((i-1)*Block_Size_x+1:i*Block_Size_x,(j-1)*Block_Size_y+1:j*Block_Size_y)=Inverse(:,:,ii);
            end             
         end
    end        
end

if (0)
    clear PCACOEFs;
    clear PCACOEFsCompressed;
    clear Inverse_PCACOEFs;
end



PSNR =[];
rate=[];
for stepsize = all_stepsizes,

stepsize
DCTCOEFs = cell(sizex/Block_Size_x,sizey/Block_Size_y,endframe/Block_Size_t);
DCTCOEFsCompressed = cell(sizex/Block_Size_x,sizey/Block_Size_y,endframe/Block_Size_t);
Inverse_DCTCOEFs = cell(sizex/Block_Size_x,sizey/Block_Size_y,endframe/Block_Size_t);

for i=1:-Support_t_begin+Block_Size_t
    Y_encoded{1,i}=Y{1,i};
end

maxindexmat = -inf*ones(Block_Size_x, Block_Size_y, Block_Size_t);
minindexmat = inf*ones(Block_Size_x, Block_Size_y, Block_Size_t);

for k=-Support_t_begin+Block_Size_t:Block_Size_t:endframe-Block_Size_t
    
    % initialization of the current and previous matrices
    for i=1:Block_Size_t
      IC(:,:,i) = Y{1,k+i-1}-Y_encoded_PCA{1,k+i-1};
    end
    
    % DCT
    for i=1:Block_Size_t
        Y_encoded{1,k+i-1}=IC(:,:,i);
    end   
    
    for i=Support_XY_end/Block_Size_x+1:(sizex/Block_Size_x-Support_XY_end/Block_Size_x)
        for j=Support_XY_end/Block_Size_y+1:(sizey/Block_Size_y-Support_XY_end/Block_Size_y)
            DCTCOEFs{i,j,(k)/Block_Size_t+1} = DCT3(IC((i-1)*Block_Size_x+1:i*Block_Size_x,(j-1)*Block_Size_y+1:j*Block_Size_y,:));
            % Quantize
            [DCTCOEFsCompressed{i,j,(k)/Block_Size_t+1}, maxindexmat, minindexmat] = QuantizeCoefs3(DCTCOEFs{i,j,(k)/Block_Size_t+1}, stepsize, maxindexmat, minindexmat);
            Inverse_DCTCOEFs{i,j,(k)/Block_Size_t+1} = IDCT3(DCTCOEFsCompressed{i,j,(k)/Block_Size_t+1}*stepsize);
            for ii=1:Block_Size_t
                Y_encoded{1,k+ii-1}((i-1)*Block_Size_x+1:i*Block_Size_x,(j-1)*Block_Size_y+1:j*Block_Size_y)=Inverse_DCTCOEFs{i,j,(k)/Block_Size_t+1}(:,:,ii);
            end             
         end
    end  
    for i=1:Block_Size_t
        Y_encoded{1,k+i-1}=Y_encoded{1,k+i-1}+Y_encoded_PCA{1,k+i-1};
    end   

end

%PSNR calculation
PSNR = [ PSNR GetPSNR(Y, Y_encoded, [Support_XY_end+1 sizex-Support_XY_end], [1+Support_XY_end sizey-Support_XY_end], [StartSupport endframe-1])]


% entropy Calculation
EntropyPCA=0;
for kk = 1:numberofPCAcoefs
    All_X = PCACOEFsCompressed(kk,:,:,:);
    All_X = All_X(:);
    minx = min(All_X);
    maxx = max(All_X);
    XX=minx:1:maxx;
    [HH,C]=hist(All_X(:),XX);
    HH=HH/sum(HH);
    indexx=find(HH>0);
    EntropyPCA=EntropyPCA + sum(-HH(indexx).*log2(HH(indexx)));
end
    
EntropyPCA

% set up symbol frequency structures
symbolfreq = cell(Block_Size_x, Block_Size_y, Block_Size_t);
for t=1:Block_Size_t
	for i=1:Block_Size_x
	   for j=1:Block_Size_y
			symbolfreq{i,j,t} = zeros(1, maxindexmat(i,j,t)-minindexmat(i,j,t)+1);      
	   end
	end
end

% find symbol frequencies
for kk=StartSupport:Block_Size_t:endframe-Block_Size_t
    	kk    
    	for i=Support_XY_end/Block_Size_x+1:(sizex/Block_Size_x-Support_XY_end/Block_Size_x)
        for j=Support_XY_end/Block_Size_y+1:(sizey/Block_Size_y-Support_XY_end/Block_Size_y)
           symbolfreq = FindSymbolFreq3(DCTCOEFsCompressed{i,j,(kk)/Block_Size_t+1}, symbolfreq, maxindexmat, minindexmat);
        end
     end
end

% calculate probabilities
prob = cell(Block_Size_x, Block_Size_y, Block_Size_t);
H = zeros(Block_Size_x, Block_Size_y, Block_Size_t);
for t=1:Block_Size_t
	for i=1:Block_Size_x
	   for j=1:Block_Size_y
	      prob{i,j,t} = symbolfreq{i,j,t}/sum(symbolfreq{i,j,t});      
	     	nonzero = find(prob{i,j,t}~=0);
	      H(i,j,t) = sum(-prob{i,j,t}(nonzero).*log2(prob{i,j,t}(nonzero)));
	   end
	end
end
Htotal = sum(sum(sum(sum(H))));
rate = [ rate (Htotal+EntropyPCA)/(Block_Size_x*Block_Size_y*Block_Size_t)]
eval(['save ./SAVEALL/NoMC_3D_PCA_DCT_' num2str(stepsize) outputfile ' Htotal PSNR stepsize Y_encoded H']);
end





