

function [Y_encoded, PSNR, rate ] = Volume_MC_3D_DCT(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);
% Variables to hold Motion Vectors 
MotionVectorsX = zeros(sizex/Block_Size_x,sizey/Block_Size_y,endframe/Block_Size_t);
MotionVectorsY = zeros(sizex/Block_Size_x,sizey/Block_Size_y,endframe/Block_Size_t);
MotionVectorsT = zeros(sizex/Block_Size_x,sizey/Block_Size_y,endframe/Block_Size_t);

Y_encoded = cell(1,endframe);
PSNR = [];
rate=[];

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
    k
    kk=1+(k-(-Support_t_begin+Block_Size_t))/Block_Size_t;
    % initialization of the current and previous matrices
    for i=1:Block_Size_t
      IC(:,:,i) = Y{1,k+i};
    end
    for i=1:Block_Size_t-Support_t_begin
      IP(:,:,i) = Y_encoded{1,k+i-(-Support_t_begin+Block_Size_t)};
    end  
    % motion estimation
    [ MotionVectorsX(:,:,kk), MotionVectorsY(:,:,kk), MotionVectorsT(:,:,kk), IPredicted ]  = BlockMatching(IC,IP,Block_Size_x,Block_Size_y,Block_Size_t,Support_XY_begin,Support_XY_end,Support_t_begin);
    %[ MotionVectorsX(:,:,kk), MotionVectorsY(:,:,kk), MotionVectorsT(:,:,kk), IPredicted ]  = BlockMatching_STD(IC,IP,Block_Size_x,Block_Size_y,Block_Size_t,Support_XY_begin,Support_XY_end,Support_t_begin);
    % DCT
    Ierror = IC-IPredicted;
    for i=1:Block_Size_t
        Y_encoded{1,k+i}=Ierror(:,:,i);
    end   
    
    %for i=1:Block_Size_t
    %    Y_encoded{1,k+i}=IPredicted(:,:,i);
    %end   
    if(0)
    for i=1:Block_Size_t
        subplot(1,3,1);       
        imagesc(Ierror(:,:,i),[1 50]); 
        subplot(1,3,2);       
        imagesc(IPredicted(:,:,i),[1 255]); 
        subplot(1,3,3);
        imagesc(Y_encoded{1,k+i},[1 250]);
        error_sum = sum(sum(abs(Ierror(:,:,i))))
        keyboard;
    end  
    end
    %keyboard;
    
    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(Ierror((i-1)*Block_Size_x+1:i*Block_Size_x,(j-1)*Block_Size_y+1:j*Block_Size_y,1:Block_Size_t));
            % 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}((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}=Y_encoded{1,k+i}+IPredicted(:,:,i);
    end
end

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


% Entropy Calculation for motion vectors


XX=Support_XY_begin:1:Support_XY_end;
[HH,C]=hist(MotionVectorsX(:),XX);
HH=HH/sum(HH);
indexx=find(HH>0);
Entropydx=sum(-HH(indexx).*log2(HH(indexx)))

HH=hist(MotionVectorsY(:),XX);
HH=HH/sum(HH);
indexx=find(HH>0);
Entropydy=sum(-HH(indexx).*log2(HH(indexx)))


if (Support_t_begin<0)
    XX=Support_t_begin:1:0;
    [HH,C]=hist(MotionVectorsT(:),XX);
    HH=HH/sum(HH);
    indexx=find(HH>0);
    Entropydt=sum(-HH(indexx).*log2(HH(indexx)));
else
    Entropydt=0;
end


% calculate the entropy

% 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-1)/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)))); 
Hmotion = Entropydx + Entropydy + Entropydt
rate = (Htotal+Hmotion)/(Block_Size_x*Block_Size_y*Block_Size_t);


eval(['save ./SAVEALL/Volume_MC_3D_DCT_' num2str(stepsize) outputfile ' Htotal Hmotion PSNR rate stepsize Y_encoded H']);




