CN112184782B - 一种骨关节的自动配准定位方法及装置 - Google Patents
一种骨关节的自动配准定位方法及装置 Download PDFInfo
- Publication number
- CN112184782B CN112184782B CN202010994656.5A CN202010994656A CN112184782B CN 112184782 B CN112184782 B CN 112184782B CN 202010994656 A CN202010994656 A CN 202010994656A CN 112184782 B CN112184782 B CN 112184782B
- Authority
- CN
- China
- Prior art keywords
- biplane
- synchronous dynamic
- ray
- image
- ray image
- 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
Links
- 210000000988 bone and bone Anatomy 0.000 title claims abstract description 107
- 238000000034 method Methods 0.000 title claims abstract description 31
- 238000002591 computed tomography Methods 0.000 claims abstract description 98
- 230000001360 synchronised effect Effects 0.000 claims abstract description 78
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 29
- 239000011159 matrix material Substances 0.000 claims abstract description 26
- 238000012937 correction Methods 0.000 claims abstract description 21
- 238000005457 optimization Methods 0.000 claims abstract description 21
- 230000000007 visual effect Effects 0.000 claims abstract description 20
- 238000012545 processing Methods 0.000 claims abstract description 6
- 239000002245 particle Substances 0.000 claims description 31
- 238000004364 calculation method Methods 0.000 claims description 14
- 230000005855 radiation Effects 0.000 claims description 13
- 230000006870 function Effects 0.000 claims description 6
- 238000012360 testing method Methods 0.000 claims description 6
- 238000009826 distribution Methods 0.000 claims description 4
- 238000005516 engineering process Methods 0.000 claims description 4
- 230000010354 integration Effects 0.000 claims description 4
- 230000003247 decreasing effect Effects 0.000 claims description 3
- 238000003384 imaging method Methods 0.000 claims description 3
- 238000013519 translation Methods 0.000 claims description 3
- 230000005484 gravity Effects 0.000 claims description 2
- 230000002285 radioactive effect Effects 0.000 claims description 2
- 238000011156 evaluation Methods 0.000 abstract description 2
- 210000000689 upper leg Anatomy 0.000 description 7
- 210000000629 knee joint Anatomy 0.000 description 4
- HPTJABJPZMULFH-UHFFFAOYSA-N 12-[(Cyclohexylcarbamoyl)amino]dodecanoic acid Chemical compound OC(=O)CCCCCCCCCCCNC(=O)NC1CCCCC1 HPTJABJPZMULFH-UHFFFAOYSA-N 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000005266 casting Methods 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 210000004872 soft tissue Anatomy 0.000 description 2
- 101100194606 Mus musculus Rfxank gene Proteins 0.000 description 1
- 230000002159 abnormal effect Effects 0.000 description 1
- 235000013405 beer Nutrition 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 201000008482 osteoarthritis Diseases 0.000 description 1
- 230000000149 penetrating effect Effects 0.000 description 1
- 230000008092 positive effect Effects 0.000 description 1
- 230000002980 postoperative effect Effects 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 230000002035 prolonged effect Effects 0.000 description 1
- 238000010845 search algorithm Methods 0.000 description 1
- 238000002922 simulated annealing Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 210000001519 tissue Anatomy 0.000 description 1
- 239000013598 vector Substances 0.000 description 1
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
- 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
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/004—Artificial life, i.e. computing arrangements simulating life
- G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/11—Region-based segmentation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/187—Segmentation; Edge detection involving region growing; involving region merging; involving connected component labelling
-
- 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
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Health & Medical Sciences (AREA)
- Computing Systems (AREA)
- Computational Linguistics (AREA)
- Data Mining & Analysis (AREA)
- Evolutionary Computation (AREA)
- Biomedical Technology (AREA)
- Molecular Biology (AREA)
- Biophysics (AREA)
- General Engineering & Computer Science (AREA)
- Artificial Intelligence (AREA)
- Mathematical Physics (AREA)
- Software Systems (AREA)
- Life Sciences & Earth Sciences (AREA)
- Health & Medical Sciences (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
本发明公开一种骨关节的自动配准定位方法，S1获取受试者的双平面同步动态X光影像以及CT体数据，并对CT体数据进行处理得到CT骨骼数据；S2对双平面同步动态X光影像进行空间校正，得到双平面同步动态X光影像的不同拍摄视角之间的相对位置；S3对校正后的X光影像进行去噪处理；S4基于CT骨骼数据生成虚拟X光影像，即DRR影像；S5计算DRR影像与对应视角的双平面同步动态X光影像相似度，采用第一优化算法更新不同时刻不同目标骨的位置变化矩阵，得到相对应的最终位置变化矩阵，从而确定骨关节的相对运动。本发明加快了DRR影像的生成速度，采用全变分去噪模型对X光影像降噪，大大提高了定位效率，通过自动配准的算法实现了对骨关节运动的快速评估。
Description
技术领域
本发明属于影像技术处理领域，尤其涉及一种骨关节的自动配准定位方法及装置。
背景技术
关节损伤或术后疗效不佳会导致的骨关节异常接触模式，将可能造成继发性的骨性关节炎，因此确定骨关节的精准位置将有助于为诊断提供有效信息。目前，现有技术中，利用两台交叉摆置的X光系成像系统搭配计算机断层扫描(CT)系统以对目标骨关节以进行二维影像到三维体数据配准进而重现骨关节的相对位置，是一种骨关节运动影像评估装置。结合X光与CT体数据重现骨关节运动的配准过程依赖数字投影重建技术(DigitallyReconstructed Radiograph，DRR)。通过将不同组织的CT值转化为X射线线性衰减系数，可以由Beer定律模拟X射线穿透CT体数据到达探测器平板后的信号强度得到仿真X射线影像。光线投射法(Ray Casting)由于其精度相对较高，在X射线仿真中被广泛应用。利用优化算法迭代体数据空间位置，得到不同虚拟X光影像，计算X光影像与真实X射线的相似度，得到最佳相似度下的骨关节空间位置。
然而，现有技术中利用X光影像来还原骨关节空间位置的优化搜索算法，该种方法需要生成大量DRR影像，然后DRR影像的生成需要对CT影像进行大量插值运算，其时间复杂度为O(n2)，大量生成DRR影像大大延长了骨关节空间定位的时间。此外，真实X光影像与DRR影像所属模态不同，其特征具有一定差异性，且真实X光影像噪点多，往往难以实现准确的配准。
发明内容
本发明的技术目的是提供一种骨关节的自动配准定位方法及装置，以获取骨关节的精准位置信息的技术效果。
为解决上述问题，本发明的技术方案为一种骨关节的自动配准定位方法，包括以下步骤：
S1：获取受试者的双平面同步动态X光影像以及CT体数据，并对CT体数据进行分割得到目标骨的CT骨骼数据；
S2：对双平面同步动态X光影像进行空间校正，得到双平面同步动态X光影像的不同拍摄视角之间的相对位置；
S3：对校正后的双平面同步动态X光影像进行去噪处理；
S4：对CT骨骼数据依据用于拍摄双平面同步动态X光影像的设备的成像参数以生成虚拟X光影像，即DRR影像；
S5：计算DRR影像与对应视角的双平面同步动态X光影像相似度，采用第一优化算法更新不同时刻不同目标骨的位置变化矩阵，得到相对应的最终位置变化矩阵；
S6：根据不同时刻不同目标骨的最终位置变化矩阵，确定骨关节的相对运动。
其中，在步骤S1中，对CT体数据进行分割得到目标骨的CT骨骼数据进一步包括以下步骤：
A1：在CT体数据内选取目标骨的种子点；
A2：根据种子点通过区域增长的方式进行分割，得到CT骨骼数据。
其中，步骤S2进一步包括以下步骤：
通过铅点网格对双平面同步动态X光影像进行空间校正，确定双平面同步动态X光影像的不同拍摄视角之间的相对位置，其中，不同拍摄视角之间的相对位置为不同组的X光发射器的放射源与X光接收器的相对位置，通过不同组的X光发射器与X光接收器拍摄得到双平面同步动态X光影像。
其中，步骤S3进一步包括：
建立全变分去噪模型，用于对空间校正后不同时刻的双平面同步动态X光影像进行去噪，利用第二优化算法求解全变分去噪模型得到能量最小化的泛函，以获取去噪X光影像。
具体地，全变分去噪模型的公式如下：
其中，v为去噪X光影像，u为校正后不同时刻的双平面同步动态X光影像，E(u,v)为保真项，用以避免去噪后的X光影像与去噪前的X光影像差异过多，γU(v)为全变分正则项，能容忍区域不连续性；
其中，E(u,v)为
即去噪后与去噪前对应位置像素差的平方和，U(v)为∑i,j|vi+1,j-vi,j|+|vi,j+1-vi,j|，即全变分模型。
进一步优选地，步骤S3之后还包括以下步骤：
以CT骨骼数据建立解剖坐标系，解剖坐标系的设置以目标骨的解剖形态决定；
基于解剖坐标系设置CT骨骼数据的初始位置，初始位置设置于直角坐标系的原点的最近点处，直角坐标系为以不同组的X光发射器的中心与X光接收器的中心连线的两线交点为原点的坐标系。
具体地，基于初始位置得到初始变化矩阵为：
T0＝(x,y,z,Rx,Ry,Rz)
其中，x、y、z分别为沿解剖坐标系的X轴、Y轴、Z轴的平移量，Rx,Ry,Rz分别为绕解剖坐标系X轴、Y轴、Z轴的旋转角度。
其中，在步骤S5中的第一优化算法为粒子群优化算法，其包括以下步骤：
S51：设定粒子群优化算法的初始参数，包括初始种群数量设定为N，最大迭代次数设定为G，惯性权重设定为λ，学习因子设定为μ1、μ2，搜索空间维度为对应待配准目标骨的运动空间维度Ф，初始速度设定为V，初始搜索空间位置为初始位置，对于初始种群采用相似度算法评价DRR影像与双平面同步动态X光影像相似度得到最佳相似度Fbest_global与初始最佳位置Pbest_global；
S52：对于N个粒子的种群中，第g次迭代种群后的第i个粒子的位置为
此时其速度为/>
通过相似性测试函数计算得到相似度为fi g；
对于第g次迭代种群的第i个粒子所历经的最佳相似度为Fbest_local，最佳搜索空间位置为Pbest_local；
对于第g次迭代种群的第i个粒子的速度更新为
其中r1，r2为随机值，搜索空间位置更新为/>
对于第g次迭代种群的权重设定为λg遵循线性递减权值策略进行更新,更新定义为
S53：当第g次迭代种群的第i个粒子相似度fi g>Fbest_global，将Fbest_global更新为fi g，Pbest_global更新为
否则Fbest_global，Pbest_global不更新；
S54：当迭代次数达到最大种群代数G时，所获得的最终最佳位置所对应的变化矩阵为最终位置变化矩阵Tfinal。
具体地，在步骤S52中，相似性测试函数进一步包括以下步骤：
B1：在直角坐标系下，基于轴对齐包围盒的线盒相交技术判断去噪后的双平面同步动态X光影像与CT体数据是否相交，并得到去噪后的双平面同步动态X光影像穿过CT体数据路径的大小参数；
B2：对步骤B1中的路径的大小参数进行插值、积分，得到去噪后的双平面同步动态X光影像穿过CT体数据的衰减系数的积分值，从而得到CT骨骼数据转化为DRR影像的各点像素值；
B3：计算去噪后的双平面同步动态X光影像与对应视角下DRR影像的归一化互信息，其计算定义为
其中H(DRR)为计算DRR影像的熵，H(F)为计算去噪后的双平面同步动态X光影像的熵，H(DRR,F)为计算影像的联合熵，影像的熵的计算定义为
其中hi为图像中灰度值为i的像素点总数，N表示图像的灰度级数；
DRR影像与去噪后的双平面同步动态X光影像的联合熵计算定义为
其中，/>
为DRR影像与去噪后的双平面同步动态X光影像的联合概率密度分布；
B4：对于不同视角的去噪后的双平面同步动态X光影像与DRR影像计算其总体相似度F＝wF1×NMIF1+wF2×NMIF2，令fi g＝F，其中，wF1为第一视角所占比重，wF2为第二视角所占比重，NMI为对应视角下的去噪后的双平面同步动态X光影像与DRR影像的归一化互信息。
一种骨关节的自动配准定位装置包括：
数据采集模块，用于获取受试者的双平面同步动态X光影像以及CT体数据；
空间校正模块，用于对双平面同步动态X光影像进行空间校正，得到拍摄双平面同步动态X光影像的放射源与探测器之间的相对位置；
影像处理模块，用于对空间校正后的不同时刻的双平面同步动态X光影像去噪，基于计算所得的目标骨的位置生成虚拟X光影像，对双平面同步动态X光影像与虚拟X光影像相似度作比较；
初始定位模块，用于基于放射源与探测器之间的相对位置和CT体数据中骨关节的解剖坐标系得到初始位置变化矩阵；
精确定位模块，用于基于若干组X光影像和的虚拟X光影像，利用粒子群算法进行空间位置搜索，比较同一视角的X光影像与虚拟X光影像的相似度，迭代更新最佳空间位置直到得到所对应的最终位置变化矩阵；
运动定位模块，用于根据最终位置变化矩阵，确定骨关节的空间位置运动结果。
本发明由于采用以上技术方案，使其与现有技术相比具有以下的优点和积极效果：
本发明通过使用计算机图形卡并进行粒子群优化算法搜索最优空间位置，加快了DRR影像的生成速度，从而减少了骨关节空间定位所需的时间，大大提高了定位效率。利用全变分去噪的方法提取X光影像的轮廓，全变分去噪模型能在有效去噪的同时保留足够的图像边缘信息，克服了X光影像噪点多、提取轮廓难的问题。因此，本发明能够鲁棒地实现目标骨关节的空间定位，将有助于利用X光影像与CT体数据对患者进行精准的骨关节定位，方便后续治疗。
附图说明
通过阅读下文优选实施方式的详细描述，各种其他的优点和益处对于本领域普通技术人员将变得清楚明了。附图仅用于示出优选实施方式的目的，而并不认为是对本发明的限制。
图1为本发明一个实施例中的骨关节的自动配准定位方法的流程示意图；
图2为本发明一个实施例中的骨关节的自动配准定位方法的粒子群优化算法的流程示意图；
图3为本发明一个实施例中的骨关节的自动配准定位方法的相似性测试算法的流程示意图；
图4为本发明一个实施例中的骨关节的自动配准定位装置结构框图；
图5为本发明一个实施例中的骨关节的自动配准定位方法的X光机采集数据示意图；
图6为本发明一个实施例中的骨关节的自动配准定位方法的CT骨骼数据位置的模拟X光射线示意图。
具体实施方式
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对照附图说明本发明的具体实施方式。显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图，并获得其他的实施方式。
为使图面简洁，各图中只示意性地表示出了与本发明相关的部分，它们并不代表其作为产品的实际结构。另外，以使图面简洁便于理解，在有些图中具有相同结构或功能的部件，仅示意性地绘示了其中的一个，或仅标出了其中的一个。在本文中，“一个”不仅表示“仅此一个”，也可以表示“多于一个”的情形。
以下结合附图和具体实施例对本发明提出的一种骨关节的自动配准定位方法及装置作进一步详细说明。根据下面说明和权利要求书，本发明的优点和特征将更清楚。
实施例1
参看图1和图5，本实施例提供一种骨关节的自动配准定位方法，包括以下步骤，S1：获取受试者的双平面同步动态X光影像以及CT体数据，并对CT体数据进行分割得到目标骨的CT骨骼数据；S2：对双平面同步动态X光影像进行扭曲校正，得到双平面同步动态X光影像的不同拍摄视角之间的相对位置；S3：对不同时刻的双平面同步动态X光影像分别进行提取轮廓并进行膨胀操作，得到多组X光影像加权轮廓信息；S4：在CT骨骼数据的解剖坐标系下，基于不同拍摄视角之间的相对位置，得到不同目标骨的初始位置变化矩阵；S5：在CT体数据拍摄的直角坐标系下，基于CT骨骼数据和多组X光影像加权轮廓信息，分别利用模拟退火算法对初始位置变化矩阵进行空间搜索，得到不同时刻对应的最终位置变化矩阵；S6：根据不同时刻对应的最终位置变化矩阵，通过目标骨之间的相对运动确定骨关节的空间位置变化。
现对本实施例进行详细说明：
参看图1，在本实施例中，步骤S1具体为：以膝关节的自动配准定位为例，首先需要获取受试者膝关节的X光影像和CT体数据，利用两台Philps C形臂X光机拍摄受试者的膝关节部位，得到双平面同步动态X光影像，其中，双平面是指两台X光机拍摄，得到两个不同拍摄视角的X光影像；同步是指将两台X光机拍摄时间相同的X光影像组合起来，由于两台X光机启动拍摄之间存在时间误差，因此需要通过程序进行配对，使每组两张X光影像是各自X光机启动后相同时间段拍摄的；动态是指两台X光机连续拍摄多张X光影像。利用CT机拍摄受试者的膝关节部位，得到CT体数据。后续，将双平面同步动态X光影像和CT体数据以Dicom影像格式进行图像处理。
在本实施例中，还需要对上述得到的CT体数据进行预处理并建立解剖坐标系，此为准备事项，可在步骤S4及之前任意时间段操作，具体地，首先在CT体数据内手动选取目标骨的种子点，即选取股骨部分的种子点；根据该种子点进行区域增长的方式进行分割，提取出目标骨的CT骨骼数据，以剥离股骨及其周围的软组织，以减少软组织对匹配定位的影响；以股骨内、外侧髁的最突起部分连线作为Z轴，以连线的中心作为解剖坐标系的中心，拟合股骨干确定ZOY平面，垂直于平面股骨向前方作X轴，Y轴垂直于XOZ平面指向身体近端。上述只是说明了股骨部分的解剖坐标系的构建，针对不同的目标骨其所需构建的解剖坐标系各不相同。
步骤S2具体为：对任一时刻的一组双平面同步动态X光影像进行空间校正，由于拍摄出的X光影像与实际的两台X光机各自的放射源与探测器之间的X射线方向存在偏差，因此需要建立虚拟的双平面X光拍摄空间，将两台X光机的放射源与探测器的位置放入虚拟的双平面X光拍摄空间，该拍摄空间内设有标定位置的铅点网格，以确定两台X光机的放射源与探测器的空间位置从而确定多项式校正系数，通过多项式校正系数从而得到矫正后的双平面同步动态X光影像。
在本实施例中，还需要得到CT体数据的直角坐标系，即包括股骨的一整体CT体数据，基于上述在拍摄空间中得到两台X光机各自的放射源与探测器的相对位置，将两台X光机各自的放射源与探测器的中心连线，两条连线的交点位置作为原点建立直角坐标系，且直角坐标系的X、Y、Z轴与股骨的CT体数据的各边平行。
参看图1，在本实施例中，步骤S3具体为：建立全变分去噪模型，用于对空间校正后各组双平面同步动态X光影像进行去噪，利用第二优化算法即原始对偶算法求解全变分去噪模型得到能量最小化的泛函，以获取去噪X光影像。
TVL1去噪模型公式如下：
其中，v为去噪后X光影像，u为含噪声的原始X光影像F1或F2，E(u,v)定义为
即去噪后与去噪前对应位置像素差的平方和，U(v)定义为∑i,j|vi+1,j-vi,j|+|vi,j+1-vi,j|，即全变分模型；E(u,v)；为保真项以避免去噪后的图像与原始图像差异过多，γU(v)为全变分正则项，能够容忍一定区域不连续性，此去噪模型能够保证对于X光在有效去噪的同时保留足够的图像边缘信息。
参看图1，在本实施例中，在步骤S3与步骤S4之间还包括以下步骤：基于解剖坐标系设置CT骨骼数据的一初始位置，初始位置设置于与直角坐标系的原点的最近点处。其中，基于初始位置得到初始变化矩阵为：
T0＝(x,y,z,Rx,Ry,Rz)
其中，x、y、z分别为沿解剖坐标系的X轴、Y轴、Z轴的平移量，Rx,Ry,Rz分别为绕解剖坐标系X轴、Y轴、Z轴的旋转角度。
参看图1，在本实施例中，步骤S4具体为：需要得到双平面同步动态X光影像的DRR影像，将双平面同步动态X光影像导入GPU的常量内存中，利用初始变化矩阵T0将CT骨骼数据从初始位置转换到对应位置，将目标骨的CT骨骼数据传输到支持统一计算设备架构(Compute Unified Device Architecture，CUDA)计算机图形卡的纹理寄存器中，利用计算机图形卡的纹理寄存器提供的告诉缓存进行三次线性插值以获得X射线经过CT体数据时的数值，通过对每束X射线的路径上经过体数据进行积分，如图6所示，计算获得对应CT位置时的模拟X光射线，从而得到虚拟X光影像即DRR影像。
参看图1和图3，步骤S5具体为：利用粒子群优化算法对初始变化矩阵V0进行空间搜索，得到最终位置变化矩阵Tfinal。其中，
设置粒子群优化算法所需的初始参数，具体地，初始种群数量设定为N，最大迭代次数设定为G，在本实施例中种群数量设定为30，最大迭代次数设定为100；惯性权重设定为λ，学习因子设定为μ1,μ2，在本实施例中初始惯性权重λ1＝0.9，λ2＝0.4，学习因子μ1＝2，μ2＝2，搜索空间维度为对应待配准目标骨的运动空间维度Ф，在本实施例中设置为Ф＝6，计算初始种群适应度，即计算初始种群中各粒子的图像相似度，找出初始种群中的最佳初始相似度Fbest_global及其对应的最佳初始位置Pbest_global，第g次迭代种群后的第i个粒子的位置表示为
速度表示为/>
其中，当前粒子的位置/>
速度位置/>
中的下标均代表所选择的粒子i，上标则均表示第g次迭代后的种群，其中/>
也对应了位置变化矩阵/>
对于第g次迭代后的种群，其惯性权重也随迭代次数而衰减，且遵循线性递减权值策略，即
在此实例中/>
在第g次迭代种群中的每个粒子，其速度与位置也相应更新，其速度更新为
其中r1，r2为随机值，位置更新为/>
其所历经的最佳相似度为Fbest _local，最佳位置为Pbest_local。在当前第g次迭代后的种群中，对每个粒子判断其适应度，如对第i个粒子，如果其相似度fi g>Fbest_global，将Fbest_global更新为fi g，Pbest_global更新为/>
否则Fbest_global，Pbest_global不更新。如在本实例中第50次迭代时得到数量为30的粒子种群，对于每个粒子判断其适应度与Fbest_global大小，完成比较后更新最大Fbest_local且高于Fbest_global为当前的Fbest_global，且将其对应的Pbest_local为更新为当前的最佳位置Pbest_global，完成后进行下一次的种群迭代。完成100次迭代后所得目标骨空间位置Pbest_global即为最终所得到的Tfinal。
参看图1、图4和图6，在本实施例中，现对相似性测试函数S(Vi)计算进行详细说明：
在直角坐标下，基于轴对齐包围盒的线盒相交技术判断X射线与CT体数据是否相交，并得到X射线穿过CT数据的路径；利用计算机图形卡纹理寄存器对X射线穿过目标骨的CT骨骼数据的路径进行插值、积分，得到X射线穿过CT骨骼数据的衰减系数的积分值，从而得到CT骨骼数据转化为DRR影像后，DRR影像的各点像素值。
具体地，基于CT骨骼数据、两台X光机的放射源位置S1、S2与探测器位置D1、D2进行计算，分别计算放射源S1、S2到D1、D2每个像素位置X射线的单位向量dir，通过轴对齐包围盒(Axis aligned Bounding box)线盒相交技术判断是否X射线与CT骨骼数据相交，并返回穿过CT骨骼数据的路径，计算公式如下：ttop＝(CTtop-S1)/dir；tbot＝(CTbottom-S1)/dir，其中ttop、tbot为X射线与CT骨骼数据各面所在平面相交的数值，CTtop为CT骨骼数据的顶点坐标值，其值均为正，CTbottom为CT骨骼数据的底部坐标值；接着，通过公式tmin＝min(ttop,tbot)，tmax＝max(ttop,tbot)取得最值，其实际物理意义是指沿着X射线方向较近、较远的点；然后，通过公式tnear＝max(max(tmin.x,tmin.y)，max(tmin.x,tmin.z))、tfar＝min(min(tmax.x,tmax.y),min(tmax.x,tmax.z))，得到X射线进入CT骨骼数据的点tnear、X射线穿出CT骨骼数据的点tfar，当tfar大于tnear，X射线穿过CT骨骼数据，否则X射线不穿过CT骨骼数据，tfar-tnear则代表X射线经CT骨骼数据的路径。利用计算机图形卡纹理寄存器对X射线经过的CT骨骼数据的路径进行插值，并进行积分，得到X射线穿过CT的衰减系数的积分值，设定X射线步长step为0.5，CT骨骼数据转化为DRR图像的每个点的像素值计算步骤如下：若travelLength<(tfar-tnear)，travelLength为X射线在CT骨骼数据内已插值计算的路径，则对该路径上的未插值计算的一点进行如下操作，t＝tex3D(CT骨骼数据未插值计算的一点，即当前位置)，对该点进行三次线性插值，t为插值后该点的值，其中，tex3D为支持CUDA技术的计算机图形卡自带三次线性插值函数；利用插值的数据在积分路径上做积分，DRR(i)＝DRR(i-1)+step*t，得到DRR图像该点的像素值；然后更新该路径上的下一点位置，下一位置＝当前位置+dir/CT每点体素尺寸，travelLength＝当前位置的路程+step，后续进行迭代后的travelLength与(tfar-tnear)之间的判断，直至计算完该路径上所有点的像素；通过多条X射线计算的结果进行叠加，得到该DRR影像整体的各点像素值。
计算去噪后的双平面同步动态X光影像与对应视角下DRR影像的归一化互信息，其计算定义为
其中H(DRR)为计算DRR影像的熵，H(F)为计算去噪后的双平面同步动态X光影像的熵，H(DRR,F)为计算影像的联合熵，影像的熵的计算定义为
其中hi为图像中灰度值为i的像素点总数，N表示图像的灰度级数；
DRR影像与去噪后的双平面同步动态X光影像的联合熵计算定义为
其中，/>
为DRR影像与去噪后的双平面同步动态X光影像的联合概率密度分布。
对DRR影像与对应的实际的X光影像计算相似度，即其互信息，其中，相似度计算公式为S(Vi)＝wF1×NMIF1+wF2×NMIF2，在此实例中，wF1＝wF2＝0.5。
最后步骤S6：根据不同时刻对应的最终位置变化矩阵Vfinal，从而确定目标骨之间的关节空间位置。
实施例2
一种骨关节的自动配准定位装置包括：
数据采集模块，参看图4和图5，用于获取受试者的双平面同步动态X光影像以及CT体数据，在本实施例中数据采集模块包括2台用于拍摄双平面同步动态X光影像的Philps C形臂X光机，当然也可以为其它种类的X光机，以及同于拍摄CT数据的CT机。
空间校正模块，用于对双平面同步动态X光影像进行空间定位，得到拍摄双平面同步动态X光影像的放射源与探测器之间的相对位置，即为得到DRR图像的各点像素做好准备。
影像处理模块，用于对不同时刻的双平面同步动态X光影像去噪，基于计算所得的目标追踪骨骼的位置生成虚拟X光影像，对X光影像与虚拟X光影像相似度作比较，在本实施例中对双平面同步动态X光影像进行上述若干处理，最终得到X光影像与预估位置所生成虚拟X光影像的归一化互信息。利用全变分去噪的方法提取X光影像的电子噪声，降低了DRR影像与真实X光影像的差异型，提高了DRR与X光影像相似度计算的准确性。
初始定位模块，用于基于X光机放射源与探测器之间的相对位置和CT体数据中骨关节的解剖坐标系得到初始位置变化矩阵，基于初始位置变化矩阵通过精确定位模块的粒子群优化算法不断迭代以及调整，得到最终位置变化矩阵，为运动定位模块提供计算基础。
精确定位模块，用于基于若干组X光去噪影像和利用预估目标骨的空间位置所生成的DRR影像，同时利用粒子群优化算法对初始位置变化矩阵进行空间搜索，得到若干组X光影像和DRR影像的归一化互信息比较以及迭代完成后所对应的最终位置变化矩阵；由于在计算机图像软件中运用了纹理内存对目标骨的CT数据进行插值处理，加快了DRR影像的生成速度，从而减少了骨关节空间定位所需的时间，大大提高了定位效率。
运动定位模块，用于根据最终位置变化矩阵，确定骨关节的空间位置运动结果，此外，本发明能够鲁棒地实现目标骨关节的空间定位，将有助于利用X光影像与CT体数据对患者进行精准的骨关节定位，方便后续治疗。
上面结合附图对本发明的实施方式作了详细说明，但是本发明并不限于上述实施方式。即使对本发明作出各种变化，倘若这些变化属于本发明权利要求及其等同技术的范围之内，则仍落入在本发明的保护范围之中。
Claims (9)
1.一种骨关节的自动配准定位方法，其特征在于，包括以下步骤：
S1：获取受试者的双平面同步动态X光影像以及CT体数据，并对所述CT体数据进行分割得到目标骨的CT骨骼数据；
S2：对所述双平面同步动态X光影像进行空间校正，得到所述双平面同步动态X光影像的不同拍摄视角之间的相对位置；
S3：对校正后的所述双平面同步动态X光影像进行去噪处理；
S4：对所述CT骨骼数据依据用于拍摄所述双平面同步动态X光影像的设备的成像参数以生成虚拟X光影像，即DRR影像；
S5：计算若干组所述DRR影像与对应视角的所述双平面同步动态X光影像相似度，采用第一优化算法进行空间位置搜索，比较同一视角的所述双平面同步动态X光影像与所述DRR影像的相似度，迭代更新不同时刻不同目标骨的位置变化矩阵，直到得到相对应的最终位置变化矩阵，所述第一优化算法为粒子群优化算法；
S6：根据不同时刻不同目标骨的所述最终位置变化矩阵，确定骨关节的相对运动；
其中，在所述步骤S5中，通过相似性测试函数获取相似度，包括以下步骤：
B1：在直角坐标系下，基于轴对齐包围盒的线盒相交技术判断去噪后的所述双平面同步动态X光影像与所述CT体数据是否相交，并得到去噪后的所述双平面同步动态X光影像穿过所述CT体数据路径的大小参数；
B2：对所述步骤B1中的路径的大小参数进行插值、积分，得到去噪后的所述双平面同步动态X光影像穿过所述CT体数据的衰减系数的积分值，从而得到所述CT骨骼数据转化为DRR影像的各点像素值；
B3：计算去噪后的所述双平面同步动态X光影像与对应视角下所述DRR影像的归一化互信息，其计算定义为
其中H(DRR)为计算所述DRR影像的熵，H(F)为计算去噪后的所述双平面同步动态X光影像的熵，H(DRR，F)为计算影像的联合熵，影像的熵的计算定义为/>
其中hi为图像中灰度值为i的像素点总数，N表示图像的灰度级数；
所述DRR影像与去噪后的所述双平面同步动态X光影像的联合熵计算定义为
其中，/>
为所述DRR影像与所述去噪后的所述双平面同步动态X光影像的联合概率密度分布；
B4：对于不同视角的去噪后的所述双平面同步动态X光影像与所述DRR影像计算其总体相似度F＝wF1×NMIF1+wF2×NMIF2，令fi g＝F，其中，fi g为相似度，wF1为第一视角所占比重，wF2为第二视角所占比重，NMI为对应视角下的去噪后的所述双平面同步动态X光影像与所述DRR影像的归一化互信息。
2.根据权利要求1所述的骨关节的自动配准定位方法，其特征在于，在所述步骤S1中，所述对所述CT体数据进行分割得到目标骨的CT骨骼数据进一步包括以下步骤：
A1：在所述CT体数据内选取目标骨的种子点；
A2：根据所述种子点通过区域增长的方式进行分割，得到所述CT骨骼数据。
3.根据权利要求1所述的骨关节的自动配准定位方法，其特征在于，所述步骤S2进一步包括以下步骤：
通过铅点网格对所述双平面同步动态X光影像进行空间校正，确定所述双平面同步动态X光影像的不同拍摄视角之间的相对位置，其中，所述不同拍摄视角之间的相对位置为不同组的X光发射器的放射源与X光接收器的相对位置，通过不同组的所述X光发射器与所述X光接收器拍摄得到所述双平面同步动态X光影像。
4.根据权利要求1所述的骨关节的自动配准定位方法，其特征在于，所述步骤S3进一步包括：
建立全变分去噪模型，用于对空间校正后不同时刻的所述双平面同步动态X光影像进行去噪，利用第二优化算法求解所述全变分去噪模型得到能量最小化的泛函，以获取去噪X光影像。
5.根据权利要求4所述的骨关节的自动配准定位方法，其特征在于，所述全变分去噪模型的公式如下：
其中，v为所述去噪X光影像，u为校正后不同时刻的所述双平面同步动态X光影像，E(u，v)为保真项，用以避免去噪后的X光影像与去噪前的X光影像差异过多，γU(v)为全变分正则项，能容忍区域不连续性；
其中，E(u，v)为
即去噪后与去噪前对应位置像素差的平方和，U(v)为∑i，j|vi+1，j-vi，j|+|vi，j+1-vi，j|，即全变分模型。
6.根据权利要求1所述的骨关节的自动配准定位方法，其特征在于，所述步骤S3之后还包括以下步骤：
以所述CT骨骼数据建立解剖坐标系，所述解剖坐标系的设置以目标骨的解剖形态决定；
基于所述解剖坐标系设置CT骨骼数据的初始位置，所述初始位置设置于直角坐标系的原点的最近点处，所述直角坐标系为以不同组的X光发射器的中心与X光接收器的中心连线的两线交点为原点的坐标系。
7.根据权利要求6所述的骨关节的自动配准定位方法，其特征在于，基于所述初始位置得到初始变化矩阵为：
T0＝(x，y，z，Rx，Ry，Rz)
其中，x、y、z分别为沿所述解剖坐标系的X轴、Y轴、Z轴的平移量，Rx，Ry，Rz分别为绕所述解剖坐标系X轴、Y轴、Z轴的旋转角度。
8.根据权利要求7所述的骨关节的自动配准定位方法，其特征在于，在所述步骤S5中的第一优化算法为粒子群优化算法，其包括以下步骤：
S51：设定所述粒子群优化算法的初始参数，包括初始种群数量设定为N，最大迭代次数设定为G，惯性权重设定为λ，学习因子设定为μ1、μ2，搜索空间维度为对应待配准目标骨的运动空间维度Ф，初始速度设定为V，初始搜索空间位置为所述初始位置，对于所述初始种群采用相似度算法评价DRR影像与所述双平面同步动态X光影像相似度得到最佳相似度Fbest _globa与初始最佳位置Pbest_global；
S52：对于N个粒子的种群中，第g次迭代种群后的第i个粒子的位置为
此时其速度为/>
通过相似性测试函数计算得到相似度为fi g；
对于第g次迭代种群的第i个粒子所历经的最佳相似度为Fbest_local，最佳搜索空间位置为Pbest_local；
对于第g次迭代种群的第i个粒子的速度更新为
其中r1，r2为随机值，搜索空间位置更新为/>
对于第g次迭代种群的权重设定为λg遵循线性递减权值策略进行更新,更新定义为
S53：当第g次迭代种群的第i个粒子相似度fi g>Fbest_global，将Fbest_global更新为fi g，pbest _global更新为
否则Fbest_global，Pbest_global不更新；
S54：当迭代次数达到最大种群代数G时，所获得的最终最佳位置所对应的变化矩阵为最终位置变化矩阵Tfinal。
9.一种骨关节的自动配准定位装置，配置有如权利要求1至8任意一项所述的骨关节的自动配准定位方法，其特征在于，包括：
数据采集模块，用于获取受试者的双平面同步动态X光影像以及CT体数据；
空间校正模块，用于对所述双平面同步动态X光影像进行空间校正，得到拍摄所述双平面同步动态X光影像的放射源与探测器之间的相对位置；
影像处理模块，用于对空间校正后的不同时刻的所述双平面同步动态X光影像去噪，基于计算所得的目标骨的位置生成虚拟X光影像，对所述双平面同步动态X光影像与所述虚拟X光影像相似度作比较；
初始定位模块，用于基于所述放射源与所述探测器之间的相对位置和所述CT体数据中骨关节的解剖坐标系得到初始位置变化矩阵；
精确定位模块，用于基于若干组所述X光影像和所述的虚拟X光影像，利用粒子群算法进行空间位置搜索，比较同一视角的所述X光影像与所述虚拟X光影像的相似度，迭代更新最佳空间位置直到得到所对应的最终位置变化矩阵；
运动定位模块，用于根据所述最终位置变化矩阵，确定骨关节的空间位置运动结果。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010994656.5A CN112184782B (zh) | 2020-09-21 | 2020-09-21 | 一种骨关节的自动配准定位方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010994656.5A CN112184782B (zh) | 2020-09-21 | 2020-09-21 | 一种骨关节的自动配准定位方法及装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112184782A CN112184782A (zh) | 2021-01-05 |
| CN112184782B true CN112184782B (zh) | 2023-05-23 |
Family
ID=73955683
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010994656.5A Active CN112184782B (zh) | 2020-09-21 | 2020-09-21 | 一种骨关节的自动配准定位方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112184782B (zh) |
Families Citing this family (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113040908B (zh) * | 2021-02-02 | 2022-03-25 | 武汉联影智融医疗科技有限公司 | 手术导航的配准方法、装置、计算机设备和存储介质 |
| CN113421226B (zh) * | 2021-06-03 | 2022-11-01 | 山东师范大学 | 基于互信息的ct-dr多模态食管图像配准方法及系统 |
| CN113570648B (zh) * | 2021-07-30 | 2023-09-26 | 武汉联影智融医疗科技有限公司 | 多骨骼影像配准方法、电子装置以及医学导航系统 |
| CN113870331B (zh) * | 2021-10-07 | 2024-07-26 | 浙江大学 | 一种基于深度学习的胸部ct与x光实时配准算法 |
| CN114748083A (zh) * | 2022-03-18 | 2022-07-15 | 上海涛影医疗科技有限公司 | 一种多视角曝光x光影像定位方法和系统 |
| CN115068110A (zh) * | 2022-06-14 | 2022-09-20 | 中国人民解放军总医院第一医学中心 | 用于股骨颈骨折手术导航的图像配准方法及系统 |
| CN115082534B (zh) * | 2022-07-21 | 2022-12-16 | 杭州三坛医疗科技有限公司 | 双平面图像配准方法、装置及机器人 |
| CN117726752B (zh) * | 2023-12-14 | 2024-08-30 | 凝动万生医疗科技(武汉)有限公司 | 一种快速的三维头颅数据正位校准方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009043224A1 (fr) * | 2007-09-24 | 2009-04-09 | Shu Jia | Procédé de reconstruction d'image utilisant une photographie volumique par rayons x |
| CN102147919A (zh) * | 2010-02-10 | 2011-08-10 | 昆明医学院第一附属医院 | 一种校正术前三维图像的术中配准方法和装置 |
| CN104637061A (zh) * | 2015-01-30 | 2015-05-20 | 中国科学院自动化研究所 | 一种二维三维医学图像配准方法 |
-
2020
- 2020-09-21 CN CN202010994656.5A patent/CN112184782B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009043224A1 (fr) * | 2007-09-24 | 2009-04-09 | Shu Jia | Procédé de reconstruction d'image utilisant une photographie volumique par rayons x |
| CN102147919A (zh) * | 2010-02-10 | 2011-08-10 | 昆明医学院第一附属医院 | 一种校正术前三维图像的术中配准方法和装置 |
| CN104637061A (zh) * | 2015-01-30 | 2015-05-20 | 中国科学院自动化研究所 | 一种二维三维医学图像配准方法 |
Non-Patent Citations (4)
| Title |
|---|
| 基于PDE的医学图像增强技术研究;荆萍;《中国优秀博硕士学位论文全文数据库（硕士）信息科技辑》;20110815;第24-30页 * |
| 浅谈图像的全变分和去噪;TomHeaven;《https://blog.csdn.net/hanlin_tan/article/details/52448803》;20160906;第1-2页 * |
| 脊柱X光透视图像与CT体积图像的配准研究;陆丽英;《http://www.doc88.com/p-1844789005195.html》;20130307;第4-8,15-17,29-30,40-41 * |
| 陆丽英.脊柱X光透视图像与CT体积图像的配准研究.《http://www.doc88.com/p-1844789005195.html》.2013, * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112184782A (zh) | 2021-01-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112184782B (zh) | 一种骨关节的自动配准定位方法及装置 | |
| CN109949899B (zh) | 图像三维测量方法、电子设备、存储介质及程序产品 | |
| Ehlke et al. | Fast generation of virtual X-ray images for reconstruction of 3D anatomy | |
| US20170278302A1 (en) | Method and device for registering an image to a model | |
| WO2021238171A1 (zh) | 图像配准方法及其相关的模型训练方法、设备、装置 | |
| US10350434B2 (en) | Patient-specific radiation dose assessment in medical therapy | |
| CN112614169B (zh) | 基于深度学习网络的2d/3d脊椎ct层级配准方法 | |
| TWI517093B (zh) | Computer tomography reconstruction method | |
| CN112598649B (zh) | 基于生成对抗网络的2d/3d脊椎ct非刚性配准方法 | |
| US20110019791A1 (en) | Selection of optimal views for computed tomography reconstruction | |
| US20240221190A1 (en) | Methods and systems for registration | |
| Afzal et al. | Rgb-d multi-view system calibration for full 3d scene reconstruction | |
| KR102036834B1 (ko) | 이미지 처리 방법 | |
| CN110270015B (zh) | 一种基于多序列MRI的sCT生成方法 | |
| CN117159138B (zh) | 一种基于关节定位的2d-3d配准方法及系统 | |
| Bögel et al. | Respiratory Motion Compensation Using Diaphragm Tracking for Cone‐Beam C‐Arm CT: A Simulation and a Phantom Study | |
| Hatt et al. | Robust 5DOF transesophageal echo probe tracking at fluoroscopic frame rates | |
| Oulbacha et al. | MRI to C‐arm spine registration through Pseudo‐3D CycleGANs with differentiable histograms | |
| CN113729747B (zh) | 一种球形金属标记的锥束ct金属伪影去除系统及去除方法 | |
| CN112085833B (zh) | 一种锥束ct与影像融合结合的颈椎在体三维运动的分析方法 | |
| Villa-Uriol et al. | Automatic creation of three-dimensional avatars | |
| Ju et al. | Multi-view stereophotogrammetry for post-mastectomy breast reconstruction | |
| Docef et al. | Reconstruction of 4D deformed CT for moving anatomy | |
| Chen et al. | 3-D model reconstruction from C-arm images | |
| Bodensteiner et al. | Motion and positional error correction for cone beam 3D-reconstruction with mobile C-arms |
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 | ||
| TA01 | Transfer of patent application right |
Effective date of registration: 20211220 Address after: 201210 room 405-2, 4th floor, building 9, No. 1206, Zhangjiang Road, China (Shanghai) pilot Free Trade Zone, Pudong New Area, Shanghai Applicant after: SHANGHAI TAOYING MEDICAL TECHNOLOGY CO.,LTD. Address before: 200240 No. 800, Dongchuan Road, Shanghai, Minhang District Applicant before: SHANGHAI JIAO TONG University |
|
| TA01 | Transfer of patent application right | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |