本发明公开了一种本构模型与有限元结合的脉络膜新生血管生长预测方法，图像预处理；区域提取和划分：将图像分割为CNV区域、外视网膜层、内视网膜层和脉络膜层4个区域；网格化：对4个区域进行四面体网格生成；建模：运用超弹性生物力学模型与反应扩散方程建模，将脉络膜新生血管生长后的质量改变作为源项加入到方程中，使变形梯度张量根据新生血管的生长持续变化；优化模型，算出最佳准确率，进行参数检验；根据每个时间点预测的参数拟合一条参数曲线，预测最后一个时间点的生长参数，得到预测结果。本发明的方法能够更加灵活和个性化的模式生物机械模型，模型假设的组织是正交各向异性，对非线性、大变形区域也有很好的预测结果，准确度高。
1.一种本构模型与有限元结合的脉络膜新生血管生长预测方法，其特征是，包括以下步骤：步骤1.图像预处理：对三维视网膜频域光学相干断层成像SD-OCT图像进行下采样处理；
步骤2.区域提取和划分：将图像分割为脉络膜新生血管CNV区域、外视网膜层、内视网膜层和脉络膜层4个区域；
步骤3.网格化：对CNV区域、内视网膜层、外视网膜层和脉络膜层进行多尺度、自适应的四面体网格生成；
步骤4.运用超弹性生物力学模型与反应扩散方程建模：在超弹性生物力学模型结构中，整合各区域不同的体积模量、剪切模量及其对应的体积弹性响应、等容弹性响应，将脉络膜新生血管生长后的质量改变作为源项加入到方程中，使变形梯度张量根据新生血管的生长持续变化；
步骤6.根据每个时间点预测的参数拟合一条参数曲线，预测最后一个时间点的生长增值率参数，得到预测结果；
其中，t为时间，D是各向异性的扩散张量，包含三个组成部分Dx，Dy，Dz，这三个分量分别为x轴，y轴，z轴的分量；ρ代表生长增殖率，N是每单位体积的细胞数，Kmax代表最大细胞容量，当血管内皮生长因子的各向异性浓度 时，有如下公式：(4b)质量效应与超弹性生物力学模型
其中， 表示格林-拉格朗日应变张量，FT代表变形梯度张量F的转置张量，I为单位矩阵；J代表变形梯度张量F的行列式，即J＝det(F)；κ和μ分别代表体积模量和剪切模量；Tr表示矩阵的迹；I1是柯西-格林变形张量的第一主不变量；
通过公式4中的梯度力f，将柯西应力张量σ＝J-1FSFT和血管内皮生长因子的各向异性浓度c联系起来：div(σ)+f＝0；f＝-ξ▽c (4)其中，f代表由c产生的梯度力即密度梯度力，ξ代表一个需要被估计的常数，该常数取决于材料的生物学特性；S代表PK2-第二Piola-Kirchhoff应力张量。
2.根据权利要求1所述的本构模型与有限元结合的脉络膜新生血管生长预测方法，其特征是，步骤1中的图像预处理包括以下步骤：(1a)下采样处理：
(1b)使用CAVASS软件，进行OCT图像配准。
3.根据权利要求1所述的本构模型与有限元结合的脉络膜新生血管生长预测方法，其特征是，步骤3中采用一个三维曲面和体积的网格生成器对CNV区域、内视网膜层、外视网膜层和脉络膜层进行多尺度、自适应的四面体网格生成。
4.根据权利要求1所述的本构模型与有限元结合的脉络膜新生血管生长预测方法，其特征是，步骤5中，以一阶马尔科夫序列优化生长增值率ρ参数：前n个时间点的数据，每个时间点均与之前的一个时间点相关，计算ρ的值。
5.根据权利要求4所述的本构模型与有限元结合的脉络膜新生血管生长预测方法，其特征是，从之前所有的时间点进行生长增值率参数的估计，然后使用下面的公式来计算预测的准确性：是基于之前所有时间点Ii,c(i＝1,2,...，10)预测得到的结果，Ii+1,c是当前时间点的金标准；
最后，估计参数：νi＝{Dx,Dy,Dz,ρ；c}。
6.根据权利要求5所述的本构模型与有限元结合的脉络膜新生血管生长预测方法，其特征是，定义预测结果与金标准的重合像素点计算反映预测结果准确性的各项指标：
7.根据权利要求4、5或6所述的本构模型与有限元结合的脉络膜新生血管生长预测方法，其特征是，步骤6中，利用得出的前n个时间点的生长增值率ρ值，画出拟合曲线，从而预测出第n+1个时间点的生长增值率ρ；将预测出的第n+1个时间点的生长增值率ρ代入第n个时间点，得出的预测图像与真实的第n+1个图像重叠，计算预测准确率。
本构模型与有限元结合的脉络膜新生血管生长预测方法
技术领域
[0001] 本发明涉及计算机视觉、图像处理与分析技术领域，属于生长预测的建模方法，尤其是应用于SD-OCT(频域光学相干断层成像)的脉络膜新生血管(CNV)生长预测的建模方法。
背景技术
[0002] 脉络膜位于视网膜和巩膜之间，由纤维组织、小血管和毛细血管组成。脉络膜的血循环营养视网膜外层、晶状体和玻璃体等，由于脉络膜内血流量大、流速较慢、而其本身软而薄的特性，使其对外界影响或冲击不耐受，导致病原体在此处易滞留，造成脉络膜疾病及各种眼底疾病，因此，研究脉络膜新生血管(CNV)的生长可以得到很多重要的数据，为后续医学研究提供基本的支撑数据。
[0003] 由于新生血管的生长部位在对外界环境极其敏感和脆弱的眼部视网膜内，实时入侵性的对新生血管进行检测显然是不可能的。现有的新生血管生长预测方法通常仅采用单一的增长模型如有限元法进行预测，这种方法存在一定的缺点和局限性：①基于各向同性的假设，②没有考虑到不同组织结构的不同生物学特性，③只能预测线性变化，④只能预测小变形的变化，对大变形变化的预测不准确，等等。
发明内容
[0004] 本发明所要解决的技术问题是克服现有技术的缺陷，本发明首次提出了一个综合性的方法——基于超弹性材料的生物力学模型和反应扩散方程相结合的方法，进行脉络膜新生血管的生长预测分析，准确度高。
[0005] 为解决上述技术问题，本发明提供一种脉络膜新生血管三维生成预测方法。
[0006] 该方法主要包括6个步骤：
[0007] 步骤1.图像预处理：对三维视网膜频域光学相干断层成像SD-OCT图像进行下采样处理；
[0008] 步骤2.区域提取和划分：将图像分割为脉络膜新生血管CNV区域、外视网膜层、内视网膜层和脉络膜层4个区域；
[0009] 步骤3.网格化：对CNV区域、内视网膜层、外视网膜层和脉络膜层进行多尺度、自适应的四面体网格生成；
[0010] 步骤4.运用超弹性生物力学模型与反应扩散方程建模：在超弹性生物力学模型结构中，整合各区域不同的体积模量、剪切模量及其对应的体积弹性响应、等容弹性响应，将脉络膜新生血管生长后的质量改变作为源项加入到方程中，使变形梯度张量根据新生血管的生长持续变化；
[0011] 步骤5.优化模型，算出最佳准确率，进行参数检验；
[0012] 步骤6.根据每个时间点预测的参数拟合一条参数曲线，预测最后一个时间点的生长增值率参数，得到预测结果。
[0013] 步骤1中的图像预处理包括以下步骤：
[0014] (1a)下采样处理：
[0015] (1b)使用CAVASS软件，进行OCT图像配准。
[0016] 步骤3中采用一个三维曲面和体积的网格生成器对CNV区域、内视网膜层、外视网膜层和脉络膜层进行多尺度、自适应的四面体网格生成。
[0017] 步骤4包括以下步骤：
[0018] (4a)反应扩散方程的入侵和扩散
[0019] 逻辑回归型增长的反应扩散方程表示为：
[0020]
[0021] 其中，t为时间，D是各向异性的扩散张量，包含三个组成部分Dx，Dy，Dz，这三个分量分别为x轴，y轴，z轴的分量；ρ代表生长增殖率，N是每单位体积的细胞数，Kmax代表最大细胞容量，
[0022] 当血管内皮生长因子的各向异性浓度 时，有如下公式：
[0023]
[0024] (4b)质量效应与超弹性生物力学模型
[0025] 材料的非线性应力-应变行为的应变能密度函数为：
[0026]
[0027] 其中， 表示格林-拉格朗日应变张量，FT代表变形梯度张量F的转置张量，I为单位矩阵；J代表变形梯度张量F的行列式，即J＝det(F)；κ和μ分别代表体积模量和剪切模量；Tr表示矩阵的迹；I1是柯西-格林变形张量的第一主不变量；
[0028] 通过公式4中的梯度力f，将柯西应力张量σ＝J-1FSFT和血管内皮生长因子的各向异性浓度c联系起来：
[0029]
[0030] 其中，f代表由c产生的梯度力即密度梯度力，ξ代表一个需要被估计的常数，该常数取决于材料的生物学特性；S代表PK2-第二Piola-Kirchhoff应力张量。
[0031] 步骤5中，以一阶马尔科夫序列优化生长增值率ρ参数：前n个时间点的数据，每个时间点均与之前的一个时间点相关，计算ρ的值。
[0032] 从之前所有的时间点进行生长增值率参数的估计，然后使用下面的公式来计算预测的准确性：
[0033]
[0034] 是基于之前所有时间点Ii,c(i＝1,2,...，10)预测得到的结果，Ii+1,c是当前时间点的金标准；
[0035] 最后，估计参数：νi＝{Dx,Dy,Dz,ρ；c}。
[0036] 定义预测结果与金标准的重合像素点
[0037] 计算反映预测结果准确性的各项指标：
[0038] (1)召回率
[0039] (2)错误率
[0040] (3)相似度DICE系数
[0041] 步骤6中，利用得出的前n个时间点的生长增值率ρ值，画出拟合曲线，从而预测出第n+1个时间点的生长增值率ρ；将预测出的第n+1个时间点的生长增值率ρ代入第n个时间点，得出的预测图像与真实的第n+1个图像重叠，计算预测准确率。
[0042] 本发明所达到的有益效果：
[0043] 本发明首次提出了一个综合性的生长预测方法——基于超弹性材料的生物力学模型和反应扩散方程相结合的方法，进行脉络膜新生血管的生长预测分析。
[0044] 本发明的方法结合超弹性材料和反应反应扩散方程，能够更加灵活和个性化的模式生物机械模型。不同于以往的材料是各向同性假设，该模型假设的组织是正交各向异性。在疾病部位(脉络膜新生血管区域)和周围的组织(包括外层视网膜和脉络膜，视网膜内层等)模拟出不同组织结构的不同生物学特性，且基于正交各向异性的假设，对非线性、大变形区域也有很好的预测结果，准确度高。
附图说明
[0045] 图1是脉络膜新生血管预测的整体框图；
[0046] 图2a是OCT的原始图像(512×1024×128)网格化后的结果；
[0047] 图2b是下采样的图像(64×64×64)网格化后的结果；
[0048] 图2c表示原始图像的节点信息和单元信息；
[0049] 图2d表示下采样图像的节点信息和单元信息；
[0050] 图3a是从下采样数据重构的网格划分；
[0051] 图3b是脉络膜新生血管的细胞内体积(ICV)；
[0052] 图3c是ICV的等值面图像；
[0053] 图3d是应力分布的显示：每一个切片上，CNV区域的粗线条圆形和细线条椭圆形分别代表强应力(一般分布在CNV边缘)和弱应力(一般分布在CNV内部)的位置；
[0054] 图4是不同病人的CNV体积差异示意图。
具体实施方式
[0055] 下面对本发明的具体实施方式作进一步详细的描述：
[0056] 本方法的基本框图如图1所示，主要包括以下步骤：
[0057] 1)图像预处理
[0058] (a)下采样处理：为了提高处理效率，优化内存使用，本实施例中将每个像素512*1024*128的三维视网膜SD-OCT图像下采样到像素64*64*64。
[0059] (b)使用CAVASS软件(MIPG开发的软件)，进行OCT图像配准。
[0060] 2)区域提取和划分
[0061] 手动标记金标准，进行视网膜感兴趣区域(CNV区域)的提取、周围组织(内视网膜层、外视网膜层、脉络膜层)的划分；
[0062] 3)网格化，用一个三维曲面和体积的网格生成器对CNV区域、内视网膜层、外视网膜层和脉络膜层进行多尺度、自适应的四面体网格生成，见图2a、图2b、图2c、图2d；
[0063] 4)对CNV区域、内视网膜层、外视网膜层和脉络膜层这4个区域建模，其中各个区域进行不同的超弹性生物力学模型与反应扩散方程模型进行参数训练，使其根据不同的生理学特性进行生长：
[0064] (a)反应扩散方程的入侵和扩散
[0065] 逻辑回归型(Logistic)增长的反应扩散方程可以表示为：
[0066]
[0067] 其中，t为时间，D是各向异性的扩散张量，包含三个组成部分Dx，Dy，Dz，这三个分量分别为x轴，y轴，z轴的分量。ρ代表生长增殖率，N是每单位体积的细胞数，Kmax代表最大细胞容量，当血管内皮生长因子(VEGF)的各向异性浓度 时，有如下公式：
[0068]
[0069] (b)质量效应与超弹性生物力学模型
[0070] 在模型结构中，通过方程(3)比较全面地整合各区域不同的体积模量、剪切模量及其在本构方程中对应的体积弹性响应部分、等容弹性响应部分，将CNV生长后的质量改变作为源项加入到超弹性生物力学的本构方程中，使变形梯度张量根据新生血管的生长持续变化：
[0071] 超弹性的生物力学模型可用于较大变形的预测，材料的非线性应力-应变行为的应变能密度函数为：
[0072]
[0073] 其中， 表示格林-拉格朗日应变张量，FT代表变形梯度张量F的转置张量，其分量关系为(FT)ij＝Fji,I为单位矩阵；J代表变形梯度张量F的行列式，即J＝det(F)；κ和μ分别代表体积模量和剪切模量；Tr表示矩阵的迹；I1是柯西-格林(Cauchy-Green)变形张量的第一主不变量。
[0074] 数据重构的网格划分以及各个切片上的应力分布，见图3a、图3b、图3c、图3d。
[0075] 同时，使用静态平衡方程来模拟该过程：
[0076]
[0077] 该方程通过梯度力f，将柯西应力张量σ＝J-1FSFT(其中S代表PK2-第二Piola-Kirchhoff应力张量)，和血管内皮生长因子的各向异性浓度c联系起来：
[0078] 其中，f代表由c产生的梯度力(即密度梯度力)，ξ代表一个需要被估计的常数，该常数取决于材料的生物学特性。
[0079] 5)优化模型，算出最佳准确率，进行参数检验。
[0080] 以一阶马尔科夫序列优化参数ρ(生长增值率)：前n(本实施例中设n＝10)个时间点的数据，每个时间点均与之前的一个时间点相关，计算ρ的值。
[0081] 从之前所有的时间点进行模型参数(生长增值率)的估计，然后使用下面的公式来计算预测的准确性：
[0082]
[0083] 是基于之前所有时间点Ii,c(i＝1,2,...，n)预测得到的结果，Ii+1,c是当前时间点的金标准，w1、w2分别代表权重。
[0084] 根据脉络膜新生血管的生理特性，设置CNV、脉络膜(choroid)和其他部分的参数κ和μ分别为：κCNV＝5kPa,κchoroid＝6kPa,κothers＝3kPa,μCNV＝μchoroid＝5*103N/m2,μothers＝1*103N/m2。最后，我们要估计的参数是：νi＝{Dx,Dy,Dz，ρ；c}。定义预测结果与金标准的重合像素点 计算反映预测结果准确性的各项指标，这些指标就作为评价和检验模型
的标准：(1)召回率 (2)错误率 (3)相似度DICE系数
[0085] 6)根据每个时间点预测的参数拟合一条参数曲线，预测最后一个时间点的生长参数，得到预测结果。
[0086] 利用得出的前n(本实施例中设n＝10)个时间点的生长增值率ρ值，画出拟合曲线，从而预测出第n+1个时间点的生长增值率ρ。将预测出的第n+1个时间点的生长增值率ρ代入第n个时间点，得出的预测图像与真实的第n+1个图像重叠，计算预测准确率。
[0087] 7)实验结果：
[0088] 如图4中所示，为10个时间点的六位病人的CNV体积，从图中可以看出，CNV的特异性明显，波动各不相同。实际上，不同病人不同时间点的CNV不管在体积大小上，还是形态上都差异很大。
[0089] 由于视网膜对光的敏感性，还在结果中增加了应力分布的显示(见后附图3d)，有助于更加直观的显示力学对CNV区域的影响。
[0090] 根据第一时间点上的变形，预测的第二时间点参数；根据前两个时间点的变形，预测第三个时间点的形变参数。以此类推，将前n(n＝10)个时间点的生长增值率ρ拟合出的曲线，代入第n个时间点，得出预测结果 将其与真实值(n+1时间点的金标准)比较，得出预测结果,如表1。
[0091]
[0092] 表1是本方法的预测结果，可以发现，除了个别CNV前后时间点变化差异巨大的病人外，其余的预测结果都得到了比较好的准确度。
[0093] 至此，一种三维的CNV生长预测模型已经实现并得到验证。本发明提出了一种融合超弹性生物力学模型、反应扩散方程和质量源的思想，结合目标函数的优化参数自动检测，且在结果中显示全区域的应力分布的技术方案，实验结果表明，本发明的生长预测方法能够比较好的解决脉络膜新生血管(CNV)的生长预测问题，具有较高的预测准确度。
[0094] 以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。
