三维锥束CT反投影滤波重建matlab程序

FDK进行锥束CT重建程序，包括滤波和反投影两个程序

%compiled by: Fangzheng
%date: 2007/12/28
clear all
Nfile=200;%360度的图像采集数目
ANG_pace=2*pi/Nfile;%步进角
L_middle=761;%中平面位置
D=963;
rr=min（1920-DD-1）;
H=2*rr+1;

R=3248;%光源到探测器的距离
gama_max=atan（rr/R）;%最大扇角
N_shadows=Nfile%ceil（Nfile/2+2*gama_max/ANG_pace）-1%半扫描需要的张数
Begin_shadows=1;%Begin_shadows=5;%起始扫描的步进数
End_shadows=Begin_shadows+N_shadows-1;%终止扫描的步进数

beta=0: ANG_pace: ANG_pace*End_shadows;%beta角序列

for L_reconstruct=761:761%重建起止切片只能从48到977
    Z=L_reconstruct-L_middle%重建切片Z坐标
    
    L_left  = floor（min（（R*Z/（R+sqrt（2）*rr））+L_middle（R*Z/（R-sqrt（2）*rr））+L_middle））%需要读取Z上最小位置
    L_right =  ceil（max（（R*Z/（R+sqrt（2）*rr））+L_middle（R*Z/（R-sqrt（2）*rr））+L_middle））%需要读取Z上最大位置
    N_l=L_right-L_left+1;%需要Z序列数
    FDK2=zeros（N_lHN_shadows）;

    for i=1:N_shadows
        ii=i+Begin_shadows-1;
        fndat = sprintf（‘I:/10.4-mid/%d.mid‘ ii）;
        fd = fopen（fndat ‘r‘）;
        x=fread（fd ‘float‘）;%读原始中间文件
        fclose（fd）;
        x=reshape（x[15361920]）;
        FDK2（::i）=x（L_left:L_rightD-rr:D+rr）;%生成必要数据三维阵
    end 

    fbpf=zeros（H）;%生成H*H的零方阵  
% fbpf=zeros（2048）;
    for i=1:H
        x=i-rr;
        for j=1:H
            y=j-rr;
              if sqrt（x*x+y*y）<（rr-50）
                 for k=Begin_shadows:End_shadows
                    u=R+x*cos（beta（k））+y*sin（beta（k））;
                    a=R*（-x*sin（beta（k））+y*cos（beta（k）））/u;
                    aa=round（a+rr）;
                    b=Z*R/u;
                    bb=round（b+L_middle-L_left+1）;
                    kk=k-Begin_shadows+1;
                    fbpf（ij）=fbpf（ij）+FDK2（bbaakk）*R/（u*u）;
                 end
              end
        end
    end%逆投影
    fndat = sprintf（‘I:/10.4-rst/%d.rst‘L_reconstruct）
    fd = fopen（fndat ‘w‘）;
    fwrite（fd fbpf ‘float‘）;%写切片文件
    fclose（fd）;
    clear fbpf;
    clear FDK2;
end
figureimagesc（fbpf）title（int2str（L_reconstruct））

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       2062  2013-07-13 20:51  xb_fdkback.m

     文件       1418  2009-10-06 09:38  xb_fdkfilter.m

----------- ---------  ---------- -----  ----

                 3480                    2