CN116703994B - 医学图像配准的方法、计算设备和计算机可读存储介质 - Google Patents

医学图像配准的方法、计算设备和计算机可读存储介质 Download PDF

Info

Publication number
CN116703994B
CN116703994B CN202310951097.3A CN202310951097A CN116703994B CN 116703994 B CN116703994 B CN 116703994B CN 202310951097 A CN202310951097 A CN 202310951097A CN 116703994 B CN116703994 B CN 116703994B
Authority
CN
China
Prior art keywords
image
transformation
target image
transform
floating
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202310951097.3A
Other languages
English (en)
Other versions
CN116703994A (zh
Inventor
罗欢
罗园明
钱沛东
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Boyi Huixin Hangzhou Network Technology Co ltd
Original Assignee
Boyi Huixin Hangzhou Network Technology Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Boyi Huixin Hangzhou Network Technology Co ltd filed Critical Boyi Huixin Hangzhou Network Technology Co ltd
Priority to CN202310951097.3A priority Critical patent/CN116703994B/zh
Publication of CN116703994A publication Critical patent/CN116703994A/zh
Application granted granted Critical
Publication of CN116703994B publication Critical patent/CN116703994B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/38—Registration of image sequences
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/37—Determination of transform parameters for the alignment of images, i.e. image registration using transform domain methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10088—Magnetic resonance imaging [MRI]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10104—Positron emission tomography [PET]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30101—Blood vessel; Artery; Vein; Vascular

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Image Processing (AREA)

Abstract

本公开提供了一种医学图像配准的方法、计算设备和计算机可读存储介质。该方法包括：采集目标对象的多个时相的医学图像；从所述多个时相的医学图像中确定一个参考时相的医学图像作为目标图像并且选择另一时相的医学图像作为浮动图像；确定所述目标图像和所述浮动图像之间的位移信息，其中所述位移信息由第一变换和第二变换确定，并且使用所述第一变换从所述目标图像到所述浮动图像的路径与使用所述第二变换从所述浮动图像到所述目标图像的路径相等；以及基于所述参考时相与所述多个时相中的所有其他时相的医学图像的位移信息确定用于对所述多个时相的医学图像进行配准的位移场。

Description

医学图像配准的方法、计算设备和计算机可读存储介质
技术领域
本发明概括而言涉及医学图像处理领域，更具体地，涉及一种医学图像配准的方法、计算设备和计算机可读存储介质。
背景技术
医学图像配准是指将不同时间拍摄的医学图像或者不同模态下的医学图像，甚至不同个体拍摄得到的医学图像通过配准的方式，使得两幅图像的对应点达到空间位置和解剖结构的一致，进而实现信息互补，协助医生做出更好更准确的诊断。医学图像配准可以用于各种类型的医学图像，如CT（计算机断层扫描，Computed Tomography），MRI（核磁共振成像，Nuclear Magnetic Resonance Imaging），PET（正电子发射型计算机断层显像，Positron Emission Computed Tomography）等。医学图像配准是医学图像处理中重要的一个环节，能够为医生提供手术规划，辅助诊断等功能。
医学图像配准的核心问题可以描述为如下过程，设置一种特定的空间变换方式，在这种变换方式下，使得某种衡量两幅图像相似度的指标达到最小。这一过程可以用如下式子来表示：
其中，I代表两幅图像中的目标图像，J代表两幅图像中的浮动图像（即，待配准图像），M代表衡量两幅图像相似度的指标，T代表所设定的空间变换方式。医学图像配准的结果是求取使得上述相似度指标M最小时的变换T，即，浮动图像J相对于目标图像I的位移信息。
目前医学图像配准有如下几类常用的方法。第一类是基于外部特征的配准方法，通过人工标注的方式，在图像上选取特征信息，可以是特征点，也可以是特征曲线、特征质心、特征轴，特征面，通过对这些特征信息进行人工匹配来实现两幅图像的配准。第二类是基于光流场理论的非刚性配准算法，其中比较有代表性质的是Demons算法。该算法通过利用目标图像的梯度以及目标图像和浮动图像的灰度插值来计算浮动图像的位移量，并利用这个位移量对浮动图像进行变换，迭代这一过程，使得变换后的浮动图像和目标图像的相似度达到阈值或是达到迭代次数上限，以得到最终结果。
然而，对于上述第一类方法，配准的结果主要依赖于人工选取的特征信息，不仅费时费力，而且由于主观性影响，对同样两幅图像的多次配准可能在结果上有较大的差异。此外，该方式对于非刚性位移的配准实现起来比较困难。对于上述第二类方法，迭代收敛的速度较慢，可能导致整体运行时间较长，并且配准精度无法达到很高。
发明内容
针对上述问题中的至少一个，本发明提供了一种医学图像配准方法，其通过将目标图像和浮动图像之间的空间变换构建为两个对称变换，能够更加准确地确定位移信息，从而实现精确配准。
根据本发明的一个方面，提供了一种医学图像配准方法。该方法包括：采集目标对象的多个时相的医学图像；从所述多个时相的医学图像中确定一个参考时相的医学图像作为目标图像并且选择另一时相的医学图像作为浮动图像；确定所述目标图像和所述浮动图像之间的位移信息，其中所述位移信息由第一变换和第二变换确定，并且使用所述第一变换从所述目标图像到所述浮动图像的路径与使用所述第二变换从所述浮动图像到所述目标图像的路径相等；以及基于所述参考时相与所述多个时相中的所有其他时相的医学图像的位移信息确定用于对所述多个时相的医学图像进行配准的位移场。
在一些实现中，确定所述目标图像和所述浮动图像之间的位移信息包括：基于所述目标图像的速度场和所述浮动图像的速度场以及所述目标图像和所述浮动图像之间的互相关数确定所述目标图像和所述浮动图像之间的相似度矩阵；确定使得所述相似度矩阵最小的第一变换和第二变换；以及基于所述第一变换和所述第二变换确定所述目标图像和所述浮动图像之间的位移信息。
在一些实现中，确定所述目标图像和所述浮动图像之间的相似度矩阵之前还包括确定所述目标图像和所述浮动图像之间的互相关数，并且确定所述目标图像和所述浮动图像之间的互相关数还包括：初始化所述第一变换和所述第二变换，其中初始化的第一变换和第二变换分别为映射到自身的变换；利用初始化的第一变换对所述目标图像进行半时间间隔变换以获取第一变换目标图像，并且利用初始化的第二变换对所述浮动图像进行半时间间隔变换以获取第一变换浮动图像；基于所述第一变换目标图像和所述第一变换目标图像中任一位置周围的多个体素的平均值确定第二变换目标图像，基于所述第一变换浮动图像和所述第一变换浮动图像中任一位置周围的多个体素的平均值确定第二变换浮动图像；以及基于所述第二变换目标图像和所述第二变换浮动图像确定所述目标图像和所述浮动图像之间的互相关数。
在一些实现中，确定所述目标图像和所述浮动图像之间的相似度矩阵包括：对所述目标图像的速度场的2范数和所述浮动图像的速度场的2范数之和在半时间间隔进行积分以确定所述相似度矩阵的第一项；对所述目标图像和所述浮动图像之间的互相关数在整个图像域上进行积分以确定所述相似度矩阵的第二项；以及基于所述第一项和所述第二项确定所述相似度矩阵。
在一些实现中，确定使得所述相似度矩阵最小的第一变换和第二变换包括：对所述相似度矩阵求所述第一变换的梯度值和所述第二变换的梯度值；基于所述第一变换的梯度值对所述第一变换进行更新并且基于所述第二变换的梯度值对所述第二变换进行更新直至满足收敛条件，并且确定所述目标图像和所述浮动图像之间的位移信息包括：基于收敛条件下的第一变换和第二变换确定所述目标图像和所述浮动图像之间的位移信息。
在一些实现中，基于所述第一变换的梯度值对所述第一变换进行更新并且基于所述第二变换的梯度值对所述第二变换进行更新直至满足收敛条件包括：在对所述第一变换和所述第二变换进行更新之后，确定是否满足所述收敛条件；响应于确定满足所述收敛条件，确定所述收敛条件下的第一变换和第二变换；响应于确定不满足所述收敛条件，基于更新后的第一变换和第二变换确定更新后的相似度矩阵。
在一些实现中，所述收敛条件包括迭代次数或收敛阈值。
在一些实现中，确定所述目标图像和所述浮动图像之间的位移信息包括：基于所述第一变换及其反变换以及所述第二变换及其反变换确定所述目标图像和所述浮动图像之间的空间变换作为所述位移信息。
在一些实现中，从所述多个时相的医学图像中确定一个参考时相的医学图像作为目标图像并且选择另一时相的医学图像作为浮动图像包括：对所述参考时相的医学图像进行第一降采样以获取所述目标图像并且对所述另一时相的医学图像进行第一降采样以获取所述浮动图像。
在一些实现中，确定所述目标图像和所述浮动图像之间的位移信息还包括：基于所述目标图像和所述浮动图像确定第一降采样下的所述第一变换和所述第二变换之后，对所述参考时相的医学图像进行第二降采样以获取第二目标图像并且对所述另一时相的医学图像进行第二降采样以获取第二浮动图像，并且以第一降采样下的所述第一变换和所述第二变换作为初始化的第一变换和第二变换来获取第二降采样下的第一变换和第二变换，其中所述第二降采样的采样间隔小于所述第一降采样的采样间隔。
在一些实现中，确定所述目标图像和所述浮动图像之间的位移信息还包括：在获取第二降采样下的第一变换和第二变换之后，对所述参考时相的医学图像进行第三降采样以获取第三目标图像并且对所述另一时相的医学图像进行第三降采样以获取第三浮动图像，并且以第二降采样下的第一变换和第二变换作为初始化的第一变换和第二变换来进一步获取第三降采样下的第一变换和第二变换，其中所述第三降采样的采样间隔小于所述第二降采样的采样间隔。
在一些实现中，所述第一降采样为1/8采样，所述第二降采样为1/4采样，所述第三降采样为1/2采样。
根据本发明的另一个方面，提供了一种计算设备。该计算设备包括：至少一个处理器；以及至少一个存储器，该至少一个存储器被耦合到该至少一个处理器并且存储用于由该至少一个处理器执行的指令，该指令当由该至少一个处理器执行时，使得该计算设备执行根据上述方法的步骤。
根据本发明的再一个方面，提供了一种计算机可读存储介质，其上存储有计算机程序代码，该计算机程序代码在被运行时执行如上所述的方法。
附图说明
通过参考下列附图所给出的本发明的具体实施方式的描述，将更好地理解本发明，并且本发明的其他目的、细节、特点和优点将变得更加显而易见。
图1示出了用于实现根据本发明的实施例的医学图像配准方法的系统的示意图。
图2示出了根据本发明的一些实施例的用于医学图像配准的方法的流程图。
图3示出了根据本发明一些实施例的用于确定目标图像和所述浮动图像之间的位移信息的过程的进一步详细流程图。
图4示出了根据本发明一些实施例的确定目标图像和浮动图像之间的互相关数的过程的示意性流程图。
图5示出了根据本发明一些实施例的用于确定使得相似度矩阵最小的第一变换和第二变换的过程的流程图。
图6示出了位移信息的一个示例性示意图。
图7示出了适合实现本发明的实施例的计算设备的结构方框图。
具体实施方式
下面将参照附图更详细地描述本发明的优选实施方式。虽然附图中显示了本发明的优选实施方式，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本发明更加透彻和完整，并且能够将本发明的范围完整的传达给本领域的技术人员。
在下文的描述中，出于说明各种发明的实施例的目的阐述了某些具体细节以提供对各种发明实施例的透彻理解。但是，相关领域技术人员将认识到可在无这些具体细节中的一个或多个细节的情况来实践实施例。在其它情形下，与本申请相关联的熟知的装置、结构和技术可能并未详细地示出或描述从而避免不必要地混淆实施例的描述。
除非语境有其它需要，在整个说明书和权利要求中，词语“包括”和其变型，诸如“包含”和“具有”应被理解为开放的、包含的含义，即应解释为“包括，但不限于”。
在整个说明书中对“一个实施例”或“一些实施例”的提及表示结合实施例所描述的特定特点、结构或特征包括于至少一个实施例中。因此，在整个说明书的各个位置“在一个实施例中”或“在一些实施例”中的出现不一定全都指相同实施例。另外，特定特点、结构或特征可在一个或多个实施例中以任何方式组合。
此外，说明书和权利要求中所用的第一、第二、第三、第四等术语，仅仅出于描述清楚起见来区分各个对象，而并不限定其所描述的对象的大小或其他顺序等。
图1示出了用于实现根据本发明的实施例的医学图像配准方法的系统1的示意图。如图1中所示，系统1可以包括操作台10、扫描床20和射线发生器30，其例如可以是一个CT系统。在系统1工作时，患者可以躺在扫描床20上，医生或操作员可以通过操作台10控制扫描床20移动，以使得射线发生器30发出的射线对患者特定部位进行扫描，并将扫描产生的医学图像返回给操作台10。
在操作台10处，或者，在独立于操作台10的另一计算设备（如医生的计算设备，图中未示出）处，可以对上述产生的医学图像进行处理和分析以获取所需要的结果。在这种情况下，操作台10或另一计算设备（本文中也统称为计算设备）可以包括至少一个处理器和与该至少一个处理器耦合的至少一个存储器，该存储器中存储有可由该至少一个处理器执行的指令，该指令在被该至少一个处理器执行时执行如下所述的方法的至少一部分。计算设备的具体结构例如可以如下结合图7所述。
图2示出了根据本发明的一些实施例的用于医学图像配准的方法100的流程图。方法100例如可以由图1中所示的系统1中的操作台10或另一计算设备执行。以下以在操作台10中执行为例，结合图1至图6对方法100进行描述。
如图2中所示，方法100可以包括方框110，其中可以采集目标对象的多个时相的医学图像。这里，目标对象是指患者或者患者的特定部位，如头部、胸部、腹部等。医学图像的类型可以包括如上所述的CT图像、MRI图像、PET图像等。更具体地，在针对血管的具体应用中，该医学图像可以是CT血管成像图像（即CTA图像），尤其是多维动态CTA图像。
在方框110，例如可以采用心电门控技术，采集目标对象的多个时相的医学图像。心电门控是指在采集心脏周期性节律性运动带来的影像数据时，开启心动周期上的一段“时间窗”，使得采集数据与周期性节律性的心电活动同步，从而相对制动心脏运动的磁共振生理同步采集技术。这里，在每个心动周期采集的影像称为目标对象的一个时相的医学影像。取决于医学影像采集的不同目的，所需要采集的时相数可能不同。例如，在采集CTA图像以用于诸如动脉瘤之类的血管疾病诊断的应用中，所采集的时相数应当至少为21个，以确保采集数据足够满足对血管模型的建模。此外，医学影像采集的其他注意事项包括：时相间隔应当尽量均匀，例如相邻采样之间可以间隔1个时相或2个时相、采样时尽量使患者保持呼吸心率平稳、采集的图像伪影要少、血管流腔边界清晰可辩等。
在方框120，可以从方框110采集的多个时相的医学图像中确定一个参考时相的医学图像作为目标图像并且选择另一时相的医学图像作为浮动图像。
如前所述，图像配准是将不同时相拍摄的两个不同图像的对应点达到空间位置和解剖结构的一致。因此，可以以一个参考时相的医学图像（即目标图像）为基准，确定所有其他时相的医学图像（即浮动图像）相对于该目标图像的位移信息以构建所有医学图像的位移场信息。这里，参考时相可以是任一时相。例如，可以选择采集的起始时相作为参考时相。
在方框130，可以确定方框120选择的目标图像和浮动图像之间的位移信息。例如，对于目标图像I和浮动图像J，该位移信息可以由目标图像I和浮动图像J之间的一个空间变换T来表示。利用该空间变换T能够将目标图像I变换到新的坐标系中，即配准的变换过程，，将空间变换T作用于一个图像（如目标图像I）可以得到一个新的图像，将空间变换T的逆变换T-1作用于该新的图像则能够得到原来的图像。该变换T的参数包括时间t，空间坐标x，以及图像域上的速度场v(x,t)，其中速度场v(x,t)是平方可积的连续向量场，满足/>。
在本文中，该位移信息（即空间变换T）可以拆分为第一变换T1和第二变换T2，并且使用第一变换T1从目标图像I到浮动图像J的路径与使用第二变换T2从浮动图像J到目标图像I的路径相等。因此，本发明中的空间变换T也可以称为对称微分同胚变换。
对称微分同胚变换是将空间变换T(x,1)分为两部分进行计算，分别为第一变换T1(x,1)和第二变换T2(x,1)，两者分别作用于参考图像I和浮动图像J，并满足如下条件：
(1)
(2)
(3)
对称微分同胚变换保证了无论采用何种相似度矩阵和优化策略，从目标图像I到浮动图像J和从浮动图像J到目标图像I计算时的路径相等，从而使得浮动图像和目标图像在计算过程中的地位是对等的。
图3示出了根据本发明一些实施例的用于确定目标图像和浮动图像之间的位移信息的过程（方框130）的进一步详细流程图。
如图3中所示，在方框132，可以基于目标图像I的速度场和浮动图像J的速度场以及目标图像I和浮动图像J之间的互相关数确定目标图像I和浮动图像J之间的相似度矩阵。
如前所述，图像配准的过程是寻找到使得衡量两幅图像相似度的指标达到最小的空间变换T。在本文中，基于目标图像I本身的信息、浮动图像J本身的信息以及目标图像I和浮动图像J之间的互相关信息来构建衡量相似度的指标。
这里，目标图像I和浮动图像J之间的互相关数可以针对目标图像I和每个浮动图像J预先确定并与相应的空间位置x相关联地存储。图4示出了根据本发明一些实施例的确定目标图像I和浮动图像J之间的互相关数的过程（方框131）的示意性流程图。
如图4中所示，在方框1312，可以初始化第一变换T1和第二变换T2，其中初始化的第一变换T1和第二变换T2分别为映射到自身的变换，即,/>。例如，可以将第一变换T1和第二变换T2分别初始化为与目标图像I和浮动图像J大小相同的单位矩阵。
在方框1314，可以利用初始化的第一变换T1对目标图像I进行半时间间隔变换以获取第一变换目标图像，利用初始化的第二变换T2对浮动图像J进行半时间间隔变换以获取第一变换浮动图像/>。
在方框1316，可以基于第一变换目标图像和第一变换目标图像/>中任一位置x周围的多个体素的平均值确定第二变换目标图像/>，并且基于第一变换浮动图像/>和第一变换浮动图像/>中任一位置x周围的多个体素的平均值确定第二变换浮动图像/>。
例如，可以如下确定第二变换目标图像和第二变换浮动图像/>：
(4)
(5)
其中，表示I1中空间坐标x周围的多个体素（如5*5*5个5体素）的平均值，表示J1中空间坐标x周围的多个体素（如5*5*5个体素）的平均值。
在方框1318，可以基于第二变换目标图像和第二变换浮动图像/>确定目标图像I和浮动图像J之间的互相关数。
在一些实施例中，目标图像I和浮动图像J之间的互相关数可以表示为
（6）
可以预先确定并存储该互相关数，或者更一般性地，可以存储该互相关数的各个组成部分。例如，可以将该互相关数改记为：
（7）
这样，可以分别存储与每个空间位置x相对应的参数A、B和C。
在这种情况下，在方框132中，衡量目标图像I和浮动图像J之间的相似度的指标（即相似度矩阵）M可以表示为：
（8）
可以看出，该相似度矩阵M包括由目标图像I和浮动图像J本身的信息确定的第一项和由目标图像I和浮动图像J之间的互相关数确定的第二项。
更具体地，如公式（8）中所示，相似度矩阵M的第一项可以通过对目标图像I的速度场的2范数和浮动图像J的速度场/>的2范数之和在半时间间隔进行积分来确定，即第一项为/>。
这里，半时间间隔是指，将目标图像I和浮动图像J之间的时间间隔归一化为1，也称为单位时间间隔的情况下，该单位时间间隔的一半。
其中各个速度场(i=1或2)可以通过现有技术中或者未来开发的任何方法来获取。例如，速度场/>可以通过如下迭代过程得到：
(9)
更具体的来说，根据迭代前的，可以计算得到迭代前的/>，而后可以计算得到/>的梯度以对/>进行更新，再计算更新后的/>。
此外，如公式（8）中所示，相似度矩阵M的第二项可以通过对目标图像I和浮动图像J之间的互相关数在整个图像域/>上进行积分来确定，即第二项为。
上述第一项和第二项之和构成了目标图像I和浮动图像J之间的相似度矩阵M。
返回继续图3，在方框134，可以确定使得相似度矩阵M最小的第一变换T1和第二变换T2，即确定时的第一变换T1和第二变换T2。
在一些实施例中，可以利用梯度下降法来确定使得相似度矩阵M最小的第一变换T1和第二变换T2。图5示出了根据本发明一些实施例的用于确定使得相似度矩阵M最小的第一变换T1和第二变换T2的过程（方框134）的流程图。
如图5中所示，在方框1342，可以对相似度矩阵M分别求第一变换T1的梯度值和第二变换T2的梯度值。例如，对相似度矩阵M求第一变换T1的梯度值为：
(10)
对相似度矩阵M求第二变换T2的梯度值为：
（11)
其中，为变换Ti的雅可比矩阵（i=1或2）。
然后，在方框1344，可以分别基于第一变换T1的梯度值对第一变换T1进行更新并且基于第二变换T2的梯度值/>对第二变换T2进行更新直至满足收敛条件。
例如，更新后的第一变换T1可以表示为
(12)
更新后的第二变换T2可以表示为
(13）
其中，是收敛步长，其可以是常数（例如1）或者是可变的系数。进一步地，在收敛步长/>可变的情况下，其可以随着迭代次数而减小。
在对第一变换T1和第二变换T2进行更新之后，可以确定是否满足预定的收敛条件。这里，收敛条件可以包括迭代次数或收敛阈值。具体地，在收敛条件是迭代次数的情况下，可以确定上述基于梯度值对第一变换和第二变换进行更新的过程的执行次数是否等于该迭代次数，并且在确定执行次数等于该迭代次数的情况下，判断满足收敛条件。在收敛条件是收敛阈值的情况下，可以确定更新前后的第一变换或第二变换之差是否小于该收敛阈值，并且在小于该收敛阈值的情况下，判断满足收敛条件。在一些实例中，该收敛阈值例如是10-6的量级。
在确定满足收敛条件时，该收敛条件下的第一变换和第二变换即为最终所确定的使得相似度矩阵M最小的第一变换和第二变换。
在确定不满足收敛条件时，可以基于更新后的第一变换T'1和第二变换T'2重新执行上述方框132以确定更新后的相似度矩阵M，并且对更新后的相似度矩阵M执行上述方框134直至满足收敛条件。
继续图3，在方框136，可以基于第一变换T1和第二变换T2确定目标图像I和浮动图像J之间的位移信息。
这里的第一变换T1和第二变换T2是在上述方框1344中满足收敛条件下的第一变换T1和第二变换T2。
可以基于第一变换T1及其逆变换以及第二变换T2及其逆变换来确定上述空间变换T及其逆变换以作为该位移信息。例如，该位移信息可以表示为如下的空间变换T及其逆变换T-1：
(14)
(15)
可以看出，目标图像I和浮动图像J之间的空间变换T等于第一变换T1，并且等于对在半时间间隔的第一变换T1执行半时间间隔的第二变换T2逆变换；相应地，目标图像I和浮动图像J之间的空间变换T的逆变换等于第一变换T1的逆变换，并且等于第二变换T2，并且等于对在半时间间隔的第二变换T2执行半时间间隔的第一变换T1逆变换。
图6示出了位移信息的一个示例性示意图。如图6中所示，每个小箭头指示了目标图像与浮动图像之间的对应坐标点的位移，所有坐标点的位移的集合（即图中的所有小箭头）指示了两个图像之间的位移信息，即空间变换T。其中，图6中的R、L、S、I分别代表了真实世界空间的右侧、左侧、上侧和下侧。为了使得图像便于查看，对小箭头进行了降采样，即原来每个体素内都含有位移信息，而图6中仅为每一小块空间（包含了多个体素）显示一个位移信息。同时，为了使得箭头更易于被观察到，对位移距离进行了放大，箭头的长度按比例进行了放大，并不代表真实空间的位移距离，但箭头的方向对应真实的位移方向。
对于参考时相和每个其他时相都可以执行如上所述的过程以确定参考时相的目标图像与其他时相的浮动图像之间的位移信息。
在这种情况下，在方框140，可以基于所获取的参考时相与多个时相中的所有其他时相的医学图像的位移信息确定用于对多个时相的医学图像进行配准的位移场。
如上所述，在上述方框130中，在使用迭代方式确定目标图像和浮动图像之间的位移信息时，使用迭代方式对第一变换T1和第二变换T2进行更新直至达到收敛。在此过程中，第一变换T1和第二变换T2被初始化为单位矩阵，如果将迭代次数设置得过低或者将收敛阈值设置得过大，则很可能在满足收敛条件时并不能得到最佳的空间变换，而如果将迭代次数设置得过高或者将收敛阈值设置得过小，则可能会降低运行效率。
针对这种情况，为了提高迭代过程中的收敛效率，在本公开的一些实施例中，对上述方法100采用多尺度配准，即通过采样间隔依次降低的多次降采样来对目标图像I和浮动图像J进行降采样，并且将前一次降采样获得的第一变换T1和第二变换T2作为下一次降采样的第一变换T1和第二变换T2的初始值。
取决于不同的应用需求，采用的多尺度配准的尺度数量可以不同。例如，在采用2个尺度的情况下，可以先进行1/2降采样下的图像配准，再执行正常图像的图像配准；在采用3个尺度的情况下，可以先采用1/4降采样进行图像配准，再进行1/2降采样下的图像配准，最后执行正常图像的图像配准；在采用4个尺度的情况下，可以依次采用1/8降采样、1/4降采样、1/2降采样和正常图像的图像配准。
具体地，以4个尺度（即4层配准）为例，在上述方框120中，可以对参考时相的医学图像进行第一降采样以获取该目标图像I并且对另一时相的医学图像进行第一降采样以获取浮动图像J。例如，这里的第一降采样可以是1/8采样。通过这种方式，可以将目标图像I和浮动图像J缩小至各自的原图像的1/8。在这种情况下，可以基于降采样后的目标图像I和浮动图像J来获取第一降采样下的第一变换T1和第二变换T2。注意，此时的第一变换T1和第二变换T2并不是最终的第一变换T1和第二变换T2。
然后，与上述类似的，对参考时相的医学图像进行第二降采样以获取第二目标图像并且对该另一时相的医学图像进行第二降采样以获取第二浮动图像，并且以第一降采样下的第一变换T1和第二变换T2作为初始化的第一变换T1和第二变换T2来获取第二降采样下的第一变换T1和第二变换T2。这里，第二降采样的采样间隔小于第一降采样的采样间隔。例如，在第一降采样是1/8采样的情况下，第二降采样可以是1/4采样。
在获取第二降采样下的第一变换T1和第二变换T2之后，对参考时相的医学图像进行第三降采样以获取第三目标图像并且对该另一时相的医学图像进行第三降采样以获取第三浮动图像，并且以第二降采样下的第一变换T1和第二变换T2作为初始化的第一变换T1和第二变换T2来进一步获取第三降采样下的第一变换T1和第二变换T2。这里，第三降采样的采样间隔小于第二降采样的采样间隔。例如，在在第二降采样是1/4采样的情况下，第三降采样可以是1/2采样。
最后，可以以第三降采样下的第一变换T1和第二变换T2作为第一变换T1和第二变换T2的初始值来确定正常图像下的第一变换T1和第二变换T2，从而确定目标图像I和浮动图像J之间的位移信息。
利用本文公开的方案，通过采用对称微分同胚变换提高了医学图像配准的精度。此外，在一些实施例中，通过采用体现目标图像与浮动图像之间的互相关信息的相似度矩阵进一步提高了配准精度。进一步的，在一些实施例中，通过采用多尺度配准，能够使得图像配准加速收敛，从而进一步提高配准速度。
图7示出了适合实现本发明的实施例的计算设备700的结构方框图。计算设备700例如可以是如上所述的用于执行方法100的操作台10或另一计算设备。
如图7中所示，计算设备700可以包括一个或多个中央处理单元（CPU）710（图中仅示意性地示出了一个），其可以根据存储在只读存储器（ROM）720中的计算机程序指令或者从存储单元780加载到随机访问存储器（RAM）730中的计算机程序指令，来执行各种适当的动作和处理。在RAM 730中，还可存储计算设备700操作所需的各种程序和数据。CPU 710、ROM 720以及RAM 730通过总线740彼此相连。输入/输出（I/O）接口750也连接至总线740。
计算设备700中的多个部件连接至I/O接口750，包括：输入单元760，例如键盘、鼠标等；输出单元770，例如各种类型的显示器、扬声器等；存储单元780，例如磁盘、光盘等；以及通信单元790，例如网卡、调制解调器、无线通信收发机等。通信单元790允许计算设备700通过诸如因特网的计算机网络和/或各种电信网络与其他设备交换信息/数据。
上文所描述的方法100例如可由计算设备700（如操作台10或另一计算设备）的CPU710执行。例如，在一些实施例中，方法100可被实现为计算机软件程序，其被有形地包括于机器可读介质，例如存储单元780。在一些实施例中，计算机程序的部分或者全部可以经由ROM 720和/或通信单元790而被载入和/或安装到计算设备700上。当计算机程序被加载到RAM 730并由CPU 710执行时，可以执行上文描述的方法100的一个或多个操作。此外，通信单元790可以支持有线或无线通信功能。
本领域技术人员可以理解，图7所示的计算设备700仅是示意性的。在一些实施例中，计算设备700可以包含更多或更少的部件。
以上结合附图对根据本发明的用于进行医学图像配准的方法100以及可用作操作台10或另一计算设备的计算设备700进行了描述。然而本领域技术人员可以理解，方法100的步骤及其子步骤的执行并不局限于图中所示和以上所述的顺序，而是可以以任何其他合理的顺序来执行。此外，计算设备700也不必须包括图7中所示的所有组件，其可以仅仅包括执行本发明中所述的功能所必须的其中一些组件，并且这些组件的连接方式也不局限于图中所示的形式。
本发明可以是方法、装置、系统和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质，其上载有用于执行本发明的各个方面的计算机可读程序指令。
在一个或多个示例性设计中，可以用硬件、软件、固件或它们的任意组合来实现本发明所述的功能。例如，如果用软件来实现，则可以将所述功能作为一个或多个指令或代码存储在计算机可读介质上，或者作为计算机可读介质上的一个或多个指令或代码来传输。
本文公开的装置的各个单元可以使用分立硬件组件来实现，也可以集成地实现在一个硬件组件，如处理器上。例如，可以用通用处理器、数字信号处理器（DSP）、专用集成电路（ASIC）、现场可编程门阵列（FPGA）或其它可编程逻辑器件、分立门或者晶体管逻辑、分立硬件组件或用于执行本文所述的功能的任意组合来实现或执行结合本发明所描述的各种示例性的逻辑块、模块和电路。
本领域普通技术人员还应当理解，结合本发明的实施例描述的各种示例性的逻辑块、模块、电路和算法步骤可以实现成电子硬件、计算机软件或二者的组合。
本发明的以上描述用于使本领域的任何普通技术人员能够实现或使用本发明。对于本领域普通技术人员来说，本发明的各种修改都是显而易见的，并且本文定义的一般性原理也可以在不脱离本发明的精神和保护范围的情况下应用于其它变形。因此，本发明并不限于本文所述的实例和设计，而是与本文公开的原理和新颖性特性的最广范围相一致。

Claims (13)

1.一种医学图像配准的方法，包括：
采集目标对象的多个时相的医学图像；
从所述多个时相的医学图像中确定一个参考时相的医学图像作为目标图像并且选择另一时相的医学图像作为浮动图像；
确定所述目标图像和所述浮动图像之间的位移信息，其中所述位移信息由第一变换和第二变换确定，并且使用所述第一变换从所述目标图像到所述浮动图像的路径与使用所述第二变换从所述浮动图像到所述目标图像的路径相等；以及
基于所述参考时相与所述多个时相中的所有其他时相的医学图像的位移信息确定用于对所述多个时相的医学图像进行配准的位移场，
其中确定所述目标图像和所述浮动图像之间的位移信息包括：
基于所述目标图像的速度场和所述浮动图像的速度场以及所述目标图像和所述浮动图像之间的互相关数确定所述目标图像和所述浮动图像之间的相似度矩阵；
确定使得所述相似度矩阵最小的第一变换和第二变换；以及
基于所述第一变换和所述第二变换确定所述目标图像和所述浮动图像之间的位移信息。
2.如权利要求1所述的方法，其中确定所述目标图像和所述浮动图像之间的相似度矩阵之前还包括确定所述目标图像和所述浮动图像之间的互相关数，并且确定所述目标图像和所述浮动图像之间的互相关数还包括：
初始化所述第一变换和所述第二变换，其中初始化的第一变换和第二变换分别为映射到自身的变换；
利用初始化的第一变换对所述目标图像进行半时间间隔变换以获取第一变换目标图像，并且利用初始化的第二变换对所述浮动图像进行半时间间隔变换以获取第一变换浮动图像；
基于所述第一变换目标图像和所述第一变换目标图像中任一位置周围的多个体素的平均值确定第二变换目标图像，基于所述第一变换浮动图像和所述第一变换浮动图像中任一位置周围的多个体素的平均值确定第二变换浮动图像；以及
基于所述第二变换目标图像和所述第二变换浮动图像确定所述目标图像和所述浮动图像之间的互相关数。
3.如权利要求1所述的方法，其中确定所述目标图像和所述浮动图像之间的相似度矩阵包括：
对所述目标图像的速度场的2范数和所述浮动图像的速度场的2范数之和在半时间间隔进行积分以确定所述相似度矩阵的第一项；
对所述目标图像和所述浮动图像之间的互相关数在整个图像域上进行积分以确定所述相似度矩阵的第二项；以及
基于所述第一项和所述第二项确定所述相似度矩阵。
4.如权利要求1所述的方法，其中确定使得所述相似度矩阵最小的第一变换和第二变换包括：
对所述相似度矩阵求所述第一变换的梯度值和所述第二变换的梯度值；
基于所述第一变换的梯度值对所述第一变换进行更新并且基于所述第二变换的梯度值对所述第二变换进行更新直至满足收敛条件，并且
确定所述目标图像和所述浮动图像之间的位移信息包括：
基于收敛条件下的第一变换和第二变换确定所述目标图像和所述浮动图像之间的位移信息。
5.如权利要求4所述的方法，其中基于所述第一变换的梯度值对所述第一变换进行更新并且基于所述第二变换的梯度值对所述第二变换进行更新直至满足收敛条件包括：
在对所述第一变换和所述第二变换进行更新之后，确定是否满足所述收敛条件；
响应于确定满足所述收敛条件，确定所述收敛条件下的第一变换和第二变换；
响应于确定不满足所述收敛条件，基于更新后的第一变换和第二变换确定更新后的相似度矩阵。
6.如权利要求5所述的方法，其中所述收敛条件包括迭代次数或收敛阈值。
7.如权利要求1所述的方法，其中确定所述目标图像和所述浮动图像之间的位移信息包括：
基于所述第一变换及其反变换以及所述第二变换及其反变换确定所述目标图像和所述浮动图像之间的空间变换作为所述位移信息。
8.如权利要求1所述的方法，其中从所述多个时相的医学图像中确定一个参考时相的医学图像作为目标图像并且选择另一时相的医学图像作为浮动图像包括：
对所述参考时相的医学图像进行第一降采样以获取所述目标图像并且对所述另一时相的医学图像进行第一降采样以获取所述浮动图像。
9.如权利要求8所述的方法，其中确定所述目标图像和所述浮动图像之间的位移信息还包括：
基于所述目标图像和所述浮动图像确定第一降采样下的所述第一变换和所述第二变换之后，对所述参考时相的医学图像进行第二降采样以获取第二目标图像并且对所述另一时相的医学图像进行第二降采样以获取第二浮动图像，并且以第一降采样下的所述第一变换和所述第二变换作为初始化的第一变换和第二变换来获取第二降采样下的第一变换和第二变换，其中所述第二降采样的采样间隔小于所述第一降采样的采样间隔。
10.如权利要求9所述的方法，其中确定所述目标图像和所述浮动图像之间的位移信息还包括：
在获取第二降采样下的第一变换和第二变换之后，对所述参考时相的医学图像进行第三降采样以获取第三目标图像并且对所述另一时相的医学图像进行第三降采样以获取第三浮动图像，并且以第二降采样下的第一变换和第二变换作为初始化的第一变换和第二变换来进一步获取第三降采样下的第一变换和第二变换，其中所述第三降采样的采样间隔小于所述第二降采样的采样间隔。
11.如权利要求10所述的方法，其中所述第一降采样为1/8采样，所述第二降采样为1/4采样，所述第三降采样为1/2采样。
12. 一种计算设备，包括：
至少一个处理器；以及
至少一个存储器，所述至少一个存储器被耦合到所述至少一个处理器并且存储用于由所述至少一个处理器执行的指令，所述指令当由所述至少一个处理器执行时，使得所述计算设备执行根据权利要求1至11中任一项所述的方法的步骤。
13.一种计算机可读存储介质，其上存储有计算机程序代码，所述计算机程序代码在被运行时执行如权利要求1至11中任一项所述的方法。
CN202310951097.3A 2023-07-31 2023-07-31 医学图像配准的方法、计算设备和计算机可读存储介质 Active CN116703994B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310951097.3A CN116703994B (zh) 2023-07-31 2023-07-31 医学图像配准的方法、计算设备和计算机可读存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310951097.3A CN116703994B (zh) 2023-07-31 2023-07-31 医学图像配准的方法、计算设备和计算机可读存储介质

Publications (2)

Publication Number Publication Date
CN116703994A CN116703994A (zh) 2023-09-05
CN116703994B true CN116703994B (zh) 2023-10-24

Family

ID=87843587

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310951097.3A Active CN116703994B (zh) 2023-07-31 2023-07-31 医学图像配准的方法、计算设备和计算机可读存储介质

Country Status (1)

Country Link
CN (1) CN116703994B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117476241B (zh) * 2023-12-28 2024-04-19 柏意慧心(杭州)网络科技有限公司 用于确定血管的血流量的方法、计算设备和介质

Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6611615B1 (en) * 1999-06-25 2003-08-26 University Of Iowa Research Foundation Method and apparatus for generating consistent image registration
CN103400376A (zh) * 2013-07-19 2013-11-20 南方医科大学 一种乳腺动态增强磁共振图像序列的配准方法
CN104091337A (zh) * 2014-07-11 2014-10-08 北京工业大学 一种基于PCA及微分同胚Demons的变形医学图像配准方法
CN106097347A (zh) * 2016-06-14 2016-11-09 福州大学 一种多模态医学图像配准与可视化方法
CN109461140A (zh) * 2018-09-29 2019-03-12 沈阳东软医疗系统有限公司 图像处理方法及装置、设备和存储介质
CN110473233A (zh) * 2019-07-26 2019-11-19 上海联影智能医疗科技有限公司 配准方法、计算机设备和存储介质
CN110473234A (zh) * 2019-09-04 2019-11-19 中国科学院近代物理研究所 微分同胚Demons图像配准方法、系统及存储介质
CN110533641A (zh) * 2019-08-20 2019-12-03 东软医疗系统股份有限公司 一种多模态医学图像配准方法和装置
CN112116642A (zh) * 2020-09-27 2020-12-22 上海联影医疗科技股份有限公司 医学图像的配准方法、电子装置和存储介质
CN113628259A (zh) * 2021-06-09 2021-11-09 维沃移动通信（杭州）有限公司 图像的配准处理方法及装置
CN115170622A (zh) * 2022-05-11 2022-10-11 复旦大学 基于transformer的医学图像配准方法及系统

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR102273020B1 (ko) * 2014-06-18 2021-07-05 삼성전자주식회사 의료 영상 정합 방법 및 그 장치
US20210049757A1 (en) * 2019-08-14 2021-02-18 Nvidia Corporation Neural network for image registration and image segmentation trained using a registration simulator
CN113808178A (zh) * 2020-06-11 2021-12-17 通用电气精准医疗有限责任公司 图像配准方法及其模型训练方法

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6611615B1 (en) * 1999-06-25 2003-08-26 University Of Iowa Research Foundation Method and apparatus for generating consistent image registration
CN103400376A (zh) * 2013-07-19 2013-11-20 南方医科大学 一种乳腺动态增强磁共振图像序列的配准方法
CN104091337A (zh) * 2014-07-11 2014-10-08 北京工业大学 一种基于PCA及微分同胚Demons的变形医学图像配准方法
CN106097347A (zh) * 2016-06-14 2016-11-09 福州大学 一种多模态医学图像配准与可视化方法
CN109461140A (zh) * 2018-09-29 2019-03-12 沈阳东软医疗系统有限公司 图像处理方法及装置、设备和存储介质
CN110473233A (zh) * 2019-07-26 2019-11-19 上海联影智能医疗科技有限公司 配准方法、计算机设备和存储介质
CN110533641A (zh) * 2019-08-20 2019-12-03 东软医疗系统股份有限公司 一种多模态医学图像配准方法和装置
CN110473234A (zh) * 2019-09-04 2019-11-19 中国科学院近代物理研究所 微分同胚Demons图像配准方法、系统及存储介质
CN112116642A (zh) * 2020-09-27 2020-12-22 上海联影医疗科技股份有限公司 医学图像的配准方法、电子装置和存储介质
CN113628259A (zh) * 2021-06-09 2021-11-09 维沃移动通信（杭州）有限公司 图像的配准处理方法及装置
CN115170622A (zh) * 2022-05-11 2022-10-11 复旦大学 基于transformer的医学图像配准方法及系统

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Evaluation of a Learning-based Deformable Registration Method on Abdominal CT Images;R. Bhattacharjee;F. Heitz;V. Noblet;S. Sharma;N. Sharma;IRBM;第42卷(第2期);第94-105页 *
医学图像配准技术;申艳平;;中国医学物理学杂志(01);第49-53页 *
李碧草 ; 张俊峰 ; 杨冠羽 ; 舒华忠.基于结构图像表示和微分同胚Demons算法的多模态医学图像配准.东南大学学报(自然科学版).2015,(第5期),第39-43页. *

Also Published As

Publication number Publication date
CN116703994A (zh) 2023-09-05

Similar Documents

Publication Publication Date Title
CN110475505B (zh) 利用全卷积网络的自动分割
CN109074639B (zh) 医学成像系统中的图像配准系统和方法
Sermesant et al. Deformable biomechanical models: Application to 4D cardiac image analysis
Frangi et al. Model-based quantitation of 3-D magnetic resonance angiographic images
Lynch et al. Segmentation of the left ventricle of the heart in 3-D+ t MRI data using an optimized nonrigid temporal model
Ferrant et al. Serial registration of intraoperative MR images of the brain
CN107886508B (zh) 差分减影方法和医学图像处理方法及系统
Park et al. GGO nodule volume-preserving nonrigid lung registration using GLCM texture analysis
Rohkohl et al. Interventional 4D motion estimation and reconstruction of cardiac vasculature without motion periodicity assumption
US8861891B2 (en) Hierarchical atlas-based segmentation
CN112508965A (zh) 医学影像中正常器官的轮廓线自动勾画系统
CN116705330B (zh) 确定血管壁的弹性特征的方法、计算设备和介质
CN110546685B (zh) 图像分割和分割预测
CN116703994B (zh) 医学图像配准的方法、计算设备和计算机可读存储介质
US7650025B2 (en) System and method for body extraction in medical image volumes
JP2008511395A (ja) 一連の画像における動き修正のための方法およびシステム
Heyde et al. Anatomical image registration using volume conservation to assess cardiac deformation from 3D ultrasound recordings
Kim et al. Automatic segmentation of the left ventricle in echocardiographic images using convolutional neural networks
CN108898578B (zh) 一种医疗图像的处理方法、装置及计算机存储介质
Lamash et al. Strain analysis from 4-D cardiac CT image data
Chenoune et al. Segmentation of cardiac cine-MR images and myocardial deformation assessment using level set methods
CN115830016A (zh) 医学图像配准模型训练方法及设备
CN112381824B (zh) 一种对图像的几何特征进行提取的方法和相关产品
Li et al. RSU-Net: U-net based on residual and self-attention mechanism in the segmentation of cardiac magnetic resonance images
CN116664635B (zh) 构建目标对象的多维动态模型的方法、计算设备和介质

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant