一种用于锥束CT图像散射修正的阻挡光栅优化方法及装置与流程

本发明涉及锥束CT应用相关的医学成像领域，尤其涉及一种用于锥束CT图像散射修正的阻挡光栅优化方法及装置。

背景技术：

本部分旨在为权利要求书中陈述的本发明的实施方式提供背景或上下文。此处的描述不因为包括在本部分中就承认是现有技术。

影像引导放射治疗技术(image guide radiation therapy,IGRT)是目前精准放疗中的最主要手段。IGRT可以提供精确的位置信息用于精准病灶定位和肿瘤照射，已经在当代微创手术和放射治疗中广泛使用。安装在治疗床边的X射线CBCT是当代图像引导的重要实现手段。

然而，伴随着单次投影中辐照体积的扩大，CBCT(cone-beam computed tomography,CBCT)图像的质量因散射污染这一基本物理过程的限制而急剧恶化。因为传统的CT重建理论假设x射线沿着直线传播，探测到的射线强度随着路径积分指数衰减。而散射光子偏离入射束方向，不能被建模到传统CT重建理论中，因而成为了CT图像重建的误差来源。研究发现，散射信号强度随着x射线辐照体积的增加而单调增强。CBCT的散射污染严重影响了其CT值精度、低对比度目标的检测能力及剂量计算的准确性，这些缺点直接导致了CBCT无法在临床上广泛应用。在没有散射修正的CBCT系统中对人体进行扫描，散射伪影导致的CT值误差可达350HU，因此CBCT主要还是用于初步定位和摆位，在介入和放射治疗中的进一步应用受到严重限制。散射修正是提高CBCT图像质量首要解决的问题。

目前已知的散射修正的方法主要可以分为两类：预处理和后处理方法。预处理方法校正散射主要是通过附加硬件装置，阻止散射光子到达探测器，这样散射信号以及相关联的统计噪声就被抑制在投影中。两个典型的预处理方法的例子包括，增加探测器与物体之间的空气隙和使用抗散射线栅。随着空气隙加宽，扩散开的散射光子的探测率会降低，而源信号则不受影响。但该方法受CBCT设备本身物理空间的限制，同时还需要增加X线剂量来弥补距离的增加，在临床实际中并不实用。抗散射线栅使用聚焦于射线源的铅栅网格，能阻挡非聚焦入射角的散射光。该方法也存在对散射光的衰减效率不高的缺陷。目前商用线栅只能提供约3倍的SPR降低率，无法保证高散射环境下的CBCT图像质量。此外，它还需要增加病人的受照剂量来补偿被衰减掉的源射线强度。

鉴于预处理方法的局限性，目前的研究更多的集中在后处理上。后处理是指按传统方式获取散射污染的投影图像后再做散射修正。由于不可能从理论上预测随机散射时间，即使采用了完美的后处理方法，散射噪声还是会留在图像中。后处理方法有很多种，包括：解析建模法、蒙特卡罗模拟法、源调制法和测量法。解析建模法认为散射信号是源信号通过散射核后的响应，散射核一般是经过测量或模拟获得。若散射核具有线性平移空间不变特性，则其计算速度很快；但是相应的散射估计精度有限且对复杂物体需要繁琐的调整参数。蒙特卡罗模拟法通过模拟光子与受照物体之间的相互作用，为散射信号建立更加精确的统计学模型，但该方法计算量极大，十分耗时，目前的计算机运算能力限制了其在要求近实时处理的CBCT图像重建中的应用。源调制是在X射线源和物体之间加入高频调制器，根据散射和源信号不同的响应特性，在频率域上把它们分离开。该方法不增加病患曝光剂量和扫描时间，但其临床应用效果受到实际物理因素的限制，如旋臂震动和大焦点尺寸等。

基于测量的散射修正是与本发明最为相似的实现方法。该方法在CBCT射线源前方加入射束阻挡光栅用于估计散射信号，这样探测器上形成仅包含散射信号的阴影区。由于散射分布主要是低频成分且受阻挡光栅扰动很小，则全野的散射分布可以对探测器阴影区的散射抽样插值得出。基于测量的方法可以获得精确的散射估计，但这是以源信号损失为代价的。因此，人们通常需要对每个角度扫描两次(一次有阻挡光栅，另一次挪开阻挡光栅)，或者在扫描过程中移动阻挡光栅。

阻挡器的低成本与易制作性，使得基于测量的散射修正方法成为可能。在中国专利201410380731.3中，提出了一种基于互补光栅的锥束CT散射校正方法及其装置。该发明，利用互补光栅扫描及少量的计算，实现投影图像的散射校正，利用散射校正后的投影图像就可以重建出散射修正切片图像。在中国专利201010574162.8中，提出一种CT系统的散射校正方法及CT系统，通过获取亮场图像，将阻挡器放置在探测器和被扫物体之间，扫描后得到衰减投影图像，再分别对被扫描物体和散射校正器进行扫描得到投影图像集和散射校正图像。再根据亮场图像、散射校正图像以及衰减投影图像生成散射信号分布，最后通过投影图像集与散射信号分布之差得到修正后的投影。

现有利用阻挡器进行散射修正方法中，大多都需要两次扫描以补偿阻挡的原始投影信号，因其增加了病人的受照剂量，不能付诸临床应用。此外，这些方法都用在桌面CBCT实验平台，由于实验平台采用扫描物体旋转替代光源和探测器的旋转，CT系统的等中心很稳定，阻挡器在不同时刻投影的位置几乎保持一致，这样利用阻挡器进行散射修正的难度就会大大降低，这些方法只适用于等中心稳定的实验室桌面CBCT系统。然而，临床CBCT系统由于机架旋转过程中的震动以及旋转等中心偏差，导致阻挡光栅投影位置随机架旋转而变化，造成光栅区域的散射信号难以准确提取。

技术实现要素：

本发明的目的是提供一种用于锥束CT图像散射修正的阻挡光栅优化方法及装置，以解决现有的利用阻挡器进行散射修正时增加扫描次数且散射信号难以准确提取的问题。

为了达到上述目的，本发明实施例提供一种用于锥束CT图像散射修正的阻挡光栅优化方法，包括：扫描阻挡光栅，建立所述阻挡光栅的摆动模型；设定初始投影中所述阻挡光栅在探测器纵向上的初始坐标，并根据所述摆动模型，建立CBCT图像数据缺失体素值与阻挡光栅在探测器纵向上的坐标之间的目标函数；利用网格自适应直接搜索算法最小化所述目标函数，生成所述阻挡光栅在探测器纵向上的优化坐标。

进一步地，在一实施例中，所述扫描阻挡光栅，建立所述阻挡光栅的摆动模型，具体包括：在所述阻挡光栅的投影上，使用最大类间方差法确定所述阻挡光栅的一个合适阈值，对其进行基于阈值的图像分割，生成二值图像，通过生成的二值图像确定所述阻挡光栅的坐标位置，得到所述阻挡光栅的摆动模型。

进一步地，在一实施例中，所述设定初始投影中所述阻挡光栅在探测器纵向上的初始坐标，包括：设定所述阻挡光栅的片数为n，初始投影中第i片阻挡光栅在探测器纵向上的坐标为：G＝(g₁,g₂,…,g_n)^T。

进一步地，在一实施例中，所述根据所述摆动模型，建立CBCT图像数据缺失体素值与阻挡光栅在探测器纵向上的坐标之间的目标函数，具体包括：根据所述摆动模型以及初始坐标，模拟所述阻挡光栅在每幅投影的投影图像；对所述每幅投影图像进行反投影重建；根据反投影重建后的图像，建立左右半扇重建图像融合后数据缺失的体素值与所述阻挡光栅在探测器纵向上的坐标之间的目标函数。

进一步地，在一实施例中，所述利用网格自适应直接搜索算法最小化所述目标函数，生成所述阻挡光栅在探测器纵向上的优化坐标，包括：用初始设定的阻挡光栅的参数作为网格自适应直接搜索算法的初始值，生成所述阻挡光栅在探测器纵向上的优化坐标。

进一步地，在一实施例中，所述初始设定的阻挡光栅为等间隔均匀分布。

进一步地，在一实施例中，所述阻挡光栅为指交叉型光栅。

为了达到上述目的，本发明实施例还提供一种用于锥束CT图像散射修正的阻挡光栅优化装置，包括：摆动模型建立模块，用于扫描阻挡光栅，建立所述阻挡光栅的摆动模型；目标函数建立模块，用于设定初始投影中所述阻挡光栅在探测器纵向上的初始坐标，并根据所述摆动模型，建立CBCT图像数据缺失体素值与阻挡光栅在探测器纵向上的坐标之间的目标函数；最小化模块，用于利用网格自适应直接搜索算法最小化所述目标函数，生成所述阻挡光栅在探测器纵向上的优化坐标。

进一步地，在一实施例中，所述摆动模型建立模块用于扫描阻挡光栅，建立所述阻挡光栅的摆动模型，具体包括：在所述阻挡光栅的投影上，使用最大类间方差法确定所述阻挡光栅的一个合适阈值，对其进行基于阈值的图像分割，生成二值图像，通过生成的二值图像确定所述阻挡光栅的坐标位置，得到所述阻挡光栅的摆动模型。

进一步地，在一实施例中，所述目标函数建立模块具体包括：初始坐标设定单元，用于设定所述阻挡光栅的片数为n，初始投影中第i片阻挡光栅在探测器纵向上的坐标为：G＝(g₁,g₂,…,g_n)^T；投影图像模拟单元，用于根据所述摆动模型以及初始坐标，模拟所述阻挡光栅在每幅投影的投影图像；反投影单元，用于对所述每幅投影图像进行反投影重建；函数建立单元，用于根据反投影重建后的图像，建立左右半扇重建图像融合后数据缺失的体素值与所述阻挡光栅在探测器纵向上的坐标之间的目标函数。

本发明实施例的用于锥束CT图像散射修正的阻挡光栅优化方法及装置，提出了一种全新的无需源补偿而且适合于临床CBCT的散射修正方法。本发明借助图像分割方法建立阻挡光栅投影因悬臂等中心偏移与机架震动引起的摆动数学模型，成功将阻挡光栅应用于临床锥束CT散射修正；对阻挡光栅在图像域造成的数据缺失进行数学优化建模，定量评估阻挡器对重建图像的影响，借助网格自适应直接搜索法求解阻挡器的几何最优结构，为基于阻挡器测量的散射修正方法建立了坚实的理论基础，进一步揭示阻挡器设计对临床锥束CT散射修正的重要性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例的用于锥束CT图像散射修正的阻挡光栅优化方法的处理流程图；

图2为本发明实施例的用于锥束CT图像散射修正的阻挡光栅优化方法的原理示意图；

图3为图2所示实施例中得到的优化后的阻挡光栅的几何模型；

图4为本发明实施例的用于锥束CT图像散射修正的阻挡光栅优化装置的结构示意图；

图5为图4所示实施例中的目标函数建立模块102的结构示意图；

图6为依据本发明的阻挡光栅优化方法得到的实际修正效果对比示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本领域技术技术人员知道，本发明的实施方式可以实现为一种系统、装置、设备、方法或计算机程序产品。因此，本公开可以具体实现为以下形式，即：完全的硬件、完全的软件(包括固件、驻留软件、微代码等)，或者硬件和软件结合的形式。

下面参考本发明的若干代表性实施方式，详细阐释本发明的原理和精神。

本发明实施例提出“指交叉”型阻挡光栅对重建图像质量的数学模型，该数学模型将光栅抖动考虑进去，引入网格自适应直接搜索的方法求解目标函数，得到适合于临床CBCT散射修正的阻挡光栅几何设计。然后将设计的阻挡光栅放置在射线源前，准确提取散射样本后通过插值法估计每个投影的散射分布，最后通过半扇扫描重建算法精确重建出散射修正后的图像，单次扫描实现临床CBCT散射修正。

图1为本发明实施例的用于锥束CT图像散射修正的阻挡光栅优化方法的处理流程图。如图1所示，包括：

步骤S101，扫描阻挡光栅，建立所述阻挡光栅的摆动模型；

步骤S102，设定初始投影中所述阻挡光栅在探测器纵向上的初始坐标，并根据所述摆动模型，建立CBCT图像数据缺失体素值与阻挡光栅在探测器纵向上的坐标之间的目标函数；

步骤S103，利用网格自适应直接搜索算法最小化所述目标函数，生成所述阻挡光栅在探测器纵向上的优化坐标。

具体实施时，参考图2，为本发明实施例的用于锥束CT图像散射修正的阻挡光栅优化方法的原理示意图。图1的步骤S101对应图2的100(建立光栅摆动模型)，图1的步骤S102对应图2的200(阻挡光栅数学建模)，图1的步骤S103对应图2的300(模型的求解)。在图2所示实施例中，采用10根光栅为例进行原理说明。

在本实施例的步骤S101中，参看图2的100，首先需要利用经验设计一个阻挡光栅，然后对其进行扫描，用于研究铅条投影在该临床CBCT因悬臂等中心偏移与机架震动引起的摆动模型。在一较佳实施例中，所设计的阻挡光栅为“指交叉”光栅。

在扫描得到的阻挡光栅的投影上，使用最大类间方差法确定所述阻挡光栅的一个合适阈值，对其进行基于阈值的图像分割，将感兴趣区域转化为一副二值图像，可得到阻挡光栅中心点在该方向上的坐标值，获得机架角度与阻挡光栅位置变化的关系，即，通过生成的二值图像确定所述阻挡光栅的坐标位置，根据光栅的坐标位置可以得到所述阻挡光栅的摆动模型。

在本实施例的步骤S102中，参看图2的200，要建立阻挡光栅与图像数据缺失体素数量的目标函数。其原理为：当阻挡光栅间隔越宽，通过两组半扇扫描算法重建出的CBCT图像数据缺失的体素越少。在保证散射估计准确性的同时，阻挡光栅的设计需要保证图像数据缺失体素数量的最小化，因此可以设计一个关于阻挡光栅放置位置的目标函数。

首先要设定初始投影中所述阻挡光栅在探测器纵向上的初始坐标。设定所述阻挡光栅的片数为n，初始投影中第i片阻挡光栅在探测器纵向上的坐标即为：

G＝(g₁,g₂,…,g_n)^T (1)

然后根据所述摆动模型以及初始坐标，模拟所述阻挡光栅在每幅投影的投影图像。具体实施时，模拟阻挡光栅的投影图像，令阻挡区域为1，未被阻挡的区域为0。由于摆动和机架等中心偏差的原因，阻挡光栅在每幅投影的坐标位置并不相同，d_j为阻挡光栅在第j幅投影偏移量。

式(2)和式(3)中的p_i,j(u,v)分别为左、右侧单片阻挡光栅的在第j幅投影的图像，U为探测器横向像素宽度，w₀为阻挡光栅宽度，u,v分别为探测器的横、纵坐标值。

将式(2)和式(3)相加，得到：

P_j为阻挡光栅的在第j幅投影的图像。

对m幅投影P_j分别乘上加权函数f(u,v)和f(-u,v)进行反投影重建，得到：

可以注意到，M_l和M_r中受到阻挡光栅影响的体素值为非零，而未受到阻挡光栅影响的体素值为零，M_l·M_r为非零的区域即为左右半扇重建图像融合后数据缺失的体素，阻挡光栅的优化设计要求重建体积中数据缺失的体素值之和(需要二值化)的最小化。目标函数建立过程如图2所示，可以得到下述的目标函数：

其中，为约束条件，表示左右光栅必须交错分布，且光栅的纵坐标不能超出探测器纵坐标的范围。

在本实施例的步骤S103中，参看图2的300，用初始设定的阻挡光栅的参数作为网格自适应直接搜索算法的初始值，生成所述阻挡光栅在探测器纵向上的优化坐标。其中，所述初始设定的阻挡光栅的参数为等间隔均匀分布，当然也可以采用其他间隔分布。

引入网格自适应直接搜索的方法，求解目标函数。本文所提优化模型中，震动和机架等中心偏差的引起的阻挡光栅震动的偏移量d_j为目标跟踪所得。目前针对最优化问题的求解方法有数学规划法、启发式算法和直接搜索法等。阻挡光栅的优化模型具有很强的非线性，并且没有导数信息可以提供，因此不能使用数学规划法求解；而像启发式算法例如模拟退火算法和遗传算法等虽然在全局搜索能力上更胜一筹，但是在局部的搜索能力却存在不足，收敛的速度也很慢。由于每生成一次目标函数，需要进行一次反投影操作，若本发明的目标函数使用启发式算法，会出现大量的冗余迭代，耗时长。网格自适应直接搜索(mesh-adaptive direct search,MADS)算法的决策向量可以是离散的、连续的和二进制等，目标函数和它的约束条件可以是“黑盒函数”，MADS算法适合求解多变量的数学模型，因此本发明采用MADS算法对阻挡光栅优化模型进行求解。求解时，在一较佳实施例中，用均匀分布阻挡光栅的参数作为MADS目标函数的初始值，求解得出的光栅几何模型如图3所示。

在通过以上步骤得到优化的阻挡光栅结构后，将其放置在射线源和被扫描物体之间，采集物体投影后，利用基于阈值的图像分割方法确定光栅所在的位置后，利用所采集的散射信号，插值估计出散射分布，原始图像减去散射分布后，得到散射修正后的投影图像。然后再利用基于Parker函数的半扇扫描重建算法，重建出散射修正后的图像。

应当注意，尽管在附图中以特定顺序描述了本发明方法的操作，但是，这并非要求或者暗示必须按照该特定顺序来执行这些操作，或是必须执行全部所示的操作才能实现期望的结果。附加地或备选地，可以省略某些步骤，将多个步骤合并为一个步骤执行，和/或将一个步骤分解为多个步骤执行。

在介绍了本发明示例性实施方式的方法之后，接下来，参考图4对本发明示例性实施方式的用于锥束CT图像散射修正的阻挡光栅优化装置进行介绍。该装置的实施可以参见上述方法的实施，重复之处不再赘述。以下所使用的术语“模块”和“单元”，可以是实现预定功能的软件和/或硬件。尽管以下实施例所描述的模块较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。

图4为本发明实施例的用于锥束CT图像散射修正的阻挡光栅优化装置的结构示意图。如图4所示，包括：

摆动模型建立模块101，用于扫描阻挡光栅，建立所述阻挡光栅的摆动模型；目标函数建立模块102，用于设定初始投影中所述阻挡光栅在探测器纵向上的初始坐标，并根据所述摆动模型，建立CBCT图像数据缺失体素值与阻挡光栅在探测器纵向上的坐标之间的目标函数；最小化模块103，用于利用网格自适应直接搜索算法最小化所述目标函数，生成所述阻挡光栅在探测器纵向上的优化坐标。

进一步地，在本实施例中，所述摆动模型建立模块101用于扫描阻挡光栅，建立所述阻挡光栅的摆动模型，具体包括：在所述阻挡光栅的投影上，使用最大类间方差法确定所述阻挡光栅的一个合适阈值，对其进行基于阈值的图像分割，生成二值图像，通过生成的二值图像确定所述阻挡光栅的坐标位置，得到所述阻挡光栅的摆动模型。

进一步地，在本实施例中，如图5所示，所述目标函数建立模块102具体包括：

初始坐标设定单元1021，用于设定所述阻挡光栅的片数为n，初始投影中第i片阻挡光栅在探测器纵向上的坐标为：G＝(g₁,g₂,…,g_n)^T；

投影图像模拟单元1022，用于根据所述摆动模型以及初始坐标，模拟所述阻挡光栅在每幅投影的投影图像；

反投影单元1023，用于对所述每幅投影图像进行反投影重建；

函数建立单元1024，用于根据反投影重建后的图像，建立左右半扇重建图像融合后数据缺失的体素值与所述阻挡光栅在探测器纵向上的坐标之间的目标函数。

此外，尽管在上文详细描述中提及用于锥束CT图像散射修正的阻挡光栅优化装置的若干单元，但是这种划分仅仅并非强制性的。实际上，根据本发明的实施方式，上文描述的两个或更多单元的特征和功能可以在一个单元中具体化。同样，上文描述的一个单元的特征和功能也可以进一步划分为由多个单元来具体化。

通过本发明的实施方式进行设计优化的阻挡光栅已经过实验、模拟和使用，证明确实可行。在瓦里安Trilogy的CBCT上，利用Catphan504模体进行实验，采用本发明设计优化的阻挡光栅进行散射修正，感兴趣区域的CT误差由115HU下降至11HU，对比度提高了1.45倍，修正效果如图6所示，(a)列为修正前的图像，(b)列为修正后的图像。

本发明实施例的用于锥束CT图像散射修正的阻挡光栅优化方法及装置，提出了一种全新的无需源补偿而且适合于临床CBCT的散射修正方法。本发明借助图像分割方法建立阻挡光栅投影因悬臂等中心偏移与机架震动引起的摆动数学模型，成功将阻挡光栅应用于临床锥束CT散射修正；对阻挡光栅在图像域造成的数据缺失进行数学优化建模，定量评估阻挡器对重建图像的影响，借助网格自适应直接搜索法求解阻挡器的几何最优结构，为基于阻挡器测量的散射修正方法建立了坚实的理论基础，进一步揭示阻挡器设计对临床锥束CT散射修正的重要性。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。