% QTLLST2DMDL find a QT decomposition based on MDL
%   [LIST,COST]=QTLLST2DMDL(COEF,LMAX,RATE) returns the optimal level list for
%   COEF going down to LMAX and using the quantization matrices at rate
%   RATE. This program uses the GLOBAL value of the entropy.
%
%   Remark: this program uses the global variable quantization.
%
%   See also: QTLLST2DMDLORIG,QTLLST2DMDL2,QTLLST2DL1.
function [liste, coef]=qtllst2dmdl(signal, lmax, rate)

[liste, coef] = recurs_testdst2d(signal, 0, lmax,rate);


function [liste, coef, border] = recurs_testdst2d(signal, level, lmax, rate)

[m,n] = size(signal);
m2 = (m+1)/2;
n2 = (n+1)/2;

b1 = signal(:,1);
b2 = signal(end,:);
b3 = signal(:,end);
b4 = signal(1,:);


%Quantize and unquantize the borders
b1(2:end-1) = b1(2:end-1) - (b1(1) + [1:m-2]'/(m-1) * (b1(end)-b1(1)));
b1(2:end-1) = sqrt(2/(m-1)) * dst(b1(2:end-1));
b1(2:end-1) = quantizeborder(b1(2:end-1),rate,1,0);
b1(2:end-1) = quantizeborder(b1(2:end-1),rate,1,1);
b1(2:end-1) = sqrt(2/(m-1)) * dst(b1(2:end-1));
b1(2:end-1) = b1(2:end-1) + (b1(1) + [1:m-2]'/(m-1) * (b1(end)-b1(1)));

b2(2:end-1) = b2(2:end-1) - (b2(1) + [1:n-2]/(n-1) * (b2(end)-b2(1)));
b2(2:end-1) = sqrt(2/(n-1)) * dst(b2(2:end-1));
b2(2:end-1) = quantizeborder(b2(2:end-1),rate,0,0);
b2(2:end-1) = quantizeborder(b2(2:end-1),rate,0,1);
b2(2:end-1) = sqrt(2/(n-1)) * dst(b2(2:end-1));
b2(2:end-1) = b2(2:end-1) + (b2(1) + [1:n-2]/(n-1) * (b2(end)-b2(1)));

b3(2:end-1) = b3(2:end-1) - (b3(1) + [1:m-2]'/(m-1) * (b3(end)-b3(1)));
b3(2:end-1) = sqrt(2/(n-1)) * dst(b3(2:end-1));
b3(2:end-1) = quantizeborder(b3(2:end-1),rate,1,0);
b3(2:end-1) = quantizeborder(b3(2:end-1),rate,1,1);
b3(2:end-1) = sqrt(2/(m-1)) * dst(b3(2:end-1));
b3(2:end-1) = b3(2:end-1) + (b3(1) + [1:m-2]'/(m-1) * (b3(end)-b3(1)));

b4(2:end-1) = b4(2:end-1) - (b4(1) + [1:n-2]/(n-1) * (b4(end)-b4(1)));
b4(2:end-1) = sqrt(2/(n-1)) * dst(b4(2:end-1));
b4(2:end-1) = quantizeborder(b4(2:end-1),rate,0,0);
b4(2:end-1) = quantizeborder(b4(2:end-1),rate,0,1);
b4(2:end-1) = sqrt(2/(n-1)) * dst(b4(2:end-1));
b4(2:end-1) = b4(2:end-1) + (b4(1) + [1:n-2]/(n-1) * (b4(end)-b4(1)));

lap = zeros(m,n);
lap(:,1) = b1;
lap(1,:) = b4;
lap(end,:) = b2;
lap(:,end) = b3;

lap = solvelaplace(lap);


r = signal(2:end-1,2:end-1) - lap;

r = dst(dst(r)')' * sqrt(4/(m-1)/(n-1));

%
r = quantizeinside(r,rate, 0);
%
r = r(zigind(m-2));
indmax = max(find(r~=0));
r = r(1:indmax);
r(end+1) = -980124;

cost = entropy(r(:))*length(r);
coef = r(:);
border = [];

clear r; %We won't need it anymore and it's better to free the memory...
clear lap;

liste = level;
border = [];

if(level == lmax)
    return;
end;

[listec1, coefc1, borderc1] = recurs_testdst2d(signal(1:m2,1:n2), level+1, lmax, rate);
[listec2, coefc2, borderc2] = recurs_testdst2d(signal(m2:m,1:n2), level+1, lmax, rate);
[listec3, coefc3, borderc3] = recurs_testdst2d(signal(1:m2,n2:n), level+1, lmax, rate);
[listec4, coefc4, borderc4] = recurs_testdst2d(signal(m2:m,n2:n), level+1, lmax, rate);


%Go with the inside coefficients
coef2 = [coefc1(:);coefc2(:);coefc3(:);coefc4(:)];
cost2 = entropy(coef2)*length(coef2);
%Done with inside coefficients


%Deal with the borders
b1 = signal(:,n2);
b2 = signal(m2,:);

eb1 = b1(2:m2-1) - (b1(1) + [1:m2-2]' / (m2-1) * (b1(m2)-b1(1)));
eb1 = dst(eb1)*sqrt(2/(m2-1));

eb2 = b1(m2+1:m-1) - (b1(m2) + [1:m2-2]' / (m2-1) * (b1(m)-b1(m2)));
eb2 = dst(eb2)*sqrt(2/(m2-1));

eb3 = b2(2:n2-1) - (b2(1) + [1:n2-2] / (n2-1) * (b2(n2)-b2(1)));
eb3 = dst(eb3)*sqrt(2/(n2-1));

eb4 = b2(n2+1:n-1) - (b2(n2) + [1:n2-2] / (n2-1) * (b2(n)-b2(n2)));
eb4 = dst(eb4)*sqrt(2/(n2-1));



eb1 = quantizeborder(eb1,rate, 1,0);
% indmax = max(find(eb1~=0));
% eb1 = eb1(1:indmax);
% eb1(end+1) = -980124;

eb2 = quantizeborder(eb2,rate, 1,0);
% indmax = max(find(eb2~=0));
% eb2 = eb2(1:indmax);
% eb2(end+1) = -980124;

eb3 = quantizeborder(eb3,rate, 0,0);
% indmax = max(find(eb3~=0));
% eb3 = eb3(1:indmax);
% eb3(end+1) = -980124;

eb4 = quantizeborder(eb4,rate, 0,0);
% indmax = max(find(eb4~=0));
% eb4 = eb4(1:indmax);
% eb4(end+1) = -980124;



border2 = [ borderc1(:); borderc2(:); borderc3(:); borderc4(:); eb1(:); eb2(:); eb3(:); eb4(:)];

cb2 = entropy(border2(:))*length(border2(:));
%Done with the borders

disp([ 'Level ' num2str(level) ', cost parent=' num2str(cost) ', cost enfants=' num2str(cost2) ', cost borders=' num2str(cb2)]);

if(cost2+cb2 < cost)
    coef = coef2;
    liste = [listec1 listec2 listec3 listec4];
    border = border2;
end



function e=entropy(coef)
mi = round(min(coef(:)));
ma = round(max(coef(:)));

if(ma-mi < 3)
    ma = mi+3;
end

[n,v] = hist(coef(:),mi:ma);

ne = n;
nrle = n;

ne = ne(find(n~=0));
netot = sum(ne);

if(ne(1) == netot)
    e = 0;
else
    ne = ne/netot;    
    e = -sum(ne.*log2(ne));
end
