锥束计算机断层扫描图像注册与重叠方法与流程

本发明涉及一种体图像的注册与重叠方法，具体而言，涉及一种锥束计算机断层扫描图像注册与重叠方法。

背景技术：

锥束计算机断层扫描(cbct)图像是口腔临床中常用于获取特定病人的颅面形态信息的一种体图像。与传统的二维放射图像相比，cbct图像可以提供病人的完整的三维几何信息。通过对比治疗前后的cbct图像，可以对治疗所引起的三维颅面形态的变化进行量化分析。由于病人在治疗过程中的机体组织生长以及在不同的图像采集过程中的姿态变化，在利用治疗前后的cbct图像进行治疗评价以及量化分析之前，必须对治疗前后的三维cbct图像注册与重叠，即，进行体图像的注册与重叠。

医学图像的注册与重叠通常用于建立不同的体图像之间的解剖结构的对应关系。现有技术中的三维cbct图像的注册与重叠仍面临如下挑战：1)体图像属于大规模数据，仅在中等分辨率的颅面cbct图像当中就包含有百万以上的体素，因此传统医学图像注册中使用的测度函数(例如，互信息以及规范互相关等测度)相对费时并易于陷入局部极小的情况；2)在颅面cbct图像所记录的解剖结构当中，一部分结构在治疗过程中保持形态不变，而另一些结构则会随时间而发生变化，因此在体图像的注册与重叠中需要对于不同的结构区别对待，然而手工地进行解剖结构的标注耗时、费力并依赖于操作者，因此需要一种有效的标注工具用于定义感兴趣结构(soi)；3)颅面cbct图像通常由于放射剂量小而图像分辨率较低且图像质量较差，因而在图像中，部分骨性结构会由于与周围组织的灰度相似而变得模糊，难以在参考cbct图像和目标cbct图像中确定其对应关系，因此需要一 种鲁棒的对应估计算法以实现可靠的cbct图像的注册与重叠。

医学图像处理领域中的自动图像注册与重叠已经进行了多年的研究。近年来，颅面体图像的注册与重叠大多使用基于体素的各种方法。互信息(mutualinformation)是在体图像的注册与重叠方法中经常使用的测度函数，其用于估计体图像的统计灰度并试图找到使各个体图像之间信息最大的空间变换。由于互信息是一个全局度量，因此用于估计局部最优的变换相对较难。进而，引入了体图像相位以替代图像的灰度，通过计算局部相位的互信息使得结构关联最大化。通常，将颅底部分的体素看作是稳定的结构并进行刚性注册与重叠的方法大多基于互信息测度函数。在这些方法中，颅面体图像在重叠后可用于评价下颌区，或者可用于比较颌面手术中骨骼以及软组织的形态变化。一些基于互信息测度函数的体图像注册与重叠软件相继出现，其中mirit软件用于计算正颌手术病人在治疗前后的体图像的刚性注册；ondemand3d软件用于计算基于前颅底的cbct图像重叠；maxilim软件用于对颅底的三维模型进行重叠。

此外，基于标记的方法也可以用于重叠体图像并度量颅面结构由于生长以及治疗所造成的变化。标记结合普氏(procrustes)测度几乎被用于所有的形态学研究。普氏分析以及薄板样条可用于样本曲线轮廓的多变量分析，其中以手工的方式选取三维标记并计算普氏几何变换，以估计下颌由于生长与治疗所造成的变化。

方向无关的三维尺度不变特征变换(scaleinvariantfeaturetransform，sift)通常用于临床ct图像的注册与重叠。例如，stephaneallaire等人于2008年在“2012ieee计算机视觉与模式识别研讨会计算机学会会议(computersocietyconferenceoncomputervisionandpatternrecognitionworkshops)”上发表的题为“fullorientationinvarianceandimprovedfeatureselectivityof3dsiftwithapplicationtomedicalimageanalysis”的文章中，给出了一种针对医疗图像分析应用的三维sift方法。

在现有技术中，通常直接对原始体图像中的稠密体素进行操作。 由于数据量较大且测度函数复杂，因而计算效率较低。此外，在现有技术中，对于感兴趣区域的手工标注不仅操作繁复而且依赖于操作者。在现有的基于标记的注册与重叠方法中，需要手工操作来定义三维标记。另外，在三维描述子的计算中，方向确定以及梯度直方图的估计较为复杂并且计算效率较低。

技术实现要素：

提出了本发明的技术方案，以解决现有技术中的上述问题之一或其他技术问题。

根据本发明的一个方面，提供了一种体图像的注册与重叠方法，包括：分别在参考体图像和目标体图像中提取参考图像子集和目标图像子集；使用半监督学习(semi-supervisedlearning)方法分别在参考图像子集和目标图像子集中定义感兴趣结构(soi)；计算参考图像子集的soi和目标图像子集的soi中的各个体素的球灰度积分算子；基于球灰度积分算子通过协同嵌入来求解参考图像子集与目标图像子集之间的对应关系；基于得到的对应关系计算参考体图像与目标体图像之间的刚性变换参数；以及将计算得到的刚性变换参数应用于参考体图像，并将刚性变换后的参考体图像与目标体图像重叠。

根据本发明的实施例，分别在参考体图像和目标体图像中提取参考图像子集和目标图像子集的步骤可以包括：分别对参考体图像和目标体图像进行高斯平滑；针对参考体图像计算各阶高斯平滑图像之间的一阶参考差分图像；针对目标体图像计算各阶高斯平滑图像之间的一阶目标差分图像；提取一阶参考差分图像的局部极值作为参考体图像的参考图像子集；以及提取一阶目标差分图像的局部极值作为目标体图像的目标图像子集。

根据本发明的实施例，可以针对参考体图像和目标体图像分别设定阈值，并且从提取的参考图像子集和目标图像子集中去除其灰度值小于阈值的体素。

根据本发明的实施例，在提取参考图像子集和目标图像子集的步骤之前，所述方法还可以包括利用直方图匹配算子对参考体图像和 目标体图像进行处理。

根据本发明的实施例，可以基于预先输入的soi的体积先验来使用半监督学习方法分别在参考图像子集和目标图像子集中定义soi。可以通过预先对属于不同类型的soi的体素进行计数得到所述soi的体积先验。

根据本发明的实施例，计算体素的球灰度积分算子的步骤可以包括：以待计算的体素作为球心并且以预定长度作为半径形成一个球，待计算的体素为球心体素；将形成的球等分成n个体锥；根据公式(1)对n等分后的每个体锥内的体素的灰度进行积分：

公式(1)

其中，r表示球的半径，bj表示第j个体锥，1≤j≤n，v表示体锥bj所涵盖的空间中的各个体素，函数i(v)返回的值为体素v与球心体素的灰度差，uj,r表示将半径为r的球等分成n的体锥后，第j个体锥bj内的体素的灰度积分；以及将针对n个体锥计算后的n个结果排序，得到n维向量ur＝(u1,r,u2,r,……,un,r)，其中，ur表示球心体素对于半径为r的球的灰度积分。

根据本发明的实施例，计算体素的球灰度积分算子的步骤还可以包括：以待计算的体素作为球心选取m个不同的半径，得到m组n维向量1≤i≤m。

根据本发明的实施例，n＝40，m＝6，ri∈[2,10]，其中，ri的单位为体素数量。

根据本发明的实施例，基于球灰度积分算子通过协同嵌入来求解参考图像子集与目标图像子集之间的对应关系的步骤可以包括：建立与参考图像子集和目标图像子集对应的图结构ge(a,ψ)，其中，结点集a为考图像子集和目标图像子集中的体素，边连接ψ包含各图像子集内部的连接以及各图像子集之间的连接；对于图结构ge(a,ψ)，基于图像子集的各个体素之间的欧式距离和与各个体素对应的球灰度积分算子之间的距离计算在各图像子集内部连接各个体素的边权值，并且基于与各个体素对应的球灰度积分算子之间的距离计算在各图像子集之间连接各个体素的边权值；基于计算得到的各个边权 值建立协同嵌入的目标函数；通过对目标函数进行拉普拉斯矩阵特征分解，求得与最小特征值对应的特征向量，即为参考图像子集和目标图像子集在协同嵌入的空间中的坐标；以及在协同嵌入空间中定义阈值，当参考图像子集中的体素和目标图像子集中的体素在协同嵌入空间中的距离小于阈值时，认定这一对体素是对应的。

根据本发明的实施例，参考体图像与目标体图像之间的刚性变换参数可以包括从参考体图像变换至目标体图像的旋转矩阵和平移向量。

附图说明

通过参照附图对本发明构思的示例性实施例进行详细描述，本发明构思的示例实施例的上述和其他特征和优点将变得更加清楚。附图旨在描述本发明构思的示例实施例，而不应当解释为限制权利要求的预期范围。在附图中：

图1示出了根据本发明构思的实施例的体图像的注册与重叠方法的流程图；以及

图2a和图2b示出了根据本发明构思的实施例的体图像的注册与重叠方法中使用的球灰度积分算子的示意图。

具体实施方式

下面将参照附图详细描述本发明构思的示例性实施例。

然而，本发明构思可以按照许多不同的形式示例，并且不应理解为限于本文所阐述的示例性实施例。此外，提供这些实施例是为了使得本公开将是彻底而完整的，并且将向本领域的技术人员充分地传达本发明构思的范围。在附图中，相同的附图标记将始终用于指示相同或相似的元件。

应当理解的是，虽然在本文中使用了术语第一、第二等来描述各个元件，然而这些元件不应当被这些术语所限定。这些术语仅用于将一个元件与另一个元件区分开。例如，第一元件可被称作第二元件，并且类似地，第二元件可被称作第一元件，而没有背离本发明构思的 示例实施例的范围。如在本文中所使用的那样，术语“和/或”包括一个或多个所列相关项目的任意和全部组合。

本文所使用的术语仅用于描述本文所阐述的示例性实施例，而非旨在限定本发明构思。如同本文所使用的那样，除非上下文中另外明确表示，否则单数形式“一个”、“一”和“该”也旨在包括复数形式。还应当理解，当术语“包括”、“提供有”和/或“具有”用于示例性实施例中时，其指示了存在所述特征、整体、步骤、操作、元件和/或部件，但并不排除存在或增加其他一个或多个特征、整体、步骤、操作、元件、部件和/或它们的组合。

还应当注意到，在一些替代实现方式中，所示出的功能/动作会按照与图中标注的顺序不同的方式发生。例如，连续示出的两个图实际上可以基本同时执行，或者有时可以按照相反的顺序执行，这取决于所涉及的功能/动作。

除非另有说明，否则本文使用的包括说明性术语或技术术语的所有术语应该被理解为具有对于本领域普通技术人员之一明显的含义。另外，在普通词典中定义并且在以下描述中使用的术语应该被理解为具有与相关描述中使用的含义等同的含义，并且，除非本文中另有说明，否则所述术语不应被理解为理想的或过于正式的。还应当认识到，在没有进行特定声明的情况下，本文所提及的“图像”为“体图像”。

图1示出了根据本发明构思的实施例的体图像的注册与重叠方法的流程图。

参照图1，根据本发明构思的实施例的体图像的注册与重叠方法100包括：分别在参考体图像和目标体图像中提取参考图像子集和目标图像子集(s110)；使用半监督学习方法分别在参考图像子集和目标图像子集中定义感兴趣结构(s120)；计算参考图像子集的感兴趣结构和目标图像子集的感兴趣结构中的各个体素的球灰度积分算子(s130)；基于球灰度积分算子通过协同嵌入来求解参考图像子集与目标图像子集之间的对应关系(s140)；基于得到的对应关系计算参考体图像与目标体图像之间的刚性变换参数(s150)；以及将计算得 到的刚性变换参数应用于参考体图像，并将刚性变换后的参考体图像与目标体图像重叠(s160)。

在提取图像子集的步骤(s110)中，分别从参考体图像与目标图像中提取各自的图像子集。参考体图像和目标体图像可以是锥束计算机断层扫描(cbct)体图像。

具体而言，可以对每个cbct体图像进行高斯平滑，并计算各阶高斯平滑图像之间的一阶差分图像。也就是说，针对参考体图像计算各阶高斯平滑图像之间的一阶参考差分图像，并且针对目标体图像计算各阶高斯平滑图像之间的一阶目标差分图像。随后，分别提取一阶参考差分图像的局部极值和一阶目标差分图像的局部极值以构成参考图像子集和目标图像子集，即，从cbct体图像提取的图像子集可以由极值位置上的各个体素进行定义。通过从cbct体图像提取图像子集，可以极大地降低计算数据的规模。

根据本发明的实施例，在从cbct体图像提取图像子集之前，可以利用直方图匹配算子对参考cbct体图像与目标cbct体图像进行处理(s100)，如图1所示。然而，本发明构思不限于此，而是可以利用其他匹配算子对体图像进行处理，或者不对体图像进行处理。

此外，由于基于差分图像的局部极值提取的体素可能位于体图像背景区域，从而使得计算过程不稳定，因此，对于参考cbct体图像与目标cbct体图像可分别定义阈值，以从提取的图像子集中去除其灰度值小于预定阈值的体素。

在定义感兴趣结构(soi)的步骤(s120)中，分别针对参考图像子集和目标图像子集使用半监督学习方法。

在口腔临床应用中，体图像的注册与重叠通常会基于稳定的解剖结构，例如，前颅底、颧弓等。这些稳定的解剖结构构成了可以在体图像的注册与重叠方法中使用的soi。利用半监督学习方法，可以将用户给定的少量结构标记扩散到整个子集。

通过半监督学习方法，可以使用少量的已标记数据来训练大量的未标记数据。半监督学习介于无监督学习(即，不使用任何已标记的训练数据)和监督学习(即，全部使用已标记的训练数据)。研究 表明，通过结合少量的已标记数据，可以大大提高未标记数据的学习准确性。

具体而言，首先由用户对少量的切面图像(即，构成参考体图像和目标体图像的切面图像，其为平面图像)进行标记，标记的内容可以包括soi的起始切面图像、终止切面图像以及中间切面图像。随后，基于已标记的体素，对图像子集应用半监督学习方法，使得标记可以扩散到图像子集中的所有体素。

对体素进行标记的过程即为根据体素所属的解剖结构对该体素进行分类的过程。在口腔临床应用中，soi的类型可以包括：前颅底(cb)、左颧弓(zl)，右颧弓(zr)以及其它结构(ob)。因此，可以使用cb、zl、zr和ob分别对体素进行标记。

本申请的发明人注意到，在针对口腔临床应用中的cbct体图像使用现有技术中的半监督学习方法时，会出现标记扩散退化的现象，例如，由于soi与其他骨骼或组织紧密相邻，在使用现有技术中的半监督学习方法时，会使标记扩散至soi之外的其他骨骼或组织。为了避免标记扩散退化的出现，根据本发明的实施例，在传统的半监督学习方法中引入了soi的体积先验(volumepriori)作为边界约束条件，即，基于预先输入的soi的体积先验来使用半监督学习方法分别在参考cbct体图像的图像子集和目标cbct体图像的图像子集中定义soi。可以通过预先对属于不同类型的soi的体素进行计数来得到所述soi的体积先验。具体而言，任一类型的soi的体积先验可以表示为属于该类型的soi的体素数量与图像子集的总体素数量的比值。例如，经过预先统计，属于前颅底(cb)的体素数量为ncb，图像子集的总体素数量为ntotal，则前颅底的体积先验δcb可以表示为ncb/ntotal。

在对图像子集中的所有体素完成标记之后，可以选择图像子集中属于某一特定的soi(例如，前颅底、颧弓等)的体素来进行重叠参数的估计。

在计算体素的球灰度积分算子的步骤(s130)中，针对参考图像子集的soi和目标图像子集的soi中的各个体素来计算球灰度积分算子。体素的球灰度积分算子在随后建立参考图像子集与目标图像子 集之间的对应关系(步骤s140)时可以用作判定两个体素(即，参考体图像中的体素和目标体图像中的体素)是否对应的原始特征空间。使用体素的球灰度积分算子作为原始特征空间，可以描述图像子集中的体素的空间分布和纹理上下文信息。此外，与现有技术中的其他可以用作原始特征空间的算子相比，体素的球灰度积分算子的计算量相对较小。

图2a和图2b示出了根据本发明构思的实施例的体图像的注册与重叠方法中使用的球灰度积分算子的示意图。图2a以平面图的方式示意性地示出了对于不同半径的球计算球灰度积分算子，图2b是示意性地示出了将不同半径的球n等分后的体锥的透视图。

参照图2a和图2b，在计算体素的球灰度积分算子时，以该体素为球心(以下也称作“球心体素”)并以r为半径形成一个球，并且将所形成的球n(例如，n＝40)等分，如图2b所示。例如，e.b.saff和a.b.j.kuijlaars于1997年12月在“themathematicalintelligencer”第19卷、第5-11页上发表的题为“distributingmanypointsonasphere”的文章中，给出了一种将球面进行n等分的方法。

随后，参照图2a，对n等分后的每个体锥内的体素的灰度进行积分。具体的积分公式参见以下公式(2)：

公式(2)

其中，表示将半径为ri(1≤i≤m)的球n等分后，第j(1≤j≤n)个体锥内的体素的灰度积分；bj表示第j个体锥；函数i(v)返回的值为体素v与球心体素的灰度差。

由于体图像以离散的方式表示，因此在实际的计算过程中，与体锥bj对应的灰度积分可近似表示为以下公式(3)：

公式(3)

其中，v表示体锥bj所涵盖的空间中的各个体素。

分别计算与每个体锥bj(1≤j≤n)对应的灰度积分，并对计算后的结果进行排序(例如，由低到高排序或由高到低排序)，以避免旋转对体素特征描述产生的影响。由此，对于半径为ri的n等分球， 可以得到n维特征描述向量以作为球心体素的球灰度积分算子。

根据本发明的一个实施例，可以选取多个不同的ri(1≤i≤m)来得到n×m维的特征描述矩阵以作为球心体素的球灰度积分算子。根据本发明的一个实施例，m可以为6，半径ri∈[2,10]，ri的单位为体素数量。

以球灰度积分作为算子，可以描述图像子集中的体素的空间分布以及纹理上下文信息。但是，在这样的原始特征空间中，对称的体素和结构会引起映射混淆，并且这样的原始特征空间维度较高。为了避免对称体素和结构在原始特征空间中引起的映射混淆，同时降低特征空间的维度，可以通过协同嵌入(jointembedding)算法将参考图像子集和目标图像子集嵌入到维度更低的流形(manifold)，并在协同嵌入空间求解参考图像子集和目标图像子集的对应关系。

在此，为了建立有语义的对应关系，协同嵌入必须满足以下两个要求：在参考图像子集和目标图像子集中相似的体素在协同嵌入空间中必须足够接近；以及各个图像子集的全局结构在协同嵌入流形后能够保持。

具体方式为，建立与参考图像子集sref和目标图像子集star对应的图结构ge(a,ψ)，其中结点集a为考图像子集sref和目标图像子集star中的体素，边连接ψ包含各个图像子集内部的连接ψintra以及各图像子集之间的连接ψinter，并且ψ＝ψinter∪ψinter。

对于如此构造的图结构ge(a,ψ)，基于图像子集的各个体素之间的欧式距离和与各个体素对应的球灰度积分算子之间的距离计算在各图像子集内部连接各个体素的边权值w^intra，并且基于与各个体素对应的球灰度积分算子之间的距离计算在各图像子集之间连接各个体素的边权值w^inter。因此，在求解参考图像子集sref和目标图像子集star的对应关系时，一方面考虑到图像子集的各个体素之间的欧式距离，另一方面考虑到与各个体素对应的球灰度积分算子之间的距离，即，基于球灰度积分算子通过协同嵌入来求解参考图像子集与目标图像子集之间的对应关系(步骤s140)。

基于计算得到的各个边权值w^intra和w^inter建立协同嵌入的目标函数。根据本发明的实施例，目标函数e(f)可以由以下公式(4)表示：

公式(4)

其中，f是将体素从原始特征空间变换到流形的投影函数，sk和sl表示图像子集sref和目标图像子集star，即，k,l∈{ref,tar}。当并且时，oij＝1，否则oij＝0，其中soik和soil表示图像子集sref和目标图像子集star中的soi。当k和l同时表示参考图像或同时表示目标图像时(即，k＝l)，w^kl＝w^intra；当k和l分别表示参考图像和目标图像时(即，k≠l)，w^kl＝w^inter。常数系数αkl用于调节图像子集的结构保持和图像子集之间的相似体素的近邻嵌入。

通过上述方式构建的目标函数e(f)可以满足协同嵌入的两个要求，即，在参考图像子集和目标图像子集中相似的体素在协同嵌入空间中必须足够接近(即，k≠l的情况)；以及各个图像子集的全局结构在协同嵌入流形后能够保持(即，k＝l的情况)。

由于协同嵌入的目标函数与拉普拉斯特征映射一致，通过对拉普拉斯矩阵进行特征分解，得到与最小特征值对应的特征向量即为参考图像子集sref和目标图像子集star在协同嵌入的空间中的坐标，即，目标函数的最优解。

由于参考图像子集sref和目标图像子集star的定义是独立进行的，因此不能保证参考图像子集sref中的体素和目标图像子集star中的体素一一对应。可以在嵌入空间中定义阈值，只有当参考图像子集sref中的体素和目标图像子集star中的体素在嵌入空间中的距离小于该阈值时，才认为这一对体素是对应的。

基于得到的对应关系，可以计算参考体图像与目标体图像之间的刚性变换参数(s150)，即，求解从参考体图像变换至目标体图像的旋转矩阵和平移向量。可以依据由不同的soi的协同嵌入所获取的图像子集之间的对应关系来计算刚性变换参数。

随后，可以将计算得到的刚性变换参数应用于参考体图像，并将刚性变换后的参考体图像与目标体图像重叠(s160)。

可以对重叠后的稳定骨结构(例如，前颅底)表面的距离进行 计算，以对基于协同嵌入的体图像重叠的精度进行验证。经实验验证，该距离在0.5mm以下，完全可以满足口腔临床对于体图像重叠的精度要求。

根据本发明实施例的体图像的注册与重叠方法，可以快速有效地对治疗前后的cbct体图像进行刚性注册与重叠。通过提取图像子集，有效地降低了计算中的数据规模。

此外，本发明还提出了用于描述图像子集的旋转无关的球灰度积分算子，以对体素的上下文信息进行描述。由于球灰度积分算子仅包含线性计算，从而有效克服传统三维特征描述中的方向估计、梯度直方图等复杂的计算。

对于稳定结构(或感兴趣结构)的标注采用了半监督学习方法(或半监督的随机游走算法)，从而避免了全部由用户手工交互地进行标注，降低了对于操作者的依赖性并减轻了交互繁复的问题。

提出了基于集协同嵌入算法来求解图像子集的对应关系的方法，其中通过对协同嵌入空间中的少数占优分量的操作就可以获取各图像子集中的体素之间的对应关系，与在原始特征空间中进行注册的计算相比，可有效降低数据维度。同时在协同嵌入空间中保持图像子集的整体结构，从而避免了参考体图像与目标体图像中的对称体素或者结构所造成的映射混淆，进而实现可靠重叠。

根据本发明实施例的体图像的注册与重叠方法，有效地克服了依据传统方法的体图像重叠中的计算量大以及容易陷入局部最小的问题。

各个示例性实施例的各种优点和效果并不限于上述描述，并且可以通过本公开中的具体实施例的解释而易于理解这些优点和作用。

虽然上文示出并描述了各个示例性实施例，但是本领域的技术人员将显而易见，可以做出修改和变化而没有背离如权利要求所定义的本发明构思的范围。