CN114913149A - 一种基于ct影像的头部可变形统计图谱构建方法 - Google Patents

一种基于ct影像的头部可变形统计图谱构建方法 Download PDF

Info

Publication number
CN114913149A
CN114913149A CN202210512328.6A CN202210512328A CN114913149A CN 114913149 A CN114913149 A CN 114913149A CN 202210512328 A CN202210512328 A CN 202210512328A CN 114913149 A CN114913149 A CN 114913149A
Authority
CN
China
Prior art keywords
shape
training sample
head
anatomical structure
anatomical
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
Application number
CN202210512328.6A
Other languages
English (en)
Other versions
CN114913149B (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.)
Yancheng Institute of Technology
Yancheng Institute of Technology Technology Transfer Center Co Ltd
Original Assignee
Yancheng Institute of Technology
Yancheng Institute of Technology Technology Transfer Center 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 Yancheng Institute of Technology, Yancheng Institute of Technology Technology Transfer Center Co Ltd filed Critical Yancheng Institute of Technology
Priority to CN202210512328.6A priority Critical patent/CN114913149B/zh
Publication of CN114913149A publication Critical patent/CN114913149A/zh
Application granted granted Critical
Publication of CN114913149B publication Critical patent/CN114913149B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/14—Transformations for image registration, e.g. adjusting or mapping for alignment of images
    • G06T3/147—Transformations for image registration, e.g. adjusting or mapping for alignment of images using affine transformations
    • 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/136—Segmentation; Edge detection involving thresholding
    • 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/344—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods involving models
    • 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/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30008—Bone
    • 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/30016—Brain
    • 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/30088—Skin; Dermal
    • 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/30196—Human being; Person
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2210/00—Indexing scheme for image generation or computer graphics
    • G06T2210/41—Medical
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2210/00—Indexing scheme for image generation or computer graphics
    • G06T2210/44—Morphing

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Quality & Reliability (AREA)
  • Radiology & Medical Imaging (AREA)
  • Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • General Health & Medical Sciences (AREA)
  • Computer Graphics (AREA)
  • Geometry (AREA)
  • Software Systems (AREA)
  • Image Processing (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

本发明提供了一种基于CT影像的头部可变形统计图谱构建方法，包括：S1：基于最大熵阈值分割法对每个训练样本对应的头部CT影像进行分割并转化，获得每个训练样本中不同解剖结构对应的目标网格；S2：获取头部参考模板和对应的3D多边形网格；S3：基于解剖标志点约束的薄板样条鲁棒点匹配算法，采用轮廓匹配的分层式配准策略，将所述3D多边形网格与每个训练样本中不同解剖结构对应的目标网格进行配准，获得每个训练样本对应的个性化解剖结构表示；S4：基于所述个性化解剖结构表示和统计形状建模方法构建出头部可变形统计图谱；用以采用中国人群CT影像作为数据集，提出基于统计建模构建包含人群个体头部之间解剖结构差异的头部可变形统计图谱的方法。

Description

一种基于CT影像的头部可变形统计图谱构建方法
技术领域
本发明涉及医学影像处理和断层影像解剖学技术领域，特别涉及一种基于CT影像的头部可变形统计图谱构建方法。
背景技术
人体的头部是一个复杂的多功能结构系统。当前，对头部的研究和实践大都基于头部数字图谱展开。头部数字图谱是指运用信息技术，将头部的解剖结构信息数字化和可视化，建立的多层次、数字化的头部模型。头部数字图谱的出现，为头部的相关研究提供了精准、完备的解剖学先验知识和仿真平台。头部数字图谱广泛应用于神经外科、神经放射学、神经病学和神经教育学等领域。现有的头部图谱大多是基于多边形网格来表示其解剖结构，以便于控制图谱的曲面形状。Nowinski等人构建了包含头部肌肉、颅内血管系统、白质束、颅神经、腺体、颅骨和颅外血管系统的头部数字图谱。Okubo等人构建了名为“BodyParts3D”的人体多边形网格图谱，包括完整的头部解剖结构。现阶段，大多数可用的头部数字图谱都是基于一个特定的个体开发的。然而，随着个性化医疗的发展对头部数字图谱提出了个性化的需求，即要求图谱能够反映个体的解剖差异。为了满足这一要求，Lee等人利用50例不同个体的头部MR影像构建了人群头部模型库，库中包含了50例个体不同的解剖特征。但是，这个头部模型库仅仅整理了50例个体，而没有对这50例个体之间的解剖差异进行统计建模。
从以上说明来看，现有头部数字解剖图谱大多基于单一的、特定的人体解剖形态构建，在医学研究和解剖学教育领域有着广泛的应用。然而，在针对不同个体的个性化诊疗、仿真和人体工程学设计等领域，囿于现有头部数字解剖图谱缺乏个体之间解剖形态差异，难以准确表达不同个体的解剖形态特征，从而限制了其在个性化建模中的应用
因此，本发明提出了一种基于CT影像的头部可变形统计图谱构建方法。
发明内容
本发明提供一种基于CT影像的头部可变形统计图谱构建方法，用以采用中国人群CT影像作为数据集，提出基于统计建模构建包含人群个体头部之间解剖结构差异的头部可变形统计图谱的方法，可以应用于头部个性化模拟和个性化医疗领域，如模拟电磁辐射剂量评估、整形美容手术制定和模拟、放射治疗路径规划和导航等。
本发明提供一种基于CT影像的头部可变形统计图谱构建方法，包括：
S1：基于最大熵阈值分割法对每个训练样本对应的头部CT影像进行分割并转化，获得每个训练样本中不同解剖结构对应的目标网格；
S2：获取头部参考模板和对应的3D多边形网格；
S3：基于解剖标志点约束的薄板样条鲁棒点匹配算法，采用轮廓匹配的分层式配准策略，将所述3D多边形网格与每个训练样本中不同解剖结构对应的目标网格进行配准，获得每个训练样本对应的个性化解剖结构表示；
S4：基于所述个性化解剖结构表示和统计形状建模方法构建出头部可变形统计图谱。
优选的，所述的一种基于CT影像的头部可变形统计图谱构建方法，S1：基于最大熵阈值分割法对每个训练样本对应的头部CT影像进行分割并转化，获得每个训练样本中不同解剖结构对应的目标网格，包括：
S101：在海量头部CT影像中筛选出无症状个体头部CT影像作为训练样本；
S102：基于最大熵阈值分割法对所述训练样本对应的头部CT影像进行分割，获得对应的可分割解剖结构和未分割解剖结构；
S103：基于等值面提取方法将所述可分割解剖结构和所述未分割解剖结构转化为对应的三角面片网格，并将对应三角面片网格作为对应解剖结构的目标网格。
优选的，所述的一种基于CT影像的头部可变形统计图谱构建方法，S2：获取头部参考模板和对应的3D多边形网格，包括：
S201：获取通用人体数字模型的头部模型作为所述头部参考模板；
S202：基于人体MR影像提取出所述头部参考模板的解剖结构信息；
S203：基于所述解剖结构信息构建出对应的3D多边形网格。
优选的，所述的一种基于CT影像的头部可变形统计图谱构建方法，S3：基于解剖标志点约束的薄板样条鲁棒点匹配算法，采用轮廓匹配的分层式配准策略，将所述3D多边形网格与每个训练样本中不同解剖结构对应的目标网格进行配准，获得每个训练样本对应的个性化解剖结构表示，包括：
基于解剖标志点约束的薄板样条鲁棒点匹配算法，采用轮廓匹配的分层式配准策略，将所述3D多边形网格中包含的第一解剖结构与对应训练样本中包含的第二解剖结构对应的目标网格进行配准，获得对应训练样本的整体配准结果；
基于预设拓扑结构网格表示所述整体配准结果，获得每个训练样本对应的个性化解剖结构表示。
优选的，所述的一种基于CT影像的头部可变形统计图谱构建方法，基于解剖标志点约束的薄板样条鲁棒点匹配算法，采用轮廓匹配的分层式配准策略，将所述3D多边形网格中包含的第一解剖结构与对应训练样本中包含的第二解剖结构对应的目标网格进行配准，获得对应训练样本的整体配准结果，包括：
基于所述解剖标志点约束的薄板样条鲁棒点匹配算法，将所述3D多边形网格中包含的第一可分割解剖结构与所述训练样本中包含的对应第二可分割解剖结构进行点云配准，获得对应的第一配准结果；
基于所述第一配准结果确定出所述头部参考模板至对应训练样本的空间形变场；
确定出所述3D多边形网格中包含的所有第一解剖结构之间的第一相对位置信息和所述训练样本中包含的所有第二解剖结构之间的第二相对位置信息；
基于所述训练样本对应的空间形变场和对应的第二相对位置信息以及所述第一相对位置信息，将所述3D多边形网格中包含的第一未分割参考解剖结构与对应训练样本中包含的第二未分割解剖结构进行二次配准，获得对应训练样本的第二配准结果；
综合所述训练样本对应的第一配准结果和第二配准结果，获得对应训练样本的整体配准结果。
优选的，所述的一种基于CT影像的头部可变形统计图谱构建方法，基于所述解剖标志点约束的薄板样条鲁棒点匹配算法，将所述3D多边形网格中包含的第一可分割解剖结构与所述训练样本中包含的对应第二可分割解剖结构进行点云配准，获得对应的第一配准结果，包括：
人工标定出所述第一可分割解剖结构中包含的第一关键解剖标志点与所述训练样本中包含的对应第二可分割解剖结构中包含的第二关键解剖标志点；
确定出所述第一可分割解剖结构的第一点云表示和对应第二可分割解剖结构的第二点云表示；
基于所述第一点云表示和初始仿射变换矩阵以及初始非线性变形系数矩阵，确定出对应的初始非刚性变换函数；
基于所述初始非刚性变换函数确定出所述第一点云表示对应的初始非刚性变换映射；
基于预设的温度系数和所述第二点云表示以及所述初始非刚性变换映射确定出对应的权重系数，基于所述权重系数和所述第二点云表示确定出对应的权重系数矩阵；
基于所述第一关键解剖标志点、所述第二关键解剖标志点、所述第一点云表示、所述第二点云表示、所述权重系数矩阵和预设的正则化参数构建出对应的代价函数，将所述代价函数最小化，确定出新的非刚性变换函数；
基于最新确定的非刚性变换函数和所述第一点云表示确定出所述第二可分割解剖结构对应的与所述第一点云表示配准的第三点云表示，判断所述第二点云表示中每个点和所述第三点云表示中对应点之间的距离是否都小于预设阈值，若是，则将所述第三点云表示作为所述3D多边形网格中包含的第一可分割解剖结构对应的第一配准结果；
否则，基于预设减小梯度设置新的温度系数和新的正则化参数，基于新的温度系数和最新确定的非刚性变换函数确定出新的权重系数矩阵，基于新的权重系数矩阵和新的正则化参数构建出新的代价函数，将新的代价函数最小化，确定出新的非刚性变换函数，基于最新确定的非刚性变换函数确定出所述第二可分割解剖结构对应的与所述第一点云表示配准的第四点云表示，直至所述第二点云表示中每个点和所述第四点云表示中对应点之间的距离都小于预设阈值时，则将所述第四点云表示作为所述3D多边形网格中包含的第一可分割解剖结构对应的第一配准结果。
优选的，所述的一种基于CT影像的头部可变形统计图谱构建方法，基于预设拓扑结构网格表示所述整体配准结果，获得每个训练样本对应的个性化解剖结构表示，包括：
基于所述整体配准结果获得所述3D多边形网格中每个顶点至对应训练样本中对应顶点的映射关系；
基于预设拓扑结构网格表示所述整体配准结果，获得每个训练样本对应的个性化解剖结构表示。
优选的，所述的一种基于CT影像的头部可变形统计图谱构建方法，S4：基于所述个性化解剖结构表示和统计形状建模方法构建出头部可变形统计图谱，包括：
基于广义普氏分析对所述个性化解剖结构中所述训练样本对应的整体配准结果中的形状网格进行空间方位归一化，获得所述训练样本对应的归一化形状网格；
基于主成分分析方法和所述整体配准结果，提取所述归一化形状网格中的统计形变分量，基于所述统计形变分量构建出对应的可变形图谱。
优选的，所述的一种基于CT影像的头部可变形统计图谱构建方法，基于广义普氏分析对所述训练样本对应的整体配准结果中的形状网格进行空间方位归一化，获得所述训练样本对应的归一化形状网格，包括：
基于所述形状网格中包含的所有点在参考坐标系下对应的第一坐标值，计算出所述形状网格的形状中心对应的第二坐标值；
基于所述形状中心对应的第二坐标值对所述形状网格中包含的所有点进行去中心化处理，获得对应的原点对齐形状网格；
基于所述形状网格中包含的所有点在参考坐标系下对应的第一坐标值和所述形状中心对应的第二坐标值计算出对应训练样本的形状尺寸测度；
将所述原点对齐形状网格中每个点对应的第三坐标值除以对应的形状尺寸测度，获得对应的尺度归一化形状网络；
在所有训练样本对应的尺度归一化形状网络中任选一个尺度归一化形状网络作为标准形状；
将所有训练样本对应的尺度归一化形状网络中除所述标准形状以外剩余的待旋转形状与所述标准形状进行对齐，获得新的训练样本；
计算出新的训练样本对应的平均形状，计算出所述平均形状和所述参考形状之间的普氏距离平方值；
当所述普氏距离平方值大于设定阈值时，则将平均形状设置为新的参考形状，将所有训练样本对应的尺度归一化形状网络中除新的参考形状以外剩余的待旋转形状与新的参考形状进行对齐，直至最新获得的平均形状与最新的参考形状之间的普氏距离不大于设定阈值时，则将最新获得的训练样本作为对应的归一化形状网格。
优选的，所述的一种基于CT影像的头部可变形统计图谱构建方法，基于主成分分析方法和所述整体配准结果，提取所述归一化形状网格中的统计形变分量，基于所述统计形变分量构建出对应的可变形图谱，包括：
计算出所有归一化形状网格的平均形状网格和协方差矩阵；
对所述协方差矩阵进行特征值分解，获得对应的由特征向量构成的正交矩阵和由特征值构成的对角矩阵；
将所述对角矩阵中包含的特征值从大到小排序获得对应的特征值序列；
将所述特征值序列中前n个特征值当作对应所述可变形图谱的数据集，将所述数据集中包含的每个第一特征值在所述正交矩阵中对应的第一特征向量当作对应归一化形状网格的统计形变分量；
基于所有归一化形状网格对应的权重变形参数和对应的统计形变分量的乘积与所述平均形状网格的线性组合，构建出对应的可变形图谱。
本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在所写的说明书、权利要求书、以及附图中所特别指出的结构来实现和获得。
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：
图1为本发明实施例中一种基于CT影像的头部可变形统计图谱构建方法流程图；
图2为本发明实施例中又一种基于CT影像的头部可变形统计图谱构建方法流程图；
图3为本发明实施例中再一种基于CT影像的头部可变形统计图谱构建方法流程图；
图4为本发明实施例中头部可变形统计图谱的构建流程图；
图5为本发明实施例中头部可变形统计图谱包含的男性头部解剖结构图；
图6为本发明实施例中头部可变形统计图谱包含的女性头部解剖结构图；
图7为本发明实施例中头部可变形统计图谱中与脂肪量变化相关的变形模式的女性头部脂肪量变化示意图；
图8为本发明实施例中头部可变形统计图谱中与脂肪量变化相关的变形模式的男性头部脂肪量变化示意图；
图9为本发明实施例中头部可变形统计图谱头部形状变化相关的变形模式的头部形状变化示意图；
图10为本发明实施例中头部可变形统计图谱头部形状变化相关的变形模式的脑颅形状变化示意图；
图11为本发明实施例中头部可变形统计图谱骨骼测量的R指数示意图；
图12为本发明实施例中头部可变形统计图谱活体测量的Z指数示意图；
图13为本发明实施例中头部可变形统计图谱配准到男性测试个体CT影像的配准结果示意图；
图14为本发明实施例中头部可变形统计图谱配准到女性测试个体CT影像的配准结果示意图。
具体实施方式
以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。
实施例1：
本发明提供了一种基于CT影像的头部可变形统计图谱构建方法，参考图1和4，包括：
S1：基于最大熵阈值分割法对每个训练样本对应的头部CT影像进行分割并转化，获得每个训练样本中不同解剖结构对应的目标网格；
S2：获取头部参考模板和对应的3D多边形网格；
S3：基于解剖标志点约束的薄板样条鲁棒点匹配算法，采用轮廓匹配的分层式配准策略，将所述3D多边形网格与每个训练样本中不同解剖结构对应的目标网格进行配准，获得每个训练样本对应的个性化解剖结构表示；
S4：基于所述个性化解剖结构表示和统计形状建模方法构建出头部可变形统计图谱。
该实施例中，最大熵阈值分割法即为选择合适的阈值使得各个子集的信息量最大。通过阈值法对头部的皮肤和颅骨进行分割，接着对皮肤和颅骨之间的区域进行阈值分割来分离面部肌肉和皮下脂肪，经过阈值分割，得到皮肤、颅骨、面部肌肉和皮下脂肪等头部解剖结构，其中，最大熵准则的核心在于将图像的灰度直方图分成独立的几个子集，使得各个子集的总熵最大。
该实施例中，训练样本即为无症状的个体头部CT影像。
该实施例中，不同解剖结构包括皮肤、颅骨、面部肌肉和皮下脂肪等头部解剖结构。
该实施例中，目标网格即为基于最大熵阈值分割法对每个训练样本对应的头部CT影像进行分割并转化后获得的三角面片网格。
该实施例中，头部参考模板即为后续用于与训练样本对应的头部CT影像进行配准的参考的头部模板。
该实施例中，3D多边形网格即为基于人体MR影像提取头部参考模板的解剖结构信息构建出的三维多边形网格。
该实施例中，个性化解剖结构表示即为将3D多边形网格与每个训练样本中不同解剖结构对应的目标网格进行配准后获得的每个训练样本对应的解剖结构表示(用同一网格拓扑表示的不同训练个体的解剖结构)。
该实施例中，统计形状建模方法即为通过SSM(statistical shape model，统计模型)方法建模，对训练样本进行形状分析，学习头部在人群中的形变方式，最终获得训练个体之间的解剖差异。
该实施例中，头部可变形统计图谱既是标准的头部解剖数字图谱，又具有个性化变形能力的头部解剖数字图谱。
该实施例中，CT(Computed Tomography，电子计算机断层扫描)影像，即电子计算机断层扫描影像。
该实施例中，本发明首先，将不同训练样本的头部CT影像分割成对应的解剖结构，接着将人体头部解剖结构的参考模板网格配准到每个分割后的影像，获得每个训练样本的个性化解剖结构表示。最后，基于配准到训练样本的参考模板，采用统计形状建模的方法构建头部DSA(deformable statistical atlas，头部可变形统计图谱)。
该实施例中，参考图4，展示了头部DSA构建的原理流程；首先，将不同个体的头部CT影像分割成对应的解剖结构；接着将人体头部解剖结构的参考模板网格配准到每个分割后的影像，获得每个个体的个性化解剖结构表示；最后，基于配准到训练样本的参考模板，采用统计形状建模的方法构建头部可变形统计图谱。
以上技术的有益效果为：以中国人群健康成人的CT影像为基础，构建了中国人群的男、女头部DSA(deformable statistical atlas,可变形统计图谱)，头部DSA包含了从真实个体的医学影像中获取的实际的解剖差异，使用统计形状建模的方法从训练样本中学习真实的解剖变化，构建的头部DSA可以通过调整形变参数匹配到个体头部，进而得到个体头部精确的解剖结构，可以应用于个性化模拟和个性化医疗，如模拟电磁辐射剂量评估、整形美容手术制定和模拟、放射治疗路径规划和导航等。
实施例2：
在实施例1的基础上，所述的一种基于CT影像的头部可变形统计图谱构建方法，S1：基于最大熵阈值分割法对每个训练样本对应的头部CT影像进行分割并转化，获得每个训练样本中不同解剖结构对应的目标网格，参考图2和5以及6，包括：
S101：在海量头部CT影像中筛选出无症状个体头部CT影像作为训练样本；
S102：基于最大熵阈值分割法对所述训练样本对应的头部CT影像进行分割，获得对应的可分割解剖结构和未分割解剖结构；
S103：基于等值面提取方法将所述可分割解剖结构和所述未分割解剖结构转化为对应的三角面片网格，并将对应三角面片网格作为对应解剖结构的目标网格。
该实施例中，可分割解剖结构即为皮肤、颅骨、面部肌肉和皮下脂肪等可分割的解剖结构。
该实施例中，未分割解剖结构即为脑组织等未分割的解剖结构。
该实施例中，基于等值面提取方法将所述可分割解剖结构和所述未分割解剖结构转化为对应的三角面片网格即为利用移动立体(marching cubes)算法，采用间接的等值面提取方法将分割后的解剖结构转化为三角面片网格。
该实施例中，目标网格即为下一步参考模板对应的网格配准到训练样本对应的个体的目标网格。
该实施例中，首先筛选无症状的个体头部CT影像作为训练数据，对于获取的训练影像，需要进行图像分割，得到尽可能准确的头部各个解剖结构，图像分割的作用是从图像中提取出感兴趣区域或者将图像细分为构成它的子区域；由于头部CT影像是典型的灰度图像，因此本发明采用基于最大熵准则的阈值分割方法对训练图像进行分割；最大熵准则的核心在于将图像的灰度直方图分成独立的几个子集，使得各个子集的总熵最大；从信息论角度出发，就是选择合适的阈值使得各个子集的信息量最大；通过阈值法对头部的皮肤和颅骨进行分割，接着对皮肤和颅骨之间的区域进行阈值分割来分离面部肌肉和皮下脂肪；经过阈值分割，得到皮肤、颅骨、面部肌肉和皮下脂肪等头部解剖结构。最后，利用移动立体(marching cubes)算法，采用间接的等值面提取方法将分割后的解剖结构转化为三角面片网格，作为下一步模板网格配准到训练样本对应的个体的目标网格。
该实施例中，参考图5和6，展示了头部可变形统计图谱包含的男性和女性解剖结构；第一到第三列展示了头部可变形统计图谱正面、侧面和背面三个不同视角呈现的所有解剖结构的半透明效果图；第四至第七列分别展示了头部可变形统计图谱皮肤、肌肉、颅骨和大脑等结构。
以上技术的有益效果为：基于最大熵阈值分割法实现尽可能获得准确的头部各个解剖结构，基于等值面提取方法将所述可分割解剖结构和所述未分割解剖结构转化为对应的三角面片网格，为下一步与参考模板对应的网格进行配准提供了目标网格，也为生成可变形统计图谱提供了初步数据基础。
实施例3：
在实施例2的基础上，所述的一种基于CT影像的头部可变形统计图谱构建方法，S2：获取头部参考模板和对应的3D多边形网格，参考图3和图7到图12包括：
S201：获取通用人体数字模型的头部模型作为所述头部参考模板；
S202：基于人体MR影像提取出所述头部参考模板的解剖结构信息；
S203：基于所述解剖结构信息构建出对应的3D多边形网格。
该实施例中，通用人体数字模型即为BodyParts3D(是日本东京大学生命科学数据库中心的Okubo教授团队研发的关于人体详细解剖结构的3D多边形网格数字模板。这个人体模板采用通用的解剖学坐标体系)，该模板符合亚洲人特征，因此选择BodyParts3D作为头部参考模板。
该实施例中，人体MR(magnetic resonance，核磁共振)影像即为核磁共振影像。
该实施例中，解剖结构信息即为基于人体MR影像提取出头部参考模板解剖结构相关的信息。
该实施例中，参考图7到12，图7和8中展示了头部可变形统计图谱中与脂肪量变化相关的变形模式女性和男性；不同的列对应于不同的形状参数值；在前视图和侧视图中，第一行使用不透明皮肤渲染头部可变形统计图谱，第二行和第三行分别使用透明皮肤渲染头部可变形统计图谱；箭头指向增加的面部脂肪和颈部后部脂肪；图7和8展示出了男性和女性头部可变形统计图谱前三个主成分对应的典型解剖变形模式；其中，男性和女性头部可变形统计图谱显示出相似的变形模式；在图7和8中，女性头部可变形统计图谱的形状变化模式形状变化模式1和男性头部可变形统计图谱的形状变化模式3对应于脂肪量的变化。当女性头部可变形统计图谱的形状参数1或男性头部可变形统计图谱的形状参数3增加时，在皮肤不透明视图中两个头部可变形统计图谱的脸都变得更胖，而在皮肤透明视图中面部脂肪和颈部后部脂肪都变的更厚；
图9和10中展示了头部形状变化相关的变形模式头部形状变化，头部的长度和面部比例变化情况；前视图中的虚线线条说明面部比例的变化，侧视图中的箭头指向女性面部突出脑颅形状变化；图9和10展示了头部可变形统计图谱整体形状的变化情况；男性和女性头部可变形统计图谱的形状变化模式2都对应于头部长度和面部比例的变化，具体如图9和10所示。当形变参数2变大时，整个脸的长度增加，但是前额的高度却降低了；这说明脸部前额以下部分的增长是导致脸部总长度增加的原因；对于女性头部可变形统计图谱，形状变化模式2也对应于面部突出的程度(如图9和10侧视图中箭头所示)，然而，这种变化在男性头部可变形统计图谱中并不明显；图9和10显示了脑颅部分的形状变化；当女性头部可变形统计图谱的形变参数3或男性头部可变形统计图谱形变参数1增长时，脑颅部分的宽度减小(如前视图所示)，而横向比例增加(如侧视图所示)；
图11和12展示了头部可变形统计图谱骨骼测量的R指数和活体测量的Z指数；图中每个测量项的3个结果从左至右分别对应头部可变形统计图谱的前3种变形模式的测量值；圆圈表示平均形状，误差条对应于i的取值范围，测量项的名称用马丁编码表示；通过在范围内改变可变形统计图谱形状参数i的取值，计算每种变形模式相应的R指数和Z指数；男性和女性的头部DSAs的前三种变形模式分别占所有变形模式的43.60％和49.98％；图11和12中活体测量和骨骼测量的测量项名称使用马丁编码表示；为了对可变形统计图谱进行活体测量和骨骼测量，在可变形统计图谱的皮肤和颅骨上手动标定相应的解剖标志点，然后对其进行活体测量和骨骼测量；从图11和12可以看出，可变形统计图谱的骨骼测量项的R指数在[-0.2,0.2]范围内，而除了测量项M7之外，其它测量项都落在[-0.1,0.1]的范围内；结果表明可变形统计图谱与中国人群平均值的偏差小于平均值的10％；除了测量项M8之外，可变形统计图谱的活体测量项的Z指数都位于[-1.96,1.96]范围内；通常认为中国人群的统计数据服从正态分布，则位于[-1.96,1.96]范围内的Z指数表示可变形统计图谱属于中国人群总体分布的概率大于95％；图11和12所示的结果表明，可变形统计图谱的测量值与中国成人人群统计数据基本一致。
以上技术的有益效果为：本发明采用通用人体数字模型BodyParts3D的头部作为参考模板可以补偿CT影像中分辨率低的软组织区域，并且人体数字模型BodyParts3D的头部符合亚洲人特征，从而提高了生成的可变形统计图谱的准确性。
实施例4：
在实施例3的基础上，所述的一种基于CT影像的头部可变形统计图谱构建方法，S3：基于解剖标志点约束的薄板样条鲁棒点匹配算法，采用轮廓匹配的分层式配准策略，将所述3D多边形网格与每个训练样本中不同解剖结构对应的目标网格进行配准，获得每个训练样本对应的个性化解剖结构表示，包括：
基于解剖标志点约束的薄板样条鲁棒点匹配算法，采用轮廓匹配的分层式配准策略，将所述3D多边形网格中包含的第一解剖结构与对应训练样本中包含的第二解剖结构对应的目标网格进行配准，获得对应训练样本的整体配准结果；
基于预设拓扑结构网格表示所述整体配准结果，获得每个训练样本对应的个性化解剖结构表示。
该实施例中，基于轮廓匹配的分层式配准策略的核心思想是利用可分割解剖结构部分与未分割解剖结构部分的位置相对关系，借助可分割解剖结构部分配准过程获得的空间形变场，把未分割解剖结构部分的形状网格映射到训练样本对应的个体，从而得到未分割解剖结构部分的配准结果。
该实施例中，第一解剖结构即为头部参考模板对应的3D多边形网格中包含的解剖结构。
该实施例中，第二解剖结构即为训练样本中包含的解剖结构。
该实施例中，整体配准结果即为基于解剖标志点约束的薄板样条鲁棒点匹配算法，采用轮廓匹配的分层式配准策略，将3D多边形网格中包含的第一解剖结构与对应训练样本中包含的第二解剖结构对应的目标网格进行配准后，获得的对应训练样本中包含的每个第二解剖结构在3D多边形网格中对应的第一解剖结构。
该实施例中，基于预设拓扑结构网格表示所述整体配准结果，获得每个训练样本对应的个性化解剖结构表示，即为：通过把已有模板配准到训练样本对应的个体，得到个体的模板变形表征，然后将所有训练样本对应的个体都用同样拓扑结构的网格表示，具有相同的顶点数目，并且每一个顶点在不同的个体曲面中代表相同的解剖位置，这样就得到了训练样本对应的个体之间的解剖结构对应关系。
以上技术的有益效果为：基于轮廓匹配的分层式配准策略，将已有的头部参考模板对应的3D多边形网格与训练样本对应的头部CT影像对应的目标网格进行配准，从第一层可分割解剖结构配准到第二层未分割解剖结构配准，最终实现所有解剖结构配准，实现将头部参考模板中包含的解剖结构和训练样本对应的头部CT影像中包含的解剖结构进行一一对应配准，进一步提高了生成的可变形统计图谱的准确性。
实施例5：
在实施例4的基础上，所述的一种基于CT影像的头部可变形统计图谱构建方法，基于解剖标志点约束的薄板样条鲁棒点匹配算法，采用轮廓匹配的分层式配准策略，将所述3D多边形网格中包含的第一解剖结构与对应训练样本中包含的第二解剖结构对应的目标网格进行配准，获得对应训练样本的整体配准结果，参考图13和14，包括：
基于所述解剖标志点约束的薄板样条鲁棒点匹配算法，将所述3D多边形网格中包含的第一可分割解剖结构与所述训练样本中包含的对应第二可分割解剖结构进行点云配准，获得对应的第一配准结果；
基于所述第一配准结果确定出所述头部参考模板至对应训练样本的空间形变场；
确定出所述3D多边形网格中包含的所有第一解剖结构之间的第一相对位置信息和所述训练样本中包含的所有第二解剖结构之间的第二相对位置信息；
基于所述训练样本对应的空间形变场和对应的第二相对位置信息以及所述第一相对位置信息，将所述3D多边形网格中包含的第一未分割参考解剖结构与对应训练样本中包含的第二未分割解剖结构进行二次配准，获得对应训练样本的第二配准结果；
综合所述训练样本对应的第一配准结果和第二配准结果，获得对应训练样本的整体配准结果。
该实施例中，第一可分割解剖结构即为3D多边形网格中包含的可分割解剖结构。
该实施例中，第二可分割解剖结构即为训练样本中包含的可分割解剖结构。
该实施例中，第一配准结果即为解剖标志点约束的薄板样条鲁棒点匹配算法将3D多边形网格中包含的第一可分割解剖结构与训练样本中包含的对应第二可分割解剖结构进行点云配准后获得的配准结果。
该实施例中，空间形变场即为基于第一配准结果确定出的头部参考模板至对应训练样本的空间形变映射关系。
该实施例中，第一相对位置信息即为3D多边形网格中包含的所有第一解剖结构之间的相对位置信息。
该实施例中，第二相对位置信息即为训练样本中包含的所有第二解剖结构之间的相对位置信息。
该实施例中，第二配准结果即为基于训练样本对应的空间形变场和对应的第二相对位置信息以及所述第一相对位置信息，将3D多边形网格中包含的第一未分割参考解剖结构与对应训练样本中包含的第二未分割解剖结构进行二次配准后，获得的对应训练样本的第二未分割解剖结构在3D多边形网格中对应的第一未分割参考解剖结构。
该实施例中，第一关键解剖标志点即为第一解剖结构中包含的关键解剖标志点，关键解剖标志点即为手动标定的有解剖学意义的特定的解剖结构位置或关键的解剖结构特征点。
该实施例中，第二关键解剖标志点即为第二解剖结构中包含的关键解剖标志点，关键解剖标志点即为手动标定的有解剖学意义的特定的解剖结构位置或关键的解剖结构特征点。
该实施例中，已知训练样本对应的头部CT影像(即训练个体)和头部参考模板，需要把头部参考模板中所有解剖结构配准到个体中，则包括：
个体的解剖结构用ASS(anatomical structure of subject，ASS)表示，参考模板的解剖结构用AST(anatomical structure of template,AST)表示；
ASS1和ASS2代表个体中已分割的解剖结构，ASS3代表个体中未分割的解剖结构；AST1和AST2以及AST3表示头部参考模板中与个体中ASS1和ASS2以及ASS3一一对应的解剖结构；参考模板具有解剖结构的先验知识，包含AST1和AST2以及AST3的轮廓形状和解剖空间中的相对位置信息；
首先将头部参考模板中的AST1和AST2配准到个体的ASS1和ASS2，得到模板配准到个体的空间形变场；利用获得的空间形变场，对头部参考模板中的AST3进行变换，从而得到AST3到个体中的ASS3的整体配准结果。
该实施例中，参考图13和14，展示了覆盖在测试个体上的配准后头部可变形统计图谱；头部可变形统计图谱的配准结果为目标CT影像的解剖结构提供了一个合理的估计，特别是对于CT影像中对比度较低、边界不清晰的脑组织部分，同样给出了比较准确的描述；头部可变形统计图谱的这一特性可以用于放射治疗路径规划和导航、头部电磁辐射评估、整形手术等针对个性化个体的应用场景，实现对个体的准确模拟和仿真测试。
以上技术的有益效果为：在头部CT影像的配准过程中，先采用解剖标志点约束的薄板样条鲁棒点匹配算法把参考模板配准到个体可分割的解剖结构，接着利用轮廓匹配方法配准得到的空间形变场将参考模板配准到个体未分割的解剖结构，从而实现参考模板到个体整个头部的配准，整个的配准过程采用的基于轮廓形状的配准方法，从第一层可分割解剖结构配准到第二层未分割解剖结构配准，最终实现所有解剖结构配准，也进一步提高了生成的可变形统计图谱的准确性。
实施例6：
在实施例5的基础上，所述的一种基于CT影像的头部可变形统计图谱构建方法，基于所述解剖标志点约束的薄板样条鲁棒点匹配算法，将所述3D多边形网格中包含的第一可分割解剖结构与所述训练样本中包含的对应第二可分割解剖结构进行点云配准，获得对应的第一配准结果，包括：
人工标定出所述第一可分割解剖结构中包含的第一关键解剖标志点与所述训练样本中包含的对应第二可分割解剖结构中包含的第二关键解剖标志点；
确定出所述第一可分割解剖结构的第一点云表示和对应第二可分割解剖结构的第二点云表示；
基于所述第一点云表示和初始仿射变换矩阵以及初始非线性变形系数矩阵，确定出对应的初始非刚性变换函数；
基于所述初始非刚性变换函数确定出所述第一点云表示对应的初始非刚性变换映射；
基于预设的温度系数和所述第二点云表示以及所述初始非刚性变换映射确定出对应的权重系数，基于所述权重系数和所述第二点云表示确定出对应的权重系数矩阵；
基于所述第一关键解剖标志点、所述第二关键解剖标志点、所述第一点云表示、所述第二点云表示、所述权重系数矩阵和预设的正则化参数构建出对应的代价函数，将所述代价函数最小化，确定出新的非刚性变换函数；
基于最新确定的非刚性变换函数和所述第一点云表示确定出所述第二可分割解剖结构对应的与所述第一点云表示配准的第三点云表示，判断所述第二点云表示中每个点和所述第三点云表示中对应点之间的距离是否都小于预设阈值，若是，则将所述第三点云表示作为所述3D多边形网格中包含的第一可分割解剖结构对应的第一配准结果；
否则，基于预设减小梯度设置新的温度系数和新的正则化参数，基于新的温度系数和最新确定的非刚性变换函数确定出新的权重系数矩阵，基于新的权重系数矩阵和新的正则化参数构建出新的代价函数，将新的代价函数最小化，确定出新的非刚性变换函数，基于最新确定的非刚性变换函数确定出所述第二可分割解剖结构对应的与所述第一点云表示配准的第四点云表示，直至所述第二点云表示中每个点和所述第四点云表示中对应点之间的距离都小于预设阈值时，则将所述第四点云表示作为所述3D多边形网格中包含的第一可分割解剖结构对应的第一配准结果。
该实施例中，求解AL-TPS-RPM算法的过程是反复迭代求最优f(·)的过程。在算法初始化时，一般设d为单位矩阵，ω为0矩阵，那么f(·)就有了初始值；接着，通过式(4)、(5)和(6)求出yi，带入式(1)，求得新的yi；接着进行新一轮的迭代运算，更新f(·)，直至收敛。最后，得到TPS变换函数f(·)；AL-TPS-RPM算法的过程如下：
初始化参数T和λ1以及λ2；
初始化权重参数矩阵M，参数d和ω；
开始A：确定性退火；
开始B：交替更新：
第一步:使用(5)和(6)更新对应权重参数矩阵M；
第二步:使用(1)更新变换参数d和ω；
结束B；
减小T和λ；
结束A。
该实施例中，第一点云表示即为第一解剖结构的点云表示，即为目标网格中的所有顶点。
该实施例中，第二点云表示即为第二解剖结构的点云表示，在第一步中，已经对所有个体CT影像利用marching cubes算法将分割后的结构转化为三角面片网格，只要提取网格的所有顶点，得到一系列离散的点就是对应解剖结构在三维空间中的点云。
该实施例中，假设头部参考模板对应的第一点云表示为Te＝*t1,t2,···,tn}，包含n个点，ti∈Te(1≤i≤n)；训练个体对应的第二点云表示为Su＝*s1,s2,···,sm}，包含m个点，sj∈Su(1≤j≤m)；人工在第二点云表示中标定的解剖标志点表示为AL＝*al1,al2,···,alk}；用函数f(·)表示AL-TPS-RPM的非刚性变换，使得第一点云表示中的点ti被映射到一个新的位置处f(ti)；引入算子L来作为映射的平滑度约束，具体光滑度测度形式为‖Lf(ti)‖2；为了得到点云匹配的最优非刚性变换函数f(·)；最小化下面AL-TPS-RPM的代价函数：
Figure BDA0003638485180000201
Figure BDA0003638485180000211
其中，E(f)为代价函数，i为第一点云表示中包含的当前计算的点云，n为第一点云表示中包含的点云总个数，yi为ti在Su中所对应的目标点表示，f(ti)为第一点云表示中的点ti的非刚性变换映射，‖yi-f(ti)‖为yi与f(ti)之间的距离(头部参考模板中的点云和训练个体对应点云之间的距离)，人工在第二点云表示中标定的解剖标志点的点云表示，‖si-f(ti)‖为si与f(ti)之间的距离(代表了参考模板中的第一关键解剖标志点和训练个体中的第二关键解剖标志点之间的距离)，L为作为映射的平滑度约束的算子，‖Lf(ti)‖2为TPS变换函数的平滑约束项，λ1和λ2为两个正则化参数，预设正则化参数为初始设置的，d为4×4的仿射变换矩阵，
Figure BDA0003638485180000212
和IAL(ti∈AL)是指示函数，I为指示函数
Figure BDA0003638485180000213
和IAL(ti∈AL)，trace[[d-I]T[d-I]]为[d-I]T,d-I]的迹，,d-I]为d和I做差后获得的矩阵，,d-I]T为,d-I]的转置矩阵；
最小化AL-TPS-RPM的代价函数E(f)，即为寻找最优的TPS变换函数，使得模板和个体的轮廓点云最匹配；
其中，
Figure BDA0003638485180000214
和IAL(ti∈AL)的定义为：
Figure BDA0003638485180000215
Figure BDA0003638485180000216
yi具体的计算公式如下：
Figure BDA0003638485180000217
其中，mij构成的矩阵即为权重系数矩阵，权重系数矩阵为初始设置的，sj∈Su(1≤j≤m)为训练个体对应的第二点云表示Su＝*s1,s2,···,sm}中的第j个点；
式(4)根据sj与f(ti)之间的距离对Su中的点进行加权求和，得到ti在Su中所对应的目标点yi，从而将目标点云Su转化为与参考模板Te中的点一一对应的点云Ys＝*y1,y2,···,yn}，以便简化AL-TPS-RPM的代价函数，把模板准确地配准到个体；
Figure BDA0003638485180000221
其中，
Figure BDA0003638485180000222
其中，T温度系数(用来控制点云匹配的精度，温度系数越小，点云匹配的精度越高)，(si-f(ti))T(si-f(ti))用于度量si与f(ti)之间的距离，权重系数mij越小，exp()为以自然常数e为底的指数函数，j为第一点云在第二点云中对应的当前计算的目标点，m为第一点云在第二点云中对应的目标点总个数；
其中，权重系数m′ij构成的矩阵即为权重系数矩阵；
TPS变换函数是基于TPS插值的非线性空间变换，当式(1)中正则化参数λ1和λ2固定时，存在由两个矩阵d和ω确定的唯一的极小值函数：
f(ti)＝tid+φ(ti)·ω； (7)
其中，d为4×4的仿射变换矩阵，ω为n×4的非线性变形系数矩阵，φ(ti)为点ti的1×n的TPS算子向量，φ(ti)包含了第一点云表示点集内部结构关系的信息，可以通过ti和Te求出，φ(ti)当与非线性变形系数矩阵ω相结合时会产生非刚性的变形；
因此，TPS变换函数完全由仿射变换矩阵d和非刚性变形矩阵ω决定，d和ω的值确定了，TPS变换函数就确定了。
该实施例中，初始仿射变换矩阵即为4×4矩阵，在算法初始化时，一般设初始仿射变换矩阵为单位矩阵。
该实施例中，初始非线性变形系数矩阵即为0矩阵。
该实施例中，基于所述第一点云表示和初始仿射变换矩阵以及初始非线性变形系数矩阵，确定出对应的初始非刚性变换函数，包括：
将初始仿射变换矩阵和初始非线性变形系数矩阵代入至公式(7)，即确定出对应的初始非刚性变换函数。
该实施例中，基于所述初始非刚性变换函数确定出所述第一点云表示对应的初始非刚性变换映射，包括：
假设头部参考模板对应的第一点云表示为Te＝*t1,t2,···,tn}，包含n个点，ti∈Te(1≤i≤n)；用函数f(·)表示AL-TPS-RPM的非刚性变换，使得第一点云表示中的点ti被映射到一个新的位置处f(ti)，f(ti)即为第一点云表示中的点ti对应的非刚性变换映射。
该实施例中，基于预设的温度系数和所述第二点云表示以及所述初始非刚性变换映射确定出对应的权重系数，基于所述权重系数和所述第二点云表示确定出对应的权重系数矩阵，包括：
将预设的温度系数和所述第二点云表示以及所述初始非刚性变换映射代入至公式(5)和公式(6)，求得的mij即为对应的权重系数；
权重系数m′ij构成的矩阵即为权重系数矩阵。
该实施例中，基于所述第一关键解剖标志点、所述第二关键解剖标志点、所述第一点云表示、所述第二点云表示、所述权重系数矩阵和预设的正则化参数构建出对应的代价函数，将所述代价函数最小化，确定出新的非刚性变换函数，包括：
将第一关键解剖标志点、第二关键解剖标志点、第一点云表示、第二点云表示、权重系数矩阵、预设的正则化参数、最新确定的非刚性变换映射代入至公式(1)至公式(4)，进而构建出对应的代价函数，将构建出的代价函数最小化，确定出新的仿射变换矩阵和新的非线性变形系数矩阵，基于新的仿射变换矩阵和新的非线性变形系数矩阵确定出新的非刚性变换函数。
该实施例中，基于最新确定的非刚性变换函数和所述第一点云表示确定出所述第二解剖结构对应的与所述第一点云表示配准的第三点云表示，包括：
将第一点云表示代入最新确定的非刚性变换函数确定出第二解剖结构对应的与所述第一点云表示配准的第三点云表示。
该实施例中，预设阈值即为点云配准结果满足要求时对应的整体配准结果中训练个体每个点和头部参考模板中对应点之间的最大距离。
该实施例中，预设减小梯度即为在点云匹配的迭代过程中预设减小梯度和正则化参数在每一次迭代过程中预先设置的较小梯度。
该实施例中，基于新的温度系数和最新确定的非刚性变换函数确定出新的权重系数矩阵，包括：
将新的温度系数和最新确定的非刚性变换函数代入至公式(5)和公式(6)确定出新的权重系数矩阵。
该实施例中，基于新的权重系数矩阵和新的正则化参数构建出新的代价函数，包括：
将新的权重系数矩阵和新的正则化参数代入至公式(1)至公式(4)，进而构建出新的代价函数。
该实施例中，基于最新确定的非刚性变换函数确定出所述第二解剖结构对应的与所述第一点云表示配准的第四点云表示，包括：
将第一点云表示代入至最新确定的非刚性变换函数确定出第二解剖结构对应的与第一点云表示配准的第四点云表示。
该实施例中，提出在现有参考模板和训练样本的点云基础上，手动标定有解剖学意义的关键解剖标志点，通过特定的解剖结构位置或具有解剖意义的关键的解剖结构特征点来引导和约束点云轮廓匹配的过程；解剖标志点是指在解剖学上有意义的点，指示特定的解剖位置。通常解剖标志点可以保证指定的解剖结构属于同一器官或同一物种；通过这些解剖标志点，确保配准过程中参考模板形状网格与训练样本之间的关键解剖结构位置的准确对应。
该实施例中，提出基于解剖标志点约束的薄板样条鲁棒点匹配(anatomicallandmark thin-plate spline robust point matching，AL-TPS-RPM)算法实现BodyParts3D模板配准到个体CT影像；AL-TPS-RPM是在TPS-RPM算法基础上加入解剖标志点引导配准过程，通过控制点弯曲函数实现参考模板形状点云配准到训练样本的形状点云；AL-TPS-RPM是一种基于TPS的非刚性空间映射参数化自动点云匹配的配准方法，算法鲁棒性强，自动化程度高，对初始值的设定不敏感。TPS是一种通用样条函数工具，它可以生成一个通过所有控制点并且弯曲最小、平滑的函数映射。选择TPS参数化非刚性映射是因为它是可以明确地分解为全局仿射分量和局部非仿射分量的样条，同时基于空间映射的二阶导数最小化弯曲能量。
以上技术的有益效果为：通过这些解剖标志点，确保配准过程中参考模板形状网格与训练样本之间的关键解剖结构位置的准确对应，实现从头部模板配准到个体CT影像，配准到训练个体的头部模板网格具有平滑的表面弥补了前一个器官分割步骤中产生的分割不平滑的问题。
实施例7：
在实施例6的基础上，所述的一种基于CT影像的头部可变形统计图谱构建方法，基于预设拓扑结构网格表示所述整体配准结果，获得每个训练样本对应的个性化解剖结构表示，包括：
基于所述整体配准结果获得所述3D多边形网格中每个顶点至对应训练样本中对应顶点的映射关系；
基于预设拓扑结构网格表示所述整体配准结果，获得每个训练样本对应的个性化解剖结构表示。
该实施例中，对于个体CT影像中未分割的大脑结构，利用颅骨配准得到的空间变换函数将参考模板中的大脑映射到个体CT图像中；采用TPS变换将整个头部参考模板配准到每个个体的影像空间，得到参考模板到所有个体头部结构的映射；将模板配准到所有个体之后，得到模板网格的每个顶点到不同训练个体的相同解剖位置的映射。这样，就可以用同一网格拓扑表示不同训练个体的解剖结构。
以上技术的有益效果为：尽管在前一个器官分割步骤中，阈值法的分割结果会产生不平滑的边缘，但是配准到训练个体的头部模板网格具有平滑的表面弥补了分割不平滑的问题，且实现了用同一网格拓扑表示不同训练个体的解剖结构。
实施例8：
在实施例7的基础上，所述的一种基于CT影像的头部可变形统计图谱构建方法，S4：基于所述个性化解剖结构表示和统计形状建模方法构建出头部可变形统计图谱，包括：
基于广义普氏分析对所述个性化解剖结构中所述训练样本对应的整体配准结果中的形状网格进行空间方位归一化，获得所述训练样本对应的归一化形状网格；
基于主成分分析方法和所述整体配准结果，提取所述归一化形状网格中的统计形变分量，基于所述统计形变分量构建出对应的可变形图谱。
该实施例中，基于统计形状模型的方法构建头部可变形图谱过程大体分两步：首先，应用广义普氏分析(Generalized Procrustes Analysis,GPA)对训练个体配准后的形状网格规范化，对所有头部训练样本在空间位置和方向上做归一化处理。然后利用主成分分析(principal component analysis,PCA)方法提取头部的统计形变分量，构造点分布模型类型的可变形图谱。
该实施例中，不同个体的头部在空间位置和形状上差异较大，在构建统计形状模型前，需要先对所有个体点云做GPA处理，建立一个参考坐标系(参考标准形状)，把所有个体的头部的姿势对齐，去掉个体头部本身的位置和方向影响，最终实现所有个体的形状点云对齐；GPA方法和Procrustes分析方法的原理基本是一致的，区别在于确定个体的参考方向，GPA是最佳确定的，而Procrustes分析方法是任意选择的；也就是说，GPA应用Procrustes分析方法来优化叠加一组对象，而不是将它们叠加到任意选择的形状上；两种方法都以相同的方式进行缩放和转换；当只处理两种形状时，GPA相当于普通的Procrustes分析；
通常Procrustes分析的形状对齐分为四个步骤：
(1)计算每个形状的中心；
(2)对齐位置，去中心化处理；
(3)缩放每个形状，使其大小相等；
(4)通过旋转对齐方向；
而在Procrustes分析方法的基础上，GPA方法的在形状对齐做了改变，具体分为四个步骤：
(1)在训练集中任意选择参考形状；
(2)将训练集中的所有形状与参考形状对齐；
(3)计算对齐后的训练集新的平均形状；
(4)如果平均形状和参考形状之间的Procrustes距离高于设定的阈值，将平均形状设置为新的参考形状，然后继续执行步骤(2)；
在构建头部可变形图谱时，头部的尺寸差异是人群之间的典型解剖形态差异，因此本发明采用部分Procrustes分析操作，在GPA过程中不对个体做尺寸归一化，即保留个体的尺寸信息，只做空间位置和方向的归一化；经过GPA处理之后，求得所有n个个体的归一化结果。
该实施例中，归一化形状网格即为基于广义普氏分析对个性化解剖结构中训练样本对应的整体配准结果中的形状网格进行空间方位归一化后获得的训练样本对应的归一化后的形状网格。
该实施例中，统计形变分量即为基于主成分分析方法和整体配准结果在归一化形状网络中提取出的表征对应归一化形状网格与头部参考模板对应的3D多边形网格之间的形变分量。
以上技术的有益效果为：基于广义普氏分析对整体配准结果中的形状网格进行空间方位归一化，获得所述训练样本对应的归一化形状网格，为后续提取出对应个体的统计形变分量并构建出可变形图谱提供了关键基础。
实施例9：
在实施例8的基础上，所述的一种基于CT影像的头部可变形统计图谱构建方法，基于广义普氏分析对所述训练样本对应的整体配准结果中的形状网格进行空间方位归一化，获得所述训练样本对应的归一化形状网格，包括：
基于所述形状网格中包含的所有点在参考坐标系下对应的第一坐标值，计算出所述形状网格的形状中心对应的第二坐标值；
基于所述形状中心对应的第二坐标值对所述形状网格中包含的所有点进行去中心化处理，获得对应的原点对齐形状网格；
基于所述形状网格中包含的所有点在参考坐标系下对应的第一坐标值和所述形状中心对应的第二坐标值计算出对应训练样本的形状尺寸测度；
将所述原点对齐形状网格中每个点对应的第三坐标值除以对应的形状尺寸测度，获得对应的尺度归一化形状网络；
在所有训练样本对应的尺度归一化形状网络中任选一个尺度归一化形状网络作为标准形状；
将所有训练样本对应的尺度归一化形状网络中除所述标准形状以外剩余的待旋转形状与所述标准形状进行对齐，获得新的训练样本；
计算出新的训练样本对应的平均形状，计算出所述平均形状和所述参考形状之间的普氏距离平方值；
当所述普氏距离平方值大于设定阈值时，则将平均形状设置为新的参考形状，将所有训练样本对应的尺度归一化形状网络中除新的参考形状以外剩余的待旋转形状与新的参考形状进行对齐，直至最新获得的平均形状与最新的参考形状之间的普氏距离不大于设定阈值时，则将最新获得的训练样本作为对应的归一化形状网格。
该实施例中，基于所述形状网格中包含的所有点在参考坐标系下对应的第一坐标值，计算出所述形状网格的形状中心对应的第二坐标值，包括：
个体Sur1的形状中心可以通过对点云中所有点求和取平均得到：
Figure BDA0003638485180000281
式中，
Figure BDA0003638485180000282
为个体Sur1的形状中心，i为个体Sur1中包含的当前计算的点云，n为个体Sur1中包含的所有点云的总个数，Sur1i为个体Sur1中包含的第i个点云。
该实施例中，第一坐标值即为形状网格中包含的所有点在参考坐标系下对应的坐标值。
该实施例中，第二坐标值即为形状网格的形状中心对应的坐标值。
该实施例中，基于所述形状中心对应的第二坐标值对所述形状网格中包含的所有点进行去中心化处理，获得对应的原点对齐形状网格，包括：
将个体Sur1点云中所有点减去点云的中心
Figure BDA0003638485180000291
将点云中心由其均值平移至坐标原点，即可实现个体去中心化；
个体Sur1的形状尺寸测度定义为：
Figure BDA0003638485180000292
式中，Ssm(Sur1)为个体Sur1的形状尺寸测度，i为个体Sur1中包含的当前计算的点云，n为个体Sur1中包含的所有点云的总个数，Sur1i为个体中包含的第i个点云，
Figure BDA0003638485180000293
为Sur1i和
Figure BDA0003638485180000294
之间的距离；
将个体Sur1中所有点除以它的形状尺寸测度，即可实现个体尺度归一化，即获得对应的原点对齐形状网格。
该实施例中，标准形状即为在所有训练样本对应的尺度归一化形状网络中任选的一个尺度归一化形状网络。
该实施例中，新的训练样本即为将所有训练样本对应的尺度归一化形状网络中除所述标准形状以外剩余的待旋转形状与所述标准形状进行对齐后获得的形状网络。
该实施例中，计算出新的训练样本对应的平均形状，计算出所述平均形状和所述参考形状之间的普氏距离平方值，包括：
通常用Procrustes距离(Procrustes distance，普氏距离)的平方度量两个个体形状之间的差异；Procrustes距离是一种最小二乘型形状度量测度，它要求两个要对齐的形状具有一对一的点对应关系；Procrustes分析的过程，是通过最小化形状差异测度即Procrustes距离来获得最相似的空间位置、尺寸和方向；对于个体Sur1和Sur2，其Procrustes距离的平方是点云中所有点距离平方的总和：
Figure BDA0003638485180000301
其中，
Figure BDA0003638485180000302
为普氏距离平方值，Sur1为平均形状对应的个体的点云表示，Sur2为参考形状对应的个体的点云表示，i为平均形状对应的个体或参考形状对应的个体中包含的当前计算的点云，n为平均形状对应的个体或参考形状对应的个体中包含的点云总个数，Sur1i为平均形状对应的个体中包含的第i个点云，Sur2i为参考形状对应的个体中包含的第i个点云，‖Sur1-Sur2‖为Sur1和Sur2之间的距离，‖Sur1i-Sur2i‖为Sur1i和Sur2i之间的距离；
对于两个个Sur1和Sur2，假设以Sur1为标准形状，个体Sur2经过平移和尺寸变换之后，固定参考对象，然后围绕原点旋转个体Sur2，使得Sur1和Sur2的Procrustes距离的平方不大于设定阈值。
该实施例中，在Procrustes分析的形状对齐方法的基础上在形状对齐步骤进行改变，具体步骤包括：
(1)计算每个形状的中心；
(2)对齐位置，去中心化处理；
(3)缩放每个形状，使其大小相等；
(4)在训练集中任意选择参考形状；
(5)将训练集中的所有形状与参考形状对齐；
(6)计算对齐后的训练集新的平均形状；
(7)如果平均形状和参考形状之间的Procrustes距离高于设定的阈值，将平均形状设置为新的参考形状，然后继续执行步骤(2)。
以上技术的有益效果为：在构建头部可变形图谱时，头部的尺寸差异是人群之间的典型解剖形态差异，本发明采用部分Procrustes分析操作和GPA过程相结合，既保留例了个体的尺寸信息，只做空间位置和方向的归一化，又经过GPA处理之后，实现所有个体的精准归一化结果。
实施例10：
在实施例9的基础上，所述的一种基于CT影像的头部可变形统计图谱构建方法，基于主成分分析方法和所述整体配准结果，提取所述归一化形状网格中的统计形变分量，基于所述统计形变分量构建出对应的可变形图谱，包括：
计算出所有归一化形状网格的平均形状网格和协方差矩阵；
对所述协方差矩阵进行特征值分解，获得对应的由特征向量构成的正交矩阵和由特征值构成的对角矩阵；
将所述对角矩阵中包含的特征值从大到小排序获得对应的特征值序列；
将所述特征值序列中前n个特征值当作对应所述可变形图谱的数据集，将所述数据集中包含的每个第一特征值在所述正交矩阵中对应的第一特征向量当作对应归一化形状网格的统计形变分量；
基于所有归一化形状网格对应的权重变形参数和对应的统计形变分量的乘积与所述平均形状网格的线性组合，构建出对应的可变形图谱。
该实施例中，计算出所有归一化形状网格的平均形状网格和协方差矩阵，包括：
计算每个个体头部形状网格的平均形状网格：
Figure BDA0003638485180000311
式中，
Figure BDA0003638485180000312
为所有头部形状网格的平均形状向量，i为当前计算的头部形状网格，n为头部形状网格总个数，Surgi为第i个头部形状网格；
计算每个头部形状网格的协方差矩阵：
Figure BDA0003638485180000313
式中，C为当前计算的头部形状网格的协方差矩阵，i为当前计算的头部形状网格，n为头部形状网格总个数，Surgi为第i个头部形状网格，
Figure BDA0003638485180000314
为所有头部形状网格的平均形状向量，
Figure BDA0003638485180000315
为
Figure BDA0003638485180000316
的转置矩阵。
该实施例中，对所述协方差矩阵进行特征值分解，获得对应的由特征向量构成的正交矩阵和由特征值构成的对角矩阵，包括：
对协方差矩阵C进行特征值分解(eigen decomposition)：
C＝Q∑QT； (13)
其中，Q＝{φ,φ2,···，φn}是协方差矩阵C的特征向量组成的正交矩阵，QT为Q的转置矩阵，∑＝diag(σ1,σ2,···,σn)是协方差矩阵C的特征值对角矩阵；
特征值满足从大到小的顺序排列，即σ1≥σ2≥···≥σn，选择前面最大的n′(n′≤n)个特征值作为对整个数据集的表示,得到的特征向量*φ1,φ2,···，φn′}对应于形状网格的形状变化模式，而特征值σi(i∈[1,n])则是每种变化模式的对应方差，形状变化模式按其相应方差值大小排序(σ1≥σ2≥···≥σn)，变形模式1对应最大方差值，变形模式2对应第二大方差值，依此类推；因此，变形模式i在所有模式中所占的比例可以通过其方差值与所有方差值之比求得，即：
Figure BDA0003638485180000321
该实施例中，将所述对角矩阵中包含的特征值从大到小排序获得对应的特征值序列即为：特征值满足从大到小的顺序排列，即σ1≥σ2≥···≥σn。
该实施例中，数据集即为特征值序列中前n′(n′≤n)个特征值。
该实施例中，第一特征值即为数据集中包含的特征值。
该实施例中，第一特征向量即为第一特征值在正交矩阵中对应的特征向量。
该实施例中，基于所有归一化形状网格对应的权重变形参数和对应的统计形变分量的乘积与所述平均形状网格的线性组合，构建出对应的可变形图谱，包括：
该实施例中，SSM(统计模型)通常用平均形状加上不同形状变化模式的线性组合表示：
Figure BDA0003638485180000322
其中SSSM是SSM生成的形状实例，是由包含N个网格顶点三维坐标的形状向量{SSSM1,SSSM2,···,SSSMN}表示，
Figure BDA0003638485180000323
是所有形状网格的平均形状向量，ai{＝1,···,n}是形状变化模式权重的变形参数，不同的ai值对应SSM的不同形变实例，当ai值连续变化时，SSM的形状会相应的连续变化，从而引起头部解剖形状的实时变形。
以上技术的有益效果为：将每个训练个体对应的统计形变分量与所有训练个体形状网格的平均形状向量进行线性组合，构造出了点分布类型的头部可变形统计图谱。
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (10)

1.一种基于CT影像的头部可变形统计图谱构建方法，其特征在于，包括：
S1：基于最大熵阈值分割法对每个训练样本对应的头部CT影像进行分割并转化，获得每个训练样本中不同解剖结构对应的目标网格；
S2：获取头部参考模板和对应的3D多边形网格；
S3：基于解剖标志点约束的薄板样条鲁棒点匹配算法，采用轮廓匹配的分层式配准策略，将所述3D多边形网格与每个训练样本中不同解剖结构对应的目标网格进行配准，获得每个训练样本对应的个性化解剖结构表示；
S4：基于所述个性化解剖结构表示和统计形状建模方法构建出头部可变形统计图谱。
2.根据权利要求1所述的一种基于CT影像的头部可变形统计图谱构建方法，其特征在于，S1：基于最大熵阈值分割法对每个训练样本对应的头部CT影像进行分割并转化，获得每个训练样本中不同解剖结构对应的目标网格，包括：
S101：在海量头部CT影像中筛选出无症状个体头部CT影像作为训练样本；
S102：基于最大熵阈值分割法对所述训练样本对应的头部CT影像进行分割，获得对应的可分割解剖结构和未分割解剖结构；
S103：基于等值面提取方法将所述可分割解剖结构和所述未分割解剖结构转化为对应的三角面片网格，并将对应三角面片网格作为对应解剖结构的目标网格。
3.根据权利要求2所述的一种基于CT影像的头部可变形统计图谱构建方法，其特征在于，S2：获取头部参考模板和对应的3D多边形网格，包括：
S201：获取通用人体数字模型的头部模型作为所述头部参考模板；
S202：基于人体MR影像提取出所述头部参考模板的解剖结构信息；
S203：基于所述解剖结构信息构建出对应的3D多边形网格。
4.根据权利要求3所述的一种基于CT影像的头部可变形统计图谱构建方法，其特征在于，S3：基于解剖标志点约束的薄板样条鲁棒点匹配算法，采用轮廓匹配的分层式配准策略，将所述3D多边形网格与每个训练样本中不同解剖结构对应的目标网格进行配准，获得每个训练样本对应的个性化解剖结构表示，包括：
基于解剖标志点约束的薄板样条鲁棒点匹配算法，采用轮廓匹配的分层式配准策略，将所述3D多边形网格中包含的第一解剖结构与对应训练样本中包含的第二解剖结构对应的目标网格进行配准，获得对应训练样本的整体配准结果；
基于预设拓扑结构网格表示所述整体配准结果，获得每个训练样本对应的个性化解剖结构表示。
5.根据权利要求4所述的一种基于CT影像的头部可变形统计图谱构建方法，其特征在于，基于解剖标志点约束的薄板样条鲁棒点匹配算法，采用轮廓匹配的分层式配准策略，将所述3D多边形网格中包含的第一解剖结构与对应训练样本中包含的第二解剖结构对应的目标网格进行配准，获得对应训练样本的整体配准结果，包括：
基于所述解剖标志点约束的薄板样条鲁棒点匹配算法，将所述3D多边形网格中包含的第一可分割解剖结构与所述训练样本中包含的对应第二可分割解剖结构进行点云配准，获得对应的第一配准结果；
基于所述第一配准结果确定出所述头部参考模板至对应训练样本的空间形变场；
确定出所述3D多边形网格中包含的所有第一解剖结构之间的第一相对位置信息和所述训练样本中包含的所有第二解剖结构之间的第二相对位置信息；
基于所述训练样本对应的空间形变场和对应的第二相对位置信息以及所述第一相对位置信息，将所述3D多边形网格中包含的第一未分割参考解剖结构与对应训练样本中包含的第二未分割解剖结构进行二次配准，获得对应训练样本的第二配准结果；
综合所述训练样本对应的第一配准结果和第二配准结果，获得对应训练样本的整体配准结果。
6.根据权利要求5所述的一种基于CT影像的头部可变形统计图谱构建方法，其特征在于，基于所述解剖标志点约束的薄板样条鲁棒点匹配算法，将所述3D多边形网格中包含的第一可分割解剖结构与所述训练样本中包含的对应第二可分割解剖结构进行点云配准，获得对应的第一配准结果，包括：
人工标定出所述第一可分割解剖结构中包含的第一关键解剖标志点与所述训练样本中包含的对应第二可分割解剖结构中包含的第二关键解剖标志点；
确定出所述第一可分割解剖结构的第一点云表示和对应第二可分割解剖结构的第二点云表示；
基于所述第一点云表示和初始仿射变换矩阵以及初始非线性变形系数矩阵，确定出对应的初始非刚性变换函数；
基于所述初始非刚性变换函数确定出所述第一点云表示对应的初始非刚性变换映射；
基于预设的温度系数和所述第二点云表示以及所述初始非刚性变换映射确定出对应的权重系数，基于所述权重系数和所述第二点云表示确定出对应的权重系数矩阵；
基于所述第一关键解剖标志点、所述第二关键解剖标志点、所述第一点云表示、所述第二点云表示、所述权重系数矩阵和预设的正则化参数构建出对应的代价函数，将所述代价函数最小化，确定出新的非刚性变换函数；
基于最新确定的非刚性变换函数和所述第一点云表示确定出所述第二可分割解剖结构对应的与所述第一点云表示配准的第三点云表示，判断所述第二点云表示中每个点和所述第三点云表示中对应点之间的距离是否都小于预设阈值，若是，则将所述第三点云表示作为所述3D多边形网格中包含的第一可分割解剖结构对应的第一配准结果；
否则，基于预设减小梯度设置新的温度系数和新的正则化参数，基于新的温度系数和最新确定的非刚性变换函数确定出新的权重系数矩阵，基于新的权重系数矩阵和新的正则化参数构建出新的代价函数，将新的代价函数最小化，确定出新的非刚性变换函数，基于最新确定的非刚性变换函数确定出所述第二解剖结构对应的与所述第一点云表示配准的第四点云表示，直至所述第二点云表示中每个点和所述第四点云表示中对应点之间的距离都小于预设阈值时，则将所述第四点云表示作为所述3D多边形网格中包含的第一可分割解剖结构对应的第一配准结果。
7.根据权利要求6所述的一种基于CT影像的头部可变形统计图谱构建方法，其特征在于，基于预设拓扑结构网格表示所述整体配准结果，获得每个训练样本对应的个性化解剖结构表示，包括：
基于所述整体配准结果获得所述3D多边形网格中每个顶点至对应训练样本中对应顶点的映射关系；
基于预设拓扑结构网格表示所述整体配准结果，获得每个训练样本对应的个性化解剖结构表示。
8.根据权利要求7所述的一种基于CT影像的头部可变形统计图谱构建方法，其特征在于，S4：基于所述个性化解剖结构表示和统计形状建模方法构建出头部可变形统计图谱，包括：
基于广义普氏分析对所述个性化解剖结构中所述训练样本对应的整体配准结果中的形状网格进行空间方位归一化，获得所述训练样本对应的归一化形状网格；
基于主成分分析方法和所述整体配准结果，提取所述归一化形状网格中的统计形变分量，基于所述统计形变分量构建出对应的可变形图谱。
9.根据权利要求8所述的一种基于CT影像的头部可变形统计图谱构建方法，其特征在于，基于广义普氏分析对所述训练样本对应的整体配准结果中的形状网格进行空间方位归一化，获得所述训练样本对应的归一化形状网格，包括：
基于所述形状网格中包含的所有点在参考坐标系下对应的第一坐标值，计算出所述形状网格的形状中心对应的第二坐标值；
基于所述形状中心对应的第二坐标值对所述形状网格中包含的所有点进行去中心化处理，获得对应的原点对齐形状网格；
基于所述形状网格中包含的所有点在参考坐标系下对应的第一坐标值和所述形状中心对应的第二坐标值计算出对应训练样本的形状尺寸测度；
将所述原点对齐形状网格中每个点对应的第三坐标值除以对应的形状尺寸测度，获得对应的尺度归一化形状网络；
在所有训练样本对应的尺度归一化形状网络中任选一个尺度归一化形状网络作为标准形状；
将所有训练样本对应的尺度归一化形状网络中除所述标准形状以外剩余的待旋转形状与所述标准形状进行对齐，获得新的训练样本；
计算出新的训练样本对应的平均形状，计算出所述平均形状和所述参考形状之间的普氏距离平方值；
当所述普氏距离平方值大于设定阈值时，则将平均形状设置为新的参考形状，将所有训练样本对应的尺度归一化形状网络中除新的参考形状以外剩余的待旋转形状与新的参考形状进行对齐，直至最新获得的平均形状与最新的参考形状之间的普氏距离不大于设定阈值时，则将最新获得的训练样本作为对应的归一化形状网格。
10.根据权利要求9所述的一种基于CT影像的头部可变形统计图谱构建方法，其特征在于，基于主成分分析方法和所述整体配准结果，提取所述归一化形状网格中的统计形变分量，基于所述统计形变分量构建出对应的可变形图谱，包括：
计算出所有归一化形状网格的平均形状网格和协方差矩阵；
对所述协方差矩阵进行特征值分解，获得对应的由特征向量构成的正交矩阵和由特征值构成的对角矩阵；
将所述对角矩阵中包含的特征值从大到小排序获得对应的特征值序列；
将所述特征值序列中前n个特征值当作对应所述可变形图谱的数据集，将所述数据集中包含的每个第一特征值在所述正交矩阵中对应的第一特征向量当作对应归一化形状网格的统计形变分量；
基于所有归一化形状网格对应的权重变形参数和对应的统计形变分量的乘积与所述平均形状网格的线性组合，构建出对应的可变形图谱。
CN202210512328.6A 2022-05-11 2022-05-11 一种基于ct影像的头部可变形统计图谱构建方法 Active CN114913149B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210512328.6A CN114913149B (zh) 2022-05-11 2022-05-11 一种基于ct影像的头部可变形统计图谱构建方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210512328.6A CN114913149B (zh) 2022-05-11 2022-05-11 一种基于ct影像的头部可变形统计图谱构建方法

Publications (2)

Publication Number Publication Date
CN114913149A true CN114913149A (zh) 2022-08-16
CN114913149B CN114913149B (zh) 2023-03-10

Family

ID=82766994

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210512328.6A Active CN114913149B (zh) 2022-05-11 2022-05-11 一种基于ct影像的头部可变形统计图谱构建方法

Country Status (1)

Country Link
CN (1) CN114913149B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150213646A1 (en) * 2014-01-28 2015-07-30 Siemens Aktiengesellschaft Method and System for Constructing Personalized Avatars Using a Parameterized Deformable Mesh
CN110335358A (zh) * 2019-06-18 2019-10-15 大连理工大学 可变形数字人解剖学模型的个性化变形方法
CN114359309A (zh) * 2022-01-12 2022-04-15 大连理工大学 基于标定点检测和形状灰度模型匹配的医学图像分割方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150213646A1 (en) * 2014-01-28 2015-07-30 Siemens Aktiengesellschaft Method and System for Constructing Personalized Avatars Using a Parameterized Deformable Mesh
CN110335358A (zh) * 2019-06-18 2019-10-15 大连理工大学 可变形数字人解剖学模型的个性化变形方法
CN114359309A (zh) * 2022-01-12 2022-04-15 大连理工大学 基于标定点检测和形状灰度模型匹配的医学图像分割方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
HONGKAI WANG ET AL.: "Deformable torso phantoms of Chinese adults for personalized anatomy modelling", 《ANATOMICAL SOCIETY》 *
ZHAOFENG CHEN ET AL.: "Deformable Head Atlas of Chinese Adults Incorporating Inter-Subject Anatomical Variations", 《DIGITAL OBJECT IDENTIFIER》 *

Also Published As

Publication number Publication date
CN114913149B (zh) 2023-03-10

Similar Documents

Publication Publication Date Title
Kelemen et al. Three-dimensional model-based segmentation of brain MRI
Montagnat et al. Volumetric medical images segmentation using shape constrained deformable models
Thompson et al. High-resolution random mesh algorithms for creating a probabilistic 3D surface atlas of the human brain
Kelemen et al. Elastic model-based segmentation of 3-D neuroradiological data sets
Yang et al. Automatic segmentation of parotids from CT scans using multiple atlases
Deng et al. A novel skull registration based on global and local deformations for craniofacial reconstruction
CN109118455B (zh) 一种基于现代人软组织分布的古人类头骨颅面交互复原方法
Desvignes et al. 3D semi-landmarks based statistical face reconstruction
CN104867104A (zh) 基于xct图像非刚度配准的目标鼠解剖结构图谱获取方法
CN115830016B (zh) 医学图像配准模型训练方法及设备
CN109978998B (zh) 一种基于面部软组织和形状空间的古人类颅面复原方法
CN115116586A (zh) 一种基于联合配准的可变形统计图谱构建方法
Berar et al. 3D statistical facial reconstruction
CN111127488B (zh) 一种基于统计形状模型自动构建患者解剖结构模型的方法
WO2002003304A2 (en) Predicting changes in characteristics of an object
CN114913149B (zh) 一种基于ct影像的头部可变形统计图谱构建方法
Xie et al. Tissue feature-based and segmented deformable image registration for improved modeling of shear movement of lungs
Kim et al. Organ shape modeling based on the laplacian deformation framework for surface-based morphometry studies
Xue et al. Lung respiratory motion estimation based on fast Kalman filtering and 4D CT image registration
CN110322491B (zh) 一种可变形小鼠全身图谱与小鼠影像配准的算法
CN114359309A (zh) 基于标定点检测和形状灰度模型匹配的医学图像分割方法
Lötjönen et al. Four-chamber 3-D statistical shape model from cardiac short-axis and long-axis MR images
Chen et al. Inter-subject shape correspondence computation from medical images without organ segmentation
Wu et al. Multi-organ Statistical Shape Model Building Using a Non-rigid ICP Based Surface Registration
Dalal et al. 3D open-surface shape correspondence for statistical shape modeling: Identifying topologically consistent landmarks

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