CN104598909A - X射线图像中的标志识别和处理 - Google Patents

Abstract

本发明涉及X射线图像中的标志识别和处理。一种用于在投影图像中快速地检测和定位放置于患者上的外部标志的机制。本发明的实施例即使在存在周围密集的解剖结构和大量的患者运动时也能够检测标志。一旦投影图像上的标志点的位置被提取出，这些标志点就可以用以执行基于标志的患者运动检测。本发明的实施例也可以用于运动修正之外，例如用于扫描器校准、自动头部测量、以及质量控制评估。

Description

X射线图像中的标志识别和处理

本申请是申请日为2010年8月26日、申请号为201080062012.7且题为“X射线图像中的标志识别和处理”的中国专利申请的分案申请。

技术领域

本发明的实施例涉及图像处理。具体而言，实施例提供用于对位于X射线图像中的标志关联点进行检测、定位及分类的方法和系统。

背景技术

现代成像系统在单个过程或扫描期间获取许多图像。使用已知为重建的过程来组合图像，以制作表示穿过患者的有限宽度“切片”的一个或多个图像。因为在为医疗或牙科过程做准备时或者为了适当地校准成像设备，经常要对所得图像加以分析，所以重要的是图像尽可能地清楚。然而，扫描期间患者移动会导致图像中的误差或模糊。

在成像过程之前常将标志放置于患者上。这些标志在所得图像上显现为标志点，并用以识别图像中的特定特征。标志点也可以用以针对患者移动来修正图像。在典型的运动修正情形中，必须于在扫描期间捕捉的数百个(例如，600个)或更多图像中的每个中识别相对较小数量的标志(例如，7个)。因而，使用标志来检测及修正患者运动仅在能够在图像中自动识别数千个标志点(例如，使用以上的示例为4200个标志点)的时候才是可行的。标志点识别由于X射线的非线性特性以及X射线图像的相对较大的动态范围而进一步被复杂化。

又一复杂之处在于，标志点识别和所得运动修正需要快速执行，这是因为其必须在图像的重建可以开始之前完成。现有的用于识别图像上的标志相似点的方法依赖于这些点的一般为圆形的特性，并采用诸如图像梯度处理或者霍夫(Hough)变换的一些变型的点识别技术，这些技术每图像要花费30秒之久才能完成。即使最简单的方法也要花费每图像1/10秒或更多，且经常要花费每图像5秒之久。

这些运动修正方法的准确度与提取标志点位置的准确度直接相关。例如，在大多数修正方法中，标志必须定位至像素的1/4至1/8的精度。这个程度的准确度必须在存在宽的背景动态范围值、对图像背景的非线性依赖性、以及噪声的情况下实现。许多现有的修正方法不能正确处理图像背景的非线性特性，导致图像边缘处的位置误差。这些现有的方法当图像中的标志彼此相对较为靠近时也不能很好地工作。

一旦标志点在图像中被识别出，这些标志点也就需要被映射回其所源自于的物理标志。物理标志通过它们的三维(“3D”)位置坐标而被唯一识别。因而，为了制作与特定标志相关联的每个投影图像的点的列表，一个步骤是识别每个标志的物理位置。这可以通过识别与由标志所产生的点的位置如何随着台架(gantry)转动而在投影图像之间的变化最为一致的每个标志的3D位置来进行。然而，因为患者移动或者计划外的系统或设备移动，所以每个标志的关联点位置会与其预期的位置明显改变。因此，识别标志的三维坐标的准确度在存在噪声或患者运动或者当多个标志位置彼此靠近时会受到损害。现有的用于识别标志的三维位置的技术可使用依赖于位置的电磁场，其由远程传感器检测并转换成位置。除该发射器和传感器系统的复杂度以外，系统所使用的技术也依赖于多个传感器，以确定标志的三维位置。其他定位技术包括光学技术，这类技术依赖于对多个被检测部位之间的固定关系的了解，从而增加了技术的复杂度。

即使在与图像点相关联的物理标志被识别，并且位置被分配给每个之后，将每个图像上的每个点分配给其原标志时仍然存在问题。在此完成之前，所有图像中的所有点仅作为单个列表存在。点分配(标志跟踪)试图将每个图像上的每个标志与实际生成它的标志相关联。在完成点分配之后，每个标志将处理与该标志相关联的标志点的分离列表。每个标志的分离列表可具有每个投影图像上的仅一个或零个标志点。重要的是，标志与其关联的标志点之间的正确对应关系即使在存在运动和测量误差的情况下也要被维持。

发明内容

实施例提供用于快速识别图像中的标志关联点的系统和方法，其然后可以用于运动修正系统中。一个方法包括处理每个投影图像，以减小除圆形物体(即，标志点的形状)之外所有物体的强度，识别每个图像中的局部极大值，生长局部极大值以生成候选区域，并对每个候选区域应用圆度和投影连续性标准。该方法以与其他方法相比将总图像处理时间减小至少一个数量级的方式，使用一系列应用的简单均值滤波器、差值以及除法/减法操作，来减小除圆形物体以外所有物体的强度。一个方法包括获得图像、应用2D高通滤波器以生成背景抑制后的图像、以及从所述背景图像中提取标志关联点。

实施例还提供用于对图像中的标志点位置进行准确的子像素定位的系统和方法。一个方法包括在两个维度中的每个上对图像中所识别的每个标志点的位置独立地进行评价。该评价通过对标志点的一侧或两侧采样，从而确定每个标志点的背景轮廓来进行。采用基线轮廓和标志点轮廓之比来修正背景影响，同时解决X射线衰减的非线性性质。残留的背景偏移通过减去标志点轮廓的拖影的平均值来去除。最后，对所得标志点轮廓进行评价，以确定其一阶矩(moment)，该一阶矩表示标志点的定位。替代地，可以通过使用位置拟合参数，将理论标志点轮廓拟合到该标志点轮廓来处理轮廓。该方法以明显比投影图像像素大小的分辨率高的分辨率，确定在计算机断层摄影扫描的投影图像中的标志点的位置。此外，该方法在存在明显不均质的背景和相邻结构时执行标志识别或定位。一个方法包括获得包含标志点的图像、得出背景修正后的标志点轮廓、得出背景和基线修正后的标志点轮廓、以及定义标志点的子像素位置以作为背景修正后的标志点轮廓的尖峰的中心。

实施例也提供用于快速地识别并分配近似的3D位置给与图像中的标志点相关联的标志的系统和方法。一个方法包括获取投影图像之间的标志点的位置变化，并将位置变化拟合到通过使用Z和R柱坐标的不同的可能值所定义的理论投影轨迹。该步骤包括为柱坐标维度中的一个或多个来定义一组可能值，并从该组中选择这些值中的一个。该方法然后在每个投影图像中确定落在与由从所述组中选择的值所定义的轨迹相距一距离范围内的标志点。然后通过找出每个标志点距轨迹范围的边界的均方距离来给落在该范围内的每个标志点分配度量。对于落在该范围内的所有标志点对度量求和，以得出与特定实验坐标值集合相关联的单个度量。对可能坐标值的集合中的每个值重复这些步骤，并评价度量的所得值(作为坐标值的函数)以确定函数尖峰，这些函数尖峰均假定为表示一个或多个标志的群集(cluster)。将与尖峰相关联的一个或多个柱坐标的这些坐标值分配给对应的得出的物理标志。对其余的柱坐标重复这些步骤，直到给标志分配了所有三个坐标为止。因为度量拟合的性质，所以在存在噪音和患者运动时可以找出平均标志定位，并且提取并分开位置彼此靠近的多个标志。一个方法包括获得图像序列，各个图像表示扫描器的一个转动角度，并具有第一值和第二值。该方法还包括分析贯穿该图像序列的标志点位置的第一值的性能，以确定标志的径向距离以及标志的角度位置。该方法还包括分析贯穿该图像序列的标志点位置的第二值的性能，以确定标志的轴向位置。

实施例也提供用于通过称为标志跟踪的处理而将每个投影图像中的每个标志点分配给最可能与该标志点相关联的单个物理标志的系统和方法。一个方法包括获得带有多个标志点的图像序列，并计算每个图像中的物理标志的预期U和V位置及范围。该方法还包括选择物理标志，并生成候选点列表，所述候选点列表由落在为该物理标志从理论上确定的每个投影图像中的预期标志点定位的范围内的标志点所构成。该方法进一步包括，通过将该列表中每个标志点的实际图像定位减去物理标志的理论轨迹，而将与理论的误差分配给每个标志点，从而来处理该列表。该方法然后逐个投影地遍历(traverse)该列表，并且基于U和V方向上的候选列表中的连续点的误差之间的一阶时间导数以及其他标准，将该列表缩窄到每投影图像零或一个“优良点”。

通过对具体说明书以及附图加以考虑，本发明的其他方面将会变得明显。

附图说明

图1示意性示出锥形束计算机断层摄影系统。

图2示意性示出图1的计算机。

图3示出包括大量标志点的图像。

图4是示出标志点识别方法的流程图。

图5是示出图4的标志点识别方法的背景减小过程的流程图。

图5A是示出1D累积均值滤波器的应用的流程图。

图6是作为图5的背景减小过程的一部分的、用以在一个维度上对图像进行快速均值滤波的过程的图解说明。

图7示出在应用了图5的背景减小过程之后的图3的图像。

图8是示出图4的标志点识别方法的标志点提取过程的流程图。

图9和图10示出在执行了图8的标志点提取过程的步骤之后的图7的图像。

图11是示出标志点子像素定位方法的流程图。

图12是示出图11的标志点子像素定位方法的标志点轮廓得出过程的流程图。

图13是在图12的标志点轮廓得出过程中进行分析的围绕标志点的示例性像素列的图解说明。

图14是示出在图12的标志点轮廓得出过程中所使用的标志点轮廓和基线轮廓的示例性值的曲线图。

图15是示出图11的标志点子像素定位方法的标志点位置确定过程的流程图。

图16是示出使用图15的积分(integral)标志中心确定过程的示例性标志点轮廓曲线的子像素标志点中心位置的图示。

图17是示出使用图15的曲线拟合标志点中心确定过程，对标志点轮廓曲线的示例性最佳拟合高斯曲线的子像素标志点中心位置的图示。

图18是示出物理标志三维定位方法的流程图。

图19A-B是示出对标志点在U维度上的性能(behavior)进行分析的过程的流程图。

图19C-D是示出对标志点在V维度上的性能进行分析的过程的流程图。

图20是示出在U方向上存在运动时使的试验值拟合到通过投影图像序列的大量标志点的多个U位置的过程的示例的曲线图。

图21是示出R取固定值作为函数的试验R和值对的累积拟合度量的示例性一维绘图的曲线图。

图22是示出使Z的试验值拟合到贯穿投影图像序列的三个标志的V轨迹的过程的示例的曲线图。

图23是示出作为函数Z的累积拟合度量的示例性一维绘图的曲线图。

图24A-E是示出标志点到物理标志映射方法的流程图。

图25是示出关于贯穿投影图像序列而识别的所有标志点的示例性物理标志和示例性U位置的预期U范围的曲线图。

图26是示出具有在图25中所示U范围内的U位置的图25中所示标志点的曲线图。

图27是示例性候选标志点列表的图解说明。

图28是示例性候选标志点列表的变平过程的图解说明。

图29是示例性候选标志点列表的同步化过程的图解说明。

图30是当“下一个”投影图像包含两个或更多候选标志点时处理“下一个”投影图像的图解说明。

图31是当“下一个”投影图像包含单个候选标志点时处理“下一个”投影图像的图解说明。

图32是当“下一个”投影图像不包含候选标志点时处理“下一个”投影图像的图解说明。

图33是当优良点列表包含优良点和可疑点时处理“下一个”投影图像的图解说明。

具体实施方式

在详细解释本发明的任何实施例之前，要理解的是，本发明在其应用时不限于在下面的描述中所阐述的或者在附图中示出的构造的细节以及部件的布置。本发明能够具有其他实施例，并且能够以各种方式实践或实施。

还应予以注意的是，可利用多个基于硬件和软件的装置、以及多个不同结构的部件来实现本发明。此外，如随后的篇章中所描述的，图中所示特定配置旨在例示本发明的实施例。替代配置也是可能的。

X射线是一种形式的射线，其类似于光线，但可以穿透人体的射线。然而，当X射线穿过人体时，X射线的强度减弱。强度的减弱量与X射线穿过的组织的密度有关。X射线的这一特性被用以产生人体以及其他物体内部的结构的图像。

X射线系统传统上包括产生X射线的发散束的X射线管。人或物体置于X射线管与一张膜之间，这张膜对在每个定位处射在其上的X射线的强度敏感。开关X射线管，对该膜进行显影，并且所得图像呈现X射线路径中的结构的“阴影(shadow)图像”，且其中，结构越密集则在图像中显得越白。

现代X射线过程使用一个或多个检测器，而非使用膜。检测器以电子方式检测并量化到达检测器的X射线的量(即，密度)。使用检测器，建立X射线计算机断层摄影(“CT”)系统，其使X射线源绕着患者转动，并且利用在患者的与X射线源相反的一侧的单个宽条状(single-wide strip)检测器元件，以电子方式检测所得X射线。针对所有不同的X射线位置收集来自检测器的数据，然后在已知为重建的过程中对其进行组合。组合后的图像表示穿过患者的单个有限宽度“切片”，其中每个点处的图像密度表示特定物理定位处的组织的X射线密度。X射线源、检测器、以及允许X射线源转动的机械结构的组合已知为“CT台架”。

通过在图像采集或扫描之间移动患者、X射线源、或者两者的同时重复上述过程，可以获取多个切片。例如，移动支承患者的平台且同时也转动X射线源，代替数据切片而产生“螺线(helix)”。另外，将检测器条或环的大小或宽度从单行检测器增加到多行(例如增加到256行)，允许对更多数据进行更加快速的采集。此外，用较大的二维检测器替换检测器条，获取每个X射线源位置处的整个检测器面板图像，而非仅单条数据。对在数量上达到600或者更多的这些的图像的收集已知为投影图像。每个投影图像表示随着X射线源和检测器绕着患者同步地转动从不同的透视或角度对患者的X射线快照。因为锥形X射线束需要覆盖二维检测器，所以这种CT成像已知为锥形束(“CB”)CT成像。图1示意性示出了根据本发明的一个实施例的CB CT成像系统10。

CB CT成像系统10包括扫描器12和计算机14。扫描器12包括CT台架13，其包括X射线源16、检测器18、和转动载台20。X射线源16和检测器18在转动载台20上在彼此对面对准，该转动载台20使X射线源16和检测器18绕着患者22移动。患者22由座位24所支承。图1中所示的成像系统10为牙科成像系统。因此，患者22坐在座位24中并将他或她的下颚放置到托架(rest)26中。当台架13转动以完成对患者的头部的扫描时，托架26保持患者的头部相对静止。

扫描器12将图像数据输出至计算机14。图像数据表示在扫描期间由检测器18所检测到的X射线的强度水平。计算机14连接至控制台30，其包括显示器32和一个或多个输入和/或输出装置，例如键盘34。使用者使用控制台30以与计算机14交互。例如，用户可以使用控制台30以从计算机14请求图像或其他数据。计算机14将所请求的信息提供至控制台30，并且控制台30将该信息显示在显示器32上、将该信息打印至打印机(未示出)，和/或将该信息保存至计算机可读存储器模块(未示出)。

图2示意性示出了根据本发明的一个实施例的图1的计算机14。计算机14包括输入/输出接口40、电子处理单元(“EPU”)42、以及一个或多个存储器模块，诸如随机访问存储器(“RAM”)模块44和只读存储器(“ROM”)模块46。输入/输出接口40从扫描器12接收图像数据，并将该图像数据提供给EPU 42。在一些实施例中，输入/输出接口40将图像数据储存至存储器模块，诸如RAM模块44，并且EPU 42从存储器模块获得图像数据用于处理。输入/输出接口40也可将数据传输至扫描器12，诸如校准和操作数据或指令。

EPU 42接收图像数据，并通过执行一个或多个应用或模块来处理信息。在一些实施例中，应用或模块储存在存储器中，诸如ROM模块46。在其他实施例中，应用或模块可储存在硬盘存储单元上并被传送至RAM用于执行。如图2所示，ROM模块46储存标志点识别模块50、标志点子像素定位模块52、物理标志三维(“3D”)定位模块54、标志点映射模块56、患者运动识别模块58、和图像重建模块59。EPU42检索并执行标志点识别模块50，以识别所接收的图像数据中的标志点(参见图4-10)。EPU 42检索并执行标志点子像素定位模块52，以提供对标志点位置的准确的子像素定位(参见图11-17)。EPU 42检索并执行物理标志3D定位模块54，以确定的通过物理标志在投影图像中的标志点识别的所述物理标志的三维位置(参见图18-23)。EPU42检索并执行标志点映射模块56，以将在投影图像中识别的标志点映射到特定的物理标志(参见图24-33)。EPU 42检索并执行患者运动识别模块58，以得出患者运动信息。EPU 42检索并执行图像重建模块59，以将投影图像与运动信息结合并生成一组运动修正轴图像。

完成台架旋转需要大约8到40秒。在此期间，患者可能移动，导致所得的图像的模糊。典型的图像分辨率为0.25毫米的量级。因此，这个量级的患者移动通常导致图像模糊，并且大范围的患者移动可以使得所得图像对于其所期望的临床目的而言变得不可接受。

计算机14(参见图12)可以通过跟踪所得图像中刚性物体的移动来针对患者移动对图像进行修正。例如，在其中没有患者移动的理想条件下，被成像的刚性物体随着台架绕着患者转动而以明确定义的方式改变在投影图像的二维上的位置。患者移动导致图像中物体的期望位置之间的偏离。通过测量明确定义的物体与其期望位置的偏离，可以测量并修正患者运动量。具体而言，如果图像中有至少三个物体，则可以组合这些物体与其期望位置的所测量的偏离来确定六参数患者位置误差值(例如，患者运动矢量)，其可以应用于图像以修正患者移动。

为了确保在图像中存在所期望的数量的明确定义的刚性物体，可以将一个或多个基准标志置于患者上。这些标志一般包括铅或钢BB，其较为致密并可防止或限制X射线透过。这些标志也可以由其他材料制成，并且构造成其他的形式或形状，其在扫描期间生成的投影图像中以高比例可见。每个标志或多个标志可位于多层粘合剂之间，并且该粘合剂可以涂覆于患者，以确保标志在手术期间不移动。

将标志置于患者上，使得由CB CT成像系统所生成的每个视野或图像均包括至少三个标志点。例如，在患者上可放置七到九个标志，以确保至少三个标志点在由CB CT成像系统所生成的每个图像内，以减少位置测量噪声，并允许结果的统计组合。在一些实施例中，标志均匀地分隔，以相对于患者具有最大的空间分布，并避免对于正常的图像判读(interpretation)的干扰。

识别图像中的标志

在将标志置于患者上之后，对患者进行扫描，诸如，对患者的头部进行扫描，并将表示投影图像序列的所得图像数据传输至计算机14(参见图2)。图3示出了包括八个标志点102的示例性投影图像100。作为针对患者移动而修正所生成的投影图像时的初始步骤，计算机14处理所生成的投影图像，以识别每个图像中的标志点。图4示出了根据本发明的一个实施例的标志点识别方法60。当EPU 42执行标志点识别模块50时，由计算机14的EPU 42来执行标志点识别方法60。

如图4所示，方法60的第一步通常包括减小每个投影图像的背景的强度(步骤62)。此背景抑制步骤使图像的背景部分最小化。图像的背景部分包括除了标志点之外的图像中的所有事物。因为基准标志具有明确定义的形状(即，球形或圆形)，所以步骤62使图像中未具有此明确定义的形状的所有事物最小化。

在一些实施例中，通过对每个图像应用线性高通滤波器来实现背景抑制。为了获得可接受的速度，执行三个优化。首先，通过将原始图像与低通滤波后的图像组合来实现高通操作。其次，代替在两个图像方向上同时执行滤波器，而在两个或更多个维度上(即，以0°、45°、90°和135°)对图像依次应用一维高通操作。最后，利用通过下述的累积像素方法所实施的简单均值滤波器来实现低通滤波器操作。

一维高通滤波器操作将趋于加强在滤波器的方向上快速变化的区域。在不同方向上执行的一系列这类操作将趋于相对于所有其他形状而加强圆形物体。由于标志所产生的标志点为圆形，所以这一系列滤波将加强标志点，且减弱任何周围的结构。

图5示出了背景减小步骤62中所涉及的过程或步骤。如图5所示，通过首先对原始图像A进行简单的一维的、累积均值低通图像滤波(步骤69)，以生成滤波后的图像B来执行多维高通滤波操作。然后，如下所述，将滤波后的图像B与原始图像A组合以生成图像C(步骤70)。使用来自于先前的滤波步骤的输出图像作为对下一个滤波步骤的输入，对每个滤波方向重复这两个步骤(步骤71-73)。图5所示过程允许快速的、多维高通滤波，以减小图像的背景。图6示出了一个方向上的一维低通滤波中所涉及的这些步骤的示例。

代替使用更典型的基于内核的滤波方法，通过使用累积均值滤波器(图5A中所绘)能更加快速地实现步骤69和71的滤波。此快速滤波通过生成新图像(图像B1)来实现(步骤74)，其中，新图像中的各个像素均被设定为在滤波器的方向上从该特定像素向下直至第一像素的、来自原始图像A的所有像素值的累积和。例如，如图6所示，当计算机14一开始从扫描器12获得原始图像(图像A)时，该方法考虑包括七个像素的图像A的一行(其中，处理方向为水平的)，然后通过将该行中的各个像素设定为该像素左侧的像素值们加上该像素的值之和来生成累积图像(图像B1)。

一经在步骤74生成累积图像，图像B1就移位与低通滤波因子相等的量(例如，2)，以生成移位后的图像(图像C1)(步骤75)。通过对移位后的边缘像素简单地调零来处理图像边缘。例如，如图6所示，累积图像(图像B1)中的每个像素向左移位两个位置，以生成移位后的图像(图像C1)。因为前两个像素(即，边缘像素)的左侧两个位置没有像素，所以这些像素设定为0。

然后将累积图像(图像B1)减去移位后的图像(图像C1)以生成差值图像(图像D1)(步骤76)。如图6所示，为了减去图像，将未移位的累积图像(图像B1)中的各个像素值减去移位后的图像(图像C1)中的对应像素值。最后的步骤是通过将差值图像中的每个像素除以滤波器因子(移位值)以归一化差值图像(图像D1)，从而产生低通滤波后的图像B(步骤77)。

步骤74-77的计算时间不依赖于低通滤波器的内核大小。这些步骤中的计算粗略涉及2×N个加法和减法，而不是(内核大小)×N个操作，其一般需要普通的基于内核的滤波器。对于512×512个像素的图像和典型的平滑因子15，这表示计算时间相对于传统的1-D内核方法实现了7.5:1的减少。与2-D内核相比，(相对于两个1D累积均值滤波器操作)对于滤波器因子15，速度提高为56:1。

在生成一维低通滤波后图像之后，通过将低通图像(图像B)与用于生成其的图像(图像A)组合而生成高通滤波后图像(图像C)(步骤73)。这可通过两个图像之间的简单差值操作来进行。然而，除法操作用于第一方向滤波步骤70中的处理。使用除法以补偿X射线衰减的指数性质和归一化图像(图像B)，以生成图像C(步骤70)。对差值图像执行除法是有效的，这是因为X射线衰减是以材料密度的负指数的形式出现的。因此，由于除法等同于对数相减然后求指数，所以对差值图像执行除法避免了计算开销高的取对数和求指数操作。

在顺序的维度中的每个上重复应用一维的、累积均值低通滤波器，并组合输入和输出图像的过程(步骤71-73)。然而，在处理第一维度之后，在步骤72可以使用减法代替除法，这是因为对X射线衰减的指数性质的修正通过在处理第一维度时执行除法(步骤70)来完成。一般来讲，0、45、90和135度处的四个通道(pass)足以将球状标志点从宽广范围的背景中适当地隔离。图7示出了在应用图5的背景减小方法之后的图3的图像100。如图7所示，背景的大部分(即，图像中除标志点102外的部分)已被减小或暗化，允许标志点102更容易识别。

回到图4，一旦背景通过上述步骤减小，则从高通滤波后的图像中提取标志点(步骤80)。图8示出了从高通图像中提取标志点的过程或步骤。如图8所示，该过程开始于在每个高通图像中定义局部极大值(步骤82)。此步骤包括识别每个图像中的每个其所有8个相邻像素都具有比其小的值的像素。接着，使用所识别的局部极大值作为种子点以生长区域(步骤84)。此过程将特定区域内所包含的局部极大值组合为单个候选物体或区域。在一些实施例中，所产生的区域的大小或范围受根据经验指定的下阈值的限制。图9示出了在区域生成局部极大值(步骤82和84)之后的图7的图像100，其生成候选区域104。

在步骤86，对每个候选区域应用形状标准，以识别哪个区域最有可能代表标志点。在一些实施例中，此过程通过取在(1)候选区域(a)的像素面积与(2)通过候选区域的平均半径的平方乘以π((a的平均半径)²×π)限定的面积之间的比率来执行。该比值用于确定候选区域是否具有合适的形状以作为标志点(例如，候选区域是否具有圆形形状)。

最后，在步骤88，通过比较相邻投影图像中的候选区域的定位并应用接近度(nearness)标准来进一步处理候选区域，以丢弃不可能的候选者。例如，如果候选区域在顺序的投影图像之间其位置移动太多，则该候选区域不可能代表标志点。类似地，如果两个候选区域太靠近在一起，则这两个区域不可能代表标志点。图10示出了在对候选区域104应用形状标准并丢弃不可能的候选区域104(步骤86和88)之后的图9的图像100。如图10所示，在这些步骤之后，仅留下与八个标志点102(参见图7和8)相对应的候选区域104。

在一些实施例中，对图像应用二维高通滤波器包括首先对该图像应用二维低通滤波器以生成低通图像。之后，将原始图像与低通图像组合。

在一些实施例中，标志点识别方法60可以在时间上以每图像少于50毫秒的量级来执行。另外，当使用球形标志时，标志点识别方法60表现出对标志的高选择率，以及对非标志相关联结构的高丢弃率。例如，在一些实施例中，标志点识别方法60可以具有100％的灵敏度，并且在实际情形中具有比40:1更好的特异性(低于2.5％的错误检测)。

确定标志的子像素中心点

即使在使用标志点识别方法60在图像中识别标志点之后，在这些标志点可以用于执行运动修正之前也需要对其进一步地或者更准确地限定。具体地，必需以比由标志点识别方法60所提供的像素分辨率(子像素)好得多的分辨率来知晓标志点的有效二维位置。可能需要小到像素分辨率的1/8的分辨率。再者，这必须在图像内以及图像之间都存在极大变化的背景中完成。图11示出了根据本发明一个实施例的标志点子像素定位方法110。当计算机14的EPU 42执行标志点子像素定位模块52时，由EPU 42执行标志点子像素定位方法110。在一些实施例中，标志点子像素定位方法110如下所述使用由标志点识别方法60所识别的标志点作为起始点，来确定准确的子像素标志定位。

如图11所示，标志点子像素定位方法110通常包括标志点轮廓(profile)得出过程(步骤112)和标志点中心确定过程(步骤114)，对投影图像中的两个维度中的每个重复执行这两个步骤(步骤116)。如以下更详细描述的(参见图12-14)，标志点轮廓得出过程112应用一系列步骤，以去除背景中的大部分，并在每个期望的维度上生成多像素标志点轮廓。在标志点中心确定过程114中，将所生成的标志点轮廓与标志的预期形状相比较，并得出一维子像素标志点中心定位(参见图15-16)。

图12示出了根据本发明一个实施例的标志点轮廓得出过程112。考虑得出标志点的垂直位置的操作，其中，垂直方向上的像素序列称为列，且水平方向上的像素序列称为行。用于得出标志点的水平位置的过程是相似的，但所有对行的提及都由列替换，而对列的提及都由行替换。像素的行和像素的列均可称为像素线(pixel line)。

该过程的第一步涉及通过在围绕标志的背景列的范围内对像素值求平均值来得出垂直背景轮廓。图13示出了围绕标志点122并以其为中心的像素的示例性列121。步骤120对围绕标志点122的背景列121a的像素值求平均值，并排除以标志点122为中心的标志列121b的像素值。对背景列121a求平均值包括对跨列121a的每一行的像素值求和并将和除以列数。对背景列121a求平均值生成单列背景值，其为标志点122的背景轮廓。

接着，通过对标志列121b求平均值来得出标志点轮廓，其中，所述标志列121b是以标志点122为中心的列(步骤124)。再者，对标志列121b求平均值包括对跨列121b的每行的像素值求和并将和除以列数。对标志列121b求平均值生成单列标志值，其为标志点122的标志点轮廓。

在步骤126，以逐行为基础，将背景轮廓除以标志点轮廓。例如，将组成背景轮廓的平均值的列中的每一行中的值被除以组成标志点轮廓的平均值的列中的对应的行的值。按此顺序对行值执行除法确保出现图像反转(例如，暗区变亮而亮区变暗)。另外，步骤126使用除法(例如，代替减法)以补偿X射线衰减的负指数性质(这是因为，如所指出的，除法在数学上等同于对两个轮廓的对数之差求指数)。

作为替代实施例，像素值的对数可以在找出标志点轮廓和背景轮廓之前计算。然后，通过以逐个元素为基础将两个轮廓相减并对每个元素求指数来实现对标志点轮廓的背景补偿。此方法较之于以上方法略微精确一些，这是因为其能更准确地处理在求平均值的行上的图像变化，但其比简单的求平均值和除法的方法更耗时间。

此时，已获得标志点的合理的基线修正后的轮廓。然而，该轮廓的拖影(tail)与零相比仍会存在一些残留的偏移，如果拖影是非对称的则尤其如此。该轮廓通过下列程序而达到零。平均基线偏移值通过对离开标志点的尖峰(例如，具有最高值的标志点内的像素)的拖影中的列上的像素值求平均值来确定，例如在其中像素值不再随着其从标志点中心移动离开而减小的列上对像素值求平均值来确定(步骤128)。关注于从此平均像素中排除与和此标志点相邻的任何结构，例如图像中的其他标志或其他结构相关联的拖影。图14用曲线图示出了示例性标志点轮廓129和示例性基线轮廓130的值，以及对标志点轮廓133进行完全基线修正之后的结果。

回到图11，在生成基线修正后的标志点轮廓之后，定位方法110的第二步包括标志点位置确定过程114。图15更详细地示出了此过程。如图15所示，存在至少两个替代项或选项以用于执行此过程。第一选项包括确定标志点轮廓的积分(步骤140)，以及将标志点位置定义为表示该积分的一半的沿着积分的位置(步骤142)。此标志点位置有效地为一阶矩的积分，或者是在其中积分在该点左侧的部分等于积分在该点右侧的部分的位置。图16示出了示例性标志点轮廓曲线145的示例性子像素标志点位置定位144。

回到图15，用于确定标志点位置的第二选项使用最小二乘法或其他拟合将适当的曲线函数拟合到标志点轮廓(步骤146)，并且将标志点位置定义为表示曲线函数的中心的参数(步骤148)。图17示出了拟合到示例性标志点轮廓曲线151的最佳拟合高斯曲线150的示例性子像素标志点位置定位149。

步骤146中所使用的曲线函数可以是标志点的高斯曲线或预期曲线，该曲线将与成像物理学有关的问题，例如部分体积效应(partialvolume effect)和指数衰减特性纳入考虑。通常，步骤146所使用的曲线函数可以是其中曲线定义参数之一表示曲线中心的任何曲线。曲线函数也可以处理非常量的基线轮廓。在一些实施例中，如果选择处理非常量基线轮廓的曲线函数，则标志点轮廓得出步骤112中的步骤128和132可以去除。此外，可以基于可测量的图像特性，例如已知的标志维度、像素大小、以及其他维度上的标志的近似子像素定位，来针对每个标志、图像、或者图像组而动态地修改曲线函数。

确定标志的3D柱坐标

即使在图像中识别标志点和标志点中心(例如，使用标志点识别方法60和标志点子像素定位方法110)之后，标志点的积累也需要分类成使得每个标志点均与负责生成该点的物理标志适当地相关联。为了做到这一点，第一步是唯一地识别每个物理标志点。物理标志点均通过其三维(3D)坐标而唯一识别。因此，将标志点与物理标志相关联的第一步是识别物理标志并且给每个分配其3D坐标。然而，此过程由于患者移动以及在投影图像中确定标志点位置的误差而会复杂化。另外，由于存在多个标志，所以在一些投影图像中，一些标志点会重叠在其他标志点之上。

图18示出了用于识别物理标志中的每个并在3D空间中给每个物理标志分配位置或定位的物理标志3D定位方法160。当计算机14的EPU 42执行物理标志3D定位模块54时，由EPU 42来执行物理标志3D定位方法160。物理标志3D定位方法160通常贯穿(through)随着CT台架转动而生成的图像序列来跟踪由标志产生的点的位置，以将标志映射到3D空间中的位置。

现有的跟踪技术在多个被检部位之间使用多个传感器和多个固定关系来跟踪患者在3D空间中的位置。然而，方法160中所应用的跟踪技术或方法使用隔离的单点检测部位。因此，代替比较同时收集的多个部位的特性，该跟踪方法收集单个点随时间的特性，并将观察到的随时间的位置变化匹配到3D位置，其将产生相同的随时间的位置变化。此跟踪方法除CT成像台架外无需额外的设备，并且容许相对较大量的图像噪音和患者移动。

如图18所示，方法160开始于得到每个标志点在每个投影图像上在标为“U”和“V”的两个图像维度上的位置(例如，通过标志子像素定位方法110而确定的每个标志点的中心)(步骤162)。U是与CT台架的转动圆形相切的维度，且V是与CT台架的转动轴线相平行的维度。标志点的这些维度位置可以通过标志点识别方法60、标志点子像素定位方法110、或其组合来确定。

如果标志在3D空间中具有固定的位置，则其在每个投影图像上的U和V位置(或值)都可由下列公式表示：

$U=\frac{\mathrm{DSD}*R*\sin (\mathrm{\θ}-\mathrm{\φ})}{\mathrm{DSO}-R*\cos (\mathrm{\θ}-\mathrm{\φ})}$

$V=\frac{\mathrm{DSD}*Z}{\mathrm{DSO}-R*\cos (\mathrm{\θ}-\mathrm{\φ})}$

其中，3D柱坐标(R＝径向距离，且Z＝高度或轴线位置)表示标志的3D位置，θ是CT台架自其起始位置的转动角度，DSD是从X射线源到检测器的距离，且DSO是从X射线源到物体或患者的距离。

为了确定标志的3D位置，必须使用上述U和V公式以及标志的已知U和V值来确定标志的柱坐标或参数(R、和Z)。然而，特定投影图像中的U和V位置由于测量误差而会与它们的正确值不同。而且，R、和Z由于患者运动而成为时间的函数(例如，R(t)、和Z(t))。

由于这些困难，方法160为R、和Z找出“平均”值(步骤164)。这些平均值是产生每个投影图像中标志的实际位置与如上述U和V公式所指定的每个投影图像中标志的预期位置的最小和平方误差(minimum sum squared error)的值。通常，下列误差公式提供这些平均值：

$\mathrm{Err}=\underset{\mathrm{\θ}=1}{\overset{n}{\mathrm{\Σ}}}({[U_{\mathrm{meas}}-\frac{\mathrm{DSD}*R_{\mathrm{avg}}*\sin (\mathrm{\θ}-{\mathrm{\φ}}_{\mathrm{avg}})}{\mathrm{DSO}-R_{\mathrm{avg}}*\cos (\mathrm{\θ}-{\mathrm{\φ}}_{\mathrm{avg}})}]}^2+{[V_{\mathrm{meas}}-\frac{\mathrm{DSD}*Z_{\mathrm{avg}}}{\mathrm{DSO}-R_{\mathrm{avg}}*\cos (\mathrm{\θ}-{\mathrm{\φ}}_{\mathrm{avg}})}]}^2)$

然而，为R_avg、Z_avg和找出期望值的过程由于U维度不依赖于Z而简化。这样，该过程可以分成两个过程或步骤：(1)分析标志在U维度上的性能(参见图19A-B中的方法166)和(2)分析标志在V维度上的性能(参见图19C-D中的方法168)。

图19A-B示出了根据本发明一个实施例的分析标志在U维度上的性能的方法166。如图19A所示，方法166开始于选择表示R和参数的所有可能值对的子采样的R和值对的数集(步骤170)。子采样集合的大小表现为在R、的可获得准确度与执行速度之间的折衷。通常，对R和值加以选择，以对物理上可能的R和值的范围进行均一采样。即是说，R在0和最大扫描器半径视场之间，且在0°和360°之间。替代地，可以假定R的值并且仅选择的可能值的子采样。因此，尽管参照R和值的对来描述方法166的其余步骤，但应予以理解的是，这些步骤也可以用R的假定值和的变化值来执行。

如图19A所示，在选择子采样集合之后，从该子采样集合选择R和值对以进行处理(步骤172)。接着，使用所选的R和值对以及上述U公式，为每个投影图像确定“理论U”值(步骤174)。因为每个投影图像均与CT台架的特定转动角度相关联，所以在U公式中使用与每个投影图像相关联的θ值来确定每个“理论U”。然后围绕每个“理论U”值定义U范围(步骤176)。U范围以“理论U”值为中心，并且在每个方向上距中心均具有如下的幅度，该幅度等于由患者运动而导致的标志与其平均或者预期U位置的最大的预期U变量。

接着，将每个投影图像中识别的标志点的U位置与每个投影图像的U范围相比较。在步骤178中，选择识别的标志点。如果所选择的投影图像中识别的标志点具有落在投影图像的U范围内的U值(如在步骤180中确定的)，则将该标志点识别为标志点候选者并且方法推进到步骤182。忽略落在该范围之外的标志点。应予以理解的是，多于一个的标志点可具有在特定投影图像的U范围内的U值。实际上，如果存在位于患者上的多个相邻的物理标志，则投影图像中所识别的多于一个的标志将可能具有在特定投影图像的U范围内的U值(即，θ的特定值)。

图20示出了贯穿300个投影图像的序列的大量标志的U位置。如图20所示，执行值的三个试验或变化，而同时假定了对应的R值并且使其保持不变。具体地，第一个试验使用了60度的值，第二个试验使用了100度的值，且第三个试验使用了140度的值。这些试验中每个的U范围都用虚线示出在曲线图中。对于每个投影图像(参见x轴)，具有U范围内的U位置的每个标志点将被标记为候选标志。

接着，为每个投影图像中的每个候选标志点确定拟合度量(fittingmetric)(步骤182)。在一个实施例中，拟合度量为标志候选者的U值与U范围的较近边界之间的距离的平方。例如，如果标志具有20的U值，且U范围为0到30，则拟合度量将设定为10(其为20与较近的U范围边界30之间的距离)的平方，或100。如果标志的U值与“理论U”值相同，则拟合度量最大化。然后对来自所有候选点的拟合度量值求和，以生成累积拟合度量，该累积拟合度量被分配给当前选择的R和值对(步骤184)。如图19A的步骤186中所示，对所有投影图像中的所有点重复步骤178-184。

如图19A所示，对子采样集合中的每个R和值对重复步骤172-186(步骤188)，以将累积拟合度量分配给每个值对。在给所有R和值对都分配了累积拟合度量之后，方法166评价作为R和的函数的、用于所有R和值对的累积拟合度量(图19B中的步骤190)。方法166使用尖峰检测，以确定此函数中的R和的尖峰(步骤192)。函数的尖峰是尖峰R和值对，其每个均表示一个或多个物理标志。方法166将每个尖峰R和值对分配给与该尖峰相关联的对应的标志(步骤194)。图21示出了作为的函数(假定R的值并使其保持不变)的用于R和值对的累积拟合度量的一维绘图(例如，直方图)的示例。如图21所示，绘制的函数的尖峰表示物理标志的值。

在分析了标志点在U维度上的性能之后，尽管候选标志已定位并分配了R和值，但对于相同的R和值可能存在多于一个的标志。为了应对这种情形，对标志点在V维度上的性能进行分析。图19C-D示出了根据本发明一个实施例的分析标志点在V维度上的性能的方法168。如图19C所示，方法168开始于利用在方法166期间分配给该标志候选者的尖峰R和值对(参见图19B中的步骤194)。其使用这些值来识别候选点，以进一步在Z方面进行评价(步骤200)。然后选择Z的所有可能值的子采样(步骤202)。对子采样中的值的数量加以选择，作为在执行速度和期望的Z值分辨率之间的折衷。对子采样的范围加以选择，以覆盖Z的物理上可能的值(可以处于扫描视场内某些点处的点)的范围。

然后从子采样集合选择Z值(步骤204)。使用先前为这一个或多个标志确定的R和值对、所选择的Z值以及上述的V公式，为每个投影图像确定“理论V”值(步骤206)。由于每个投影图像均与CT台架的特定转动角度相关联，所以在V公式中使用与每个投影图像相关联的θ值以确定每个理论V。然后围绕每个理论V值定义V范围(步骤208)。V范围以“理论V”值为中心，并且在每个方向上距中心均具有以下的幅度，该幅度等于由患者运动而导致的标志与其平均或者预期V位置的最大的预期V变量。

接着，作为候选者列表的一部分的每个标志点的V位置将其V值与此Z试验的V范围比较。在步骤210中，选择投影图像中的点。如果该点具有落在投影图像的V范围内的V值(步骤212)，则标志可识别为进一步限制的标志点候选者。为了阐释清楚，点必须落在与其R和值相关联的U值的范围内，以及落在与其R、和Z值相关联的V值范围内，以便被视为候选者。如上所述，在特定投影图像中所识别的多于一个的标志点可具有在与该投影图像相关联的当前V和U范围内的V和U值。如在图19C的步骤216中所示，对所有投影图像中的所有点重复步骤210-215。

图22示出了贯穿300个投影图像的序列的三个标志的V位置。图22中，示出了Z值的三个试验或变化。具体地，第一个试验使用了70毫米的Z值，第二个试验使用了30毫米的Z值，且第三个试验使用了-15毫米的Z值。这些试验中每个的V范围都用虚线示出在曲线图中。对于每个投影图像(参见x轴)，具有V范围内的V位置的每个标志点将被标记为标志候选。

回到图19C，为每个投影图像中的每个标志候选确定拟合度量(步骤214)。拟合度量为标志候选者的V值与V范围的较近边界之间的距离的平方。例如，如果标志具有20的V值，且V范围为0到30，则拟合度量将设定为10(其为10——20与较近的V范围边界30之间的距离——的平方)。然后对所有投影图像的度量拟合求和，以生成累积拟合度量，该累积拟合度量被分配给当前选择的Z值(步骤215)。

如图19C所示，对子采样集合中的每个Z值重复步骤204-218(步骤217)，以将累积拟合度量分配给每个值。在给子采样集合中的所有Z值对都分配了累积拟合度量之后，方法168评价Z和所得累积拟合度量之间的关系(图19D中的步骤218)。方法168然后使用尖峰检测，以确定该关系中的尖峰(步骤220)。

累积拟合度量-Z关系中的每个尖峰与唯一物理标志相关联。尖峰的Z值表示对应的物理标志的Z值。所分配的R和值是先前为一个或多个标志所找出的。这些标志中的每个目前具有用于所有三个柱坐标(R、Z)的已知值。图23示出了作为Z的函数的、累积拟合度量的一维绘图(例如，直方图)的示例。图23中，存在对应于三个单独的物理标志的三个尖峰，每个均在R和的相同值处。

对通过方法166所确定的其他R和对中的每个重复(步骤224)方法168(步骤200-222)，以隔离所有唯一物理标志并为每个确定Z值。方法168当所有标志候选者都已分离成唯一物理标志且每个物理标志都已分配了所有三个柱坐标时终止。所分配的柱坐标表示标志的平均R、和Z值，这些值表示标志的3D位置。

以上参照图19A-D所描述的方法166和168表示两步过程，该两步过程使用两个分离的方法将3D柱坐标分配给标志(即，第一方法分配R和值，且第二方法分配Z值)。标志3D定位方法160的第一替代实施方式将这些方法组合成单个操作，其中，选择R、和Z的值的子采样集合，并使用这些值的每个组合以基于组合后的U和V误差来确定单个拟合度量。该单个过程的结果是具有与每个隔离的物体标志相一致的尖峰值的三维阵列。然而，此实施方式较之于分离的方法实施会更耗时。具体地，如果R、和Z分别具有m、n和q个子采样，则第一替代实施方式相对于用于分离实施方式的(m*n)+q而言需要m*n*q次迭代。

也可使用标志3D定位方法160的第二替代实施方式来加速分开的实施方式。在第二替代实施方式中，按原来描述的两步过程找出R、和Z值，而第一方法166中处理的R值的数量(m)将通过以较粗略的频率对R子采样而减少。这由于拟合对R值不太敏感因而是合理的。然后，可以增加第三方法，该第三方法保持和Z固定并使用上述的拟合度量过程和R的较精细采样的子采样集合来优化R。

标志3D定位方法160的第三替代实施方式也可以使子采样对于一个参数、两个参数或者所有三个参数而言较粗略。此实施方式也可以包括用于每个标志候选者的最终拟合方法，其对覆盖较窄范围的每个参数子采样集合使用较精细的采样速率。此实施方式既可提高方法160的速度又可提高标志位置坐标的最终准确度。

对标志3D定位方法160的一个通常的挑战是处理其中标志的位置彼此非常靠近但对于每个标志仍允许合理的运动范围的情形。如果标志被分开少于最大预期运动，则方法160可能将两个标志标注为好似它们是一个标志，并给这一个标志分配在两个标志的实际位置之间的位置。为了防止这种情形，应当以足够的物理分隔将标志放置于患者上。但方法160需要处理其中没有遵循放置指导原则的情形。可以对方法160加以修改，以识别其中很大百分比的投影图像具有多于一个的在候选点接受范围内的标志的情形，或者识别其中累计度量尖峰大于单个标志的预期的情形。如果这些情形中任一个被检测到，则方法160可以包括另外的处理以更好地隔离这两个标志。一种方法将是，当交叠的标志可疑时，以较精细的子采样以及缩窄到仅覆盖该可疑标志的范围重复定位过程。

第二个挑战是在于当标志运动导致在累积拟合度量中出现两个尖峰时。在此情况下，单个标志将显现为两个。为了应对这一可能性，方法160可以将θ范围分成子范围，分离地在每个子范围内进行拟合，然后合并结果。另外，方法160可以适应性地改变用于三个坐标参数R、和Z的候选点接受范围。例如，方法160可以小接受范围开始，以确保紧密相邻的标志充分隔离，以及假定很少运动。如果在这些范围中未包括足够数量的点，则大概是由于运动，所以可以增大范围，直到找出足够的候选标志点为止。

通常，标志3D定位方法160定位标志，以在扫描过程中存在多达两厘米的患者运动时将分辨率调节到好于一毫米。方法160在存在24个或更多标志时也能有效地操作，并且容许测量噪音。最后，方法160可以用少于大约5秒的时间处理300个投影图像上的24个标志。

将点映射到患者上的物理标志

即使在识别并定位图像标志点，以及识别对应的物理标志并给其分配3D坐标之后，也必须将每个标志点分配给置于患者上的特定物理标志。在600个图像和九个标志的典型的CB CT扫描中，在投影图像中将会存在超过5000个被识别的点，其每个必须分配给九个物理标志之一。

就单个投影图像中所识别的标志点来看，没有(至少通常没有)足够的信息将在图像中所识别的特定标志点分配给位于患者上的特定物理标志。在理想情况下，如果这些标志在物理上是不动的，且其物理标志是确切知道的，则可以为投影图像中所识别的任何标志点定义从一个投影图像到连续图像的标志点位置变化的唯一的U和V轨迹。仅一个物理标志将与特定标志从一个投影图像到下一个构成的点位置轨迹相一致。因而，两个投影图像理论上足以将这两个标志点分配给其对应的物理标志。

然而，当标志在投影图像中的物理定位并非确切知道或者如果在扫描期间存在明显的患者运动量(或非预期的台架运动)的时候，会出现另外的挑战。在这些情况下，轨迹不确定性的程度增加，令投影图像中所识别的标志点对置于患者上的物理标志的适当分配复杂化。

图24A-E示出了根据本发明一个实施例的用于将投影图像中所识别的标志点分配给置于患者上的特定物理标志的物理标志映射方法250。通常，方法250在存在患者运动时跟踪连续投影图像中的标志点位置，并以逐个投影图像为基础评价每个标志点，以识别最可能与特定物理标志相关联的标志点。方法250使用大量不同标准，包括为标志点所测量的U和V以及与物理标志相关联的理论U和V之间的距离，以及这些距离的时间导数。

如图24A所示，方法250的初始步骤包括识别与扫描相关联的所有物理标志并得到每个的近似3D定位(步骤252)。这可以使用物理测量、来自重建的图像数据集合中的测量值、或者对投影图像的直接分析来进行。在一些实施例中，方法250使用上述物理标志定位方法160来获得用于每个物理标志的近似3D定位。

接着，方法250选择一个物理标志，并为每个投影图像计算预期的U位置(U_image)和预期的V位置(V_image)(步骤254)。方法250使用下列等式，这些等式描述预期在空间中的物理点是如何基于CT台架的特定几何形状而在图像之间移动的。

例如，假定CT台架绕着转动中心平滑转动，R是转动平面中一点与该中心的距离，是标志的角位置，且Z是标志与转动的中心平面(即，包括定义X射线源的定位的点的平面)的距离。DSD是从X射线源到检测器平面的距离，且DSO是从X射线源到CT台架的转动中心的距离。每个投影图像在台架转动角θ处获得。使用这些等式，方法250用每个投影图像角θ处的3D位置R、和Z来定义物理标志的预期U和V投影图像位置。

方法250然后使用预期U和V标志图像位置来定义用于物理标志的轨迹的预期U范围和预期V范围(步骤250)。这些范围足够大，以包括与该物理标志相关联的所有点，即使存在患者运动也如此。例如，如果预期患者在图像获取过程内移动20毫米之多，则用于特定物理标志的预期U和V范围将为：

和

方法250然后生成该物理标志的候选点列表(步骤258)，并添加在每个投影图像中所识别的、既落在预期U范围又落在预期V范围内的每个标志点(步骤260)。图25是示出了特定物理标志的示例性预期U范围(点划线之间)以及贯穿投影图像序列所识别的标志点的示例性U位置的曲线图。图26示出了图25中所示具有在图25所示预期U范围内的U位置的标志。这些标志中也具有在预期V范围内的V位置的任何标志都将添加到物理标志的候选点列表。例如，图27图示出了候选点列表的一部分。如图27所示，候选点列表包括来自每个投影图像的具有落在预期U和V范围内的U位置和V位置的每个标志。将从该候选点列表中选择与物理标志相关联的标志点的最终列表。如图27所示，候选点列表可以包括特定投影图像内的多于一个的标志(本示例中在投影1和投影6中都出现了两个候选点)。这出现在两个标志贯穿投影图像序列与彼此的路径相交时。如果两个物理标志在患者上布置为彼此靠近，则在特定投影图像的候选点列表中也可以存在多于一个的标志点。

如图24A和24B所示，在生成并填入(populate)候选点列表(步骤258和260)之后，方法250通过对于U和V维度，将候选点列表中所包含的每个点的实际标志点的定位减去物理标志的“预期”标志定位来处理候选点列表(步骤262)。此过程使标志的定位变平，以生成“与预期的误差”候选点列表。该“与预期的误差”候选点列表因而表示在投影图像中反映的物理标志的预期位置与实际标志点位置之间的差异(discrepancy)或误差的量(例如，由患者移动或测量误差而导致的)。图28图示出了应用于示例性候选点列表的此变平过程。

一经为物理标志生成“与预期的误差”候选点列表，就发生点跟踪过程。此过程的目的是从一个投影进行到下一个投影，并在下一个投影图像中选择最“可能”是正确点的一个点，即是说，该点最可能是由相关联的物理标志真实生成了的。应用大量的启发法(heuristic)以便实现这一点。整个过程通常具有两部分：同步化和跟踪。在同步化部分，对点列表进行处理，每次一个投影，直至发现特定点具有与物理标志相关联的高度可能性为止。跟踪部分始于该点并逐个投影继续，选择与预期由该标志形成的轨迹最为一致的一个(或零个)点(从该投影处的零个至多个点的列表中)。如果对于一定设置数量的图像没有找到一致的点，则同步丧失且将不得不重新获得。跟踪过程可以在投影编号递增的方向或投影编号递减的方向上进行。

在跟踪过程中，在每个新投影处，存在三种可能情况：不存在点，存在一个点，以及存在多于一个的点。运动和测量噪音使过程复杂化的原因在于，“正确点”可能不在其预期位置处。这通过实施“可疑点”的构思来解决。如果最佳候选点落在一些标准之外，则将其标为可疑。通过在将来的投影中的点位置的情境下查找该点的位置来去除此标记，并将该点包含或丢弃。以下提供此过程的具体细节。

在一些实施例中，方法250的跟踪部分，步骤266-352，可以在投影编号递增的方向或投影编号递减的方向上执行。在后一种情况下，应予以理解的是，对“下一个投影”和“前一个投影”以及与次序有关的其他术语的提及是相对于处理方向的。

回到图24B，在步骤262处理候选点列表之后，方法250生成物理标志的空的优良点列表(步骤264)。如下面将更详细描述的，用在投影图像中所识别的具有与物理标志相关联的高度可能性的标志点填充优良点列表。第一步是同步化(步骤266-270)。方法250检验候选者列表，每次一个投影(例如，从第一投影图像开始)，直至找到所定义的数量的顺序的投影图像，例如3个，其中每个图像具有候选点列表中的单个候选点(步骤266)。将与所检测的投影图像序列中的第一投影图像相关联的点标为“基准优良点”或基准点(步骤268)。另外，方法250将该点添加到优良点列表(步骤270)。图29图示出示例性候选点列表上的同步化过程。如图29所示，候选点列表包括三个顺序的投影图像的序列271，其中该序列271中的每个图像具有候选点列表中的仅一个点。在该同步化过程中，所检测的序列中的第一投影图像(标记为271a)包含单个候选点271b，其被标记为基准优良点(步骤268)并添加到优良点列表(步骤270)。

在取得同步化之后，方法250然后移动到该序列中的“下一个”投影图像(步骤272)。在一些实施例中，方法250将使用在同步化过程中找到的投影图像，并且从该点先向后然后向前进行。向后的处理将进行至到达具有在优良点列表中的点的投影、或者出现首次与第一投影同步的情况为止。向前的过程将进行至到达最后的投影或者直到失去同步为止。再者，当处理顺序的投影时，“下一个”的定义将取决于处理是向前还是向后。

在移至下一个投影图像之后，方法250确定该投影图像中有多少候选点(步骤274)。在一些实施例中，对于步骤274存在两类结果。下一个投影图像可以包括一个或多个候选点，或者零个候选点。图30示出了其中下一个投影图像275具有两个候选点的情形。图31示出了其中下一个投影图像275仅有一个候选点的情形，且图32示出了其中下一个投影图像275具有零个候选点的情形。

回到图24B和24C，如果下一个投影图像具有一个或多个候选点，则方法250确定在下一个投影图像中的每个候选点与优良点列表中的最近添加的“优良点”之间的、U和V位置误差相对于θ(投影角度)的导数(步骤276)。该导数有效地表示如图27中连接两点之间的线的平坦度。零导数意味着实际和理论之间的误差量在两个图像之间没有变化。

接着，方法250选择与最接近零的导数相关联的候选点(步骤278)并将该点添加到优良点列表(步骤280)。如果仅存在单个候选点，则选择该点。因为候选点列表被变平以去除误差以外的所有影响(例如，由患者运动或测量噪音导致的)，所以两个相邻的投影中的标志点之间的位置误差导数的大的变化可表明(1)图像中的结构或物体已被错误定义为标志或者(2)该标志点实际上与不同的物理标志相关联。

因此，候选点列表中与相邻的投影图像相关联并且具有接近零的导数的两个标志可能是与同一物理标志相关联的相邻标志。例如，如图30所示，下一个投影图像275包括两个候选点292a、292b。基准“优良点”293与候选点292a之间的导数比基准优良点293与候选点292b之间的导数更接近零。因此，在步骤278将选择候选点292a。作为另一个示例，如图31所示，下一个投影图像275包括一个候选点292a。在步骤280，计算基准优良点293与候选点292b之间的误差导数，并将候选点292a添加到优良点列表作为优良点。

在选择与最接近零的导数相关联的候选点之后，方法250确定该候选点的导数是否通过导数测试(步骤282)。导数测试可包括确定导数是否具有与零相差超过预定量，例如0.02mm的值。如果不是，则候选点通过导数测试，且方法250将该点标为基准优良点(步骤284)。如果该点没有通过导数测试，则方法250将该点标为“可疑点”(步骤286)。该新添加的点用于处理后续投影图像中的候选点。

可以对用于导数测试的预定量加以设定和改变，以修改方法250的灵敏度。例如，可以将预定量设定为较大值，以便允许较少的可疑点被添加到优良点列表，这可增大最终映射到特定物理标志的标志点的量。在一些实施例中，导数测试可以基于下一个投影图像中的候选点的数量、优良点列表是否包含可疑点等而变化。在一些实施例中，可不考虑对标为可疑的点进一步处理，而是总是将其视为不良的并且不将其添加到优良点列表。在这些情况下，相关联的投影将具有零个点。

再次回到图24B的步骤274，如果下一个投影图像不包含任何候选点，则方法250跳过下一个投影图像并继续使用最近添加的基准优良点(以及最近添加的可疑点，如果可用的话)以处理后续投影图像中的候选点(步骤292)。例如，如图32所示，下一个投影图像275不包含任何候选点。因此，方法250跳过下一个投影图像275并继续使用最近添加的基准优良点293(即，来自当前投影图像295)以处理后续投影图像297中的(一个或多个)候选点(例如，候选点296a)。

可疑点可以由于测量噪音导致特定点的错误放置或是由于患者或计划外的台架运动导致了物理标志的实际移动而存在。为了适当的运动修正，期望丢弃与前者有关的点，并将与后者有关的点包含在内。对可疑点的处理达到了这些运动修正目标。基本上，对处理可疑点的方案的实施是基于由于运动而可疑的点将可能具有在多个相邻投影中反映的性能这一观点。因而，一旦标记了可疑点，最终分类就取决于对将来的以及可能的另外的过去的投影性能的检验。

如上所述，如果方法250在步骤288或290中确定在优良点列表中存在先前的可疑点(也称为旧的可疑点)，并且下一个投影图像具有一个或多个候选点，则需要特殊处理(参见图24D中的步骤300和302)。如果是这样，则在该可疑点和优良点列表中最近添加的点之间形成误差导数(步骤300和302)。在形成了与可疑点相关联的导数之后，方法250确定该导数是否通过导数测试(步骤304和306)。导数测试可以包括确定导数是否具有与零相差超过预定量，例如0.02mm的值。如果导数相差少于预定值，则称可疑点通过导数测试。如果导数相差与预定值相同或者比其大，则称可疑点未通过导数测试。如果该可疑点通过导数测试，则从该点去除可疑标记(步骤308和310)。如果可疑点未通过导数测试，则从优良点列表中去除该点(步骤312和314)，这是因为其未通过两个导数测试，一个是对前一点，一个是对后一点。

如果先前的可疑点通过导数测试，并且最近添加到优良点列表中的点未通过它的导数测试，则该先前的可疑点现在成为后续投影图像中的处理时的基准优良点(步骤316)。在此情况下，假定已经发生运动并且先前的可疑点现在表示新的标志位置的最佳指示。

例如，如图33所示，可疑点321和投影图像322的候选点323之间的一阶导数最接近零，并通过导数测试。另外，优良点320和候选点323之间的对应的一阶导数未通过导数测试。因此，在步骤310，可疑点321的可疑标记被去除。先前标为可疑的点，即点321，将用作用于处理下一个投影图像的基准点(步骤316)。

存在大量关于点跟踪的可替代实施例，存在许多关于处理可疑点的可替代处理方式。在一个中，方法250被配置为，只有在对于超过一个的预定数量的顺序投影图像，可疑点与候选点的误差导数比关于基准优良点的对应导数更接近零时，才接受或标记其为优良点。例如，方法250可确定，对于一行中的三个投影图像，可疑点与候选点之间的导数是否比基准优良点与候选点之间的对应的导数更接近零。如果是，则方法250将该可疑点重新标为优良点。可对用于接受可疑点所需要的投影图像的数量加以设定和修改，以改变方法250的灵敏度(即，最终接受为优良点的可疑点的数量)。在另一个替代实施例中，可以根据任何进一步的考虑而立即移除任何可疑点。在这类情况中，可能存在更多依赖性或者重同步化。在另一个实施例中，可疑点标记将代替将单独的点标注为可疑的。在此情况下，设定可疑点标记所进行的处理可能依据候选点列表是否具有单个点或者多于一个的点而不同。

如在步骤294中确定的，如果多于特定数量的投影，例如五个，没有候选点，则断言该处理已经失去同步(步骤299)。在这些情况下，重复上述同步化处理，以重新获得同步。

如图24C和24D所示，在处理了下一个投影图像中的(一个或多个)候选点之后，方法250将下一个投影图像定义为新的当前投影图像(步骤350)，并对新的当前投影图像重复步骤272-350(即，那些适用的步骤)。该过程被重复直到从初始的当前投影图像的在一个方向上的所有投影图像都被作为“下一个投影图像”而处理(例如，到达第一个投影图像并将其作为下一个投影图像处理)(步骤352)。

在一些实施例中，方法250的投影图像步骤266-352可以完全执行两次，一次是从最低的投影编号到最高的，一次是从最高的投影编号到最低的，其中，每个方向均生成单独的优良点列表(步骤354)。一旦生成点列表，则逐个投影地将其组合成新的列表，并排除任何未存在于两个原始的列表中的点(步骤355)。该组合步骤确保最终的优良点列表中所包含的任何点都具有较高的可能性是正确的。图30-32示出了向后的或反方向上处理投影图像。图33示出了在向前的方向上处理投影图像。

在投影图像在任一方向上都进行了处理之后，方法250对每个物理标志重复步骤254-354(步骤356)。例如，如果在扫描前将九个物理标志放置于患者上面，则方法250对这九个物理标志中的每个均重复步骤254-354。在对每个物理标志执行了步骤254-354之后，方法250为每个物理标志生成优良点列表，并且每个优良点列表均包含表示在投影图像中所识别的具有与特定物理标志相关联的最高可能性的那些标志点的“优良点”(步骤360)。每个优良点列表可以包含仅零个或一个来自每个投影图像的与特定物理标志相关联的标志点。可以使用优良点列表中的“优良点”以提取患者运动信息。例如，可以将贯穿投影图像序列反映为标志点的物理标志的预期定位和贯穿该投影图像序列与物理标志相关联的标志的实际定位相比较，以确定患者移动量。

一个通常的挑战在于，如果物理标志的位置彼此靠近(例如，比患者运动所允许的范围更靠近在一起)，则难以将标志适当地分配给特定物理标志。这可以通过确保物理标志在被放置于患者上面时使其分开适当距离来加以避免。尽管如此，但方法250设计为至少部分地处理未遵循放置指导原则的情形。如果该情形出现，则跟踪过程可能不正确地分配标志点，即，被称为“轨迹跳跃(track jump)”的过程。一旦一个标志点“跳跃”到错误的轨迹，则不会存在跳回到正确的轨迹的良好基础，除非新的(以及错误的)轨迹中的一系列点丢失。尽管方法250可以处理标志轨迹上的丢失的点，但过多的丢失的点会导致轨迹跳跃。这种可能性可以通过为物理标志的候选点列表选择小范围的候选点(例如，由于患者移动量的更低的预期而导致)来最小化。

另外，在两个方向(例如，向前以及向后)上都对投影图像进行处理可以识别并排除轨迹跳跃。例如，通过将在向前方向上处理投影图像的结果与在向后方向上处理投影图像的结果相比较，方法250可以识别扩大的长度差异，这些差异是跟踪问题的指示。识别跟踪问题的另一个方式是比较用于不同物理标志的优良点列表。如果两个优良点列表都包括多于某一些(例如，一个或两个)的相同标志，则方法250可以检测并修正不正确的轨迹问题。然而，利用合理的标志点检测和分开大约一厘米或更多的物理标志，对于方法250而言通常不会存在轨迹跳跃。

虽然已经结合CB CT成像，对标志点识别方法60、标志点子像素定位方法110、物理标志3D定位方法160、以及标志点映射方法250进行了描述，但这些方法也可以用于CT、MRI、超声、其他形式的医学成像，以及非医学成像，例如摄影。例如，标志点识别方法60可以用于其中具有明确定义的形状的物体(例如，小的物体)需要与复杂的背景快速隔离的情形。另外，虽然标志点识别方法60描述为使用圆形或球形标志，但方法60也可以用于具有其他明确定义的形状的标志点。例如，可以修改步骤80中的步骤86和88(方法60的一部分)以辨识其他形状。方法60也可以使用患者内部的明确定义的解剖学上的物体作为标志。此外，方法60可以用作用于通过搜索特定标志物(landmark)的通常形状来提取解剖学上的标志物的起始点。

类似地，标志点子像素定位方法110可以应用于任何任意形状的标志，只要该形状的特点是已知的并且该标志在特定图像中的通常方位可以确定即可。各种方法都可以用于确定标志的通常方位。例如，可以使用标志点在一个(或两个)维度上的延伸以帮助限定该标志点轮廓在其他维度上的放置，从而在通常意义上对其进行评价。然后可以用简单的一阶矩计算或者细节更为复杂的曲线拟合来处理标志点轮廓。另外，可以将图像形状与在不同的方位处所预期的形状的图谱(atlas)相比较，并可以选择合适的一个。此外，可以使用CT图像数据的重建来产生标志的3D图像，以在定义和拟合适当的曲线函数之前确定其方位。

标志点子像素定位方法110在医学成像以及非医学成像中也有其他潜在的应用。例如，在牙科成像中，方法110可以用于在正牙治疗时放置托架(bracket)，其中托架的位置可以评价到大约100微米或更好的精确度。方法100也可用于准确配准在不同时间采集但利用相同的外部或内部的标志的图像。类似地，标志3D定位方法160的基本方法可用于无标志CT扫描以及其他的非CT扫描应用，例如天文学。

标志点映射方法250也有其他的潜在的应用，而不只是患者运动跟踪。例如，方法250可以用于创建用于经过非理想转动的台架的动态校准。在此情况下，方法250将应用于包含嵌入式标志的固定人体模型(phantom)。这种应用对于带有X射线图像增强器(有时称为C-Arm扫描设置)的CT系统尤其有用。此外，方法250通常可用于其中物体具有在实际实践中可能并不理想的“理想”行进的系统。

另外，上述方法可以使用各种数据处理技术。例如，标志点映射方法250可以使用各种列表以及列表和表格的组合，以用于跟踪候选点、优良点和可疑点。例如，方法250可以生成候选点列表并去除非优良点而非生成单独的优良点列表。另外，方法250可以生成一个或多个点列表，并使用标记或位，以在候选点、优良点、不良点以及可疑点之间进行区分。

此外，虽然图2将ROM模块46示出为存储有分离的模块(例如，52、54、56等)，但在一些实施例中，这些模块可以组合并分布到ROM模块46中所储存的一个或多个模块、其他存储器模块、或其组合中。例如，在一些实施例中，相同的模块执行标志识别方法和标志子像素定位方法。另外，ROM模块46可以包含用于执行除上述的那些函数外的函数的其他模块。在一些实施例中，模块可以存储在硬盘、闪存或其他半永久性介质上，并传送至RAM进行执行。

因此，本发明在其他事物当中提供用于识别图像中的标志点、确定图像中所识别的标志点的子像素中心点、为图像中所识别的标志确定3D坐标、以及将标志点映射到放置于患者上面的特定物理标志的方法和系统。本发明的各种特征和优点将在随附权利要求书中阐明。

Claims (17)

1.一种用于识别图像中的标志关联点的方法，该方法由包括扫描器、带有电子处理单元的计算机、和存储有可由所述电子处理单元执行的标志点识别模块的存储器模块的成像系统执行，该方法包括：

在所述计算机处，基于由所述扫描器产生的图像数据获得图像；

利用所述电子处理单元，对所述图像应用具有X射线补偿的二维高通滤波器以抑制背景，从而生成背景抑制后的图像；

利用所述电子处理单元，对所述图像应用二维低通滤波器以生成低通图像，并且其中，应用所述二维高通滤波器包括将所述图像与所述低通图像组合；以及

利用所述电子处理单元，从所述背景抑制后的图像中提取标志关联点。

2.根据权利要求1所述的方法，其中，应用所述二维高通滤波器包括对所述图像顺序地应用两个或更多个一维高通滤波器，所述两个或更多个一维高通滤波器中的各个是在不同的二维方向上对所述图像执行的。

3.根据权利要求1所述的方法，其中，应用所述二维高通滤波器包括对所述图像应用两个或更多个一维高通滤波器，以及其中，应用各个一维高通滤波器包括组合下列各项中的至少一项：

所述图像的一维低通滤波后的表示与所述图像，以及

来自所述一维高通滤波器的较早应用的临时图像的一维低通滤波后的表示与所述临时图像。

4.根据权利要求3所述的方法，其中，所述图像的所述一维低通滤波后的表示通过应用累积均值滤波来生成，所述滤波包括：

使用所述电子处理单元，从所述图像生成累积图像，其中所述累积图像中的各个像素表示在所述滤波器方向上的较低像素的累积和；

利用所述电子处理单元，将所述累积图像以滤波因子移位，以生成移位后的图像；

利用所述电子处理单元，将所述移位后的图像减去所述累积图像，以生成差值图像；

利用所述电子处理单元，将所述差值图像除以所述滤波因子。

5.根据权利要求3所述的方法，其中，将所述图像与所述图像的所述一维低通滤波后的表示组合考虑了X射线的非线性衰减，并且包括将所述图像逐个像素地除以所述图像的所述一维低通滤波后的表示。

6.根据权利要求5所述的方法，其中，逐个像素的除法使用所述图像的所述一维低通滤波后的表示的像素作为被除数来进行。

7.根据权利要求3所述的方法，其中，将所述图像与所述图像的所述一维低通滤波后的表示组合包括减法。

8.根据权利要求1所述的方法，其中，从所述背景抑制后的图像中提取所述标志关联点包括在所述背景抑制后的图像中定义局部极大值。

9.根据权利要求8所述的方法，其中，在所述背景抑制后的图像中定义所述局部极大值包括在所述背景抑制后的图像中识别含有具有较小的像素值的八个相邻像素的像素。

10.根据权利要求8所述的方法，其中，从所述背景抑制后的图像中提取所述标志关联点进一步包括基于所述局部极大值来生长候选区域。

11.根据权利要求10所述的方法，其中，从所述背景抑制后的图像中提取所述标志关联点进一步包括对每个候选区域应用形状标准。

12.根据权利要求11所述的方法，其中，对每个候选区域应用所述形状标准包括确定每个候选区域的像素面积和由相应候选区域的平均半径的平方乘以π所定义的面积之比。

13.根据权利要求10所述的方法，其中，从所述背景抑制后的图像中提取所述标志关联点进一步包括通过比较相邻投影图像中的候选区域并且在候选区域之间应用接近度标准来排除候选区域。

14.一种成像系统，包括扫描器、带有电子处理单元的计算机、和存储有可由所述电子处理单元执行的标志点识别模块的存储器模块，其中，所述标志点识别模块当执行时被配置为：

在所述计算机处，基于由所述扫描器产生的图像数据获得图像；

对所述图像应用具有x射线补偿的二维高通滤波器以抑制背景，从而生成背景抑制后的图像；

从所述背景抑制后的图像中提取标志关联点；

对所述图像应用二维低通滤波器以生成低通图像；以及

通过将所述图像与所述低通图像组合来应用所述二维高通滤波器。

15.根据权利要求14所述的成像系统，其中所述标志点识别模块当执行时进一步被配置为：通过对所述图像顺序地应用两个或更多个一维高通滤波器来应用所述二维高通滤波器，所述两个或更多个一维高通滤波器的各个是在不同的二维方向上对所述图像执行的。

16.根据权利要求14所述的成像系统，其中，所述标志点识别模块当执行时进一步被配置为通过在所述背景抑制后的图像中定义局部极大值来从所述背景抑制后的图像中提取所述标志关联点，其中定义所述局部极大值包括在所述背景抑制后的图像中识别含有具有较小的像素值的八个相邻像素的像素。

17.根据权利要求16所述的成像系统，其中所述标志点识别模块当执行时进一步被配置为通过下列步骤来从所述背景抑制后的图像中提取所述标志关联点：

基于所述局部极大值来生长候选区域，

对每个候选区域应用形状标准，以及

通过比较相邻投影图像中的候选区域并且在候选区域之间应用接近度标准来排除候选区域。