CN101926648B

CN101926648B - Cardiac wall stress-strain measuring method based on four-dimensional medical image

Info

Publication number: CN101926648B
Application number: CN2009101010724A
Authority: CN
Inventors: 管秋; 杜雅慧; 陈胜勇; 刘盛; 蒋超; 王芳; 王万良
Original assignee: Zhejiang University of Technology ZJUT
Current assignee: Zhejiang University of Technology ZJUT
Priority date: 2009-08-03
Filing date: 2009-08-03
Publication date: 2012-01-25
Anticipated expiration: 2029-08-03
Also published as: CN101926648A

Abstract

A heart wall stress-strain measurement method based on a four-dimensional medical image, comprising the following steps: 1) establishing a B-spline surface model P(u, v) of the heart to be detected, and there are n marker points ai( i=1, 2,...n), the normalized coordinates corresponding to a _i are [u _i , v _i ], the displacement of a _i is [rθ _i rψ _i rR _i ], and a continuous 2) Connect the points in the displacement field with the position of the moving point, and fit the three displacements respectively; 3) Take the displacement vector of any point on the central wall in step 2) and obtain the global coordinates according to d ₀ position and the corresponding displacement, and then apply the rotation matrix T, that is, the formula (9), to rotate the position and displacement to the local coordinates; then use the formulas (8) and (6) to calculate the strain and stress respectively. The invention has good deformation matching, high precision and greatly reduces computational complexity.

Description

Measuring method of heart wall stress and strain based on 4D medical image

技术领域 technical field

本发明涉及图像处理、计算机视觉、医学图像分析、心血管医学、机械动力学和运动学领域的一种心壁应力应变测量方法。 The invention relates to a method for measuring heart wall stress and strain in the fields of image processing, computer vision, medical image analysis, cardiovascular medicine, mechanical dynamics and kinematics. the

背景技术 Background technique

心脏是人体血液的动力中枢，因此心脏功能的分析研究也备受关注，而左心室在整个心动过程中提供主要的动能，因此，对心脏功能的分析研究可以转化为对左心室的分析研究。左心室的功能参数有很多，比如全局参数：左心室容积(LVV)，射血份数(EF)，心室壁厚(WT)等，这些参数可以在一定程度上反映整个心室的功能强弱，但却无法定位具体的病变部位，因而近年来学界加大了对局部参数的研究力度，其中最主要的局部参数是应力与应变。 The heart is the power center of human blood, so the analysis and research of heart function has also attracted much attention, and the left ventricle provides the main kinetic energy during the whole cardiac process. Therefore, the analysis and research of heart function can be transformed into the analysis and research of left ventricle. There are many functional parameters of the left ventricle, such as global parameters: left ventricular volume (LVV), ejection fraction (EF), ventricular wall thickness (WT), etc. These parameters can reflect the function of the entire ventricle to a certain extent, However, it is impossible to locate the specific lesion site. Therefore, in recent years, the academic circle has intensified the research on local parameters, and the most important local parameters are stress and strain. the

应变与应力是指柔性体在受到一定外力作用的情况下发生的形变以及由这种形变引起的抵御外力的自身作用力，应力与应变可以很好地反应柔性体在不同部位所受到的力以其形变情况。就左心室而言，如果某一部位发生了病变从而导致其心室壁的性质发生变化，必然会体现在应力应变图上。在医学研究方面，也有很多文献开始以应力应变作为指标研究不同的心脏疾病，有的用左心室舒张期应力的变化率来评价心室的整体功能，也有通过对比心力衰竭患者与正常人左心室应变值研究分析了应变与心脏功能分析的关系，或者研究了冠心病与高血压患者在应力应变参数上的变化。 Strain and stress refer to the deformation of the flexible body when it is subjected to a certain external force and the self-acting force caused by this deformation to resist the external force. Stress and strain can well reflect the force received by the flexible body in different parts. its deformation. As far as the left ventricle is concerned, if a lesion occurs in a certain part, which leads to changes in the properties of the ventricular wall, it will inevitably be reflected in the stress-strain diagram. In terms of medical research, there are also many literatures that use stress and strain as indicators to study different heart diseases. Some use the rate of change of left ventricular diastolic stress to evaluate the overall function of the ventricle, and some compare the left ventricular strain of patients with heart failure and normal people. Value studies analyze the relationship between strain and cardiac function analysis, or study changes in stress-strain parameters in patients with coronary heart disease and hypertension. the

目前对应力与应变的计算主要有两种方式，即从应力到应变和从应变到应力。前一种方式通常以心室运动的电位变化为依据建立电位模型，由模型得到心室的应力，再利用有限元方法计算对应位置的应变。后一种方式通常是由心室壁的位移得到应变，再由应变计算出应力。本文采用的是从应变计算应力的方式，因为该方式可以以心室图片为数据源，无需接触，且对于不同个体的适应性较好。由于应力应变之间的关系遵循一定力学对应关系(弹性力学已经给出)，因而该方式的关键在于如何获得应变值，而应变是位移相对于位置的导数，所以求取应变的关键在于如何求解左心室壁各处的位移。 At present, there are mainly two ways to calculate stress and strain, that is, from stress to strain and from strain to stress. The former method usually establishes a potential model based on the potential change of ventricular movement, obtains the stress of the ventricle from the model, and then uses the finite element method to calculate the strain at the corresponding position. The latter method usually obtains the strain from the displacement of the ventricular wall, and then calculates the stress from the strain. This paper adopts the method of calculating stress from strain, because this method can use ventricle pictures as the data source, without contact, and has better adaptability to different individuals. Since the relationship between stress and strain follows a certain mechanical correspondence (elastic mechanics has been given), the key to this method lies in how to obtain the strain value, and strain is the derivative of displacement with respect to position, so the key to obtaining strain is how to solve Displacement throughout the left ventricular wall. the

通过各种成像和标记方法，我们可以很方便地得到心动周期不同时刻的左心室模型，要得到位移就需要在不同时刻的模型上找到同一点的位置。目前的主流方法是将不同时刻模型分别划分为有限个单元，在模型之间寻找最相似的单元从而进行位置匹配最终计算出位移和应变。这种通过匹配单元计算位移的方法有其无法克服的问题，首先它的前提是前后时刻的形变是微小的，形状是大致相同的，事实上左心室壁的位移并不严格满足以上假设；其次是它虽然能计算单元结点处的位移，但对于单元内部的点依然要采用近似插值计算，不够准确；第三是为了保证匹配精度每次计算必须对整个左心室室壁的所有单元都进行计算，从而剧增了运算量。 Through various imaging and marking methods, we can easily obtain the left ventricle model at different moments of the cardiac cycle. To obtain the displacement, we need to find the same point on the model at different moments. The current mainstream method is to divide the models at different times into finite units, and find the most similar units between the models to perform position matching and finally calculate the displacement and strain. This method of calculating displacement by matching units has its insurmountable problems. First, its premise is that the deformation at the front and rear moments is small and the shape is roughly the same. In fact, the displacement of the left ventricular wall does not strictly meet the above assumptions; secondly Although it can calculate the displacement at the node of the unit, it still needs to use approximate interpolation calculation for the points inside the unit, which is not accurate enough; the third is that in order to ensure the matching accuracy, each calculation must be performed on all units of the entire left ventricle wall Calculation, thus dramatically increasing the amount of computation. the

发明内容 Contents of the invention

为了克服已有的心壁应力应变测量方法的形变匹配性差、精度低、计算复杂的不足，本发明提供一种形变匹配性良好、精度高、大大减少计算复杂度的基于四维医学图像的心壁应力应变测量方法。 In order to overcome the deficiencies of poor deformation matching, low precision, and complicated calculation of the existing heart wall stress and strain measurement methods, the present invention provides a heart wall based on four-dimensional medical images with good deformation matching, high precision, and greatly reduced computational complexity. Stress-strain measurement method. the

本发明解决其技术问题所采用的技术方案是： The technical solution adopted by the present invention to solve its technical problems is:

一种基于四维医学图像的心壁应力应变测量方法，所述心壁应力应变测量方法包括以下步骤： A method for measuring heart wall stress-strain based on four-dimensional medical images, said method for measuring heart wall stress-strain comprises the following steps:

1)、将待检测的心脏建立B样条曲面模型P(u，v)，其中u，v是B样条曲面定义时的规范化参数，参数域均为[0，1]；设模型上存在n个位移已知的标记点a_i(i＝1，2，...n)，a_i的在球面坐标系下的空间位置为[θ_i ψ_i R_i]，在P(u，v)上的规范化坐标为[u_i，v_i]，且a_i的位移为[rθ_i rψ_i rR_i]，通过离散的n个点的位移矢量建立整个模型的rθ、rψ和rR的连续位移场； 1), establish the B-spline surface model P(u, v) of the heart to be detected, where u, v are the normalization parameters when the B-spline surface is defined, and the parameter fields are all [0, 1]; n marker points a _i (i=1, 2,...n) whose displacement is known, the spatial position of a _i in the spherical coordinate system is [θ _i ψ _i R _i ], in P(u, v The normalized coordinates on ) are [u _i , v _i ], and the displacement of a _i is [rθ _i rψ _i rR _i ], the continuous displacement of rθ, rψ and rR of the whole model is established by discrete displacement vectors of n points field;

2)、利用规范化坐标u，v拟合位移场并将位移场中的点与实际的空间点建立对应关系，对三个位移分别拟合，其中： 2), use the standardized coordinates u, v to fit the displacement field and establish a corresponding relationship between the points in the displacement field and the actual space points, and respectively fit the three displacements, wherein:

rθ位移的拟合过程：使用[u_i v_i rθ_i]拟合出了场曲面Q(s，t)，根据B样条曲面特点Q(s，t)实际包括3个方程，即Q^u(s，t)＝u，Q^v(s，t)＝v，Q^θ(s，t)＝rθ；已知[θ ψR]＝P(u，v)，[u v rθ]＝Q(s，t)，定义P^-1和Q^-1分别为P(u，v)与Q(s，t)的逆运算；对模型上任意点a₀[θ₀ ψ₀ R₀]，计算其坐标值θ的位移量rθ₀的值，过程如下： The fitting process of rθ displacement: using [u _i v _i rθ _i ] to fit the field surface Q(s, t), according to the characteristics of the B-spline surface Q(s, t) actually includes three equations, that is, Q ^u (s, t)=u, Q ^v (s, t)=v, Q ^θ (s, t)=rθ; known [θ ψR]=P(u,v), [u v rθ]=Q(s , t), define P ^-1 and Q ^-1 as the inverse operation of P(u, v) and Q(s, t) respectively; for any point a ₀ [θ ₀ ψ ₀ R ₀ ] on the model, calculate its coordinates The value of the displacement rθ ₀ of the value θ, the process is as follows:

Step1：将坐标转换至规范坐标系上，[u₀ v₀]＝P^-1(θ₀ ψ₀ R₀)； Step1: Convert the coordinates to the standard coordinate system, [u ₀ v ₀ ]=P ^-1 (θ ₀ ψ ₀ R ₀ );

Step2：计算该点在位移场中坐标s₀，t₀，[s₀ t₀]＝Q^-1(u₀ v₀)； Step2: Calculate the coordinates s ₀ , t ₀ of the point in the displacement field, [s ₀ t ₀ ]=Q ^-1 (u ₀ v ₀ );

Step3：得到位移rθ，rθ＝Q^θ(s₀，t₀)； Step3: Obtain the displacement rθ, rθ=Q ^θ (s ₀ , t ₀ );

rψ位移的拟合过程：使用[u_i v_i rψ_i]拟合出了场曲面Q(s，t)，且有Q^u(s，t)＝u，Q^v(s，t)＝v，Q^ψ(s，t)＝rψ；[θψ R]＝P(u，v)，[u v rψ]＝Q(s，t)；设计算点a₀[θ₀ ψ₀ R₀]的rψ位移值，过程如下： The fitting process of rψ displacement: use [u _i v _i rψ _i ] to fit the field surface Q(s, t), and Q ^u (s, t)=u, Q ^v (s, t)=v , Q ^ψ (s, t) = rψ; [θψ R] = P(u, v), [u v rψ] = Q(s, t); design calculation point a ₀ [θ ₀ ψ ₀ R ₀ ] rψ Displacement value, the process is as follows:

Step1：将坐标转换至规范坐标系上，[u₀ v₀]＝P^-1(ψ₀ ψ₀ R₀)； Step1: Convert the coordinates to the standard coordinate system, [u ₀ v ₀ ]=P ^-1 (ψ ₀ ψ ₀ R ₀ );

Step3：得到位移rψ，rψ＝Q^ψ(s₀，t₀)； Step3: Obtain the displacement rψ, rψ=Q ^ψ (s ₀ , t ₀ );

rR位移的拟合过程：使用[u_i v_i rR_i]拟合出了场曲面Q(s，t)，且有Q^u(s，t)＝u，Q^v(s，t)＝v，Q^R(s，t)＝rR；[θ ψ R]]＝P(u，v)，[u v rR]＝Q(s，t)；设计算点a₀[θ₀ ψ₀ R₀]的rR位移值，过程如下： Fitting process of rR displacement: use [u _i v _i rR _i ] to fit the field surface Q(s, t), and Q ^u (s, t)=u, Q ^v (s, t)=v , Q ^R (s, t) = rR; [θ ψ R]] = P (u, v), [u v rR] = Q (s, t); design calculation point a ₀ [θ ₀ ψ ₀ R ₀ ] rR displacement value, the process is as follows:

Step1：将坐标转换至规范坐标系上，[u₀ v₀]＝P^-1(R₀ ψ₀ R₀)； Step1: Convert the coordinates to the standard coordinate system, [u ₀ v ₀ ]=P ^-1 (R ₀ ψ ₀ R ₀ );

Step3：得到位移rR，rR＝Q^R(s₀，t₀)； Step3: Obtain displacement rR, rR=Q ^R (s ₀ , t ₀ );

3)、将步骤2)中心壁上任意点a₀[θ₀ ψ₀ R₀]的位移量分别记做N^θ(θ₀ ψ₀ R₀)、N^ψ(θ₀ψ₀ R₀)和N^R(θ₀ ψ₀ R₀)；由于应变是在笛卡尔坐标系x，y，z下定义的，我们需将该定义转化至球面坐标系下，设a₀在笛卡尔系下坐标为[x₀ y₀ z₀]，其位移量分别为M^x(x₀，y₀，z₀)，M^y(x₀，y₀，z₀)，M^z(x₀，y₀，z₀)，则根据公式(8)可直接由球面坐标系下位移求出应变，并根据(6)求出应力。 3) Record the displacement of any point a ₀ [θ ₀ ψ ₀ R ₀ ] on the central wall in step 2) as N ^θ (θ ₀ ψ ₀ R ₀ ), N ^ψ (θ ₀ ψ ₀ R ₀ ) and N ^R (θ ₀ ψ ₀ R ₀ ); since the strain is defined in the Cartesian coordinate system x, y, z, we need to convert this definition to the spherical coordinate system, let the coordinates of a ₀ in the Cartesian system be [x ₀ y ₀ z ₀ ], the displacements are M ^x (x ₀ , y ₀ , z ₀ ), M ^y (x ₀ , y ₀ , z ₀ ), M ^z (x ₀ , y ₀ , z ₀ ), the strain can be obtained directly from the displacement in the spherical coordinate system according to formula (8), and the stress can be obtained according to (6).

$\{\begin{matrix} {ϵ ϵ}_{x x} = = \frac{&PartialD; &PartialD; u u}{&PartialD; &PartialD; x x} = = \frac{{M m}^{x x} (({x x}_{00} + + dx dx,, {y the y}_{00},, {z z}_{00})) - - {M m}^{x x} (({x x}_{00},, {y the y}_{00},, {z z}_{00}))}{dx dx} = = \frac{{N N}^{θ θ} (({θ θ}_{00} + + dθ dθ,, {ψ ψ}_{00},, {R R}_{00})) - - {N N}^{θ θ} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00}))}{{R R}_{00} dθ dθ} \\ {ϵ ϵ}_{y the y} = = \frac{&PartialD; &PartialD; v v}{&PartialD; &PartialD; y the y} = = \frac{{M m}^{y the y} (({x x}_{00},, {y the y}_{00} + + dy dy,, {z z}_{00})) - - {M m}^{y the y} (({x x}_{00},, {y the y}_{00},, {z z}_{00}))}{dy dy} = = \frac{{N N}^{ψ ψ} (({θ θ}_{00},, {ψ ψ}_{00} + + dψ dψ,, {R R}_{00})) - - {N N}^{ψ ψ} (({θ θ}_{00} {,, ψ ψ}_{00},, {R R}_{00}))}{{R R}_{00} dψ dψ} \\ {ϵ ϵ}_{z z} = = \frac{&PartialD; &PartialD; w w}{&PartialD; &PartialD; z z} = = \frac{{M m}^{z z} (({x x}_{00},, {y the y}_{00},, {z z}_{00} + + dz dz)) - - {M m}^{z z} (({x x}_{00},, {y the y}_{00},, {z z}_{00}))}{dz dz} = = \frac{{N N}^{R R} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00} + + dR d)) - - {N N}^{R R} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00}))}{dR d} \\ {γ γ}_{xy xy} = = \frac{&PartialD; &PartialD; u u}{&PartialD; &PartialD; y the y} + + \frac{&PartialD; &PartialD; v v}{&PartialD; &PartialD; x x} = = \frac{{N N}^{θ θ} (({θ θ}_{00},, {ψ ψ}_{00} + + dψ dψ,, {R R}_{00})) - - {N N}^{θ θ} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00}))}{{R R}_{00} dψ dψ} + + \frac{{N N}^{ψ ψ} (({θ θ}_{00} + + dθ dθ,, {ψ ψ}_{00},, {R R}_{00})) - - {N N}^{ψ ψ} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00}))}{{R R}_{00} dθ dθ} \\ {γ γ}_{yz yz} = = \frac{&PartialD; &PartialD; v v}{&PartialD; &PartialD; z z} + + \frac{&PartialD; &PartialD; w w}{&PartialD; &PartialD; y the y} = = \frac{{N N}^{ψ ψ} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00} + + dR d)) - - {N N}^{ψ ψ} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00}))}{dR d} + + \frac{{N N}^{R R} (({θ θ}_{00},, {ψ ψ}_{00} + + dψ dψ,, {R R}_{00})) - - {N N}^{R R} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00}))}{{R R}_{00} dψ dψ} \\ {γ γ}_{xz xz} = = \frac{&PartialD; &PartialD; u u}{&PartialD; &PartialD; z z} + + \frac{&PartialD; &PartialD; w w}{&PartialD; &PartialD; x x} = = \frac{{N N}^{θ θ} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00} + + dR d)) - - {N N}^{θ θ} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00}))}{dR d} + + \frac{{N N}^{R R} (({θ θ}_{00} + + dθ dθ,, {ψ ψ}_{00},, {R R}_{00})) - - {N N}^{R R} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00}))}{{R R}_{00} dθ dθ} \end{matrix} - - - - - - ((88))$

笛卡尔坐标系下的应变定义：物体位移向量δ由x，y，z轴上的位移u，v，w组成。应变向量ε有3个主应变和三个剪应变，与位移的关系为： Definition of strain in the Cartesian coordinate system: The object displacement vector δ consists of displacements u, v, and w on the x, y, and z axes. The strain vector ε has three principal strains and three shear strains, and the relationship with displacement is:

$\{\begin{matrix} {ϵ ϵ}_{x x} = = \frac{&PartialD; &PartialD; u u}{&PartialD; &PartialD; x x} \\ {ϵ ϵ}_{y the y} = = \frac{&PartialD; &PartialD; v v}{&PartialD; &PartialD; y the y} \\ {ϵ ϵ}_{z z} = = \frac{&PartialD; &PartialD; w w}{&PartialD; &PartialD; z z} \\ {γ γ}_{xy xy} = = \frac{&PartialD; &PartialD; u u}{&PartialD; &PartialD; y the y} + + \frac{&PartialD; &PartialD; v v}{&PartialD; &PartialD; x x} \\ {γ γ}_{yz yz} = = \frac{&PartialD; &PartialD; v v}{&PartialD; &PartialD; z z} + + \frac{&PartialD; &PartialD; w w}{&PartialD; &PartialD; y the y} \\ {γ γ}_{xz xz} = = \frac{&PartialD; &PartialD; u u}{&PartialD; &PartialD; z z} + + \frac{&PartialD; &PartialD; w w}{&PartialD; &PartialD; x x} \end{matrix} - - - - - - ((55))$

应力向量σ定义如下： The stress vector σ is defined as follows:

{σ}＝[σ_x σ_y σ_z τ_xy τ_yz τ_xz]＝Dε (6) {σ}＝[σ _x σ _y σ _z τ _xy τ _yz τ _xz ]＝Dε (6)

其中，D为弹性矩阵，由研究对象的弹性模量E及泊松比μ确定，对于确定的材料可视为已知常量。 Among them, D is the elastic matrix, which is determined by the elastic modulus E and Poisson's ratio μ of the research object, and can be regarded as a known constant for a certain material. the

作为优选的一种方案：在所述步骤1)中，B样条曲面模型建立过程：由一组控制点di(i＝0，1，2，...，n)生成的k次B样条曲线定义为： As a preferred solution: in the step 1), the B-spline surface model building process: k times of B-sample generated by a group of control points di (i=0, 1, 2, ..., n) A curve is defined as:

$P P ((u u)) = = {Σ Σ}_{i i = = 00}^{n no} {d d}_{i i} {N N}_{i i,, k k} ((u u)),, {u u}_{i i} \leq \leq u u < < {u u}_{i i + + 11} - - - - - - ((11))$

其中，N_i，k(u)为B样条基函数，其经典的de-Boor Cox定义式为一组递推公式，即： Among them, N _{i, k} (u) is the B-spline basis function, and its classic de-Boor Cox definition formula is a set of recursive formulas, namely:

式中u_i为一组非递减规范化序列，称作节点序列，通常取i＝0，1，2，...，n+k+1；同理由(m+1)×(n+1)个控制点生成的三维B样条曲面定义为： In the formula, u _i is a group of non-decreasing normalized sequences, called node sequences, usually take i=0, 1, 2,..., n+k+1; the same reason (m+1)×(n+1) The three-dimensional B-spline surface generated by control points is defined as:

$P P ((u u,, v v)) = = {Σ Σ}_{i i = = 00}^{m m} {Σ Σ}_{j j = = 00}^{n no} {d d}_{i i,, j j} {N N}_{i i,, k k 11} ((u u,, v v)) {N N}_{j j,, k k 22} ((u u,, v v)),, {u u}_{i i} \leq \leq u u < < {u u}_{i i + + 11},, {v v}_{j j} \leq \leq v v < < {v v}_{j j + + 11} - - - - - - ((33))$

其中，k1、k2分别为u和v方向上的曲线阶数，样条曲面就是样条的推广，它将一组空间坐标下的点转换到了一个1×1的空间u-v内，将u-v空间定义为规范化空间。 Among them, k1 and k2 are the curve orders in the u and v directions respectively, and the spline surface is the extension of the spline, which transforms the points under a set of space coordinates into a 1×1 space u-v, and defines the u-v space for the normalized space. the

本发明的有益效果主要表现在：形变匹配性良好、精度高、大大减少计算复杂度。 The beneficial effects of the invention are mainly manifested in: good deformation matching, high precision and greatly reduced computational complexity. the

附图说明 Description of drawings

图1是NURBS曲面逆运算流程图。 Figure 1 is a flowchart of the inverse operation of NURBS surfaces. the

图2是应力应变定义及壳状单元的示意图。 Figure 2 is a schematic diagram of the definition of stress and strain and shell elements. the

图3是位移场与空间位置对应关系的示意图。 Fig. 3 is a schematic diagram of the corresponding relationship between the displacement field and the spatial position. the

图4是左心室7个NURBS模型的示意图。 Figure 4 is a schematic diagram of 7 NURBS models of the left ventricle. the

图5是左心室内壁第一相位标记点的位移矢量图。 Fig. 5 is a vector diagram of the displacement of the first phase marker point on the inner wall of the left ventricle. the

图6是左心室内壁第一相位标记点的对比图。 Fig. 6 is a comparison diagram of the first phase marker point on the inner wall of the left ventricle. the

图7是左心室内壁第一相位x方向位移场曲面图。 Fig. 7 is a curved surface diagram of the displacement field in the x direction of the first phase of the inner wall of the left ventricle. the

图8是左心室内壁第一相位y方向位移场曲面图。 Fig. 8 is a curved surface diagram of the displacement field in the y direction of the first phase of the inner wall of the left ventricle. the

图9是左心室内壁第一相位z方向位移场曲面图。 Fig. 9 is a curved surface diagram of the displacement field in the z direction of the first phase of the inner wall of the left ventricle. the

图10是左心室内壁在7个相位处的应力场分布色图。 Figure 10 is a color map of the stress field distribution of the inner wall of the left ventricle at 7 phases. the

具体实施方式Detailed ways

下面结合附图对本发明作进一步描述。 The present invention will be further described below in conjunction with the accompanying drawings. the

参照图1～图10，一种基于四维医学图像的心壁应力应变测量方法，所述心壁应力应变测量方法包括以下步骤： Referring to Figures 1 to 10, a method for measuring heart wall stress and strain based on four-dimensional medical images, the method for measuring heart wall stress and strain includes the following steps:

rθ位移的拟合过程：使用[u_i v_i rθ_i]拟合出了场曲面Q(s，t)，根据B样条曲面特点Q(s，t)实际包括3个方程，即Q^u(s，t)＝u，Q^v(s，t)＝v，Q^θ(s，t)＝rθ；已知[θψR]＝P(u，v)，[u v rθ]＝Q(s，t)，定义P^-1和Q^-1分别为P(u，v)与Q(s，t)的逆运算；对模型上任意点a₀[θ₀ ψ₀ R₀]，计算其坐标值θ的位移量rθ₀的值，过程如下： The fitting process of rθ displacement: using [u _i v _i rθ _i ] to fit the field surface Q(s, t), according to the characteristics of the B-spline surface Q(s, t) actually includes three equations, that is, Q ^u (s, t)=u, Q ^v (s, t)=v, Q ^θ (s, t)=rθ; known [θψR]=P(u,v), [u v rθ]=Q(s, t), define P ^-1 and Q ^-1 as the inverse operation of P(u, v) and Q(s, t) respectively; for any point a ₀ [θ ₀ ψ ₀ R ₀ ] on the model, calculate its coordinate value The value of displacement rθ ₀ of θ, the process is as follows:

rψ位移的拟合过程：使用[u_i v_i rψ_i]拟合出了场曲面Q(s，t)，且有Q^u(s，t)＝u，^ψ(s，t)＝rψ；[θψ R ]＝P(u，v)，[u v rψ]＝Q(s，t)；设计算点a₀[θ₀ ψ₀ R₀]的rψ位移值，过程如下： The fitting process of rψ displacement: use [u _i v _i rψ _i ] to fit the field surface Q(s, t), and Q ^u (s, t)=u, ^ψ (s, t)=rψ; [θψ R ]=P(u,v), [u v rψ]=Q(s,t); design and calculate the rψ displacement value of point a ₀ [θ ₀ ψ ₀ R ₀ ], the process is as follows:

$\{\begin{matrix} {ϵ ϵ}_{x x} = = \frac{&PartialD; &PartialD; u u}{&PartialD; &PartialD; x x} = = \frac{{M m}^{x x} (({x x}_{00} + + dx dx,, {y the y}_{00},, {z z}_{00})) - - {M m}^{x x} (({x x}_{00},, {y the y}_{00},, {z z}_{00}))}{dx dx} = = \frac{{N N}^{θ θ} (({θ θ}_{00} + + dθ dθ,, {ψ ψ}_{00},, {R R}_{00})) - - {N N}^{θ θ} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00}))}{{R R}_{00} dθ dθ} \\ {ϵ ϵ}_{y the y} = = \frac{&PartialD; &PartialD; v v}{&PartialD; &PartialD; y the y} = = \frac{{M m}^{y the y} (({x x}_{00},, {y the y}_{00} + + dy dy,, {z z}_{00})) - - {M m}^{y the y} (({x x}_{00},, {y the y}_{00},, {z z}_{00}))}{dy dy} = = \frac{{N N}^{ψ ψ} (({θ θ}_{00},, {ψ ψ}_{00} + + dψ dψ,, {R R}_{00})) - - {N N}^{ψ ψ} (({θ θ}_{00} {,, ψ ψ}_{00},, {R R}_{00}))}{{R R}_{00} dψ dψ} \\ {ϵ ϵ}_{z z} = = \frac{&PartialD; &PartialD; w w}{&PartialD; &PartialD; z z} = = \frac{{M m}^{z z} (({x x}_{00},, {y the y}_{00},, {z z}_{00} + + dz dz)) - - {M m}^{z z} (({x x}_{00},, {y the y}_{00},, {z z}_{00}))}{dz dz} = = \frac{{N N}^{R R} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00} + + dR d)) - - {N N}^{R R} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00}))}{dR d} \\ {γ γ}_{xy xy} = = \frac{&PartialD; &PartialD; u u}{&PartialD; &PartialD; y the y} + + \frac{&PartialD; &PartialD; v v}{&PartialD; &PartialD; x x} = = \frac{{N N}^{θ θ} (({θ θ}_{00},, {ψ ψ}_{00} + + dψ dψ,, {R R}_{00})) - - {N N}^{θ θ} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00}))}{{R R}_{00} dψ dψ} + + \frac{{N N}^{ψ ψ} (({θ θ}_{00} + + dθ dθ,, {ψ ψ}_{00},, {R R}_{00})) - - {N N}^{ψ ψ} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00}))}{{R R}_{00} dθ dθ} \\ {γ γ}_{yz yz} = = \frac{&PartialD; &PartialD; v v}{&PartialD; &PartialD; z z} + + \frac{&PartialD; &PartialD; w w}{&PartialD; &PartialD; y the y} = = \frac{{N N}^{ψ ψ} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00} + + dR d)) - - {N N}^{ψ ψ} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00}))}{dR d} + + \frac{{N N}^{R R} (({θ θ}_{00},, {ψ ψ}_{00} + + dψ dψ,, {R R}_{00})) - - {N N}^{R R} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00}))}{{R R}_{00} dψ dψ} \\ {γ γ}_{xz xz} = = \frac{&PartialD; &PartialD; u u}{&PartialD; &PartialD; z z} + + \frac{&PartialD; &PartialD; w w}{&PartialD; &PartialD; x x} = = \frac{{N N}^{θ θ} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00} + + dR d)) - - {N N}^{θ θ} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00}))}{dR d} + + \frac{{N N}^{R R} (({θ θ}_{00} + + dθ dθ,, {ψ ψ}_{00},, {R R}_{00})) - - {N N}^{R R} (({θ θ}_{00},, {ψ ψ}_{00},, {R R}_{00}))}{{R R}_{00} dθ dθ} \end{matrix}$

(8) (8)

笛卡尔坐标系下的应变定义：物体位移向量δ由x，y，z轴上的位移u，v，w组成。应变向量ε有3个主应变和三个剪应变，与位移的关系为： Definition of strain in the Cartesian coordinate system: The object displacement vector δ is composed of displacements u, v, and w on the x, y, and z axes. The strain vector ε has three principal strains and three shear strains, and the relationship with displacement is:

应力向量σ定义如下： The stress vector σ is defined as follows:

本实施例的模型的建立过程为： The establishment process of the model of this embodiment is:

B样条方法及其逆运算：由一组控制点di(i＝0，1，2，...，n)生成的k次B样条曲线定义为： B-spline method and its inverse operation: the k-degree B-spline curve generated by a set of control points di (i=0, 1, 2, ..., n) is defined as:

式中u_i为一组非递减规范化序列，称作节点序列，通常取i＝0，1，2，...，n+k+1。同理由(m+1)×(n+1)个控制点生成的三维B样条曲面定义为： In the formula, u _i is a group of non-decreasing normalized sequences, called node sequences, usually i=0, 1, 2, ..., n+k+1. For the same reason, the three-dimensional B-spline surface generated by (m+1)×(n+1) control points is defined as:

其中，基函数的递推同(2)，k1、k2分别为u和v方向上的曲线阶数。样条曲面就是样条的推广，它将一组空间坐标下的点转换到了一个边长为1的平面空间u-v内，我们将u-v空间称为规范化空间。以上的定义是B样条方法的正向计算，即已知u，v求p(u，v)，在很多时候，我们需要在建立好的B样条模型P(u，v)上求的一点d[x，y，z]的参数u，v，我们称这样的计算为B样条曲面的逆运算，记为P^-1(x，y，z)，可知[u，v]＝P^-1(x，y，z)。P^-1(x，y，z)实际包含3个函数，即Px^-1(x)，Py^-1(y)和Pz^-1(z)，它们分别是取d_i，j的x，y，z轴分量为控制点P(u，v)的逆函数。对于逆函数的计算我们采用迭代的方法在u-v空间中进行搜索得到，迭代过程如图1。 Among them, the recursion of the basis function is the same as (2), and k1 and k2 are the curve orders in the u and v directions respectively. The spline surface is the extension of the spline, which transforms the points under a set of space coordinates into a plane space uv with a side length of 1. We call the uv space a normalized space. The above definition is the forward calculation of the B-spline method, that is, knowing u and v to find p(u, v), in many cases, we need to find it on the established B-spline model P(u, v) The parameters u, v of a point d[x, y, z], we call this calculation the inverse operation of the B-spline surface, denoted as P ^-1 (x, y, z), we know that [u, v] = P ^-1 (x, y, z). P ^-1 (x, y, z) actually contains 3 functions, namely Px ^-1 (x), Py ^-1 (y) and Pz ^-1 (z), which are x _{, y of d i, j} respectively , the z-axis component is the inverse function of the control point P(u, v). For the calculation of the inverse function, we use an iterative method to search in the uv space, and the iterative process is shown in Figure 1.

图1中ε为求解精度，m定义了逼近速度，Pd(i，j)定义如下： In Figure 1, ε is the solution accuracy, m defines the approach speed, and Pd(i, j) is defined as follows:

$Pd PD ((i i,, j j)) < < 00 &DoubleLeftRightArrow; &DoubleLeftRightArrow; \{\begin{matrix} [[{P P}^{x x} ((mu mu + + id id,, mv mv + + jd jd)) - - x x]] * * [[{P P}^{x x} ((mu mu + + ((i i + + 11)) d d,, mv mv + + ((j j + + 11)) d d)) - - x x]] < < 00 \\ [[{P P}^{y the y} ((mu mu + + id id,, mv mv + + jd jd)) - - y the y]] * * [[{P P}^{y the y} ((mu mu + + ((i i + + 11)) d d,, mv mv + + ((j j + + 11)) d d)) - - y the y]] < < 00 \\ [[{P P}^{z z} ((mu mu + + id id,, mv mv + + jd jd)) - - z z]] * * [[{P P}^{z z} ((mu mu + + ((i i + + 11)) d d,, mv mv + + ((j j + + 11)) d d)) - - z z]] < < 00 \end{matrix} - - - - - - ((44))$

结束时的mu+id和mu+jd即为所求u，v。 The mu+id and mu+jd at the end are the requested u and v. the

弹性力学方程及单元坐标变换：根据弹性力学知识在笛卡尔坐标系下定义物体位移向量δ由x，y，z轴上的位移u，v，w组成。应变向量ε有3个主应变和三个剪应变，与位移的关系为： Elastic mechanics equation and unit coordinate transformation: According to the knowledge of elastic mechanics, the displacement vector δ of the object is defined in the Cartesian coordinate system, which is composed of the displacement u, v, and w on the x, y, and z axes. The strain vector ε has three principal strains and three shear strains, and the relationship with displacement is:

应力向量σ定义如下： The stress vector σ is defined as follows:

其中，D为弹性矩阵，由研究对象的弹性模量E及泊松比μ确定，对于确定的材料可视为已知常量。由公式(5)和(6)可知求解应力应变的关键可归结于如何根据位移计算应变。 Among them, D is the elastic matrix, which is determined by the elastic modulus E and Poisson's ratio μ of the research object, and can be regarded as a known constant for a certain material. It can be seen from formulas (5) and (6) that the key to solving stress and strain can be attributed to how to calculate strain according to displacement. the

公式(5)和(6)的定义是在笛卡尔坐标系下的，也即研究对象的局部坐标，如图2(a)所示。由于我们的研究模型为壳状，不同位置点的局部坐标和全局坐标具有一定旋转和形变关系，如图2(b)。为了克服形状干扰，我们采用球面坐标系进行计算。笛卡尔坐标系与球面坐标系转化关系如下： The definitions of formulas (5) and (6) are in the Cartesian coordinate system, that is, the local coordinates of the research object, as shown in Figure 2(a). Since our research model is shell-shaped, the local coordinates and global coordinates of different positions have certain rotation and deformation relationships, as shown in Figure 2(b). In order to overcome shape interference, we use a spherical coordinate system for calculation. The conversion relationship between the Cartesian coordinate system and the spherical coordinate system is as follows:

$\{\begin{matrix} R R = = \sqrt{{x x}^{22} + + {y the y}^{22} + + {z z}^{22}} \\ θ θ = = arccos arccos ((\frac{z z}{R R})) \\ ψ ψ = = arccos arccos ((\frac{x x}{\sqrt{{x x}^{22} + + {y the y}^{22}}})) \end{matrix} - - - - - - ((77))$

球面坐标系对于曲面与平面的变形有很好的适应性。如图2(b)中所示，已知点a₀[θ₀ψ₀ R₀]的位移量分别记做N^θ(θ₀ ψ₀ R₀)、N^ψ(θ₀ ψ₀ R₀)和N^R(θ₀ ψ₀ R₀)，a₀在笛卡尔系下坐标为[x₀ y₀ z₀]，其位移量分别为M^x(x₀，y₀，z₀)，M^y(x₀，y₀，z₀)，M^z(x₀，y₀，z₀)，则根据公式(8)可直接由球面坐标系下位移求出应变 The spherical coordinate system has good adaptability to the deformation of curved surface and plane. As shown in Figure 2(b), the displacements of the known point a ₀ [θ ₀ ψ ₀ R ₀ ] are recorded as N ^θ (θ ₀ ψ ₀ R ₀ ), N ^ψ (θ ₀ ψ ₀ R ₀ ) and N ^R (θ ₀ ψ ₀ R ₀ ), the coordinates of a ₀ in the Cartesian system are [x ₀ y ₀ z ₀ ], and their displacements are M ^x (x ₀ , y ₀ , z ₀ ), M ^y (x ₀ , y ₀ , z ₀ ), M ^z (x ₀ , y ₀ , z ₀ ), then according to the formula (8), the strain can be obtained directly from the displacement in the spherical coordinate system

且γ_xy＝γ_yx，γ_yz＝γ_zy，γ_xz＝γ_zx。 And γ _xy =γ _yx , γ _yz =γ _zy , γ _xz =γ _zx .

根据位移场计算应变：已知B样条曲面模型P(u，v)，模型上存在n个位移已知标记点a_i(i＝1，2，...n)，a_i对应的规范化坐标为[u_i，v_i]，a_i的位移为[rθ_i rψ_i rR_i]，我们的目标是要通过离散的n个位移建立一个连续的位移场。对于一个没有出现断裂的变形体，连续位置的形变也必然是连续的，因此对于位移场我们也通过对离散点进行B样条拟合来获得。为了将位移场中的点与运动点的位置建立联系，我们对三个位移分别拟合。以rθ场为例，我们使用[u_i v_i rθ_i]拟合出了场曲面Q(s，t)，且有Q^u(s，t)＝u，Q^v(s，t)＝v，Q^θ(s，t)＝rθ。各个空间下对应点的坐标关系如图3。 Calculate the strain according to the displacement field: the B-spline surface model P(u, v) is known, and there are n mark points a _i (i=1, 2,...n) with known displacements on the model, and the normalization corresponding to a _i The coordinates are [u _i , v _i ], the displacement of a _i is [rθ _i rψ _i rR _i ], our goal is to establish a continuous displacement field through discrete n displacements. For a deformable body without fractures, the deformation at continuous positions must also be continuous, so we also obtain the displacement field by B-spline fitting to discrete points. To relate the points in the displacement field to the positions of the moving points, we fit the three displacements separately. Taking the rθ field as an example, we use [u _i v _i rθ _i ] to fit the field surface Q(s, t), and Q ^u (s, t)=u, Q ^v (s, t)=v , Q ^θ (s, t)=rθ. The coordinate relationship of corresponding points in each space is shown in Figure 3.

由图3可以清楚的看到[θψ R]＝P(u，v)，[u v rθ]＝Q(s，t)。假设我们要计算点d₀[θ₀ψ₀ R₀]的rθ位移值，过程如下： It can be clearly seen from Fig. 3 that [θψ R]=P(u,v), [u v rθ]=Q(s,t). Suppose we want to calculate the rθ displacement value of point d ₀ [θ ₀ ψ ₀ R ₀ ], the process is as follows:

Step1：将坐标转换至规范坐标系上，[u₀ v₀]＝P^-1(θ₀ ψ₀ R₀) ； Step1: Convert the coordinates to the standard coordinate system, [u ₀ v ₀ ]=P ^-1 (θ ₀ ψ ₀ R ₀ );

Step3：得到位移rθ，rθ＝Q^θ(s₀，t₀)。 Step3: Obtain the displacement rθ, rθ=Q ^θ (s ₀ , t ₀ ).

通过这种变换我们可以得到心壁上任意点的位移向量，我们将这种计算记做Trans^θ(θ ψ R)，相应的还有Trans^ψ(θ ψ R)，Trans^R(θ ψ R)。通过以上计算可以得到全局坐标下d₀的位置及对应位移，接着应用公式(9)将位置及位移旋转至局部坐标下，最后用公式(8)和(6)分别计算应变和应力。 Through this transformation, we can get the displacement vector of any point on the heart wall. We record this calculation as Trans ^θ (θ ψ R), and there are correspondingly Trans ^ψ (θ ψ R), Trans ^R (θ ψ R) . Through the above calculations, the position and corresponding displacement of d ₀ in the global coordinates can be obtained, and then the position and displacement are rotated to the local coordinates by applying formula (9), and finally the strain and stress are calculated by formulas (8) and (6), respectively.

本实施例的模型及位移场：本文实验的数据源选择左心室的SPECT图像，这是由于在SPECT图像中物理点和标记点很清楚的可以看到，大量的物理点能在很短的时间间隔内被标记，而且可以在心脏周期内被全程跟踪。选择左心室内外壁共1165个标记点(其中内壁355个，外壁810个)，每隔100ms跟踪记录它们的位置，获得一个心动周期内7个相位的标记点位置，记做D_i(i＝1，2，...，7)。根据公式(7)，将D_i转化至球坐标下记做D_i’，并以D_i’为源应用NURBS方法拟合左心室7个时刻的内外面模型，模型如图4所示。 The model and displacement field of this embodiment: The data source of the experiment in this paper is the SPECT image of the left ventricle. This is because the physical points and marker points can be clearly seen in the SPECT image, and a large number of physical points can be displayed in a short time. Intervals are marked and can be tracked throughout the cardiac cycle. Select a total of 1165 marker points on the inner and outer walls of the left ventricle (including 355 on the inner wall and 810 on the outer wall), track and record their positions every 100 ms, and obtain the positions of the marker points of 7 phases in one cardiac cycle, which are recorded as D _i (i= 1, 2, ..., 7). According to the formula (7), D _i was transformed into spherical coordinates and denoted as D _i ', and using D _i ' as the source, the NURBS method was used to fit the inner and outer models of the left ventricle at 7 moments. The model is shown in Figure 4.

通过做差运算可以找到标记点在相邻时刻的位移为D_i+1’-D_i’，图5为第一相位到第二相位的位移矢量图。该矢量图直观的反应了左心室心室壁的运动方向及大小，与图6的图形进行比较发现位移的总体方向呈螺旋扭曲形，非常的相近也验证了标记点位移的正确性。根据位移及规范化坐标我们对位移进行拟合得到了各个时刻的位移场，图7、图8和图9为第一相位左心室内壁的位移场曲面图。从位移场中我们可以明显的看出在左心室运动时，心尖处的位移大于心底，而左心室后部的位移大于前部。 The displacement of the marking point at adjacent moments can be found by performing difference operation as D _i+1 '-D _i ', and Fig. 5 is the displacement vector diagram from the first phase to the second phase. The vector diagram intuitively reflects the direction and size of the movement of the left ventricle wall. Compared with the graph in Figure 6, it is found that the overall direction of the displacement is spiral twisted, which is very similar and verifies the correctness of the displacement of the marker points. According to the displacement and normalized coordinates, we fitted the displacement to obtain the displacement field at each moment. Figure 7, Figure 8 and Figure 9 are the surface diagrams of the displacement field of the inner wall of the left ventricle in the first phase. From the displacement field, we can clearly see that when the left ventricle moves, the displacement at the apex is greater than that at the bottom of the heart, and the displacement at the rear of the left ventricle is greater than that at the front.

众所周知对于心脏的研究无法以实体实验的方式进行验证，因此我们将结果与其他方法获得的结果进行比较。现有技术通过有限元方法建立模型对左心室形变进行研究得到了结论，即心室扭曲运动，心尖处绕心底旋转，心室内壁扭曲程度大于外壁，这都与我们得到的左心室位移场相吻合，从而验证了位移场的正确性。 Heart studies are notoriously incapable of being validated experimentally, so we compared the results with those obtained by other methods. In the prior art, the research on the deformation of the left ventricle by establishing a model with the finite element method has drawn conclusions, that is, the twisting movement of the ventricle, the rotation of the apex around the bottom of the heart, and the twisting degree of the inner wall of the ventricle is greater than that of the outer wall, which is consistent with the displacement field of the left ventricle we obtained. Thus the correctness of the displacement field is verified. the

位移场同时体现出了左心室室壁的形变特点，对它进行定量分析可以的到一些标志性的数据。我们通过位移场对心动一周期内心室长轴与内壁半径进行分析，得到图7、8、9所示结果，与现有技术得到的结果进行比较，可以看到变化方向与幅度都大体相同。 The displacement field also reflects the deformation characteristics of the left ventricle wall, and some landmark data can be obtained by quantitative analysis of it. We analyzed the long axis and inner wall radius of the cardiac chamber through the displacement field, and obtained the results shown in Figures 7, 8, and 9. Compared with the results obtained by the prior art, we can see that the direction and magnitude of the changes are roughly the same. the

应力应变计算：分析左心室的应变量具有十分重要的意义，因为它可以定量来研究心肌性能，它可以将心肌的局部变形分离出来。由于采用不同的标记成像技术，标记结果会有一定得差异，导致数据结果也存在一定得差异，在本文中我们将所得结果与tagged MRI技术所得得研究结果进行对比，该结果得到了广泛的认可。通过对比，我们得到的应变场与现有技术得到的结论一致；即心尖到心底应变逐渐减小，心后应变大于心前。 Stress-strain calculation: It is of great significance to analyze the strain of the left ventricle, because it can be used to quantitatively study the performance of the myocardium, and it can separate the local deformation of the myocardium. Due to the use of different tagging imaging techniques, there will be some differences in the tagging results, resulting in some differences in the data results. In this paper, we compared the results obtained with the research results obtained by tagged MRI technology, and the results have been widely recognized. . By comparison, the strain field obtained by us is consistent with the conclusion obtained by the prior art; that is, the strain from the apex to the bottom of the heart gradually decreases, and the strain behind the heart is greater than that in the front of the heart. the

左心室壁的应力场计算具有重要意义，心壁应力是心肌冠脉血流和耗氧情况的主要决定因素之一。根据应变场，弹性模量与泊松比分别取E＝，λ＝，计算心室内外壁的应力值。图10所示为左心室内壁在7个相位处的应力场分布色图，这里的应力为根据米责思约束条件计算的到得统一值，色彩的变化根据应力增大而逐渐变暖，到红色处达到最大，该结果与现有技术得到的结果相一致，即总体分布不均匀，心室内壁应力大于外壁，在心尖和心底部应力较大，心尖处还出现应力集中现象。 The calculation of the stress field of the left ventricular wall is of great significance, and the stress of the heart wall is one of the main determinants of myocardial coronary blood flow and oxygen consumption. According to the strain field, the elastic modulus and Poisson's ratio are respectively E=, λ=, and the stress values of the inner and outer walls of the ventricle are calculated. Figure 10 shows the color map of the stress field distribution of the inner wall of the left ventricle at seven phases. The stress here is a uniform value calculated according to the Mies constraints. The red point reaches the maximum. This result is consistent with the results obtained by the prior art, that is, the overall distribution is uneven, the stress on the inner wall of the ventricle is greater than that on the outer wall, the stress on the apex and the bottom of the heart is greater, and stress concentration also occurs at the apex. the

Claims

1. cardiac wall stress-strain measuring method based on four-dimensional medical image, said cardiac wall stress-strain measuring method may further comprise the steps:

1), with heart to be detected set up B-spline surface model P (u, v), the standardization parameter when u wherein, v are the B-spline surface definition, parameter field is [0,1], establishes to have n the known gauge point a of displacement on the model _i(i=1,2 ... n), a _iThe locus under spheric coordinate system be [θ _iψ _iR _i], (u, the standardization coordinate of v) going up is [u at P _i, v _i], and a _iDisplacement be [r ^{θ i}r ^{Ψ irRi}], set up the r of The model through n the discrete displacement vector of putting ^θ, r ^ψAnd r ^RThe continuous dislocation field;

2), utilize standardization coordinate u, v match displacement field is also set up corresponding relation with point in the displacement field and actual spatial point, to three displacements match respectively, wherein:

The fit procedure of r θ displacement: use [u _iv _ir ^{θ i}] (s, t), (s t) actually comprises 3 equations, i.e. Q according to B-spline surface characteristics Q to have simulated curved surface Q ^u(s, t)=u, Q ^v(s, t)=v, Q ^θ(s, t)=r θ; Known [θ Ψ R]=P (u, v), [u v r θ]=Q (s, t), definition P ^-1And Q ^-1(u is v) with Q (s, inverse operation t) to be respectively P; To arbitrfary point a on the model ₀[θ ₀ψ ₀R ₀], calculate the displacement r θ of its coordinate figure θ ₀Value, process is following:

Step1: on Coordinate Conversion to standard coordinate system, [u ₀v ₀]=P ^-1(θ ₀ψ ₀R ₀);

Step2: calculate this coordinate s in displacement field ₀, t ₀, [s ₀t ₀]=Q ^-1(u ₀v ₀);

Step3: obtain displacement r θ, r θ=Q ^θ(s ₀, t ₀);

The fit procedure of r ψ displacement: use [u _iv _iR Ψ _i] (s t), and has Q to have simulated a curved surface Q ^u(s, t)=u, Q ^v(s, t)=v, Q ^Ψ(s, t)=r Ψ; [θ Ψ R]=P (u, v), [uvr ψ]=Q (s, t); If calculation level a ₀[θ ₀ψ ₀R ₀] r Ψ shift value, process is following:

Step1: on Coordinate Conversion to standard coordinate system, [u ₀v ₀]=P ^-1(ψ ₀ψ ₀R ₀);

Step3: obtain displacement r Ψ, r Ψ=Q ^ψ(s ₀, t ₀);

r ^RThe fit procedure of displacement: use [u _iv _ir ^Ri] (s t), and has Q to have simulated a curved surface Q ^u(s, t)=u, Q ^v(s, t)=v, Q ^R(s, t)=r ^R[θ Ψ R]=P (u, v), [u v rR]=Q (s, t); If calculation level α ₀[θ ₀ψ ₀R ₀] r ^RShift value, process is following:

Step1: on Coordinate Conversion to standard coordinate system, [u ₀v ₀]=P ^-1(R ₀ψ ₀R ₀);

Step3: obtain displacement r ^R, r ^R=Q ^R(s ₀, t ₀);

3), with step 2) arbitrfary point a on the center wall ₀[θ ₀ψ ₀R ₀] displacement remember respectively and be N ^θ(θ ₀ψ ₀R ₀), N ^ψ(θ ₀ψ ₀R ₀) and N ^R(θ ₀ψ ₀R ₀); Because strain is at cartesian coordinate system x, y, z is undefined, and we need this definition is converted under the spheric coordinate system, establish a ₀Descartes is that coordinate is [x down ₀y ₀z ₀], its displacement is respectively M ^x(x ₀, y ₀, z ₀), M ^y(x ₀, y ₀, z ₀), M ^z(x ₀, y ₀, z ₀), then directly obtain strain, and obtain stress according to (6) by the spheric coordinate system bottom offset according to formula (8):

\{\begin{matrix} ϵ_{x} = \frac{&PartialD; u}{&PartialD; x} = \frac{M^{x} (x_{0} + dx, y_{0}, z_{0}) - M^{x} (x_{0}, y_{0}, z_{0})}{dx} = \frac{N^{θ} (θ_{0} + dθ, ψ_{0}, R_{0}) - N^{θ} (θ_{0}, ψ_{0}, R_{0})}{R_{0} dθ} \\ ϵ_{y} = \frac{&PartialD; v}{&PartialD; y} = \frac{M^{y} (x_{0}, y_{0}, + dy, z_{0}) - M^{y} (x_{0}, y_{0}, z_{0})}{dy} = \frac{N^{ψ} (θ_{0}, ψ_{0} + dψ, R_{0}) - N^{ψ} (θ_{0}, ψ_{0}, R_{0})}{R_{0} dψ} \\ ϵ_{z} = \frac{&PartialD; w}{&PartialD; z} = \frac{M^{z} (x_{0}, y_{0}, z_{0} + dz) - M^{z} (x_{0}, y_{0}, z_{0})}{dz} = \frac{N^{R} (θ_{0}, ψ_{0}, R_{0} + dR) - N^{R} (θ_{0}, ψ_{0}, R_{0})}{dR} \\ γ_{xy} = \frac{&PartialD; u}{&PartialD; y} + \frac{&PartialD; v}{&PartialD; x} = \frac{N^{θ} (θ_{0}, ψ_{0} + dψ, R_{0}) - N^{θ} (θ_{0}, ψ_{0}, R_{0})}{R_{0} dψ} + \frac{N^{ψ} (θ_{0} + dθ, ψ_{0}, R_{0}) - N^{ψ} (θ_{0}, ψ_{0}, R_{0})}{R_{0} dθ} \\ γ_{yz} = \frac{&PartialD; v}{&PartialD; z} + \frac{&PartialD; w}{&PartialD; y} = \frac{N^{ψ} (θ_{0}, ψ_{0}, R_{0} + dR) - N^{ψ} (θ_{0}, ψ_{0}, R_{0})}{dR} + \frac{N^{R} (θ_{0}, ψ_{0} + dψ, R_{0}) - N^{R} (θ_{0}, ψ_{0}, R_{0})}{R_{0} dψ} \\ γ_{xz} = \frac{&PartialD; u}{&PartialD; z} + \frac{&PartialD; w}{&PartialD; x} = \frac{N^{θ} (θ_{0}, ψ_{0}, R_{0} + dR) - N^{θ} (θ_{0}, ψ_{0}, R_{0})}{dR} + \frac{N^{R} (θ_{0} + dθ, ψ_{0}, R_{0}) - N^{R} (θ_{0}, ψ_{0}, R_{0})}{R_{0} dθ} \end{matrix} - - - (8)

Strain under cartesian coordinate system definition: ohject displacement vector δ is by x, y, and the displacement components u on the z axle, v, w forms.Answer variable vector ε that 3 principal strains and three shearing strains are arranged, be with the relation of displacement:

\{\begin{matrix} ϵ_{x} = \frac{&PartialD; u}{&PartialD; x} \\ ϵ_{y} = \frac{&PartialD; v}{&PartialD; y} \\ ϵ_{z} = \frac{&PartialD; w}{&PartialD; z} \\ γ_{xy} = \frac{&PartialD; u}{&PartialD; y} + \frac{&PartialD; v}{&PartialD; x} \\ γ_{yz} = \frac{&PartialD; v}{&PartialD; z} + \frac{&PartialD; w}{&PartialD; y} \\ γ_{xz} = \frac{&PartialD; u}{&PartialD; z} + \frac{&PartialD; w}{&PartialD; x} \end{matrix} - - - (5)

Stress vector σ defines as follows:

{σ}＝[σ _xσ _yσ _zτ _xyτ _yzτ _xz]＝Dε(6)

Wherein, D is an elastic matrix, and is definite by the elastic modulus E and the Poisson's ratio μ of object of study, can be considered known constant for the material of confirming.

2. the cardiac wall stress-strain measuring method based on four-dimensional medical image as claimed in claim 1 is characterized in that:

In said step 1), B-spline surface modelling process: by one group of control point di (i=0,1,2 ..., k B-spline curves that n) generate are defined as:

Wherein, N _{I, k}(u) be the B spline base function, its classical de-Boor Cox definition is one group of recurrence formula, that is:

U in the formula _iBe one group of non-decreasing standardization sequence, be called sequence node, get i=0 usually, 1,2 ..., n+k+1; The three-dimensional B-spline surface that the individual control point of same reason (m+1) * (n+1) generates is defined as:

P (u, v) = Σ_{i = 0}^{m} Σ_{j = 0}^{n} d_{i, j} N_{i, k 1} (u, v) N_{j, k 2} (u, v), u_{i} \leq u < u_{i + 1}, v_{j} \leq v < v_{j + 1} - - - (3)

Wherein, k1, k2 are respectively the curve exponent number on u and the v direction, and spline surface is exactly the popularization of batten, and it has been transformed into the point under one group of space coordinates in one 1 * 1 the space u-v, is normalized space with the u-v definition space.