滤波反投影matlab仿真程序

本程序用matlab实现滤波反投影，有助于我们理解滤波反投影算法的具体原理

clear;

disp（‘Display Original‘）
M=256;
N=256;
I = phantom（MN）;
figure（1）;
imshow（I）;

coordinateSource=[];
projMatrix=[];
detector=[];
proj=load（char（‘projection.mat‘））;
phyRatoDig=proj.phyRatoDig;
projMatrix=proj.projection;
yDetector=proj.yDetector;
nDetectors=proj.nDetectors;

figure（2）
showimge（projMatrix3605120max（max（projMatrix）））;

% coordinateSource=proj.coordinateSource;
% coordinateDetector=proj.coordinateDetector;
% 
D_dig=proj.focalDistance_dig;
sourceToDetector_dig=proj.focalDistance_dig+proj.detecDistance_dig;
s=[];
s=D_dig/sourceToDetector_dig*yDetector（1:）*phyRatoDig;
Detector=yDetector（1:）*phyRatoDig;
% 
 pe=[];
 M=D_dig./sqrt（D_dig.^2+s.^2）;
 nViews=proj.nViews;
% 
 for i=1:nViews
    pe（i:）=projMatrix（i:）.*M;%加权投影
     
 end
 figure（3）;
 showimge（pe3605120max（max（pe）））;


disp（‘Filtering‘）
filternum=128;
filter_ramp=zeros（filternum1）;
for j=1:filternum   % 16 point ramp filter
    i=j-1-filternum/2;
    if（i==0）
     filter_ramp（j1）=1/（8.0）;
    elseif （mod（i2）==0）
            filter_ramp（j1）=0;
    elseif （mod（i2）==1）
            filter_ramp（j1）=-0.5/（i*i*pi*pi）;
    end
end
m=1;
figure（4）;

plot（filter_ramp）;


pfilter=[];
length_conv=filternum+nDetectors-1;%639
pPro=zeros（length_conv1）;
temp_pro=zeros（nDetectors1）;
h_filter=filternum/2;%64
ii=length_conv-h_filter-nDetectors;%63
 for s=1:nViews % sample-loop
    % for  pp=1:h_filter
    pro_left =（pe（s1）+pe（s2））/2.0;
    pro_right=（pe（snDetectors）+pe（snDetectors-1））/2.0;
    
    for pp=1:h_filter;               %left part
        pPro（pp1）=pro_left;  %左边补零      
    end
%    
    for pp=1:nDetectors                      %middle part
     pPro（h_filter+pp1）=pe（spp）;%中间部分不变
    end
%    
   for pp=h_filter+nDetectors+1:length_conv
    pPro（pp）=pro_right;%右边补零
    end
%   result_conv    
   for n=1:nDetectors 
       result_conv=0;
   	for jj=1:filternum%1：128
        pPmove=pPro（n+jj-11）;
        result_conv=result_conv+pPmove*filter_ramp（jj1）;
    end
    pfilter（sn）=result_conv;
   end
 
    
 end
figure（5）;
showimge（pfilter360512min（min（pfilter））max（max（pfilter）））;%输出卷积滤波结果

% %%%%%%%%%%%%%%%%%%%%%%reArrange%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% %%%%%%%%%%%%%%%%%%%%%%Back projection%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  deltaBeta=proj.deltaTheta;
%  xCoodinate=proj.xCoodinate;
%  yCoodinate=proj.yCoodinate;
%  
%  unitDis= proj.unitDis;
% Detector =coordinateDetector（2:）;
disp（‘BackProjection‘）
detecLength=proj.detecLength;
unitDis=detecLength/（nDetectors-1）;
unitDis_dig=unitDis*phyRatoDig;

deltaBeta=2*pi/nViews;
 M=proj.M;
 N=proj.N;
 fReconstruct=[];
 for i=1:M
      x=i-（M+1）/2;%待重建点的横坐标
     for j=1:N
        y=（N+1）/2-j;%待重建点的纵坐标
       r=sqrt（x^2+y^2）;
        theta=atan2（yx）;%极坐标（rtheta）
%     
        result=0; 
     for s=1:nViews
     beta= （s-1）*deltaBeta;%中心线与y轴夹角
     s1=D_dig*r*cos（beta-theta）/（D_dig+

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

     文件       3946  2011-08-12 14:41  FBPRec\FBPRec\FanRec.asv

     文件       3946  2011-08-12 14:42  FBPRec\FBPRec\FanRec.m

     文件       6367  2009-12-16 18:21  FBPRec\FBPRec\limited_Views.m

     文件    1120183  2009-12-21 18:51  FBPRec\FBPRec\projection.mat

     文件       3622  2011-08-12 12:41  FBPRec\FBPRec\projecton.asv

     文件       3622  2011-08-12 12:42  FBPRec\FBPRec\projecton.m

     文件        365  2009-11-22 21:14  FBPRec\FBPRec\showimge.m

     目录          0  2011-08-12 13:41  FBPRec\FBPRec

     目录          0  2011-08-10 14:01  FBPRec

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

              1142051                    9