CN112053330B - 基于pca和tssm模型的膈肌预测系统及方法 - Google Patents

基于pca和tssm模型的膈肌预测系统及方法 Download PDF

Info

Publication number
CN112053330B
CN112053330B CN202010884814.1A CN202010884814A CN112053330B CN 112053330 B CN112053330 B CN 112053330B CN 202010884814 A CN202010884814 A CN 202010884814A CN 112053330 B CN112053330 B CN 112053330B
Authority
CN
China
Prior art keywords
diaphragm
dimensional
displacement
chest
patient
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
CN202010884814.1A
Other languages
English (en)
Other versions
CN112053330A (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.)
Xidian University
Original Assignee
Xidian University
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 Xidian University filed Critical Xidian University
Priority to CN202010884814.1A priority Critical patent/CN112053330B/zh
Publication of CN112053330A publication Critical patent/CN112053330A/zh
Application granted granted Critical
Publication of CN112053330B publication Critical patent/CN112053330B/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
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G06F18/2135—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on approximation criteria, e.g. principal component analysis
    • 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/60—Analysis of geometric attributes
    • G06T7/62—Analysis of geometric attributes of area, perimeter, diameter or volume
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/70—Determining position or orientation of objects or cameras
    • G06T7/73—Determining position or orientation of objects or cameras 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/20—Special algorithmic details
    • G06T2207/20081—Training; Learning
    • 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/30061—Lung

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • General Health & Medical Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Computation (AREA)
  • Quality & Reliability (AREA)
  • Radiology & Medical Imaging (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Medical Informatics (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Health & Medical Sciences (AREA)
  • Algebra (AREA)
  • Geometry (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

本发明公开了一种基于PCA和TSSM模型的膈肌预测系统及方法，其实现的思路是：将4DCT扫描获取的呼吸时相数据重建成每个相位的人体三维图像；从人体三维图像中分割出膈肌图像，计算膈肌位移；从人体三维图像中分割出胸腹表面图像，计算胸腹表面位移；对膈肌位移和胸腹表面位移分别进行主成分分析，使其降维并映射到各自的d维子空间中；利用岭回归优化计算出从d维子空间的胸腹表面位移映射到d维子空间的膈肌位移的变换矩阵。本发明具有对患者的伤害更少，选取的膈肌质心更加准确，抗干扰性更好的优点。

Description

基于PCA和TSSM模型的膈肌预测系统及方法
技术领域
本发明属于图像处理技术领域，更进一步涉及医学影像辅助干预技术领域中的一种基于主成分分析PCA(Principal Components Analysis)和分步子空间映射TSSM(Two-Step Subspace Mapping)模型的膈肌预测系统及方法。本发明在不使用标记物的情况下，通过体外胸腹部表面的位移预测人体膈肌的运动。
背景技术
动态目标跟踪法是目前最为有效的呼气管理技术，直接方式是植入体内的标志物进行成像，进而进行目标跟踪。但实验对象需要进行昂贵的、创伤性手术。植入标志物的替代方法，是使用呼吸信息对内部解剖结构的运动进行建模。通过肺活量测定和实时位置跟踪系统，研究肺肿瘤运动与呼吸运动的相关性。相比于其它解剖结构，膈肌与呼吸运动具有更强的相关性。另一方面，人体内部结构的运动情况，可以通过机体外部来进行观察和预测。
在Malinowski等人在其发表的论文"Mitigating Errors in ExternalRespiratory Surrogate-Based Models of Tumor Position"(J International Journalof Radiation Oncology,Biology,Physics,"vol.82,no.5,2012)中公开公布了一种减轻外呼吸替代肿瘤位置模型的误差的方法。该方法的主要实现步骤是：在患者身上贴附了多个标志物，并建立了标志物-肿瘤的位移关系模型。该方法通过评估肿瘤位置预测模型的潜在误差来源，最终选择了偏最小二乘回归(PLS)进行模型拟合。该方法存在的不足之处是，模型的准确性与患者数量显著相关，采用外标志物建立的肿瘤位置预测模型的准确性受患者间差异、肿瘤-替代物相关性、训练数据选择的影响。
在南方医科大学申请的专利“一种胸部数字合成X射线体层成像中呼吸运动分析方法”(公开号：CN 109146842 A，申请号：CN 201810745977.4，申请日：2018年07月09日)中公开了一种胸部数字合成X射线体层成像中呼吸运动分析方法，其制备方法包括以下步骤：S1，通过胸部数字合成体层成像，获取CTS投影数据；S2,从S1中获取的CTS投影数据提取膈肌的运动轨迹；S3,从S2获取的膈肌运动进行进行正弦二次模型拟合。该方法存在的不足之处是：部分患者在整个数据采集过程中可能无法屏住呼吸，在采集结束时开始呼吸。且质心仅以最大振幅的点来表示，忽略了膈肌在左右方向上的位移，导致振幅最大的点正下方振幅最小时的点并不是同一个点从而使得简化对象不准确的问题。
发明内容
本发明的目的是针对上述已有技术的不足，提出了一种基于PCA和TSSM模型的膈肌预测系统及方法，用于解决PLS模型外标志物建立的肿瘤位置预测模型的准确性受患者个体差异、测量精度、肿瘤-替代物相关性、训练数据选择和正弦二次模型对患者数据采集时可能无法屏住呼吸而产生数据采集不全的问题。
实现本发明目的的思路是：对于PLS模型和正弦二次模型都存在的初始采集的数据就会受到患者差异影响的问题，本发明采用了主成分分析法对患者的数据进行降维处理，保留了不同患者之间的共性，极大的消除个体差异性；PLS模型中标志物的移动会导致模型误差随时间越来越大。本发明采用了ICP算法，通过对不同相位的胸腹表面图像的像素点进行配准计算出胸腹表面位移，而人体本身的胸腹表面的相对位置是相对固定的，使得模型不会随时间变差，也减少了患者植入标志物的痛苦。
本发明的系统包括图像分割模块，膈肌位移计算模块、胸腹表面位移模块、图像输出模块、PCA模块、模型拟合模块；其中，
所述膈肌图像分割模块，用于选取所选患者的人体三维图像，将所选相位的人体三维图像以像素值差值法处理为仅有右肺下膈肌的三维图像；
所述胸腹表面图像分割模块，用于选取所选患者的所选相位的人体三维图像，将该图像以掩膜法处理为仅有胸腹表面的三维图像；
所述膈肌位移计算模块，将膈肌三维图像的全部膈肌像素点的空间三维坐标值取均值，将该均值作为右肺下膈肌的质心坐标值；将该坐标值与所选患者的右肺下膈肌的初始质心坐标值相减，得到所选相位的膈肌位移；
所述胸腹表面位移计算模块，用于使用ICP算法将所选患者的所选相位的胸腹部三维图像的空间三维坐标与呼吸周期的0％相位的胸腹表面图的空间三维坐标进行配准，将配准后得到的平移向量作为所选相位的胸腹表面位移；
所述PCA模块，用于将每个患者全部相位的膈肌位移连接为一个10个x轴上的位移值+10个y轴上的位移值+10个z轴上的位移的30维的向量，并进行归一化；用于将每个患者的归一化后的30维的向量，进行主成分分析法处理，降维并映射到d维子空间Rv中；用于将每个患者每个相位的胸腹表面位移归一化为一个10个x轴上的位移值+10个y轴上的位移值+10个z轴上的位移的30维的向量，并进行归一化；用于将每个患者归一化后的30维向量，进行主成分分析法处理，降维并映射到d维子空间Ru中；用于将经过TSSM模型计算后得到的膈肌位移
Figure GDA0004188921550000031
按照PCA特征子空间的系数和基，即/>
Figure GDA0004188921550000032
计算得到归一化的膈肌位移数据/>
Figure GDA0004188921550000033
按照膈肌位移的归一化公式，将归一化的膈肌位移数据/>
Figure GDA0004188921550000034
逆归一化，得到最终的膈肌位移vout；
所述模型拟合模块，用于利用线性岭回归优化法将映射到Rv中的膈肌位移数据与映射到Ru中的胸腹部表面位移数据连接起来，得到变换矩阵β，求出β的最优化表达式βopt，
Figure GDA0004188921550000035
即为TSSM模型。
所述图像输出模块在计算机上显示出横坐标为相位，纵坐标为膈肌位移的曲线图。
本发明方法的具体步骤包括如下：
(1)呼吸时相数据采集：
(1a)选取患有肺部癌症或胰腺癌的至少20个患者，利用4DCT扫描仪获取每个患者连续的至少5个呼吸周期的时相数据，计算每个患者的平均呼吸周期，从5个呼吸周期中选取最接近平均呼吸周期的一个时相数据，作为该患者的呼吸周期；
(1b)用与4DCT扫描仪连接的计算机打开每个患者的平均呼吸周期时相数据，将该患者呼吸周期从0％开始，每隔10％处选取一次时相数据重建成10个相位的人体三维图像；
(1c)从20个患者中随机选取16个患者，将所选16个患者的人体三维图像包含的全部时相数据集组成初始训练数据，剩余4个患者的人体三维图像作为测试数据；
(2)计算每个患者的4DCT图像的膈肌位移和胸腹表面位移：
(2a)从16个患者中选取一个未选过的患者：
(2b)膈肌图像分割模块选取所选患者呼吸周期0％相位的患者体三维图像，利用像素值差值法，将所选相位的人体三维图像处理为仅有右肺下膈肌的三维图像；
(2c)胸腹表面图像分割模块选取所选患者的呼吸周期的0％相位的人体三维图像，将该图像以掩膜法处理为仅有胸腹表面的三维图像；
(2d)选取所选患者的一个未选过相位的人体三维图像，采用与步骤(2b)相同的处理方法，对所选相位的人体三维图像进行处理，得到所选相位的人体三维图像中仅有右肺下膈肌的膈肌三维图像，膈肌位移计算模块将该膈肌三维图像的全部膈肌像素点的空间三维坐标取均值，将该均值作为右肺下膈肌的质心坐标；将质心坐标值与所选患者的右肺下膈肌的初始质心坐标相减，得到所选相位的膈肌位移；
(2e)采用与步骤(2c)中相同的处理方法，对所选相位的人体三维图像进行处理，得到所选相位的人体三维图像中仅有的胸腹表面的胸腹部三维图像，胸腹表面位移计算模块用ICP算法将所选患者的所选相位的胸腹部三维图像的空间三维坐标值与呼吸周期的0％相位的胸腹表面图的空间三维坐标值进行配准，将配准后得到的平移向量作为所选相位的胸腹表面位移；
(2f)判断是否选完所选患者的所有相位的人体三维图像，若是，则执行步骤(2g)，否则，执行步骤(2d)；
(2g)判断是否选完16个患者，若是，则执行步骤(3)，否则，执行步骤(2a)；
(3)将膈肌位移数据和胸腹表面位移数据映射到各自的子空间中：
(3a)模块将每个患者全部相位的膈肌位移连接为一个10个x轴上的位移值+10个y轴上的位移值+10个z轴上的位移的30维的向量，并进行归一化；
(3b)PCA模块将每个患者的归一化后的30维的向量，进行主成分分析法处理，降维并映射到d维子空间Rv中；
(3c)PCA模块将每个患者每个相位的胸腹表面位移归一化为一个10个x轴上的位移值+10个y轴上的位移值+10个z轴上的位移的30维的向量，并进行归一化；
(3d)PCA模块将每个患者归一化后的30维向量，进行主成分分析法处理，降维并映射到d维子空间Ru中；
(4)模型拟合模块利用线性岭回归优化法将映射到Rv中的膈肌位移数据与映射到Ru中的胸腹部表面位移数据连接起来，得到变换矩阵β，求出β的最优化表达式βopt，
Figure GDA0004188921550000051
即为TSSM模型。
(5)实时输出测试患者的相位-膈肌位移的曲线图：
(5a)用4DCT扫描仪获取测试患者的呼吸时相数据；
(5b)测试患者的呼吸时相数据依次经过：胸腹表面图像分割模块、胸腹表面位移计算模块、PCA模块、TSSM模型、PCA模块的逆归一化处理、最后输出测试患者的膈肌位移；
(5c)图像输出模块在计算机上显示出横坐标为膈肌位移相位，纵坐标为膈肌位移的曲线图。
(6)进一步将TSSM推广到基于核映射的非线性子空间：
将胸腹表面位移数据用经过核函数变换后的线性组合来表示，基于核函数的非线性岭回归优化变换kTSSM模型。
本发明与现有技术相比具有以下优点：
第一，由于本发明系统中的胸腹表面图像分割模块用于从每个相位的人体三维中分割出胸腹表面图像，通过计算胸腹表面位移预测出膈肌位移，克服了现有技术在偏最小二乘回归模型中，通过计算外标志物位移计算出肿瘤位移的技术中需要在患者胸腹表面植入标志物从而对患者造成创伤的问题，使得本发明具有对患者的伤害更少的优点。
第二，由于本发明中采用了膈肌图像分割，并以所有膈肌像素点的平均值来表示膈肌质心的方法，克服了现有技术在正弦二次模型中，仅以膈肌上呼吸振幅最大的点表示膈肌质心，但膈肌的运动在左右方向上有一定的偏移，导致在振幅最高处和振幅最低处的点实际上并不是膈肌上同一个点的缺点，使得本发明具有选取的膈肌质心更加准确的优点。
第三，由于本发明方法中采用了主成分分析法，将患者的膈肌位移和胸腹表面位移降维并映射到仅包涵数据本质特征的子空间中，克服了现有技术中采集呼吸时相数据时就易受患者的个体差异、噪声的影响的问题，使得本发明具有更抗干扰的优点。
附图说明
图1是本发明系统的模块示意图；
图2是本发明方法的流程图；
图3是本发明仿真图；
图4是本发明的性能评价箱线图。
具体实施方式
下面结合附图对本发明做进一步的详细描述。
参照附图1，对本发明的系统做进一步的详细描述。
所述膈肌图像分割模块，用于选取所选患者呼吸周期0％相位的患者体三维图像，将所选相位的人体三维图像以像素值差值法处理为仅有右肺下膈肌的三维图像。
所述胸腹表面图像分割模块，用于选取所选患者的所选相位的人体三维图像，将该图像以掩膜法处理为仅有胸腹表面的三维图像。
所述膈肌位移计算模块，将膈肌三维图像的全部膈肌像素点的空间三维坐标取均值，将该均值作为右肺下膈肌的质心坐标；将该坐标与所选患者的右肺下膈肌的初始质心坐标相减，得到所选相位的膈肌位移。
所述胸腹表面位移计算模块，用于用ICP算法将所选患者的所选相位的胸腹部三维图像的空间三维坐标与呼吸周期的0％相位的胸腹表面图的空间三维坐标进行配准，将配准后得到的平移向量作为所选相位的胸腹表面位移。
所述PCA模块，用于将每个患者的10个相位归一化后的向量，进行主成分分析法处理，得到降维并映射到d维子空间Rv中的；用于将每个患者的10个相位，归一化后的向量，进行主成分分析法处理，得到降维并映射到d维子空间Rv中的；用于将经过TSSM模型计算得到的膈肌位移
Figure GDA0004188921550000071
根据PCA特征子空间的系数和基来计算膈肌位移的预测数据，即/>
Figure GDA0004188921550000072
得到归一化的膈肌位移数据/>
Figure GDA0004188921550000073
根据膈肌位移的归一化公式，将归一化的膈肌位移数据/>
Figure GDA0004188921550000074
逆归一化，得到最终的膈肌位移vout。
所述模型拟合模块，用于利用线性岭回归优化法将映射到Rv中的膈肌位移数据与映射到Ru中的胸腹部表面位移数据连接起来，得到变换矩阵β，求出β的最优化表达式βopt，
Figure GDA0004188921550000075
即为TSSM模型。
所述图像输出模块在计算机上显示出横坐标为相位，纵坐标为膈肌位移的曲线图。
参照附图2，对本发明的方法做进一步的详细描述。
步骤1，呼吸时相数据采集。
选取患有肺部癌或胰腺癌或肝癌的至少20个患者，利用4DCT扫描仪获取每个患者连续的至少5个呼吸周期的时相数据，计算每个患者的平均呼吸周期，从5个呼吸周期中选取最接近平均呼吸周期的一个时相数据，作为该患者的呼吸周期。
用与4DCT扫描仪连接的计算机打开每个患者的平均呼吸周期时相数据，将该患者呼吸周期从0％开始，每隔10％处选取一次时相数据重建成10个相位的人体三维图像。
从20个患者中随机选取16个患者，将所选16个患者的人体三维图像包含的全部时相数据集组成初始训练集；剩余4个患者的全部时相数据集组成测试集。
步骤2，计算出每个患者的4DCT图像的膈肌位移和胸腹表面位移：
(2.1)从16个患者中选取一个未选过的患者：
(2.2)膈肌图像分割模块选取所选患者呼吸周期0％相位的患者体三维图像，将所选相位的人体三维图像以像素值差值法处理为仅有右肺下膈肌的三维图像，对右肺下膈肌的三维图像中全部膈肌像素点的空间三维坐标取均值，将该均值作为右肺下膈肌的初始质心坐标。
所述的像素值差值法的步骤如下：
(2.2.1)将所选相位的三维人体图像中的最大目标区域作为身体区域；通过形态学方法中的开运算，去除每个相位的三维人体图像中背景中不属于人体，但像素不为0的孤立区域；利用形态学方法中的闭运算填充身体区域中的孔，得到分割后的身体区域。
(2.2.2)从分割后的身体区域中，分割出肺部区域；对三维人体图像，将背景像素值均置为0，分割后的身体区域像素值均置为1，肺部区域像素值均置为3；得到三维人体分割图像。
(2.2.3)用三维ROI掩膜对三维人体分割图像进行分割，得到从前向后看的肺部三维图像。
(2.2.4)对肺部三维图像，从左到右在矢状平面上计算肺部区域的各个切片的面积，生成横坐标为左右方向的三维ROI掩膜切面，纵坐标为肺部面积的曲线图；选择曲线两个峰之间的最低点，作为中位点，定位出该点横坐标的坐标值。
(2.2.5)在三维人体分割图像中，将中位点以左的人体和肺部的像素点全部置为0；再从下到上找到像素值差值为2的分界面，该边界面上的全部像素组成右肺下表面三维图像。
(2.2.6)对该右肺下表面，从上方，先通过形态学运算中的腐蚀去除掉一些属于肺部，但不属于膈肌部分的像素点，再在经过腐蚀后的右肺下表面上填充被过度腐蚀掉的膈肌像素点。
(2.3)胸腹表面图像分割模块选取所选患者的呼吸周期的0％相位的人体三维图像，将该图像以掩膜法处理为仅有胸腹表面的三维图像。
所述的掩膜法的步骤如下：
(2.3.1)对人体三维图像中背景区域的所有像素值均设置为0、身体区域中的所有像素值均设置为1；沿x轴负方向找出像素值差值为1的边界面，将该边界面称为初始胸腹部表面。
(2.3.2)在初始胸腹部表面上沿y轴正方向找到最宽的位置，最宽位置对应人体的肩膀，去除肩膀之上的区域，进行腐蚀操作，得到了每个相位的胸腹表面图。
(2.3.3)将所选患者的所有相位的胸腹表面图取交集，得到胸腹部掩膜，将该掩膜与初始胸腹部表面取交集，得到最终胸腹部表面。
(2.4)选取所选患者的一个未选过相位的人体三维图像，采用像素值差值法，对所选相位的人体三维图像进行处理，得到所选相位的人体三维图像中仅有右肺下膈肌的膈肌三维图像，膈肌位移计算模块将该膈肌三维图像的全部膈肌像素点的空间三维坐标取均值，将该均值作为右肺下膈肌的质心坐标。
(2.5)将所选患者所选相位的右肺下膈肌的质心坐标与所选患者的右肺下膈肌的初始质心坐标相减，得到所选相位的膈肌位移，其表达式如下：
Figure GDA0004188921550000091
其中，i＝1,2,……16表示训练集内的患者的序号，k＝1,2,……,10表示4DCT图像的相位的序号，
Figure GDA0004188921550000092
表示第i个患者的第k个相位的右肺下膈肌全部的像素点的质心，/>
Figure GDA0004188921550000093
表示第i个患者呼吸周期的0％相位的右肺下膈肌全部的像素点的质心；/>
Figure GDA0004188921550000094
表示第i个患者第k个相位右肺下膈肌全部的像素点在x轴上的坐标值的平均值，/>
Figure GDA0004188921550000095
表示第i个患者呼吸周期的0％相位的右肺下膈肌全部的像素点在x轴上坐标值的平均值；/>
Figure GDA0004188921550000096
表示第i个患者第k个相位右肺下膈肌全部的像素点在y轴上的坐标值的平均值，/>
Figure GDA0004188921550000097
表示第i个患者的呼吸周期0％相位右肺下膈肌全部的像素点在y轴上的坐标值的平均值，/>
Figure GDA0004188921550000098
表示第i个患者第k个相位右肺下膈肌全部的像素点在z轴上的坐标值的平均值，/>
Figure GDA0004188921550000099
表示第i个患者的呼吸周期0％相位右肺下膈肌全部的像素点在z轴上的坐标值的平均值；/>
Figure GDA00041889215500000910
表示第i个患者第k个相位右肺下膈肌质心的位移；/>
Figure GDA00041889215500000911
表示第i个患者第k个相位右肺下膈肌的质心在x轴上的位移，/>
Figure GDA00041889215500000912
表示第i个患者第k个相位右肺下膈肌的质心在y轴上的位移，/>
Figure GDA00041889215500000913
表示第i个患者第k个相位右肺下膈肌的质心在z轴上的位移。
为了便于表达，本发明后续将忽略平均符号和差值符号，使用vi(k)＝{oi(k),pi(k),qi(k)},k＝1,2,……,10表示第i个病人第k个相位的膈肌位移。
(2.6)采用掩膜法，对所选相位的人体三维图像进行处理，得到所选相位的人体三维图像中仅有的胸腹表面的胸腹部三维图像，胸腹表面位移计算模块用ICP算法将所选患者的所选相位的胸腹部三维图像的空间三维坐标与呼吸周期的0％相位的胸腹表面图的空间三维坐标进行配准，将配准后得到的平移向量作为所选相位的胸腹表面位移，其表达式如下：
Figure GDA0004188921550000101
其中，
Figure GDA0004188921550000102
表示第i个患者第k个相位的胸腹表面的位移，/>
Figure GDA0004188921550000103
表示第i个患者第k个相位的沿x轴方向上的位移，/>
Figure GDA0004188921550000104
表示第i个患者第k个相位的沿y轴方向上的位移，
Figure GDA0004188921550000105
表示第i个患者第k个相位的沿z轴方向上的位移。
为了便于表达，本发明后续将忽略平均符号和插值符号，使用ui(k)＝{ri(k),si(k),ti(k)},k＝1,2,……,10表示第i个病人第k个相位的胸腹部表面位移。
所述ICP算法是指：
(2.6.1)将呼吸周期0％相位的胸腹表面图的全部胸腹表面像素点的三维坐标作为源点云M，取点集mj∈M。
(2.6.2)将所选相位的胸腹表面图的全部胸腹表面像素点的三维坐标作为目标点云N，取点集nj∈N。
(2.6.3)计算映射矩阵H，使误差函数E(G,l)最小，映射矩阵H和误差函数E(G,l)如下：
Figure GDA0004188921550000106
Figure GDA0004188921550000111
Figure GDA0004188921550000112
其中，j表示点集mj和nj中每个点的序列号，J表示j的最大值，GJ×J表示J阶旋转矩阵，lJ×1表示J维平移向量，O1×J表示J维零向量，
Figure GDA0004188921550000113
表示使{·}取得最小值时的G和l的取值。
(2.6.4)对点集nj进行旋转和平移，计算点集nj'＝{nj'＝Gnj+l,nj∈N}。
(2.6.5)计算映射矩阵H'，使误差函数E(G',l')最小，映射矩阵H'和E(G',l')如下：
Figure GDA0004188921550000114
Figure GDA0004188921550000115
Figure GDA0004188921550000116
第6步，判断误差函数E(G',l')是否小于阈值或者大于预设的最大迭代次数，若是，则停止迭代更新，对所有平移向量求和，即为将所选相位的胸腹表面图的全部胸腹表面像素点与呼吸周期0％相位的胸腹表面图的全部胸腹表面像素点，进行配准后的平移向量，并执行(2.7)；否则，执行本步骤的(2.6.4)。
(2.7)判断是否选完所选患者的所有相位的人体三维图像，若是，执行本步骤的(2.8)，否则，执行本步骤的(2.4)。
(2.8)判断是否选完16个患者，若是，执行步骤3，否则，执行本步骤的(2.1)。
步骤3，将膈肌位移数据和胸腹表面位移数据映射到各自的子空间中。
PCA模块利用主成分分析法，将膈肌位移vi和胸腹表面位移ui，进行归一化，降维并映射到各自的d维子空间Rv，Ru中。
所述的主成分分析法是指：
(3.1)分别将每个患者的10个相位的膈肌位移值中每个相位的x轴、y轴、z轴的位移值，分别组成一个10维向量，再将这3个向量，按照x轴向量+y轴向量+z轴向量，组成一个30维的向量Vi；
Vi＝{oi(1),oi(2),...,oi(10),pi(1),pi(2),...,pi(10),qi(1),qi(2),...qi(10)}
采用与本步骤的(3.1)相同的处理方法，将每个患者的10个相位的胸腹表面位移值组成一个30维的向量Ui；
Ui＝{ri(1),ri(2),...,ri(10),si(1),si(2),...,si(10),ti(1),ti(2),...ti(10)}
(3.2)对向量Vi，Ui分别进行归一化处理，得到
Figure GDA0004188921550000121
和/>
Figure GDA0004188921550000122
归一化公式如下：
Figure GDA0004188921550000123
Figure GDA0004188921550000124
其中，
Figure GDA0004188921550000125
表示第i个患者膈肌位移归一化后的30维向量，i表示所选患者的序号，Vi(n)表示向量Vi内每一个元素，/>
Figure GDA0004188921550000126
表示第i个患者膈肌位移向量Vi(n)中30个元素中最小的元素，/>
Figure GDA0004188921550000127
表示第i个患者膈肌位移向量Vi(n)中30个元素中最大的元素；/>
Figure GDA0004188921550000128
表示第i个患者胸腹表面位移归一化后的30维向量，Ui(n)表示向量Ui内每一个元素，/>
Figure GDA0004188921550000129
第i个患者胸腹表面位移3向量Ui(n)30个元素中最小的元素，/>
Figure GDA00041889215500001210
表示第i个患者胸腹表面位移向量Ui(n)30个元素中最大的元素。
(3.3)按照下式，计算协方差矩阵C和B：
Figure GDA00041889215500001211
Figure GDA00041889215500001212
其中，
Figure GDA0004188921550000131
表示第i个患者膈肌位移归一化后的30维向量，/>
Figure GDA0004188921550000132
表示/>
Figure GDA0004188921550000133
的转置，/>
Figure GDA0004188921550000134
表示第i个患者胸腹表面位移归一化后的30维向量，/>
Figure GDA0004188921550000135
表示/>
Figure GDA0004188921550000136
的转置。
(3.4)按照下式，计算协方差矩阵B的特征值及对应的特征向量：
det(λI-B)＝0
(λI-B)e＝0
其中，det(·)表示行列式符号，λ表示协方差矩阵B的特征值，I表示和协方差矩阵B阶数相同的单位矩阵，e表示特征值λ对应的特征向量。
按照与本步骤的(3.4)相同的处理方式，计算协方差矩阵C的特征值及对应的特征向量。
(3.5)对协方差矩阵B的全部特征值，从大到小选择大于阈值ε的d个特征值，阈值ε在构造域的过程中，需要在数据完整性和数据冗余之间进行平衡。
按照与本步骤的(3.5)相同的处理方式，选择协方差矩阵C的d个特征值。
(3.6)根据d个特征值，得到d个对应的特征向量，该d个特征向量分别作为膈肌位移数据和胸腹部表面位移数据的基，表示为ωv和ωu，ωu,ωv∈R30×d，分别构建出d维子空间Rv和Ru。
(3.7)
Figure GDA0004188921550000137
和/>
Figure GDA0004188921550000138
是正交化的(/>
Figure GDA0004188921550000139
且/>
Figure GDA00041889215500001310
其中Id是d的单位矩阵)；分别通过运算/>
Figure GDA00041889215500001311
和/>
Figure GDA00041889215500001312
将膈肌位移数据和胸腹表面位移数据映射到各自的子空间中，其中，/>
Figure GDA00041889215500001313
是/>
Figure GDA00041889215500001314
在子空间Ru中的对应点，/>
Figure GDA00041889215500001315
是/>
Figure GDA00041889215500001316
在子空间Rv中的对应点。
步骤4，模型拟合模块利用线性岭回归优化法将映射到Rv中的膈肌位移数据与映射到Ru中的胸腹部表面位移数据连接起来，得到变换矩阵β，求出β的最优化表达式βopt，
Figure GDA00041889215500001317
即为TSSM模型。
所述的变换矩阵β的公式如下：
Figure GDA00041889215500001318
其中，
Figure GDA0004188921550000141
表示当{·}取最小值时的β的值，||·||2表示2范数，λ≥0是控制收缩量的规则化参数。
所述的最优表达式βopt操作是指：对β求导，并使等式等于0，则可以获得最优表达式，如下所示：
Figure GDA0004188921550000142
其中，opt表示最优，I表示单位矩阵。
TSSM模型表达式如下：
Figure GDA0004188921550000143
其中，
Figure GDA0004188921550000144
表示测试数据经过TSSM模型的输出结果，/>
Figure GDA0004188921550000145
表示测试患者的胸腹表面位移在d维子空间Ru的表示。
步骤5，实时输出测试患者的相位-膈肌位移的坐标图。
用4DCT扫描仪获取测试患者的呼吸时相数据。
将测试患者的呼吸时相数据依次经过：胸腹表面图像分割模块、胸腹表面位移计算模块、PCA模块、TSSM模型、PCA模块的逆归一化处理，最后输出测试患者的膈肌位移。
在计算机上显示出横坐标为相位，纵坐标为膈肌位移的曲线图。
步骤6，推导出基于核函数的非线性子空间的kTSSM模型。
用
Figure GDA0004188921550000146
替换掉TSSM中的/>
Figure GDA0004188921550000147
用K(·,·)替换掉两个高维子空间Wu中的两个向量的内积，得到非线性岭回归优化的变换矩阵β'opt和kTSSM模型。
Figure GDA0004188921550000148
其中，K(·,·)表示映射的核函数，
Figure GDA0004188921550000149
表示测试数据中d维子空间Σ中的胸腹表面位移，/>
Figure GDA00041889215500001410
表示训练数据中d维子空间Ru中的胸腹表面位移，φ表示非线性子空间Ru到高维子空间Wu的映射。
利用
Figure GDA00041889215500001411
替换掉非线性子空间Ru的胸腹表面位移，得到基于核函数的非线性岭回归优化变换表达式如下：
Figure GDA0004188921550000151
其中，β'表示非线性岭回归优化变换矩阵。
岭回归优化方程的最优解表示为：
Figure GDA0004188921550000152
kTSSM模型表达式如下：
Figure GDA0004188921550000153
Figure GDA0004188921550000154
其中，
Figure GDA0004188921550000155
表示使用kTSSM模型计算后得到的d维子空间Rv中的膈肌位移，/>
Figure GDA0004188921550000156
表示训练数据中的d维子空间Rv中的膈肌位移。
本发明的效果可通过以下仿真进一步说明。
1.仿真实验条件：
本发明的仿真实验的硬件平台为：处理器为Inter Xeon 3.6GHz，32G内存。
本发明的仿真实验的软件配置为Windows 10操作系统和Matlab R2019b。
2.仿真内容：
本发明的仿真是用本发明的系统和本发明的方法，采用了来自20位实验对象的4DCT数据集，每个实验对象包括10个相位(具体为平均呼吸周期的0％、10％、…、90％)。随机选出16个患者的呼吸时相数据作为训练数据，剩余4名患者的呼吸时相数据作为测试数据。
本发明仿真实验是采用本发明的TSSM模型和2个kTSSM模型和一个现有技术(偏最小二乘回归模型)分别对输入的训练患者的呼吸时相数据和测试患者的呼吸时相数据进行计算，获得相位-膈肌位移的坐标图。
在仿真实验中采用的两个kTSSM模型是指：
两个kTSSM模型是指本发明步骤6中推导出基于核函数的非线性子空间的kTSSM模型：多项式核的kTSSM模型和高斯核的kTSSM模型。
在仿真实验中采用的一个现有技术偏最小二乘法模型是指：
Malinowski等人在“Mitigating Errors in External Respiratory Surrogate-Based Models of Tumor Position，International Journal of Radiation Oncology,Biology,Physics,vol.82,no.5,2012”中提到的偏最小二乘回归模型，简称PLS模型。
下面结合图3仿真图对本发明的效果做进一步的描述。
图3对应四种算法的膈肌位移的预测结果和实际结果的曲线图，图3中有4个子图，分别对应4种模型的预测结果曲线和实际结果的曲线图，在每个子图中，又有4张图，每张图对应一名测试患者的预测结果和实际结果的曲线图。图3中的x轴表示十个相位，y轴表示膈肌位移值，单位是毫米(mm)。图3中以
Figure GDA0004188921550000161
标示的曲线表示实际结果的曲线，即通过膈肌图像分割，膈肌位移计算得到的值，图3中以“米”字标示的曲线表示预测结果的曲线，即将胸腹表面位移输入预测模型后得到的值。图3(a)为对4名患者采用PLS模型计算出的预测结果和实际结果的曲线图，图3(b)为对4名患者采用线性TSSM模型计算出的预测结果和实际结果的曲线图，图3(c)为4名患者采用多项式核的kTSSM模型计算出的预测结果和实际结果的曲线图，图3(d)为4名患者采用高斯核的kTSSM模型计算出的预测结果和实际结果的曲线图。两条曲线重合的越好，表明该模型的预测效果越好。图3(a)中的4名患者的预测值与实际值之间的误差都明显大于图3(b)，说明TSSM模型相比PLS模型具有更好的预测效果。
3.仿真结果分析。
图4表示在膈肌位移在z轴方向上每个相位的预测性能评价的箱线图，每个相位都包含100次独立运行的结果。图4中，在每个方框上，中心的标记表示中值，而方框的底部和顶部边缘分别表示100次绝对误差结果中的从小到大排列的第25％和第75％的绝对误差结果，“+”符号表示异常值。图4中,x轴表示十个相位，y轴表示预测值和实际值之间的估计误差绝对值，单位是毫米(mm)。图4(a)为对PLS模型进行预测性能评价的箱线图，图4(b)为对线性TSSM模型进行预测性能评价的箱线图，图4(c)为对多项式核的TSSM模型进行预测性能评价的箱线图，图4(d)为对高斯核的TSSM模型进行预测性能评价的箱线图。在图4(a)中的预测误差值比图4(b)、(c)、(d)中的误差值大得多，说明PLS模型相比TSSM模型的预测误差更大。
利用三个指标进行定量评估四种模型的性能：均方误差MSE、R2误差和平均绝对百分比误差MAPE，利用下面公式计算MSE、R2、MAPE，将所有计算结果绘制成表1：
Figure GDA0004188921550000171
Figure GDA0004188921550000172
Figure GDA0004188921550000173
其中，f＝1,2,3,4表示测试集内患者序号，
Figure GDA0004188921550000174
是膈肌实际的位移值，即将膈肌图像从人体三维图像中分割出来计算得到的膈肌位移值，/>
Figure GDA0004188921550000175
是模型的预测值，即将胸腹表面位移值输入进模型后经过计算得到的膈肌位移值，Σ表示求和。MSE衡量模型的预测均方误差，MSE的值越低，预测结果越好。R2比较模型的优劣，消除数据值域对实验结果的影响。R2的范围为[-∞,1]，值越大预测效果越好；MAPE衡量模型的预测值与实际值的平均偏差，平均偏差值越低预测结果越好。
表1主方向预测结果的统计结果一览表
Figure GDA0004188921550000176
Figure GDA0004188921550000181
如表1所示，每个表格的前一个值代表100次独立运行的平均值，(·)中的值代表100次独立运行的标准偏差，本发明仿真实验记录MSE和MAPE的最小值，记录R2的最大值。从表1中可见由顶部到底部，噪声系数依次增加，TSSM模型相比PLS模型，MSE值更小，R2值更大，MAPE值更小，显示了TSSM算法优于PLS模型；通过TSSM算法获得的MSE值和R2值非常令人满意；由于向量值的不平衡性，MAPE值大于100％。

Claims (9)

1.一种基于PCA和TSSM模型的膈肌预测系统，包括膈肌图像分割模块，胸腹表面图像分割模块，膈肌位移计算模块，胸腹表面位移模块，图像输出模块；其特征在于，还包括PCA模块，模型拟合模块；其中，
所述膈肌图像分割模块，用于选取所选患者呼吸周期0％相位的患者体三维图像，将所选相位的人体三维图像以像素值差值法处理为仅有右肺下膈肌的三维图像；
所述胸腹表面图像分割模块，用于选取所选患者的呼吸周期的0％相位的人体三维图像，将该图像以掩膜法处理为仅有胸腹表面的三维图像；
所述膈肌位移计算模块，将膈肌三维图像的全部膈肌像素点的空间三维坐标值取均值，将该均值作为右肺下膈肌的质心坐标值；将该坐标值与所选患者的右肺下膈肌的初始质心坐标值相减，得到所选相位的膈肌位移；
所述胸腹表面位移计算模块，用于用ICP算法将所选患者的所选相位的胸腹部三维图像的空间三维坐标与呼吸周期的0％相位的胸腹表面图的空间三维坐标进行配准，将配准后得到的平移向量作为所选相位的胸腹表面位移；
所述PCA模块，用于将每个患者全部相位的膈肌位移连接为一个10个x轴上的位移值+10个y轴上的位移值+10个z轴上的位移的30维的向量，并进行归一化；用于将每个患者的归一化后的30维的向量，进行主成分分析法处理，降维并映射到d维子空间Rv中；用于将每个患者每个相位的胸腹表面位移归一化为一个10个x轴上的位移值+10个y轴上的位移值+10个z轴上的位移的30维的向量，并进行归一化；用于将每个患者归一化后的30维向量，进行主成分分析法处理，降维并映射到d维子空间Ru中；用于将从TSSM模型中获得的膈肌位移
Figure FDA0004188921540000011
按照PCA特征子空间的系数和基来计算膈肌位移的预测数据，即/>
Figure FDA0004188921540000012
得到归一化的膈肌位移数据/>
Figure FDA0004188921540000013
按照膈肌位移的归一化公式，将归一化的膈肌位移数据/>
Figure FDA0004188921540000014
逆归一化，得到最终的膈肌位移vout；
所述模型拟合模块，用于利用线性岭回归优化法将映射到Rv中的膈肌位移数据与映射到Ru中的胸腹部表面位移数据连接起来，得到变换矩阵β，求出β的最优化表达式βopt，
Figure FDA0004188921540000021
即为TSSM模型；
所述图像输出模块在计算机上显示出横坐标为相位，纵坐标为膈肌位移的曲线图。
2.根据权利要求1所述膈肌预测系统的一种基于PCA和TSSM模型的膈肌预测方法，其特征在于，利用主成分分析PCA为每组数据构建特征子空间，通过线性岭回归优化过程将膈肌数据与胸腹部表面数据连接起来，得到一个子空间映射模型；该膈肌预测方法的步骤包括如下：
(1)呼吸时相数据采集：
(1a)选取患有肺部癌症或胰腺癌的至少20个患者，利用4DCT扫描仪获取每个患者连续的至少5个呼吸周期的时相数据，计算每个患者的平均呼吸周期，从5个呼吸周期中选取最接近平均呼吸周期的一个时相数据，作为该患者的呼吸周期；
(1b)用与4DCT扫描仪连接的计算机打开每个患者的平均呼吸周期时相数据，将该患者呼吸周期从0％开始，每隔10％处选取一次时相数据重建成10个相位的人体三维图像；
(1c)从20个患者中随机选取16个患者，将所选16个患者的人体三维图像包含的全部时相数据集组成初始训练数据，剩余4个患者的人体三维图像作为测试数据；
(2)计算每个患者的4DCT图像的膈肌位移和胸腹表面位移：
(2a)从16个患者中选取一个未选过的患者：
(2b)膈肌图像分割模块选取所选患者呼吸周期0％相位的患者体三维图像，利用像素值差值法，将所选相位的人体三维图像处理为仅有右肺下膈肌的三维图像；
(2c)胸腹表面图像分割模块选取所选患者的呼吸周期的0％相位的人体三维图像，将该图像以掩膜法处理为仅有胸腹表面的三维图像；
(2d)选取所选患者的一个未选过相位的人体三维图像，采用与步骤(2b)相同的处理方法，对所选相位的人体三维图像进行处理，得到所选相位的人体三维图像中仅有右肺下膈肌的膈肌三维图像，膈肌位移计算模块将该膈肌三维图像的全部膈肌像素点的空间三维坐标取均值，将该均值作为右肺下膈肌的质心坐标；将质心坐标值与所选患者的右肺下膈肌的初始质心坐标相减，得到所选相位的膈肌位移；
(2e)采用与步骤(2c)中相同的处理方法，对所选相位的人体三维图像进行处理，得到所选相位的人体三维图像中仅有的胸腹表面的胸腹部三维图像，胸腹表面位移计算模块用ICP算法将所选患者的所选相位的胸腹部三维图像的空间三维坐标值与呼吸周期的0％相位的胸腹表面图的空间三维坐标值进行配准，将配准后得到的平移向量作为所选相位的胸腹表面位移；
(2f)判断是否选完所选患者的所有相位的人体三维图像，若是，则执行步骤(2g)，否则，执行步骤(2d)；
(2g)判断是否选完16个患者，若是，则执行步骤(3)，否则，执行步骤(2a)；
(3)将膈肌位移数据和胸腹表面位移数据映射到各自的子空间中：
(3a)PCA模块将每个患者全部相位的膈肌位移连接为一个10个x轴上的位移值+10个y轴上的位移值+10个z轴上的位移的30维的向量，并进行归一化；
(3b)PCA模块将每个患者的归一化后的30维的向量，进行主成分分析法处理，降维并映射到d维子空间Rv中；
(3c)PCA模块将每个患者每个相位的胸腹表面位移归一化为一个10个x轴上的位移值+10个y轴上的位移值+10个z轴上的位移的30维的向量，并进行归一化；
(3d)PCA模块将每个患者归一化后的30维向量，进行主成分分析法处理，降维并映射到d维子空间Ru中；
(4)模型拟合模块利用线性岭回归优化法将映射到Rv中的膈肌位移数据与映射到Ru中的胸腹部表面位移数据连接起来，得到变换矩阵β，求出β的最优化表达式βopt，
Figure FDA0004188921540000041
即为TSSM模型；
(5)实时输出测试患者的相位-膈肌位移的曲线图：
(5a)用4DCT扫描仪获取测试患者的呼吸时相数据；
(5b)测试患者的呼吸时相数据依次经过：胸腹表面图像分割模块、胸腹表面位移计算模块、PCA模块、TSSM模型、PCA模块的逆归一化处理、最后输出测试患者的膈肌位移值；
(5c)图像输出模块在计算机上显示出横坐标为膈肌位移相位，纵坐标为膈肌位移值的曲线图；
(6)推导出基于核函数的非线性子空间的kTSSM模型：
(6a)用φ表示d维子空间Ru的
Figure FDA0004188921540000042
到高维子空间Wu的映射；用核函数K表示用d维子空间Ru的两个向量内积计算出两个高维子空间Wu中的两个向量的内积；
(6b)用
Figure FDA0004188921540000043
和K(·,·)替换掉TSSM中的/>
Figure FDA0004188921540000044
和两个高维子空间Wu中的两个向量的内积，得到非线性岭回归优化的变换矩阵β'opt和kTSSM模型。
3.根据权利要求2所述的基于PCA和TSSM模型的膈肌预测方法，其特征在于，步骤(2b)中所述的像素值差值法的步骤如下：
第一步，将所选相位的三维人体图像中的最大目标区域作为身体区域；通过形态学方法中的开运算，去除每个相位的三维人体图像中背景中不属于人体，但像素不为0的孤立区域；利用形态学方法中的闭运算填充身体区域中的孔，得到分割后的身体区域；
第二步，从分割后的身体区域中，分割出肺部区域；对三维人体图像，将背景像素值均置为0，分割后的身体区域像素值均置为1，肺部区域像素值均置为3；得到三维人体分割图像；
第三步，用三维ROI掩膜对三维人体分割图像进行分割，得到从前向后看的肺部三维图像；
第四步，对肺部三维图像，从左到右在矢状平面上计算肺部区域的各个切片的面积，生成横坐标为左右方向的三维ROI掩膜切面，纵坐标为肺部面积的曲线图；选择曲线两个峰之间的最低点，作为中位点，定位出该点横坐标的坐标值；
第五步，在三维人体分割图像中，将中位点以左的人体和肺部的像素点全部置为0；再从下到上找到像素值差值为2的分界面，该边界面上的全部像素组成右肺下表面三维图像；
第六步，对该右肺下表面，从上方，先通过形态学运算中的腐蚀去除掉一些属于肺部，但不属于膈肌部分的像素点，再在经过腐蚀后的右肺下表面上填充被过度腐蚀掉的膈肌像素点。
4.根据权利要求2所述的基于PCA和TSSM模型的膈肌预测方法，其特征在于，步骤(2c)中所述的掩膜法的步骤如下：
第1步，对人体三维图像中背景区域的所有像素值均设置为0、身体区域中的所有像素值均设置为1；沿x轴负方向找出像素值差值为1的边界面，将该边界面称为初始胸腹部表面；
第2步，在初始胸腹部表面上沿y轴正方向找到最宽的位置，最宽位置对应人体的肩膀，去除肩膀之上的区域，进行腐蚀操作，得到了每个相位的胸腹表面图；
第3步，将所选患者的所有相位的胸腹表面图取交集，得到胸腹部掩膜，将该掩膜与初始胸腹部表面取交集，得到最终胸腹部表面。
5.根据权利要求2所述的基于PCA和TSSM模型的膈肌预测方法，其特征在于，步骤(2e)中所述ICP算法的步骤如下：
第1步，将呼吸周期0％相位的胸腹表面图的全部胸腹表面像素点的三维坐标作为源点云M，取点集mj∈M；
第2步，将所选相位的胸腹表面图的全部胸腹表面像素点的三维坐标作为目标点云N，取点集nj∈N；
第3步，计算映射矩阵H，使误差函数E(G,l)最小，映射矩阵H和误差函数E(G,l)如下：
Figure FDA0004188921540000061
Figure FDA0004188921540000062
Figure FDA0004188921540000063
其中，j表示点集mj和nj中每个点的序列号，J表示j的最大值，GJ×J表示J阶旋转矩阵，lJ×1表示J维平移向量，O1×J表示J维零向量，
Figure FDA0004188921540000064
表示使{·}取得最小值时的G和l的取值；
第4步，对点集nj进行旋转和平移，计算点集nj'＝{nj'＝Gnj+l,nj∈N}；
第5步，计算映射矩阵H'，使误差函数E(G',l')最小，映射矩阵H'和E(G',l')如下：
Figure FDA0004188921540000065
Figure FDA0004188921540000066
Figure FDA0004188921540000067
第6步，判断误差函数E(G',l')是否小于阈值或者大于预设的最大迭代次数，若是，则停止迭代更新，对所有平移向量求和，即为将所选相位的胸腹表面图的全部胸腹表面像素点与呼吸周期0％相位的胸腹表面图的全部胸腹表面像素点，进行配准后的平移向量；否则，执行第4步。
6.根据权利要求2所述的基于PCA和TSSM模型的膈肌预测方法，其特征在于，步骤(3b)中所述的主成分分析法是指：
第1步，分别将每个患者的10个相位的膈肌位移值中每个相位的x轴、y轴、z轴的位移值，分别组成一个10维向量，再将这3个向量，按照x轴向量+y轴向量+z轴向量，组成一个30维的向量Ui；
采用与第1步相同的处理方法，将每个患者的10个相位的胸腹表面位移值组成一个30维的向量Ui；
第2步，对向量Vi，Ui分别进行归一化处理，得到
Figure FDA0004188921540000071
和/>
Figure FDA0004188921540000072
归一化公式如下：
Figure FDA0004188921540000073
Figure FDA0004188921540000074
其中，
Figure FDA0004188921540000075
表示第i个患者膈肌位移归一化后的30维向量，i表示训练集内患者的序号，Vi(n)表示向量Vi内每一个元素，/>
Figure FDA0004188921540000076
表示第i个患者膈肌位移向量Vi(n)中30个元素中最小的元素，/>
Figure FDA0004188921540000077
表示第i个患者膈肌位移向量Vi(n)中30个元素中最大的元素；/>
Figure FDA0004188921540000078
表示第i个患者胸腹表面位移归一化后的30维向量，Ui(n)表示向量Ui内每一个元素，/>
Figure FDA0004188921540000079
第i个患者胸腹表面位移3向量Ui(n)30个元素中最小的元素，/>
Figure FDA00041889215400000710
表示第i个患者胸腹表面位移向量Ui(n)30个元素中最大的元素；
第3步，按照下式，分别计算协方差矩阵C和B：
Figure FDA00041889215400000711
Figure FDA00041889215400000712
其中，
Figure FDA00041889215400000713
表示第i个患者膈肌位移归一化后的30维向量，/>
Figure FDA00041889215400000714
表示/>
Figure FDA00041889215400000715
的转置，/>
Figure FDA00041889215400000716
表示第i个患者胸腹表面位移归一化后的30维向量，/>
Figure FDA00041889215400000717
表示/>
Figure FDA00041889215400000718
的转置；
第4步，按照下式，计算协方差矩阵B的特征值及对应的特征向量：
det(λI-B)＝0
(λI-B)e＝0
其中，det(·)表示行列式符号，λ表示协方差矩阵B的特征值，I表示和协方差矩阵B阶数相同的单位矩阵，e表示特征值λ对应的特征向量；
按照与第四步相同的处理方式，计算协方差矩阵C的特征值及对应的特征向量；
第5步，对协方差矩阵B的全部特征值，从大到小选择大于阈值ε的d个特征值，阈值ε在构造域的过程中，需要在数据完整性和数据冗余之间进行平衡；
按照与第五步相同的处理方式，选择协方差矩阵C的d个特征值；
第6步，根据d个特征值，得到d个对应的特征向量，该d个特征向量分别作为膈肌位移数据和胸腹部表面位移数据的基，表示为ωv和ωu，ωu,ωv∈R30×d，分别构建出d维子空间Rv和Ru；
第7步，
Figure FDA0004188921540000081
和/>
Figure FDA0004188921540000082
是正交化的/>
Figure FDA0004188921540000083
且/>
Figure FDA0004188921540000084
其中Id是d的单位矩阵；分别通过运算/>
Figure FDA0004188921540000085
和/>
Figure FDA0004188921540000086
将膈肌位移数据和胸腹表面位移数据映射到各自的子空间中，其中，
Figure FDA0004188921540000087
是/>
Figure FDA0004188921540000088
在子空间Ru中的对应点，/>
Figure FDA0004188921540000089
是/>
Figure FDA00041889215400000810
在子空间Rv中的对应点。
7.根据权利要求6所述基于PCA和TSSM模型的膈肌预测方法，其特征在于，步骤(4)中所述的最优表达式βopt操作是指：对β求导，并使等式等于0，则可以获得最优表达式如下：
Figure FDA00041889215400000811
其中，opt表示最优。
8.根据权利要求2所述基于PCA和TSSM模型的膈肌预测方法，其特征在于，步骤(6a)中所述的核函数K的表达式如下：
Figure FDA0004188921540000091
其中，K(·,·)表示映射的核函数，
Figure FDA0004188921540000092
表示测试数据中d维子空间Ru中的胸腹表面位移，/>
Figure FDA0004188921540000093
表示训练数据中d维子空间Ru中的胸腹表面位移，φ表示非线性子空间Ru到高维子空间Wu的映射。
9.根据权利要求8所述基于PCA和TSSM模型的膈肌预测方法，其特征在于，步骤(6b)中所述的β'opt和kTSSM模型表达式如下：
Figure FDA0004188921540000094
Figure FDA0004188921540000095
其中，
Figure FDA0004188921540000096
表示使用kTSSM模型计算后得到的d维子空间Rv中的膈肌位移，/>
Figure FDA0004188921540000097
表示训练数据中的d维子空间Rv中的膈肌位移。
CN202010884814.1A 2020-08-28 2020-08-28 基于pca和tssm模型的膈肌预测系统及方法 Active CN112053330B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010884814.1A CN112053330B (zh) 2020-08-28 2020-08-28 基于pca和tssm模型的膈肌预测系统及方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010884814.1A CN112053330B (zh) 2020-08-28 2020-08-28 基于pca和tssm模型的膈肌预测系统及方法

Publications (2)

Publication Number Publication Date
CN112053330A CN112053330A (zh) 2020-12-08
CN112053330B true CN112053330B (zh) 2023-06-20

Family

ID=73606503

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010884814.1A Active CN112053330B (zh) 2020-08-28 2020-08-28 基于pca和tssm模型的膈肌预测系统及方法

Country Status (1)

Country Link
CN (1) CN112053330B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101790748A (zh) * 2007-06-19 2010-07-28 爱克发医疗保健公司 对3d数字医学图像中的解剖实体进行分割的方法
CN104025119A (zh) * 2011-10-05 2014-09-03 塞伏瑞斯派恩公司 用于手术和介入性医疗过程中的成像系统和方法
CN111161333A (zh) * 2019-12-12 2020-05-15 中国科学院深圳先进技术研究院 一种肝脏呼吸运动模型的预测方法、装置及存储介质

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11612338B2 (en) * 2013-10-24 2023-03-28 Breathevision Ltd. Body motion monitor

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101790748A (zh) * 2007-06-19 2010-07-28 爱克发医疗保健公司 对3d数字医学图像中的解剖实体进行分割的方法
CN104025119A (zh) * 2011-10-05 2014-09-03 塞伏瑞斯派恩公司 用于手术和介入性医疗过程中的成像系统和方法
CN111161333A (zh) * 2019-12-12 2020-05-15 中国科学院深圳先进技术研究院 一种肝脏呼吸运动模型的预测方法、装置及存储介质

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
GND-PCA-Based Statistical Modeling of Diaphragm Motion Extracted from 4D MRI;Windra Swastika et al;《Computational and Mathematical Methods in Medicine》;20130526;1-10页 *
介入式治疗中靶区呼吸运动实时跟踪技术;孙百权等;《中国医学物理学杂志》;20180430;第35卷(第4期);374-376页 *
基于四维CT的肝脏各部分和膈肌的呼吸移动度的分析;龚敏等;《中国临床医学》;20120825(第04期);436-443页 *

Also Published As

Publication number Publication date
CN112053330A (zh) 2020-12-08

Similar Documents

Publication Publication Date Title
De Craene et al. Temporal diffeomorphic free-form deformation: Application to motion and strain estimation from 3D echocardiography
US9076201B1 (en) Volumetric deformable registration method for thoracic 4-D computed tomography images and method of determining regional lung function
Zacharaki et al. Non-diffeomorphic registration of brain tumor images by simulating tissue loss and tumor growth
Ibragimov et al. Segmentation of pathological structures by landmark-assisted deformable models
Heyde et al. Anatomical image registration using volume conservation to assess cardiac deformation from 3D ultrasound recordings
Baka et al. Statistical shape model-based femur kinematics from biplane fluoroscopy
Yu et al. Unsupervised 3D PET-CT image registration method using a metabolic constraint function and a multi-domain similarity measure
CN115116586A (zh) 一种基于联合配准的可变形统计图谱构建方法
Liu et al. Shape-correlated deformation statistics for respiratory motion prediction in 4D lung
KR101460908B1 (ko) 4차원 컴퓨터 단층촬영 영상의 폐종양 위치 추적 시스템 및 그 방법
CN112053330B (zh) 基于pca和tssm模型的膈肌预测系统及方法
Xue et al. Lung respiratory motion estimation based on fast Kalman filtering and 4D CT image registration
De Craene et al. Temporal diffeomorphic free-form deformation for strain quantification in 3D-US images
Urschler et al. Assessing breathing motion by shape matching of lung and diaphragm surfaces
Lafitte et al. Accelerating multi-modal image registration using a supervoxel-based variational framework
Taubmann et al. Prediction of respiration-induced internal 3-D deformation fields from dense external 3-D surface motion
Li et al. 3D intersubject warping and registration of pulmonary CT images for a human lung model
Cao et al. Intensity-and-landmarkdriven, inverse consistent, B-Spline registration and analysis for lung imagery
Tsuzuki et al. Animated solid model of the lung constructed from unsynchronized MR sequential images
Ouyang et al. Preliminary feasibility study of imaging registration between supine and prone breast CT in breast cancer radiotherapy using residual recursive cascaded networks
Coevoet et al. Introducing interactive inverse FEM simulation and its application for adaptive radiotherapy
Zhang et al. Temporal consistent 2D-3D registration of lateral cephalograms and cone-beam computed tomography images
Negahdar et al. Tracking planar lung motion in 4D CT with optical flow: validations and comparison of global, local, and local-global methods
Roose et al. Adaptive boundary conditions for physically based follow-up breast MR image registration
Peressutti et al. Personalising population-based respiratory motion models of the heart using neighbourhood approximation based on learnt anatomical features

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