直接反投影重建算法（matlab）

1.原理推导

将仿真得到的投影数据，或者扫描得到的真实投影数据，沿着扫描路径按照原路回抹到路径经过的像素上，然后累加；数学公式的推导可参考YouTube上国立台湾大学成像课程https://www.youtube.com/watch?v=0255KLvu75g&t=1295s

2.实验

方法一："iradon"函数

>> clear all;
>> close all;
>> I=phantom(256);
>> theta=0:179;
>> P=radon(I,theta);
>> rec=iradon(P, theta, 'none');
>> figure;
>> imshow(I, []), title("raw image");
>> figure;
>> imshow(rec, []), title("imaging");

重建效果 

方法二：解析法

>> clear all;
>> close all;
>> 
>> N=256;
>> I=phantom(N);
>> delta=pi/180;
>> theta=0:1:179;
>> theta_num=length(theta);
>> 
>> P=radon(I, theta);
>> [mm, nn]=size(P);
>> e=floor((mm-N-1)/2+1)+1;
>> P=P(e:N+e-1,:);
>> P1=reshape(P,N,theta_num);

>> figure;
>> imshow(I, []);
>> rec=medfuncBackprojection(theta_num, N, P1, delta);
>> figure;
>> imshow(rec, []);


>> edit medfuncBackprojection;
function rec= medfuncBackprojection(theta_num, N, R1, delta)

rec=zeros(N);
for m=1:theta_num
    pm=R1(:,m);
    Cm=(N/2)*(1-cos((m-1)*delta)-sin((m-1)*delta));
    for k1=1:N
        for k2=1:N
            Xrm=Cm+(k2-1)*cos((m-1)*delta)+(k1-1)*sin((m-1)*delta);
            n=floor(Xrm);
            t=Xrm-floor(Xrm);
            n=max(1,n);n=min(n,N-1);
            p=(1-t)*pm(n)+t*pm(n+1);
            rec(N+1-k1,k2)=rec(N+1-k1,k2)+p;
        end
    end
end

以上实验内容均学习整理于黄力宇老师等著《医学断层图像重建仿真实验》。

3.后续理解

1.点源的直接反投影重建

 2.头模的直接反投影重建