% Ovo je skripta koja daje primjer mnozenja matrica upotrebom wavelet
% transformacije. Matrica s kojom radim predstavlja sliku velicine 256x256.
% Matricu cu pomnoziti samu sa sobom i analizirati koja je razlika u
% mnozenju matrica na klasican nacin i upotrebom wavelet transformacije.

% Ucitavanje slike
[img, map] = imread('slika.tiff');
whos;

% Ucitana slika je indeksirana slika, cije detaljne karakteristike nisu
% bitne.
figure;
image(img);
colormap(gray(255)); % Prikaz sive slike s paletom od 256 nijansi.

% Klasicno mnozenje matrica.
a = double(img); % Matrica koja predstavlja sliku je tipa uint8, za koji nije definiran operator (*)
b = a*a;
figure;
imagesc(b); % Rezultat skaliramo i prikazujemo.
colormap(gray(255));

% Na matricu cemo primjeniti tzv. psi^n wavelet transformaciju. Koristiti cu 
% Haarov wavelet kao funkciju baze. Kako je matrica velicine NxN (N=256), razlaganje 
% cu vrsiti na n=8 razina, (2^n=N), buduci da razlaganja na vise razina ne bi imalo smisla.
a = img; 
[all, alh, ahl, ahh] = DWT2(a, 'haar');
[bll, blh, bhl, bhh] = DWT2(b, 'haar'); % Transformacija umnoska.
% all je usrednjena komponenta slike;
% alh je vertikalna komponenta slike;
% ahl je horizontalna komponenta slike;
% ahh je dijagonalna komponenta slike;
% Sve komponente zajedno slozene u matricu cine 2D wavelet transformaciju
% slike.
% Preslazemo komponente (psi^1):

apsi = [all, ahl; alh, ahh]; % Vertikalna i horizontalna komponenta se zamjenjuju.
bpsi = [bll, bhl; blh, bhh];

% apsi je psi^n transformacija od matrice a na prvoj razini. Na apsi primjenjujemo
% opet wavelet transformaciju i nakon preslagivanja komponenti dobivamo
% razlaganje na drugoj razini (psi^2):

[all, alh, ahl, ahh] = DWT2(apsi, 'haar');
apsi2 = [all, ahl; alh, ahh];
[bll, blh, bhl, bhh] = DWT2(bpsi, 'haar');
bpsi2 = [bll, bhl; blh, bhh];

% Buduci je n=8, razlazemo do psi^8.
[all, alh, ahl, ahh] = DWT2(apsi2, 'haar');
apsi3 = [all, ahl; alh, ahh];
[bll, blh, bhl, bhh] = DWT2(bpsi2, 'haar');
bpsi3 = [bll, bhl; blh, bhh];

[all, alh, ahl, ahh] = DWT2(apsi3, 'haar');
apsi4 = [all, ahl; alh, ahh];
[bll, blh, bhl, bhh] = DWT2(bpsi3, 'haar');
bpsi4 = [bll, bhl; blh, bhh];

[all, alh, ahl, ahh] = DWT2(apsi4, 'haar');
apsi5 = [all, ahl; alh, ahh];
[bll, blh, bhl, bhh] = DWT2(bpsi4, 'haar');
bpsi5 = [bll, bhl; blh, bhh];

[all, alh, ahl, ahh] = DWT2(apsi5, 'haar');
apsi6 = [all, ahl; alh, ahh];
[bll, blh, bhl, bhh] = DWT2(bpsi5, 'haar');
bpsi6 = [bll, bhl; blh, bhh];

[all, alh, ahl, ahh] = DWT2(apsi6, 'haar');
apsi7 = [all, ahl; alh, ahh];
[bll, blh, bhl, bhh] = DWT2(bpsi6, 'haar');
bpsi7 = [bll, bhl; blh, bhh];

[all, alh, ahl, ahh] = DWT2(apsi7, 'haar');
apsi8 = [all, ahl; alh, ahh];
[bll, blh, bhl, bhh] = DWT2(bpsi7, 'haar');
bpsi8 = [bll, bhl; blh, bhh];

% apsi8 je psi^8 wavelet transformacija ucitane slike img.
% bpsi8 je psi^8 wavelet transformacija umnoska.
% slijedi provjera da li vrijedi W(x*y) = W(x)*W(y) za ovakvu wavelet
% transformaciju: bpsi8 = apsi8*apsi8.

provjera = max(max(abs(bpsi8-apsi8*apsi8)));
cpsi8 = apsi8*apsi8;

% Prikaz tarnsformirane slike.
figure;
image(apsi8);
colormap(gray(255));

% Rekonstrukcija umnoska u wavelet domeni iz wavelet domene.
% Prije inverzne transformacije definiraju se komponente matrice.
cll = cpsi8(1:128, 1:128);
chl = cpsi8(1:128, 129:256);
clh = cpsi8(129:256, 1:128);
chh = cpsi8(129:256, 129:256);
cpsi7 = IDWT2(cll, clh, chl, chh, 'haar');

cll = cpsi7(1:128, 1:128);
chl = cpsi7(1:128, 129:256);
clh = cpsi7(129:256, 1:128);
chh = cpsi7(129:256, 129:256);
cpsi6 = IDWT2(cll, clh, chl, chh, 'haar');

cll = cpsi6(1:128, 1:128);
chl = cpsi6(1:128, 129:256);
clh = cpsi6(129:256, 1:128);
chh = cpsi6(129:256, 129:256);
cpsi5 = IDWT2(cll, clh, chl, chh, 'haar');

cll = cpsi5(1:128, 1:128);
chl = cpsi5(1:128, 129:256);
clh = cpsi5(129:256, 1:128);
chh = cpsi5(129:256, 129:256);
cpsi4 = IDWT2(cll, clh, chl, chh, 'haar');

cll = cpsi4(1:128, 1:128);
chl = cpsi4(1:128, 129:256);
clh = cpsi4(129:256, 1:128);
chh = cpsi4(129:256, 129:256);
cpsi3 = IDWT2(cll, clh, chl, chh, 'haar');

cll = cpsi3(1:128, 1:128);
chl = cpsi3(1:128, 129:256);
clh = cpsi3(129:256, 1:128);
chh = cpsi3(129:256, 129:256);
cpsi2 = IDWT2(cll, clh, chl, chh, 'haar');

cll = cpsi2(1:128, 1:128);
chl = cpsi2(1:128, 129:256);
clh = cpsi2(129:256, 1:128);
chh = cpsi2(129:256, 129:256);
cpsi = IDWT2(cll, clh, chl, chh, 'haar');

cll = cpsi(1:128, 1:128);
chl = cpsi(1:128, 129:256);
clh = cpsi(129:256, 1:128);
chh = cpsi(129:256, 129:256);
c = IDWT2(cll, clh, chl, chh, 'haar');

% Prikaz umnoska dobivenog iz wavelet transformacije.
figure;
imagesc(c);
colormap(gray(255));

% Pogreska mnozenja.
pogreska = max(max(abs(b-c)));

