CT二维平行光束重建

clc,clear;
%% 创建图像及投影
M=512; %模型图像的长宽 %建议先和接收传感器的个数一样，这里不一样
P = phantom(M); % 创建一个 shepp-logan 模型，原图像
theta=0:0.5:179.999; %角度
% theta=1:1:180; %角度
[R,xp] = radon(P,(theta*180)/pi);%利用radon变换获得不同方向上的投影
xp_offset = abs(min(xp) + 1);  %探测器的中心
% plot(xp,R(:,50));   
plot(xp,R);  

%% 傅立叶变换，滤波反投影重建算法，对投影数据做一维傅里叶变换，然后对傅里叶变换结果进行
%滤波，最后对滤波的结果进行反投影
size(R,1);%367
width = 2^nextpow2(size(R,1)); %傅立叶变换的宽度，nextpow2就是靠的最近的2的指数
proj_fft = fftshift(fft(R, width)); %快速傅里叶变换
filter =1 - 2*[0:(width/2-1), width/2:-1:1]'/width; %定义R-L是一种基础的滤波算法
% figure(1);
% plot(filter);

proj_filtered = zeros(width,360);
for i = 1:360
    %.*表示矩阵a中的元素与矩阵b中的元素按位置依次相乘
    %对每一列（每一个角度的投影数据的傅里叶变换值做滤波）
    proj_filtered(:,i) = proj_fft(:,i).*filter; 
end

%% 逆傅里叶变换并反投影
proj_ifft = real(ifft(fftshift(proj_filtered))); %取复数的实数部分
fbp2 = zeros(M); %创建一个和原图等大小的矩阵，且假设图像初值为0

for i = 1:360
    rad = theta(i); %遍历投影角度，%弧度， %这个rad 是投影角，不是投影线与x轴夹角，他们之间相差 pi/2
    for x = (-M/2+1):M/2
        for y = (-M/2+1):M/2
            %最近邻插值法
            t2 = round(x*cos(rad + pi/2)+y*sin(rad + pi/2)+xp_offset); %将每个元素舍X入到最接近的整数
            fbp2(x+M/2,y+M/2)=fbp2(x+M/2,y+M/2)+proj_ifft(t2 ,i);
        end
    end
end
fbp2 =(fbp2*pi)/180;%256x256 原图像每个像素位置的密度


%% 显示结果
figure(2);
subplot(1, 2, 1), imshow(P,[]), title('Original')
subplot(1, 2, 2), imshow(fbp2,[]), title('round(x*cos(rad + pi/2)+y*sin(rad + pi/2)+xp_offset)')