DRR(Digitally Reconstructured Radiograph)分类及原理

DRR(Digitally Reconstructured Radiograph)，全称为数字重建放射影像，其是通过将三维 (3D) 图像（Volume）透视投影到二维 (2D) 图像平面上而生成的射线照相图像的模拟。DRR被广泛应用于CT模拟定位、图像引导放射治疗（IGRT）及计算机辅助外科等领域。

目前，DRR重建算法主要采用光线投射法(Ray casting)实现。光线投射是一种经典体绘制方法，在计算机图形学和三维可视化方面得到了广泛的应用，其采用光线投射的原理来模拟X线穿透人体并经过人体组织吸收后衰减而产生DRR图像的过程。

Digitally reconstructed radiographs (DRR) are a simulation of radiographic images produced through a perspective projection of the three-dimensional (3D) image (volume) onto a two-dimensional (2D) image plane.

光线投射生成DRR的方法

传统的直接插值法

直线插值方法的的主要过程是：从X线出束口发出多条X射线穿过三维体数据；在每条射线上进行等间距采样，并利用三线性插值方法由距离采样点最近的8个体素计算出该采样点的CT值对应的衰减系数；然后从前至后对所有采样点的衰减系数进行累加，得到该条射线对应成像平面的像素点的灰度值。针对每条射线重复上述过程后，将计算得到的所有像素点合成为一幅完整的DRR图像。基于光线投射方法的DRR生成如图所示：

SAD 为 X 射线源点到体数据等中心的距离，SFD 为 X 射线源点到 DRR 成像平面的距离。从 X 射线源点出发的射线条数与 DRR 图像的像素点个数相同，两者之间是一一对应的关系。为了简化采样计算过程从而提高采样效率本文选择等步长和三线性插值来计算采样点的 CT 值。该方法也是基于探测器的(Detector-based)的一种DRR生成方式。

换一种方法说，传统直接插值的光线投射算法，其具体实现过程可以描述为：首先利用计算机模拟一个虚拟点光源，用来表示传统 X 射线光源，然后从该点发射出若干条虚拟 X 射线，射线穿过三维数据集，照射到垂直与射线中线的投影面板上。投影面板上所有射线的投影点即是 DRR 图像的所有像素点。在投射过程中， X 射线按照既定的步长穿过三维数据集，当射线每经过一个 CT 切片，算法都会模拟真实 X 射线穿过人体时的衰减而得到一个 CT 值，如果射线与切片的交点恰巧不是像素点的时候就需要通过插值算法来估算出此时的 CT 值。最终将每条射线上所有交点的 CT 值累加，通过一定的转换方式将累加值转换成图像灰度，就获得了 DRR 图像。

总之，在传统直接插值的光线投射算法的实现过程中，需要获得每条射线与三位体数据的交点的 CT 值，若交点恰巧没有位于体素点上，则需要通过插值的方式来估算出该点的 CT 值，这个的过程需要大量的插值和求整运算，所以大大降低了 DRR 图像的生成速度，因此会有多种方式来加速该射线生成方式，例如通过三维 Bresenham 直线生成算法改进来快速的生成一条射线

Bresenham 直线生成算法

Bresenham 算法是一种高效的直线生成算法，传统的 Bresenham 算法应用在二维空间内，其原理如图所示：图中点 A， D， E 为实际像素点，点 B 为实际直线与 AD 连线的交点。该算法所要完成的任务就是判断将点 B 绘制在 D 点，还是绘制在 A 点。最为简洁的方法就判断 AB 与 BD 的距离，如果 AB>BD，则将 B 绘制在 D点，反之将其绘制在 A 点。另外一种辨别方法就是判断点 B 与点 C 的坐标大小，若B>C，则绘制在 A 点，若B<C ，则绘制在 D 点。

知乎上面也有个专栏很形象地介绍了该过程：https://zhuanlan.zhihu.com/p/106155534

使用Bresenham 算法，就能够快速的获得这条直线与三维数据集的交点序列，这个序列中的每一个像素点都是最接近真实直线与三维数据集的每一个交点，这样就省去了“传统直接插值的光线追踪算法”的插值和求整运算，只需要进行少量的浮点运算，就能够获得每一条射线上交点的总的 CT 值，大大加快了 DRR 图像的生成速度。

SiddonJacobs的光线追踪算法

该方法在itk remote module中的TwoProjectionRegistration有进行实现(链接：https://github.com/InsightSoftwareConsortium/ITKTwoProjectionRegistration)，该库也可以直接在编译安装itk的过程中直接build。

该方法最初由R. L. Siddon在文章"Fast calculation of the exact radiological path for a 3-dimensional CT array"中介绍，然后由F. Jacobs, E. Sundermann, B. D. Sutter, and I. Lemahieu在文章"A fast algorithm to calculate the exact radiological path through a pixel or voxel space"进行了改进。如果要看得很细节，可以参考后者的文章"A fast algorithm to calculate the exact radiological path through a pixel or voxel space"

光线追踪算法究其细节，就是要获得该条射线穿过了哪儿些体素，然后同时要获得这些体素的CT值以及光线穿过该体素时候经过的长度。

SiddonJacobs算法的优势之处就在于，使用α来表示该条射线归一化的长度，其值为0则表示该点是在光源点，其值为1则表示该点是在Detector上的像素点，使用该α可以同时追踪光线和在CT上的像素位置。

下面来详细说明该图，每个方块代表一个体素（图上为了方便只写2D情况，3D下可以很简单地类比过去），每一条横线代表一个x平面，每一条纵线代表一个y平面，p1为光源，p2为Detector上的像素点。在最开始，只知道CT上每个x y平面的物理空间位置以及p1和p2点的物理空间位置，并不清楚入射点出射点的坐标，甚至不清楚光线是否与CT相交。

首先，通过文章"A fast algorithm to calculate the exact radiological path through a pixel or voxel space"里的公式，将所有的α求出来，α由两个数组构成，分别是“光线与全部x平面的交点” 以及 “光线与全部y平面的交点”，利用该公式，α又同时能反过来算出是与哪儿个x平面或者哪儿个y平面相交（这点很重要，因为最后求和的时候，需要累加光线与体素相交的长度及其对应的CT值，CT值就需要知道是哪儿个像素才行）。
其次，将全部的α升序排列，然后邻位相减，这样就可以得到每小段长度的比例。
最后，将上述邻位相减的α值，依次乘以光线总长度，依次得到每小段光线的真实长度；再通过第一步骤计算α的公式，依次反过来拿回每小段光线经过像素的坐标，从而得到每小段光线经过像素的CT值，将两个数组对应相乘，再求和，即可得到该条光线的模拟x光射线的数值。