matlab提供大量函数，可以方便的完成fbp算法

1)fbp算法原理：

中心切片定理 (CST) :  原数据投影的一维傅立叶变换等于原数据的二维傅立叶变换

投影  -->  一维傅立叶变换 --> 滤波 --> 二维傅立叶反变换

经过上述过程应该得到原始数据

2)投影相关知识

2.1)正投影：对投影线经过的像素做线积分，积分得到的值保存为该角度下的权值

对一组数据 P 做 Radon 变换，即做正投影，会得到两个数据

[R, xp] = radon(P,theta);

xp是投影线条数

R是theta角下第 xp 条投影线得到的线积分，即权值

2.2) 反投影：反投影是利用上面投影得到的 R 权值，把R值投回到 x y 坐标中

x y 满足 x*cos(theta) + y*sin(theta) = R 就表明点(x, y)在投影线上

注：这个theta 是投影角，不是投影线与x轴夹角，他们之间相差 pi/2

此时把R值全部赋给投影线上的每一个点，最后每个角度的R值叠加到这个点上，在除以总投影线数，就得到了原始数据

分析：有的像素点的原始数据是0，但是这样一反投影就会有数据

产生星状伪迹，所以要滤波

3) matlab 源代码

P = phantom(128); % create a shepp-logan

theta = 1:180;

% 1. projection using radon function

[R,xp] = radon(P,theta);

width = 2^nextpow2(size(R,1)); % 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(width,180);

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 = theta(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/2,y+128/2)=fbp(x+128/2,y+128/2)+proj_ifft(t+round(size(R,1)/2),i);

end

end

end

fbp = fbp/180;

% 6. show the result

subplot(1, 2, 1), imshow(P), title('Original')

subplot(1, 2, 2), imshow(fbp), title('FBP')