%
% Data una sequenza bi-dimensionale, la trasformata wavelet puo' essere prodotta in due passi,
% utilizzando l'algoritmo 1D ed agendo al primo passo sulle righe (o sulle colonne) ed al secondo
% sulle colonne (o sulle righe). La sequenza trasformata possiede dunque una suddivisione a blocchi
% del tipo:                  ThTv    FhTv
%                            ThFv    FhFv
% in cui ogni blocco è dunque una combinazione di trend e fluttuazioni orrizzontali e verticali.
%
% I livelli successivi della trasformata procedono considerando la sequenza trend "ThTv".
%
%

%
% costruiamo un'immagine a ottagono:
%
Nx=128;
Ny=Nx;
otta = ones(Nx,Ny);
for x=Nx/2+(-Nx/4:Nx/4)
  for y=Ny/2+(-Ny/4:Ny/4)
    if ((x-Nx/4+y-Ny/4 > Nx/6) & (3*Nx/4-x+y-Ny/4 > Nx/6) & (x-Nx/4+3*Ny/4-y > Nx/6) & (3*Nx/4-x+3*Ny/4-y > Nx/6))
      otta(x,y) = 255;
    end;
  end;
end;
figure(1); clf, colormap('gray'); image(otta); pause,
%
figure(11); clf, mesh(otta); pause,

%
% calcolo del primo livello della trasformata wavelet 2D
%
x = otta;
% scaling numbers e wavelet numbers Coif6 :
scalnum = [1-sqrt(7) 5+sqrt(7) 14+2*sqrt(7) 14-2*sqrt(7) 1-sqrt(7) -3+sqrt(7)]/(16*sqrt(2));
wavenum = [-3+sqrt(7) -(1-sqrt(7)) 14-2*sqrt(7) -(14+2*sqrt(7)) 5+sqrt(7) -(1-sqrt(7))]/(16*sqrt(2));
for ix=1:Nx
  totflutt = [];
  trend = x(ix,:);
  trend = filter(scalnum,[1],x(ix,:));
  trend = trend(2:2:end);
  flutt = filter(wavenum,[1],x(ix,:));
  flutt = flutt(2:2:end);
  x(ix,:) = [trend flutt];
end;
for iy=1:Ny
  totflutt = [];
  trend = x(:,iy);
  trend = filter(scalnum,[1],x(:,iy));
  trend = trend(2:2:end);
  flutt = filter(wavenum,[1],x(:,iy));
  flutt = flutt(2:2:end);
  x(:,iy) = [trend; flutt];
end;
figure(2); clf, colormap('gray'); image(x); title('Trasformata wavelet 2D'); pause,


%
% ora, azzeriamo la sotto-sequenza trend "ThTv" e calcoliamo la trasformata inversa
% 
trasf_w = x;    %  memorizzo la trasformata wavelet
x(1:Nx/2,1:Ny/2) = zeros(Nx/2,Ny/2);
figure(3); clf, colormap('gray'); image(x); pause,
%
for ix=1:Nx
  ttemp = filter(scalnum(end:-1:1),[1], upsample(x(ix,1:Ny/2),2) );
  ttemp = ttemp + filter(wavenum(end:-1:1),[1], upsample(x(ix,Ny/2+1:Ny),2) );
  x(ix,:) = ttemp;
end;
for iy=1:Ny
  ttemp = filter(scalnum(end:-1:1),[1], upsample(x(1:Nx/2,iy)',2) );
  ttemp = ttemp + filter(wavenum(end:-1:1),[1], upsample(x(Nx/2+1:Nx,iy)',2) );
  x(:,iy) = ttemp';
end;
figure(4); clf, colormap('gray'); image(x); title('edge detection'), 
%
% si può notare come le fluttuazioni mettano bene in evidenza i bordi (edges); 
%

