滤波反投影fbp算法matlab

用matlab实现ct滤波反投影fbp算法，其中投影可以直接用radon变化实现，但代码中手动编写投影部分的代码，使用的滤波器为R-L滤波器

%%投影
I = phantom（‘Modified Shepp-Logan‘256）;
nmax=ceil（length（I）*（sqrt（2）））;
I1=zeros（512）;
I1（128:383128:383）=I;
proj=zeros（512180）;
for theta=1:180
    I2=imrotate（I1-theta‘bilinear‘）;%原图旋转theta度，得到I2
    temp=length（I2）;
    temp=round（temp/2-255）;
    I3=I2（temp:temp+511temp:temp+511）;
    b=sum（I3）;%b是投影向量
    proj（:theta）=b‘;%赋值投影矩阵
end
imshow（proj[]）
%%投影

% width = 2^nextpow2（size（I1））;
width=512;
H = 2*[0:（width/2-1） width/2:-1:1]‘/width;%R-L滤波
proj_filtered = zeros（length（proj）180）; 
% plot（H）
proj_fft = fft（proj 512）;
for i = 1:180    
    proj_filtered（:i） = proj_fft（:i）.*H;  
end  

% M=256;
% fbp = zeros（M）; % 假设初值为0
proj_ifft = real（ifft（proj_filtered））; 
% imshow（proj_ifft[]）
fbp = zeros（256）; 

% 5. back-projection to the x- and y- axis  

for i = 1:180  
    % rad is the angle of the projection line  not projection angle  
    rad = i*pi/180;  
    for x = （-256/2+1）:256/2  
        for y = （-256/2+1）:256/2  
            t = round（x*cos（rad+pi/2）+y*sin（rad+pi/2））;  
            fbp（x+256/2y+256/2）=fbp（x+256/2y+256/2）+proj_ifft（t+256i）/180;  
            
        end  
    end  
%     imshow（fbp[]）
end  
fbp = fbp/360;  

subplot（1 2 1） imshow（I） title（‘Original‘）
subplot（1 2 2） imshow（fbp[]） title（‘FBP‘）

% I=dicomread（‘C:\Users\56326\Desktop\dcm\100.dcm‘）;
% imshow（I[]）
% info= dicominfo（‘C:\Users\56326\Desktop\dcm\100.dcm‘）;
% I1=imread（‘C:\Users\56326\Desktop\sonoscape\12345\20180425_145827\20180425_150454_1.jpg‘）;
% I2=imread（‘C:\Users\56326\Desktop\sonoscape\12345\20180425_145827\20180425_150545_2.jpg‘）;
% I3=imread（‘C:\Users\56326\Desktop\sonoscape\12345\20180425_145827\20180425_150603_3.jpg‘）;
% 
% figure
%%
theta = 1:180;  
  P=phantom（128）;
% 1. projection  using radon function  
[Rxp] = radon（Ptheta）;  
width = 2^nextpow2（size（R1））;  % set width for fft transformation  
  
% 2. do fft to the projection  
proj_fft = fft（R width）;  
  
% 3. filter  
% Ramp filter function  from 0 to width then to  0  
filter = 2*[0:（width/2-1） width/2:-1:1]‘/width;  
proj_filtered = zeros（width180）;  
for i = 1:180  
    proj_filtered（:i） = proj_fft（:i）.*filter;  
end  
  
% 4. do ifft to the filtered projection  
proj_ifft = real（ifft（proj_filtered））; % get the real part of the result  
  
% 5. back-projection to the x- and y- axis  
fbp = zeros（128）; % asign the original value 0  
for i = 1:180  
    % rad is the angle of the projection line  not projection angle  
    rad = i*pi/180;  
    for x = （-128/2+1）:128/2  
        for y = （-128/2+1）:128/2  
            t = round（x*cos（rad+pi/2）+y*sin（rad+pi/2））;  
            fbp（x+128/2y+128/2）=fbp（x+128/2y+128/2）+proj_ifft（t+round（size（R1）/2）i）;  
        end  
    end  
    imshow（fbp[]）
end  
fbp = fbp/180;  
  
% 6. show the result  
subplot（1 2 1） imshow（P） title（‘Original‘）  
subplot（1 2 2） imshow（fbp） title（‘FBP‘）  
%%

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

     文件       3022  2018-05-16 14:41  fbp-.m

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

                 3022                    1