%This computes 2-D Haar transform of Shepp-Logan phantom.
%and depicts it as in Fig. 7-16 (thumbnail at upper left).
%EXTRA: Reconstruct Shepp-Logan phantom from its wavelet transform.
%Note how fast this program runs (even 2-D wavelets are fast).
clear;load shepp256.mat;%256X256 Shepp-Logan phantom.
figure,imagesc(X),colormap(gray),axis off,title('Shepp-Logan Phantom')
%Compute 2-D Haar wavelet transform:
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;
%Display 2-D Haar transform as in Fig. 7-16.
%NOTE: Fig. 7-16 was enhanced for display.
ZZ=[ZLL4 ZLH4;ZHL4 ZHH4];ZZ=[ZZ ZLH3;ZHL3 ZHH3];
ZZ=[ZZ ZLH2;ZHL2 ZHH2];ZZ=[ZZ ZLH1;ZHL1 ZHH1];
figure,imagesc(ZZ),colormap(gray),axis off,title('2-D Haar transform of Shepp-Logan')
%EXTRA: Reconstruct Shepp-Logan phantom from its 2-D Haar transform.
%Inverse Haar transform (note the double spatial reversals):
%3rd stage reconstruction:
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 reconstruction:
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 reconstruction:
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),axis off,title('Reconstructed Shepp-Logan phantom')
