CN114155286B - 一种骨骼ct图像的解剖形态和材料力学特性模板库的个性化配准方法 - Google Patents

一种骨骼ct图像的解剖形态和材料力学特性模板库的个性化配准方法 Download PDF

Info

Publication number
CN114155286B
CN114155286B CN202111317598.3A CN202111317598A CN114155286B CN 114155286 B CN114155286 B CN 114155286B CN 202111317598 A CN202111317598 A CN 202111317598A CN 114155286 B CN114155286 B CN 114155286B
Authority
CN
China
Prior art keywords
image
shape
bone
patient
model
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
CN202111317598.3A
Other languages
English (en)
Other versions
CN114155286A (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.)
Dalian University of Technology
Original Assignee
Dalian University of Technology
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 Dalian University of Technology filed Critical Dalian University of Technology
Priority to CN202111317598.3A priority Critical patent/CN114155286B/zh
Publication of CN114155286A publication Critical patent/CN114155286A/zh
Application granted granted Critical
Publication of CN114155286B publication Critical patent/CN114155286B/zh
Priority to US17/841,222 priority patent/US20230146649A1/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/0002Inspection of images, e.g. flaw detection
    • G06T7/0012Biomedical image inspection
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/30Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/33Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/30Determination of transform parameters for the alignment of images, i.e. image registration
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/50Information retrieval; Database structures therefor; File system structures therefor of still image data
    • G06F16/53Querying
    • G06F16/532Query formulation, e.g. graphical querying
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/30Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/38Registration of image sequences
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/70Determining position or orientation of objects or cameras
    • G06T7/73Determining position or orientation of objects or cameras using feature-based methods
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/90Determination of colour characteristics
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00Arrangements for image or video recognition or understanding
    • G06V10/70Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/74Image or video pattern matching; Proximity measures in feature spaces
    • G06V10/75Organisation of the matching processes, e.g. simultaneous or sequential comparisons of image or video features; Coarse-fine approaches, e.g. multi-scale approaches; using context analysis; Selection of dictionaries
    • G06V10/758Involving statistics of pixels or of feature values, e.g. histogram matching
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00Arrangements for image or video recognition or understanding
    • G06V10/70Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/77Processing 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/7715Feature extraction, e.g. by transforming the feature space, e.g. multi-dimensional scaling [MDS]; Mappings, e.g. subspace methods
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00Scenes; Scene-specific elements
    • G06V20/60Type of objects
    • G06V20/64Three-dimensional objects
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10081Computed x-ray tomography [CT]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30004Biomedical image processing
    • G06T2207/30008Bone
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V2201/00Indexing scheme relating to image or video recognition or understanding
    • G06V2201/03Recognition of patterns in medical or anatomical images
    • G06V2201/033Recognition of patterns in medical or anatomical images of skeletal patterns

Abstract

本发明属于数字医疗、图像处理技术领域,提出了一种骨骼CT图像的解剖形态和材料力学特性模板库的个性化配准方法。该方法以大量健康人的骨骼CT图像构建能包含骨骼解剖形态与材料力学特性的统计模型,通过统计模型与骨骼CT图像之间个性化的配准方法实现对患者骨骼的参数化描述,使用患者图像与配准参数构建假体模板库,模板库通过个性化配准方法匹配患者CT图像,检索模板库中与患者骨骼情况相近的模板图像与假体模型,作为个性化假体植入物设计的初始参考模板。通过本发明提出的模板库个性化配准方法,将为个性化骨骼植入物的设计、骨骼的材料力学特性分析、骨骼假体库的构建等实际临床应用场景提供更多的科学保障。

Description

一种骨骼CT图像的解剖形态和材料力学特性模板库的个性化 配准方法
技术领域
本发明涉及数字医疗、图像处理技术领域,具体涉及一种骨骼CT图像的解剖形态和材料力学特性模板库的个性化配准方法图像,主要适用于患者受损骨骼假体植入物的形状设计与材料力学特性分析。
背景技术
人体骨骼系统经常会因为受到病变、创伤、骨折等因素的影响而出现受损的情况,在多数情况下需要使用假体植入物进行外科手术来对受损骨骼进行修复。为了确保手术干预的成功,需要医生制定良好的术前计划、设计选择合适的个性化假体植入物。近年来随着个性化数字医疗与3D打印技术的发展,专业的骨科医生会采用患者的计算机断层扫描(Computed Tomography,CT)图像来构建三维骨骼模型,分析个体骨骼的解剖形态并为患者定制个性化的假体植入物,使用3D打印技术进行假体制作。假体植入物的个性化设计是一项步骤繁琐且耗时的工作,如果建立一个包含解剖形态和材料力学特性的假体模板库,实现模板库与病人骨骼CT图像的个性化配准,则可以检索匹配到模板库中与病人骨骼的形态与灰度相似的模板图像以及假体模型,为医生进行个性化的假体植入物设计提供一个初始参考模板,节省设计的时间。另外,对于病人骨骼形状与灰度的参数化表示形式,也可用于虚拟手术计划、骨骼的缺陷分析、个性化植入物的材料设计等方面的研究。
为实现基于解剖形态与力学特性模板库的个性化配准方法,首先需要获取大量健康人的骨骼CT图像进行形状与灰度的统计建模。图像配准得到的变形场描述了图像间对应像素的位置变化,使用骨骼所在区域的变形场可以描述骨骼形态的变化进而构建形状模型;CT图像的灰度值反映了骨骼对X射线的衰减系数,与骨骼表观密度具有近似的线性关系,同时表征了骨骼的材料力学特性,使用CT图像的灰度值对骨骼的灰度(密度)进行灰度建模。通过统计模型与患者CT图像之间个性化的配准方法可以得到描述骨骼信息的形状参数与灰度参数,与患者的假体模型一同构建包含解剖形态和材料力学特性的假体模板库。使用相似度公式评估病人CT图像与模板库中样本图像形状与灰度的相似程度,获得与病人骨骼相似度最高的样本图像与假体模型,作为病人的参考骨骼假体,方便医生进行假体植入物的设计与骨骼的材料力学特性分析。
发明内容
本发明提出一种骨骼CT图像的解剖形态和材料力学特性模板库的个性化配准方法,以大量健康人的骨骼CT图像构建能包含骨骼解剖形态与材料力学特性的统计模型,通过统计模型与骨骼CT图像之间个性化的配准方法实现对患者骨骼的参数化描述,使用患者图像与配准参数构建假体模板库,假体模板库通过个性化配准方法匹配病人骨骼CT图像,检索模板库中与病人骨骼情况相近的模板图像与假体模型,作为医生个性化假体植入物设计的初始参考模板。
为实现上述目的,本发明通过以下技术方案予以实现:一种骨骼CT图像的解剖形态和材料力学特性模板库的个性化配准方法,包括以下步骤:
S1.构建包含骨骼解剖形态的形状模型与材料力学特性的灰度模型的统计模型;S11、数据预处理;整理大量健康群体的骨骼CT图像为样本集,对样本集的像素分辨率和灰度范围进行归一化预处理;重采样处理统一空间分辨率,分辨率选取最精细的空间分辨率;采用统一的CT值范围对图像的灰度进行归一化处理以消除异常值的影响。
S12、解剖信息特征标定;所谓解剖特征信息,是指骨骼CT图像中能够显著反应骨骼形态结构变化、密度变化或具有医学意义的标志性信息,比如解剖标定点或骨骼分割的遮罩等。对样本集的所有样本图像进行解剖标定点的标注,记为L={(x1,y1,z1),…,(xi,yi,zi),…,(xk,yk,zk)},其中k为解剖标定点的个数,(xi,yi,zi)表示第i个点的三维坐标;
S13、构建平均模板图像;选取样本集的平均图像或任意样本图像作为初始模板,将所有样本通过图像配准的方式变形到平均图像的空间,计算出新的平均图像并替代初始模板图像,继续重复配准操作,直至平均模板图像不再发生变化;
S14、形状建模;通过图像配准得到图像间对应像素的位置变化,骨骼所在区域的变形场描述骨骼形态的变化,使用变形场构建骨骼的形状模型;具体步骤如下:
S141、计算形状向量;采用图像配准方法计算样本图像与平均模板图像之间的非线性空间变换,对每个样本构建形状向量,设si=[pi,li]∈R3N为样本i的形状向量,其中pi∈R3N为平均模板配准到样本i后得到的N个像素的三维坐标集合,li∈R3K为样本图像中K个解剖标定点的三维坐标集合;
S142、使用广义普氏分析(Generalized Procrustes Analysis,GPA)对所有样本的形状向量si进行对齐处理;
S143、使用主成分分析(Principal components analysis,PCA)方法对骨骼形状变化进行建模;形状模型如下所示:
Figure GDA0003556203860000031
其中,
Figure GDA0003556203860000032
表示样本集的平均形状,Φs∈R3(N+K)×M为对样本集做主成分分析处理后得到的特征向量矩阵,其中M个特征向量代表形状模型学习到的M种骨骼变形模式;bs∈RM为形状参数,bs中的M个元素分别代表Φs中的M个变形模式与平均形状叠加的权重;s∈R3(N+K)表示模型的形状,s受bs的控制,通过调节bs的值达到控制模型变形的目的。
S15、灰度建模;CT图像的灰度值反映器官对X射线的衰减系数,骨骼CT图像与周围软组织的高对比度精确描述骨骼的几何形态;
利用CT图像的灰度值对骨骼的灰度进行统计建模;对CT图像进行降采样后使用主成分分析方法对所有骨骼CT图像进行灰度建模,其公式如下:
Figure GDA0003556203860000041
其中,
Figure GDA0003556203860000042
表示样本集的平均灰度向量;Φg∈RN×M为对样本集做主成分分析处理后得到的统计模型中灰度模型的灰度特征矩阵,包含从样本集中学习到的M种灰度变化模式,下标g表示灰度;bg∈RM为灰度参数,g∈RN表示模型的灰度,g的值受灰度参数bg的控制,通过调节灰度参数bg的值达到控制模型灰度变化的目的;
S2.构建假体模板库;
S21、数据收集与处理根据经历骨骼假体植入样本前与植入后骨骼CT图像样本,建立植入前数据集D1与植入后数据集D2;D1包含了患者受损骨骼的形状与灰度信息,植入后数据集D2包含了个性化骨骼假体的形状与材质信息;对D1、D2进行像素分辨率和灰度范围归一化预处理,对D1中的图像进行解剖标定点的标注;
S22、形状模型匹配;将统计模型中的形状模型配准拟合到患者骨骼假体植入前的骨骼CT图像上,计算患者骨骼的形状参数bs;具体步骤如下:
S221、采样;对形状模型中平均形状的顶点进行等距采样,选择少量的采样点进行后续计算;
S222、个性化配准;通过图像的非线性配准方法得到平均模板到患者骨骼图像的映射关系即变形场,使用变形场对形状模型中平均形状的顶点进行变形,获得描述患者骨骼的形状向量;
非线性配准的损失函数使用图像间的灰度相似性测度与解剖标志点的距离测度作为约束条件,目标函数公式如下:
E=ωiEilEl
其中,E为模板与患者图像之间的加权相似测度,El为解剖标定点的距离测度;ωi和ωl为两个测度的权重,二者之和为1;
S223、将配准后的患者骨骼形状向量与平均形状向量进行对齐;其中,相似变换由缩放系数和旋转矩阵两部分组成,缩放系数通过计算配准后的骨骼形状向量与平均形状向量的顶点到质心之间距离的比值得出,旋转矩阵通过对两个形状向量的顶点矩阵进行奇异值分解计算得到;
S224、设置形变参数bs的个数;直接定义形变参数bs的个数或选择能够保留数据信息95%以上的形变参数bs个数;根据形变参数bs的个数从形状模型中取出对应的形状特征矩阵Φs
S225、计算患者骨骼的形状参数bs;根据从形状模型中取出的形状特征矩阵Φs以及配准后患者骨骼的形状向量,使用最小二乘法进行线性回归拟合,计算线性方程组Φsbs=sp得到患者骨骼的形状参数bs;其中,sp表示患者骨骼形状向量与平均形状向量之间的位移偏差,Φd表示统计模型中形状模型包含的形状特征矩阵。
S23、灰度模型匹配;将统计模型中的灰度模型拟合到患者骨骼假体植入前的骨骼图像上,计算患者骨骼的灰度参数bg
Φgbg=gp
其中,gp表示采样患者骨骼的灰度向量与平均灰度向量的标准差,Φg表示统计模型中灰度模型的灰度特征矩阵;
S24、获取骨骼假体模型;术后带有假体的CT图像使用阈值分割方法提取假体植入物的标签图,使用Marching Cubes算法将假体的标签图转化为三维曲面模型;金属假体植入物在CT图像中通常会呈现出较高的CT值,如金属铝的CT值在2000HU左右,不锈钢与钛合金在普通CT图像中会达到最高CT值,即3071HU。与人体内部器官的CT值相比有明显的区分,对术后带有假体的CT图像使用阈值分割方法提取假体植入物的标签图,然后使用MarchingCubes算法将假体的标签图转化为三维曲面模型。
S25、构建假体模板库;以统计模型、患者CT图像数据、患者骨骼的形状参数和灰度参数、假体模型为基础构建假体模板库;
假体模板库中存储的内容包含下表所示的部分:
假体模板库表格
Figure GDA0003556203860000061
由于形状特征矩阵、形状参数、平均形状向量、灰度特征矩阵、灰度参数、平均灰度向量的数据量太大无法直接存储在数据库列表中,可选择保存为二进制文件后将其存储路径保存在数据库中。
S3.模板库的个性化配准;通过假体模板库与CT图像的个性化配准得到骨骼的形状参数与灰度参数,检索假体模板库中与CT图像中骨骼最佳匹配的假体模型;
S31、数据预处理;获取需要进行骨骼假体植入手术的病人骨骼CT图像,对CT图像进行统一的空间与灰度归一化处理,对病人骨骼CT图像进行解剖标定点的标注;
S32、个性化配准与参数化表示;使用统计模型与病人骨骼CT图像进行个性化配准,通过形状模型与灰度模型的匹配计算得到病人骨骼的形状参数b′s与灰度参数b′g
S33、数据匹配;根据配准得到描述病人骨骼的形状参数与灰度参数,在假体模板库中使用相似度计算公式检索库中已存储的与病人骨骼近似的模板图像数据,相似度计算公式如下式:
Figure GDA0003556203860000071
其中,
Figure GDA0003556203860000072
为病人骨骼的形状参数与假体模板库中模板图像的形状参数之间的方差,衡量两幅图像间骨骼形状的相似程度;
Figure GDA0003556203860000073
为病人骨骼的灰度参数与假体模板库中模板图像的灰度参数之间的方差,衡量两幅图像之间骨骼灰度的相似程度;λ1和λ2分别为形状相似度与灰度相似度的权重,二者之和为1,根据不同的权重分配λ1和λ2的值,实现以形状相似度为主或以灰度相似度为主的检索匹配方式;
对假体模板库中的数据计算与病人数据间的相似度D,比较后得到与病人骨骼最相似的模板CT图像,获取该图像对应的骨骼假体模型作为病人骨骼的个性化骨骼假体的初始参考模板,用于后续骨骼假体的分析与设计。
图像图像图像图像图像图像图像
图像图像图像图像图像本发明的有益效果是:借助健康人骨骼CT图像构建统计模型,通过个性化的配准方法实现对骨骼CT图像的解剖形状与材料力学特征的参数化描述,依据现有数据构建假体模板库,通过模板库与病人图像的个性化配准找到模板库中与病人骨骼最佳匹配的假体模型,达到病人骨骼假体初始化设计的目的。另外,对于病人骨骼形状与灰度参数化表示形式,也可用于虚拟手术计划、骨骼的缺陷分析、个性化植入物的材料设计等方面。
附图说明
图1是模板库的个性化配准的总体工作流程图。
图2是统计模型的构建流程图。
图3是非线性配准方法的流程图。
图4(a)个性化配准平均模板图像;
图4(b)个性化配准患者图像;
图5(a)是患者一的CT图像;
图5(b)是患者一的假体模型图;
图5(c)是患者二的CT图像;
图5(d)是患者二的假体模型图;
图5(e)是患者三的CT图像;
图5(f)是患者三的假体模型图。
具体实施方式
以下髋骨为例,结合具体实施步骤对本发明做进一步说明,如图1所示,解剖形态和力学特性模板库的个性化配准方法,包括如下步骤:
步骤一、构建包含骨骼解剖形态与材料力学特性的统计模型。
S11、数据预处理。整理大量健康人的骨骼CT图像为样本集,对样本集进行归一化预处理,包括像素分辨率和灰度范围归一化。
对于来自不同医院的CT图像,采用统一的空间分辨率做重采样以消除空间分辨率的不一致,为确保不丢失图像信息,重采样的分辨率选取多家医院中最精细的空间分辨率;采用统一的CT值范围对图像的灰度进行归一化处理以消除异常值的影响。
S12、解剖信息特征标定。所谓解剖特征信息,是指骨骼CT图像中能够显著反应骨骼形态结构变化、密度变化或具有医学意义的标志性信息,比如解剖标定点或骨骼分割的遮罩等。
通过专家对样本集的所有样本图像进行解剖标定点的标注,记为L={(x1,y1,z1),…,(xi,yi,zi),…,(xk,yk,zk)},其中k为解剖标定点的个数,(xi,yi,zi)表示第i个点的三维坐标。
S13、构建平均模板图像。首先选取样本集的平均图像或任意样本图像作为初始模板,然后将所有样本通过图像配准的方式变形到平均图像的空间,计算出新的平均图像并替代初始模板图像,继续重复配准操作,循环迭代直至平均模板图像不再发生变化。
S14、形状建模。图像配准得到的变形场描述了图像间对应像素的位置变化,骨骼所在区域的变形场可以描述骨骼形态的变化,使用变形场构建骨骼的形状模型。
具体做法为:
S141、计算形状向量。采用图像配准方法计算样本图像与平均模板图像之间的非线性空间变换,对每个样本构建形状向量,设si=[pi,li]∈R3N为样本i的形状向量,其中pi∈R3N为平均模板配准到样本i后得到的所有像素(共N个)的三维坐标集合,li∈R3K为样本图像中所有解剖标定点(共K个)的三维坐标集合。
S142、广义普氏分析(Generalized Procrustes Analysis,GPA)。在使用形状向量si进行统计建模之前,需要消除其他因素引起的形状变化。所有的形状向量si需要在公共坐标系中经过旋转与平移进行对齐,使用GPA方法来对所有样本的形状向量si进行对齐处理。
S143、主成分分析(Principal components analysis,PCA)。对所有的形状向量si,使用PCA方法对骨骼形状变化进行建模,得到的形状模型如下所示:
Figure GDA0003556203860000101
其中,
Figure GDA0003556203860000102
表示样本集的平均形状,Φs∈R3(N+K)×M为对样本集做PCA处理后得到的特征向量矩阵,s表示形状,其中M个特征向量,M为样本数,代表形状模型学习到的M种骨骼变形模式;bs∈RM为形状参数,bs中的M个元素分别代表Φs中的M个变形模式与平均形状叠加的权重。s∈R3(N+K)表示模型的形状,s受bs的控制,通过调节bs的值达到控制模型变形的目的。
S15、灰度建模。CT图像的灰度值反映了器官对X射线的衰减系数,骨骼CT图像与周围软组织的高对比度能够精确描述骨骼的几何形态,同时,CT值与骨骼表观密度具有近似的线性关系,而骨骼表观密度与骨骼材料特性存在幂指数关系的经验公式,能够精确地描述骨骼的材料力学特性。
使用CT图像的灰度值可以对骨骼的灰度进行统计建模。采用相同的尺度标准对样本图像的像素进行降采样处理,可以减少后续步骤的计算量,提高灰度建模的速度。设gi∈RN为在样本i图像中采样得到的灰度向量,它代表了平均模板中各像素在样本i中对应位置处的图像灰度。以样本集所有的灰度向量进行建模,构建的模型表示如下:
Figure GDA0003556203860000111
其中
Figure GDA0003556203860000112
表示样本集的平均灰度向量,Φg∈RN×M为对样本集做PCA处理后得到的统计模型中灰度模型的灰度特征矩阵,g表示灰度,包含从样本集中学习到的M中灰度变化模式;bg∈RM为灰度参数,g∈RN表示模型的灰度,g的值受到灰度参数bg的控制,通过调节灰度参数bg的值达到控制模型灰度变化的目的。统计模型的构建流程如图2所示。
步骤二、个性化配准骨骼CT图像构建假体模板库。
S21、数据收集与处理。收集整理经历了骨骼假体植入手术患者的术前与术后的骨骼CT图像,建立术前数据集D1与术后数据集D2,其中术前数据集D1包含了患者受损骨骼的形态与灰度(密度)信息,术后数据集D2包含了个性化骨骼假体的形状与材质信息。参照步骤一中的S11方法,对数据集D1、D2进行像素分辨率和灰度范围归一化预处理,通过专家对术前数据集D1中的图像进行解剖标定点的标注。
S22、形状模型匹配。将统计模型中的形状模型配准拟合到患者术前的骨骼CT图像上,计算患者骨骼的形状参数bs
具体步骤包括:
S221、采样。对形状模型中平均形状的顶点进行等距采样,选择少量的采样点,如10000个顶点进行后续计算,减少计算的成本和时间;
S222、个性化配准。通过图像的非线性配准方法得到平均模板到患者骨骼图像的映射关系即变,使用变形场对形状模型中平均形状的顶点量
非线性配准的损失函数使用图像间的灰度相似性测度与解剖标志点的距离测度作为约束条件,使得模板能够正确匹配患者图像,目标函数公式如下:
Ei=ωiEilEl
患者图像间的灰度相似性测度(如互信息或互相关),基于灰度的配准帮助模板与患者骨骼CT图像的未受损部分完成对齐。El为解剖标定点的距离测度,用于确保模板与患者图像中解剖标志点之间的对齐,对骨骼受损部分的配准,其周围的解剖标志点可以为其提供引导,确保受损区域的配准不会发生过度的形变。ωi和ωl为两个测度的权重,其和为1。非线性配准方法的流程如图3所示。
S223、将配准后的患者骨骼形状向量与平均形状向量进行对齐,消除由相似变换带来的形状影响。其中,相似变换由缩放系数和旋转矩阵两部分组成,缩放系数通过计算配准后的骨骼形状向量与平均形状向量的顶点到质心之间距离的比值得出,旋转矩阵通过对两个形状向量的顶点矩阵进行奇异值分解计算得到。
S224、设置形变参数bs的个数。可以直接定义形变参数bs的个数,也可以选择能够保留数据信息95%以上的形变参数bs个数,根据形变参数bs的个数从形状模型中取出对应的形状特征矩阵Φg
S225、计算患者骨骼的形状参数bs。根据从形状模型中取出的形状特征矩阵Φs以及配准后患者骨骼的形状向量,使用最小二乘法进行线性回归拟合,计算线性方程组Φsbs=sp得到患者骨骼的形状参数bs,其中,sp表示患者骨骼形状向量与平均形状向量之间的位移偏差,Φs表示统计模型中形状模型包含的形状特征矩阵。
S23、灰度模型匹配。将统计模型中的灰度模型拟合到患者的骨骼图像上,计算患者骨骼的灰度参数bg。将描述患者骨骼形状的点云从物理空间坐标系转移到矩阵空间坐标系中,该点云同时表述了平均模板与患者图像像素间的映射关系,可以从患者图像中计算出平均模板配准到患者图像后的结果。
根据配准后图像的灰度向量与灰度模型中的灰度特征矩阵Φg,使用最小二乘法进行线性回归拟合,计算线性方程组Φgbg=gp得到患者骨骼的灰度参数bg,其中,gp表示采样后的患者骨骼的灰度向量与平均灰度向量的标准差,Φg表示统计模型中灰度模型的灰度特征矩阵。个性化配置准后患者骨骼的参数化表示如图4所示。
S24、获取骨骼假体模型。金属假体植入物在CT图像中通常会呈现出较高的CT值,如金属铝的CT值在2000HU左右,不锈钢与钛合金在普通CT图像中会达到最高CT值,即3071HU。与人体内部器官的CT值相比有明显的区分,对术后带有假体的CT图像使用阈值分割方法可以提取假体植入物的标签图,然后使用Marching Cubes算法将假体的标签图转化为三维曲面模型。术后患者的CT图像与假体模型如图5所示。
S25、构建假体模板库。以统计模型、患者CT图像数据、患者骨骼的形状参数和灰度参数、假体模型为基础构建假体模板库。假体模板库中存储的内容包含下表所示的部分:
假体模板库表格
Figure GDA0003556203860000131
Figure GDA0003556203860000141
由于形状特征矩阵、形状参数、平均形状向量、灰度特征矩阵、灰度参数、平均灰度向量的数据量太大无法直接存储在数据库列表中,可选择保存为二进制文件后将其存储路径保存在数据库中。
步骤三、假体模板库的个性化配准方法。
S31、数据预处理。获取需要进行骨骼假体植入手术的病人骨骼CT图像,参照步骤一中S11、S12方法,对CT图像进行统一的空间与灰度归一化处理,通过专家对病人骨骼CT图像进行解剖标定点的标注。
S32、个性化配准与参数化表示。参照步骤二中S22、S23方法,使用统计模型与病人骨骼CT图像进行个性化配准,通过形状模型与灰度模型的匹配计算得到病人骨骼的形状参数b′s与灰度参数b′g,。
S33、数据匹配。根据配准得到描述病人骨骼的形状参数与灰度参数,在假体模板库中使用相似度计算公式检索库中已存储的与病人骨骼近似的模板图像数据,相似度计算公式如下式:
Figure GDA0003556203860000142
其中,
Figure GDA0003556203860000143
为病人骨骼的形状参数与假体模板库中模板图像的形状参数之间的方差,衡量两幅图像间骨骼形状的相似程度;
Figure GDA0003556203860000144
为病人骨骼的灰度参数与假体模板库中模板图像的灰度参数之间的方差,衡量两幅图像之间骨骼灰度的相似程度;λ1和λ2分别为形状相似度与灰度相似度的权重,其和为1,根据不同的权重分配λ1和λ2的值,实现以形状相似度为主或以灰度相似度为主的检索匹配方式。
对假体模板库中的每例数据计算与病人数据间的相似度D值,比较后得到与病人骨骼最相似的模板CT图像,获取该图像对应的骨骼假体模型作为病人骨骼的个性化骨骼假体的初始参考模板,提供给医生来进行后续骨骼假体的分析与设计。
此处仅以髋骨配准为例加以说明,虽然人体其他骨骼的解剖形态与灰度变化情况各有不同,但本发明的思想和方法同样适用于其他骨骼。尽管已经示出和描述了本发明的实施例,对于本领域的普通技术人员而言,在不脱离本发明的原理和基本思想的情况下所做出的若干修改或改进均应视为本发明的保护范围。

Claims (3)

1.一种骨骼CT图像的解剖形态和材料力学特性模板库的个性化配准方法,其特征在于,包括以下步骤:
S1、构建包含骨骼解剖形态的形状模型与材料力学特性的灰度模型的统计模型;
S11、数据预处理;整理大量健康群体的骨骼CT图像为样本集,对样本集的像素分辨率和灰度范围进行归一化预处理;重采样处理统一空间分辨率,分辨率选取最精细的空间分辨率;
S12、解剖信息特征标定;对样本集的所有样本图像进行解剖标定点的标注,记为L={(x1,y1,z1),…,(xi,yi,zi),…,(xk,yk,zk)},其中k为解剖标定点的个数,(xi,yi,zi)表示第i个点的三维坐标;
S13、构建平均模板图像;选取样本集的平均图像或任意样本图像作为初始模板,将所有样本通过图像配准的方式变形到平均图像的空间,计算出新的平均图像并替代初始模板图像,继续重复配准操作,直至平均模板图像不再发生变化;
S14、形状建模;通过图像配准得到图像间对应像素的位置变化,骨骼所在区域的变形场描述骨骼形态的变化,使用变形场构建骨骼的形状模型;
S15、灰度建模;CT图像的灰度值反映器官对X射线的衰减系数,骨骼CT图像与周围软组织的高对比度精确描述骨骼的几何形态;
利用CT图像的灰度值对骨骼的灰度进行统计建模;对CT图像进行降采样后使用主成分分析方法对所有骨骼CT影像进行灰度建模,其公式如下:
Figure FDA0003556203850000011
其中,
Figure FDA0003556203850000012
表示样本集的平均灰度向量;Φg∈RN×M为对样本集做主成分分析处理后得到的统计模型中灰度模型的灰度特征矩阵,包含从样本集中学习到的M种灰度变化模式,下标g表示灰度;bg∈RM为灰度参数,g∈RN表示模型的灰度,g的值受灰度参数bg的控制,通过调节灰度参数bg的值达到控制模型灰度变化的目的;
S2.构建假体模板库;
S21、数据收集与处理根据经历骨骼假体植入样本前与植入后骨骼CT图像样本,建立植入前数据集D1与植入后数据集D2;D1包含了患者受损骨骼的形状与灰度信息,植入后数据集D2包含了个性化骨骼假体的形状与材质信息;对D1、D2进行像素分辨率和灰度范围归一化预处理,对D1中的图像进行解剖标定点的标注;
S22、形状模型匹配;将统计模型中的形状模型配准拟合到患者骨骼假体植入前的骨骼CT图像上,计算患者骨骼的形状参数bs
S23、灰度模型匹配;将统计模型中的灰度模型拟合到患者骨骼假体植入前的骨骼图像上,计算患者骨骼的灰度参数bg
Φgbg=gp,其中,gp表示采样患者骨骼的灰度向量与平均灰度向量的标准差,Φg表示统计模型中灰度模型的灰度特征矩阵;
S24、获取骨骼假体模型;术后带有假体的CT图像使用阈值分割方法提取假体植入物的标签图,使用Marching Cubes算法将假体的标签图转化为三维曲面模型;
S25、构建假体模板库;以统计模型、患者CT图像数据、患者骨骼的形状参数和灰度参数、假体模型为基础构建假体模板库;
S3.模板库的个性化配准;通过假体模板库与CT图像的个性化配准得到骨骼的形状参数与灰度参数,检索假体模板库中与CT图像中骨骼最佳匹配的假体模型;
S31、数据预处理;获取需要进行骨骼假体植入手术的病人骨骼CT图像,对CT图像进行统一的空间与灰度归一化处理,对病人骨骼CT图像进行解剖标定点的标注;
S32、个性化配准与参数化表示;使用统计模型与病人骨骼CT图像进行个性化配准,通过形状模型与灰度模型的匹配计算得到病人骨骼的形状参数b′s与灰度参数b′g
S33、数据匹配;根据配准得到描述病人骨骼的形状参数与灰度参数,在假体模板库中使用相似度计算公式检索库中已存储的与病人骨骼近似的模板图像数据,相似度计算公式如下式:
Figure FDA0003556203850000031
其中,
Figure FDA0003556203850000032
为病人骨骼的形状参数与假体模板库中模板图像的形状参数之间的方差,衡量两幅图像间骨骼形状的相似程度;
Figure FDA0003556203850000033
为病人骨骼的灰度参数与假体模板库中模板图像的灰度参数之间的方差,衡量两幅图像之间骨骼灰度的相似程度;λ1和λ2分别为形状相似度与灰度相似度的权重,二者之和为1,根据不同的权重分配λ1和λ2的值,实现以形状相似度为主或以灰度相似度为主的检索匹配方式;
对假体模板库中的数据计算与病人数据间的相似度D,比较后得到与病人骨骼最相似的模板CT图像,获取该图像对应的骨骼假体模型作为病人骨骼的个性化骨骼假体的初始参考模板,用于后续骨骼假体的分析与设计。
2.根据权利要求1所述的骨骼CT图像的解剖形态和材料力学特性模板库的个性化配准方法,其特征在于,S14中形状建模包括步骤如下:
S141、计算形状向量;采用图像配准方法计算样本图像与平均模板图像之间的非线性空间变换,对每个样本构建形状向量,设si=[pi,li]∈R3N为样本i的形状向量,其中pi∈R3N为平均模板配准到样本i后得到的N个像素的三维坐标集合,li∈R3K为样本图像中K个解剖标定点的三维坐标集合;
S142、使用广义普氏分析对所有样本的形状向量si进行对齐处理;
S143、使用主成分分析方法对骨骼形状变化进行建模;形状模型如下所示:
Figure FDA0003556203850000041
其中,
Figure FDA0003556203850000042
表示样本集的平均形状,Φs∈R3(N+K)×M为对样本集做主成分分析处理后得到的特征向量矩阵,其中M个特征向量代表形状模型学习到的M种骨骼变形模式;bs∈RM为形状参数,bs中的M个元素分别代表Φs中的M个变形模式与平均形状叠加的权重;s∈R3(N+K)表示模型的形状,s受bs的控制,通过调节bs的值达到控制模型变形的目的。
3.根据权利要求1或2所述的骨骼CT图像的解剖形态和材料力学特性模板库的个性化配准方法,其特征在于,S22形状模型匹配中,包括步骤如下:
S221、采样;对形状模型中平均形状的顶点进行等距采样,选择少量的采样点进行后续计算;
S222、个性化配准;通过图像的非线性配准方法得到平均模板到患者骨骼图像的映射关系即变形场,使用变形场对形状模型中平均形状的顶点进行变形,获得描述患者骨骼的形状向量;
非线性配准的损失函数使用图像间的灰度相似性测度与解剖标志点的距离测度作为约束条件,目标函数公式如下:
E=ωiEilEl
其中,E为模板与患者图像之间的加权相似测度,Ei为模板与患者图像间的灰度相似性测度,El为解剖标定点的距离测度;ωi和ωl为两个测度的权重,二者之和为1;
S223、将配准后的患者骨骼形状向量与平均形状向量进行对齐;其中,相似变换由缩放系数和旋转矩阵两部分组成,缩放系数通过计算配准后的骨骼形状向量与平均形状向量的顶点到质心之间距离的比值得出,旋转矩阵通过对两个形状向量的顶点矩阵进行奇异值分解计算得到;
S224、设置形变参数bs的个数;直接定义形变参数bs的个数或选择能够保留数据信息95%以上的形变参数bs个数;根据形变参数bs的个数从形状模型中取出对应的形状特征矩阵Φs
S225、计算患者骨骼的形状参数bs;根据从形状模型中取出的形状特征矩阵Φs以及配准后患者骨骼的形状向量,使用最小二乘法进行线性回归拟合,计算线性方程组Φsbs=sp得到患者骨骼的形状参数bs;其中,sp表示患者骨骼形状向量与平均形状向量之间的位移偏差,Φs表示统计模型中形状模型包含的形状特征矩阵。
CN202111317598.3A 2021-11-09 2021-11-09 一种骨骼ct图像的解剖形态和材料力学特性模板库的个性化配准方法 Active CN114155286B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN202111317598.3A CN114155286B (zh) 2021-11-09 2021-11-09 一种骨骼ct图像的解剖形态和材料力学特性模板库的个性化配准方法
US17/841,222 US20230146649A1 (en) 2021-11-09 2022-06-15 Personalized registration method for template library of anatomical morphology and mechanical properties of materials of bone ct images

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111317598.3A CN114155286B (zh) 2021-11-09 2021-11-09 一种骨骼ct图像的解剖形态和材料力学特性模板库的个性化配准方法

Publications (2)

Publication Number Publication Date
CN114155286A CN114155286A (zh) 2022-03-08
CN114155286B true CN114155286B (zh) 2022-05-10

Family

ID=80459617

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111317598.3A Active CN114155286B (zh) 2021-11-09 2021-11-09 一种骨骼ct图像的解剖形态和材料力学特性模板库的个性化配准方法

Country Status (2)

Country Link
US (1) US20230146649A1 (zh)
CN (1) CN114155286B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115292828B (zh) * 2022-07-12 2023-04-18 中国人民解放军总医院第四医学中心 一种治疗干骺端骨折解剖型钢板及其形态设计方法及装置
CN116029161B (zh) * 2023-03-27 2023-09-08 北京爱康宜诚医疗器材有限公司 一种具有渐变孔隙率的假体设计方法
CN116452755B (zh) * 2023-06-15 2023-09-22 成就医学科技(天津)有限公司 一种骨骼模型构建方法、系统、介质及设备
CN116503453B (zh) * 2023-06-21 2023-09-26 福建自贸试验区厦门片区Manteia数据科技有限公司 图像配准方法、装置、计算机可读存储介质及电子设备
CN117152256B (zh) * 2023-10-30 2024-02-13 中国人民解放军总医院第一医学中心 一种基于模板的骨盆模型通道定位方法及装置

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006099538A (ja) * 2004-09-30 2006-04-13 Hitachi Ltd 形状モデル作成方法、それを実装したコンピュータプログラムおよび形状モデル作成システム
EP2901419A1 (en) * 2012-09-27 2015-08-05 Siemens Product Lifecycle Management Software Inc. Multi-bone segmentation for 3d computed tomography
CN109003283A (zh) * 2018-03-26 2018-12-14 天津工业大学 一种基于主动形状模型的主动脉轮廓分割算法
CN110148128A (zh) * 2019-05-23 2019-08-20 中南大学 一种补全病变骨骼以获得骨骼预期参考模型的方法
CN110335358A (zh) * 2019-06-18 2019-10-15 大连理工大学 可变形数字人解剖学模型的个性化变形方法
CN112233801A (zh) * 2020-12-17 2021-01-15 季华实验室 内置假体拓扑优化数学模型构建方法及拓扑优化设计方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006099538A (ja) * 2004-09-30 2006-04-13 Hitachi Ltd 形状モデル作成方法、それを実装したコンピュータプログラムおよび形状モデル作成システム
EP2901419A1 (en) * 2012-09-27 2015-08-05 Siemens Product Lifecycle Management Software Inc. Multi-bone segmentation for 3d computed tomography
CN109003283A (zh) * 2018-03-26 2018-12-14 天津工业大学 一种基于主动形状模型的主动脉轮廓分割算法
CN110148128A (zh) * 2019-05-23 2019-08-20 中南大学 一种补全病变骨骼以获得骨骼预期参考模型的方法
CN110335358A (zh) * 2019-06-18 2019-10-15 大连理工大学 可变形数字人解剖学模型的个性化变形方法
CN112233801A (zh) * 2020-12-17 2021-01-15 季华实验室 内置假体拓扑优化数学模型构建方法及拓扑优化设计方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Planning of skull reconstruction based on a statistical shape model combined with geometric morphometrics;Marc Anton Fuessinger等;《Springer》;20171028;第519-529页 *
基于分形维数特征的肺结节形状建模;赵海等;《东北大学学报(自然科学版)》;20181130;第39卷(第11期);第1545-1550页 *
面向模板影像学习的医学影像数据库构建研究;翟浩宇;《万方数据库》;20201208;第1-60页 *

Also Published As

Publication number Publication date
US20230146649A1 (en) 2023-05-11
CN114155286A (zh) 2022-03-08

Similar Documents

Publication Publication Date Title
CN114155286B (zh) 一种骨骼ct图像的解剖形态和材料力学特性模板库的个性化配准方法
US10068671B2 (en) Methods and systems for producing an implant
Lamecker et al. Atlas-based 3D-shape reconstruction from X-ray images
Claes et al. Computerized craniofacial reconstruction: conceptual framework and review
Bryan et al. Statistical modelling of the whole human femur incorporating geometric and material properties
CN113674281B (zh) 一种基于深度形状学习的肝脏ct自动分割方法
Sigal et al. Mesh-morphing algorithms for specimen-specific finite element modeling
CN109493308A (zh) 基于条件多判别生成对抗网络的医疗图像合成与分类方法
Ibragimov et al. Segmentation of pathological structures by landmark-assisted deformable models
Kainmueller et al. An articulated statistical shape model for accurate hip joint segmentation
CN114494183B (zh) 一种基于人工智能的髋臼半径自动测量方法及系统
AU2020101836A4 (en) A method for generating femoral x-ray films based on deep learning and digital reconstruction of radiological image
Pereañez et al. Accurate segmentation of vertebral bodies and processes using statistical shape decomposition and conditional models
CN112037200A (zh) 一种医学影像中解剖特征自动识别与模型重建方法
WO2024001140A1 (zh) 一种椎体亚区域分割方法、装置及存储介质
CN106327479A (zh) 血管造影中介下先心病术中血管辨识的装置及方法
de Buhan et al. A facial reconstruction method based on new mesh deformation techniques
CN107980149A (zh) 用于脊椎标记的方法、装置和系统
CN113838048B (zh) 一种十字交叉韧带术前止点中心定位及韧带长度计算方法
CN113515875A (zh) 基于多模态图像的骨生物力学建模方法、系统及装置
Berar et al. 3D statistical facial reconstruction
CN111080676B (zh) 一种通过在线分类跟踪内窥镜图像序列特征点的方法
CN117350143A (zh) 一种心脏房室瓣置换术后力学环境的评估方法
CN115252233B (zh) 基于深度学习的全髋关节置换术中髋臼杯的自动化规划方法
Staring et al. Nonrigid registration using a rigidity constraint

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