专利名称：：基于四维医学图像的心壁应力应变测量方法
技术领域：
：本发明涉及图像处理、计算机视觉、医学图像分析、心血管医学、机械动力学和运动学领域的一种心壁应力应变测量方法。
背景技术：
：心脏是人体血液的动力中枢，因此心脏功能的分析研究也备受关注，而左心室在整个心动过程中提供主要的动能，因此，对心脏功能的分析研究可以转化为对左心室的分析研究。左心室的功能参数有很多，比如全局参数左心室容积(LVV)，射血份数(EF)，心室壁厚(WT)等，这些参数可以在一定程度上反映整个心室的功能强弱，但却无法定位具体的病变部位，因而近年来学界加大了对局部参数的研究力度，其中最主要的局部参数是应力与应变。应变与应力是指柔性体在受到一定外力作用的情况下发生的形变以及由这种形变引起的抵御外力的自身作用力，应力与应变可以很好地反应柔性体在不同部位所受到的力以其形变情况。就左心室而言，如果某一部位发生了病变从而导致其心室壁的性质发生变化，必然会体现在应力应变图上。在医学研究方面，也有很多文献开始以应力应变作为指标研究不同的心脏疾病，有的用左心室舒张期应力的变化率来评价心室的整体功能，也有通过对比心力衰竭患者与正常人左心室应变值研究分析了应变与心脏功能分析的关系，或者研究了冠心病与高血压患者在应力应变参数上的变化。目前对应力与应变的计算主要有两种方式，即从应力到应变和从应变到应力。前一种方式通常以心室运动的电位变化为依据建立电位模型，由模型得到心室的应力，再利用有限元方法计算对应位置的应变。后一种方式通常是由心室壁的位移得到应变，再由应变计算出应力。本文采用的是从应变计算应力的方式，因为该方式可以以心室图片为数据源，无需接触，且对于不同个体的适应性较好。由于应力应变之间的关系遵循一定力学对应关系(弹性力学已经给出)，因而该方式的关键在于如何获得应变值，而应变是位移相对于位置的导数，所以求取应变的关键在于如何求解左心室壁各处的位移。通过各种成像和标记方法，我们可以很方便地得到心动周期不同时刻的左心室模型，要得到位移就需要在不同时刻的模型上找到同一点的位置。目前的主流方法是将不同时刻模型分别划分为有限个单元，在模型之间寻找最相似的单元从而进行位置匹配最终计算出位移和应变。这种通过匹配单元计算位移的方法有其无法克服的问题，首先它的前提是前后时刻的形变是微小的，形状是大致相同的，事实上左心室壁的位移并不严格满足以上假设；其次是它虽然能计算单元结点处的位移，但对于单元内部的点依然要采用近似插值计算，不够准确；第三是为了保证匹配精度每次计算必须对整个左心室室壁的所有单元都进行计算，从而剧增了运算量。
发明内容为了克服已有的心壁应力应变测量方法的形变匹配性差、精度低、计算复杂的不5足，本发明提供一种形变匹配性良好、精度高、大大减少计算复杂度的基于四维医学图像的心壁应力应变测量方法。本发明解决其技术问题所采用的技术方案是一种基于四维医学图像的心壁应力应变测量方法，所述心壁应力应变测量方法包括以下步骤1)、将待检测的心脏建立B样条曲面模型P(U，ν)，其中U，ν是B样条曲面定义时的规范化参数，参数域均为W，l]；设模型上存在η个位移已知的标记点(i=1,2,...η),Bi的在球面坐标系下的空间位置为[QiUTiRi],在P(u，v)上的规范化坐标为[ui;Vi]，且Bi的位移为[rθirVirRj，通过离散的η个点的位移矢量建立整个模型的rθ、rψ和rR的连续位移场；2)、利用规范化坐标U，ν拟合位移场并将位移场中的点与实际的空间点建立对应关系，对三个位移分别拟合，其中Γθ位移的拟合过程使用[UiVirej拟合出了场曲面Q(s，t)，根据B样条曲面特点Q(s，t)实际包括3个方程，即Qu(s，t)=u，Qv(s，t)=v,Qe(s,t)=γΘ；已知[θ￥R]=P(u,ν),[uνrθ]=Q(s，t)，定义P—1和Q-1分别为P(u，ν)与Q(s，t)的逆运算；对模型上任意点aje^RJ，计算其坐标值θ的位移量rθ^的值，过程如下Stepl将坐标转换至规范坐标系上，[U0ν0]=P—1(θQψ0R0)；St印2计算该点在位移场中坐标s0,t0,[s0t0]=Q-1Oi0v0)；St印3得到位移rθ，rθ=Q0(s0,t0)；rΨ位移的拟合过程使用[UiVirVi]拟合出了场曲面Q(s，t)，且有Qu(s，t)=u，Qv(s，t)=ν,Qv(s,t)=ΓΨ；[θψR]=P(u，ν),[uνγψ]=Q(s，t)；设计算点a0[θ0Ψ。RJ的rΨ位移值，过程如下Stepl将坐标转换至规范坐标系上，[U0ν0]=Ρ^(ψ0V0R0)；St印2计算该点在位移场中坐标s0,t0,[s0t0]=Q-1Oi0v0)；St印3得到位移rΨ，rΨ=Qv(s0,t0)；rR位移的拟合过程使用[UiVirRj拟合出了场曲面Q(s，t)，且有Qu(s，t)=u,Qv(s，t)=v,QE(s,t)=rR；[θψR]]=P(u,ν),[uνrR]=Q(s，t)；设计算点a0[eoΨ0RJ的rR位移值，过程如下Stepl将坐标转换至规范坐标系上，[u0V0]=P-1(R0Ψ0R0)；St印2计算该点在位移场中坐标s0,t0,[s0t0]=Q-1Oi0v0)；St印3得到位移rR,rR=Qe(s0,t0)；3)、将步骤2)中心壁上任意点aje。ψ0R0]的位移量分别记做N0(θ。ψ0R0)、Νν(θ。ψ。Rq)和炉(9。ψ。Rq)；由于应变是在笛卡尔坐标系x，y，ζ下定义的，我们需将该定义转化至球面坐标系下，设％在笛卡尔系下坐标为[XtlY0、!，其位移量分别为化丨^，^，z0),My(x0,y0,z0),Mz(x0,y0,Z(1)，则根据公式⑶可直接由球面坐标系下位移求出应变，并根据(6)求出应力。(8)笛卡尔坐标系下的应变定义物体位移向量δ由x，y，z轴上的位移u，v，w组成t应变向量ε有3个主应变和三个剪应变，与位移的关系为(6)其中，D为弹性矩阵，由研究对象的弹性模量E及泊松比μ确定，对于确定的材料可视为已知常量。作为优选的一种方案在所述步骤1)中，B样条曲面模型建立过程由一组控制点di(i=0,1,2,...,η)生成的k次B样条曲线定义为其中，Niik(U)为B样条基函数，其经典的de-BoorCox定义式为一组递推公式，即式中Ui为一组非递减规范化序列，称作节点序列，通常取i=0，1，2，...，n+k+1同理由(m+1)X(n+1)个控制点生成的三维B样条曲面定义为(3)其中，kl、k2分别为u和ν方向上的曲线阶数，样条曲面就是样条的推广，它将-组空间坐标下的点转换到了一个IX1的空间U-V内，将U-V空间定义为规范化空间。本发明的有益效果主要表现在形变匹配性良好、精度高、大大减少计算复杂度。图1是NURBS曲面逆运算流程图。图2是应力应变定义及壳状单元的示意图。图3是位移场与空间位置对应关系的示意图。图4是左心室7个NURBS模型的示意图。图5是左心室内壁第一相位标记点的位移矢量图。图6是左心室内壁第一相位标记点的对比图。图7是左心室内壁第一相位χ方向位移场曲面图。图8是左心室内壁第一相位y方向位移场曲面图。图9是左心室内壁第一相位ζ方向位移场曲面图。图10是左心室内壁在7个相位处的应力场分布色图。具体实施例方式下面结合附图对本发明作进一步描述。参照图1图10，一种基于四维医学图像的心壁应力应变测量方法，所述心壁应力应变测量方法包括以下步骤1)、将待检测的心脏建立B样条曲面模型P(u，ν)，其中U，ν是B样条曲面定义时的规范化参数，参数域均为W，l]；设模型上存在η个位移已知的标记点(i=1,2,...η),Bi的在球面坐标系下的空间位置为[QiUTiRi]，在P(u，ν)上的规范化坐标为[Ui,Vi]，且Bi的位移为[rθirVirRj，通过离散的η个点的位移矢量建立整个模型的rθ、rψ禾口rR的连续位移场；2)、利用规范化坐标U，ν拟合位移场并将位移场中的点与实际的空间点建立对应关系，对三个位移分别拟合，其中ΓΘ位移的拟合过程使用[UiViΓθ,]拟合出了场曲面Q(s，t)，根据B样条曲面特点Q(s，t)实际包括3个方程，即Qu(s，t)=u,Qv(s，t)=ν,Q0(s,t)=γΘ；已知[ΘvR]=P(u,v),[uνΓθ]=Q(s，t)，定义P—1和Q-1分别为P(u，ν)与Q(s，t)的逆运算；对模型上任意点a0[θ^RJ，计算其坐标值θ的位移量I^tl的值，过程如下Stepl将坐标转换至规范坐标系上，[U0ν0]=P—1(θQψ0R0)；St印2计算该点在位移场中坐标s0,t0,[s0t0]=Q-1Oi0v0)；0056]St印3得到位移rθ，rθ=Q0(s0,t0)；ru/位移的拟合过程使用[UiViri^]拟合出了场曲面Q(s，t)，且有Qu(s，t)=u，Qv(s，t)=ν,Qv(s,t)=γψ；[θψR]=P(u，ν),[uνγψ]=Q(s，t)；设计算点a0[θ0Ψ。RJ的rΨ位移值，过程如下Stepl将坐标转换至规范坐标系上，[U0ν0]=Ρ—Ηψ。Ψ。R0)；St印2计算该点在位移场中坐标s0,t0,[s0t0]=Q-1Oi0v0)；St印3得到位移rΨ，rΨ=Qv(s0,t0)；rR位移的拟合过程使用[UiVirRj拟合出了场曲面Q(s，t)，且有Qu(s，t)=u,Qv(s，t)=v,QE(s,t)=rR；[θψR]]=P(u,ν),[uνrR]=Q(s，t)；设计算点a0[eoΨ0RJ的rR位移值，过程如下Stepl将坐标转换至规范坐标系上，[u0V0]=P-1(R0Ψ0R0)；St印2计算该点在位移场中坐标s0,t0,[s0t0]=Q-1Oi0v0)；St印3得到位移rR,rR=Qe(s0,t0)；3)、将步骤2)中心壁上任意点&0[90ψ0RJ的位移量分别记做N0(θ0ψ0R0)、Νν(θ。ψ。Rq)和炉(9。ψ。Rq)；由于应变是在笛卡尔坐标系x，y，ζ下定义的，我们需将该定义转化至球面坐标系下，设％在笛卡尔系下坐标为[XtlY0、!，其位移量分别为化丨^，^，z0),My(x0,y0,z0),Mz(x0,y0,Z(1)，则根据公式⑶可直接由球面坐标系下位移求出应变，并根据(6)求出应力。(8)笛卡尔坐标系下的应变定义物体位移向量δ由X，y，ζ轴上的位移u，v，w组成。应变向量ε有3个主应变和三个剪应变，与位移的关系为应力向量ο定义如下{σ}=[σχσyσζτxyτyzτχζ]=Dε(6)其中，D为弹性矩阵，由研究对象的弹性模量E及泊松比μ确定，对于确定的材料可视为已知常量。本实施例的模型的建立过程为B样条方法及其逆运算由一组控制点di(i=0,1,2,...,η)生成的k次B样条曲线定义为其中，Niik(U)为B样条基函数，其经典的de-BoorCox定义式为一组递推公式，即式中Ui为一组非递减规范化序列，称作节点序列，通常取i=0，1，2，...，n+k+L同理由(m+1)X(n+1)个控制点生成的三维B样条曲面定义为其中，基函数的递推同(2)，kl、k2分别为11和￥方向上的曲线阶数。样条曲面就是样条的推广，它将一组空间坐标下的点转换到了一个边长为1的平面空间u-v内，我们将u-v空间称为规范化空间。以上的定义是B样条方法的正向计算，即已知u，ν求ρ(u，ν)，在很多时候，我们需要在建立好的B样条模型P(u，ν)上求的一点d[x，y，z]的参数u，v，我们称这样的计算为B样条曲面的逆运算，记为F1O^y,z)，可知[U，ν]ζΓΥχ,γ,ζ)。F1OC,y，ζ)实际包含3个函数，即Pf1(χ)，Py-1(y)和Ρζ—1(ζ)，它们分别是取C^j的x，y，ζ轴分量为控制点P(u，ν)的逆函数。对于逆函数的计算我们采用迭代的方法在u-v空间中进行搜索得到，迭代过程如图1。图1中ε为求解精度，m定义了逼近速度，Pd(i，j)定义如下*[P1{mu+{i+\)d,mv+{j+\)d)-x]<0Pd(i,;')<0ο][Py(mu+id,mv+jd)[Py{mu+(/+\)d,mv+(j+\)d)-少]<0*[P2(mu+(J+l)d,mv+(j+\)d)-z]<0(4)结束时的mu+id和mu+jd即为所求u，ν。弹性力学方程及单元坐标变换根据弹性力学知识在笛卡尔坐标系下定义物体位移向量δ由x，y，ζ轴上的位移u，ν，w组成。应变向量ε有3个主应变和三个剪应变，与位移的关系为其中，D为弹性矩阵，由研究对象的弹性模量E及泊松比μ确定，对于确定的材料可视为已知常量。由公式(5)和(6)可知求解应力应变的关键可归结于如何根据位移计算应变。公式(5)和(6)的定义是在笛卡尔坐标系下的，也即研究对象的局部坐标，如图2(a)所示。由于我们的研究模型为壳状，不同位置点的局部坐标和全局坐标具有一定旋转和形变关系，如图2(b)。为了克服形状干扰，我们采用球面坐标系进行计算。笛卡尔坐标系与球面坐标系转化关系如下球面坐标系对于曲面与平面的变形有很好的适应性。如图2(b)中所示，已知点B0[θOV0R0]的位移量分别记做N0(θ0V0R0)、ΝΨ(θ0Ψ0R0)和Νκ(θ0Ψ0R0)，a0在笛卡尔系下坐标为[χ。y。ζ。]，其位移量分别为Mx(X。，y。，z0)，My(x0,y。，ζ。)，Mz(χ0,y。，ζ。)，则根据公式(8)可直接由球面坐标系下位移求出应变且Yxy=Yyx,Yyz=Yzy,Yxz=Yzxo(8)根据位移场计算应变已知B样条曲面模型P(u，ν)，模型上存在η个位移已知标记点Bi(i=1,2,...n),Bi对应的规范化坐标为[ui;Vi]，Bi的位移为[rθirVirRj，我们的目标是要通过离散的η个位移建立一个连续的位移场。对于一个没有出现断裂的变形体，连续位置的形变也必然是连续的，因此对于位移场我们也通过对离散点进行B样条拟合来获得。为了将位移场中的点与运动点的位置建立联系，我们对三个位移分别拟合。以ΓΘ场为例，我们使用[UiVirej拟合出了场曲面Q(s，t)，且有Qu(s，t)=u，Qv(s，t)=ν,Qe(s，t)=rθ。各个空间下对应点的坐标关系如图3。由图3可以清楚的看到[θψR]=P(u，ν)，[uνrθ]=Q(s，t)。假设我们要计算点d0[θοΨοR0]的rθ位移值，过程如下Stepl将坐标转换至规范坐标系上，[u0V0]=P—1(θQΨ(1R。)；St印2计算该点在位移场中坐标s0)t0,[s0t0]=Q-1Oi0v0)；St印3得到位移rθ，rθ=Q0(s0,t0)。通过这种变换我们可以得到心壁上任意点的位移向量，我们将这种计算记做Trans0(θψR)，相应的还有Transv(θψR),TransE(6ψR)。通过以上计算可以得到全局坐标下Cltl的位置及对应位移，接着应用公式(9)将位置及位移旋转至局部坐标下，最后用公式(8)和(6)分别计算应变和应力。本实施例的模型及位移场本文实验的数据源选择左心室的SPECT图像，这是由于在SPECT图像中物理点和标记点很清楚的可以看到，大量的物理点能在很短的时间间隔内被标记，而且可以在心脏周期内被全程跟踪。选择左心室内外壁共1165个标记点(其中内壁355个，外壁810个)，每隔IOOms跟踪记录它们的位置，获得一个心动周期内7个相位的标记点位置，记做Di(i=1，2，...，7)。根据公式(7)，将Di转化至球坐标下记做D/，并以Di'为源应用NURBS方法拟合左心室7个时刻的内外面模型，模型如图4所示。通过做差运算可以找到标记点在相邻时刻的位移为Di+1’-D/，图5为第一相位到第二相位的位移矢量图。该矢量图直观的反应了左心室心室壁的运动方向及大小，与图6的图形进行比较发现位移的总体方向呈螺旋扭曲形，非常的相近也验证了标记点位移的正确性。根据位移及规范化坐标我们对位移进行拟合得到了各个时刻的位移场，图7、图8和图9为第一相位左心室内壁的位移场曲面图。从位移场中我们可以明显的看出在左心室运动时，心尖处的位移大于心底，而左心室后部的位移大于前部。众所周知对于心脏的研究无法以实体实验的方式进行验证，因此我们将结果与其他方法获得的结果进行比较。现有技术通过有限元方法建立模型对左心室形变进行研究得到了结论，即心室扭曲运动，心尖处绕心底旋转，心室内壁扭曲程度大于外壁，这都与我们得到的左心室位移场相吻合，从而验证了位移场的正确性。位移场同时体现出了左心室室壁的形变特点，对它进行定量分析可以的到一些标志性的数据。我们通过位移场对心动一周期内心室长轴与内壁半径进行分析，得到图7、8、9所示结果，与现有技术得到的结果进行比较，可以看到变化方向与幅度都大体相同。应力应变计算分析左心室的应变量具有十分重要的意义，因为它可以定量来研究心肌性能，它可以将心肌的局部变形分离出来。由于采用不同的标记成像技术，标记结果会有一定得差异，导致数据结果也存在一定得差异，在本文中我们将所得结果与taggedMRI技术所得得研究结果进行对比，该结果得到了广泛的认可。通过对比，我们得到的应变场与现有技术得到的结论一致；即心尖到心底应变逐渐减小，心后应变大于心前。左心室壁的应力场计算具有重要意义，心壁应力是心肌冠脉血流和耗氧情况的主要决定因素之一。根据应变场，弹性模量与泊松比分别取E=，λ=，计算心室内外壁的应力值。图10所示为左心室内壁在7个相位处的应力场分布色图，这里的应力为根据米责思约束条件计算的到得统一值，色彩的变化根据应力增大而逐渐变暖，到红色处达到最大，该结果与现有技术得到的结果相一致，即总体分布不均勻，心室内壁应力大于外壁，在心尖和心底部应力较大，心尖处还出现应力集中现象。权利要求一种基于四维医学图像的心壁应力应变测量方法，所述心壁应力应变测量方法包括以下步骤1)、将待检测的心脏建立B样条曲面模型P(u，v)，其中u，v是B样条曲面定义时的规范化参数，参数域均为，设模型上存在n个位移已知的标记点ai(i＝1，2，...n)，ai的在球面坐标系下的空间位置为[θiψiRi]，在P(u，v)上的规范化坐标为[ui，vi]，且ai的位移为[rθirψirRi]，通过离散的n个点的位移矢量建立整个模型的rθ、rψ和rR的连续位移场；2)、利用规范化坐标u，v拟合位移场并将位移场中的点与实际的空间点建立对应关系，对三个位移分别拟合，其中rθ位移的拟合过程使用[uivirθi]拟合出了场曲面Q(s，t)，根据B样条曲面特点Q(s，t)实际包括3个方程，即Qu(s，t)＝u，Qv(s，t)＝v，Qθ(s，t)＝rθ；已知[θψR]＝P(u，v)，[uvrθ]＝Q(s，t)，定义P1和Q1分别为P(u，v)与Q(s，t)的逆运算；对模型上任意点a0[θ0ψ0R0]，计算其坐标值θ的位移量rθ0的值，过程如下Step1将坐标转换至规范坐标系上，[u0v0]＝P1(θ0ψ0R0)；Step2计算该点在位移场中坐标s0，t0，[s0t0]＝Q1(u0v0)；Step3得到位移rθ，rθ＝Qθ(s0，t0)；rψ位移的拟合过程使用[uivirψi]拟合出了场曲面Q(s，t)，且有Qu(s，t)＝u，Qv(s，t)＝v，Qψ(s，t)＝rψ；[θψR]＝P(u，v)，[uvrψ]＝Q(s，t)；设计算点a0[θ0ψ0R0]的rψ位移值，过程如下Step1将坐标转换至规范坐标系上，[u0v0]＝P1(ψ0ψ0R0)；Step2计算该点在位移场中坐标s0，t0，[s0t0]＝Q1(u0v0)；Step3得到位移rψ，rψ＝Qψ(s0t0)；rR位移的拟合过程使用[uivirRi]拟合出了场曲面Q(s，t)，且有Qu(s，t)＝u，Qv(s，t)＝v，QR(s，t)＝rR；[θψR]]＝P(u，v)，[uvrR]＝Q(s，t)；设计算点a0[θ0ψ0R0]的rR位移值，过程如下Step1将坐标转换至规范坐标系上，[u0v0]＝P1(R0ψ0R0)；Step2计算该点在位移场中坐标s0，t0，[s0t0]＝Q1(u0v0)；Step3得到位移rR，rR＝QR(s0，t0)；3)、将步骤2)中心壁上任意点a0[θ0ψ0R0]的位移量分别记做Nθ(θ0ψ0R0)、Nψ(θ0ψ0R0)和NR(θ0ψ0R0)；由于应变是在笛卡尔坐标系x，y，z下定义的，我们需将该定义转化至球面坐标系下，设a0在笛卡尔系下坐标为[x0y0z0]，其位移量分别为Mx(x0，y0，z0)，My(x0，y0，z0)，Mz(x0，y0，z0)，则根据公式(8)直接由球面坐标系下位移求出应变，并根据(6)求出应力<mrow><mfencedopen='{'close=''><mtable><mtr><mtd><msub><mi>ϵ</mi><mi>x</mi></msub><mo>=</mo><mfrac><mrow><mo>∂</mo><mi>u</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><msup><mi>M</mi><mi>x</mi></msup><mrow><mo>(</mo><msub><mi>x</mi><mn>0</mn></msub><mo>+</mo><mi>dx</mi><mo>,</mo><msub><mi>y</mi><mn>0</mn></msub><mo>,</mo><msub><mi>z</mi><mn>0</mn></msub><mo>)</mo></mrow><mo>-</mo><msup><mi>M</mi><mi>x</mi></msup><mrow><mo>(</mo><msub><mi>x</mi><mn>0</mn></msub><mo>,</mo><msub><mi>y</mi><mn>0</mn></msub><mo>,</mo><msub><mi>z</mi><mn>0</mn></msub><mo>)</mo></mrow></mrow><mi>dx</mi></mfrac><mo>=</mo><mfrac><mrow><msup><mi>N</mi><mi>θ</mi></msup><mrow><mo>(</mo><msub><mi>θ</mi><mn>0</mn></msub><mo>+</mo><mi>dθ</mi><mo>,</mo><msub><mi>ψ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>R</mi><mn>0</mn></msub><mo>)</mo></mrow><mo>-</mo><msup><mi>N</mi><mi>θ</mi></msup><mrow><mo>(</mo><msub><mi>θ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>ψ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>R</mi><mn>0</mn></msub><mo>)</mo></mrow></mrow><mrow><msub><mi>R</mi><mn>0</mn></msub><mi>dθ</mi></mrow></mfrac></mtd></mtr><mtr><mtd><msub><mi>ϵ</mi><mi>y</mi></msub><mo>=</mo><mfrac><mrow><mo>∂</mo><mi>v</mi></mrow><mrow><mo>∂</mo><mi>y</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><msup><mi>M</mi><mi>y</mi></msup><mrow><mo>(</mo><msub><mi>x</mi><mn>0</mn></msub><mo>,</mo><msub><mi>y</mi><mn>0</mn></msub><mo>+</mo><mi>dy</mi><mo>,</mo><msub><mi>z</mi><mn>0</mn></msub><mo>)</mo></mrow><mo>-</mo><msup><mi>M</mi><mi>y</mi></msup><mrow><mo>(</mo><msub><mi>x</mi><mn>0</mn></msub><mo>,</mo><msub><mi>y</mi><mn>0</mn></msub><mo>,</mo><msub><mi>z</mi><mn>0</mn></msub><mo>)</mo></mrow></mrow><mi>dy</mi></mfrac><mo>=</mo><mfrac><mrow><msup><mi>N</mi><mi>ψ</mi></msup><mrow><mo>(</mo><msub><mi>θ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>ψ</mi><mn>0</mn></msub><mo>+</mo><mi>dψ</mi><mo>,</mo><msub><mi>R</mi><mn>0</mn></msub><mo>)</mo></mrow><mo>-</mo><msup><mi>N</mi><mi>ψ</mi></msup><mrow><mo>(</mo><msub><mi>θ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>ψ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>R</mi><mn>0</mn></msub><mo>)</mo></mrow></mrow><mrow><msub><mi>R</mi><mn>0</mn></msub><mi>dψ</mi></mrow></mfrac></mtd></mtr><mtr><mtd><msub><mi>ϵ</mi><mi>z</mi></msub><mo>=</mo><mfrac><mrow><mo>∂</mo><mi>w</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><msup><mi>M</mi><mi>z</mi></msup><mrow><mo>(</mo><msub><mi>x</mi><mn>0</mn></msub><mo>,</mo><msub><mi>y</mi><mn>0</mn></msub><mo>,</mo><msub><mi>z</mi><mn>0</mn></msub><mo>+</mo><mi>dz</mi><mo>)</mo></mrow><mo>-</mo><msup><mi>M</mi><mi>z</mi></msup><mrow><mo>(</mo><msub><mi>x</mi><mn>0</mn></msub><mo>,</mo><msub><mi>y</mi><mn>0</mn></msub><mo>,</mo><msub><mi>z</mi><mn>0</mn></msub><mo>)</mo></mrow></mrow><mi>dz</mi></mfrac><mo>=</mo><mfrac><mrow><msup><mi>N</mi><mi>R</mi></msup><mrow><mo>(</mo><msub><mi>θ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>ψ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>R</mi><mn>0</mn></msub><mo>+</mo><mi>dR</mi><mo>)</mo></mrow><mo>-</mo><msup><mi>N</mi><mi>R</mi></msup><mrow><mo>(</mo><msub><mi>θ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>ψ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>R</mi><mn>0</mn></msub><mo>)</mo></mrow></mrow><mi>dR</mi></mfrac></mtd></mtr><mtr><mtd><msub><mi>γ</mi><mi>xy</mi></msub><mo>=</mo><mfrac><mrow><mo>∂</mo><mi>u</mi></mrow><mrow><mo>∂</mo><mi>y</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mi>v</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><msup><mi>N</mi><mi>θ</mi></msup><mrow><mo>(</mo><msub><mi>θ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>ψ</mi><mn>0</mn></msub><mo>+</mo><mi>dψ</mi><mo>,</mo><msub><mi>R</mi><mn>0</mn></msub><mo>)</mo></mrow><mo>-</mo><msup><mi>N</mi><mi>θ</mi></msup><mrow><mo>(</mo><msub><mi>θ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>ψ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>R</mi><mn>0</mn></msub><mo>)</mo></mrow></mrow><mrow><msub><mi>R</mi><mn>0</mn></msub><mi>dψ</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><msup><mi>N</mi><mi>ψ</mi></msup><mrow><mo>(</mo><msub><mi>θ</mi><mn>0</mn></msub><mo>+</mo><mi>dθ</mi><mo>,</mo><msub><mi>ψ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>R</mi><mn>0</mn></msub><mo>)</mo></mrow><mo>-</mo><msup><mi>N</mi><mi>ψ</mi></msup><mrow><mo>(</mo><msub><mi>θ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>ψ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>R</mi><mn>0</mn></msub><mo>)</mo></mrow></mrow><mrow><msub><mi>R</mi><mn>0</mn></msub><mi>dθ</mi></mrow></mfrac></mtd></mtr><mtr><mtd><msub><mi>γ</mi><mi>yz</mi></msub><mo>=</mo><mfrac><mrow><mo>∂</mo><mi>v</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mi>w</mi></mrow><mrow><mo>∂</mo><mi>y</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><msup><mi>N</mi><mi>ψ</mi></msup><mrow><mo>(</mo><msub><mi>θ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>ψ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>R</mi><mn>0</mn></msub><mo>+</mo><mi>dR</mi><mo>)</mo></mrow><mo>-</mo><msup><mi>N</mi><mi>ψ</mi></msup><mrow><mo>(</mo><msub><mi>θ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>ψ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>R</mi><mn>0</mn></msub><mo>)</mo></mrow></mrow><mi>dR</mi></mfrac><mo>+</mo><mfrac><mrow><msup><mi>N</mi><mi>R</mi></msup><mrow><mo>(</mo><msub><mi>θ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>ψ</mi><mn>0</mn></msub><mo>+</mo><mi>dψ</mi><mo>,</mo><msub><mi>R</mi><mn>0</mn></msub><mo>)</mo></mrow><mo>-</mo><msup><mi>N</mi><mi>R</mi></msup><mrow><mo>(</mo><msub><mi>θ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>ψ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>R</mi><mn>0</mn></msub><mo>)</mo></mrow></mrow><mrow><msub><mi>R</mi><mn>0</mn></msub><mi>dψ</mi></mrow></mfrac></mtd></mtr><mtr><mtd><msub><mi>γ</mi><mi>xz</mi></msub><mo>=</mo><mfrac><mrow><mo>∂</mo><mi>u</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mi>w</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><msup><mi>N</mi><mi>θ</mi></msup><mrow><mo>(</mo><msub><mi>θ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>ψ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>R</mi><mn>0</mn></msub><mo>+</mo><mi>dR</mi><mo>)</mo></mrow><mo>-</mo><msup><mi>N</mi><mi>θ</mi></msup><mrow><mo>(</mo><msub><mi>θ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>ψ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>R</mi><mn>0</mn></msub><mo>)</mo></mrow></mrow><mi>dR</mi></mfrac><mo>+</mo><mfrac><mrow><msup><mi>N</mi><mi>R</mi></msup><mrow><mo>(</mo><msub><mi>θ</mi><mn>0</mn></msub><mo>+</mo><mi>dθ</mi><mo>,</mo><msub><mi>ψ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>R</mi><mn>0</mn></msub><mo>)</mo></mrow><mo>-</mo><msup><mi>N</mi><mi>R</mi></msup><mrow><mo>(</mo><msub><mi>θ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>ψ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>R</mi><mn>0</mn></msub><mo>)</mo></mrow></mrow><mrow><msub><mi>R</mi><mn>0</mn></msub><mi>dθ</mi></mrow></mfrac></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>8</mn><mo>)</mo></mrow></mrow>笛卡尔坐标系下的应变定义物体位移向量δ由x，y，z轴上的位移u，v，w组成。应变向量ε有3个主应变和三个剪应变，与位移的关系为<mrow><mfencedopen='{'close=''><mtable><mtr><mtd><msub><mi>ϵ</mi><mi>x</mi></msub><mo>=</mo><mfrac><mrow><mo>∂</mo><mi>u</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac></mtd></mtr><mtr><mtd><msub><mi>ϵ</mi><mi>y</mi></msub><mo>=</mo><mfrac><mrow><mo>∂</mo><mi>v</mi></mrow><mrow><mo>∂</mo><mi>y</mi></mrow></mfrac></mtd></mtr><mtr><mtd><msub><mi>ϵ</mi><mi>z</mi></msub><mo>=</mo><mfrac><mrow><mo>∂</mo><mi>w</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac></mtd></mtr><mtr><mtd><msub><mi>γ</mi><mi>xy</mi></msub><mo>=</mo><mfrac><mrow><mo>∂</mo><mi>u</mi></mrow><mrow><mo>∂</mo><mi>y</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mi>v</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac></mtd></mtr><mtr><mtd><msub><mi>γ</mi><mi>yz</mi></msub><mo>=</mo><mfrac><mrow><mo>∂</mo><mi>v</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mi>w</mi></mrow><mrow><mo>∂</mo><mi>y</mi></mrow></mfrac></mtd></mtr><mtr><mtd><msub><mi>γ</mi><mi>xz</mi></msub><mo>=</mo><mfrac><mrow><mo>∂</mo><mi>u</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mi>w</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>5</mn><mo>)</mo></mrow></mrow>应力向量σ定义如下{σ}＝[σxσyσzτxyτyzτxz]＝Dε(6)其中，D为弹性矩阵，由研究对象的弹性模量E及泊松比μ确定，对于确定的材料可视为已知常量。2.如权利要求1所述的基于四维医学图像的心壁应力应变测量方法，其特征在于在所述步骤1)中，B样条曲面模型建立过程由一组控制点di(i=0，1，2，...，η)生成的k次B样条曲线定义为其中，Niik(U)为B样条基函数，其经典的de-BoorCox定义式为一组递推公式，即N⑷-J1'若u,+1式中Ui为一组非递减规范化序列,称作节点序列,通常取i=0,1,2，·.·,η+k+l；同理由(m+1)X(n+1)个控制点生成的三维B样条曲面定义为其中，kl、k2分别为u和ν方向上的曲线阶数，样条曲面就是样条的推广，它将一组空间坐标下的点转换到了一个1X1的空间u-v内，将u-v空间定义为规范化空间。全文摘要一种基于四维医学图像的心壁应力应变测量方法，包括以下步骤1)将待检测的心脏建立B样条曲面模型P(u，v)，模型上存在n个位移已知标记点ai(i＝1，2，...n)，ai对应的规范化坐标为[ui，vi]，ai的位移为[rθirψirRi]，通过离散的n个位移建立一个连续的位移场；2)将位移场中的点与运动点的位置建立联系，对三个位移分别拟合；3)将步骤2)中心壁上任意点的位移向量，并根据得到全局坐标下d0的位置及对应位移，接着应用旋转矩阵T，即公式(9)，将位置及位移旋转至局部坐标下；再用公式(8)和(6)分别计算应变和应力。本发明形变匹配性良好、精度高、大大减少计算复杂度。文档编号A61B5/11GK101926648SQ20091010107公开日2010年12月29日申请日期2009年8月3日优先权日2009年8月3日发明者刘盛,杜雅慧,王万良,王芳,管秋,蒋超,陈胜勇申请人:浙江工业大学