CT图像重建算法------射线驱动投影模型

在模拟CT扫描过程中被测物体保持静止不动，射线源与探测器绕Z轴做逆时针旋转，每隔一定角度进行一次投影数据的计算，而计算方法则取决于采用什么样的投影算法。本文主要介绍了投影算法的分类，并详细描述了射线驱动算法中的joseph算法原理。

一、图像重建算法分类

CT图像重建算法主要有3类：1、反投影法；2、迭代重建算法；3、解析法:包括滤波反投影法和傅里叶变换法

二、迭代重建算法分类

迭代重建算法在求解的过程中需要不断地求解矩阵元素，目前系统矩阵的建模方式主要分为４种，分别为
（1）像素驱动模型（2）射线驱动模型（3）距离驱动模型（4）面积积分模型。

（1）像素驱动模型（PDM）通常应用在FBP重建反投影的计算；

（2）射线驱动模型（RDM）多用于迭代算法中的正向投影算法，这种模型下，射线不再看作是
有宽度的，它的原理是通过对射线经过的像素进行加权求和来模拟线积分；
包括：Joseph算法（使用线性插值）、Siddon算法（使用临近插值，也是射线驱动模型中最经典算法）
其中Joseph算法使用线性插值模型，其重建质量比Siddon算法要理想，但由于计算量大导致重建速度较慢。

（3）距离驱动模型（DDM）相较于前两种模型较为先进。这种模型下考虑到了像素与探测器单
元的宽度。另还有改进的距离驱动模型（IDDM）；

（4）面积积分模型（AIM）将射线穿过像素时与像素单元边界形成的图形面积作为权
因子。
本文主要介绍迭代重建算法中射线驱动模型中

三、Joseph算法原理

Joseph’s method即约瑟夫方法，是射线驱动模型下的一种前向投影算法，通过对射线经过的像素进行加权求和来模拟线积分。

1、二维Joseph算法

为了方便理解我们可以先从二维角度观察算法原理，如下图所示为joseph算法的二维原理图。

图1 Joseph算法基本原理图
算法主要根据焦点（也就是光源）与探测器每一个像素的中心两点所建立的光线斜率判断插值方式。Joseph利用线性插值，对于三维情况就是双线性插值。线性插值方式分为在X轴上插值和在Y轴插值两种。判断方法为：
若射线更平行于Y轴，则在Y轴上求和，在X轴上插值。也就是在X轴上投影，最后在Y轴方向求和。
反之，若射线更平行于X轴，则在X轴上求和，在Y轴上插值。也就是在Y轴上投影，最后在X轴方向求和。如图2

图2 Joseph算法在Y轴上投影的情况

数学上判断方式有很多，比如：方法1，计算焦点和当前探测器像素中心点两点的斜率；方法2，计算两点X和Y坐标的插值，比较他们的大小关系，若|xs-xd|<=|ys-yd|，则在X轴插值，若|xs-xd|>|ys-yd|，则在Y轴插值。图1就是在X轴插值的情况。根据插值方式，也就是插值坐标，就开始进行插值计算了。
二维情况下，重建图像被看作一个给定大小的网格像素图像（m*n）。像素大小为p，当射线穿过重建区域时，穿过图像的射线附近像素对其有影响，只考虑这些像素的衰减作用，其他像素对该射线的衰减无贡献。这个贡献称为权重因子，也就是线性插值过程所得的结果。具体过程为（以X轴插值为例）：
（1）判断相交，寻找交点
首先确定采样角以及所要使用的探测器像素确定射线L，而后判断L是否与重建区域相交，若相交，相交在哪一层像素的中心点所在平面上，并求出交点坐标。
判断相交的方式为遍历重建区域每一层像素中心点平面，也就是每一层中心点的纵坐标，判断射线是否与当前遍历的层相交，若相交则求出交点坐标，若不想交则继续遍历下一层。
（2）计算权因子
对于有交点的那一层，例如图1中最后一层，计算交点与相邻两个像素中点的距离占总距离的比值，分别为wl/deta和(deta-wl)/deta，即为公式3中的两个权因子：

其中f(Xk(l)，Yl)×较小的距离，而右侧的像素衰减值却乘以较大值，这是因为交点越接近一个像素点时，认为该像素对该路径的贡献越大，所以是反向相乘的。
式中前面的为比例因子，是一个雅可比矩阵，其中的deta表示相邻两个像素中心之间的距离，cos值为采样角的余弦值
（3）
当某一采样角度下，当前探测器像素对应的射线在某一层的插值结算结束后，即寻找下一个相交的层继续计算；
当所有层遍历结束后，转换到下一个探测器像素，继续以上步骤；
当所有探测器像素均计算完毕后，就继续下一个采样角度的所有计算。
至此，一个简单的二位Joseph算法过程执行完毕。

2、三维Joseph算法

三维算法基本的原理与二维是保持一致，只是需要加上一个Z方向产生的分量。
与二维不同的需要修改的部分：
（1）以在X轴上投影的一个情况作为例子，如图3所示。

图2 三维Joseph算法原理
对于三维情况，一条射线产生的权重因子需要由两个角度的因子共同决定，例如图2所示的情况下，分别需要XOY平面和YOZ平面求得的插值，如图3所示，而后将两个平面得到的插值进行乘积得到最终的三维插值结果

图4 X轴插值情况的两个视角
（2）对于最后的投影值的公式
需要将最前面的比例因子修改为deta/(cos(xoy）*cos（yoz）)，三维公式大致为：

注意：
1、由于我们模拟的是锥束投影，因此所说的投影角是每一个采样角度下焦点与探测器中心连线的旋转角度，而当前采样角度下整个探测器的其他像素块所对应的射线角度并不是投影角，所以最后求得的权因子需要再除以探测器每一个像素对应射线的角度的余弦值。
2、对于三维Joseph算法，在X轴投影时，需要XOY和YOZ两个平面；而在Y轴上投影时，需要YOX和ZOX两个平面的插值结果。
3、边界问题处理一定分清楚情况，确定好极限位置
在判断射线是否与重建区域相交时，需要考虑边界的极端情况，比如图1中第一行，交点位置已经超出重建区域，此类极端情况根据项目需求可以自行设定，一般的处理方式是，若交点超出重建区域的边界，但超出区域在半个像素内，仍然保留交点，继续计算。

四、Siddon算法

Siddon算法整体计算结构与Joseph算法基本一致，仅仅是在计算权重因子时有略微的区别。
在Siddon算法中，每个像素对于投影的贡献由像素内射线长度决定，例如下图射线与重建区域的相交情况：

此时，投影值P可以表示为：

五、坐标系建立