CN103049663A - 磁共振弹性成像中的弹性模量重建方法和系统 - Google Patents
磁共振弹性成像中的弹性模量重建方法和系统 Download PDFInfo
- Publication number
- CN103049663A CN103049663A CN2012105724496A CN201210572449A CN103049663A CN 103049663 A CN103049663 A CN 103049663A CN 2012105724496 A CN2012105724496 A CN 2012105724496A CN 201210572449 A CN201210572449 A CN 201210572449A CN 103049663 A CN103049663 A CN 103049663A
- Authority
- CN
- China
- Prior art keywords
- displacement
- elastic modulus
- subregion
- unit
- volume element
- 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
Links
Images
Landscapes
- Magnetic Resonance Imaging Apparatus (AREA)
Abstract
本发明提供了一种磁共振弹性成像中的弹性模量重建方法和系统。所述方法包括:将成像组织的表面假设为平面域,通过有限体积元算法求解得到所述平面域中的位移初值,以及所述平面域中子区域的位移、所述位移对未知弹性模量的导数;通过所述有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值,直至所述弹性模量值对应的最优化问题平方差小于预设的容忍误差,且达到预设的迭代次数时停止牛顿迭代;将所述最终迭代得到的弹性模量值组成所述成像组织的弹性模量分布。采用本发明能降低计算量。
Description
技术领域
本发明涉及磁共振成像技术,特别是涉及一种磁共振弹性成像中的弹性模量重建方法和系统。
背景技术
弹性是人体组织物理性质中一种重要的机械力学参数,生物组织的弹性变化通常是与一定的病理现象紧密相关的,也就是说,病变组织和正常组织往往存在着弹性模量的差异,这一差异为临床上疾病的诊断提供了重要的参考信息。
磁共振弹性成像(Magnetic Resonance Elastography,简称MRE)作为一种无创成像方法,能够直观地显示和量化人体内部组织弹性,实现对人体内部组织的弹性成像,使得“影像触诊”成为了可能,在乳腺癌检测、肝硬化分期、动脉粥样硬化斑块、肌肉损伤、大脑疾病检测和射频消融等治疗和监控方面具有重要意义。
磁共振弹性成像中弹性模量重建的方法是一个由质点位移图反推弹性分布的逆问题求解,因此,其本质上是不稳定的。为了避免该问题的病态性,弹性弹性模块重建方法将根据应用范围进行假设和简化。目前提出的弹性模量重建方法包括:(1)局部频率估计(Local Frequency Estimation,简称LFE)算法及其变种,该算法将假设介质是均匀的和不可压缩的,并忽略波动中的衰减,机械波在介质中的传播方程因而简化为亥姆霍兹方程,以该方程为模型进行直接逆问题代数求解,但是局部频率估计算法存在着分辨率低、精度有限的缺陷,对尖锐的边界无法估计出精确的弹性系数,其假设也不适用于某些临床中;(2)基于有限元分析的弹性模量重建算法,计算出一幅质点位移图,通过最小化该质点位移图和磁共振质点位移图得到弹性系数分布图,与局部频率估计算法及其变种相比较,该方法对介质等没有做特定假设,对噪声的干扰不敏感,可产生较高分辨率的图像,但计算量非常庞大。
发明内容
基于此,提供一种能降低计算量的磁共振弹性成像中的弹性模量重建方法。
此外,还有必要提供一种能降低计算量的磁共振弹性成像中的弹性模量重建系统。
一种磁共振弹性成像中的弹性模量重建方法,包括如下步骤:
将成像组织的表面假设为平面域,通过有限体积元算法求解得到所述平面域中的位移初值,以及所述平面域中子区域的位移、所述位移对未知弹性模量的导数;
通过所述有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值,直至所述弹性模量值对应的最优化问题平方差小于预设的容忍误差,且达到预设的迭代次数时停止牛顿迭代;
将所述最终迭代得到的弹性模量值组成所述成像组织的弹性模量分布。
在其中一个实施例中,所述通过有限体积元算法求解得到所述平面域中的位移初值的步骤为:
根据设定的弹性模量初值通过有限体积元算法计算得到平面域中的位移初值。
在其中一个实施例中,所述通过有限体积元算法求解得到所述平面域中子区域的位移、所述位移对未知弹性模量的导数的步骤包括:
将所述平面域划分为若干个子区域,并划分所述子区域为若干个单元,所述子区域之间和单元之间不存在重叠且任一单元的顶点均不在其它单元的边上,平面域边界的顶点为单元的顶点;
在所述子区域中,以构成单元的顶点作为所述子区域的节点,并构建所述节点对应的有限体积元方程,并以所述位移初值作为有限体积元方程中的初值进行求解得到所述子区域对应的位移;
通过所述求解得到的位移对包含了所述位移对未知弹性模量的导数的方程进行求解得到所述位移对未知弹性模量的导数,所述位移对未知弹性模量的导数是与位移所在的子区域相对应的。
在其中一个实施例中,所述通过所述有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值的步骤包括:
通过每一子区域所对应的位移和导数进行运算得到所述子区域对应的弹性模量改进值;
根据所述弹性模量改进值以设定的弹性模量初值为起始进行牛顿迭代得到与当前迭代次数对应的弹性模量值。
在其中一个实施例中,所述根据所述弹性模量改进值以设定的弹性模量初值为起始进行牛顿迭代得到与当前迭代次数对应的弹性模量值的步骤之后还包括:
获取通过对所述成像组织进行磁共振成像所得到的位移图;
根据当前迭代得到的弹性模量值和所述位移图得到最优化问题平方差;
判断所述最优化问题平方差是否小于预设的容忍误差,若是,则进一步判断所述平面域中的每一单元是否均包含于至少一个子区域中,若是,则
判断所述单元中的最小单元对应的迭代次数是否达至预设的迭代次数,若是,则
停止进行牛顿迭代。
在其中一个实施例中,还包括:
若判断到所述单元中的最小单元对应的迭代次数未达到预设的迭代次数,则
返回所述通过有限体积元算法求解得到所述平面域中的位移初值,以及所述平面域中子区域的位移、所述位移对未知弹性模量的导数的步骤。
一种磁共振弹性成像中的弹性模量重建系统,包括:
有限体积元运算模块,用于将成像组织的表面假设为平面域,通过有限体积元算法求解得到所述平面域中的位移初值,以及所述平面域中子区域的位移、所述位移对未知弹性模量的导数;
迭代模块,用于通过所述有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值,直至所述弹性模量值对应的最优化问题平方差小于预设的容忍误差,且达到预设的迭代次数时停止牛顿迭代;
分布形成模块,用于将所述最终迭代得到的弹性模量值组成所述成像组织的弹性模量分布。
在其中一个实施例中,所述有限体积元运算模块还用于根据设定的弹性模量初值通过有限体积元算法计算得到平面域中的位移初值。
在其中一个实施例中,所述有限体积元运算模块包括:
划分单元,用于将所述平面域划分为若干个子区域,并划分所述子区域为若干个单元,所述子区域之间和单元之间不存在重叠且任一单元的顶点均不在其它单元的边上,平面域边界的顶点为单元的顶点;
子区域位移求解单元,用于在所述子区域中,以构成单元的顶点作为所述子区域的节点,并构建所述节点对应的有限体积元方程,并以所述位移初值作为有限体积元方程中的初值进行求解得到所述子区域对应的位移;
子区域导数求解单元,用于通过所述求解得到的位移对包含了所述位移对未知弹性模量的导数的方程进行求解得到所述位移对未知弹性模量的导数,所述位移对未知弹性模量的导数是与位移所在的子区域相对应的。
在其中一个实施例中,所述迭代模块包括:
改进值运算单元,用于通过每一子区域所对应的位移和导数进行运算得到所述子区域对应的弹性模量改进值;
弹性模量迭代单元,用于根据所述弹性模量改进值以设定的弹性模量初值为起始进行牛顿迭代得到与当前迭代次数对应的弹性模量值。
在其中一个实施例中,所述迭代模块还包括:
位移图处理单元,用于获取通过对所述成像组织进行磁共振成像所得到的位移图,并根据当前迭代得到的弹性模量值和所述位移图得到最优化问题平方差;
判断单元,用于判断所述最优化问题平方差是否小于预设的容忍误差,若是,则进一步判断所述平面域中的每一单元是否均包含于至少一个子区域中,若是,则判断所述单元中的最小单元对应的迭代次数是否达到预设的迭代次数,若是,则停止进行牛顿迭代。
在其中一个实施例中,所述判断单元还用于若判断到所述单元中的最小单元对应的迭代次数未达到预设的迭代次数,则通知所述有限体积元运算模块。
上述磁共振弹性成像中的弹性模量重建方法和系统,将成像组织的表面假设为平面域,引入有限体积元算法以求解得到位移和该位移对未知弹性模量的导数,进而通过求解得到的位移和导数进行牛顿迭代得到弹性模量值,并在弹性模量值对应的最优化问题平方差小于预设的容忍误差,且达到预设的迭代次数时停止牛顿迭代,进而将最终迭代得到的弹性模量值组成成像组的弹性模量分布,通过引入有限体积元算法,使得弹性模量的重建在保证了计算精度的前提下降低了计算量,提高了计算速度。
附图说明
图1为一个实施例中磁共振弹性成像中的弹性模量重建方法的流程图;
图2为图1中通过有限体积元算法求解得到平面域中子区域的位移、位移对未知弹性模量的导数的方法流程图;
图3为一个实施例中子区域的原始剖分的示意图;
图4为一个实施例中对子区域的某一节点的外心对偶剖分单元的示意图;
图5为一个实施例中通过有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值的方法流程图;
图6为另一个实施例中通过有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值的方法流程图;
图7为一个实施例中平面域的示意图;
图8为一个实施例中磁共振弹性成像中的弹性模量重建系统的结构示意图;
图9为图8中有限体积元运行模块的结构示意图;
图10为一个实施例中迭代模块的结构示意图;
图11为另一个实施例中迭代模块的结构示意图。
具体实施方式
如图1所示,一种磁共振弹性成像中的弹性模量重建方法,包括如下步骤:
步骤S10,将成像组织的表面假设为平面域,通过有限体积元算法求解得到平面域中的位移初值,以及平面域中子区域的位移、位移对未知弹性模量的导数。
本实施例中,在一定条件下,可将一个在外力作用下处于平衡状态的弹性体视为平面弹性问题,而在磁共振弹性成像中,成像组织被施加的外力是非常小的,可将该成像组织视为处于平衡状态的,因此,可将磁共振弹性成像中应力和应变之间的关系视为平面弹性问题。
将成像组织的表面视为平面域Ω,获取通过平面域对应的平衡方程所构建得到的二阶椭圆微分方程组,并在构建的二阶椭圆微分方程组中对未知弹性模量求导得到包含了位移对未知弹性模量的导数的方程。
进一步的,对于平面域Ω,其边界为描述成像组织平衡的状态变量包括三组,即,应力张量σ=(σ11,σ22,σ12)T、应变张量ε=(ε11,ε22,ε12)T和位移向量u=(u1,u2)T。假设Ω是均匀各向同性的弹性体,令
其中,正数λ和μ是Lam é常数:
其中,ν是Poisson系数(Poisson’s Ratio,是已知的常数),E是弹性模量(Young’sModulus)。σ,ε,u满足如下三组方程:
σ=Aε(应力-应变关系)(3)
其中,f是体积力。
由Green公式可推出:
在求解时,利用公式(1)-(3)消去部分未知量,留下的是待求量,即弹性模量E,得到包含了位移向量u的二阶椭圆微分方程组:
其中
在边界Γ1上有
其中,H1(Ω)是平面域Ω上的Hilbert空间.
进一步的,为求解出位移对未知弹性模量的导数,将对(4)式两端对未知弹性模量求导,以得到包含了位移对未知弹性模量的导数的方程。
在构建得到了二阶椭圆微分方程组和包含了位移对未知弹性模量的导数的方程之后,将应用有限体积元算法求解位移和位移对未知弹性模量的导数。
有限体积元法是吸取了有限差分法和有限元法的优点发展起来的,又称广义差分法,可以处理复杂的边值条件及不规则区域,并且由于离散化得到的线性方程组是稀疏的,在不降低数值解的收敛阶的前提下求解方程的计算量小,计算速度快,因此,将有限体积元法应用于磁共振弹性成像中将大大地降低了弹性模量重建的计算量。
在一个实施例中,上述通过有限体积元算法求解得到平面域中的位移初值的具体过程为:根据设定的弹性模量初值通过有限体积元算法计算得到平面域中的位移初值。
本实施例中,将首先进行位移初值的求解,以方便后续的位移求解过程。在位移初值的求解过程中,将进行弹性模量初值的假设,并根据这一假设通过有限体积元算法得到平面域所对应的位移初值,其中,该位移初值计算时所使用的边界Γ上的初值取为对成像组织进行磁共振成像所得到的位移图中对应的位移值。
如图2所示,在一个实施例中,上述通过有限体积元算法求解得到平面域中子区域的位移、位移对未知弹性模量的导数的步骤包括:
步骤S110,将平面域划分为若干个子区域,并划分子区域为若干个单元,子区域之间和单元之间不存在重叠且任一单元的顶点均不在其它单元的边上,平面域边界的顶点为单元的顶点。
本实施例中,将平面域Ω划分为若干个子区域,然后分别划分每一子区域为若干个单元之和,在优选的实施例中,划分得到的单元为三角形的形状,进而使得子区域之间、单元之间互相不存在重叠区域,并且任一单元的顶点不会在其它单元的边上,边界Γ的每一顶点均为单元的顶点,从而每一子区域将对应了一个原始剖分Th,如图3所示,其中h是所有三角形单元边的最长边。记KQ为Th内的三角形单元,Q是三角形单元的外心。
进一步的,如图4所示,记是原始剖分Th的对偶剖分,包括是重心对偶剖分和外心对偶剖分。以下设是外心对偶剖分。图4表示的是原始剖分Th中以P0为顶点的所有三角单元和围绕P0的对偶单元(图中阴影剖分)。外心对偶剖分做法:设Th的任一单元的内角不大于90°,取三角单元□P0PiPi+1(i=1,2,…6,P7=P1)的外心Qi为对偶剖分的节点,即依次连结各个外心点可得到外心对偶剖分,此时,是的中垂线,分别过边中点Mi。
步骤S130,在子区域中,以构成单元的顶点作为子区域的节点,并构建节点对应的有限体积元方程,并以位移初值作为有限体积元方程中的初值进行求解得到子区域对应的位移。
本实施例中,对于每一子区域,设是定义在原始分剖Th上的试探函数空间,它是分片线性函数空间,即为一次多项式,完全由三角形单元中三个顶点上的值所确定,将Th的内节点编号为1,2,…,N0。其中,节点分两类,给定力条件的节点编号为N0+1,…,N1,给定位移条件的节点编号为N1+1,…,N。用表示节点i∈{1,2,…,N1}的基函数,则对uh∈Uh,可表示为:(Γ0h是Γ0的近似)
设节点基函数为ψj=(ψ1j,ψ2j)T。则得到有限体积元法方程是:
其中,Γ1h是Γ1的近似。
其中
总之可列出2N1个形如(9)式的有限体积元方程和Γ0h上的位移边界条件。也就是说,假设已知给定的表面力,和正数λ、μLam é常数,就可以计算出位移所对应的数值。
步骤S150,通过求解得到的位移对包含了位移对未知弹性模量的导数的方程进行求解得到位移对未知弹性模量的导数,位移对未知弹性模量的导数是与位移所在的子区域相对应的。
本实施例中,在式(4)中对每一子区域的未知弹性模量求导数,即可得到慎勿使形谍了位移对未知弹性模量的导数的方程,即:
其中z用于子区域z,(13)式的形式如同(4),这里是关于u′j的方程,方程右端可看为是(4)式中的f,因此可以用解(4)式的有限体积元法解(13)式,得到u′j的解,在此不再赘述。
步骤S30,通过有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值,直至弹性模量值对应的最优化问题平方差小于预设的容忍误差,且达到预设的迭代次数时停止牛顿迭代。
其中ΔEz是弹性模量Ez的改进值,它是如下正则化方程组的解
因此,由上述推导可知,在进行牛顿迭代计算弹性模量值时,将通过有限体积元算法求解得到的位移和导数经由式(15)和式(16)计算得到改进值,进而方可进行下一步的迭代。
如图5所示,在一个实施例中,上述通过有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值的步骤包括:
步骤S310,通过每一子区域对应的位移和导数进行运算得到子区域对应的弹性模量改进值。
步骤S330,根据弹性模量改进值以设定的弹性模量初值为起始进行牛顿迭代得到与当前迭代次数对应的弹性模量值。
本实施例中,通过计算得到的弹性模量改进值和上一次迭代所得到的弹性模量值通过式(14)进行计算得到当前迭代次数所对应的弹性模量值。
如图6所示,在一个实施例中,上述步骤S330之后还包括如下步骤:
步骤S301,获取通过对成像组织进行磁共振成像所得到的位移图。
步骤S302,根据当前迭代得到的弹性模量值和位移图得到最优化问题平方差。
本实施例中,由于磁共振成像中的弹性模量重建问题被归结为一个带约束条件的最优化问题,其目标函数即为求解得到磁共振成像所测得的位移与计算得到的位移之间的最小平方和差。
该最优化问题为:min F(E)
其中
这里和分别是在成像组织位置i处通过磁共振测得的x方向和y方向的位移,和在成像组织位置i处用有限体积元法计算出的位移,一共有N个不同的位置。E是一个M维的弹性参数向量,它由一个连续基底φ张成,用它来定义整个感兴趣区域的弹性组织分布。
如图7所示,假设整个感兴趣区域为Ω,将Ω分为若干个子区域Ωsub,可将最优化问题重写为子区域上的最优化问题。假设一共有Q个子区域,有
其中,Fz(Ez)是第z个子区域上的最优化函数,对所有的子区域求和可得最优化问题为:
这里
这里和分别是在成像组织第z个子区域上iz处磁共振测得的x方向和y方向的位移,和在成像组织z个子区域上iz处用有限体积元法计算出的位移,一共有Nz个不同的位置。Ez是Mz的向量,有Nz<N,Mz<M。对于在方程两边分别对Ez求一阶导数,并令导数等于0,可得如下非线性方程组:
进而通过如上所述的牛顿迭代法求解上述非线性方程组。
步骤S303,判断最优化问题平方差是否小于预设的容忍误差,若是,则进入步骤S304,若否,则返回步骤S310。
本实施例中,预设的容忍误差是根据实际需要进行灵活设定的。
步骤S304,进一步判断平面域中的每一单元是否均包含于至少一个子区域中,若是,则进入步骤S305,若否,则返回步骤S110。
步骤S305,判断单元中的最小单元对应的迭代次数是否达至预设的迭代次数,若是,则进入步骤S306,若否,则返回步骤S10。
步骤S306,停止进行牛顿迭代。
步骤S50,将最终迭代得到的弹性模量值组成成像组织的弹性模量分布。
如图8所示,在一个实施例中,一种磁共振弹性成像中的弹性模量重建系统,包括有限体积元运算模块10、迭代模块30和分布形成模块50。
有限体积元运算模块10,用于将成像组织的表面假设为平面域,通过有限体积元算法求解得到平面域中的位移初值,以及平面域中子区域的位移、位移对未知弹性模量的导数。
本实施例中,在一定条件下,可将一个在外力作用下处于平衡状态的弹性体视为平面弹性问题,而在磁共振弹性成像中,成像组织被施加的外力是非常小的,可将该成像组织视为处于平衡状态的,因此,有限体积元运算模块10可将磁共振弹性成像中应力和应变之间的关系视为平面弹性问题。
有限体积元运算模块10将成像组织的表面视为平面域Ω,获取通过平面域对应的平衡方程所构建得到的二阶椭圆微分方程组,并在构建的二阶椭圆微分方程组中对未知弹性模量求导得到包含了位移对未知弹性模量的导数的方程。
进一步的,有限体积元运算模块10对于平面域Ω,其边界为描述成像组织平衡的状态变量包括三组,即,应力张量σ=(σ11,σ22,σ12)T、应变张量ε=(ε11,ε22,ε12)T和位移向量u=(u1,u2)T。假设Ω是均匀各向同性的弹性体,令
其中,正数λ和μ是Lam é常数:
其中,ν是Poisson系数(Poisson’s Ratio,是已知的常数),E是弹性模量(Young’sModulus)。σ,ε,u满足如下三组方程:
σ=Aε(应力-应变关系)(3)
其中,f是体积力。
由Green公式可推出:
在求解时,利用公式(1)-(3)消去部分未知量,留下的是待求量,即弹性模量E,得到包含了位移向量u的二阶椭圆微分方程组:
其中,以分别乘(4)式两端,关于x∈Ω积分,然后利用Green公式,得积分形式:
其中
在边界Γ1上有
故方程(4)的变分形式是:求u∈(H1(Ω))2,使得
其中,H1(Ω)是平面域Ω上的Hilbert空间.
进一步的,为求解出位移对未知弹性模量的导数,有限体积元运算模块10将对(4)式两端对未知弹性模量求导,以得到包含了位移对未知弹性模量的导数的方程。
在构建得到了二阶椭圆微分方程组和包含了位移对未知弹性模量的导数的方程之后,有限体积元运算模块10将应用有限体积元算法求解位移和位移对未知弹性模量的导数。
有限体积元法是吸取了有限差分法和有限元法的优点发展起来的,又称广义差分法,可以处理复杂的边值条件及不规则区域,并且由于离散化得到的线性方程组是稀疏的,在不降低数值解的收敛阶的前提下求解方程的计算量小,计算速度快,因此,将有限体积元法应用于磁共振弹性成像中将大大地降低了弹性模量重建的计算量。
在一个实施例中,上述有限体积元运算模块10还用于根据设定的弹性模量初值通过有限体积元算法计算得到平面域中的位移初值。
本实施例中,有限体积元运算模块10将首先进行位移初值的求解,以方便后续的位移求解过程。在位移初值的求解过程中,将进行弹性模量初值的假设,并根据这一假设通过有限体积元算法得到平面域所对应的位移初值,其中,该位移初值计算时所使用的边界Γ上的初值取为对成像组织进行磁共振成像所得到的位移图中对应的位移值。
如图9所示,在一个实施例中,上述有限体积元运算模块10包括划分单元110、子区域位移求解单元130和子区域导数求解单元150。
划分单元110,用于将平面域划分为若干个子区域,并划分子区域为若干个单元,该子区域之间和单元之间不存在重叠且任一单元的顶点均不在其它单元的边上,平面域边界的顶点为单元的顶点。
本实施例中,划分单元110将平面域Ω划分为若干个子区域,然后分别划分每一子区域为若干个单元之和,在优选的实施例中,划分单元110划分得到的单元为三角形的形状,进而使得子区域之间、单元之间互相不存在重叠区域,并且任一单元的顶点不会在其它单元的边上,边界Γ的每一顶点均为单元的顶点,从而每一子区域将对应了一个原始剖分Th,其中h是所有三角形单元边的最长边。记KQ为Th内的三角形单元,Q是三角形单元的外心。
进一步的,划分单元110记是原始剖分Th的对偶剖分,包括是重心对偶剖分和外心对偶剖分。以下设是外心对偶剖分。外心对偶剖分做法:设Th的任一单元的内角不大于90°,取三角单元□P0PiPi+1(i=1,2,…6,P7=P1)的外心Qi为对偶剖分的节点,即依次连结各个外心点可得到外心对偶剖分,此时,是的中垂线,分别过边中点Mi。
子区域位移求解单元130,用于在子区域中,以构成单元的顶点作为子区域的节点,并构建节点对应的有限体积元方程,并以位移初值作为有限体积元方程中的初值进行求解得到子区域对应的位移。
本实施例中,对于每一子区域,子区域位移求解单元130设是定义在原始分剖Th上的试探函数空间,它是分片线性函数空间,即
子区域位移求解单元130将Th的内节点编号为1,2,…,N0。其中,节点分两类,给定力条件的节点编号为N0+1,…,N1,给定位移条件的节点编号为N1+1,…,N。用表示节点i∈{1,2,…,N1}的基函数,则对uh∈Uh,可表示为:(Γ0h是Γ0的近似)
设节点基函数为ψj=(ψ1j,ψ2j)T。则得到有限体积元法方程是:
求uh∈Uh,使得
其中,Γ1h是Γ1的近似。
将式(4)两端在上积分,应用Green公式并用uh代替u,得
其中
其中,(x1(Pi),x2(Pi))是点Pi的坐标,是含外心Q1的单元面积。其余折线段上的积分类推。
总之子区域位移求解单元130可列出2N1个形如(9)式的有限体积元方程和Γ0h上的位移边界条件。也就是说,假设已知给定的表面力,和正数λ、μLamé常数,就可以计算出位移所对应的数值。
子区域导数求解单元150,用于通过求解得到的位移对包含了位移对未知弹性模量的导数的方程进行求解得到位移对未知弹性模量的导数,位移对未知弹性模量的导数是与位移所在的子区域相对应的。
迭代模块30,用于通过有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值,直至弹性模量值对应的最优化问题平方差小于预设的容忍误差,且达到预设的迭代次数时停止牛顿迭代。
本实施例中,迭代模块30应用牛顿迭代法进行弹性模量值的,假设迭代第n步的解为则第n+1步的解为
其中ΔEz是弹性模量Ez的改进值,它是如下正则化方程组的解
因此,由上述推导可知,迭代模块30在进行牛顿迭代计算弹性模量值时,将通过有限体积元算法求解得到的位移和导数经由式(15)和式(16)计算得到改进值,进而方可进行下一步的迭代。
如图10所示,在一个实施例中,上述迭代模块30包括改进值运算单元310和弹性模量迭代单元330。
改进值运算单元310,用于通过每一子区域所对应的位移和导数进行运算得到子区域对应的弹性模量改进值。
弹性模量迭代单元330,用于根据弹性模量改进值以设定的弹性模量初值为起始进行牛顿迭代得到与当前失效次数对应的弹性模量值。
本实施例中,弹性模量迭代单元330通过计算得到的弹性模量改进值和上一次迭代所得到的弹性模量值通过式(14)进行计算得到当前迭代次数所对应的弹性模量值。
如图11所示,在一个实施例中,上述迭代模块30还包括位移图像处理单元301和判断单元303。
位移图处理单元301,用于获取通过对成像组织进行磁共振成像所得到的位移图,并根据当前迭代得到的弹性模量值和位移图得到最优化问题方差。
本实施例中,由于磁共振成像中的弹性模量重建问题被归结为一个带约束条件的最优化问题,其目标函数即为求解得到磁共振成像所测得的位移与计算得到的位移之间的最小平方和差。
该最优化问题为:min F(E)
其中
这里和分别是在成像组织位置i处通过磁共振测得的x方向和y方向的位移,和在成像组织位置i处用有限体积元法计算出的位移,一共有N个不同的位置。E是一个M维的弹性参数向量,它由一个连续基底φ张成,用它来定义整个感兴趣区域的弹性组织分布。
位移图处理单元301假设整个感兴趣区域为Ω,将Ω分为若干个子区域Ωsub,可将最优化问题重写为子区域上的最优化问题。假设一共有Q个子区域,有
其中,Fz(Ez)是第z个子区域上的最优化函数,对所有的子区域求和可得最优化问题为:
这里
这里和分别是在成像组织第z个子区域上iz处磁共振测得的x方向和y方向的位移,和在成像组织z个子区域上iz处用有限体积元法计算出的位移,一共有Nz个不同的位置。Ez是Mz的向量,有Nz<N,Mz<M。对于在方程两边分别对Ez求一阶导数,并令导数等于0,可得如下非线性方程组:
进而位移图处理单元301通过如上所述的牛顿迭代法求解上述非线性方程组。
判断单元303,用于判断最优化问题平方差是否小于预设的容忍误差,若是,则进一步判断单元中的最小单元对应的迭代次数是否达到预设的迭代次数,若是,则停止进行牛顿迭代。
本实施例中,预设的容忍误差是根据实际需要进行灵活设定的。
在另一个实施例中,上述判断单元303还用于若判断到单元中的最小单元对应的迭代次数未达到预设的迭代次数,则通知上述有限体积元运算模块10。
分布形成模块50,用于将最终迭代得到的弹性模量值组成成像组织的弹性模量分布。
上述磁共振弹性成像中的弹性模量重建方法和系统,将成像组织的表面假设为平面域,引入有限体积元算法以求解得到位移和该位移对未知弹性模量的导数,进而通过求解得到的位移和导数进行牛顿迭代得到弹性模量值,并在弹性模量值对应的最优化问题平方差小于预设的容忍误差,且达到预设的迭代次数时停止牛顿迭代,进而将最终迭代得到的弹性模量值组成成像组的弹性模量分布,通过引入有限体积元算法,使得弹性模量的重建在保证了计算精度的前提下降低了计算量,提高了计算速度。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程,是可以通过计算机程序来指令相关的硬件来完成,所述的程序可存储于一计算机可读取存储介质中,该程序在执行时,可包括如上述各方法的实施例的流程。其中,所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory,ROM)或随机存储记忆体(Random Access Memory,RAM)等。
以上所述实施例仅表达了本发明的几种实施方式,其描述较为具体和详细,但并不能因此而理解为对本发明专利范围的限制。应当指出的是,对于本领域的普通技术人员来说,在不脱离本发明构思的前提下,还可以做出若干变形和改进,这些都属于本发明的保护范围。因此,本发明专利的保护范围应以所附权利要求为准。
Claims (12)
1.一种磁共振弹性成像中的弹性模量重建方法,包括如下步骤:
将成像组织的表面假设为平面域,通过有限体积元算法求解得到所述平面域中的位移初值,以及所述平面域中子区域的位移、所述位移对未知弹性模量的导数;
通过所述有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值,直至所述弹性模量值对应的最优化问题平方差小于预设的容忍误差,且达到预设的迭代次数时停止牛顿迭代;
将所述最终迭代得到的弹性模量值组成所述成像组织的弹性模量分布。
2.根据权利要求1所述的磁共振弹性成像中的弹性模量重建方法,其特征在于,所述通过有限体积元算法求解得到所述平面域中的位移初值的步骤为:
根据设定的弹性模量初值通过有限体积元算法计算得到平面域中的位移初值。
3.根据权利要求1所述的磁共振弹性成像中的弹性模量重建方法,其特征在于,所述通过有限体积元算法求解得到所述平面域中子区域的位移、所述位移对未知弹性模量的导数的步骤包括:
将所述平面域划分为若干个子区域,并划分所述子区域为若干个单元,所述子区域之间和单元之间不存在重叠且任一单元的顶点均不在其它单元的边上,平面域边界的顶点为单元的顶点;
在所述子区域中,以构成单元的顶点作为所述子区域的节点,并构建所述节点对应的有限体积元方程,并以所述位移初值作为有限体积元方程中的初值进行求解得到所述子区域对应的位移;
通过所述求解得到的位移对包含了所述位移对未知弹性模量的导数的方程进行求解得到所述位移对未知弹性模量的导数,所述位移对未知弹性模量的导数是与位移所在的子区域相对应的。
4.根据权利要求3所述的磁共振弹性成像中的弹性模量重建方法,其特征在于,所述通过所述有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值的步骤包括:
通过每一子区域所对应的位移和导数进行运算得到所述子区域对应的弹性模量改进值;
根据所述弹性模量改进值以设定的弹性模量初值为起始进行牛顿迭代得到与当前迭代次数对应的弹性模量值。
5.根据权利要求4所述的磁共振弹性成像中的弹性模量重建方法,其特征在于,所述根据所述弹性模量改进值以设定的弹性模量初值为起始进行牛顿迭代得到与当前迭代次数对应的弹性模量值的步骤之后还包括:
获取通过对所述成像组织进行磁共振成像所得到的位移图;
根据当前迭代得到的弹性模量值和所述位移图得到最优化问题平方差;
判断所述最优化问题平方差是否小于预设的容忍误差,若是,则进一步判断所述平面域中的每一单元是否均包含于至少一个子区域中,若是,则
判断所述单元中的最小单元对应的迭代次数是否达至预设的迭代次数,若是,则
停止进行牛顿迭代。
6.根据权利要求5所述的磁共振弹性成像中的弹性模量重建方法,其特征在于,还包括:
若判断到所述单元中的最小单元对应的迭代次数未达到预设的迭代次数,则
返回所述通过有限体积元算法求解得到所述平面域中的位移初值,以及所述平面域中子区域的位移、所述位移对未知弹性模量的导数的步骤。
7.一种磁共振弹性成像中的弹性模量重建系统,其特征在于,包括:
有限体积元运算模块,用于将成像组织的表面假设为平面域,通过有限体积元算法求解得到所述平面域中的位移初值,以及所述平面域中子区域的位移、所述位移对未知弹性模量的导数;
迭代模块,用于通过所述有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值,直至所述弹性模量值对应的最优化问题平方差小于预设的容忍误差,且达到预设的迭代次数时停止牛顿迭代;
分布形成模块,用于将所述最终迭代得到的弹性模量值组成所述成像组织的弹性模量分布。
8.根据权利要求7所述的磁共振弹性成像中的弹性模量重建系统,其特征在于,所述有限体积元运算模块还用于根据设定的弹性模量初值通过有限体积元算法计算得到平面域中的位移初值。
9.根据权利要求7所述的磁共振弹性成像中的弹性模量重建系统,其特征在于,所述有限体积元运算模块包括:
划分单元,用于将所述平面域划分为若干个子区域,并划分所述子区域为若干个单元,所述子区域之间和单元之间不存在重叠且任一单元的顶点均不在其它单元的边上,平面域边界的顶点为单元的顶点;
子区域位移求解单元,用于在所述子区域中,以构成单元的顶点作为所述子区域的节点,并构建所述节点对应的有限体积元方程,并以所述位移初值作为有限体积元方程中的初值进行求解得到所述子区域对应的位移;
子区域导数求解单元,用于通过所述求解得到的位移对包含了所述位移对未知弹性模量的导数的方程进行求解得到所述位移对未知弹性模量的导数,所述位移对未知弹性模量的导数是与位移所在的子区域相对应的。
10.根据权利要求9所述的磁共振弹性成像中的弹性模量重建系统,其特征在于,所述迭代模块包括:
改进值运算单元,用于通过每一子区域所对应的位移和导数进行运算得到所述子区域对应的弹性模量改进值;
弹性模量迭代单元,用于根据所述弹性模量改进值以设定的弹性模量初值为起始进行牛顿迭代得到与当前迭代次数对应的弹性模量值。
11.根据权利要求10所述的磁共振弹性成像中的弹性模量重建系统,其特征在于,所述迭代模块还包括:
位移图处理单元,用于获取通过对所述成像组织进行磁共振成像所得到的位移图,并根据当前迭代得到的弹性模量值和所述位移图得到最优化问题平方差;
判断单元,用于判断所述最优化问题平方差是否小于预设的容忍误差,若是,则进一步判断所述平面域中的每一单元是否均包含于至少一个子区域中,若是,则判断所述单元中的最小单元对应的迭代次数是否达到预设的迭代次数,若是,则停止进行牛顿迭代。
12.根据权利要求11所述的磁共振弹性成像中的弹性模量重建系统,其特征在于,所述判断单元还用于若判断到所述单元中的最小单元对应的迭代次数未达到预设的迭代次数,则通知所述有限体积元运算模块。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201210572449.6A CN103049663B (zh) | 2012-12-25 | 2012-12-25 | 磁共振弹性成像中的弹性模量重建方法和系统 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201210572449.6A CN103049663B (zh) | 2012-12-25 | 2012-12-25 | 磁共振弹性成像中的弹性模量重建方法和系统 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN103049663A true CN103049663A (zh) | 2013-04-17 |
CN103049663B CN103049663B (zh) | 2016-06-08 |
Family
ID=48062297
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201210572449.6A Active CN103049663B (zh) | 2012-12-25 | 2012-12-25 | 磁共振弹性成像中的弹性模量重建方法和系统 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN103049663B (zh) |
Cited By (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102920457A (zh) * | 2011-12-12 | 2013-02-13 | 中国科学院深圳先进技术研究院 | 磁共振弹性成像精确度检测方法 |
CN105068785A (zh) * | 2015-04-22 | 2015-11-18 | 清华大学 | 一种并行计算方法及系统 |
CN110916662A (zh) * | 2019-12-05 | 2020-03-27 | 无锡鸣石峻致医疗科技有限公司 | 一种便携式核磁共振器官弹性无创定量检测系统 |
CN112327233A (zh) * | 2020-11-02 | 2021-02-05 | 上海交通大学 | 多相位快速磁共振弹性成像采集与重建方法及系统 |
CN116026682A (zh) * | 2023-03-30 | 2023-04-28 | 浙江大学 | 基于qme的快速弹性成像计算方法 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2012145312A2 (en) * | 2011-04-22 | 2012-10-26 | Mayo Foundation For Medical Education And Research | System and method for magnetic resonance elastography of the breast |
CN102764141A (zh) * | 2012-07-20 | 2012-11-07 | 中国科学院深圳先进技术研究院 | 弹性成像方法和系统及其中的生物组织位移估计方法和系统 |
-
2012
- 2012-12-25 CN CN201210572449.6A patent/CN103049663B/zh active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2012145312A2 (en) * | 2011-04-22 | 2012-10-26 | Mayo Foundation For Medical Education And Research | System and method for magnetic resonance elastography of the breast |
CN102764141A (zh) * | 2012-07-20 | 2012-11-07 | 中国科学院深圳先进技术研究院 | 弹性成像方法和系统及其中的生物组织位移估计方法和系统 |
Non-Patent Citations (1)
Title |
---|
ELIJAH E.W. VAN HOUTEN等: "Three-Dimensional Subzone-Based Reconstruction Algorithm for MR Elastography", 《MAGNETIC RESONANCE IN MEDICINE》 * |
Cited By (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102920457A (zh) * | 2011-12-12 | 2013-02-13 | 中国科学院深圳先进技术研究院 | 磁共振弹性成像精确度检测方法 |
CN102920457B (zh) * | 2011-12-12 | 2015-03-11 | 中国科学院深圳先进技术研究院 | 磁共振弹性成像精确度检测方法 |
CN105068785A (zh) * | 2015-04-22 | 2015-11-18 | 清华大学 | 一种并行计算方法及系统 |
CN105068785B (zh) * | 2015-04-22 | 2018-04-10 | 清华大学 | 一种并行计算方法及系统 |
CN110916662A (zh) * | 2019-12-05 | 2020-03-27 | 无锡鸣石峻致医疗科技有限公司 | 一种便携式核磁共振器官弹性无创定量检测系统 |
CN112327233A (zh) * | 2020-11-02 | 2021-02-05 | 上海交通大学 | 多相位快速磁共振弹性成像采集与重建方法及系统 |
CN116026682A (zh) * | 2023-03-30 | 2023-04-28 | 浙江大学 | 基于qme的快速弹性成像计算方法 |
Also Published As
Publication number | Publication date |
---|---|
CN103049663B (zh) | 2016-06-08 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Mohamed et al. | Deformable registration of brain tumor images via a statistical model of tumor-induced deformation | |
Barnhill et al. | Nonlinear multiscale regularisation in MR elastography: Towards fine feature mapping | |
CN103049663A (zh) | 磁共振弹性成像中的弹性模量重建方法和系统 | |
Buttinger-Kreuzhuber et al. | A fast second-order shallow water scheme on two-dimensional structured grids over abrupt topography | |
Jeong et al. | Interactive visualization of volumetric white matter connectivity in DT-MRI using a parallel-hardware Hamilton-Jacobi solver | |
McGarry et al. | A heterogenous, time harmonic, nearly incompressible transverse isotropic finite element brain simulation platform for MR elastography | |
Lombardi | Inverse problems in 1D hemodynamics on systemic networks: A sequential approach | |
Kamali et al. | Elasticity imaging using physics-informed neural networks: Spatial discovery of elastic modulus and Poisson's ratio | |
Karaçali et al. | Simulation of tissue atrophy using a topology preserving transformation model | |
Tan et al. | Gradient-based optimization for poroelastic and viscoelastic MR elastography | |
Babaniyi et al. | Recovering vector displacement estimates in quasistatic elastography using sparse relaxation of the momentum equation | |
Wang et al. | A novel cortical thickness estimation method based on volumetric Laplace–Beltrami operator and heat kernel | |
Miller et al. | Estimation of transversely isotropic material properties from magnetic resonance elastography using the optimised virtual fields method | |
Zhao et al. | Cerebral vascular strains in dynamic head impact using an upgraded model with brain material property heterogeneity | |
Tripura et al. | A wavelet neural operator based elastography for localization and quantification of tumors | |
Nguyen et al. | A Fourier‐series‐based virtual fields method for the identification of three‐dimensional stiffness distributions and its application to incompressible materials | |
Mohammed et al. | Model-based quantitative elasticity reconstruction using ADMM | |
Vigneron et al. | On extended finite element method (xfem) for modelling of organ deformations associated with surgical cuts | |
Feng et al. | An adjoint-based method for a linear mechanically-coupled tumor model: application to estimate the spatial variation of murine glioma growth based on diffusion weighted magnetic resonance imaging | |
Galarce et al. | Displacement and pressure reconstruction from magnetic resonance elastography images: Application to an in silico brain model | |
Rodríguez-Arós et al. | Numerical analysis of a frictional contact problem for viscoelastic materials with long-term memory | |
Iffrig et al. | A new method for quantifying abdominal aortic wall shear stress using phase contrast magnetic resonance imaging and the Womersley solution | |
Hajhashemkhani et al. | A novel method for the identification of the unloaded configuration of a deformed hyperelastic body | |
CN102920457B (zh) | 磁共振弹性成像精确度检测方法 | |
EP3438989A1 (en) | Method and apparatus for predicting fluid flow through a subject conduit |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant |