基于MATLAB的CT平行束投影数据的仿真

实验二 平行束投影数据的仿真

实验任务：用两种方法实现对SL头模型平行束投影数据的仿真计算，头模型图像大小为256*256，投影角度为0，1，。。。，179。

方法一：利用Randon变换产生投影数据
I=phantom(256);
theta=0:179;
P=radon(I,theta);
figure;
imshow(I,[]),title('256*256头模型图像');
figure;
imagesc(P),colormap(gray),colorbar,title('180°平行束投影图像');

推导过程参考教材：医学断层图像重建仿真实验

方法二：利用解析法产生投影数据
%%--=-=仿真参数设置-=---%%
N = 256;% 重建图像大小，探测器通道个数
theta =0:1:179;% 投影角度
%%=====产生仿真数据- ====%%
I= phantom(N);% 产生 shepp_logan 头模型
N_d = 2 * ceil(norm(size(I)-floor((size(I)-1)/2)-1))+1;% 设定探测器通道个数
P = medfuncParallelBeamForwardProjection(theta,N,N_d);% 产生投影数据
%%-====仿真结果显示-====-%%
figure;% 显示原始图像
imshow(I,[]),title('256×256 头模型图像');
figure; imagesc(P),colormap(gray),colorbar,title('180°平行束投影图像');

子程序：
function P = medfuncParallelBeamForwardProjection(theta,N,N_d)
    shep=[   1   .69     .92     0       0       0
        -.8 .6624   08740   0       -.0184  0
        -.2 .1100   .3100   .22     0       -10
        -.2 .1600   .4100   -.22    0       15
        .1  .2100   .2500   0       .35     5
        .1  .0460   .0560   .06       .1      15
        .1  .0260   .0460   .13       -.1     20
        .1  .0460   .0130   -.08    -.605   -5
        .1  .0430   .0230   .03       -.606   3
        .1  .0230   .0460   .06     -.605   0];
   theta_num = length(theta);
    P = zeros(N_d,theta_num); % 存放投影数据
    rho = shep(:,1).';  % 圆对应密度
    ae =0.5*N*shep(:,2).'; %椭圆短半轴
    be = 0.5*N*shep(:,3).'; % 椭圆长半轴
    xe =0.5*N*shep(:,4).'; %圆中心x坐标
    ye = 0.5*N*shep(:,5).'; %圆中心y坐标
    alpha = shep(:,6).'; % 圆旋转角度
    alpha = alpha * pi/180;
    theta = theta * pi/180; % 角度换算成弧度
    TT =-(N_d-1)/2:(N_d-1)/2;% 探测器坐标
    for k1 = 1 : theta_num
        P_theta = zeros(1,N_d);% 存放每一角度投影数据
        for k2 = 1: max(size(xe))
            a =(ae(k2) * cos(theta(k1)-alpha(k2)))^2+(be(k2) * sin(theta(k1)-alpha(k2)))^2;
            temp = a- (TT-xe(k2) * cos(theta(k1))-ye(k2) * sin(theta(k1))).^2;
            ind = temp>0;% 根号内值需为非负
            P_theta(ind)= P_theta(ind) +rho(k2) * (2 * ae(k2) * be(k2) * sqrt(temp(ind)))./a;
        end
        P(:,k1) = P_theta.';
    end
end