function [codec_intra, codec_inter, decpsnr, Bytes0, Bytes1] = MDLTPC_codec_2DES_new_core(img, cfgfile, M, qstep0, qstep1, MDindex)
% Test performance of image coding with Wiener filter residual
% [outimg, outimg] = MDLTPC_codec_2DES_new_core(imgname, cfgfile, M, qstep0, qstep1,MDindex)
%   imgname:    image file that can be recognized by imread( ) function.
%   cfgfile:    pre/postfilter file,
%   M:          block size,
%   qstep0:     I quantization step,
%   qstep1:     P quantization step,
%   MDindex:    Index of the description
% example:
%  MDLTPC_codec_2DES_new_core('..\images\lena.pgm','MDpredict.dat',8,10,10,0);
%
% Jie Liang, June 10, 2007

if nargin < 6
    help MDLTPC_codec_2DES_new_core;
    return;
end

%output file name format: %MDLTPC_Zip_00, _01, _10, _11. 
%00 and 10: base layer. 01 and 11: enhanced layer
zfile0 =   ['MDLTPC_Zip_', num2str(MDindex), '0.out'];   
zfile1 =   ['MDLTPC_Zip_', num2str(MDindex), '1.out'];   
pgmfile =  ['MDLTPC_Dec_', num2str(MDindex), '.pgm'];   

DCT = dct(eye(M));

%read filters
[row, col]=size(img);
rowblk = row / M; colblk = col / M;
[preflt, postflt, ErrorPostFilter, DCTwiener,Wiener1Nbr] = ReadPrePostFile(cfgfile);
%Wiener1Nbr is for left/top nbr only, Wiener1Nbr2 is for right/bot nbr.
%Guoqian's code has problem here. Check.
Wiener1Nbr2 = fliplr(flipud(Wiener1Nbr));  

blkmap = kron([0, 1; 1 0], ones(M)); 
mask = kron(ones(rowblk/2, colblk/2), blkmap); 
intraidx = find(mask == MDindex);
interidx = find(mask == 1 - MDindex);

% Get reconstrucuted intra blocks
preimg = PrePostLT(img, preflt, M);     

%orgintra = preimg(intraidx);
%orgintra = reshape(orgintra, [row/2, col]);
%figure; imshow(orgintra, []);

outimg = imgdct(preimg, M);
outimg = dzoneq(outimg,qstep0);
refimg = outimg*qstep0;
refimg = imgdct(refimg,M,1);

%Get prediction residual for each inter block
for r = 1 : row/M
    for c = 1 : col/M
        curidx = mod(r+c, 2);   %MD index of the current block
        if curidx == MDindex
            %skip intra block 
            continue;
        end
        
        %Average hz and vt pred, get residual
        ridx = (r-1)*M + 1 : r*M;
        cidx = (c-1)*M + 1 : c*M;

        nbrcnt = 0;
        %horizontal prediction
        if (c == 1)    %left boundary
            nbrcnt = nbrcnt + 1;
            right = refimg(ridx, cidx + M);
            hpred = (Wiener1Nbr2 * right')';
        elseif (c == col/M)   %right boundary
            nbrcnt = nbrcnt + 1;
            left = refimg(ridx, cidx - M);
            hpred = (Wiener1Nbr * left')';
        else %other middle blocks
            nbrcnt = nbrcnt + 2;
            hnbr = [refimg(ridx,cidx - M), refimg(ridx,cidx + M)];
            hpred = (DCTwiener * hnbr')' * 2;
        end

        %vertical prediction
        if (r == 1)    %top boundary
            nbrcnt = nbrcnt + 1;
            bot = refimg(ridx + M, cidx);
            vpred = (Wiener1Nbr2 * bot);
        elseif (r == row/M)   %bot boundary
            nbrcnt = nbrcnt + 1;
            top = refimg(ridx - M, cidx);
            vpred = (Wiener1Nbr * top);
        else %other middle blocks
            nbrcnt = nbrcnt + 2;
            vnbr = [refimg(ridx - M, cidx); refimg(ridx + M,cidx)];
            vpred = (DCTwiener * vnbr) * 2;
        end

        %avg hz and vt prediction, more weight if 2 nbrs are available
        pred = (hpred + vpred) / nbrcnt;
        
       outimg(ridx, cidx) = preimg(ridx, cidx) - pred;       
    end
end

%save intra and inter into FAT subimage (better entropy coding than TALL image)
intra = outimg(intraidx);
intra = reshape(intra, [row/2, col]);   

fid = fopen('tmp0.raw', 'wb');
fwrite(fid, col, 'int');    %w
fwrite(fid, row/2, 'int');    %h
fwrite(fid, intra', 'float');
fclose(fid);

%DCT and Q for inter blocks, save for entropy coding
inter = outimg(interidx);
inter = reshape(inter, [row/2, col]);

%figure; imshow(inter, []);
if 0
	fid = fopen('mdc_lena_inter.raw','wb');
    fwrite(fid, row/2, 'int');
    fwrite(fid, col, 'int');
    fwrite(fid, inter, 'float');
    fclose(fid);
end
orginter = inter;

inter = imgdct(inter, M);      
inter = dzoneq(inter, qstep1); 

fid = fopen('tmp1.raw', 'wb');
fwrite(fid, col, 'int');    %w
fwrite(fid, row/2, 'int');    %h
fwrite(fid, inter', 'float');   
fclose(fid);

%-c 1: entropy only, -f 1: ECEB
[e,out] = dos(['iencoder.exe tmp0.raw ', zfile0, ' ', cfgfile,' 1 -c 1 -f 1']);
point = strfind(out, 'bytes');
Bytes0 = str2double(out(point(1)+5:point(2)-1));

[e,out] = dos(['iencoder.exe tmp1.raw ', zfile1, ' ', cfgfile,' 1 -c 1 -f 1']);
point = strfind(out, 'bytes');
Bytes1 = str2double(out(point(1)+5:point(2)-1));

%decoding
[e,out] = dos(['idecoder.exe ', zfile0, ' dtmp0.raw ', cfgfile, ' -c 1 -f 1']);
[e,out] = dos(['idecoder.exe ', zfile1, ' dtmp1.raw ', cfgfile, ' -c 1 -f 1']);

fid = fopen('dtmp0.raw', 'rb');
codec_col = fread(fid, 1, 'int');
codec_row = fread(fid, 1, 'int');
codec_intra = fread(fid, [codec_col, codec_row], 'float32');
codec_intra = codec_intra';
fclose(fid);

fid = fopen('dtmp1.raw', 'rb');
codec_col = fread(fid, 1, 'int');
codec_row = fread(fid, 1, 'int');
codec_inter = fread(fid, [codec_col, codec_row], 'float32');
codec_inter = codec_inter';
fclose(fid);

%DeQ, IDCT, fill up into recimg
codec_intra = codec_intra * qstep0;
codec_intra = imgdct(codec_intra, M, 1);    %IDCT
codec_inter = codec_inter * qstep1;
codec_inter = imgdct(codec_inter, M, 1);    %IDCT

recimg = zeros(row, col);
recimg(intraidx) = reshape(codec_intra, [row/2 * col, 1]);
recimg(interidx) = reshape(codec_inter, [row/2 * col, 1]);

%Reconstruct inter blocks
refimg = [];    %clear old refimg, to ensure using decoded result
for r = 1 : row/M
    for c = 1 : col/M
        curidx = mod(r+c, 2);   %MD index of the current block
        if curidx == MDindex
            %skip intra block 
            continue;
        end
        
        %Average hz and vt pred, get residual
        ridx = (r-1)*M + 1 : r*M;
        cidx = (c-1)*M + 1 : c*M;

        nbrcnt = 0;
        %horizontal prediction
        if (c == 1)    %left boundary
            nbrcnt = nbrcnt + 1;
            right = recimg(ridx, cidx + M);
            hpred = (Wiener1Nbr2 * right')';
        elseif (c == col/M)   %right boundary
            nbrcnt = nbrcnt + 1;
            left = recimg(ridx, cidx - M);
            hpred = (Wiener1Nbr * left')';
        else %other middle blocks
            nbrcnt = nbrcnt + 2;
            hnbr = [recimg(ridx,cidx - M), recimg(ridx,cidx + M)];
            hpred = (DCTwiener * hnbr')' * 2;
        end

        %vertical prediction
        if (r == 1)    %top boundary
            nbrcnt = nbrcnt + 1;
            bot = recimg(ridx + M, cidx);
            vpred = (Wiener1Nbr2 * bot);
        elseif (r == row/M)   %bot boundary
            nbrcnt = nbrcnt + 1;
            top = recimg(ridx - M, cidx);
            vpred = (Wiener1Nbr * top);
        else %other middle blocks
            nbrcnt = nbrcnt + 2;
            vnbr = [recimg(ridx - M, cidx); recimg(ridx + M,cidx)];
            vpred = (DCTwiener * vnbr) * 2;
        end

        %avg hz and vt prediction, more weight if 2 nbrs are available
        pred = (hpred + vpred) / nbrcnt;
        
       recimg(ridx, cidx) = recimg(ridx, cidx) + pred;       
    end
end

recimg = PrePostLT(recimg, postflt, M);
recimg = round(recimg);
recimg(recimg > 255) = 255;
recimg(recimg < 0) = 0;
recimg = uint8(recimg);
imwrite(recimg, pgmfile, 'pgm');
decpsnr = psnr(img, pgmfile);

dos('del tmp0.raw');
dos('del tmp1.raw');
dos('del dtmp0.raw');
dos('del dtmp1.raw');
