平行束投影数据的仿真（matlab）

1.Shepp-Logan头模型显示仿真

clc;
clear all;
close all;
I = phantom(256);

%phantom函数用于生成一个Shepp-Logan头模型;
%方法一：phantom(DEF,N); DEF为‘Shepp-Logan'和'Modified Shepp-Logan'(默认）；
%方法二：phantom(E,N); E为自定义头模型参数，可不填；N为生成图像大小，默认256；

figure;
imshow(I, []); 

%imshow函数用于显示图像（主要为灰度图像）；
%imshow(I); 显示灰度图像I；
%imshow(I, [low high]);  显示介于low和high之间的灰度图像（0-255），大于high显白，小于low显黑；  

2.平行束投影原理

3.仿真实验

方法一：利用radon变换产生投影数据 

RADON变换：是一种原始灰度图像到（ρ,θ）二维矩阵的映射，映射关系是灰度图像的像素值对参数（ρ,θ）的直线的线积分，或者叫投影，可理解为垂直于这个直线方向的像素值累计求和；

clc;
clear all;
close all;
I = phantom(256);
theta = 0: 179; %投影角度从0到179°;
P = radon(I, theta); 
%radon函数:返回灰度图像I在theta角度下的radon变换;

figure;
imshow(I, []);

figure;
imagesc(P), colormap(gray), colorbar;
%imagesc:自动缩放读入的图像；
%colormap：定义图像显示的颜色；
%colorbar：用于显示色彩条；

 

方法二：利用第二章内容计算投影数据

%edit medfuncParallelBeamForwardProjection；创建子程序
function P = medfuncParallelBeamForwardProjection(theta, N, N_d)
theta_num = length(theta);
P = zeros(N_d, theta_num);

%     椭圆密度rho  长轴ae   短轴be     xe     ye   旋转角度alpha
shep = [1         .69       .92       0      0        0
        -.8       .6624     .8740     0     -.0184    0
        -.2       .1100     .3100     .22    0        -18
        .1        .2100     .2500     0      .35      0
        .1        .0460     .0460     0      .1       0
        .1        .0460     .0230    -.08   -.605     0
        .1        .0230     .0230     0     -.606     0
        .1        .0230     .0460     .06   -.605     0];

rho = shep(: , 1).';
ae = 0.5 * N * shep(:, 2).';
be = 0.5 * N * shep(:, 3).';
xe = 0.5 * N * shep(:, 4).';
ye = 0.5 * N * shep(:, 5).';
alpha = shep(: , 6).';
alpha = alpha * pi/180;
theta = theta * pi/180;
TT = -(N_d-1)/2: (N_d-1)/2;      %探测器坐标-183到183，总计367个
for k1 = 1: theta_num            %theta_num(k1)代表t
    P_theta = zeros(1, N_d);     %一维向量[1, 367]
    for k2 = 1: max(size(xe))    %k2代表第几个椭圆
        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;          %保证根号>0,temp为根号内结果
        P_theta(ind) = P_theta(ind) + rho(k2) * (2*ae(k2)*be(k2)*sqrt(temp(ind)))./a;
    end
    P(: , k1) = P_theta.';
end

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

4.后续理解

4.1平行束投影公式理解

 4.2Sinogram举例

一个点在sinogram中展现应该为一条平行线：因为该点在中心，所以无论探测器从何角度投影，都会是中心位置的探测器探测得到同一投影值，所以是一样平行线，且灰度值一样。如下图：

而不在原点的一点在sinogram上展现应为正弦图，d=r*sinΘ，如下图：

 

一条线在sinogram中展现应该为一个亮点与很多条暗线： 亮点为投影角度刚好与该线重合，此时投影值最大，且仅有一个探测器探测到，反应为一个亮点；而在其他位置探测器探测到的投影值一样大，所以在sinogram中，在同一投影角度下，展现为一条竖直暗线。如下图：