%Computation of 2-D Haar transform of Shepp-Logan Phantom:
clear;load shepp256.mat;imagesc(X),colormap(gray)
ZLL0=X;%1st stage decomposition
ZLL1=ZLL0(1:2:256,1:2:256)+ZLL0(1:2:256,2:2:256);
ZLL1=ZLL1+ZLL0(2:2:256,1:2:256)+ZLL0(2:2:256,2:2:256);
ZLH1=ZLL0(1:2:256,2:2:256)-ZLL0(1:2:256,1:2:256);
ZLH1=ZLH1+ZLL0(2:2:256,2:2:256)-ZLL0(2:2:256,1:2:256);
ZHL1=ZLL0(2:2:256,2:2:256)-ZLL0(1:2:256,2:2:256);
ZHL1=ZHL1+ZLL0(2:2:256,1:2:256)-ZLL0(1:2:256,1:2:256);
ZHH1=ZLL0(2:2:256,2:2:256)-ZLL0(1:2:256,2:2:256);
ZHH1=ZHH1-ZLL0(2:2:256,1:2:256)+ZLL0(1:2:256,1:2:256);
ZLL1=ZLL1/2;ZLH1=ZLH1/2;ZHL1=ZHL1/2;ZHH1=ZHH1/2;
%2nd stage decomposition
ZLL2=ZLL1(1:2:128,1:2:128)+ZLL1(1:2:128,2:2:128);
ZLL2=ZLL2+ZLL1(2:2:128,1:2:128)+ZLL1(2:2:128,2:2:128);
ZLH2=ZLL1(1:2:128,2:2:128)-ZLL1(1:2:128,1:2:128);
ZLH2=ZLH2+ZLL1(2:2:128,2:2:128)-ZLL1(2:2:128,1:2:128);
ZHL2=ZLL1(2:2:128,2:2:128)-ZLL1(1:2:128,2:2:128);
ZHL2=ZHL2+ZLL1(2:2:128,1:2:128)-ZLL1(1:2:128,1:2:128);
ZHH2=ZLL1(2:2:128,2:2:128)-ZLL1(1:2:128,2:2:128);
ZHH2=ZHH2-ZLL1(2:2:128,1:2:128)+ZLL1(1:2:128,1:2:128);
ZLL2=ZLL2/2;ZLH2=ZLH2/2;ZHL2=ZHL2/2;ZHH2=ZHH2/2;
%3rd stage decomposition
ZLL3=ZLL2(1:2:64,1:2:64)+ZLL2(1:2:64,2:2:64);
ZLL3=ZLL3+ZLL2(2:2:64,1:2:64)+ZLL2(2:2:64,2:2:64);
ZLH3=ZLL2(1:2:64,2:2:64)-ZLL2(1:2:64,1:2:64);
ZLH3=ZLH3+ZLL2(2:2:64,2:2:64)-ZLL2(2:2:64,1:2:64);
ZHL3=ZLL2(2:2:64,2:2:64)-ZLL2(1:2:64,2:2:64);
ZHL3=ZHL3+ZLL2(2:2:64,1:2:64)-ZLL2(1:2:64,1:2:64);
ZHH3=ZLL2(2:2:64,2:2:64)-ZLL2(1:2:64,2:2:64);
ZHH3=ZHH3-ZLL2(2:2:64,1:2:64)+ZLL2(1:2:64,1:2:64);
ZLL3=ZLL3/2;ZLH3=ZLH3/2;ZHL3=ZHL3/2;ZHH3=ZHH3/2;
%4th stage decomposition
ZLL4=ZLL3(1:2:32,1:2:32)+ZLL3(1:2:32,2:2:32);
ZLL4=ZLL4+ZLL3(2:2:32,1:2:32)+ZLL3(2:2:32,2:2:32);
ZLH4=ZLL3(1:2:32,2:2:32)-ZLL3(1:2:32,1:2:32);
ZLH4=ZLH4+ZLL3(2:2:32,2:2:32)-ZLL3(2:2:32,1:2:32);
ZHL4=ZLL3(2:2:32,2:2:32)-ZLL3(1:2:32,2:2:32);
ZHL4=ZHL4+ZLL3(2:2:32,1:2:32)-ZLL3(1:2:32,1:2:32);
ZHH4=ZLL3(2:2:32,2:2:32)-ZLL3(1:2:32,2:2:32);
ZHH4=ZHH4-ZLL3(2:2:32,1:2:32)+ZLL3(1:2:32,1:2:32);
ZLL4=ZLL4/2;ZLH4=ZLH4/2;ZHL4=ZHL4/2;ZHH4=ZHH4/2;
ZZ=[ZLL4 ZLH4;ZHL4 ZHH4];%Display as image:
ZZ=[ZZ ZLH3;ZHL3 ZHH3];
ZZ=[ZZ ZLH2;ZHL2 ZHH2];
ZZ=[ZZ ZLH1;ZHL1 ZHH1];
figure,imagesc(ZZ),colormap(gray)
%3764 nonzero XX;3713 nonzero ZZ. 
%Inverse Haar (note double reversals, & FAST):
WLL2(2:2:64,2:2:64)=(ZLL3+ZLH3+ZHL3+ZHH3)/2;
WLL2(2:2:64,1:2:64)=(ZLL3-ZLH3+ZHL3-ZHH3)/2;
WLL2(1:2:64,2:2:64)=(ZLL3+ZLH3-ZHL3-ZHH3)/2;
WLL2(1:2:64,1:2:64)=(ZLL3-ZLH3-ZHL3+ZHH3)/2;
%2nd stage:
WLL1(2:2:128,2:2:128)=(WLL2+ZLH2+ZHL2+ZHH2)/2;
WLL1(2:2:128,1:2:128)=(WLL2-ZLH2+ZHL2-ZHH2)/2;
WLL1(1:2:128,2:2:128)=(WLL2+ZLH2-ZHL2-ZHH2)/2;
WLL1(1:2:128,1:2:128)=(WLL2-ZLH2-ZHL2+ZHH2)/2;
%1st stage:
W(2:2:256,2:2:256)=(WLL1+ZLH1+ZHL1+ZHH1)/2;
W(2:2:256,1:2:256)=(WLL1-ZLH1+ZHL1-ZHH1)/2;
W(1:2:256,2:2:256)=(WLL1+ZLH1-ZHL1-ZHH1)/2;
W(1:2:256,1:2:256)=(WLL1-ZLH1-ZHL1+ZHH1)/2;
figure,imagesc(W),colormap(gray)
