0 引言

传统的放射治疗采用解剖影像引导，即观察肿瘤的形态进行追踪和治疗，虽然具有良好的空间分辨率，但却不能很好地反映组织的生理学信息。然而肿瘤的生物学变化比形态学变化早^[1]，其中灌注是指血流通过毛细血管网将氧和营养物质传送给组织细胞的过程，在一定程度上能反映器官组织的血流动力学状态和功能情况，功能成像(或称灌注成像)即是通过影像学手段来直观显示活体灌注过程并进行定量或半定量分析^[2]，从而达到对肿瘤病变、复发以及癌细胞转移的早期检测、发现，进而实施治疗。

锥形束电子计算机断层扫描(Cone-Beam Computed Tomography, CBCT)系统^[3]作为放射治疗的重要设备，使用实时有源面阵探测器对检测物体进行一系列不同角度的射线测量，能够在高质量的诊断图像提供亚毫米的分辨率，相比常规CT扫描剂量降低为原来的1/15，由于空间分辨率、射线利用率高在医学应用领域有旷阔的应用前景^[4-5]。在同样的扫描参数与体层厚度上，CBCT扫描一周的时间为数十秒，相比动态对比增强CT(Dynamic Contrast Enhanced CT, DCE-CT)图像扫描一周只用1 s时间，CBCT的每个探测器单元能得到更多的信息量，即CBCT在重建图像的密度分辨率上占有更大优势，便于医学发现和诊断组织功能情况。

在CBCT灌注成像的研究中，文献^[6]针对C-臂CBCT数据采集速度相对缓慢，使用了一种时间恢复方法得到灌注增强曲线；文献^[7]研究了影响C-臂CBCT灌注成像的硬件因素；文献^[8]在经动脉化疗栓塞过程中，使用超快锥束成像数据，通过横断面成像、血容量灌注以及影像融合来进行肿瘤检测；文献^[9-10]分别考虑了射线穿过体素的像素值、射线覆盖体素的射线长度，设计了点、线积分模型获得投影数据，使用伽马函数来描述时间-密度曲线(Time-Density Curve, TDC)，利用半卷积模型获取灌注参数。本文在之前研究的基础上，探究使用面积积分模型获得投影数据，验证CBCT灌注成像模型的有效性。

1 功能CT与CBCT采样的问题

功能CT(functional CT, fCT)^[11-12]又叫CT灌注，它通过计算组织灌注参数来评价组织器官的灌注状态，进而在活体状态下，反映出活体内肿瘤微血管的变化和微循环功能。通过动态CT进行单位置扫描可以得到时间-密度曲线，通过去卷积^[13]方法，从主动脉增强的峰值分离出组织的时间-密度曲线的最大梯度，提供量化的组织灌注信息^[14]。由于CBCT扫描一周的时间为数十秒，无法获取某个位置下体素随着时间的变化，但是依据血流动力学理论，在2~3 min的时间内血流动力学参数值不会有明显变化，所以只要能得到每一体素的近似时间-密度曲线，仍可运用CBCT进行灌注成像。

为了验证CBCT可以进行功能成像，本研究采用每秒获得的DCE-CT图像来模拟CBCT投影数据。DCE-CT可1 s对物体完成一周扫描，获得投影数据，从而可求得每秒的实际体素时间-密度值；使用DCE-CT图像作为物体，分析CBCT投影的物理过程建立面积积分模型，得到模拟CBCT投影数据，使用CBCT体素时间模型获得模拟体素时间-密度曲线；通过对比实际体素时间-密度值与模拟体素时间-密度曲线，就可以验证CBCT是否可以进行功能成像。

2 CBCT功能成像方法 2.1 获取模拟投影数据

投影矩阵反映探测器上的投影与物体的关系，是影响投影图像的重要因素，在建立CT系统模型过程中，射线束与物体的体素之间的关系一般是用射线与体素的覆盖情况来描述的，常见的射线覆盖模型有点积分模型、线积分模型和面积积分模型等^[15]。面积积分模型(后简称面积模型)中将探测器的采样范围看作单位线段，将射线是看作射线源到探测器线段之间的面积，相对于点模型和线模型，面积模型有更好的图像质量，但是计算也更加复杂和耗时。

可以利用几何模型求解面积权值矩阵，如图 1^[16]所示，射线源到探测器形成一个狭窄的扇形束，每个体素的面积权值表示为：

${{w}_{ij}}={{S}_{ij}}/{{\delta }^{2}}$

(1)

其中：S_ij表示探测器i穿过在索引为j的体素的面积；δ²表示体素的面积值。某个体素最多可能被两条射线穿过，先看体素被单射线穿过的情况。将射线看作一个y＝mx+b的方程，根据m的值可将射线穿过的情况分为如图 2的六种。假设b=0，当-1≤m＜0时，射线从左到右穿过x轴，穿过的下一个体素索引的横坐标(j值，代表列)和纵坐标(i值，代表行)都递增；当0 ＜m≤1时，射线从左到右穿过x轴，横坐标递增而纵坐标递减；对于m＜-1和m>1的情况，将方程中的x和y变换位置，即转换为x＝m′y+b′，则m′为0＜m′≤1和-1≤m′＜0，射线从左到右穿过y轴，它们的纵坐标都是递增，对于m＜-1纵坐标递增而m>1纵坐标递减。当x＝0时，纵坐标递增而横坐标保持不变；当m＝0时，横坐标递增而纵坐标保持不变。

现在，我们的关注点在于计算一个射线穿过时的面积计算，观察发现，根据射线最后穿出体素的位置，可以将一条射线穿过的情况分为三种，以0＜m≤1为例，如图 3所示。

如图 3所示，射线穿过第一个体素k时，射线与左边缘相交的长度为d₀，与右边缘相交的长度为d，d＝d₀+Δ(Δ＝m×δ)。当d＜δ时(如图 3(a) )，L从体素k的右方边缘穿出，然后使d₀=d，下一个被穿过的体素为k+1；当d＝δ时(如图 3(b) )，L从体素k的右上方点穿出，然后使d₀=0，下一个被穿过的体素为k+1；当d>δ时(如图 3(c) )，L从体素k的上边缘穿出，穿过k-n，然后使d₀=d，下一个被穿过的体素为k-n+1。各部分的面积公式如下：

${{S}_{1}}=\left( 2d-\Delta \right)\times \delta /2$

(2)

${{S}_{2}}=\left( 2\delta -\Delta \right)\times \delta /2$

(3)

$\left\{ \begin{array}{*{35}{l}} {{S}_{3}}={{\delta }^{2}}-{{\left( \Delta -d \right)}^{2}}/\left( 2m \right) \\ {{S}_{4}}={{d}^{2}}/2m \\ \end{array} \right.$

(4)

根据以上总结的规律，使用伪代码呈现一条射线穿过体素的面积、索引的求解，如图 4所示，以0 ＜m≤1的情况为例，Ω代表一个n × n＝N的被扫描物体图像(Ω ={(x，y)|-D ≤ x ≤D；-D ≤ y ≤ D} )，它的每个体素k的索引由i和j表示(0≤i ＜n，0≤j＜n，k= i×n+j)，射线L₁从左边穿过Ω，首先应计算出被最先穿过的体素k₁、i₁、j₁，以及穿过k₁的左边缘长度d₁。经过计算得出：h=-m×D+b+D，i₁=n-1-(h/δ)，j₁=0，k₁=i₁×n+j₁，d₁=D×(1-m)+b-(n-1-i₁)×δ。现在使用{index}、{index_i}、{index_ j}、{area}分别储存射线穿过的k、i、j，以及所占的面积，伪代码1内容如下：

输入 射线穿过的初始体素k₁及其i₁、 j₁、d₁；

输出 数组{index}、{index_i}、{index_j}、{area} 。

Δ←m×δ, num←0, k←k₁, i←i₁, j←j₁, d←d₁

while i≥0 and j ＜n do

{ d←d+Δ

if d>δ then

d←d-δ, index←k, k←k-n

area←δ²-(Δ-d)²∕(2m)

index_i←i, index_ j←j

i←i-1, num←num+1

if i≥0 then

index←k, k←k+1

area←d²/(2m)

index_i←i, index_ j←j

num←num+1, j←j+1

else return num

else if d＜δ then

index←k, k←k+1

index_i←i, index_ j←j

area←(2d-Δ)×δ∕2

j←j+1, num←num+1

else

d←0, index←k, k←k-n+1

index_i←i, index_ j←j

area←(2δ-Δ)×δ∕2

j←j+1, i←i-1, num←num+1;

}

return num

当射线从底部穿过，如图 4的射线L₂所示，设定它穿过第一个像素k₂的面积为d₂²/(2m)，接下来的算法和上述伪代码1一致。对于-1≤m＜0的情况，除了坐标变化规律有差异，面积的计算方法是一致的，此处不再赘述。而m＜-1、m>1的情况，须将m变为m′，并注意坐标规律的变化。另需注意m＝0和x＝0的情况，以m＝0为例，每个被像素覆盖的像素的面积都设定为S=d₁×δ，伪代码2如下，将坐标变化规律进行修改可转变为x＝0的情况。

S←d₁×δ, num←0, k←k₁

for j←0 to n-1 step j←j+1

index←k, k←k+1

area←S

index_i←i₁, index_j←j

num←num+1

return num

当射线被两条射线穿过时，可以通过计算它们覆盖的面积的差，图 5为当L₁为-1≤m＜0，L₂为m＜-1时的情况，面积S的值为S₁-S₂。由此可见，两条射线穿过的情况在于计算一个射线的覆盖面积。经统计，两条射线穿过的情况一共有八种，如表 1所示。

由以上理论，对于一个探测器单元的投影，可能射线覆盖体素时会有三种情况：单射线覆盖、双射线覆盖、全覆盖(即图 7的黑色阴影处的体素)，方法是从左到右以列的方式进行计算，如图 6所示，列的索引为j(0≤j＜n)，L₁为处于上方的射线，L₂位于下方，设定U、V为L₁、L₂经过的体素数量。以列被穿过L₁、L₂的情况，可以将其分为三组：设j₀=min{j₁[0], j₂[0]}，j₁=max{j₁[U-1], j₂[V-1]}，j_min=max{j₁[0], j₂[0]}，j_max=min{j₁[U-1], j₂[V-1]}}。第一组C₁为j₀ ＜j_min时，列的体素仅有上方射线穿过；从j_min至j_max为第二组C₂，即列的体素由两条射线都穿过；若j_max ＜j₁，从j_max+1至j₁，列的体素仅有下方射线穿过，这为第三组C₃。

由于C₁、C₃组是只有一条射线穿过的情况，使用伪代码1讨论的方法可以容易地进行计算，故此处讨论C₂的情况。如图 7，设物体Ω被L₁、L₂穿过，根据伪代码1的算法，使用{index₁}、{i₁}、{j₁}、{area₁}和{index₂}、{i₂}、{j₂}、{area₂}分别储存射线L₁、L₂穿过的体素索引以及所占的面积。当j_min≤j≤j_max时，设定n₁和n₂分别为{j₁}和{j₂}的第一个索引index值，即有j=j₁[n₁]=j₂[n₂]，L₁将穿过k₁，k₁＝index[n₁]，s₁＝area[n₁]，然后使n₁=n₁+1，L₁将穿过下一个j=j₁[n₁]的体素即k₂，L₁将陆续穿过k₃至k₅；L₂将穿过k₆，k₆＝index[n₂]，s₆＝area[n₂]，同样使n₂=n₂+1，L₂将穿过下一个j=j₂[n₂]的体素即k₇。在j_min≤j≤j_max的区域中，可以看出有四种覆盖情况：

1) 只被L₁穿过，即如图 7的k₁至k₃、k₅(浅灰色区域)，它们的面积权值为s/δ²；

2) 只被L₂穿过，即如图 7的k₆至k₁₁(深灰色区域)，它们的面积权值为1-s/δ²；

3) 同时被L₁、L₂穿过，即如图 7的k₄＝k₁₂，它的面积权值为(s₄-s₁₂)/δ²(斜底纹区域)；

4) 都不会被L₁、L₂穿过，即如图 7的黑色阴影区域。

由此给出伪代码3的算法，使S=δ²，p、q分别代表同一列中，被L₁、L₂穿过的体素数量(如在j_max列中，k₄、k₅、k₁₀、k₁₁位于同一列，p和q的值为2) ，将单射线覆盖、双射线覆盖、全覆盖三种情况下的面积权值存至{weight}，索引存至{pixel_index}，在计算过程中它们目前的项的个数为num。

//S为体素的单位面积，j为穿过的体素列索引，j_min≤j≤j_max

//n为图像总列数

//当前j下L₁、L₂首次穿过的体素列索引j₁、j₂、索引n₁、n₂

//U、V为L₁、L₂经过的体素数量

// j_min=max{j₁, j₂}，j_max=min{j₁, j₂}

//p、q为同一列中L₁、L₂穿过的体素数量

//index_max为L₁穿过的索引最大值，index_min为L₂穿过的索

//引最小值

//面积权值数组{weight}，索引数组{pixel_index}，{weight}或//{pixel_index}的个数num

for j←j_min to j_max step j←j+1

p←0, q←0

//当L₁经过当前j的首体素

if j₁=j then

pixel_index←index₁

weight←area₁∕S

k₁←num, q←q+1, num←num+1, n₁←n₁+1

//继续计算当前j下L₁经过的非首体素

if n₁＜U and j₁=j then

pixel_index←index₁

weight←area₁∕S

k₂←num, q←q+1, num←num+1, n₁←n₁+1

//L₁只穿过一个体素

if q=1 then index_max←pixel_index

else index_max←max{pixel_index, pixel_index}

//当L₂经过当前j的首体素

if j₂=j then

//当L₁的首体素也是此体素

if index₂=pixel_index then

weight←weight-area₂∕S

k₃←k₁, p←p+1, n₂←n₂+1

else

pixel_index←index₂

weight←1-area₂∕S

k₃←num, p←p+1, num←num+1, n₂←n₂+1

//继续计算当前j下L₂经过的非首体素

if n₂＜V and j₂=j then

//当L₁的首体素也是此体素

if index₂=pixel_index then

weight←weight-area₂∕S

k₄←k₁, p←p+1, n₂←n₂+1

//当L₁的非首体素也是此体素

else if index₂=pixel_index then

weight←weight-area₂∕S

k₄←k₂, p←p+1, n₂←n₂+1

else

pixel_index←index₂

weight←1-area₂∕S

k₄←num, p←p+1, num←num+1, n₂←n₂+1

//L₂只穿过一个体素

if p=1 then index_min←pixel_index

else index_min←min{pixel_index, pixel_index}

//当L₁、L₂同时穿过的列超过1

if index_min-index_max>n then

for k←index_max+n to index_min-n step k←k+n

pixel_index←k

weight←1, num←num+1

2.2 体素时间-密度模型

放射治疗过程中，组织的血流动力学特征一般不会出现明显的变化。功能CT中对比剂团注下获得的时间-密度曲线通常可用伽马函数来描述^[17]。利用获取的CBCT的投影数据，对每一个体素的时间-密度变化进行数学建模，利用最优化参数求解技术拟合求解模型参数，描绘出各个体素的TDC曲线，用于观察组织血流灌注过程。

文献^[9]中引入反正切函数来模拟TDC曲线到达峰值后的变化趋势，以表征生物体实际的血流动力学情况，模型表达式为式(5) ：

$X(t)=\left\{ \begin{matrix} P \\ {{(t-T)}^{a}}\times {{e}^{-b(t-T)}}+c\times \arctan (t-T)+P \\ \end{matrix}\begin{matrix} {} \\ {} \\ \end{matrix}\begin{matrix} t ＜T \\ t\ge \text{T} \\ \end{matrix} \right.$

(5)

其中：X(t)是体素随时间变化的密度值；P表示未注入对比剂时物体的密度；t是注入对比剂后经过的时间；T是对比剂进入扫描层的时间；a、b、c是模型参数。

2.3 去卷积模型求灌注参数

利用现有的成熟血流动力学模型--去卷积模型，计算出血流(Blood Flow, BF)、血容(Blood Volume, BV)、平均通过时间(Mean Transit Time, MTT)、表面通透性(Permeability of capillary vessel Surface, PS)等功能参数值。

如图 8的模型将毛细血管看作两个同心圆柱体，代表着血管外空间(即细胞间隙)和血管内空间。血管外空间是对比剂充分混合的区域，对比剂浓度C_e(t)随着时间变化，血管内空间的对比剂浓度C_b(x, t)代表对比剂从动脉到毛细血管过程中的减少量。C_a(t)是动脉TDC，即动脉的对比剂浓度，C_v(t)为静脉出口处的局部对比剂浓度;F代表它们的体积流量;V_b和V_e分别为血管内空间和血管外空间空间的容积;PS为毛细血管表面通透性，与血流大小决定了对比剂由血管外空间向血管内空间的流量。

较短时间内，血管外空间内的浓度保持不变，C_e(t)可看为常数，此时血管内空间的对比剂浓度C_b(x, t)可看作功能参数和时间的函数。因此在去卷积模型^[18]中，使用时间域可将每单位组织对比剂的量表示如下：

$Q\left( t \right)=BF\cdot {{C}_{a}}\left( t \right)\otimes R\left( t \right)={{C}_{a}}\left( t \right)\otimes BF\cdot R\left( t \right)$

(6)

其中：Q(t)是组织TDC，即组织驻留函数；BF·R(t)是血流量与脉冲驻留函数的乘积。R(t)的表达式为式(7) ：

$R(t)=\left\{ \begin{align} & 1, \\ & {{E}^{2}}{{e}^{-\frac{BF\cdot E}{{{V}_{e}}}(2t-3MTT)}}, \\ \end{align} \right.\begin{matrix} 0＜t＜MTT \\ t\ge MTT \\ \end{matrix}$

(7)

其中，E与PS和BF相关，公式为：

$E=1-{{e}^{-PS/BF}}$

(8)

若R(t)和C_a(t)已知，则可以通过式(1) 计算Q(t)。CT能够检测的是Q(t)和C_a(t)，R(t)常常不易检测，使用Q(t)和C_a(t)来计算R(t)的过程便是卷积的逆运算，所以称为去卷积。经过去卷积运算得到BF·R(t)，它的平台高度即为BF，R(t)曲线下面部分的面积和即为BV，BV/BF的值便是MTT。

3 实验结果及分析 3.1 实验材料与准备

实验图像为多伦多玛格丽特公主医院提供的在腿部植入VX2型肿瘤的新西兰大白兔，肿瘤植入两周后进行DEC-CT扫描获得的二维序列图像。CT机型：GE Lightspeed Plus 4层螺旋CT，使用GE MEDICAL SYSTEMS，管电压80kVp，管电流120 mA，20 cm VOF，扫描速度为每周1 s，轴向4层，层厚5 mm，重建图像分辨率为0.39 mm，像素矩阵大小为512×512。使用Isovue300对比剂，按1.5 mL/kg体质量(总量不超过4.5 mL)10 s以上经耳静脉连续注入。扫描周期为1 s，整个扫描持续120 s。共得到四组120幅DICOM格式DCE-CT序列图像。如图 9所示，虚线椭圆圈住区域是肿瘤区域，实线椭圆圈住的是动脉血管。由于计算速度的原因，选择如图 10的肿瘤区域(9×9像素大小，第二层切面中)为研究对象。

3.2 投影数据

使用Matlab R2014b分别实现点^[9]、线^[10]、面积模型^[16]的投影数据模拟方法，进行基于DCE-CT序列图像的CBCT投影仿真，得到如图 11所示的CBCT投影正弦图，可见面积模型呈现出更多细节，且减少了大部分锯齿状伪影。

3.3 模拟TDC与实际TDC

以Matlab R2014b为工具，利用最优化求解技术，在所用材料数据的合理参数范围：a(2.50～3.00) 、b(0.180～0.210) 、c(15～35) 、T(10～20) 内变化各体素的模型参数值。在上述参数范围内以线性遍历的方式计算出不同的X(t)值，当在各投影方向上的体素点的X(t)值和与投影到探测器上的投影值最接近时，确定此参数组为各体素点的最佳参数组。下面对实验最终得到的模型计算出的TDC与实测的TDC的偏差进行评估。

以第二个切面的图 10的9×9区域为例，由模型获得的TDC(黑色虚线)与DET-CT扫描获得TDC(黑色实线)对比，如图 12(a)体素索引为37~45(肿瘤中心)的TDC，图 12(b)体素索引为64～72(正常组织)的TDC。

用相关系数来评价两个曲线之间的相似度。以第二个切面的图 10的9×9区域81个体素点为例，计算出三种模型下得到的投影数据，它们的真实TDC与模拟TDC间的相关系数值，如图 13所示，77～81的误差可能是原图像测量误差所致。其中，点模型、线模型、面积模型的相关系数的平均值分别为0.804024983、0.825983271、0.829104662，可见面积模型所求的投影数据建立的TDC曲线更加符合实际图像。

3.4 灌注参数与伪彩灌注图

使用图 10的肿瘤区域，通过面积模型下的投影数据得到动脉TDC与模拟TDC，计算组织的灌注参数。使用2.3节的模型，利用Matlab工具，计算两条TDC的残差和来求解最优化灌注参数：设置生理参数的初始值，得到R(t)后与C_a作卷积，再与BF相乘，将结果与组织实测的TDC比较。当两者最吻合时，提取此时参数计算血流(BF)、血容(BV)、表面通透性(PS)、平均通过时间(MTT)，作为此区域的功能参数。

图 14为分别为对MTT、PS、BV和BF的第一、四个切片的值作伪彩图得到的功能图像。选用的第一切片更靠近肿瘤区域，组织的体素间区域功能参数值差异较大，说明肿瘤生长活跃，在放射治疗中可用大剂量进行照射。第四切片靠近正常组织，在解剖图像中形态上保持肿瘤区域，但是由图可见其功能性十分低，可能其中肿瘤已经坏死，对这些区域可少剂量照射观察或不照射。功能CT成像技术可以监测到组织的血流动力学变化，实施符合实际的放疗计划。

4 结语

本文利用DEC-CT图像，使用投影矩阵模拟得到投影数据，对比投影矩阵中点模型、线模型、面积模型的功能成像效果，以验证锥形束CT功能成像的可行性。通过血动力学建模分析，使用最优化数值计算方法求得基于投影数据的各体素密度随时间变化的模型参数，仿真出组织的模拟TDC与实际TDC值，利用相关系数对曲线的拟合程度进行定量分析，发现经面积模型得到的模拟数据能更贴近实际值，同时验证了TDC模型的有效性；使用去卷积模型以及计算机最优化技术获取组织血动力学的灌注参数，并使用伪彩图展示，呈现出锥形束CT功能成像效果，发现其符合图像实际情况。实验结果表明本文所提出的方法以及所用模型的正确性，利用模拟投影数据证明所用功能模型可以实现在线的锥形束CT功能成像以及功能图像引导的放射治疗。