
%Tomography from Partial Fourier Data using 
%2-D Haar wavelet transform and IST algorithm.
clear;load shepp256.mat;N=size(X,1);
M=6;JMAX=3;lambda=0.01;NMAX=500;
%Set up rays of known DFT values:
F(N/2,N/2)=0;FY(N,N)=0;FWK(N,N)=0;
for J1=1:JMAX;for J2=1:JMAX;F(J1,:)=1;F(:,J2)=1;
for theta=[1:M-1]/M*pi/2;for I=1:N/2;
F(round(I*cos(theta))+J1,round(I*sin(theta))+J2)=1;
end;end;end;end;F=F(1:N/2,1:N/2);IF=[N/2:-1:2];
FF=[F zeros(N/2,1) F(:,IF);zeros(1,N)];
FF=[FF;F(IF,:) zeros(N/2-1,1) F(IF,IF)];
imagesc(fftshift(FF)),sum(sum(FF))%Locations of known DFT
%max(max(imag(ifft(FF))))%confirm conjugate symmetric DFT
K=find(FF>0.5);FX=fft2(X);FY(K)=FX(K);%Known 2-D DFT values.
IFY=real(ifft2(FY));%Least-Squares Reconstruction
figure,imagesc(IFY),colormap(gray),title('LEAST-SQUARES')
%Initialize Landweber with all zeros:
%(Also preallocates variable spaces)
ZLL3=zeros(32,32);ZLH3=zeros(32,32);ZHL3=zeros(32,32);ZHH3=zeros(32,32);
ZLL2=zeros(64,64);ZLH2=zeros(64,64);ZHL2=zeros(64,64);ZHH2=zeros(64,64);
ZLL1=zeros(128,128);ZLH1=zeros(128,128);ZHL1=zeros(128,128);ZHH1=zeros(128,128);
for I=1:NMAX;
%Landweber:z^{k+1}=z^{k}+H'F'(y-FH'z^{k}).
%H=Haar. F=DFT. W=H'z^{k}=XHAT.
%Inverse Haar H'=H^{-1}:
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;
%FW=DFT(W) at known K. FD=(y-FH'z^{k}).
FW=fft2(W);FWK(K)=FW(K);FD=FY-FWK;
%HF'(y-FH'z^{k})=HF'(FD)=HD=H(DLL0).
DLL0=real(ifft2(FD,256,256));%Now Haar H of DLL0:
%1st stage decomposition
DLL1=DLL0(1:2:256,1:2:256)+DLL0(1:2:256,2:2:256);
DLL1=DLL1+DLL0(2:2:256,1:2:256)+DLL0(2:2:256,2:2:256);
DLH1=DLL0(1:2:256,2:2:256)-DLL0(1:2:256,1:2:256);
DLH1=DLH1+DLL0(2:2:256,2:2:256)-DLL0(2:2:256,1:2:256);
DHL1=DLL0(2:2:256,2:2:256)-DLL0(1:2:256,2:2:256);
DHL1=DHL1+DLL0(2:2:256,1:2:256)-DLL0(1:2:256,1:2:256);
DHH1=DLL0(2:2:256,2:2:256)-DLL0(1:2:256,2:2:256);
DHH1=DHH1-DLL0(2:2:256,1:2:256)+DLL0(1:2:256,1:2:256);
DLL1=DLL1/2;DLH1=DLH1/2;DHL1=DHL1/2;DHH1=DHH1/2;
%2nd stage decomposition
DLL2=DLL1(1:2:128,1:2:128)+DLL1(1:2:128,2:2:128);
DLL2=DLL2+DLL1(2:2:128,1:2:128)+DLL1(2:2:128,2:2:128);
DLH2=DLL1(1:2:128,2:2:128)-DLL1(1:2:128,1:2:128);
DLH2=DLH2+DLL1(2:2:128,2:2:128)-DLL1(2:2:128,1:2:128);
DHL2=DLL1(2:2:128,2:2:128)-DLL1(1:2:128,2:2:128);
DHL2=DHL2+DLL1(2:2:128,1:2:128)-DLL1(1:2:128,1:2:128);
DHH2=DLL1(2:2:128,2:2:128)-DLL1(1:2:128,2:2:128);
DHH2=DHH2-DLL1(2:2:128,1:2:128)+DLL1(1:2:128,1:2:128);
DLL2=DLL2/2;DLH2=DLH2/2;DHL2=DHL2/2;DHH2=DHH2/2;
%3rd stage decomposition
DLL3=DLL2(1:2:64,1:2:64)+DLL2(1:2:64,2:2:64);
DLL3=DLL3+DLL2(2:2:64,1:2:64)+DLL2(2:2:64,2:2:64);
DLH3=DLL2(1:2:64,2:2:64)-DLL2(1:2:64,1:2:64);
DLH3=DLH3+DLL2(2:2:64,2:2:64)-DLL2(2:2:64,1:2:64);
DHL3=DLL2(2:2:64,2:2:64)-DLL2(1:2:64,2:2:64);
DHL3=DHL3+DLL2(2:2:64,1:2:64)-DLL2(1:2:64,1:2:64);
DHH3=DLL2(2:2:64,2:2:64)-DLL2(1:2:64,2:2:64);
DHH3=DHH3-DLL2(2:2:64,1:2:64)+DLL2(1:2:64,1:2:64);
DLL3=DLL3/2;DLH3=DLH3/2;DHL3=DHL3/2;DHH3=DHH3/2;
%z^{k+1}=z^{k}+HF'(y-FH'z^{k}):Now add z^{k}:
ZLH1=ZLH1+DLH1;ZHL1=ZHL1+DHL1;ZHH1=ZHH1+DHH1;
ZLH2=ZLH2+DLH2;ZHL2=ZHL2+DHL2;ZHH2=ZHH2+DHH2;
ZLH3=ZLH3+DLH3;ZHL3=ZHL3+DHL3;ZHH3=ZHH3+DHH3;
ZLL3=ZLL3+DLL3;%Only use finest-scale average.
%Threshold & Shrink each wavelet component:
%Threshold small values to 0.
ZLH1(abs(ZLH1)<lambda)=0;ZHL1(abs(ZHL1)<lambda)=0;
ZLH2(abs(ZLH2)<lambda)=0;ZHL2(abs(ZHL2)<lambda)=0;
ZLH3(abs(ZLH3)<lambda)=0;ZHL3(abs(ZHL3)<lambda)=0;
ZHH1(abs(ZHH1)<lambda)=0;ZHH2(abs(ZHH2)<lambda)=0;
ZHH3(abs(ZHH3)<lambda)=0;ZLL3(abs(ZLL3)<lambda)=0;
%Shrink larger values by lambda:
ZLH1(abs(ZLH1)>lambda)=ZLH1(abs(ZLH1)>lambda)-lambda*sign(ZLH1(abs(ZLH1)>lambda));
ZLH2(abs(ZLH2)>lambda)=ZLH2(abs(ZLH2)>lambda)-lambda*sign(ZLH2(abs(ZLH2)>lambda));
ZLH3(abs(ZLH3)>lambda)=ZLH3(abs(ZLH3)>lambda)-lambda*sign(ZLH3(abs(ZLH3)>lambda));
ZHL1(abs(ZHL1)>lambda)=ZHL1(abs(ZHL1)>lambda)-lambda*sign(ZHL1(abs(ZHL1)>lambda));
ZHL2(abs(ZHL2)>lambda)=ZHL2(abs(ZHL2)>lambda)-lambda*sign(ZHL2(abs(ZHL2)>lambda));
ZHL3(abs(ZHL3)>lambda)=ZHL3(abs(ZHL3)>lambda)-lambda*sign(ZHL3(abs(ZHL3)>lambda));
ZHH1(abs(ZHH1)>lambda)=ZHH1(abs(ZHH1)>lambda)-lambda*sign(ZHH1(abs(ZHH1)>lambda));
ZHH2(abs(ZHH2)>lambda)=ZHH2(abs(ZHH2)>lambda)-lambda*sign(ZHH2(abs(ZHH2)>lambda));
ZHH3(abs(ZHH3)>lambda)=ZHH3(abs(ZHH3)>lambda)-lambda*sign(ZHH3(abs(ZHH3)>lambda));
ZLL3(abs(ZLL3)>lambda)=ZLL3(abs(ZLL3)>lambda)-lambda*sign(ZLL3(abs(ZLL3)>lambda));
end;%Final Inverse Haar to get XHAT from Z:
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),title('IST RESULT')
