CN112233155A - 一种2d-3d影像配准算法 - Google Patents
一种2d-3d影像配准算法 Download PDFInfo
- Publication number
- CN112233155A CN112233155A CN202011415198.1A CN202011415198A CN112233155A CN 112233155 A CN112233155 A CN 112233155A CN 202011415198 A CN202011415198 A CN 202011415198A CN 112233155 A CN112233155 A CN 112233155A
- Authority
- CN
- China
- Prior art keywords
- image
- preoperative
- matching
- gradient
- pose
- 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.)
- Granted
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
- 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/10116—X-ray image
-
- 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
Landscapes
- Engineering & Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
本发明公开了一种2D‑3D影像配准算法,包括步骤:(1)在术前扫描患者患处得到术前3D影像;在术中分别对患者相对应患处进行正侧位透视扫描得到正侧位2D透视影像;(2)构建术中拍摄正侧位透视影像设备的空间模型,并根据步骤(1)计算得到术前3D影像在空间模型中的初始位姿;(3)采用步骤(2)构建的空间模型生成术前3D影像的正侧位DRR影像,并将其与其对应的术中2D透视影像进行模板匹配,得到术前3D影像的粗匹配变换位姿;(4)使用梯度下降优化算法进行精准优化匹配,得到术前3D影像的最终变换位姿。本发明通过采用模板匹配在进行粗匹配,然后使用梯度下降优化算法得到最优结果,可以在保证配准速度的同时,得到非常精准的配准结果。
Description
技术领域
本发明涉及手术导航技术领域,尤其涉及一种3D CT影像和2D透视影像的配准算法。
背景技术
3D CT影像和2D透视影像配准算法计算3D CT影像和2D透视影像之间的位姿变换值,在影像导航置入手术中被用于将术前信息(基于CT的手术计划数据)整合到术中X射线投影影像中,以提高到达目标位置的精度。当前的算法实现通常使用3D CT影像生成2D DRR(数字重建放射)影像,计算DRR影像和透视影像之间的相似度,找到CT的某个位姿使生成的DRR影像和透视影像具有最佳的相似度。由于目标变量有6个自由度(其中3个表示旋转,另外3个表示平移),传统的基于模板匹配的相似度评价的算法或基于随机采样搜索的算法并不能得到精准的结果,且非常耗时,影响算法在实际场景的应用。
发明内容
发明目的:本发明针对上述不足,提出了一种快速、准确的2D-3D影像配准算法。
技术方案:
一种2D-3D影像配准算法,包括步骤:
(1)在术前扫描患者患处得到术前3D影像;在术中分别对患者相对应患处进行正侧位透视扫描得到患者患处的正侧位2D透视影像;
(2)构建术中拍摄正侧位透视影像设备的空间模型,并根据步骤(1)计算得到术前3D影像在空间模型中的初始位姿;
(3)采用步骤(2)构建的空间模型生成术前3D影像的正侧位DRR影像,并将其与其对应的2D透视影像进行模板匹配,得到术前3D影像的粗匹配变换位姿;
(4)使用梯度下降优化算法进行精准优化匹配,得到术前3D影像的最终变换位姿。
所述步骤(3)中,位姿包括三维欧拉角和三维平移向量;具体计算得到术前3D影像的粗匹配变换位姿如下:
通过模板匹配计算生成的DRR影像和其对应的2D透视影像的归一化相关匹配值,得到最佳的匹配位置,据此计算得到术前3D影像粗匹配后的平移向量;
对于初始欧拉角,分别以生成的DRR影像的每个分量角度以预设的偏差取值,并分别计算其与对应的2D透视影像之间的归一化相关匹配值,得到最高匹配值对应的欧拉角;
综上计算得到术前3D影像的粗匹配变换位姿。
所述步骤(4)中:
以术前3D影像的粗匹配变换位姿作为待优化变量;
分别计算正侧位DRR影像的梯度图及经过粗匹配变换位姿变换后的3D影像生成的对应DRR影像的梯度图,并以二者之间的相似度值作为评价函数;
使用梯度下降优化算法计算术前3D影像的最终变换位姿。
使用Scharr算子分别计算所述DRR影像的梯度图和术中2D透视影像的梯度图。
所述DRR影像的梯度图和2D透视影像的梯度图之间的相似度值metric计算如下:
使用梯度下降优化算法计算术前3D影像的最终变换位姿具体为:
(41)计算当前值T的梯度,梯度由每个分量T[i]的导数组成;
(42)设定偏差Δ,则T[i]的导数为:
(43)沿梯度方向优化当前值T,计算优化后的变量T[i]'为:
(44)计算震荡系数shakeFactor:
dPT[i]为上次迭代计算得到的PT每一分量的导数,初始值设为0;若shakeFactor<0,则优化步长stepLength调整为之前的一半,并重新确定偏差,重复步骤(42);否则转至步骤(45);
(45)判断优化步长stepLength是否满足终止条件,若满足则停止迭代,输出最终的T;若不满足,则重新确定偏差Δ,返回步骤(42)。
所述步骤(45)中的终止条件为stepLength<0.03。
有益效果:本发明的算法首先采用模板匹配在一定范围内进行粗匹配,然后使用梯度下降优化算法得到最优结果,可以在保证配准速度的同时,得到非常精准的配准结果。
附图说明
图1为本发明的2D-3D影像配准算法框架流程图。
图2为术前CT影像中模拟X光机生成DRR影像的原理图。
图3为在图像x和y方向上进行模板匹配计算的归一化相关匹配值。
具体实施方式
下面结合附图和具体实施例,进一步阐明本发明。
本发明主要应用在人体椎骨的CT透视影像配准中,图1为本发明的流程图,如图1所示,本发明的2D-3D影像配准算法包括步骤:
(1)通过C臂机或CT机在术前扫描患者患处(本实施例为脊椎)得到术前3D影像;对术前3D影像进行自动分割得到椎骨图像,进而得到椎骨的中心点位置(即椎骨中心点在影像坐标系中的坐标)和包围盒(即椎骨的外接长方体);
(2)通过C臂机或CT机在术中分别对患者患处进行正侧位透视扫描,得到患者患处的正侧位术中2D透视影像,并使用透视变换矩阵得到单个椎骨影像,在术中2D透视影像中的椎骨中心点加上标记;
(3)构建术中拍摄正侧位透视影像的C臂机或CT机的空间模型(以X射线发射管中心作为空间模型光心,C臂机成像平面作为空间模型成像平面),根据构建的空间模型中成像平面位姿及步骤(2)标记的术中2D透视影像的椎骨中心点的坐标可以计算出术中2D透视影像中的椎骨中心点在构建的空间模型中的位置坐标;
(4)根据步骤(1)得到的椎骨中心点在影像坐标系中的坐标,可以计算出3D影像在C臂机或CT机空间模型中的初始位姿T 0[r 0 ,t 0],其中r 0是表示旋转的3维欧拉角,t 0是表示平移的3维平移向量;
(5)采用步骤(3)构建的C臂机或CT机空间模型生成术前3D影像正侧位2D DRR影像,生成方法如图2所示:
从X射线源向2D DRR成像平面每个像素的方向采样,由以下公式计算最终的像素值I xy :
式中i是术前3D影像的采样值,作为衰减系数。生成算法需要大量计算,并且是可以并行计算的,所以使用CUDA加速计算;
(6)利用模板匹配进行粗配准:
通过模板匹配计算生成的2D DRR影像和其对应的术中2D透视影像的归一化相关匹配值,结果如图3 所示,最亮的像素点对应最佳的匹配位置,这样就可以得到术前3D影像粗匹配后的平移向t 1;对于初始欧拉角r 1,分别以生成的2D DRR影像的每个分量角度以预设的偏差取值,并分别计算其与对应的术中2D透视影像之间的归一化相关匹配值,得到最高匹配值对应的欧拉角r 1,最终得到粗匹配后的术前3D影像的变换位姿T 1[r 1 ,t 1];
(7)粗匹配可以将平移向量的误差控制在0.5毫米以内(透视影像一个像素对应的长度),下面需要使用梯度下降优化算法进行精准优化匹配,具体优化步骤为:
(71)将粗匹配后的变换位姿T 1[r 1 ,t 1]传递给梯度下降优化算法的初始值T,T就是优化算法的待优化变量;
(72)由当前的变换位姿T得到变换后的3D影像,并采用步骤(5)生成DRR影像,并使用透视变换矩阵得到单个椎骨影像(脊柱是非刚体,所以算法进行单椎骨配准),使用Scharr算子分别计算该DRR影像的梯度图和透视影像的梯度图,再计算两幅梯度图的相似度值metric;
其中,表示正位2D透视影像的梯度图,表示正位DRR影像的梯度图,是正位透视影像的方差值。该方法计算每个对应像素的相似度,然后将所有像素的相似度值累加。显然,透视图和DRR图越相似,越接近0,像素相似度值越接近1,的值越接近梯度的像素数,它们是正相关的;侧位相似度的计算方法类似,最后计算总的相似度值。
本算法框架在计算DRR、梯度图和相似度值都使用CUDA完成,这样不仅可以利用GPU加速计算,还可以减少数据在主机内存和设备内存(GPU显存)间传递所花费的时间,最终只需要将相似度值从显存中拷贝至主机内存。
(73)计算当前6自由度T的梯度,梯度由每个分量T[i]的导数组成;具体为:定义偏差Δ,则T[i]的导数为:
(74)调整优化步长stepLength,具体为计算震荡系数shakeFactor:
dPT[i]为上次迭代计算得到的T(定义为PT)每一分量的导数,初始值设为0。如果shakeFactor<0,则stepLength调整为之前的一半,并重新确定偏差,重复步骤(73);否则转至步骤(75);
(75)沿梯度方向优化当前值T,计算优化后的变量T[i]'为:
(76)判断优化步长stepLength是否满足终止条件(即stepLength<0.03),若满足则停止迭代,输出最终的变换位姿T;若不满足,则重新确定偏差Δ,返回步骤(73)。
以上详细描述了本发明的优选实施方式,但是本发明并不限于上述实施方式中的具体细节,在本发明的技术构思范围内,可以对本发明的技术方案进行多种等同变换(如数量、形状、位置等),这些等同变换均属于本发明的保护范围。
Claims (7)
1.一种2D-3D影像配准算法,其特征在于:包括步骤:
(1)在术前扫描患者患处得到术前3D影像;在术中分别对患者相对应患处进行正侧位透视扫描得到患者患处的正侧位2D透视影像;
(2)构建术中拍摄正侧位透视影像设备的空间模型,并根据步骤(1)计算得到术前3D影像在空间模型中的初始位姿;
(3)采用步骤(2)构建的空间模型生成术前3D影像的正侧位DRR影像,并将其与其对应的2D透视影像进行模板匹配,得到术前3D影像的粗匹配变换位姿;
(4)使用梯度下降优化算法进行精准优化匹配,得到术前3D影像的最终变换位姿。
2.根据权利要求1所述的2D-3D影像配准算法,其特征在于:所述步骤(3)中,位姿包括三维欧拉角和三维平移向量;具体计算得到术前3D影像的粗匹配变换位姿如下:
通过模板匹配计算生成的DRR影像和其对应的2D透视影像的归一化相关匹配值,得到最佳的匹配位置,据此计算得到术前3D影像粗匹配后的平移向量;
对于初始欧拉角,分别以生成的DRR影像的每个分量角度以预设的偏差取值,并分别计算其与对应的2D透视影像之间的归一化相关匹配值,得到最高匹配值对应的欧拉角;
综上计算得到术前3D影像的粗匹配变换位姿。
3.根据权利要求1所述的2D-3D影像配准算法,其特征在于:所述步骤(4)中:
以术前3D影像的粗匹配变换位姿作为待优化变量;
分别计算正侧位DRR影像的梯度图及经过粗匹配变换位姿变换后的3D影像生成的对应DRR影像的梯度图,并以二者之间的相似度值作为评价函数;
使用梯度下降优化算法计算术前3D影像的最终变换位姿。
4.根据权利要求3所述的2D-3D影像配准算法,其特征在于:使用Scharr算子分别计算所述DRR影像的梯度图和2D透视影像的梯度图。
6.根据权利要求3所述的2D-3D影像配准算法,其特征在于:使用梯度下降优化算法计算术前3D影像的最终变换位姿具体为:
(41)计算当前值T的梯度,梯度由每个分量T[i]的导数组成;
(42)设定偏差Δ,则T[i]的导数为:
(43)沿梯度方向优化当前值T,计算优化后的变量T[i]'为:
(44)计算震荡系数shakeFactor:
dPT[i]为上次迭代计算得到的PT每一分量的导数,初始值设为0;若shakeFactor<0,则优化步长stepLength调整为之前的一半,并重新确定偏差,重复步骤(42);否则转至步骤(45);
(45)判断优化步长stepLength是否满足终止条件,若满足则停止迭代,输出最终的T;若不满足,则重新确定偏差Δ,返回步骤(42)。
7.根据权利要求6所述的2D-3D影像配准算法,其特征在于:所述步骤(45)中的终止条件为stepLength<0.03。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011415198.1A CN112233155B (zh) | 2020-12-07 | 2020-12-07 | 一种2d-3d影像配准算法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011415198.1A CN112233155B (zh) | 2020-12-07 | 2020-12-07 | 一种2d-3d影像配准算法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN112233155A true CN112233155A (zh) | 2021-01-15 |
CN112233155B CN112233155B (zh) | 2021-02-26 |
Family
ID=74123566
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202011415198.1A Active CN112233155B (zh) | 2020-12-07 | 2020-12-07 | 一种2d-3d影像配准算法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112233155B (zh) |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113040908A (zh) * | 2021-02-02 | 2021-06-29 | 武汉联影智融医疗科技有限公司 | 手术导航的配准方法、装置、计算机设备和存储介质 |
CN115082534A (zh) * | 2022-07-21 | 2022-09-20 | 杭州三坛医疗科技有限公司 | 双平面图像配准方法、装置及机器人 |
WO2023006021A1 (zh) * | 2021-07-30 | 2023-02-02 | 武汉联影智融医疗科技有限公司 | 一种配准方法和系统 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20130094745A1 (en) * | 2011-09-28 | 2013-04-18 | Siemens Corporation | Non-rigid 2d/3d registration of coronary artery models with live fluoroscopy images |
CN111429491A (zh) * | 2020-03-11 | 2020-07-17 | 上海嘉奥信息科技发展有限公司 | 脊柱术前三维影像与术中二维影像配准的方法与系统 |
-
2020
- 2020-12-07 CN CN202011415198.1A patent/CN112233155B/zh active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20130094745A1 (en) * | 2011-09-28 | 2013-04-18 | Siemens Corporation | Non-rigid 2d/3d registration of coronary artery models with live fluoroscopy images |
CN111429491A (zh) * | 2020-03-11 | 2020-07-17 | 上海嘉奥信息科技发展有限公司 | 脊柱术前三维影像与术中二维影像配准的方法与系统 |
Non-Patent Citations (2)
Title |
---|
共和国之辉: "《最优化五:梯度法(梯度下降法、最优梯度法、共轭梯度法)》", 《HTTPS://BLOG.CSDN.NET/LITTLEEMPEROR/ARTICLE/DETAILS/105073755》 * |
张冉等: "《2D/3D图像配准中的相似性测度和优化算法》", 《激光与红外》 * |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113040908A (zh) * | 2021-02-02 | 2021-06-29 | 武汉联影智融医疗科技有限公司 | 手术导航的配准方法、装置、计算机设备和存储介质 |
WO2023006021A1 (zh) * | 2021-07-30 | 2023-02-02 | 武汉联影智融医疗科技有限公司 | 一种配准方法和系统 |
CN115082534A (zh) * | 2022-07-21 | 2022-09-20 | 杭州三坛医疗科技有限公司 | 双平面图像配准方法、装置及机器人 |
Also Published As
Publication number | Publication date |
---|---|
CN112233155B (zh) | 2021-02-26 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN112233155B (zh) | 一种2d-3d影像配准算法 | |
US11605185B2 (en) | System and method for generating partial surface from volumetric data for registration to surface topology image data | |
JP2022141792A5 (zh) | ||
US6381302B1 (en) | Computer assisted 2D adjustment of stereo X-ray images | |
RU2568635C2 (ru) | Регистрация двумерных/трехмерных изображений на основе признаков | |
US8184886B2 (en) | Deformable 2D-3D registration | |
US20180189966A1 (en) | System and method for guidance of laparoscopic surgical procedures through anatomical model augmentation | |
US20050047544A1 (en) | Apparatus and method for registering 2D radiographic images with images reconstructed from 3D scan data | |
US20150115178A1 (en) | Methods and systems for guiding an emission to a target | |
CN109925054B (zh) | 确定靶点路径的辅助方法、装置和系统、可读存储介质 | |
CN104268849A (zh) | 图像核对装置以及患者定位装置 | |
CN107680688B (zh) | 一种基于3d打印的盆腔仿真微创手术视觉导航验证方法 | |
JP2022534123A (ja) | 画像レジストレーション方法及びそれに関係するモデルトレーニング方法、デバイス、装置 | |
CN112289416B (zh) | 一种引导针置入精度评价方法 | |
US20170270678A1 (en) | Device and method for image registration, and non-transitory recording medium | |
CN109767458A (zh) | 一种半自动分段的顺序优化配准方法 | |
CN109925053B (zh) | 手术路径的确定方法、装置和系统、可读存储介质 | |
Maharjan et al. | A novel visualization system of using augmented reality in knee replacement surgery: Enhanced bidirectional maximum correntropy algorithm | |
CN113963057B (zh) | 成像几何关系标定方法、装置、电子设备以及存储介质 | |
CN112472293B (zh) | 一种术前三维影像与术中透视图像的配准方法 | |
KR102469141B1 (ko) | 의료용 화상 처리 장치, 의료용 화상 처리 방법, 및 프로그램 | |
CN114511597A (zh) | X光图像与ct图像的配准方法 | |
JP7513980B2 (ja) | 医用画像処理装置、治療システム、医用画像処理方法、およびプログラム | |
EP4144298A1 (en) | Object visualisation in x-ray imaging | |
JP7407831B2 (ja) | 介入装置追跡 |
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 | ||
CP03 | Change of name, title or address |
Address after: 210000 building 3, No. 34, Dazhou Road, Yuhuatai District, Nanjing, Jiangsu Province Patentee after: Tuodao Medical Technology Co.,Ltd. Address before: Room 102-86, building 6, 57 Andemen street, Yuhuatai District, Nanjing, Jiangsu 210000 Patentee before: Nanjing Tuodao Medical Technology Co.,Ltd. |
|
CP03 | Change of name, title or address |