利用射线驱动方法中的等距采样法获取二维投影的matlab代码

之前一直看不懂利用射线驱动方法中的等距采样法获取二维投影的matlab代码，今天终于看懂了，现在记录一下。

由于苦恼很多网上现成的代码运行不了，或者却文件，我这里声明一下，我的博客，以后所有的博客，我能力范围内的，我都是保证绝对能运行的(不会缺少包),我真的是烦死了，所以我自己发布的博客会尽可能避免这个问题。

这是我用的自己搞得体膜：

这是结果：

代码如下（有详细的注释）：

%%首先创建一个二维矩阵，为了能够便于之后更清楚的看到投影图像，这里采用了偏心圆
% 创建一个 256x256 的矩阵
[x, y] = meshgrid(1:256, 1:256);

% 定义新的圆心和半径
new_center = [180, 100];
radius = 100;

% 计算每个点到新圆心的距离
distances = sqrt((x - new_center(1)).^2 + (y - new_center(2)).^2);

% 创建一个 256x256 的矩阵，圆内像素值为 1，圆外像素值为 0
img = zeros(256);
img(distances <= radius) = 1;

% 显示图像
imshow(img);

%%%======================== paramsetting ======================%
%% 参数设置(注：本段程序引入了"param"，初步的面向对象思想)
param.Geom = '2D'; % 射束几何，包括"2D","3D"
% 被测二维图像的参数设置
param.ny = 256;  param.nx = 256;    % 待重建二维图像像素数[nx, nx]
param.sy = 500;  param.sx = 500; % 被测图像尺[sx, sx]，mm
param.dy = param.sy/ param.ny; param.dx = param.sx / param.nx; % 待重建图像像素尺寸
param.dh = 0.1* param.dx; %等距采样距离，即X向(射线方向)划分出2*nx-1=1023个
% 待重建二维图像X向尺寸坐标
img = phantom('Modified Shepp-Logan', param.nx);
% 直线探测器的参数设置
param.nu = 256;     % 直线探测器探元数量
param.su = 1000;    % 直线探测器尺寸，mm
param.du = param.su /param.nu;  % 直线探测器探元尺寸，mm
% 直线探测器存在的位置偏差
param.off_u = 0;    % 中心位置偏移 mm
param.rotate = 0;    % 角度偏差(度)
% 射线源和探测器距离设置
param.SOD = 750;   % 射线源到旋转中心的距，mm
param.SDD = 1500;  % 射线源到直线探测器的距离，mm
% 被测图像、探测器边界尺寸
param.img_origin = (param.sx- param.dx)/2;
param.dete_origin = (param.su- param.du)/2;
% 直线探测器坐标序列(负->正，原点在中心)
param.us = (-param.dete_origin: param.du : param.dete_origin)'; 
param.xhrzt = (-param.img_origin:param.dh:param.img_origin)';
% 二维被测图像的实际坐标矩阵(负->正，原点在中心)
[param.xs, param.ys] = meshgrid(-param.img_origin: param.dx : param.img_origin);
% 建立存放在GPU上的直线探测器实际尺寸坐标，中心为原点 
%param.us=gpuArray.colon(-param.dete_origin, param.du, param.dete_origin)'+ param.off_u;
% 被测图像扫描角度设置
param.dir = 1;      % 旋转方向(正，被测图逆时针旋转，等效源-探测器顺时针旋转)
param.dbeta = 1;    % 旋转角度步长(度)
%========================%
param.parker = 0;   % 若短扫描，使能parker加权(暂先忽略，详见下文4.5小节)
param.parker_accelerate = 1; % 采用加速版(矩阵真值索引方式)parker加权
if param.parker == 1  
    param.betas= 0:param.dbeta:180+2*atand(param.su/2/param.SDD);%短扫描 0 :π+2γm
else
    param.betas = 0:param.dbeta:360-1;   % 360度全扫描    
end
%========================%
param.betas = param.betas *param.dir;  % 旋转角度序列(度)
param.nProj = length(param.betas);     % 旋转扫描的角度个数
% 被测体旋转角度设置
param.dir   = 1;   % 被测体旋转方向，1：逆时针，-1：顺时针（注：机架旋转方向与被测体相反）
param.dang  = 1;    % 角度增加的步长（度）
% 需要进行全扫描还是扫描角度序列，自定义扫描角度范围
param.deg   =  (0:param.dang:359)-90;    
param.deg   =  param.deg*param.dir; % 旋转角度(度)序列
param.nProj =  length(param.deg);   % 角度序列长度

% 若为全扫描（360度）     ->参数parker=0 (还未知该功能作用)
% 若为短扫描（小于360度） ->参数parker=1
param.parker = 0;   

%%
%获取投影图像
[xx, yy]= meshgrid(param.xhrzt, param.us); %细化x网格，建立旋转等距采样网格，param.xhrzt是图像坐标序列，param.us是直线探测器坐标序列
proj = zeros(param.nu, param.nProj);  %预设空间存放投影数据
for i = 1:1: param.nProj   
    %已知图像中心为旋转轴，旋转的实质就是将图中的各个向量进行旋转，然
    %后将旋转前的位置的像素值赋值给旋转后的位置的像素值（图像旋转之后，
    %会出现许多空洞点。对这些空洞点必须进行填充处理，即内插/插值）。
    %旋转图像（与机架旋转方向相反）
    rx = xx.*cosd(param.deg(i)) - yy.*sind(param.deg(i)); 
    ry = xx.*sind(param.deg(i)) + yy.*cosd(param.deg(i)); % cosd计算的是角度
    img2 = interp2( param.xs, param.ys, img, rx,ry,'linear',0 ); %采用双线性内插（详见下文5.4.2小节“FDK算法的计算机实现”尾部）进行对应像素值的赋值。
    %将每行（x为射线积分方向）灰度值累加，乘采样距离dh，获得投影数据
    proj(:, i)=param.dh.*sum(img2,2);  
end

%%
%展示由等距采样法获取到的投影图像的结果
imshow(proj);