CN104021547B - 肺部 ct 的三维配准方法 - Google Patents

Abstract

一种肺部CT的三维配准方法，其包括以下步骤：载入第一肺部CT体数据和第二肺部CT体数据；以互信息为配准测度，对第一肺部CT体数据和第二肺部CT体数据进行三维刚性配准，使二者中的肺部区域整体对齐，记录X轴、Y轴和Z轴三个方向的位移参数和旋转参数，作为第一配准参数；在刚性配准的基础上进行非刚性配准，记录非刚性配准的最优位移场，作为第二配准参数；用第一配准参数和第二配准参数对第一肺部CT体数据和第二肺部CT体数据中的结节进行匹配，建立对应关系。本发明方法能够实现精确配准，并利用配准结果将同一位置的结节病灶在不同CT数据中的成像位置进行自动对应，从而能够方便医生对结节病情发展进行对比。

Description

肺部CT的三维配准方法

技术领域

本发明属于计算机视觉及图像处理技术领域，具体涉及一种肺部CT的三维配准方法。

背景技术

根据WHO2008年的全球肿瘤流行病统计数据(GLOBOCAN2008)，仅2008年全球肿瘤病例为1270万，死亡760万例[1]，其中，中国所有的肿瘤新病例为280万，是美国新发肿瘤病例的2倍(美国为140万)；同年中国肿瘤死亡病例为196万，而美国的死亡病例为57万，二者的差距为3.4倍。恶性肿瘤占我国居民死亡原因的20％以上。据WHO国际肿瘤研究理事会负责人Bernard Stewart报道，到2020年全球肿瘤发病率还将上升50％。

提高肿瘤病人的五年存活率的关键在于如何尽早的发现并确诊肿瘤病情，并进行治疗。随着CT技术的发展，高分辨率CT因为其良好的成像质量，成为肿瘤筛查的重要工具，被越来越广泛地在临床上采用，但是随之而来的问题是医生阅片工作量的增加，尤其是在对比阅片时在大量CT数据中手动进行相同解剖位置图像匹配。

因此利用计算机辅助诊断工具帮助医生高效准确地对肿瘤病情进行诊断有着巨大的现实意义。目前相关的计算机辅助诊断工具的研究都集中如何从一套CT数据中进行肺结节的识别和分割，以及根据特征对结节的良恶性进行判断。在肿瘤的诊断和治疗中，对比同一位置病灶特征的变化有着重要的临床意义，而在如何提高医生在进行病情对比时的阅片效率方面的研究很少。

在一些PACS(Picture Archiving and Communication Systems)系统中，虽然提供了自动定位相近断层图像的功能，但是这种功能精确度不高，并且只是进行了断层与断层之间的匹配，如果一个断层中出现多个结节，这种方法并不能提供有效的对应结果。

发明内容

本发明的目的是提供一种肺部CT的三维配准方法，该三维配准方法能够将同一位置的结节在不同CT数据中的成像位置进行自动对应。

本发明中提出的配准技术应用的对象是高分辨率肺部CT图像，其技术思路为：首先进行刚性配准，通过互信息将CT体数据中的肺部区域进行整体对齐；然后在刚性配准的基础上利用提取的匹配特征点对结合局部单点互信息(SMI)最大化的方式进行非刚性配准，得到最终的精确配准结果；然后利用配准结果提供的空间变换信息，进行相同结节成像位置的对应，从而方便结节的特征变化进行对比。

本发明的具体技术方案为：

一种肺部CT的三维配准方法，包括以下步骤：

载入第一肺部CT体数据和第二肺部CT体数据；

以互信息为配准测度，对第一肺部CT体数据和第二肺部CT体数据进行三维刚性配准，使二者中的肺部区域整体对齐，记录X轴、Y轴和Z轴三个方向的位移参数和旋转参数，作为第一配准参数；

在刚性配准的基础上进行非刚性配准，记录非刚性配准的最优位移场，作为第二配准参数；

用第一配准参数和第二配准参数对第一肺部CT体数据和第二肺部CT体数据中的结节进行匹配，建立对应关系。

在上述的肺部CT的三维配准方法中，优选地，所述非刚性配准包括特征点提取、特征点匹配和位移场搜索。

在上述的肺部CT的三维配准方法中，优选地，采用三维的算子提取特征点，特征点匹配中被匹配为特征点对的两个特征点r、f满足以下两个条件：一、参考图像中的特征点r在浮动图像中最邻近的点是特征点f，且特征点f在参考图像中的最邻近的点也是特征点r；二、最邻近特征点与次临近特征点间的距离比小于预设的第一阈值。

在上述的肺部CT的三维配准方法中，优选地，采用BBF算法确定最邻近的特征点。

在上述的肺部CT的三维配准方法中，优选地，所述非刚性配准用互信息结合特征点限制作为相似性测度。

在上述的肺部CT的三维配准方法中，优选地，所述位移场搜索包括以下步骤：计算位移场变化量；将位移场变化量与位移场相加，更新位移场；利用特征点对位移场进行限制，包括计算特征点x处的位移场Uⁿ⁺¹(x)，计算该位移场Uⁿ⁺¹(x)与V(x)的差值，该差值小于预设的第二阈值时，使位移场Uⁿ⁺¹(x))V(x)，否则Uⁿ⁺¹(x)保持不变，V(x)为构成特征点对的两个特征点之间的位移，该特征点对包括特征点x；计算浮动图像插值结果，计算该插值结果与参考图像之间的互信息，当互信息停止增大或迭代次数达到限定次数时，将当前位移场作为最优位移场，否则重复上述步骤。

在上述的肺部CT的三维配准方法中，优选地，更新位移场前还包括平滑位移场变化量的步骤，在利用特征点对位移场进行限制后还包括对位移场进行加权高斯平滑的步骤，所述加权高斯平滑在特征点处采用第一权值，在其它位置采用第二权值。

在上述的肺部CT的三维配准方法中，优选地，采用对称搜索的方式计算位移场变化量。

在上述的肺部CT的三维配准方法中，优选地，该三维配准方法还包括在进行配准前将第一肺部CT体数据和第二肺部CT体数据每个体素代表的物理尺寸降采样到相同尺度上。

在上述的肺部CT的三维配准方法中，优选地，采用浮动图像与参考图像之间的联合灰度分布直方图计算互信息，其中采用PV插值的方式来计算联合灰度分布直方图。

本发明方法能够实现精确配准，并利用配准结果将同一位置的结节病灶在不同CT数据中的成像位置进行自动对应，能够方便医生对结节病情发展进行对比，从而能够大大提高医生的诊断效率。

附图说明

图1为一些实施例肺部CT的三维配准方法的流程图；

图2为采用算子提取特征点匹配后的结果；

图3为浮动图像配准后的插值图像与参考图像相减图。

具体实施方式

下面结合附图和实施例对本发明进一步说明。这些更详细的描述旨在帮助理解本发明，而不应被用于限制本发明。根据本发明公开的内容，本领域技术人员明白，可以不需要一些或者所有这些特定细节即可实施本发明。而在其它情况下，为了避免将发明创造淡化，未详细描述众所周知的操作过程。

如图1所示，本肺部CT的三维配准方法包括以下步骤：

步骤S1、载入第一肺部CT体数据和第二肺部CT体数据。第一肺部CT体数据和第二肺部CT体数据是同一被测者在不同时间或不同检测条件下的肺部CT体数据。

步骤S2、以互信息为配准测度，对第一肺部CT体数据和第二肺部CT体数据进行三维刚性配准，使二者中的肺部区域整体对齐，记录X轴、Y轴和Z轴三个方向的位移参数和旋转参数，作为第一配准参数。

三维刚性配准采用互信息作为相似性测度，因此不要对肺实质进行预分割处理。互信息的定义式如下：

其中F表示浮动图像，R表示参考图像，T是将浮动图像映射到参考图像的变换矩阵，hTCg·O表示在变换矩阵T下，浮动图像与参考图像之间的联合灰度分布直方图，并由此估计浮动图像与参考图像间灰度的联合概率分布p_FR，T(f，r)及边缘概率分布p_F，T(f)和p_R，T(r)。

在计算直方图分布时，采用了PV插值的方式，PV插值是根据浮动图像映射到目标图像上的对应坐标值，利用线性加权方式计算联合直方图的分布。因此PV插值并不是真正意义上的插值，没有生成真正的插值图像，从而避免了生成插值图像时的舍入误差对直方图的估算精度造成影响。

刚性变换矩阵是具有三个位移参数和三个旋转参数的刚性变换矩阵，具体的变换方程可表示为：

V_R(P_R-C_R)＝R_x(φ_x)R_y(φ_y)R_z(φ_z)V_F(P_F-C_F)+t(t_x，t_y，t_z) (5)

其中V_R和V_F是3×3的对角阵，每个对角线元素代表分别表示一个体素在x、y、z三个方向上表示的物理尺寸；P_R、P_F分别代表参考图像和浮动图像中对应的两个体素；C_R、C_F分别表示参考图像和浮动图像的旋转中心；R_x(φ_x)、R_y(φ_y)、R_z(φ_z)分别表示绕x轴、y轴、z轴三个方向旋转的3×3旋转矩阵。

步骤S3、在刚性配准的基础上进行非刚性配准，记录非刚性配准的最优位移场，作为第二配准参数。

在较佳实施例中，非刚性配准包括三个步骤：特征点提取、特征点匹配、位移场搜索。下面分别对特征点提取、特征点匹配、位移场搜索的一些较佳实施方式进行详细说明。

(1)特征点提取

特征点提取采用了三维的角点提取算子，其定义如式(6)所示：

其中

为体数据中在x点的梯度向量，T为矩阵转置，C为3×3的矩阵，是x的3×3×3邻域范围内的的局部平均值，C^adj为矩阵C的共轭矩阵，det表示求方阵的行列式，trace表示求方阵的迹,，若C为零矩阵，那么F取0。利用式(6)求出各点的F值后，在其5×5×5的邻域范围内得到一个严格的极大值点记为特征点。

(2)特征点匹配

在特征点匹配时，采用SURF算法中的特征向量来表征每个特征点，因为在非刚性配准前已经将体数据降采样到了同一尺度空间下，并且已经用刚性配准进行了矫正，因此不需要SURF算法中的尺度空间变换和保证旋转不变性，仅仅借用其特征向量的表达形式。

SURF特征向量的计算方式为：将特征点周围16*16*16的区域分为4个4*4*4子区域，在每个子区域计算式(8)的特征向量

v＝(Σg_X，Σg_Y，Σg_Z，Σ|g_X|，Σ|g_Y|，Σ|g_Z|) (8)

其中g表示各个方向上的梯度分量。在每个子空间形成一个6维的特征向量，在一个特征点处将每个子空间的特征向量排列生成一个4*4*4*6＝384为的特征向量。

因为利用算子提取到的特征点比较多，而且在一个体数据中出现的特征点在另一个体数据中不一定会出现，因此仅仅利用特征点之间特征向量距离最短并不能很好的获得匹配的特征点对。为了解决这个问题，本发明提出了匹配为特征点对的两个特征点r、f应满足以下两个条件：一、参考序列中的特征点r在浮动序列中最邻近的点是特征点f，且特征点f在参考序列中的最邻近的点也是特征点r；二、最邻近特征点与次临近特征点间的距离比小于预设的第一阈值。一些实施例中第一阈值设置为0.6。

用暴力求解确定最小距离的特征点计算量会很大，因此较佳实施例中采用BBF(Best Bin First)算法确定具有最小距离的点(即最邻近的特征点)。BBF算法并不保证找到的点一定是特征向量之间距离最小的点，而是只考察n个距离最近可能性最大的特征点，从而减少了计算量。

图2示出了采用算子提取特征点匹配后的结果，其中子图a和b、子图c和d、子图e和f分别为三个断层的匹配结果。可以看到，最终得到的匹配特征点对会比找到的特征点少很多，每层平均只有一两个左右。

(3)位移场搜索

在非刚性配准中的一个难点是变换矩阵参数多，最优参数搜索过程繁琐。为了使搜索过程简化，避免对梯度的计算，用互信息(MI)结合特征点限制作为相似性测度，利用单点互信息(SMI)与互信息之间的关系，采用邻域搜索的方式进行位移场迭代，从而得到最优位移场。

一些实施例的非刚性配准的相似性测度A可以表示为式(9)

其中ΔU表示位移场U的梯度，x是匹配特征点的坐标，|U(x)-V(x)|²表示在特征点处，迭代得到的位移场U(x)与匹配特征点对之间的位移V(x)之间的二范数。λ和μ为权值。

单点互信息与互信息之间存在以下关系：两个序列图像重合部分各点向SMI增大的方向进行迭代，那么MI就会向增大的方向发展。该关系可通过以下推导得到证明

MI(I，J)＝H(I)+H(J)-H(I，J) (10)

上述关系中p(i)、p(j)、p(i，j)分别表示图像I，J中重合部分分别出现灰度值i和j的边缘概率分布和联合概率分布；x表示图像I，J重合部分的坐标；i(x)，j(x)表示点x处的灰度值。

在式(14)中采用了灰度直方图对灰度概率p(i)、p(j)与p(i，j)进行估计，因为

所以可以看做是函数在所有灰度值对为(i,j)的x点处的函数值和，因此可以得到式(15)。

定义单点互信息

可以得到

式(16)反映了这样一个直觉上的认知——两个图像数据I,J之间的互信息是由他们重合部分的各点共同贡献得到的。

为了体现位移场对应的空间变换U在配准中的作用，将式(16)改写为式(18)所示的形式

可见如果需要对SMI进行估计，必须知道变形后的图像J的灰度分布。

由泰勒公式可以得到

所以当ΔU很小时，有

p(j(U(x)+ΔU))≈p(j(U(x))) (20)

同理可得

p(i(x)，j(U(x)+ΔU))≈p(i(x)，j(U(x))) (21)

这意味着如果每次位移场U变化很小，则可以用上一步迭代得到插值图像的灰度分布来估计SMI。

一些实施例中，位移场搜索的迭代公式可以用式(22)表示

Uⁿ＝U^n-1ο(Id+uⁿ)+uⁿ (22)

其中Uⁿ表示第n次迭代后得到的位移场，Id表示恒等变换，uⁿ表示第n次迭代位移场的变化量，ο表示变换之间的组合。

在较佳实施例中，位移场搜索的具体迭代步骤如下：

a.初始化。利用刚性配准可以将浮动图像插值成与目标图像(即参考图像)大小相同的结果，记为J⁰，目标图像记为T。初始化位移场U^ο＝0。

b.计算位移场变化量uⁿ。这里我们将浮动图像第n次迭代后变换的插值结果记为Jⁿ，目标图像记为T。为了提高收敛速度，并进一步保证配准结果的单应性，采用对称搜索的方式。即，首先计算T中点x与其在Jⁿ中对应位置的3×3×3的邻域内各点x'的SMI(I(x)，Jⁿ(x′))，找出令SMI(x)最大的点x′，记uⁿ ₁(x)＝x′-x，最大的单点互信息为SMI₁；然后计算在Jⁿ中点x在T对应位置的3×3×3的邻域内最大单点互信息SMI₂,记uⁿ ₂(x)＝x-x'，如果SMI₁＞SMI₂，则取uⁿ(x)＝uⁿ ₁(x)，否则取uⁿ(x)＝uⁿ ₂(x)。

c.平滑位移场变化量uⁿ。利用半径为1个体素的高斯核对uⁿ进行平滑。

d.将位移场变化量与位移场相加，更新位移场。即，令Uⁿ⁺¹＝Uⁿ+uⁿ。

e.利用特征点对位移场进行限制。具体包括计算特征点x处的位移场Uⁿ⁺¹(x)，计算该位移场Uⁿ⁺¹(x)与V(x)的差值，该差值小于预设的第二阈值时，使位移场Uⁿ⁺¹(x)＝V(x)，否则Uⁿ⁺¹(x)保持不变。其中V(x)为构成特征点对的两个特征点之间的位移，该特征点对包括特征点x。一些实施例中设置为50mm。特征点的约束可以防止位移场因为自由度过高在uⁿ的临近搜索过程中U的发散。

f.对Uⁿ⁺¹进行加权高斯平滑。具体包括用半径为2个体素的高斯和对Uⁿ⁺¹进行平滑，所述加权高斯平滑在特征点处采用第一权值，在其它位置采用第二权值。为了增强特征点对整个位移场的约束，在特征点处取权值为1.0，其他位置取权值为0.01。加权平滑可用式(23)表示

g.判断是否停止迭代。具体包括计算浮动图像插值结果Jⁿ⁺¹，计算该插值结果Jⁿ⁺¹与参考图像T之间的互信息，当互信息停止增大或迭代次数达到限定次数时停止迭代，将当前位移场作为最优位移场，否则转到步骤b重复上述步骤。

图3示出了浮动图像配准后的插值图像与参考图像相减图。其中子图a-d分别为从头部到腹部的方向四个断层在降采样级别上得到的配准插值图像与目标图像相减的结果。当灰度为0时，在图中显示为灰色，可以看到在相减图中，绝大部分区域的灰度值在0附近，一些重要血管分支在两幅图像中相差不大(从图中相近的黑色与白色区域可以看出)，肺部边缘处比较饱满但仍有没对齐的部分(见子图c、d)。从图3可以直观地看到本发明配准方法的有效性。

步骤S4、用第一配准参数和第二配准参数对第一肺部CT体数据和第二肺部CT体数据中的结节进行匹配，建立对应关系。

因为配准结果是在降采样级别上先进行刚性变换再通过位移场进行非刚性变换，因此将浮动图像空间中的体素坐标映射到参考图像空间的变换方程可表示为

T(X)＝k(ΦX+t+u(X)) (24)

其中X是浮动图像中体素的坐标，Φ是刚性变换中的旋转矩阵，t是位移向量，u(X)是位移场中体素X对应的位移，k是从降采样级别变换为原数据尺度的比例系数，T(X)是浮动图像映射到目标图像中的空间位置。

为了进一步验证本发明方法的效果，进行了验证实验。具体如下：

采用了9个病人不同时间拍摄的肺部CT数据进行实验，每个病人CT数据拍摄的平均间隔时间为10天，断层厚度为1mm，每个体素在同一断层的x和y方向上代表的物理尺寸在0.6～1.0mm之间，具体的实验数据见表1。

表1.实验数据信息

用本发明方法对表1中9个病例进行配准，检测配准结果。因为每组测试数据的灰度范围一致，所以理想情况下浮动图像的插值图像与目标图像的重合部分的灰度差为0，因此可以用两者之间的均方误差(NME)来对配准结果进行度量，NME值越小，说明配准效果越好。NME的计算公如式25所示。

其中x表示浮动图像和目标图像间的重合区域，N是重合区域像素点的个数。具体的度量结果见表2。

表2.NME的统计结果

结节位置的配准偏差。将表2中的病例1中的对应结节的位置进行了手工标记，并将作本发明方法自动配准的结果与手工标记的实际结果进行比较，比较结果见表3。

表3.相同位置结节配准偏差统计

Claims (6)

1.一种肺部CT的三维配准方法，其特征在于，该三维配准方法包括以下步骤：

载入第一肺部CT体数据和第二肺部CT体数据；

以互信息为配准测度，对第一肺部CT体数据和第二肺部CT体数据进行三维刚性配准，使二者中的肺部区域整体对齐，记录X轴、Y轴和Z轴三个方向的位移参数和旋转参数，作为第一配准参数；

在刚性配准的基础上进行非刚性配准，记录非刚性配准的最优位移场，作为第二配准参数；

用第一配准参数和第二配准参数对第一肺部CT体数据和第二肺部CT体数据中的结节进行匹配，建立对应关系；

所述非刚性配准包括特征点提取、特征点匹配和位移场搜索；其中，采用三维的算子提取特征点，特征点匹配中被匹配为特征点对的两个特征点r、f满足以下两个条件：一、参考图像中的特征点r在浮动图像中最邻近的点是特征点f，且特征点f在参考图像中的最邻近的点也是特征点r；二、最邻近特征点与次临近特征点间的距离比小于预设的第一阈值；位移场搜索中，用互信息结合特征点限制作为相似性测度，利用单点互信息与互信息之间的关系，采用邻域搜索的方式进行位移场迭代，从而得到最优位移场，位移场搜索包括以下步骤：

计算位移场变化量；

将位移场变化量与位移场相加，更新位移场；

利用特征点对位移场进行限制，包括计算特征点x处的位移场Uⁿ⁺¹(x)，计算该位移场U^{n +1}(x)与V(x)的差值，该差值小于预设的第二阈值时，使位移场Uⁿ⁺¹(x)＝V(x)，否则Uⁿ⁺¹(x)保持不变，V(x)为构成特征点对的两个特征点之间的位移，该特征点对包括特征点x；

计算浮动图像插值结果，计算该插值结果与参考图像之间的互信息，当互信息停止增大或迭代次数达到限定次数时，将当前位移场作为最优位移场，否则重复上述步骤。

2.根据权利要求1所述的肺部CT的三维配准方法，其特征在于，采用BBF算法确定最邻近的特征点。

3.根据权利要求1所述的肺部CT的三维配准方法，其特征在于，更新位移场前还包括平滑位移场变化量的步骤，在利用特征点对位移场进行限制后还包括对位移场进行加权高斯平滑的步骤，所述加权高斯平滑在特征点处采用第一权值，在其它位置采用第二权值。

4.根据权利要求1所述的肺部CT的三维配准方法，其特征在于，采用对称搜索的方式计算位移场变化量。

5.根据权利要求1所述的肺部CT的三维配准方法，其特征在于，该三维配准方法还包括在进行配准前将第一肺部CT体数据和第二肺部CT体数据每个体素代表的物理尺寸降采样到相同尺度上。

6.根据权利要求1所述的肺部CT的三维配准方法，其特征在于，采用浮动图像与参考图像之间的联合灰度分布直方图计算互信息，其中采用PV插值的方式来计算联合灰度分布直方图。