matlab - 在MATLAB中实现动态Radon变换的Landweber迭代

Question

我正在研究动态逆问题和氡变换。术语“动态”是指对象（最终我将通过测量来重建其图像）随着时间的推移而变形。你可以想象人体内的心脏或肺。因此，需要修改通常的图像重建工具（如过滤后的反投影）来处理变形。就我而言，变形由中的平滑函数给出R^2，

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

其中它表示哪个粒子在(x,y)时间 instance的位置t。这里的系数c和d由下式给出

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

现在考虑以下动态逆算子（广义 Radon 变换）：

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

这\iint_{R^2}是双重积分R^2，\delta是通常的delta函数。使用变量的变化，上面的运算符可以写成

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

其中反函数φ^{-1}由下式给出

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

和雅可比

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

这t是一个决定时间和光线方向（角度）的参数。现在该函数f(φ^t (x,y))显示时间实例的真实对象（我们称之为 Shepp-Logan SL 的原始图像）的状态t。例如当t=0，那么函数φ^t就是恒等式，所以我们得到f(x,y)哪个是原始图像 SL。

我想用 Landweber 迭代进行图像重建，因为与其他方法相比，运行它似乎并不昂贵。Landweber 迭代采用以下形式

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

哪里\gamma是步长参数让我们说10^{-3}。在这里A^T，我的意思是算子的伴随，A称为反投影。因此对于一个函数h，A^Th将采用以下形式

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

因此，如果我们在上面的伴随方程中替换h为g-Af^k并代入 Landweber 迭代，我们得到

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.

这是我的目标。假设我们将原始对象 SL 从变形t=0为t=5。然后我们取变形 SL 的氡气，在处t=5给出测量值t=5。我和我的朋友编写了一个 MATLAB 代码，我们可以在其中成功地对原始对象 SL 进行变形，然后进行 Radon 变换（正向问题）。t=0现在我的目标是通过执行 Landweber 迭代公式（如上所示）从变形 SL 处的测量中重建原始对象 SL处t=5。我们也为这部分编写了代码，但我们无法重建原始 SL。以下是代码：

最初的 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;

前向变形：

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(:));

逆变形：

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

跑：

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

我想知道这里是否有人熟悉这个主题并且可以提供任何帮助。

如果不清楚，请发表评论，我将回答您的问题。先感谢您。

matlab - 在MATLAB中实现动态Radon变换的Landweber迭代

0 回答 0

Related

Reference