CT图像重建算法——等角扇束重排算法重

blue的Whale鲸鱼

于 2024-06-08 23:12:09 发布

阅读量648

点赞数 12

分类专栏： CT图像重建算法文章标签：算法 matlab 重构计算机视觉

本文链接：https://blog.csdn.net/2201_76025250/article/details/139552564

版权

CT图像重建算法专栏收录该内容

4 篇文章 0 订阅

订阅专栏

扇束投影重建算法大致分为两类:一类是重排算法，即将采集到的扇形数据重新排列成平行的射线投影数据，再用实验2所介绍的平行束算法重建;另一类算法是扇束投影直接重建算法，即根据扇形束投影数据本身的特点直接进行图像重建的方法。
根据扇束扫描检测器的布置形式，扇形射线的分布形式有两种:等角度分布的扇形束(等角扇束)和等间距分布的扇形束(等距扇束)。在等角扇束的模式下,检测器单元可以是直线分布，也可以是圆弧分布，如图3.1(a)和(b)所示。在等距扇束的模式下,检测器单元在直线上是作等距分布的，但在这种情况下实现的分布是不等角度的，如图3.2所示。

所所谓重排，是把扇束情况下得到的全部投影数据 $p(\gamma ,\displaystyle \beta )$ 重新排列整理成为不同视角下的平行射线投影数据，再利用平行束投影算法进行重建。前面已知，一条射线可以由 $(\gamma ,\displaystyle \beta )$ 表示，同时也可以看作平行射线中的一条,用（t,0）表示。将扇形射线重组成平行射线的示意图如图3.20所示。

算法的具体计算方法见《医学断层图像重建仿真实验》书籍黄力宇等著

MATLAB实现代码如下：

%扇束等角重排算法
clc;
clear all;
close all;
delta_beta=1;%旋转角度增量
beta=0:delta_beta:359;%旋转角度
N=256;%图像大小
N_d=257;%测器通道个数
SOD=250;%焦距
delta_gamma=0.25;%扇束张角增量
I=phantom(N);%建立Shepp-logan 头模型
P=FanbeamAFProj(N,beta,SOD,N_d,delta_gamma);%投影数据仿真
rec_RL=FanBeamAngleResorting(P,N,SOD,delta_beta,delta_gamma);%调用等角扇束重建函数进行重建
figure;
subplot(1,2,1),imshow(I,[]),title('原始图像');
subplot(1,2,2),imshow(rec_RL,[]),title('重排算法重建的图像');


function P=FanbeamAFProj(N,beta,SOD,N_d,delta_gamma)
shep=[1 .69 .92 0 0 0
    -.8 .6624 .8740 0 -.0184 0
    -.2 .1100 .3100 .22 0 -18
    -.2 .1600 .4100 -.22 0 18
    .1 .2100 .2500 0 .35 0
    .1 .0460 .0460 0 .1 0
    .1 .0460 .0460 0 -.1 0
    .1 .0460 .0230 -.08 -.605 0
    .1 .0230 .0230 0 -.606 0
    .1 .0230 .0460 .06 -.605 0];
gamma=delta_gamma*(-N_d/2+0.5:N_d/2-0.5);%扇束角度

rho=shep(:,1); %椭圆对应密度
ae=0.5*N*shep(:,2);%椭圆短半轴
be=0.5*N*shep(:,3);%椭圆长半轴
xe=0.5*N*shep(:,4);%椭圆中心x坐标
ye=0.5*N*shep(:,5);%椭圆中心y坐标
%投影数据生成%
beta=beta*pi/180;
alpha=shep(:,6);%椭圆旋转角度
alpha=alpha*pi/180;
gamma=gamma*pi/180;%角度换算成弧度
beta_num=length(beta);
P=zeros(N_d,beta_num); %存放投影数据
for k1=1:beta_num
    theta=beta(k1)+gamma;
    P_beta=zeros(1,N_d);%存放每一旋转角度下的投影数据
    for k2=1:length(xe)
        rsq=(ae(k2)*cos(theta-alpha(k2))).^2+(be(k2)*sin(theta-alpha(k2))).^2;
        dsq=(SOD*sin(gamma)-xe(k2)*cos(theta)-ye(k2)*sin(theta)).^2;
        temp=rsq-dsq;
        ind=temp>0;
        P_beta(ind)=P_beta(ind)+rho(k2)*(2*ae(k2)*be(k2)*sqrt(temp(ind)))./rsq(ind);
    end
    P(:,k1)=P_beta;
end

    
function rec_RL=FanBeamAngleResorting(P,N,SOD,delta_beta,delta_gamma)
[N_d,beta_num]=size(P);
delta_theta=0.5;
theta=0:delta_theta:359;
Np=257;
Mp=length(theta);
delta_gamma=delta_gamma*pi/180;
delta_theta=delta_theta*pi/180;
delta_beta=delta_beta*pi/180;
d=SOD*sin((N_d-1)/2*delta_gamma)/((Np-1)/2);
pp=zeros(N_d,Mp);PP=zeros(Np,Mp);
%第一内插%
m1=zeros(N_d,Mp);
for k1=1:N_d
    for k2=1:Mp
    t=k2*(delta_theta/delta_beta)-(k1-(N_d-1)/2-1)*(delta_gamma/delta_beta);
    n=floor(t);
    m1(k1,k2)=n;
    u=t-n;
    if n>=1 && n<beta_num
        pp(k1,k2)=(1-u)*P(k1,n)+u*P(k1,n+1);
    end
    end
end
%第二步内插%
for k1=1:Mp
    for k2=1:Np
       tt=1/delta_gamma*asin((k2-(Np-1)/2-1)*d/SOD)+(Np-1)/2+1;
       m=floor(tt);
       uu=tt-m;
       if m>=1 && m<N_d
           PP(k2,k1)=(1-uu)*pp(m,k1)+uu*pp(m+1,k1);
       end
    end
end
rec_RL=iradon(PP,theta,'linear','Ram-Lak',N);
end