CN115409879A - 图像配准的数据处理方法及装置、存储介质、电子设备 - Google Patents

Abstract

本申请公开了一种图像配准的数据处理方法及装置、存储介质、电子设备。该方法包括：根据已配准的图像对，生成多个形变场，其中，图像对包括固定图像和浮动图像；根据多个形变场，确定形变场统计特征和物理属性特征；根据配准前后的图像，确定图像相似性特征；根据形变场统计特征，物理属性特征和图像相似性特征，输入预先训练好的分类器，由分类器确定图像对的配准误差。解决了相关技术中图像配准精确度的检测，需要依赖于人工判定的标准图像，存在准确率低的问题。

Description

图像配准的数据处理方法及装置、存储介质、电子设备

技术领域

本申请涉及图像配准领域，具体而言，涉及一种图像配准的数据处理方法及装置、存储介质、电子设备。

背景技术

医学图像配准几乎存在于所有的放疗及辅助诊疗系统中，是将不同模态、不同时间序列的两幅或者多幅图像进行空间匹配的过程。弹性图像配准算法是自适应放疗的关键组成部分，在放射治疗领域的应用越来越多，常运用于通过将图像向解剖图谱配准从而对目标区域进行分割，将不同模态(如CT和MR)的医学图像配准融合以区分肿瘤组织来辅助诊断，在放射治疗时通过图像配准辅助定位照射靶区，提供剂量估计等。

医学图像配准结果的验证是一项很有意义却又复杂的问题，而目前大部分研究都集中在配准算法本身，而不是配准质量的评估，因此设计一种评估配准结果误差的方法是很重要的。配准错误容易发生在图像缺乏对比度的组织区域，在配准算法方面，自由度过多或过少、优化器选用不当、配准模型选用不当都可能导致配准错误。

临床上，图像配准产生的误差会影响医生判断，导致手术中的错误，严重威胁患者的生命健康。现大多配准方法由于缺少金标准，配准结果大部分被移交给专家进行人工判别，常用方法如分屏显示和棋盘格显示进行定性的配准结果可视化，或将结果图像叠加显示直观比较，虽然在高对比度组织结构区域识别配准与否很有效，但耗时费力且不能量化误差大小，具有主观性，易产生误差。配准的质量经常用替代度量去评估，如利用局部解剖地区的重叠度Dice去评估配准质量并不可靠，如果计算Dice 的解剖区域占了很大的体积，则它未必能真实反映组织内部一些细节性的比如气管、血管的配准误差。最常用的去评价弹性图像配准质量的方法是使用大量专家确定的解剖标记点计算目标配准误差，对配准算法精确度进行评估，真实反映局部结构上的配准误差，但这些解剖标记点的标注是一项艰巨的工作，现实中很难获取，尽管存在着解剖标记点的半自动标注方法，但十分耗时，工作量大，不适合日常上使用。此外，可变形的仿体可辅助评估弹性配准算法精确度，但是构建符合复杂临床场景的仿体是很困难的。

针对相关技术中图像配准精确度的检测，需要依赖于人工判定的标准图像，存在准确率低的问题，目前尚未提出有效的解决方案。

发明内容

本申请的主要目的在于提供一种图像配准的数据处理方法及装置、存储介质、电子设备，以解决相关技术中图像配准精确度的检测，需要依赖于人工判定的标准图像，存在准确率低的问题。

为了实现上述目的，根据本申请的一个方面，提供了一种图像配准的数据处理方法，包括：根据已配准的图像对，生成多个形变场，其中，所述图像对包括固定图像和浮动图像；根据多个形变场，确定形变场统计特征和物理属性特征；根据配准前后的图像，确定图像相似性特征；根据所述形变场统计特征，所述物理属性特征和所述图像相似性特征，输入预先训练好的分类器，由所述分类器确定所述图像对的配准误差。

可选的，根据已配准的图像对，生成多个形变场包括：对所述图像对配准过程中的网格参数进行多次扰动，得到多个扰动网格参数；根据每个所述扰动网格参数进行配准计算，生成对应的形变场；通过遍历所述多个扰动网格参数，生成多个形变场。

可选的，根据已配准的图像对，生成多个形变场包括：根据所述图像对配准过程中的网格参数进行配准计算，得到原始形变场；对所述原始形变场的体速进行多次扰动，得到多个形变场。

可选的，根据多个形变场，确定形变场统计特征和物理属性特征包括：根据多个形变场确定所述形变场统计特征，其中，所述形变场统计特征包括下列至少之一：形变场偏差，标准差，偏度，峰度；确定每个形变场的物理属性特征，其中，所述物理属性特征包括下列至少之一：拓扑特征，一致性，物理保真度。

可选的，根据配准前后的图像，确定图像相似性特征包括：根据配准前后的图像进行相似性运算，确定所述图像相似性特征，其中，所述图像相似性特征包括下列至少之一：差分指标，相关性特征，结构相似性特征，模态独立邻域描述子特征。

可选的，根据所述形变场统计特征，所述物理属性特征和所述图像相似性特征，输入预先训练好的分类器，由所述分类器确定所述图像对的配准误差之前，所述方法还包括：对所述形变场统计特征，所述物理属性特征和所述图像相似性特征，通过不同半径的统计窗进行特征扩充，得到多个扩充特征；根据所述形变场统计特征，所述物理属性特征和所述图像相似性特征，输入预先训练好的分类器，由所述分类器确定所述图像对的配准误差包括：将多个所述扩充特征输入所述分类器，由所述分类器确定所述图像对的处于多个配准误差范围的概率；将概率最高的配准误差范围，作为所述图像对的配准误差。

可选的，根据所述形变场统计特征，所述物理属性特征和所述图像相似性特征，输入预先训练好的分类器，由所述分类器确定所述图像对的配准误差之前，所述方法还包括：确定训练样本集中训练样本的形变场统计特征，物理属性特征，图像相似性特征，以及训练样本的目标配准误差；根据所述训练样本的形变场统计特征，物理属性特征，图像相似性特征，以及所述训练样本的目标配准误差，对预设的分类器进行训练，直至所述分类器的分类误差达到满足预设条件。

为了实现上述目的，根据本申请的另一方面，提供了一种图像配准的数据处理装置。该装置包括：生成模块，用于根据已配准的图像对，生成多个形变场，其中，所述图像对包括固定图像和浮动图像；第一确定模块，用于根据多个形变场，确定形变场统计特征和物理属性特征；第二确定模块，用于根据配准前后的图像，确定图像相似性特征；分类模块，用于根据所述形变场统计特征，所述物理属性特征和所述图像相似性特征，输入预先训练好的分类器，由所述分类器确定所述图像对的配准误差。

根据本申请的另一方面，还提供了一种计算机可读存储介质，所述处存储介质用于存储程序，其中，所述程序执行上述中任意一项所述的图像配准的数据处理方法。

根据本申请的另一方面，还提供了一种电子设备，包括一个或多个处理器和存储器，所述存储器用于存储一个或多个程序，其中，当所述一个或多个程序被所述一个或多个处理器执行时，使得所述一个或多个处理器实现上述中任意一项所述的图像配准的数据处理方法。

通过本申请的上述步骤，根据已配准的图像对，生成多个形变场，其中，所述图像对包括固定图像和浮动图像；根据多个形变场，确定形变场统计特征和物理属性特征；根据配准前后的图像，确定图像相似性特征；根据所述形变场统计特征，所述物理属性特征和所述图像相似性特征，输入预先训练好的分类器，由所述分类器确定所述图像对的配准误差。

通过对已经配准的固定图像和浮动图像配准过程中进行扰动，生成多个形变场，根据多个形变场确定已配准的图像对的形变场统计特征，物理属性特征和图像相似性特征，来根据分类器对图像对的配准误差进行检测，实现了提高对已配准图像配准误差确定的准确性和效率的技术效果，进而解决了相关技术中图像配准精确度的检测，需要依赖于人工判定的标准图像，存在准确率低的问题。

附图说明

构成本申请的一部分的附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。在附图中：

图1是根据本申请实施例提供的一种图像配准的数据处理方法的流程图；

图2是根据本申请实施方式提供的图像配准方法的流程图；

图3是根据本申请实施方式提供的对网格参数进行随机扰动的示意图；

图4是根据本申请实施方式提供的生成多个形变场的两种不同的方式的流程图；

图5(a)是根据本申请实施方式提供的原始形变场的示意图；

图5(b)是根据本申请实施方式提供的[-2,2]mm尺度扰动下得到的形变场的示意图；

图5(c)是根据本申请实施方式提供的[-3,3]mm尺度扰动下得到的形变场的示意图；

图5(d)是根据本申请实施方式提供的[-5,5]mm尺度扰动下得到的形变场的示意图；

图6是根据本申请实施方式提供的目标配准误差的计算流程的示意图；

图7是根据本申请实施例提供的一种图像配准的数据处理装置的示意图；

图8是根据本申请实施例提供的一种电子设备的示意图。

具体实施方式

需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本申请。

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

需要说明的是，本申请的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别类似的对象，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本申请的实施例。此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。

需要说明的是，本公开所涉及的相关信息(包括但不限于用户设备信息、用户个人信息等)和数据(包括但不限于用于展示的数据、分析的数据等)，均为经用户授权或者经过各方充分授权的信息和数据。例如，本系统和相关用户或机构间设置有接口，在获取相关信息之前，需要通过接口向前述的用户或机构发送获取请求，并在接收到前述的用户或机构反馈的同意信息后，获取相关信息。

实施例

下面结合优选的实施步骤对本发明进行说明，图1是根据本申请实施例提供的一种图像配准的数据处理方法的流程图，如图1所示，该方法包括如下步骤：

步骤S101，根据已配准的图像对，生成多个形变场，其中，图像对包括固定图像和浮动图像；

步骤S102，根据多个形变场，确定形变场统计特征和物理属性特征；

步骤S103，根据配准前后的图像，确定图像相似性特征；

步骤S104，根据形变场统计特征，物理属性特征和图像相似性特征，输入预先训练好的分类器，由分类器确定图像对的配准误差。

通过上述步骤，通过对已经配准的固定图像和浮动图像配准过程中进行扰动，生成多个形变场，根据多个形变场确定已配准的图像对的形变场统计特征，物理属性特征和图像相似性特征，来根据分类器对图像对的配准误差进行检测，实现了提高对已配准图像配准误差确定的准确性和效率的技术效果，进而解决了相关技术中图像配准精确度的检测，需要依赖于人工判定的标准图像，存在准确率低的问题。

上述步骤的执行主体可以为计算器，处理器，服务器等具有数据处理运算能力的设备。上述计算器，处理器，服务器可以与图像配准的设备具有数据连接，在图像配准设备对图像进行配准后，得到上述图像对，将上述图像对发送给上述执行主体，由执行主体进行配准误差的确定。因此，上述执行主体也可以时具有运算和数据处理能力的功能模块，该功能模块可以设置在某一硬件设备内部，在需要时，响应调用请求提供上述配准误差的服务。

上述步骤S101中，根据已配准的图像对，生成多个形变场。已配准的图像对包括固定图像和浮动图像，是经过图像配准技术处理后得到的。相关技术中图像配准精确度的检测，需要依赖于人工判定的标准图像，也即是金标准，其存在准确率低的问题。

图像配准问题通常被转化为寻找最优空间变换的过程，即一个损失函数的最优化问题。衡量两幅图像之间差异性或相似性的损失函数，与图像的最优的形变场，和图像的变形过程有关。配准是个迭代优化的问题，并且可以由原始配准形变场初始化。好的配准算法可以看作一个稳定的系统，如果施加随机形变扰动，稳定系统是可以经过迭代优化回到稳定状态。

本实施例提供了一种基于B样条构建自由形变模型进行非刚性图像配准，B样条函数是非刚性图像配准中常用的技术手段，具有最小局部支撑特性，B样条曲线是一系列控制点的加权和，控制点对B样条曲线的影响是局部的，曲线上任意点的位置可由其周围相邻的控制点的位置计算得到。三维图像中每个体素的位移矢量都可以由其周围的控制点的位移来表示，即通过某一点周围的控制点可以插值得到该点的位移矢量。

配准的优化过程中，图像同质区域的部分可能会缺乏图像强度梯度，因而不能很好驱动损失函数的优化，造成损失函数的最小值的取值范围较宽，达不到全局最优。即使在非同质区域，配准的错误也会由于优化迭代次数不足，或损失函数处于局部最优发生。在损失测度函数处于局部最优的情况，或者图像中错位的结构和同质区域初始状态的微小变化将会很大影响配准的最终解，施加一定随机扰动作为微小形变，可能反而改善了最终的配准结果，该配准可以看作是不稳定的，结果受不同初始状态的影响很大，在这些区域产生的配准结果可以认为是不可重现的。

如果配准产生了完全正确的形变场，并且图像拥有足够的结构信息，使得相似性损失测度找到了全局最优解，配准产生的形变场结果是可以复制的，施加一定额外的变形就会使测度函数变得更差。因此施加一定范围随机形变扰动经配准得到数个形变场，利用这些形变场之间的统计特性判别配准的结果能否重现，可以评估配准的几何不确定性，形变场之间的差异越小代表着配准越稳定，结果可重现，配准误差越小，进而表征误差。

因此，本实施例，通过对已配准的图像对，进行多次扰动，生成多个形变场，来提取图像特征，进行配准前图像和配准后图像的比对，进而确定出配准误差。理论上，配准精度越高，其抗干扰的能力越强，也即是经过多次干扰后，其配准前图像和配后图像的差别仍然不大。从而根据最终配准前图像和配准后图像的相似度来确定，图像配准的误差等级，也即是误差数值范围。

可选的，根据已配准的图像对，生成多个形变场包括：对图像对配准过程中的网格参数进行多次扰动，得到多个扰动网格参数；根据每个扰动网格参数进行配准计算，生成对应的形变场；通过遍历多个扰动网格参数，生成多个形变场。

需要说明的是，上述根据网格参数生成网格，然后根据网格去进行配准是图像配准的流程。本实施例的一种实施方式，可以根据对图像对配准过程中的网格参数进行多次扰动，得到多个扰动网格参数。然后，根据每个扰动网格参数进行配准计算，生成对应的形变场。通过遍历多个扰动网格参数，生成多个形变场。

如图4上部的流程所示，提供了第一种方案：给定一种B样条网格点参数(也即是上述网格参数)后，对其施加N次随机形变扰动后，得到N个初始参数T^P，作为初始状态进行配准，得到N个形变场T^P-R，根据以下公式对形变场位移量大小进行统计，观察位移量的数值分布，表征医学图像弹性配准误差。由此提出了形变场偏差 bias(T^P-R)，标准差std(T^P-R)，偏度Skewness(T^P-R)，峰度Kurtosis(T^P-R)作为特征。

可选的，根据已配准的图像对，生成多个形变场包括：根据图像对配准过程中的网格参数进行配准计算，得到原始形变场；对原始形变场的体速进行多次扰动，得到多个形变场。

需要说明的是，上述根据网格参数生成网格，然后根据网格去进行配准得到形变场是图像配准的流程。本实施例的另一种实施方式，可以根据图像对配准过程中的网格参数进行配准计算，得到原始形变场。然后对原始形变场的体速进行多次扰动，得到多个形变场。

如图4下部的流程所示，提供了第二种方案：不改变B样条网格点初始参数，而是先进行一次配准，得到原始形变场T^R并保存下来，在TR基础上对每个体素施加随机形变扰动，重复N次得到N个形变场T^R-P，受到扰动的形变场T^R-P作为下面配准的起始点，继续配准N次得到N个形变场T^R-P-R，实际上，得到的T^R-P-R应该与原始形变场T^R相似但不完全相同，因此不同T^R-P-R偏移量之间的分布是否稳定，T^R-P-R与原始形变场T^R之间的差异都逐体素反映了配准结果的不确定性。

可选的，根据多个形变场，确定形变场统计特征和物理属性特征包括：根据多个形变场确定形变场统计特征，其中，形变场统计特征包括下列至少之一：形变场偏差，标准差，偏度，峰度；确定每个形变场的物理属性特征，其中，物理属性特征包括下列至少之一：拓扑特征，一致性，物理保真度。

上述形变场统计特征，也即是形变场的不同参数的统计特征，例如，形变场偏差，标准差，偏度，峰度等。在上述第一种方案下，经形变扰动得到的N个形变场T^P-R的均值为

与T^R的偏差为bias(T^P-R)，T^P-R的标准差std(T^P-R)，std(T^P-R)和 bias(T^P-R)表达了形变场位移量偏离均值的程度，反映了位移量彼此之间波动的程度。一个完美的配准算法在施加随机形变扰动之后产生的形变场应该是一样的，理论上 std(T^P-R)和bias(T^P-R)都为0，当其越大时意味着随机形变扰动后配准得到的形变场位移量之间差异较大，配准优化过程愈发不稳定。偏度为Skewness(T^P-R)，衡量形变场位移量分布的偏斜方向和程度，表征围绕平均值的位移量分布的不对称性，偏度接近0则可认为分布是对称的，不同形变场位移量分布均匀，配准算法是稳定的。峰度为Kurtosis(T^P-R)表征了形变场位移量分布形态陡缓程度，峰度的绝对值数值越大代表不同形变场位移量分布的陡缓程度与正态分布之间差异越大。理论上Skewness(T^P-R) 和Kurtosis(T^P-R)都为0，当其值越大意味着配准结果不是全局最优解，即配准误差也越大，std(T^P-R)、bias(T^P-R)和Kurtosis(T^P-R)较大，配准误差较大。

在上述第二种方案下，经扰动、配准得到的N个形变场T^R-P-R的均值

由此求得标准差std(T^R-P-R)以及与T^R的偏差bias(TR-P-R)，偏度为Skewness(T^R-P-R)，峰度为Kurtosis(T^R-P-R)，std(T^R-P-R)和bias(T^R-P-R)与配准误差正相关，td(T^R-P-R) 和bias′值较大，配准误差较大。

上述物理属性特征可以为形变场的物理属性特征，包括拓扑特征，一致性，物理保真度。上述拓扑特征可以包：雅可比行列式JAC，微分同胚是形变场在物理上可行的必要条件，即图像中器官组织在形变场作用下只能被压缩和舒展，不能经历不可逆的空间变换，该性质与形变场的雅可比行列式有关。雅可比行列式表示了配准前后的局部体积变化，是评估配准形变场拓扑特性的一种常用指标。负值区域表示物理上无法实现的变形如结构折叠，当其值等于1时，在该点处的体积没有发生变化；大于1 时，体积发生膨胀；在0-1之间时，发生缩小。如图所示，红框中Jac值小于0，意味着形变场在该点出现奇点，导致相应器官组织发生了折叠，这是不符合实际情况的。

上述一致性可以包括：逆一致性ICE，一个好的配准算法应该具有良好的逆一致性(ICE)，逆一致性误差表征了从浮动图像IM到固定图像IF的形变场

与从固定图像IF到浮动图像IM的形变场

之间的差异，若配准算法完美，则利用该算法将固定图像IF向浮动图像IM进行配准得到

，即为

的逆。正反向配准得到的形变场经合成，得到合成形变场，ICE定义为浮动图像上某点在合成形变场下的位移，如果逆一致性完美，合成的形变场应为0，在该点的位移为0，ICE也为0。

上述物理保真度可以包括：空间平滑度SS，弹性配准中形变场的拓扑特性十分重要，深度学习中常用形变场的正则化对形变场平滑度进行约束，因此在这使用空间平滑度SS来表征形变场的物理保真度，谐波能量HE与形变场的平滑度成反比，二者数学原理近乎一致，所以在这只使用了空间平滑度，形变场相邻体素方向的变化越大，则在该区域往往产生非平滑的表现。

可选的，根据配准前后的图像，确定图像相似性特征包括：根据配准前后的图像进行相似性运算，确定图像相似性特征，其中，图像相似性特征包括下列至少之一：差分指标，相关性特征，结构相似性特征，模态独立邻域描述子特征。

上述差分指标可以包括配准完成后得到的结果图像与固定图像间的差分σ，如配准算法完美，理论上差分σ应为0。还可以包括：图像强度的平方差之和(SSD)，也可以作为配准效果的评价指标，SSD越低，两幅图像越相似。还可以包括：归一化互信息NMI：NMI常作为多模图像之间配准的测度，当两幅图像相似度越高或重合部分越大时，其相关性也越大，归一化互信息(NMI)越大。

上述模态独立邻域描述子特征可以包括：模态独立邻域描述子(MIND)，MIND是一个描述每个体素周围局部模式的特征描述子，通过观察中心体素和固定距离外的体素之间的相似性来计算定义，图像中一个体素周围的局部模式也应该是相似的，计算配准前后图像的MIND特征图，若配准效果越好，则结果图像与固定图像的MIND特征应该越相似，因此可以拿MIND之间的差异进行配准误差的表征。

上述相关性特征可以包括：互相关系数NCC：NCC常作为单模图像配准之间的测度，评判结果图像与固定图像之间的相似程度，因此可以使用NCC对配准误差进行表征，配准效果越好，NCC应该越大，逼近1。

上述结构相似性特征可以包括：结构相似性指数SSIM：一种用于量化两幅图像间的结构相似性的指标，考虑了人眼的视觉特性，从亮度、对比度以及结构量化图像的属性，SSIM值越大代表图像越相似。还可以包括：峰值信噪比PSNR：常用于衡量两张图像之间差异，是一种评价图像的客观标准，配准结果图像与固定图像越相似，则 PSNR值越大。

可选的，根据形变场统计特征，物理属性特征和图像相似性特征，输入预先训练好的分类器，由分类器确定图像对的配准误差之前，方法还包括：对形变场统计特征，物理属性特征和图像相似性特征，通过不同半径的统计窗进行特征扩充，得到多个扩充特征；根据形变场统计特征，物理属性特征和图像相似性特征，输入预先训练好的分类器，由分类器确定图像对的配准误差包括：将多个扩充特征输入分类器，由分类器确定图像对的处于多个配准误差范围的概率；将概率最高的配准误差范围，作为图像对的配准误差。

上述特征扩充可以获得更多的特征，作为训练样本对分类器进行训练。上述特征都是在图像上每个体素点计算的，为了减少特征的不连续性并改善与其他特征的交互，通过使用不同半径的窗(统计窗半径：[1,2,3,4,5,6,7,8,9]voxel)计算窗内每个点特征的平均值和最大值，从而增加特征的总集合。一个点上的特征可以通过不同大小的窗池化成18个特征，9个来自最大值，9个来自平均值。

可选的，根据形变场统计特征，物理属性特征和图像相似性特征，输入预先训练好的分类器，由分类器确定图像对的配准误差之前，方法还包括：确定训练样本集中训练样本的形变场统计特征，物理属性特征，图像相似性特征，以及训练样本的目标配准误差；根据训练样本的形变场统计特征，物理属性特征，图像相似性特征，以及训练样本的目标配准误差，对预设的分类器进行训练，直至分类器的分类误差达到满足预设条件。

在对上述分类器进行训练之前，还需要对各个图像对的特征数据对应的配准误差进行确定，作为样本标签，对分类器进行训练。本实施例中根据配准前后对应标记点的欧式距离，作为医学图像弹性配准的真实误差，也即是目标配准误差TRE，也称真实配准误差。

然后将上述特征数据和对应的样本标签也即是上述目标配准误差，输入分类器进行训练，直至分类误差满足条件。当然，可以将多个样本分为训练样本和测试样本，在训练完成后，通过测试样本进行测试，确定分类器的分类误差。

需要说明的是，在附图的流程图示出的步骤可以在诸如一组计算机可执行指令的计算机系统中执行，并且，虽然在流程图中示出了逻辑顺序，但是在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤。

需要说明的是，本申请还提供了一种可选的实施方式，下面对该实施方式进行详细说明。

本实施方式提供了一种医学图像配准精度自动评估方法，对图像弹性配准算法的精确度进行自动定量评估，减少分析时间并提供配准不确定度信息定量分析误差，帮助医生有效识别关键区域的不真实形变，做出准确的临床判断，保证患者安全。关键点在于：医学图像配准自动评估流程：基于形变场统计特征、物理属性特征和图像相似性特征结合的医学图像配准误差自动评估；基于形变场两次扰动方式的变场统计特征计算方法；适用于肺部呼吸运动的图像形变场多尺度扰动方法。

首先对本实施方式涉及到的部分字母的含义进行说明，具体如表1所示，表1为公式字母含义表。

表1公式字母含义表

图2是根据本申请实施方式提供的图像配准方法的流程图，如图2所示，该系统由固定图像IF和浮动图像IM输入，经配准和形变场扰动后，通过计算生成三大类特征，包括形变场统计特征、物理属性特征和图像相似性特征，然后所有特征进行池化扩充，输入随机森林回归器中训练，预测医学图像弹性配准误差，根据TG-132标准，将回归得出的误差值进行三分类并在图像中用伪彩色图显示，更易于医师参考判断。

所提出的系统以㎜级别逐点预测医学图像目标配准误差TRE，并在独立拥有对应标记点的数据集上进行方法的验证，该系统主要包括四部分，详细如下。

医学图像弹性配准问题通常被转化为寻找最优空间变换的过程，即一个损失函数的最优化问题。可以使用IM表示待配准的浮动图像，IF表示参考固定图像，T为浮动图像到固定图像之间的形变场，则医学图像弹性配准问题可以用公式表示成：

其中loss为衡量两幅图像之间差异性或相似性的损失函数，

表示最优的形变场，

表示图像的变形过程。

配准是个迭代优化的问题，并且可以由原始配准形变场T初始化。好的配准算法可以看作一个稳定的系统，如果施加随机形变扰动，稳定系统是可以经过迭代优化回到稳定状态。

(3)形变场统计特征(DVFs’Statisticalfeatures)

图3是根据本申请实施方式提供的对网格参数进行随机扰动的示意图，施加随机形变扰动如图3所示，随机形变扰动是针对B样条控制点参数{c₁,...,c_3M}来说的，也即是上述网格参数，其中M是B样条网格点的数量，3表示空间维度，根据控制点网格的尺寸，随机生成一定范围下的系数{p₁,...,p_3M}，使其与原先的参数相加，得到新的修改系数{c₁+p₁,...,c_3M+p_3M}，作为新的初始状态，之后进行配准计算多个形变场。

图4是根据本申请实施方式提供的生成多个形变场的两种不同的方式的流程图，如图4所示，共有两种不同的方式对配准进行随机形变扰动，产生多个形变场计算统计特征。

第一种方案：给定一种B样条网格点参数后，对其施加N次随机形变扰动后，得到N个初始参数T^P，作为初始状态进行配准，得到N个形变场T^P-R，根据以下公式对形变场位移量大小进行统计，观察位移量的数值分布，表征医学图像弹性配准误差。由此提出了形变场偏差bias(T^P-R)，形变场标准差std(T^P-R)，形变场偏度 Skewness(T^P-R)，形变场峰度Kurtosis(T^P-R)作为特征。

经形变扰动得到的N个形变场T^P-R的均值为

与T^R的偏差为bias(T^P-R)， T^P-R的标准差std(T^P-R)，std(T^P-R)和bias(T^P-R)表达了形变场位移量偏离均值的程度，反映了位移量彼此之间波动的程度。一个完美的配准算法在施加随机形变扰动之后产生的形变场应该是一样的，理论上std(T^P-R)和bias(T^P-R)都为0，当其越大时意味着随机形变扰动后配准得到的形变场位移量之间差异较大，配准优化过程愈发不稳定。偏度为Skewness(T^P-R)，衡量形变场位移量分布的偏斜方向和程度，表征围绕平均值的位移量分布的不对称性，偏度接近0则可认为分布是对称的，不同形变场位移量分布均匀，配准算法是稳定的。峰度为Kurtosis(T^P-R)表征了形变场位移量分布形态陡缓程度，峰度的绝对值数值越大代表不同形变场位移量分布的陡缓程度与正态分布之间差异越大。理论上Skewness(T^P-R)和Kurtosis(T^P-R)都为0，当其值越大意味着配准结果不是全局最优解，即配准误差也越大。std(T^P-R)、bias(T^P-R)和Kurtosis(T^P-R)越大，配准误差越大。

第二种方案：另一种方式不改变B样条网格点初始参数，而是先进行一次配准，得到原始形变场T^R并保存下来，在T^R基础上对每个体素施加随机形变扰动，重复N次得到N个形变场T^R-P，受到扰动的形变场T^R-P作为下面配准的起始点，继续配准N 次得到N个形变场T^{R -P-R}，实际上，得到的T^R-P-R应该与原始形变场T^R相似但不完全相同，因此不同T^R-P-R偏移量之间的分布是否稳定，T^R-P-R与原始形变场T^R之间的差异都逐体素反映了配准结果的不确定性。

根据以下公式计算统计特征，经扰动、配准得到的N个形变场T^R-P-R的均值

由此求得标准差std(T^R-P-R)以及与T^R的偏差bias(T^R-P-R)，偏度为 Skewness(T^R-P-R)，峰度为Kurtosis(T^R-P-R)，std(T^R-P-R)和bias(T^R-P-R)与配准误差正相关，std(T^R-P-R)和bias′值越大，配准误差越大。

图4中示出的已配准的图像对是肺部的CT图像，考虑到肺组织的不同区域在呼吸等运动过程时发生的形变幅度是不一样的，下肺区域形变较大，而肺尖运动幅度较小，靠近胸膜的肺部区域的TRE是近纵隔区域的三倍。若对B样条网格点参数只使用同一范围内的扰动形变，统计特征将只能良好地对对应范围内的配准误差进行表征，而这是不符合上述先验知识的，不能实现对可能产生的配准误差的全覆盖，所以为提升误差预测的适用性，在N次扰动中使用不同的尺度的随机形变，扰动范围分别为 [-2，2]，[-3，3]，[-5，5]mm，比例为2∶2∶1，不同扰动范围下得到的形变场如图5(a-d)所示，图5(a)是根据本申请实施方式提供的原始形变场的示意图；图5(b)是根据本申请实施方式提供的[-2，2]mm尺度扰动下得到的形变场的示意图；图5(c)是根据本申请实施方式提供的[-3，3]mm尺度扰动下得到的形变场的示意图；图5(d)是根据本申请实施方式提供的[-5，5]mm尺度扰动下得到的形变场的示意图。

(4)物理属性特征(DVF’s physiologically realistic features)

由于人体组织结构在发生形态变化时是具有一定规律的，所以医学图像上相邻像素点的位移在大小和方向上都是相近的，将两个不同呼吸相位的图像进行配准，产生的形变场理论上是连续且平滑的，所以有必要通过配准后得到的形变场的物理保真度来判断配准结果的好坏，接下来对形变场物理属性进行特征定义。

雅可比行列式JAC：微分同胚是形变场在物理上可行的必要条件，即图像中器官组织在形变场作用下只能被压缩和舒展，不能经历不可逆的空间变换，该性质与形变场的雅可比行列式有关。雅可比行列式表示了配准前后的局部体积变化，是评估配准形变场拓扑特性的一种常用指标。负值区域表示物理上无法实现的变形如结构折叠，当其值等于1时，在该点处的体积没有发生变化；大于1时，体积发生膨胀；在0-1 之间时，发生缩小。Jac值小于0，意味着形变场在该点出现奇点，导致相应器官组织发生了折叠，这是不符合实际情况的。

式中，J(i′，j′，k′)代表坐标位于(i′，j′，k′)的像素的雅克比矩阵，det()代表行列式；此处x′，y′，z′表示图像的xyz轴方向。

逆一致性ICE：一个好的配准算法应该具有良好的逆一致性(ICE)，逆一致性误差表征了从浮动图像IM到固定图像IF的形变场

与从固定图像IF到浮动图像IM 的形变场

之间的差异，若配准算法完美，则利用该算法将固定图像IF向浮动图像IM进行配准得到

即为

的逆。

正反向配准得到的形变场经合成，得到合成形变场，ICE定义为浮动图像上某点在合成形变场下的位移，如果逆一致性完美，合成的形变场应为0，在该点的位移为0， ICE也为0，点位移量越大，逆一致性越不好，配准误差越大。

空间平滑度SS：弹性配准中形变场的拓扑特性十分重要，深度学习中常用形变场的正则化对形变场平滑度进行约束，因此在这使用空间平滑度SS来表征形变场的物理保真度，谐波能量HE与形变场的平滑度成反比，二者数学原理近乎一致，所以在这只使用了空间平滑度，公式如下，形变场相邻体素方向的变化越大，则在该区域往往产生非平滑的表现。

式中，

表示梯度。

(5)图像相似性特征(Imagesimilarityfeatures)

通过比较配准前后的图像的相似性，可以直观地判断出图像的轮廓以及组织结构是否配准，接下来对配准前后图像相似性进行特征定义。

结果图像与固定图像的差分σ：配准完成后得到的结果图像与固定图像间的差分σ，如配准算法完美，理论上差分σ应为0，差分图像中配准越精确的部分纹理越不明显。此外图像强度的平方差之和(SSD)也常作为配准效果的评价指标，SSD越低，两幅图像越相似，配准算法效果越好，σ与SSD在数学原理上近似。

其中

表示图像上的像素点，

表示图像的变形过程，Ω为图像空间。

归一化互信息NMI：NMI常作为多模图像之间配准的测度，当两幅图像相似度越高或重合部分越大时，其相关性也越大，归一化互信息(NMI)越大，因此使用NMI 进行配准误差表征。NMI特征在半径为[1，2，3，4，5，6，7，8，9]voxel的图像patch上进行计算，一个点上可以得到9个NMI特征。

其中H()表示()内图像的熵。

模态独立邻域描述子(MIND)：MIND是一个描述每个体素周围局部模式的特征描述子，通过观察中心体素和固定距离外的体素之间的相似性来计算定义，图像中一个体素周围的局部模式也应该是相似的，计算配准前后图像的MIND特征图，若配准效果越好，则结果图像与固定图像的MIND特征应该越相似，因此可以拿MIND之间的差异进行配准误差的表征。

互相关系数NCC：NCC常作为单模图像配准之间的测度，评判结果图像与固定图像之间的相似程度，因此可以使用NCC对配准误差进行表征，配准效果越好，NCC 应该越大，逼近1。

式中，xi为结果图像的第j个像素点，Ω_F为固定图像的图像空间。

结构相似性指数SSIM：一种用于量化两幅图像间的结构相似性的指标，考虑了人眼的视觉特性，从亮度、对比度以及结构量化图像的属性，SSIM值越大代表图像越相似。

峰值信噪比PSNR：常用于衡量两张图像之间差异，是一种评价图像的客观标准，配准结果图像与固定图像越相似，则PSNR值越大。

(6)特征扩充

形变场偏差bias、形变场标准差std、形变场偏度Skewness、形变场峰度Kurtosis、雅可比行列式Jac、结果图像与固定图像的差分σ、逆一致性ICE、模态独立邻域描述子MIND、结构相似性指数SSIM特征都是在图像上每个体素点计算的，为了减少特征的不连续性并改善与其他特征的交互，通过使用不同半径的窗(统计窗半径： [1，2，3，4，5，6，7，8，9]voxel)计算窗内每个点特征的平均值和最大值，从而增加特征的总集合。一个点上的特征可以通过不同大小的窗池化成18个特征，9个来自最大值，9个来自平均值。

平均值

最大值

归一化互信息NMI、互相关系数NCC、峰值信噪比PSNR、空间平滑度SS特征则在[1，2，3，4，5，6，7，8，9]voxel的图像的子图像块patch上直接进行计算，一个点上可以得到9个NMI特征。相关图像信息不能在整个图像行进行计算，一般在子图像块上计算；在不同大小的子图像块上计算，得到不同的结果，也起到了数据扩充的效果。

回归标签值计算：

回归器的标签值使用目标配准误差TRE，TRE定义为配准前后对应标记点的欧式距离，作为医学图像弹性配准的真实误差，如图6所示。图6是根据本申请实施方式提供的目标配准误差的计算流程的示意图

TRE＝||T(x_F)-x_M||₂

其中T为形变场，X_F为固定图像上的标记点，X_M为浮动图像上对应的标记点。

(7)随机森林回归

使用随机森林进行回归预测。随机森林是由多颗决策树用随机的方式结合在一起，每一颗决策树之间是没有关联的，降低了错误率，训练速度比较快，能够处理很高维度的数据，即特征数量很多，并且不用做特征选择。随机森林具有计算每个特征重要性的能力，只需进行少量的额外计算，就可以显示每个特征对森林的贡献度。

随机森林由上文提及的三大类特征输入进行训练，预测得到的配准误差是mm级别的，根据TG-132标准，配准TRE从0-3mm之间被认为是结果较好，3-6mm之间被认为是结果一般，大于6mm被认为是结果较差。将真实TRE与预测得到的TRE根据 TG-132和回归结果分为三类。

本申请实施例还提供了一种图像配准的数据处理装置，需要说明的是，本申请实施例的图像配准的数据处理装置可以用于执行本申请实施例所提供的用于图像配准的数据处理方法。以下对本申请实施例提供的图像配准的数据处理装置进行介绍。

图7是根据本申请实施例提供的一种图像配准的数据处理装置的示意图，如图7所示，该装置包括：生成模块72，第一确定模块74，第二确定模块76，分类模块78，下面对该装置进行详细说明。

生成模块72，用于根据已配准的图像对，生成多个形变场，其中，所述图像对包括固定图像和浮动图像；第一确定模块74，与上述生成模块72相连，用于根据多个形变场，确定形变场统计特征和物理属性特征；第二确定模块76，与上述第一确定模块74相连，用于根据配准前后的图像，确定图像相似性特征；分类模块78，与上述第二确定模块76相连，用于根据所述形变场统计特征，所述物理属性特征和所述图像相似性特征，输入预先训练好的分类器，由所述分类器确定所述图像对的配准误差。

本申请实施例提供的图像配准的数据处理装置，通过对已经配准的固定图像和浮动图像配准过程中进行扰动，生成多个形变场，根据多个形变场确定已配准的图像对的形变场统计特征，物理属性特征和图像相似性特征，来根据分类器对图像对的配准误差进行检测，实现了提高对已配准图像配准误差确定的准确性和效率的技术效果，进而解决了相关技术中图像配准精确度的检测，需要依赖于人工判定的标准图像，存在准确率低的问题。

图像配准的数据处理装置包括处理器和存储器，上述生成模块72，第一确定模块74，第二确定模块76，分类模块78等均作为程序单元存储在存储器中，由处理器执行存储在存储器中的上述程序单元来实现相应的功能。

处理器中包含内核，由内核去存储器中调取相应的程序单元。内核可以设置一个或以上，通过调整内核参数来解决了相关技术中图像配准精确度的检测，需要依赖于人工判定的标准图像，存在准确率低的问题。

存储器可能包括计算机可读介质中的非永久性存储器，随机存取存储器(RAM)和/或非易失性内存等形式，如只读存储器(ROM)或闪存(flash RAM)，存储器包括至少一个存储芯片。

本发明实施例提供了一种计算机可读存储介质，其上存储有程序，该程序被处理器执行时实现图像配准的数据处理方法。

本发明实施例提供了一种处理器，处理器用于运行程序，其中，程序运行时执行图像配准的数据处理方法。

图8是根据本申请实施例提供的一种电子设备的示意图，如图8所示，本申请实施例提供了一种电子设备80，设备包括处理器、存储器及存储在存储器上并可在处理器上运行的程序，处理器执行程序时实现上述任一项方法的步骤。

本申请中的设备可以是服务器、PC、PAD、手机等。

本申请还提供了一种计算机程序产品，当在图像配准的数据处理设备上执行时，适于执行初始化有上述任一方法步骤的程序。

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

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

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

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

在一个典型的配置中，计算设备包括一个或多个处理器(CPU)、输入/输出接口、网络接口和内存。

存储器可能包括计算机可读介质中的非永久性存储器，随机存取存储器(RAM)和/或非易失性内存等形式，如只读存储器(ROM)或闪存(flash RAM)。存储器是计算机可读介质的示例。

计算机可读介质包括永久性和非永久性、可移动和非可移动媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。计算机的存储介质的例子包括，但不限于相变内存(PRAM)、静态随机存取存储器(SRAM)、动态随机存取存储器(DRAM)、其他类型的随机存取存储器(RAM)、只读存储器(ROM)、电可擦除可编程只读存储器(EEPROM)、快闪记忆体或其他内存技术、只读光盘只读存储器(CD-ROM)、数字多功能光盘(DVD)或其他光学存储、磁盒式磁带，磁带磁盘存储或其他磁性存储设备或任何其他非传输介质，可用于存储可以被计算设备访问的信息。按照本文中的界定，计算机可读介质不包括暂存电脑可读媒体(transitory media)，如调制的数据信号和载波。

还需要说明的是，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、商品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、商品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括要素的过程、方法、商品或者设备中还存在另外的相同要素。

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

以上仅为本申请的实施例而已，并不用于限制本申请。对于本领域技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本申请的权利要求范围之内。

Claims (10)

1.一种图像配准的数据处理方法，其特征在于，包括：

根据已配准的图像对，生成多个形变场，其中，所述图像对包括固定图像和浮动图像；

根据多个形变场，确定形变场统计特征和物理属性特征；

根据配准前后的图像，确定图像相似性特征；

根据所述形变场统计特征，所述物理属性特征和所述图像相似性特征，输入预先训练好的分类器，由所述分类器确定所述图像对的配准误差。

2.根据权利要求1所述的方法，其特征在于，根据已配准的图像对，生成多个形变场包括：

对所述图像对配准过程中的网格参数进行多次扰动，得到多个扰动网格参数；

根据每个所述扰动网格参数进行配准计算，生成对应的形变场；

通过遍历所述多个扰动网格参数，生成多个形变场。

3.根据权利要求1所述的方法，其特征在于，根据已配准的图像对，生成多个形变场包括：

根据所述图像对配准过程中的网格参数进行配准计算，得到原始形变场；

对所述原始形变场的体速进行多次扰动，得到多个形变场。

4.根据权利要求2或3所述的方法，其特征在于，根据多个形变场，确定形变场统计特征和物理属性特征包括：

根据多个形变场确定所述形变场统计特征，其中，所述形变场统计特征包括下列至少之一：形变场偏差，标准差，偏度，峰度；

确定每个形变场的物理属性特征，其中，所述物理属性特征包括下列至少之一：拓扑特征，一致性，物理保真度。

5.根据权利要求1所述的方法，其特征在于，根据配准前后的图像，确定图像相似性特征包括：

根据配准前后的图像进行相似性运算，确定所述图像相似性特征，其中，所述图像相似性特征包括下列至少之一：差分指标，相关性特征，结构相似性特征，模态独立邻域描述子特征。

6.根据权利要求1所述的方法，其特征在于，根据所述形变场统计特征，所述物理属性特征和所述图像相似性特征，输入预先训练好的分类器，由所述分类器确定所述图像对的配准误差之前，所述方法还包括：

对所述形变场统计特征，所述物理属性特征和所述图像相似性特征，通过不同半径的统计窗进行特征扩充，得到多个扩充特征；

根据所述形变场统计特征，所述物理属性特征和所述图像相似性特征，输入预先训练好的分类器，由所述分类器确定所述图像对的配准误差包括：

将多个所述扩充特征输入所述分类器，由所述分类器确定所述图像对的处于多个配准误差范围的概率；

将概率最高的配准误差范围，作为所述图像对的配准误差。

7.根据权利要求1所述的方法，其特征在于，根据所述形变场统计特征，所述物理属性特征和所述图像相似性特征，输入预先训练好的分类器，由所述分类器确定所述图像对的配准误差之前，所述方法还包括：

确定训练样本集中训练样本的形变场统计特征，物理属性特征，图像相似性特征，以及训练样本的目标配准误差；

根据所述训练样本的形变场统计特征，物理属性特征，图像相似性特征，以及所述训练样本的目标配准误差，对预设的分类器进行训练，直至所述分类器的分类误差达到满足预设条件。

8.一种图像配准的数据处理装置，其特征在于，包括：

生成模块，用于根据已配准的图像对，生成多个形变场，其中，所述图像对包括固定图像和浮动图像；

第一确定模块，用于根据多个形变场，确定形变场统计特征和物理属性特征；

第二确定模块，用于根据配准前后的图像，确定图像相似性特征；

分类模块，用于根据所述形变场统计特征，所述物理属性特征和所述图像相似性特征，输入预先训练好的分类器，由所述分类器确定所述图像对的配准误差。

9.一种计算机可读存储介质，其特征在于，所述存储介质用于存储程序，其中，所述程序执行权利要求1至7中任意一项所述的图像配准的数据处理方法。

10.一种电子设备，其特征在于，包括一个或多个处理器和存储器，所述存储器用于存储一个或多个程序，其中，当所述一个或多个程序被所述一个或多个处理器执行时，使得所述一个或多个处理器实现权利要求1至7中任意一项所述的图像配准的数据处理方法。

Images