CN116612166A - 一种多模态影像的配准融合算法 - Google Patents
一种多模态影像的配准融合算法 Download PDFInfo
- Publication number
- CN116612166A CN116612166A CN202310541933.0A CN202310541933A CN116612166A CN 116612166 A CN116612166 A CN 116612166A CN 202310541933 A CN202310541933 A CN 202310541933A CN 116612166 A CN116612166 A CN 116612166A
- Authority
- CN
- China
- Prior art keywords
- dimensional
- point
- points
- image
- vector
- 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
Links
- 230000004927 fusion Effects 0.000 title claims abstract description 27
- 210000004204 blood vessel Anatomy 0.000 claims abstract description 30
- 239000013598 vector Substances 0.000 claims description 54
- 239000011159 matrix material Substances 0.000 claims description 38
- 238000000034 method Methods 0.000 claims description 21
- 230000009466 transformation Effects 0.000 claims description 10
- 238000000605 extraction Methods 0.000 claims description 9
- 238000013519 translation Methods 0.000 claims description 8
- 238000004364 calculation method Methods 0.000 claims description 7
- 230000008569 process Effects 0.000 claims description 4
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 3
- 101001022170 Homo sapiens FYN-binding protein 2 Proteins 0.000 claims description 3
- 101000990990 Homo sapiens Midkine Proteins 0.000 claims description 3
- 102100030335 Midkine Human genes 0.000 claims description 3
- 238000000354 decomposition reaction Methods 0.000 claims description 3
- 238000009792 diffusion process Methods 0.000 claims description 3
- 238000006073 displacement reaction Methods 0.000 claims description 3
- 230000000694 effects Effects 0.000 claims description 3
- 238000001914 filtration Methods 0.000 claims description 3
- 238000003709 image segmentation Methods 0.000 claims description 3
- 239000000203 mixture Substances 0.000 claims description 3
- 238000012847 principal component analysis method Methods 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 238000012216 screening Methods 0.000 claims description 3
- 238000010845 search algorithm Methods 0.000 claims description 3
- 238000003384 imaging method Methods 0.000 abstract description 7
- 238000003745 diagnosis Methods 0.000 abstract description 6
- 230000002159 abnormal effect Effects 0.000 abstract description 3
- 208000014644 Brain disease Diseases 0.000 abstract description 2
- 238000002059 diagnostic imaging Methods 0.000 abstract description 2
- 210000003484 anatomy Anatomy 0.000 description 3
- 230000005855 radiation Effects 0.000 description 2
- 230000002792 vascular Effects 0.000 description 2
- 239000000654 additive Substances 0.000 description 1
- 230000000996 additive effect Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000002675 image-guided surgery Methods 0.000 description 1
- 238000003780 insertion Methods 0.000 description 1
- 230000037431 insertion Effects 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 210000000056 organ Anatomy 0.000 description 1
- 238000007500 overflow downdraw method Methods 0.000 description 1
- 238000010882 preoperative diagnosis Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000002604 ultrasonography Methods 0.000 description 1
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
- G06T7/337—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods involving reference images or patches
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
- G06V10/77—Processing image or video features in feature spaces; using data integration or data reduction, e.g. principal component analysis [PCA] or independent component analysis [ICA] or self-organising maps [SOM]; Blind source separation
- G06V10/80—Fusion, i.e. combining data from various sources at the sensor level, preprocessing level, feature extraction level or classification level
- G06V10/809—Fusion, i.e. combining data from various sources at the sensor level, preprocessing level, feature extraction level or classification level of classification results, e.g. where the classifiers operate on the same input data
- G06V10/811—Fusion, i.e. combining data from various sources at the sensor level, preprocessing level, feature extraction level or classification level of classification results, e.g. where the classifiers operate on the same input data the classifiers operating on different input data, e.g. multi-modal recognition
-
- 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/10028—Range image; Depth image; 3D point clouds
-
- 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/20—Special algorithmic details
- G06T2207/20212—Image combination
- G06T2207/20221—Image fusion; Image merging
-
- 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)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Artificial Intelligence (AREA)
- Health & Medical Sciences (AREA)
- Computing Systems (AREA)
- Databases & Information Systems (AREA)
- Evolutionary Computation (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Software Systems (AREA)
- Multimedia (AREA)
- Magnetic Resonance Imaging Apparatus (AREA)
Abstract
本发明适用于医学成像技术领域,提供了一种多模态影像的配准融合算法,包括以下步骤:步骤一、将术前三维HR‑MRI影像和术中获取的术中三维DSA影像完成多模态配准,并得到配准后的三维HR‑MRI影像;步骤二、配准后的三维HR‑MRI影像与术中二维DSA影像进行实时融合,并完成配准后的三维HR‑MRI影像与术中二维DSA影像融合;步骤三、依据术中形变,再对术前三维HR‑MRI影像进行变形。应用于神经介入医疗器械治疗脑部疾病的场景下的多模态的血管的配准和融合是具有切实的实际意义,不仅可以解决降低单一成像方式不足以提供诊断所需的所有正常和异常结构的解剖学和功能信息的问题,而且极大的减少了医生手术的时长。
Description
技术领域
本发明属于医学成像技术领域,尤其涉及一种多模态影像的配准融合算法。
背景技术
介入手术涉及计算机视觉、生物医学、影像学、自动控制等多个学科,是一个跨学科的研究方向,它通过多种医学影像信息的综合应用进行术前诊断、病情分析、规划手术路径,术中定位病灶空间位置、实时跟踪手术器械和调整手术器械空间位置,达到精准诊断和精准治疗的日的。
通常,术前通过诸如MRI、CT和MRA等高分辨率三维扫描,获得所需的感兴趣解剖区域的图像,这一类图像具有分辨率高、空间信息多的特点,能够更好地反映人体组织结构和生理信息,但同时成像时间长、成本高、辐射剂量相对较高,不利于在手术环境下成像。
而术中采用的数据是二维超声、X光图像和光学图像这一类二维图像,具有成像快、辐射低的特点,能够适应手术环境,但分辨率较低,想要获取精确完备的病灶位置信息和纹理信息相对困难。
术中单个图像的信息不足以进行准确的诊断,因此需要术前三维图像以显示更高维的信息,即:在图像引导手术中,通过在相同的组织或器官中对相应信息进行比对,将术前图像和术中图像映射到相同的坐标系中,使得解剖结构保持一致,也就是术前3D图像与术中2D图像的融合。但在融合中,实时二维DSA影像会因导丝导管的插入出现血管变形,导致融合精度出现下降,为解决这一问题,将术前影像进行变形预测使其可以为术中导丝导管的插入提供导航。
发明内容
本发明实施例的目的在于提供一种多模态影像的配准融合算法,旨在解决在神经介入手术中,针对单一成像方式不足以提供诊断所需的所有正常和异常结构的解剖学和功能信息以及实时二维的DSA图像作为不足以作为导航来辅助插入导丝导管的问题。
本发明是这样实现的,一种多模态影像的配准融合算法,包括以下步骤:
步骤一、将术前三维HR-MRI影像和术中获取的术中三维DSA影像完成多模态配准,并得到配准后的三维HR-MRI影像;
步骤二、配准后的三维HR-MRI影像与术中二维DSA影像进行实时融合,并完成配准后的三维HR-MRI影像与术中二维DSA影像融合;
步骤三、依据术中形变,再对术前三维HR-MRI影像进行变形。
进一步的技术方案,所述步骤一中多模态配准具体操作为:
1)提取术前三维HR-MRI影像和术中三维DSA影像的点的特征
求解法向量,提取目标点周围半径为R的球体范围内所有点的高斯分布的协方差矩阵,目标点领域内的点云表面是良定义的话,它们就可以近似为一个平面,并且协方差的特征值中仅有一个特征值远小于另外两个,这个最小的特征值对应的特征向量则被认为是这个点领域点云平面的法向量,通过数学公式表达如下:
2)寻找匹配点对
在查找匹配点对的时候利用KDtree查找近邻,匹配关系的确定需要经过如下条件的筛选,如果满足以下条件中的至少一个就去掉组点对;
a.两个点的距离超过阈值,
b.点对中的两个点至少有一个没有良定义的法向量;
c.两个点曲率对数之差绝对值超过阈值,||logσi-logσj||>εσ;
d.两点的法向量之间的角度超过阈值,
3)根据匹配点对计算两组点之间的变换
在得到了两组点云之间的匹配点对之后,通过迭代最小二乘法求取它们之间的变换,给定一个点,得到其坐标、表面曲率、法向量以及协方差,记是带法向量的点/>T代表两组点之间的旋转矩阵,其中旋转部分为矩阵R,平移部分为向量t,/>运算作用在一个带法向量的点/>如式(3):
当给定一对匹配点和当前的变换T,则误差函数eij(T)是一个6维向量,如式(4)所示:
所有点对组成的目标函数表示为式(5):
其中,是一个6X6的信息矩阵,理想情况下,要一个信息矩阵来旋转匹配的两组点云,使得点对的法向量对齐,而且惩罚法向量方向上的投影距离忽略切向的投影误差,最后将位移和法向量部分分开,选择构造/>如式(6):
其中,是点Pi c周围表面信息矩阵,/>是法向量的信息矩阵,为表面曲率足够小,设/>如式(7):
利用nicp算法所使用增量扰动的局部参数化来最小化目标函数,增量表示为ΔT=(Δtx,Δty,Δtz,Δqx,Δqy,Δqz)T,其包含平移向量和旋转单位四元数的虚部,通过在每次迭代计算中使用阻高斯牛顿算法(8):
其中H是近似海森矩阵,Jij雅克比计算;通过上式求出增量则更新,T←ΔT+T;
进一步的技术方案,所述步骤二中配准后的三维HR-MRI影像与术中二维DSA影像进行实时融合,并完成配准后的三维HR-MRI影像与术中二维DSA影像融合的具体操作为:
1)按照血管的曲率完成三维特征点的提取:
a.对于每一个点pi建立一个局部坐标系,并对所有的点设立一个搜索的半径r;
b.对任意点pj的权值为式(10):
c.计算每一个点pi的协方差矩阵并进行特征值分解:
其中,为了排除线和面的可能性,特征值应满足
d.描述三维的空间的点的弯曲程度,满足式(12)则被选为特征点:
2)自适应的给出初始三维模型的投影矩阵;
给出血管的大致走向并旋转三维模型,利用PCA主元分析法获得点云的三个主方向,获取质心,计算协方差,获得协方差矩阵,求取协方差矩阵的特征值和特长向量,特征向量即为主方向,并利用之前获得的主方向和质心,将输入点云转换至原点,且主方向与坐标系方向重回,建立变换到原点的点云的大致走向;
3)将术中二维DSA影像进行图像分割,提取二维特征点
使用传统的AKAZE算法提取二维图像的特征点,它是KAZE特征提取算法的加速版本;通过正则化PM方程与AOS方法来求解非线性扩散,从而得到尺度空间的每一层;采样的方法通过对每一层实现候选点的定位与过滤以实现关键点的提取;
然后再使用与SURF求解方向角度类似的方法实现旋转不变性特征,最终生成其二维特征点;
4)完成三维-二维DSA影像特征点的匹配;
使用全局搜索算法完成匹配关系的确定,并惩罚旋转矩阵完成匹配点关系的确定;
进一步的技术方案,所述步骤三中依据术中形变,再对术前三维HR-MRI影像进行变形具体为将已经配准融合后的术前血管模型进行变形,具体操作如下:
按照步骤二已经确定的匹配关系,计算三维和二维特征点的最小二乘公式,并依据其的梯度方向,进行血管的变形;
该算法驱动整个变形配准过程的差分度量惩罚了来自变形前的三维顶点V和来自变形后的对应的三维点V’之间的距离,也可以被描述为ARAP能量函数Earap相对于血管节点的旋转变化Ri的最小化,以便得到p’的估计,即:
其中,表示顶点vi的单环邻域,权值wij是通过曲率进行计算的;
同时算法还考虑到导丝导管对血管的拉直作用,使用血管中心线作为限制条件,限制变形公式:
本发明实施例提供的一种多模态影像的配准融合算法,应用于神经介入医疗器械治疗脑部疾病的场景下的多模态的血管的配准和融合是具有切实的实际意义,不仅可以解决降低单一成像方式不足以提供诊断所需的所有正常和异常结构的解剖学和功能信息的问题,而且极大的减少了医生手术的时长,并且能准确配准术前和术中的血管,神经介入手术是借助放射学及先进的影像引导技术的微创手术;本方法是考虑血管变形的多模态医疗图像配准融合方法,将术前体积数据与术中透视数据配准融合。通过将术前血管解剖结构叠加到透视图像上,将变形后的术前图像可以作为支持手术过程中导航的路线图,进一步提高多模态医疗图像配准融合的精确度。
附图说明
图1为本发明实施例提供的一种多模态影像的配准融合算法的整体流程图;
图2为本发明实施例提供的图1中多模态配准的流程图;
图3为本发明实施例提供的图1中三维HR-MRI影像与术中二维DSA影像融合的流程图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白,以下结合附图及实施例,对本发明进行进一步详细说明。应当理解,此处所描述的具体实施例仅仅用以解释本发明,并不用于限定本发明。
以下结合具体实施例对本发明的具体实现进行详细描述。
如图1-图3所示,为本发明一个实施例提供的一种多模态影像的配准融合算法,为保证配准的精度,将选取对血管壁厚、血管曲面凸凹特征不敏感的血管特征值作为配准特征值,后对不同影像下的同段血管进行三维重建,利用血管特征值以及局部点云配准算法对不同模态下的同段血管进行配准与融合,为术前以及术中的医生提供更完整、更全面、更精准的信息,提高诊断效率与治疗精准度实现方案如下,包括以下步骤:;
步骤一、将术前三维HR-MRI影像和术中获取的术中三维DSA影像完成多模态配准,并得到配准后的三维HR-MRI影像;将更多考虑血管模型的特点,将每个点的欧几里得坐标都可以用曲面法线来增强,采用最小化每个对应点对及其法线的马氏距离寻找最合适的匹配关系;根据得到得两点云的匹配情况,构建最小二乘方程,求解旋转和平移矩阵;
多模态配准具体操作为:
1)提取术前三维HR-MRI影像和术中三维DSA影像的点的特征
求解法向量,提取目标点周围半径为R的球体范围内所有点的高斯分布的协方差矩阵,目标点领域内的点云表面是良定义的话,它们就可以近似为一个平面,并且协方差的特征值中仅有一个特征值远小于另外两个,这个最小的特征值对应的特征向量则被认为是这个点领域点云平面的法向量,通过数学公式表达如下:
2)寻找匹配点对
在查找匹配点对的时候利用KDtree查找近邻,匹配关系的确定需要经过如下条件的筛选,如果满足以下条件中的至少一个就去掉组点对;
a.两个点的距离超过阈值,
b.点对中的两个点至少有一个没有良定义的法向量;
c.两个点曲率对数之差绝对值超过阈值,||logσi-logσj||>εσ;
d.两点的法向量之间的角度超过阈值,
3)根据匹配点对计算两组点之间的变换
在得到了两组点云之间的匹配点对之后,通过迭代最小二乘法求取它们之间的变换,给定一个点,得到其坐标、表面曲率、法向量以及协方差,记是带法向量的点/>T代表两组点之间的旋转矩阵,其中旋转部分为矩阵R,平移部分为向量t,/>运算作用在一个带法向量的点/>如式(3):
当给定一对匹配点和当前的变换T,则误差函数eij(T)是一个6维向量,如式(4)所示:
所有点对组成的目标函数表示为式(5):
其中,是一个6X6的信息矩阵,理想情况下,要一个信息矩阵来旋转匹配的两组点云,使得点对的法向量对齐,而且惩罚法向量方向上的投影距离忽略切向的投影误差(类似于point to plane ICP),最后将位移和法向量部分分开,选择构造/>如式(6):
其中,是点Pi c周围表面信息矩阵,/>是法向量的信息矩阵,为表面曲率足够小,设/>如式(7):
利用nicp算法所使用增量扰动的局部参数化来最小化目标函数,增量表示为ΔT=(Δtx,Δty,Δtz,Δqx,Δqy,Δqz)T,其包含平移向量和旋转单位四元数的虚部,通过在每次迭代计算中使用阻高斯牛顿算法(8):
其中H是近似海森矩阵,Jij雅克比计算;通过上式求出增量则更新,T←ΔT+T;
步骤二、配准后的三维HR-MRI影像与术中二维DSA影像进行实时融合,并完成配准后的三维HR-MRI影像与术中二维DSA影像融合;首先提取三维-DSA血管中心线作为三维点集;在一些情况下,也可以手动编辑中心线,以去除与细血管相对应的分支,对二维DSA图像上的血管提取中心线点,寻找三维点集的旋转和平移,使其投影能够与二维点集相匹配;
配准后的三维HR-MRI影像与术中二维DSA影像进行实时融合,并完成配准后的三维HR-MRI影像与术中二维DSA影像融合的具体操作为:
1)按照血管的曲率完成三维特征点的提取:
a.对于每一个点pi建立一个局部坐标系,并对所有的点设立一个搜索的半径r;
b.对任意点pj的权值为式(10):
c.计算每一个点pi的协方差矩阵并进行特征值分解:
其中,为了排除线和面的可能性,特征值应满足
d.描述三维的空间的点的弯曲程度,满足式(12)则被选为特征点:
2)自适应的给出初始三维模型的投影矩阵:
给出血管的大致走向并旋转三维模型,利用PCA主元分析法获得点云的三个主方向,获取质心,计算协方差,获得协方差矩阵,求取协方差矩阵的特征值和特长向量,特征向量即为主方向,并利用之前获得的主方向和质心,将输入点云转换至原点,且主方向与坐标系方向重回,建立变换到原点的点云的大致走向;
3)将术中二维DSA影像进行图像分割,提取二维特征点:
使用传统的AKAZE算法提取二维图像的特征点,该算法具备很好的具有尺度不变性与旋转不变性,它是KAZE特征提取算法的加速版本;通过正则化PM方程与AOS(加性算子分裂)方法来求解非线性扩散,从而得到尺度空间的每一层;采样的方法通过对每一层实现候选点的定位与过滤以实现关键点的提取;
然后再使用与SURF求解方向角度类似的方法实现旋转不变性特征,最终生成其二维特征点;
4)完成三维-二维DSA影像特征点的匹配:
使用全局搜索算法完成匹配关系的确定,并惩罚旋转矩阵完成匹配点关系的确定:
步骤三、依据术中形变,再对术前三维HR-MRI影像进行变形具体为将已经配准融合后的术前血管模型进行变形,具体操作如下:
按照步骤二已经确定的匹配关系,计算三维和二维特征点的最小二乘公式,并依据其的梯度方向,进行血管的变形;
该算法驱动整个变形配准过程的差分度量惩罚了来自变形前的三维顶点V和来自变形后的对应的三维点V’之间的距离,也可以被描述为ARAP能量函数Earap相对于血管节点的旋转变化Ri的最小化,以便得到p’的估计,即:
其中,表示顶点vi的单环邻域,权值wij是通过曲率进行计算的;
同时算法还考虑到导丝导管对血管的拉直作用,使用血管中心线作为限制条件,限制变形公式:
因为曲率是用来反映几何体的弯曲程度,定性的看,弯曲的越厉害,该部分的曲率越大。我们认为血管弯曲程度上越大的地方越容易被矫直,那么这个面片更应该被关注。类似注意力机制,只不过我们搜索的是面片的曲率,而其考虑的是图像特征。
利用最小二乘法寻找每个三角面尽可能刚性的变换,并求出这个其变形的位置。最小化能量公式有两个代求量(R和P’),使用全局与局部的优化方法,先利用SVD求解单个三角面的旋转矩阵代入能量公式对其求导得出预测的位置。
以上所述仅为本发明的较佳实施例而已,并不用以限制本发明,凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等,均应包含在本发明的保护范围之内。
Claims (4)
1.一种多模态影像的配准融合算法,其特征在于,包括以下步骤:
步骤一、将术前三维HR-MRI影像和术中获取的术中三维DSA影像完成多模态配准,并得到配准后的三维HR-MRI影像;
步骤二、配准后的三维HR-MRI影像与术中二维DSA影像进行实时融合,并完成配准后的三维HR-MRI影像与术中二维DSA影像融合;
步骤三、依据术中形变,再对术前三维HR-MRI影像进行变形。
2.根据权利要求1所述的多模态影像的配准融合算法,其特征在于,所述步骤一中多模态配准具体操作为:
1)提取术前三维HR-MRI影像和术中三维DSA影像的点的特征
求解法向量,提取目标点周围半径为R的球体范围内所有点的高斯分布的协方差矩阵,目标点领域内的点云表面是良定义的话,它们就可以近似为一个平面,并且协方差的特征值中仅有一个特征值远小于另外两个,这个最小的特征值对应的特征向量则被认为是这个点领域点云平面的法向量,通过数学公式表达如下:
2)寻找匹配点对
在查找匹配点对的时候利用KDtree查找近邻,匹配关系的确定需要经过如下条件的筛选,如果满足以下条件中的至少一个就去掉组点对;
a.两个点的距离超过阈值,
b.点对中的两个点至少有一个没有良定义的法向量;
c.两个点曲率对数之差绝对值超过阈值,||logσi-logσj||>εσ;
d.两点的法向量之间的角度超过阈值,
3)根据匹配点对计算两组点之间的变换
在得到了两组点云之间的匹配点对之后,通过迭代最小二乘法求取它们之间的变换,给定一个点,得到其坐标、表面曲率、法向量以及协方差,记是带法向量的点/>T代表两组点之间的旋转矩阵,其中旋转部分为矩阵R,平移部分为向量t,/>运算作用在一个带法向量的点/>如式(3):
当给定一对匹配点和当前的变换T,则误差函数eij(T)是一个6维向量,如式(4)所示:
所有点对组成的目标函数表示为式(5):
其中,是一个6X6的信息矩阵,理想情况下,要一个信息矩阵来旋转匹配的两组点云,使得点对的法向量对齐,而且惩罚法向量方向上的投影距离忽略切向的投影误差,最后将位移和法向量部分分开,选择构造/>如式(6):
其中,是点Pi c周围表面信息矩阵,/>是法向量的信息矩阵,为表面曲率足够小,设如式(7):
利用nicp算法所使用增量扰动的局部参数化来最小化目标函数,增量表示为ΔT=(Δtx,Δty,Δtz,Δqx,Δqy,Δqz)T,其包含平移向量和旋转单位四元数的虚部,通过在每次迭代计算中使用阻高斯牛顿算法(8):
其中H是近似海森矩阵,Jij雅克比计算;通过上式求出增量则更新,T←ΔT+T;
3.根据权利要求1所述的多模态影像的配准融合算法,其特征在于,所述步骤二中配准后的三维HR-MRI影像与术中二维DSA影像进行实时融合,并完成配准后的三维HR-MRI影像与术中二维DSA影像融合的具体操作为:
1)按照血管的曲率完成三维特征点的提取:
a.对于每一个点pi建立一个局部坐标系,并对所有的点设立一个搜索的半径r;
b.对任意点pj的权值为式(10):
c.计算每一个点pi的协方差矩阵并进行特征值分解:
其中,为了排除线和面的可能性,特征值应满足
d.描述三维的空间的点的弯曲程度,满足式(12)则被选为特征点:
2)自适应的给出初始三维模型的投影矩阵;
给出血管的大致走向并旋转三维模型,利用PCA主元分析法获得点云的三个主方向,获取质心,计算协方差,获得协方差矩阵,求取协方差矩阵的特征值和特长向量,特征向量即为主方向,并利用之前获得的主方向和质心,将输入点云转换至原点,且主方向与坐标系方向重回,建立变换到原点的点云的大致走向;
3)将术中二维DSA影像进行图像分割,提取二维特征点;
使用传统的AKAZE算法提取二维图像的特征点,它是KAZE特征提取算法的加速版本;通过正则化PM方程与AOS方法来求解非线性扩散,从而得到尺度空间的每一层;采样的方法通过对每一层实现候选点的定位与过滤以实现关键点的提取;
然后再使用与SURF求解方向角度类似的方法实现旋转不变性特征,最终生成其二维特征点;
4)完成三维-二维DSA影像特征点的匹配;
使用全局搜索算法完成匹配关系的确定,并惩罚旋转矩阵完成匹配点关系的确定;
4.根据权利要求3所述的多模态影像的配准融合算法,其特征在于,所述步骤三中依据术中形变,再对术前三维HR-MRI影像进行变形具体为将已经配准融合后的术前血管模型进行变形,具体操作如下:
按照步骤二已经确定的匹配关系,计算三维和二维特征点的最小二乘公式,并依据其的梯度方向,进行血管的变形;
该算法驱动整个变形配准过程的差分度量惩罚了来自变形前的三维顶点V和来自变形后的对应的三维点V’之间的距离,也可以被描述为ARAP能量函数Earap相对于血管节点的旋转变化Ri的最小化,以便得到p’的估计,即:
其中,表示顶点wi的单环邻域,权值wij是通过曲率进行计算的;
同时算法还考虑到导丝导管对血管的拉直作用,使用血管中心线作为限制条件,限制变形公式:
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202310541933.0A CN116612166A (zh) | 2023-05-15 | 2023-05-15 | 一种多模态影像的配准融合算法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202310541933.0A CN116612166A (zh) | 2023-05-15 | 2023-05-15 | 一种多模态影像的配准融合算法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN116612166A true CN116612166A (zh) | 2023-08-18 |
Family
ID=87681091
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202310541933.0A Pending CN116612166A (zh) | 2023-05-15 | 2023-05-15 | 一种多模态影像的配准融合算法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN116612166A (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116862930A (zh) * | 2023-09-04 | 2023-10-10 | 首都医科大学附属北京天坛医院 | 适用于多模态的脑血管分割方法、装置、设备和存储介质 |
-
2023
- 2023-05-15 CN CN202310541933.0A patent/CN116612166A/zh active Pending
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116862930A (zh) * | 2023-09-04 | 2023-10-10 | 首都医科大学附属北京天坛医院 | 适用于多模态的脑血管分割方法、装置、设备和存储介质 |
CN116862930B (zh) * | 2023-09-04 | 2023-11-28 | 首都医科大学附属北京天坛医院 | 适用于多模态的脑血管分割方法、装置、设备和存储介质 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Ferrante et al. | Slice-to-volume medical image registration: A survey | |
Alam et al. | Medical image registration in image guided surgery: Issues, challenges and research opportunities | |
US11890063B2 (en) | System and methods for a trackerless navigation system | |
Guan et al. | A review of point feature based medical image registration | |
Blackall et al. | Alignment of sparse freehand 3-D ultrasound with preoperative images of the liver using models of respiratory motion and deformation | |
US9687204B2 (en) | Method and system for registration of ultrasound and physiological models to X-ray fluoroscopic images | |
US8271068B2 (en) | Method for dynamic road mapping | |
US9155470B2 (en) | Method and system for model based fusion on pre-operative computed tomography and intra-operative fluoroscopy using transesophageal echocardiography | |
JP7046553B2 (ja) | 撮像装置を備える磁気追跡システムの重ね合わせ法 | |
Che et al. | Ultrasound registration: A review | |
US20060269165A1 (en) | Registration of three dimensional image data with patient in a projection imaging system | |
US20080095422A1 (en) | Alignment method for registering medical images | |
Wein et al. | Automatic bone detection and soft tissue aware ultrasound–CT registration for computer-aided orthopedic surgery | |
Song et al. | Locally rigid, vessel-based registration for laparoscopic liver surgery | |
JP2008546441A (ja) | 第1及び第2画像を比較するためのモデルに基づく弾性画像登録方法 | |
WO2001043070A2 (en) | Method and apparatus for cross modality image registration | |
Wu et al. | Three-dimensional modeling from endoscopic video using geometric constraints via feature positioning | |
CN110301883B (zh) | 用于导航管状网络的基于图像的向导 | |
CN111640143A (zh) | 一种基于PointNet的神经导航快速面配准方法及系统 | |
Cheema et al. | Image-aligned dynamic liver reconstruction using intra-operative field of views for minimal invasive surgery | |
CN115578320A (zh) | 一种骨科手术机器人全自动空间注册方法及系统 | |
CN116612166A (zh) | 一种多模态影像的配准融合算法 | |
CN115830016A (zh) | 医学图像配准模型训练方法及设备 | |
CN111260704A (zh) | 基于启发式树搜索的血管结构3d/2d刚性配准方法及装置 | |
De Silva et al. | Robust 2-D–3-D registration optimization for motion compensation during 3-D TRUS-guided biopsy using learned prostate motion data |
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 |