扇形束CT重建算法

扇形束CT成像几何中，滤波反投影算法，Matlab版本。

clear all
close all
load head.mat;
% head =phantom（256）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%获得扇形束等距离投影数据
D = 512;                                                                                %%%%射线源到图象中心的距离
[shade] = fanbeam（headD‘FanSensorGeometry‘‘line‘‘FanSensorSpacing‘1）;    %%%%获得投影数据 线形检测器 通道数sp  旋转角度360。
figure;
imshow（shade[]）;title（‘360度投影数据‘）;
clear Xangle;
%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%修正投影数据
[XsizeYsize] = size（shade）;
for n = 1:Ysize
    for m = 1:Xsize
        shade_temp（mn） = shade（mn）*D/sqrt（D.^2+（m-185）.^2）;
% shade_temp（mn） = shade（mn）;
    end
end;
% shade_temp = （D/sqrt（D.^2+（m-185）.^2））*shade;
clear Xsize;
clear Ysize;
figure
imshow（shade_temp[]）;
title（‘修正后的投影数据‘）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%进行R_L卷积滤波

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%生成R_L卷积核
t=linspace（-128127256）;
filt=0.0085*（sinc（t）/2-sinc（t/2）.^2/4）;
figure;
plot（filt）;
title（‘R__L 滤波函数 256点‘）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%卷积滤波
shade = shade_temp;
clear shade_temp;
shade_temp = zeros（881360）;
shade_temp（257:625:） = shade;
for n = 1:360
    temp（:n） = conv（shade_temp（:n）filt）;
end
clear shade_temp;
% % % % % % % % % % shade_temp =abs（temp（385:753:））;
shade_temp =temp（385:753:）;
size（shade_temp）;
figure
imshow（shade_temp[]）;
title（‘滤波后投影数据‘）;
clear filt;
clear shade;
clear temp;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%进行反投影
re = zeros（256256）;
d_th = pi/180;
for view = 1:360
    for x = 1:256
        for y = 1:256
            u = （D-（x-128.5）*cos（（view-1）*d_th）-（y-128.5）*sin（（view-1）*d_th））;
            n = D*（-（x-128.5）*sin（（view-1）*d_th）+（y-128.5）*cos（（view-1）*d_th））/u+185;
            if （1                nn = fix（n）; n0 = n-nn;
%                 re（xy） = （（1-n0）*shade_temp（nnview）+n0*shade_temp（nn+1view））*D^2/u^2;
                re（xy） = re（xy） + （（1-n0）*shade_temp（nnview）+n0*shade_temp（nn+1view））*D^2/u^2;

            else re（xy） = re（xy）;
            end
        end
    end
%     view
%     figure
%     imshow（re[]）
end
figure
imshow（re[]）
title（‘重建后的图象‘）;
toc

figure 
imshow（head[]）
title（‘原始图象‘）;
p1 = re（128:）;
p2 = head（128:）;
figure
plot（p1

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

     文件       3360  2009-06-07 20:58  fb_t.m

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

                 3360                    1