Implementation of Landweber iteration for dynamic Radon Transform in MATLAB

Question

I am working on Dynamic Inverse Problems and Radon transforms. By the term "Dynamic," I mean the object (which eventually I am going to reconstruct its image by having the its measurements), deforms as time evolves. You may imagine the Heart or the Lungs in the human body. Therefore, the usual image reconstruction tools (like filtered back-projections) need to be modified to handle the deformation. In my case, the deformation is given by a smooth function in R^2,

φ^t(x,y)=(φ(t,x),φ(t,y))=(x*(1+c*x), y*(1+d*y))

where it denotes which particle is at position(x,y)at the time instancet. Here the coefficients canddare given by

c=(5.*sin(5e-5.*t*p0./pi))^(1/4), d= (5.*sin(7e-5.*t*p0./pi))^(1/4)

Now consider the following dynamic inverse operator (generalized Radon transform):

Af(s,t)=\iint_{R^2} f(φ^t(x,y)) \delta(x*cos t + y*sin t -s ) dx dy.

Here\iint_{R^2}is double integral overR^2 and \delta is the usual delta function. Using a change of variable, above operator can be written as

g(s,t)=Af(s,t)=\iint_{R^2} J(t,x,y)f(x,y)*\delta(φ^{-1}(t,x)*cos t+φ^{-1}(t,y)*sin t-s)dxdy

where the inverse functionsφ^{-1}are given by

φ^{-1}(t,x) = (\sqrt{1+4c*x}-1)/(2c), φ^{-1}(t,y)= (\sqrt{1+4d*y}-1)/(2d)

and the Jacobian

J(t,x,y)=1/(\sqrt{(1+4c*x)(1+4d*y)).

Heret is a parameter which determines both time and the direction of ray (the angle). Now the functionf(φ^t (x,y))shows the state of true object (Original Image which we call it Shepp-Logan SL) at the time instancet. For example whent=0, then the functionφ^t is just identity, so we get f(x,y) which is the original image SL.

I want to do the image reconstruction with Landweber iterations as it seems it is not expensive to run compared to other methods. A Landweber iterations takes the following form

f^{k+1}(x,y)= f^k(x,y)+ \gamma*A^T(g-Af^k)(x,y),

where\gammais a step length parameter let say 10^{-3}. Here by A^T, I mean the adjoint of the operatorA which is known as back projection. Therefore for a functionh, A^Thwill take the following form

A^Th(x,y)=\int_{R} J(t,x,y)*h(φ^{-1}(t,x)*cos t+φ^{-1}(t,y)*sin t,t)d t.

Therefore, if we replace h by g-Af^k in above adjoint equation and plug into the Landweber iteration, we get

f^{k+1}(x,y)= f^k(x,y)+\gamma*\int_{R} J(t,x,y)*(g-Af^k)(φ^{-1}(t,x)*cost+φ^{-1}(t,y)*sint,t)dt. 

Here is my goal. Let say we deform the original object SL from t=0 to t=5. Then we take the Radon of the deformed SL at t=5 which gives us the measurement at t=5. My friend and I have written a MATLAB code where we can deform the original object SL successfully and then we take the Radon transform (Forward problem). Now my goal is to reconstruct the original object SL at t=0 by performing the Landweber iteration formula (given above) from the measurement of deformed SL at t=5. We also have written the code for this part, but we cannot reconstruct the original SL. Here are the codes:

The original Object Shepp Logan SL

function SLAmir = SLAmirFunc(Xsize)
% this function creates a Shepp-Logan like image
% xsize: is the desired size of the image
imageSizeX = Xsize;%100;
imageSizeY = Xsize;
[columnsInImage rowsInImage] = meshgrid(1:imageSizeX, 1:imageSizeY);
centerX = floor(Xsize/2);%50;
centerY = floor(Xsize/2);%50;
radius1 = floor(Xsize/2)-floor(Xsize/10);%40;
radius2 = floor(Xsize/2)-floor(Xsize/5.5);%32;
circlePixels_1 = (rowsInImage - centerY).^2 ./radius1.^2 ...
+ (columnsInImage - centerX).^2 ./radius2.^2<= 1;
centerX = floor(Xsize/2);
centerY = floor(Xsize/2);
radius1 = floor(Xsize/2)-floor(Xsize/7);%36;
radius2 = floor(Xsize/2)-floor(Xsize/5);%30;
circlePixels_2 = (rowsInImage - centerY).^2 ./radius1.^2 ...
+ (columnsInImage - centerX).^2 ./radius2.^2<= 1;
centerX = floor(Xsize/2);%50;
centerY = floor(Xsize/2)-floor(Xsize/7);%36;
radius = floor(centerY/3);%12;
circlePixels_3 = (rowsInImage - centerY).^2 ...
+ (columnsInImage - centerX).^2 <= radius.^2;
centerX = floor(Xsize/2)+floor(Xsize/7);%64;
centerY = floor(Xsize/2)+floor(Xsize/16);%56;
radius1 = floor(centerY/4.5);%12;
radius2 = floor(centerY/7);%8;
circlePixels_4 = (rowsInImage - centerY).^2 ./radius1.^2 ...
+ (columnsInImage - centerX).^2 ./radius2.^2<= 1;
circlePixels_4=imrotate(circlePixels_4,-30,'crop','bilinear');
centerX = floor(Xsize/2)-floor(Xsize/12.5);%42;
centerY = floor(Xsize/2)-floor(Xsize/10);%40;
radius1 = floor(centerY/2.5);%16;
radius2 = floor(centerY/4);%10;
circlePixels_5 = (rowsInImage - centerY).^2 ./radius1.^2 ...
+ (columnsInImage - centerX).^2 ./radius2.^2<= 1;
circlePixels_5=imrotate(circlePixels_5,45,'crop','bilinear');
circlePixels_6=circlePixels_1-.8.*circlePixels_2 + 0.2.*circlePixels_3 
+.4.*circlePixels_4 -.1.*circlePixels_5;
sqIy=floor(Xsize/2)-floor(Xsize/10);
sqIx=floor(Xsize/2)+floor(Xsize/5.6);
circlePixels_6 (sqIx:sqIx+4,sqIy:sqIy+4)=circlePixels_6 
(sqIx:sqIx+4,sqIy:sqIy+4)+.15;
SLAmir=circlePixels_6;

Forward deformation:

function PPN3 = NonAffineTransform_v4(P,t,p0)
%performs the non-affine transform
%t:angle
%size_p:size of the main image +50
size_p=max(size(P));
xinit=(size_p-1)/200;
xx=-xinit:0.01:xinit;
yy=xx;
k=1;
for i=1:max(size(P))
for j=1:max(size(P))
    xxO(k,:)=[xx(i),yy(j)];
    k=k+1;
end
end
k=1;
cx = (5.*sin(5e-5.*t*p0./pi))^(1/4);
cy = (5.*sin(7e-5.*t*p0./pi))^(1/4);
for i=1:max(size(P))
for j=1:max(size(P))
  xxt(k,:)= [xx(i).*((cx.*xx(i))^0+(cx.*xx(i))^1 +(cx.*xx(i))^2 +
  (cx.*xx(i))^3+(cx.*xx(i))^4), ...%+(cx.*xx(i))^2+(cx.*xx(i))^3+
  (cx.*xx(i))^4
             yy(j).*((cy.*yy(j))^0+(cy.*yy(j))^1 +(cy.*yy(j))^2 +
  (cy.*yy(j))^3+(cy.*yy(j))^4)];%+(cy.*yy(j))^2+(cy.*yy(j))^3+(cy.*yy(j))^4
  k=k+1;
end
end
PPN=zeros(size_p,size_p);
mxx1=(max(size(xx))-1);
myy1=(max(size(yy))-1);
xx_idx=round(xxO*100+floor(mxx1/2)+1);% idx correspond to original image
xxt_idx=round(xxt*100+floor(mxx1/2)+1);%idx correspond to the deformed image
for i=1:max(size(xx))*max(size(yy))
if (xxt_idx(i,1)>0 && xxt_idx(i,2)>0) && (xxt_idx(i,1)<max(size(xx))+1 && 
xxt_idx(i,2)<max(size(xx))+1)
PPN(xxt_idx(i,1),xxt_idx(i,2))=P(xx_idx(i,1),xx_idx(i,2));%.*NormPhi(i)
end
end
PPN1=PPN;
kk=0;
loopCnt=1;
while kk ~=max(size(PPN1)) 
 for j=1:max(size(xx))
     aa=find(PPN1(:,j));
     aams=max(size(aa));
         if (isempty(aa)~=1)
             AA=PPN1(aa(1):aa(aams),j);
             bb = find(AA==0);
             if (isempty(bb)~=1)
                 for i=1:max(size(bb))
                     PPN1(bb(i)+aa(1)-1,j)=PPN1(bb(i)+aa(1)-2,j);
                 end
             else
                 kk=kk+1;
             end
         else
             kk=kk+1;
         end
 end
 if j==max(size(xx)) && kk ~=max(size(PPN1))
     kk=0;
     loopCnt=loopCnt+1;
 end
 end
 % 
 PPN2=PPN1;
 kk=0;
 loopCnt=1;
 while kk ~=max(size(PPN2)) 
 for j=1:max(size(xx))
     aa=find(PPN2(j,:));
     aams=max(size(aa));
         if (isempty(aa)~=1)
             AA=PPN2(j,aa(1):aa(aams));
             bb = find(AA==0);
             if (isempty(bb)~=1)
                 for i=1:max(size(bb))
                     PPN2(j,bb(i)+aa(1)-1)=PPN2(j,bb(i)+aa(1)-2);
                 end
             else
                 kk=kk+1;
             end
         else
             kk=kk+1;
         end
 end
 if j==max(size(xx)) && kk ~=max(size(PPN2))
     kk=0;
     loopCnt=loopCnt+1;
 end
 end 
 PPN3=PPN2;%.*NormPhi;%./sum(NormPhi(:));

Inverse deformation:

function fxy = invNonAffineTransform_v41(g,sp,t,tt,x,y,dxx)
p0=300;
ts=max(size(t));
xs=max(size(x));
ys=max(size(y));
temp=0;
mxx1=(max(size(x))-1);
myy1=(max(size(y))-1);
fxy=zeros(xs,ys);
for i=1:xs
for j=1:ys
    for k=1:ts
        cx = (5.*sin(5e-5.*t(k)*p0./pi))^(1/4);%(t(k)+pi/2)
        cy = (5.*sin(7e-5.*t(k)*p0./pi))^(1/4);
        invphix = (sqrt(1+4.*cx.*x(i))-1)./(2.*cx);
        invphiy = (sqrt(1+4.*cy.*y(j))-1)./(2.*cy);
        Stemp = invphix.*cos(t(k)+pi/2)+ invphiy.*sin(t(k)+pi/2);
        %StempIdx= floor(Stemp./dxx);
        StempIdx = floor(Stemp.*100);
        spIdx=find(sp==StempIdx);
        J=sqrt(1/((1+2.*cx.*x(i)).*(1+2.*cy.*y(j))));
        if isempty(spIdx)
            temp=temp+0;
        else
            temp = temp +J.*g(spIdx,k);
        end
    end
     fxy(i,j) = fxy(i,j)+ tt.*temp;
     temp =0;
end
end

Run:

clear all
clc
close all
size_p=151;
P=zeros(size_p,size_p);
DisSide=40;
P1 = SLAmirFunc(size_p-DisSide);%
P(DisSide/2+1:size_p-DisSide/2,DisSide/2+1:size_p-DisSide/2)=P1;
p0=300;
size_p=max(size(P));
xinit=(size_p-1)/200;
xx=-xinit:0.01:xinit;
yy=xx;
j=1;
tt=[25:5:300].*pi./p0;
fxy=P;
for i=1:max(size(tt))
PPN3p = NonAffineTransform_v4(fxy,tt(i),p0);
PPN3 = PPN3p;
figure;imagesc(PPN3), colormap(bone), colorbar;
end
ttdeg=tt.*180./pi;
[g,sp] = radon(PPN3,ttdeg);
I1 = iradon(g,ttdeg,size_p);
Afpre=zeros(max(size(sp)),max(size(tt)));
fold=zeros(size_p,size_p);
stepsize=0.0001;
ssp=max(size(sp));
dxx=(xx(size_p)-xx(1)+1)/ssp;
Dt=tt(2)-tt(1);
for it=1:10
g1=g-Afpre;
fnew = fold + stepsize.*invNonAffineTransform_v41(g1,sp,tt,Dt,xx,yy,dxx);
for i=1:max(size(tt))
PPN3p = NonAffineTransform_v4(fnew,tt(i),p0);
clear PPN3
PPN3 = PPN3p;   
end
[Afpre,sp] = radon(PPN3,ttdeg);
fold = fnew;
end

I was wondering if anyone here is familiar with this topic and can be any help.

If this is not clear please comment, and I will be answering your questions. Thank you in advance.