CN113870324A - 多模态图像的配准方法及其配准装置和计算机可读存储介质 - Google Patents

多模态图像的配准方法及其配准装置和计算机可读存储介质 Download PDF

Info

Publication number
CN113870324A
CN113870324A CN202010609691.0A CN202010609691A CN113870324A CN 113870324 A CN113870324 A CN 113870324A CN 202010609691 A CN202010609691 A CN 202010609691A CN 113870324 A CN113870324 A CN 113870324A
Authority
CN
China
Prior art keywords
image
registered
registration
reference image
local
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.)
Pending
Application number
CN202010609691.0A
Other languages
English (en)
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.)
Shanghai Weiwei Medical Technology Co ltd
Original Assignee
Shanghai Weiwei Medical 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 Shanghai Weiwei Medical Technology Co ltd filed Critical Shanghai Weiwei Medical Technology Co ltd
Priority to CN202010609691.0A priority Critical patent/CN113870324A/zh
Publication of CN113870324A publication Critical patent/CN113870324A/zh
Pending legal-status Critical Current

Links

Images

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/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/40—Image enhancement or restoration using histogram techniques
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • 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/10132—Ultrasound image

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Quality & Reliability (AREA)
  • Computing Systems (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Measuring And Recording Apparatus For Diagnosis (AREA)

Abstract

本发明提供了一种多模态图像的配准方法及其配准装置和计算机可读存储介质，所述配准方法包括步骤：获取参考图像和待配准图像；对所述参考图像和所述待配准图像进行预处理；基于互信息作为配准测度，得到第一步的旋转平移矩阵，采用第一步的旋转平移矩阵进行初步配准；分别提取所述参考图像和所述待配准图像的目标局部对象的关键点特征向量；以及将提取到关键点特征向量进行配准，得到第二步的旋转平移矩阵，采用第二步的旋转平移矩阵进行目标局部对象配准。本发明通过基于互信息作为配准测度的初步配准和关键点特征向量配准，较好解决了局部组织形变的问题，有助于提高配准精度，辅助医生提高诊断的准确性。

Description

多模态图像的配准方法及其配准装置和计算机可读存储介质
技术领域
本发明涉及医学图像配准技术领域，具体涉及一种多模态图像的配准方法及其配准装置和计算机可读存储介质。
背景技术
图像配准在医学图像处理上具有重要意义，不同模态的图像能够提供相同组织不同层面的信息，便于医生提高疾病诊断率。医学图像配准的常见模态包括但不限于：磁共振成像(MRI)、计算机断层摄影(CT)以及超声成像(US)。
例如在肺部病灶检测中，CT图像是目前肺病最为理想的功能分子显像方法，可用于肺癌原发病灶的评估、治疗计划的设计、愈后预测、疗效评价和复发转移预测；而支气管内US图像可用于指导肺外周病灶经支气管肺活检，或指导肺门纵隔淋巴结肿块进行气道内超声引导下的针吸活检。将CT图像和支气管内US图像的两种不同模态的医学影像结合，能够更精确的对患者进行诊断。又例如在深脑部刺激手术中(简称DBS)，将电极植入脑深部特定区域，通过电脉冲发生器对脑内相关神经核团进行电刺激，从而控制异常脑活动，可显著改善遗传性的肌张力障碍(帕金森氏病)。但植入电极有误时会导致很多副作用，因此，术前医生常对患者进行MRI扫描，以获取具有较高淸晰度的软组织图像，且在术后还对颅骨和电极等进行密度敏感的CT扫描，最后将MRI图像和CT图像的两种不同模态的医学影像结合，从而获得脑软组织、颅骨和电极等结构的完整信息。
因此，需要利用图像配准技术来实现不同模态的医学影像的结合。目前已经有很多的图像配准技术，但多模态图像的配准精度一直没有得到较好的解决。目前多模态图像配准方法主要有半自动配准和自动配准；半自动配准主要基于人为标记点进行配准，实用性不高；自动配准完全基于图像处理进行配准，对图像处理技术要求较高。现有主要以自动配准为主。在自动配准过程中，患者生理状态、手术姿势、呼吸、心跳等变化都会引起器官或组织发生局部形变，导致配准精度不足。例如，现有技术中公开的一种图像配准方法，通过建立一具有多条异面直线的配准模型，根据多条异面直线中每两条异面直线之间最短距离，确定第一点分别在参考图像和待配准图像间的映射数据，利用两组映射数据计算旋转平移矩阵，最终实现多模态图像之间的配准，但未考虑器官或组织发生局部形变的问题，导致配准出现误差。再如，现有技术中还公开了一种图像配准方法，采用B样条曲面作为变形模型实现多模态图像的配准，并不断增大采样比例进行循环迭代来提高精度，但该方法较为耗时，仅考虑组织整体形变，未考虑组织的局部形变问题。再如，现有技术中还公开了一种多模态可变形配准，通过计算不同模态图像的变形场，在变形场向量的基础上进行配准，但针对不同的组织，需要提前准备组织不同形变状态的高清晰度模板图像，实用性不强。
因此，现有的多模态图像配准大多存在未考虑组织局部形变导致的配准偏差问题，配准精度较低。
发明内容
为了解决上述技术问题，本发明的目的是提供一种多模态图像的配准方法及其配准装置和计算机可读存储介质，以解决医学多模态图像中组织局部形变的问题，提高配准精度。
为了实现上述目的，根据本发明的第一方面提供了一种多模态图像的配准方法，所述配准方法包括如下步骤：
获取参考图像和待配准图像；
对所述参考图像和所述待配准图像进行预处理；
基于互信息作为配准测度，得到第一步的旋转平移矩阵，采用第一步的旋转平移矩阵对预处理后的所述参考图像和所述待配准图像进行初步配准；
分别提取所述参考图像和所述待配准图像的目标局部对象的关键点特征向量；以及
将提取到的所述参考图像和所述待配准图像的目标局部对象的关键点特征向量进行配准，得到第二步的旋转平移矩阵，采用第二步的旋转平移矩阵对初步配准后的所述参考图像和所述待配准图像进行目标局部对象配准。
可选的，对所述参考图像和/或所述待配准图像进行预处理的步骤包括：
采用高斯平滑滤波器对所述参考图像和/或所述待配准图像进行平滑处理。
可选的，基于互信息作为配准测度，得到第一步的旋转平移矩阵，采用第一步的旋转平移矩阵对预处理后的所述参考图像和所述待配准图像进行初步配准的步骤包括：
分别将所述参考图像和所述待配准图像所处的像素坐标系变换为物理坐标系，所述物理坐标系为当前物理坐标系；
基于图像灰度，分别计算所述参考图像和所述待配准图像的质心，并计算所述参考图像和所述待配准图像的质心之间的平移矩阵，且将所述平移矩阵作为初始矩阵；
基于初始矩阵，通过插值运算将所述待配准图像从所述当前物理坐标系变换到新物理坐标系中；
统计在所述当前物理坐标系下的所述参考图像和在所述新物理坐标系下的所述待配准图像的重叠区域的对应像素对的不同灰度组合出现的频度，得到联合直方图，并得到所述参考图像和所述待配准图像的联合熵；以及
根据联合熵计算互信息度量，进行随机梯度迭代，得到第一步的旋转平移矩阵，采用第一步的旋转平移矩阵对所述参考图像和所述待配准图像进行初步配准。
可选的，在计算互信息度量前，对联合熵进行归一化处理。
可选的，在将所述参考图像和所述待配准图像所处的像素坐标系变换为物理坐标系的步骤中，利用标签信息将所述参考图像和所述待配准图像所处的像素坐标系变换为物理坐标系。
可选的，所述标签信息包括：单个像素与毫米转换的间距值、像素坐标系原点与物理坐标系原点间的偏移值、和像素坐标系的方向值。
可选的，分别提取所述参考图像和所述待配准图像的目标局部对象的关键点特征向量的步骤包括：
对所述参考图像和所述待配准图像的目标局部对象选取感兴趣区域，分别得到局部参考图像和局部待配准图像；
寻找所述局部参考图像和所述局部待配准图像中的离散极值点；
将每个离散极值点置于多象限坐标系中，分别以每个离散极值点为中心，在所述局部参考图像和所述局部待配准图像中划分出对应的图像块，并将图像块分割为多个图像单元；
计算每个图像单元内各像素点在物理坐标系各坐标轴的梯度，并将每个图像单元内所有像素点的方向和梯度幅值分别累加到多象限坐标系的各象限上，形成梯度直方图；以及
根据梯度直方图，得到每个离散极值点所对应的特征向量作为关键点特征向量。
可选的，寻找所述局部参考图像和所述局部待配准图像中的离散极值点的步骤包括：
选取具有不同标准偏差的高斯分布作为不同尺度，分别对所述局部参考图像和所述局部待配准图像进行高斯平滑，并计算出高斯差分金字塔；以及
在高斯差分金字塔中遍历所有像素点的特定邻域，寻找特定邻域内灰度值最大的点作为离散极值点。
可选的，所述配准方法用于二维图像的配准，所述特定邻域为像素点的上下左右前后六邻域或所有二十六邻域，所述多象限坐标系为由二维平面的360度均分为8块～32块象限的平面坐标系，所述图像块为以离散极值点为中心，以对应高斯分布的标准偏差的3-6倍的值为半径形成的圆，所述图像单元为所述图像块的内接最大正方形均分成的M块正方形中的除去离散极值点所处的中间一块正方形的其余M-1块正方形，其中，M为大于一的正整数的平方数。
可选的，所述配准方法用于三维图像的配准，所述特定邻域为像素点的上下左右前后六邻域或所有二十六邻域，所述多象限坐标系为32球面的球坐标系，所述图像块为以离散极值点为中心，以对应高斯分布的标准偏差的3-6倍的值为半径形成的球，所述图像单元为所述图像块的内接最大正方体均分成的M块正方体中的除去离散极值点所处的中间一块正方体的其余M-1块正方体，其中，M为大于一的正整数的立方数。
可选的，将提取到的所述参考图像和所述待配准图像的目标局部对象的关键点特征向量进行配准，得到第二步的旋转平移矩阵，采用第二步的旋转平移矩阵对初步配准后的所述参考图像和所述待配准图像进行目标局部对象配准的步骤包括：
从局部参考图像和局部待配准图像的所有关键点特征向量中随机采样多对关键点特征向量；
根据采样的多对关键点特征向量计算局部参考图像和局部待配准图像之间的仿射变换矩阵；
将仿射变换矩阵应用到局部待配准图像的所有关键点特征向量，计算误差值；
当误差值符合预定条件时，将当前仿射变换矩阵作为有效的仿射变换矩阵并作为第二步的旋转平移矩阵；
当误差值不符合预定条件时，继续进行随机采样，并重新计算仿射变换矩阵，再次计算误差值，直至误差值符合预定条件后将重新计算的仿射变换矩阵作为第二步的旋转平移矩阵。
根据本发明的第二方面提供了一种多模态图像的配准装置，所述配准装置包括：
图像获取单元，用于获取参考图像和待配准图像；
图像预处理单元，用于对所述参考图像和所述待配准图像进行预处理；
初步配准单元，基于互信息作为配准测度，得到第一步的旋转平移矩阵，采用第一步的旋转平移矩阵对预处理后的所述参考图像和所述待配准图像进行初步配准；以及
目标局部对象配准单元，用于分别提取所述参考图像和所述待配准图像的目标局部对象的关键点特征向量，以及将提取到的所述参考图像和所述待配准图像的目标局部对象的关键点特征向量进行配准，得到第二步的旋转平移矩阵，采用第二步的旋转平移矩阵对初步配准后的所述参考图像和所述待配准图像进行目标局部对象配准。
根据本发明的第三方面提供了一种计算机可读存储介质，其上存储有计算机程序，当所述计算机程序被处理器执行时执行如本发明任一技术方案所述的多模态图像的配准方法的步骤。
本发明的有益效果为：
在局部组织形变后，原本坐标上的组织发生偏移，变成了其他组织，现有图像配准技术未考虑局部组织形变问题，当医生结合图像进行诊断时，会带来干扰甚至诊断错误，与现有技术相比，本发明的多模态图像配准方法及其配准装置和计算机可读存储介质通过基于互信息作为配准测度的初步配准和关键点特征向量配准，在多模态医学图像配准中，可选取感兴趣的局部组织作为目标局部对象，提取目标局部对象的关键点特征向量进行配准，较好解决了局部组织形变的问题，有助于提高配准精度，辅助医生提高诊断的准确性。
附图说明
图1为本发明一实施例的多模态图像的配准方法的流程图；
图2为步骤S3的流程图；
图3为步骤S4的流程图；
图4为步骤S42的流程图；
图5为步骤S5的流程图；
图6为术前人脑轴状面的MRI图像；
图7为术前人脑冠状面的MRI图像；
图8为术前人脑矢状面的MRI图像；
图9为术后人脑轴状面的CT图像；
图10为术后人脑冠状面的CT图像；
图11为术后人脑矢状面的CT图像；
图12为配准前同一物理坐标系下的人脑轴状面MRI图像和CT图像；
图13为配准前同一物理坐标系下的人脑冠状面MRI图像和CT图像；
图14为配准前同一物理坐标系下的人脑矢状面MRI图像和CT图像；
图15为配准后同一物理坐标系下的人脑轴状面MRI图像和CT图像；
图16为配准后同一物理坐标系下的人脑冠状面MRI图像和CT图像；
图17为配准后同一物理坐标系下的人脑矢状面MRI图像和CT图像。
具体实施方式
为使本发明的内容更加清楚易懂，以下结合说明书附图来说明本发明一实施例的多模态图像的配准方法及其配准装置和计算机可读存储介质。但可以理解，本发明并不局限于下面所描述的具体实施例，本领域的技术人员所熟知的一般替换也涵盖在本发明的保护范围内。需说明的是，附图均采用非常简化的形式且均使用非精准的比例，仅用以方便、明晰地辅助说明本发明实施例的目的。
如图1所示，本发明一实施例提供了一种多模态图像的配准方法，所述配准方法包括如下步骤：
步骤S1：获取参考图像和待配准图像；
在步骤S1中，本实施例以人脑为例进行说明，术前拍摄MRI影像，作为参考图像(如图6-图8所示)，术后拍摄CT影像(如图9-图11所示)，作为待配准图像，本实施例提供的配准方法即为将所述CT图像配准至MRI图像，当然本发明并不以此为限。所述图像可以是磁共振成像(MRI)、计算机断层摄影(CT)或超声成像(US)等获得的医学影像，例如MRI人脑体数据图像，CTA人脑体数据图像等。其中，CTA图像是计算机断层摄影(CT)获得的一种特殊的CT图像，可获得组织的三维数据结构。
步骤S2：对所述参考图像和所述待配准图像进行预处理；对所述参考图像和所述待配准图像进行预处理指的是消除图像中无关的信息，增强有关信息的可检测性，从而改善后续图像处理步骤的可靠性的图像预处理操作，包括去除图像噪声和增强图像等。
步骤S3：基于互信息(Mutual Information)作为配准测度，得到第一步的旋转平移矩阵，采用第一步的旋转平移矩阵对预处理后的所述参考图像和所述待配准图像进行初步配准；
步骤S4：分别提取所述参考图像和所述待配准图像的目标局部对象的关键点特征向量；
在步骤S4中，以人脑为例，人脑包括颅骨和颅骨内的神经核团两部分。在步骤S3基于互信息进行初步配准的步骤中，容易出现颅骨大部分区域对齐，而神经核团未对齐的情况。在DBS手术中，需要在术前通过MRI图像获取脑软组织的信息，以及在术后通过CT图像获取颅骨和电极的信息，而电极通常被植入脑深部神经核团附近，因此，本实施例中选取神经核团所处的区域为所述参考图像和所述待配准图像的目标局部对象，并对目标局部对象进行配准，以提高配准精度，辅助医生在DBS手术中进行诊断。
步骤S5：将提取到的所述参考图像和所述待配准图像的目标局部对象的关键点特征向量进行配准，得到第二步的旋转平移矩阵，采用第二步的旋转平移矩阵对初步配准后的所述参考图像和所述待配准图像进行目标局部对象配准。
在步骤S2中，对所述参考图像和/或所述待配准图像进行预处理的步骤包括采用高斯平滑滤波器对所述参考图像和/或所述待配准图像进行处理。高斯平滑滤波器可滤除图像中的噪声信息，并同时在一定程度上可模糊图像，以在步骤S3中基于互信息作为配准测度时(互信息为非光滑单峰函数，存在许多局部最优值)，避免陷入局部最优。
如图2所示，在步骤S3中，基于互信息作为配准测度，，得到第一步的旋转平移矩阵，采用第一步的旋转平移矩阵对预处理后的所述参考图像和所述待配准图像进行初步配准的步骤包括：
步骤S31：分别将所述参考图像和所述待配准图像所处的像素坐标系变换为物理坐标系。所述物理坐标系指的是图像变换后的最终的输出空间，通常以实际长度单位衡量各像素点的位置。
例如，利用所述参考图像和所述待配准图像的标签信息(Tag信息)，所述标签信息包括各图像单个像素与毫米转换的间距值(spacing值)、像素坐标系原点与物理坐标系之间的原点偏移值(origin值)和像素坐标系的方向值(direction值)等，将所述参考图像和所述待配准图像所处的像素坐标系变换为物理坐标系，以实现所述参考图像和所述待配准图像的尺寸维度信息的统一。
如图12-图14所示，分别为配准前在同一物理坐标系下的人脑轴状面、冠状面和矢状面的MRI图像和CT图像，MRI图像和CT图像未配准前不能完全重叠，人脑组织在旋转角度和平移距离上都具有偏移，需要进行配准，还不能作为医生诊断的依据。
步骤S32：基于图像灰度，分别计算所述参考图像和所述待配准图像的质心，并计算所述参考图像和所述待配准图像的质心之间的平移矩阵，且将所述平移矩阵作为初始矩阵。这里，通过计算所述初始矩阵，可减少在后续步骤计算中不必要的迭代。
步骤S33：基于初始矩阵，通过插值运算将所述待配准图像从所处的物理坐标系变换到新物理坐标系中。
所述插值运算用于解决待配准图像在变换到新坐标系后，为新物理坐标系处的像素值赋灰度值。本实施例中，插值运算可采用最近邻法、双线性法、部分体积分布插值法(Partial Volume Distribution)等方法，相比其它一些高阶的插值运算的方法，计算量较小。优选采用部分体积分布插值法，该方法利用线性插值的权重分配原则进行贡献分配，相对最近邻法和双线性法光滑性较好，有利于在后续步骤中优化离散极值点的搜索。
步骤34：统计在当前物理坐标系下的所述参考图像和在新物理坐标系下的所述待配准图像的重叠区域的对应像素对的不同灰度组合出现的频度，得到联合直方图，并得到所述参考图像和所述待配准图像的联合熵。
具体的，统计所述参考图像和所述待配准图像的重叠区域对应像素对的不同灰度组合(用(i，j)表示，其中i和j分别表示所述参考图像和所述待配准图像的灰度)出现次数，得到两者之间的联合直方图，并得到HAB(i，j)值。其中，A和B分别表示所述参考图像和所述待配准图像，HAB(i，j)为所述参考图像和所述待配准图像的联合熵。
步骤35：根据联合熵计算互信息度量，进行随机梯度迭代，得到第一步的旋转平移矩阵，采用第一步的旋转平移矩阵对所述参考图像和所述待配准图像进行初步配准。所述联合熵优选采用归一化处理，再以归一化的联合熵计算互信息度量，即采用归一化的互信息作为配准测度。采用互信息作为配准测度进行初步配准的本质是由于当所述参考图像和所述待配准图像的空间位置达到一致时，其互信息度量值应该达到最大值。但若单以基于最大互信息的方法进行配准，精确性和鲁棒性不足，原因在于，基于最大互信息的方法将图像中的所有像素点平等对待，容易产生误配，而归一化互信息方法可以较好解决该问题。
根据联合熵或归一化的联合熵计算互信息度量的计算过程具体为：将得到的联合熵HAB(i，j)值或归一化后的值依次代入下述式1至式4，计算得到互信息度量I(A，B)，并进行随机梯度迭代，最终得到第一步的旋转平移矩阵。
Figure BDA0002560551700000091
Figure BDA0002560551700000092
PA(i)＝∑jHAB(i，j) (式3)
PB(j)＝∑iHAB(i，j) (式4)
其中，i和j分别表示所述参考图像和所述待配准图像的灰度，A和B分别表示所述参考图像和所述待配准图像，I(A，B)为互信息度量，PAB(i，j)为联合概率分布，PA(i)为所述参考图像的边缘概率分布，PB(j)为所述待配准图像的边缘概率分布。
如图3所示，在步骤S4中，分别提取所述参考图像和所述待配准图像的目标局部对象的关键点特征向量的步骤包括：
步骤S41：对所述参考图像和所述待配准图像的目标局部对象选取感兴趣区域，分别得到局部参考图像和局部待配准图像。
本实施例中，选取神经核团所在的区域作为感兴趣区域，感兴趣区域可以基于先验知识选取，例如可基于神经核团与颅骨的相对位置进行选取，较优的，可将感兴趣区域制成配准文件形式直接读取，以减少人工干预成本。
步骤S42：寻找所述局部参考图像和所述局部待配准图像中的离散极值点。
步骤S43：将每个离散极值点置于多象限坐标系中，分别以每个离散极值点为中心，在所述局部参考图像和所述局部待配准图像中划分出对应的图像块，并将图像块分割为多个图像单元。
步骤S44：计算每个图像单元内各像素点在物理坐标系各坐标轴的梯度，并将每个图像单元内所有像素点的方向和梯度幅值分别累加到多象限坐标系的各象限上，形成梯度直方图。
步骤S45：根据梯度直方图，得到每个离散极值点所对应的特征向量作为关键点特征向量。
本实施例的分别提取所述参考图像和所述待配准图像的目标局部对象的关键点特征向量的方法结合了实际图像，构建了不同的梯度直方图以更好地描述关键点信息。
如图4所示，在步骤S42中，寻找所述局部参考图像和所述局部待配准图像中的离散极值点的步骤包括：
步骤S421：选取具有不同标准偏差的高斯分布作为不同尺度，分别对所述局部参考图像和所述局部待配准图像进行高斯平滑，并计算出高斯差分金字塔。例如，可通过如下步骤获取所述高斯差分金字塔：选取5种σ(σ为高斯分布的标准偏差)的高斯分布作为不同尺度，分别对所述局部参考图像和所述局部待配准图像进行高斯平滑，各得到5组图像，分别对其两两计算差值，分别得到4组图像形成高斯差分金字塔。当然，可选择其它尺度设置的高斯平滑方式，本发明对此不作限定。
步骤S422：在高斯差分金字塔中遍历所有像素点的特定邻域，寻找特定邻域内灰度值最大的点作为离散极值点。所述特定邻域为像素点的上下左右前后六邻域或所有二十六邻域，此处的所述二十六邻域是指将一块正方体均分为3*3*3块小正方体，一像素点置于中间一块小正方体处，则该像素点的周围共有26块小正方体，即二十六邻域。
在所述配准方法用于二维图形的配准时，所述特定邻域为像素点的上下左右前后六邻域或所有二十六邻域。所述多象限坐标系为由二维平面的360度均分为8块～32块象限的平面坐标系，优选为32块象限。所述图像块为以离散极值点为中心，以对应高斯分布的标准偏差的3-6倍的值为半径形成的圆，优选以对应高斯分布的标准偏差的4.5倍的值为半径(即r＝4.5*σ，σ为离散极值点对应的高斯分布的标准偏差)。所述图像单元为所述图像块的内接最大正方形均分成的M块正方形中的除去离散极值点所处的中间一块正方形的其余M-1块正方形，其中，M为大于一的正整数的平方数，例如为4块、8块、9块或16块等。
以下以所述多象限坐标系选用32象限的平面坐标系，所述图像块选用对应高斯分布的标准偏差的4.5倍的值为半径，所述图像单元选用8块正方形(即M值为9)的分割方式举例，说明所述配准方法用于二维图形的配准时的具体步骤：
所述步骤S43即：将每个离散极值点置于由二维平面的360度均分为32块的32象限的平面坐标系中，且以离散极值点为特征点中心，以对应高斯分布的标准偏差的4.5倍的值为半径(即r＝4.5*σ，σ为离散极值点对应的高斯分布的标准偏差)形成一个圆，将圆的内接最大正方形均分为9块正方形，除去离散极值点所处的中间一块正方形，共得到8块正方形。
所述步骤S44即：分别计算8块正方形内各像素点在物理坐标系的x轴和y轴方向的梯度，得到像素点的梯度，将各个像素点的梯度置于32象限的平面坐标系中，计算像素点的方向和梯度幅值，将8块正方形内的所有像素点的方向与梯度幅值分别累加到平面坐标系的32象限上，形成二维梯度直方图。
所述步骤S45即：根据二维梯度直方图，得到每个离散极值点所对应的8*32个特征向量，若存在N个离散极值点，则共得到N*8*32个特征向量作为关键点特征向量；N为大于1的正整数。
采用上述步骤，通过一种获得离散极值点的方法可较好地得到二维图像的目标局部对象的关键点特征向量，精确性较高，可普适于各种二维医学影像的图像对准，以辅助医生提高诊断的准确性。
参照上述二维图像提取关键点特征向量的方法，本实施例还提供了一种可用于三维图像的关键点特征向量的方法，在所述配准方法用于三维图像的配准时，所述特定邻域为像素点的上下左右前后六邻域或所有二十六邻域，所述多象限坐标系为32球面的球坐标系。所述图像块为以离散极值点为中心，以对应高斯分布的标准偏差的3-6倍的值为半径形成的球，优选以对应高斯分布的标准偏差的4.5倍的值为半径(即r＝4.5*σ，σ为离散极值点对应的高斯分布的标准偏差)。所述图像单元为所述图像块的内接最大正方体均分成的M块正方体中的除去离散极值点所处的中间一块正方体的其余M-1块正方体，其中，M为大于一的正整数的立方数，例如为27块。
以下以所述多象限坐标系选用32象限的球面坐标系，所述图像块选用对应高斯分布的标准偏差的4.5倍的值为半径，所述图像单元选用26块正方体(即M值为27)的分割方式举例，说明所述配准方法用于三维图形的配准时的具体步骤：
所述步骤S43即：将每个离散极值点置于32球面的球坐标系(即32象限的球坐标系)中，且以离散极值点为特征点中心，以对应高斯分布的标准偏差的4.5倍的值为半径形成一个球，将球的内接最大正方体均分为27块正方体，除去离散极值点所处的中间一块正方体，共得到26块正方体。
所述步骤S44即：分别计算26块正方体内各像素点在物理坐标系的x轴、y轴和z轴方向的梯度得到像素点的梯度，将各个像素点的梯度置于32球面的球坐标系中，计算像素点的方向和梯度幅值，将26块正方体内的所有像素点的方向与梯度幅值分别累加到球坐标系的32球面上，形成三维梯度直方图。
其中，像素点的梯度幅值即为主梯度的模，像素点在球坐标系中的主方向可按下式(式5)计算：
Figure BDA0002560551700000131
其中，θ为像素点的主方向，
Figure BDA0002560551700000132
为像素点的主梯度，
Figure BDA0002560551700000133
为像素点的主梯度落入的球面的单位法向量。
所述步骤S45即：根据三维梯度直方图，得到每个离散极值点所对应的26*32个特征向量，若存在N个离散极值点，则共得到N*26*32个特征向量作为关键点特征向量；N为大于1的正整数。
采用上述步骤，参照二维图像的目标局部对象的关键点特征向量的提取方法，本实施例提供的适用于三维图像的目标局部对象的关键点特征向量的提取方法，同样以离散极值点作为关键点特征向量，可以较好地寻获离散极值点，精确性较高，可普适于各种三维医学影像的图像对准，以辅助医生提高诊断的准确性。
如图5所示，在步骤S5中，将提取到的所述参考图像和所述待配准图像的目标局部对象的关键点特征向量进行配准，得到第二步的旋转平移矩阵，采用第二步的旋转平移矩阵对初步配准后的所述参考图像和所述待配准图像进行目标局部对象配准的步骤包括：
步骤S51：从局部参考图像和局部待配准图像的所有关键点特征向量中随机采样多对关键点特征向量。
步骤S52：根据采样的多对关键点特征向量计算局部参考图像和局部待配准图像之间的仿射变换矩阵。
步骤S53：将仿射变换矩阵应用到局部待配准图像的所有关键点特征向量，计算误差值。
步骤S54：当误差值符合预定条件时，将当前仿射变换矩阵作为有效的仿射变换矩阵并作为第二步的旋转平移矩阵。设定的预定条件例如：当误差值小于阈值或达到设定的迭代次数。
步骤S55：当误差值不符合预定条件时，继续进行随机采样，并重新计算仿射变换矩阵，再次计算误差值，直至误差值符合预定条件后将重新计算的仿射变换矩阵作为第二步的旋转平移矩阵。
如图15-图17所示，分别为在同一物理坐标系下配准后的人脑轴状面、冠状面、矢状面的CT图像和MRI图像，配准后的图像的人脑三视图均能较好地对齐，尤其是目标局部对象(神经核团)对齐效果较好，有效避免由于局部组织形变而产生的误配，配准精度较高。
本实施例还提供了一种多模态图像的配准装置，用于执行上述任一技术方案所述的多模态图像的配准方法，所述配准装置包括：
图像获取单元，用于获取参考图像和待配准图像；
图像预处理单元，用于对所述参考图像和所述待配准图像进行预处理；
初步配准单元，基于互信息作为配准测度，得到第一步的旋转平移矩阵，采用第一步的旋转平移矩阵对预处理后的所述参考图像和所述待配准图像进行初步配准；以及
目标局部对象配准单元，用于分别提取所述参考图像和所述待配准图像的目标局部对象的关键点特征向量，以及将提取到的所述参考图像和所述待配准图像的目标局部对象的关键点特征向量进行配准，得到第二步的旋转平移矩阵，采用第二步的旋转平移矩阵对初步配准后的所述参考图像和所述待配准图像进行目标局部对象配准。
本实施例还提供了一种计算机可读存储介质，其上存储有计算机程序，当所述计算机程序被处理器执行时执行上述任一技术方案所述的多模态图像的配准方法的步骤。
综上所述，本发明提供了一种多模态图像的配准方法及其配准装置和计算机可读存储介质，在基于互信息作为配准测度的初步配准后，提取目标局部对象的关键点特征向量进行目标局部对象配准，提高了目标局部对象的配准精度。尤其在医学图像配准中，不同模态的图像会因患者的生理状态变化而发生局部组织形变，在原来坐标上的组织发生偏移，变成其它组织，给医生诊断带来干扰，本发明可选取感兴趣的局部组织作为目标局部对象，提取其关键点特征向量进行配准，从而较好地解决了局部组织形变的问题，提高了配准精度和医生诊断的准确性。
上述描述仅是对本发明一些实施例的描述，并非对本发明范围的任何限定，本发明领域的普通技术人员根据上述揭示内容做的任何变更、修饰，均属于本发明的保护范围。

Claims (13)

1.一种多模态图像的配准方法，其特征在于，包括：
获取参考图像和待配准图像；
对所述参考图像和所述待配准图像进行预处理；
基于互信息作为配准测度，得到第一步的旋转平移矩阵，采用第一步的旋转平移矩阵对预处理后的所述参考图像和所述待配准图像进行初步配准；
分别提取所述参考图像和所述待配准图像的目标局部对象的关键点特征向量；以及
将提取到的所述参考图像和所述待配准图像的目标局部对象的关键点特征向量进行配准，得到第二步的旋转平移矩阵，采用第二步的旋转平移矩阵对初步配准后的所述参考图像和所述待配准图像进行目标局部对象配准。
2.如权利要求1所述的多模态图像的配准方法，其特征在于，对所述参考图像和/或所述待配准图像进行预处理的步骤包括：
采用高斯平滑滤波器对所述参考图像和/或所述待配准图像进行平滑处理。
3.如权利要求1所述的多模态图像的配准方法，其特征在于，基于互信息作为配准测度，得到第一步的旋转平移矩阵，采用第一步的旋转平移矩阵对预处理后的所述参考图像和所述待配准图像进行初步配准的步骤包括：
分别将所述参考图像和所述待配准图像所处的像素坐标系变换为物理坐标系，所述物理坐标系为当前物理坐标系；
基于图像灰度，分别计算所述参考图像和所述待配准图像的质心，并计算所述参考图像和所述待配准图像的质心之间的平移矩阵，且将所述平移矩阵作为初始矩阵；
基于初始矩阵，通过插值运算将所述待配准图像从所述当前物理坐标系变换到新物理坐标系中；
统计在所述当前物理坐标系下的所述参考图像和在所述新物理坐标系下的所述待配准图像的重叠区域的对应像素对的不同灰度组合出现的频度，得到联合直方图，并得到所述参考图像和所述待配准图像的联合熵；以及
根据联合熵计算互信息度量，进行随机梯度迭代，得到第一步的旋转平移矩阵，采用第一步的旋转平移矩阵对所述参考图像和所述待配准图像进行初步配准。
4.如权利要求3所述的多模态图像的配准方法，其特征在于，在计算互信息度量前，对联合熵进行归一化处理。
5.如权利要求3或4所述的多模态图像的配准方法，其特征在于，在将所述参考图像和所述待配准图像所处的像素坐标系变换为物理坐标系的步骤中，利用标签信息将所述参考图像和所述待配准图像所处的像素坐标系变换为物理坐标系。
6.如权利要求5所述的多模态图像的配准方法，其特征在于，所述标签信息包括：单个像素与毫米转换的间距值、像素坐标系原点与物理坐标系原点间的偏移值、和像素坐标系的方向值。
7.如权利要求1所述的多模态图像的配准方法，其特征在于，分别提取所述参考图像和所述待配准图像的目标局部对象的关键点特征向量的步骤包括：
对所述参考图像和所述待配准图像的目标局部对象选取感兴趣区域，分别得到局部参考图像和局部待配准图像；
寻找所述局部参考图像和所述局部待配准图像中的离散极值点；
将每个离散极值点置于多象限坐标系中，分别以每个离散极值点为中心，在所述局部参考图像和所述局部待配准图像中划分出对应的图像块，并将图像块分割为多个图像单元；
计算每个图像单元内各像素点在物理坐标系各坐标轴的梯度，并将每个图像单元内所有像素点的方向和梯度幅值分别累加到多象限坐标系的各象限上，形成梯度直方图；以及
根据梯度直方图，得到每个离散极值点所对应的特征向量作为关键点特征向量。
8.如权利要求7所述的多模态图像的配准方法，其特征在于，寻找所述局部参考图像和所述局部待配准图像中的离散极值点的步骤包括：
选取具有不同标准偏差的高斯分布作为不同尺度，分别对所述局部参考图像和所述局部待配准图像进行高斯平滑，并计算出高斯差分金字塔；以及
在高斯差分金字塔中遍历所有像素点的特定邻域，寻找特定邻域内灰度值最大的点作为离散极值点。
9.如权利要求8所述的多模态图像的配准方法，其特征在于，所述配准方法用于二维图像的配准，所述特定邻域为像素点的上下左右前后六邻域或所有二十六邻域，所述多象限坐标系为由二维平面的360度均分为8块～32块象限的平面坐标系，所述图像块为以离散极值点为中心，以对应高斯分布的标准偏差的3-6倍的值为半径形成的圆，所述图像单元为所述图像块的内接最大正方形均分成的M块正方形中的除去离散极值点所处的中间一块正方形的其余M-1块正方形，其中，M为大于一的正整数的平方数。
10.如权利要求8所述的多模态图像的配准方法，其特征在于，所述配准方法用于三维图像的配准，所述特定邻域为像素点的上下左右前后六邻域或所有二十六邻域，所述多象限坐标系为32球面的球坐标系，所述图像块为以离散极值点为中心，以对应高斯分布的标准偏差的3-6倍的值为半径形成的球，所述图像单元为所述图像块的内接最大正方体均分成的M块正方体中的除去离散极值点所处的中间一块正方体的其余M-1块正方体，其中，M为大于一的正整数的立方数。
11.如权利要求7-10中任一项所述的多模态图像的配准方法，其特征在于，将提取到的所述参考图像和所述待配准图像的目标局部对象的关键点特征向量进行配准，得到第二步的旋转平移矩阵，采用第二步的旋转平移矩阵对初步配准后的所述参考图像和所述待配准图像进行目标局部对象配准的步骤包括：
从局部参考图像和局部待配准图像的所有关键点特征向量中随机采样多对关键点特征向量；
根据采样的多对关键点特征向量计算局部参考图像和局部待配准图像之间的仿射变换矩阵；
将仿射变换矩阵应用到局部待配准图像的所有关键点特征向量，计算误差值；
当误差值符合预定条件时，将当前仿射变换矩阵作为有效的仿射变换矩阵并作为第二步的旋转平移矩阵；
当误差值不符合预定条件时，继续进行随机采样，并重新计算仿射变换矩阵，再次计算误差值，直至误差值符合预定条件后将重新计算的仿射变换矩阵作为第二步的旋转平移矩阵。
12.一种多模态图像的配准装置，其特征在于，包括：
图像获取单元，用于获取参考图像和待配准图像；
图像预处理单元，用于对所述参考图像和所述待配准图像进行预处理；
初步配准单元，基于互信息作为配准测度，得到第一步的旋转平移矩阵，采用第一步的旋转平移矩阵对预处理后的所述参考图像和所述待配准图像进行初步配准；以及
目标局部对象配准单元，用于分别提取所述参考图像和所述待配准图像的目标局部对象的关键点特征向量，以及将提取到的所述参考图像和所述待配准图像的目标局部对象的关键点特征向量进行配准，得到第二步的旋转平移矩阵，采用第二步的旋转平移矩阵对初步配准后的所述参考图像和所述待配准图像进行目标局部对象配准。
13.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，当所述计算机程序被处理器执行时执行如权利要求1-11中任一项所述的多模态图像的配准方法的步骤。
CN202010609691.0A 2020-06-29 2020-06-29 多模态图像的配准方法及其配准装置和计算机可读存储介质 Pending CN113870324A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010609691.0A CN113870324A (zh) 2020-06-29 2020-06-29 多模态图像的配准方法及其配准装置和计算机可读存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010609691.0A CN113870324A (zh) 2020-06-29 2020-06-29 多模态图像的配准方法及其配准装置和计算机可读存储介质

Publications (1)

Publication Number Publication Date
CN113870324A true CN113870324A (zh) 2021-12-31

Family

ID=78981089

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010609691.0A Pending CN113870324A (zh) 2020-06-29 2020-06-29 多模态图像的配准方法及其配准装置和计算机可读存储介质

Country Status (1)

Country Link
CN (1) CN113870324A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114723794A (zh) * 2022-04-12 2022-07-08 南京雷电信息技术有限公司 基于lsd直线检测的sar图像配准方法

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114723794A (zh) * 2022-04-12 2022-07-08 南京雷电信息技术有限公司 基于lsd直线检测的sar图像配准方法

Similar Documents

Publication Publication Date Title
Ferrante et al. Slice-to-volume medical image registration: A survey
US11109773B2 (en) Treating patients with TTFields with the electrode positions optimized using deformable templates
EP2462560B1 (en) Apparatus and method for registering two medical images
EP3316784B1 (en) Fiducial markers, systems, and methods of registration
US8768022B2 (en) Apparatus and methods of compensating for organ deformation, registration of internal structures to images, and applications of same
US8666128B2 (en) Methods, systems, and computer readable media for mapping regions in a model of an object comprising an anatomical structure from one image data set to images used in a diagnostic or therapeutic intervention
Zhang et al. Multi‐needle localization with attention U‐net in US‐guided HDR prostate brachytherapy
Ibragimov et al. Segmentation of tongue muscles from super-resolution magnetic resonance images
Lu et al. Precise segmentation of multiple organs in CT volumes using learning-based approach and information theory
Vásquez Osorio et al. Accurate CT/MR vessel‐guided nonrigid registration of largely deformed livers
WO2008141293A2 (en) Image segmentation system and method
Yavariabdi et al. Mapping and characterizing endometrial implants by registering 2D transvaginal ultrasound to 3D pelvic magnetic resonance images
Liu et al. Automatic localization of the anterior commissure, posterior commissure, and midsagittal plane in MRI scans using regression forests
US10945709B2 (en) Systems, methods and computer readable storage media storing instructions for image-guided interventions based on patient-specific models
CN113870324A (zh) 多模态图像的配准方法及其配准装置和计算机可读存储介质
Mitra et al. A thin-plate spline based multimodal prostate registration with optimal correspondences
Alirr et al. Automatic liver segmentation from ct scans using intensity analysis and level-set active contours
Schwing et al. Reliable extraction of the mid-sagittal plane in 3D brain MRI via hierarchical landmark detection
CN115116113A (zh) 一种光学导航方法
Little et al. The registration of multiple medical images acquired from a single subject: why, how, what next?
Hu Registration of magnetic resonance and ultrasound images for guiding prostate cancer interventions
Liu et al. Multi-modal learning-based pre-operative targeting in deep brain stimulation procedures
Masoumi Inter-Contrast and Inter-Modal Medical Image Registrations: From Traditional Energy-Based to Deep Learning Methods
Kang et al. Non-invasive Face Registration for Surgical Navigation.
CN117796905A (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