%------------------------------------------------
% decomposition de l'image en 4 sous-images
%------------------------------------------------

clear all; close all;
load escher;
load map;

[nl,nc] = size(IM);
IM = double(IM); 
x = IM;
figure('Position',[50 600 560 420]); colormap(map);
image(IM);
[h0,h1,g0,g1] = litfiltres;

% calcul du retard
zk = conv(h0,g0)+conv(h1,g1);
k = find(abs(zk)>1e-12)-1;

npix = nl*nc;

%------------------------------------------------
% analyse
%------------------------------------------------

% filtrage et decimation des lignes de x
x = [x,zeros(nl,k)]; nc = nc+k;
x0 = filter(h0,1,x'); x0 = x0';
x1 = filter(h1,1,x'); x1 = x1';
x0 = x0(:,1:2:nc);
x1 = x1(:,1:2:nc);

% filtrage et decimation des colonnes de x0, x1
nc2 = size(x0,2);
x0 = [x0;zeros(k,nc2)]; nl = nl+k;
x00 = filter(h0,1,x0);
x01 = filter(h1,1,x0);
x1 = [x1;zeros(k,nc2)];
x10 = filter(h0,1,x1);
x11 = filter(h1,1,x1);

x00 = x00(1:2:nl,:);
x01 = x01(1:2:nl,:);
x10 = x10(1:2:nl,:);
x11 = x11(1:2:nl,:);

%------------------------------------------------
% visualisation du resultat intermediaire
%------------------------------------------------

figure('Position',[650 600 560 420]); colormap(map);
x = [x00,recale(x10);recale(x01),recale(x11)];
image(x); 

%------------------------------------------------
% synthese
%------------------------------------------------

% interpolation et filtrage des colonnes
z = zeros(nl,nc2); z(1:2:nl,:) = x00; x00 = z;
z = zeros(nl,nc2); z(1:2:nl,:) = x01; x01 = z;
z = zeros(nl,nc2); z(1:2:nl,:) = x10; x10 = z;
z = zeros(nl,nc2); z(1:2:nl,:) = x11; x11 = z;
x0 = filter(g0,1,x00)+filter(g1,1,x01);
x1 = filter(g0,1,x10)+filter(g1,1,x11);

% interpolation des lignes et filtrage lignes
z = zeros(nl,nc); z(:,1:2:nc) = x0; x0 = z;
z = zeros(nl,nc); z(:,1:2:nc) = x1; x1 = z;
x = filter(g0,1,x0')+filter(g1,1,x1'); 
x = x';

x = x(k+1:nl,k+1:nc);
figure('Position',[1250 600 560 420]); colormap(map);
image(x);
er1 = max(max(abs(x-IM)));
er2 = sum(sum(abs(x-IM)))/npix;
disp([er1 er2]);

