function [D_side, D_central, totalbytes, R0, R1] = MDLTPC_codec_2DES_new(imgname, cfgfile, qstep0, qstep1, errprob)
% Prediction compensated MD lapped transform with 2 descriptions:
%  [D_side,D_central, totalbytes,R0, R1] = MDLTPC_codec_2DES_new(imgname, cfgfile, qstep0, qstep1, errprob)
%   imgname:    image file that can be recognized by imread( ) function.
%   cfgfile:    pre/postfilter file,
%   qstep0:     I quantization step,
%   qstep1:     P quantization step,
%   errprob (optional): channel loss probability. If given, the expected PSNR will be calculated.
% Output:
%   D_side:     Average side PSNR.
%   D_central:  Central PSNR.
%   totalbtyes: Overall bytes of twodesriptions.
%   R0:         Average bit rate for base layer.
%   R1:         Average bit rate for enhanced layer. R0 + R1 should equal the target overall bit rate R.
% Other output files:
%   Encoded files:  
%       Description 1: base layer: MDLTPC_Zip_00.out, enhanced layer: MDLTPC_Zip_01.out
%       Description 2: base layer: MDLTPC_Zip_10.out, enhanced layer: MDLTPC_Zip_11.out
%       (Conceptually the combination of the base layer and enhanced layer is treated as one description).
%       (They are saved separately to simplify the codec and debugging)
%   Decoded images: MDLTPC_Dec_0.pgm, MDLTPC_Dec_1.pgm, MDLTPC_Dec_central.pgm
% Example: 
%   MDLTPC_codec_2DES_new('barb2.pgm','MDpredictN8_1_2.dat',19,96);
%
% Guoqian Sun, Jie Liang, June 10, 2007

if nargin < 4
	help MDLTPC_codec_2DES_new;
    return;
end

M = 8;

%read image
if ischar(imgname)
      img = imread(imgname); 
else
      img = imgname;  
end
if img == -1
    return;
end
img=double(img);

[row, col]=size(img);
rowblk = row / M; colblk = col / M;

[preflt, postflt, ErrorPostFilter, DCTwiener,Wiener1Nbr] = ReadPrePostFile(cfgfile);

%start timing
t0 = cputime;

%Encode, decode, return decoded intr and inter for eachd escription
[codec_intra0, codec_inter0, psnr0, Bytes0, Bytes1] = MDLTPC_codec_2DES_new_core(img, cfgfile, M, qstep0, qstep1, 0);
[codec_intra1, codec_inter1, psnr1, tmp0,   tmp1]   = MDLTPC_codec_2DES_new_core(img, cfgfile, M, qstep0, qstep1, 1);
Bytes0 = Bytes0 + tmp0;
Bytes1 = Bytes1 + tmp1;

%Side distortion: get avg MSE, convert to PSNR
sidemse = (255^2 * 10^(-psnr0 / 10) + 255^2 * 10^(-psnr1 / 10)) / 2;
D_side = 10*log10(255^2 / sidemse);

%Use intra blocks of two descriptions to compute central reconstruction
blkmap = kron([0, 1; 1 0], ones(M)); 
mask = kron(ones(rowblk/2, colblk/2), blkmap); 
idx0 = find(mask == 0);
idx1 = find(mask == 1);
centralimg = zeros(row, col);
centralimg(idx0) = reshape(codec_intra0, [row/2 * col, 1]);
centralimg(idx1) = reshape(codec_intra1, [row/2 * col, 1]);

centralimg = PrePostLT(centralimg, postflt, M);
centralimg = round(centralimg);
centralimg(centralimg > 255) = 255;
centralimg(centralimg < 0) = 0;
centralimg = uint8(centralimg);
imwrite(centralimg, 'MDLTPC_Dec_Central.pgm', 'pgm');
D_central = psnr(img, 'MDLTPC_Dec_Central.pgm');
centralmse = 255^2 * 10^(-D_central / 10);

t1 = cputime;
disp(sprintf('Time used: %f sec', t1 - t0));

totalbytes = Bytes0 + Bytes1;
R0 = Bytes0*8/col/row;
R1 = Bytes1*8/col/row;
disp(sprintf('1st Layer Bytes: %d, 2nd layer bytes: %d, total: %d', Bytes0, Bytes1, totalbytes));
disp(sprintf('Total rate=%f, R0= %f, R1 = %f, redun.=%f',totalbytes*8/col/row, R0, R1, R1/R0));

if nargin < 5
    disp(sprintf('Side PSNR = %6.3f dB, central PSNR = %6.3f dB',D_side,D_central));
else
    %calculate expected MSE and PSNR
    p0 = (1 - errprob)^2;
    p1 = 2 * errprob * (1 - errprob);    
    Expmse = (p1 * sidemse + p0 * centralmse) / (p0 + p1);
    ExpPSNR = 10 * log10(255^2 / Expmse);
    disp(sprintf('Side PSNR: %6.3f dB, central PSNR: %6.3f dB, Expected PSNR: %6.3f dB',D_side, D_central, ExpPSNR));
end

% data format for plotting figures
disp(sprintf('\t%8.3f, %8.3f, %8.3f, %8.3f, %8.3f, %8.3f; ...', D_side, D_central, qstep0, qstep1, R0, R1));