function recon = dft_coder(input,alloc,factor,blocksize)

% function   recon = dft_coder(input,alloc,factor,blocksize)
%
%  This function applies  dft_block.m  with the given blocksize to input.
%  It then quantizes the resulting dft coefficients with quantizers determined
%  by alloc and factor.
%
%  alloc is a vector of dimension blocksize that determines how many bits are
%  to be used to quantize each coefficient. 
%
%  The range of the quantizer for the ith coefficient is [mins, maxs].
%  mins  is the larger of the minimum value of the coef's of type i, and
%  the mean minus factor times the standard deviation of such coefficients.
%
%  maxs  is the smaller of the maximum value of the coef's of type i, and
%  the mean plus factor times the standard deviation of such coefficients.
%
%  alloc, factor and blocksize are optional arguments.  Their default values
%  are  [8 8 8 ... 8],  5  and  8.
%  
%  The script displays a table of useful values, such as the mean and std deviation
%  of the coefficients of each type, and MSE and RMSE of the quantization errors
%  for the various coefficients.
%
%  MB, DN 10/4/01


% Define useful inline functions
rms = inline('sqrt(mean(mean(x.^2,2),1))');
msv = inline('(mean(mean(x.^2,2),1))');

if ~exist('alloc','var') | isempty('alloc')
    alloc = [8 8 8 8 8 8 8 8];
end

if ~exist('factor','var')| isempty('factor')
	factor = 5 ;  %  factor that multiplies the rms value
end

factor

if ~exist('blocksize','var') | isempty(blocksize)
    blocksize = 8;
end

%calculate MSV of input
inmsv = msv(input) ;
inrmsv = sqrt(inmsv);

% Calculate the DFT coefficients
output = dft_block(input,blocksize);
outputsize = size(output) ;

%calculate their mean values
outmeann(1:blocksize) = mean(mean(output)) ;

%calculate their ms and rms values
outmsv(1:blocksize) = msv(output); 
outrmsv = sqrt(outmsv);

%calculate var and std
outvar = outmsv - outmeann.^2 ; 
outstdd = sqrt(outvar) ;

%calculate their max and min values
maxx(1:blocksize) = max(max(output,[],2),[],1) ;
minn(1:blocksize) = min(min(output,[],2),[],1) ;


%%% Quantize the coefficients

% Default quantizer range is based on RMS values and actual maximum values
if ~exist('maxs','var') | isempty(maxs)
    maxs = min(outmeann + factor*outstdd,maxx);
	mins = max(outmeann - factor*outstdd,minn) ;
end

%maxs
maxs = shiftdim(maxs)';
%mins = [0 -maxs(2:blocksize)];

%calculate the delta's
delta = (maxs - mins) .* 2.^(-alloc) ;

%do the quantization
q_output = output;
for index=2:blocksize
    q_output(:,:,index) = quantize(output(:,:,index), 2^alloc(index),mins(index),maxs(index),1);
end
q_output(:,:,1) = quantize(output(:,:,1), 2^alloc(1),mins(1),maxs(1));

%compute vector of mse and rmse for coefficeints
for index=1:blocksize
	cmse(index) = msv(output(:,:,index)-q_output(:,:,index)) ;
end
crmse = sqrt(cmse) ;
 
%compute average mse for coefficients
meancmse = mean(cmse);  %/blocksize;
rmscmse = sqrt(meancmse); 

% Reconstruct the image
recon = abs(inverse_dft_block(q_output,blocksize));

% compute the MSE
imse = msv(input - recon) ;
irmse = sqrt(imse) ;

% Display some useful data

fprintf(1,'\n');
disp('           Re{c[0]} Re{c[1]} Re{c[2]} Re{c[3]} Re{c[4]} Im{c[1]} Im{c[2]} Im{c[3]}')
disp('-------------------------------------------------------------------------------')
fprintf(1,'Coefficient Data:\n\n');
fprintf(1,'Max:     ');
fprintf(1,'%8.1f ',maxx);
fprintf(1,'\nMin:     ');
fprintf(1,'%8.1f ',minn);
fprintf(1,'\n\nMean:     ');
fprintf(1,'%8.2f ',outmeann);
fprintf(1,'\nStd:     ');
fprintf(1,'%8.1f ',outstdd);
fprintf(1,'\n\nXMax:    ');
fprintf(1,'%8.1f ',maxs);
fprintf(1,'\nXMin:    ');
fprintf(1,'%8.1f ',mins);

fprintf(1,'\n\nBits:    ');
fprintf(1,'%8d ',alloc);
fprintf(1,'\n\nDelta:   ');
fprintf(1,'%8.2f ',delta);
fprintf(1,'\n\nRMS Err: ');
fprintf(1,'%8.2f ',crmse);
fprintf(1,'\nMSE:     ');
fprintf(1,'%8.2f ',cmse);
%fprintf(1,'\ndel^2/12:');
%fprintf(1,'%8.2f ',delta.^2/12);

fprintf(1,'\n\nSummary Data:\n');

%fprintf(1,'\n\nInput MSV:       %6.1f',inmsv) ;
fprintf(1,'\nInput RMS value:    %6.1f', inrmsv);

%fprintf(1,'\n\nOutput MSV:       %6.1f', mean(outmsv)); %/blocksize);
%fprintf(1,'\nOut MSV x 8:     %6.1f', 8 * mean(outmsv) ) ;

%fprintf(1,'\n\nAvg MSE(k)''s:     %8.4f \n',meancmse) ;
%fprintf(1,'Avg MSE(k)''s*8:   %8.4f \n',meancmse*8) ;
%fprintf(1,'RMSE:             %8.4f \n', sqrt(meancmse*8)) ;

msesq = (256/2^(sum(alloc)/blocksize))^2 / 12 ;
%fprintf(1,'\nTotal RMS Error:  %8.4f        RMSE for scalar quant:  %8.4f\n',irmse, sqrt(msesq));
%fprintf(1,'Total MSE:        %8.4f        MSE for scalar quant:   %8.4f\n\n',imse, msesq) ;

fprintf(1,'\nTotal RMS Error:  %8.2f\n',irmse);
fprintf(1,'Total MSE:        %8.2f\n\n',imse) ;

fprintf(1,'Total Bits Per Pixel: %d \n\n',sum(alloc)/blocksize);


