Spatiotemporal finite element method for motion analysis with velocity data

• US 6,236,738 B1
• Filed: 10/01/1998
• Issued: 05/22/2001
• Est. Priority Date: 04/09/1998
• Status: Expired due to Term

• Current Assignee: Board of Trustees of the Leland Stanford Junior University (Stanford University)
• Original Assignee: Board of Trustees of the Leland Stanford Junior University (Stanford University)
• Inventors: Zhu, Yudong, Pelc, Norbert J.
• Primary Examiner(s): Johns, Andrew W.
• Assistant Examiner(s): NAKHJAVAN, SHERVIN K
• Application Number:
  US09/165,232
• Time in Patent Office: 964 Days
• Field of Search: 382/100, 382/103, 382/107, 348/699
• US Class Current: 382/107
• CPC Class Codes: G06T 7/215 Motion-based segmentation
• CPC Class Codes: G06T 7/246 using feature-based methods, e.g. the tracking of corners or segments

Abstract

Disclosed is a method for nonrigid cyclic motion analysis using a series of images covering the cycle, acquired, for example, from phase contrast magnetic resonance imaging. The method is based on fitting a global spatiotemporal finite element mesh model to motion data samples of an extended region at all time frames. A spatiotemporal model is composed of time-varying finite elements, with the nonrigid motion of each characterized by a set of Fourier harmonics. The model is suitable for accurately modeling the kinematics of a cyclically moving and deforming object with complex geometry, such as that of the myocardium. The model has controllable built-in smoothing in space and time for achieving satisfactory reproducibility in the presence of noise. Motion data measured, with PC MRI for example, can be used to quantify motion and deformation by fitting the model to data.

Claims

1. A method of motion analysis of a moving region within an object comprising the steps of:

a) defining a dynamic mesh within said region, said mesh composed of non-overlapping elements each specified by a set of node points,

b) acquiring motion data for a set of locations within said region for a plurality of time frames, wherein at least some of the locations for which motion data are acquired do not correspond to a node point of any mesh element,

c) defining a spatiotemporal model of the region which relates the kinematics of the node points to the acquired motion data, and

d) computing the motion and deformation of the moving region using the motion data and said model.

2. The method as defined by claim 1 wherein said elements are two dimensional.

3. The method as defined by claim 2 wherein said elements are triangles.

4. The method as defined by claim 1 wherein said elements are three-dimensional.

5. The method as defined by claim 4 wherein said elements are tetrahedrons.

6. The method as defined by claim 1 wherein in step c) the motion is modeled as being periodic with period 1/fo and the trajectory of a node point is modeled as a Fourier series:

${\stackrel{\to <br><br>}{x}}^e<br><br>\left({\stackrel{\to <br><br>}{\mathrm{\xi <br><br>}}}_a,t\right)=\underset{k=-M}{\overset{m}{\sum <br><br>}}<br><br>{\stackrel{\to <br><br>}{g}}_a^e<br><br>\left(k\right)<br><br>{\mathrm{<br><br>}}^{\mathrm{j2}<br><br>\text{\hspace{1em}<br><br>}<br><br>\mathrm{\pi <br><br>}<br><br>\text{\hspace{1em}<br><br>}<br><br>{\mathrm{kf}}_0<br><br>t}$

where the superscript e is the index of the element containing the node point, a is the index of the node point, {right arrow over (x)}^e({right arrow over (ξ

)}_a,t) is the trajectory of the node point,

t is time within the motion period,

{right arrow over (g)}_a^e is a Fourier coefficient,

and k is the index of the Fourier series.

7. The method as defined by claim 6 wherein step d) includes a least squares fit for the coefficients {right arrow over (g)}_a^e(k).

8. The method as defined by claim 7 wherein a fitting residual is used to estimate uncertainty in motion data as follows:

${\widehat{\mathrm{\sigma <br><br>}}}_{\mathrm{\upsilon <br><br>}}^2=\frac{\mathrm{RSS}}{\mathrm{degrees}<br><br>\text{\hspace{1em}<br><br>}<br><br>\mathrm{of}<br><br>\text{\hspace{1em}<br><br>}<br><br>\mathrm{freedom}<br><br>\text{\hspace{1em}<br><br>}<br><br>\text{\hspace{1em}<br><br>}<br><br>\mathrm{for}<br><br>\text{\hspace{1em}<br><br>}<br><br>\mathrm{residue}}$

where σ

_v² is the estimated variance of the motion data and RSS is the sum of squared differences between the measured motion data and the modeled motion data.

9. The method of claim 1 wherein the motion data are velocity data.

10. The method of claim 9 wherein the velocity data are measured with phase contrast magnetic resonance imaging.

11. The method as defined by claim 1 wherein steps c) and d) are executed iteratively.

12. The method as defined by claim 1 wherein in step c) the trajectory of a specific point within an element is related to the trajectories of the nodes of that element by:

${\stackrel{\to <br><br>}{x}}^e<br><br>\left(\stackrel{\to <br><br>}{\mathrm{\xi <br><br>}},t\right)=\underset{a=1}{\overset{n_{\mathrm{en}}}{\sum <br><br>}}<br><br>N_a<br><br>\left(\stackrel{\to <br><br>}{\mathrm{\xi <br><br>}}\right)<br><br>{\stackrel{\to <br><br>}{x}}^e<br><br>\left({\stackrel{\to <br><br>}{\mathrm{\xi <br><br>}}}_a,t\right)$

where

the superscript e is the index of the element, a is the index of a node point,

n_en is the number of nodes of the eth mesh element,

{right arrow over (ξ

)} is the local coordinates of the specific point,

N_a is the shape function associated with the ath node, t is time,

{right arrow over (x)}^e({right arrow over (ξ

)}, t) is the trajectory of the specific point, and

{right arrow over (x)}^e({right arrow over (ξ

)}_a, t) is the trajectory of the ath node of the eth element.

13. The method as defined by claim 12 wherein in step c) the velocity of the specific point is related to the velocity of the nodes by:

${\stackrel{\to <br><br>}{\mathrm{\upsilon <br><br>}}}^e<br><br>\left(\stackrel{\to <br><br>}{\mathrm{\xi <br><br>}},t\right)=\underset{a=1}{\overset{n_{\mathrm{en}}}{\sum <br><br>}}<br><br>N_a<br><br>\left(\stackrel{\to <br><br>}{\mathrm{\xi <br><br>}}\right)<br><br>{\stackrel{\to <br><br>}{\mathrm{\upsilon <br><br>}}}^e<br><br>\left({\stackrel{\to <br><br>}{\mathrm{\xi <br><br>}}}_a,t\right)$

where

{right arrow over (v)}^e({right arrow over (ξ

)}, t) is the velocity of the specific point,

and {right arrow over (v)}^e({right arrow over (ξ

)}_a, t) is the velocity of the ath node of the eth element.

14. The method as defined by claim 13 wherein in step c) the motion is modeled as being periodic with period 1/fo, and wherein the spatiotemporal finite element mesh model is:

${\stackrel{\to <br><br>}{\mathrm{\upsilon <br><br>}}}^e<br><br>\left(\stackrel{\to <br><br>}{\mathrm{\xi <br><br>}},t\right)=\underset{a=1}{\overset{n_{\mathrm{en}}}{\sum <br><br>}}<br><br>N_a<br><br>\left(\stackrel{\to <br><br>}{\mathrm{\xi <br><br>}}\right)<br><br>\underset{k=-M,k\ne <br><br>0}{\overset{M}{\sum <br><br>}}<br><br>\left(\mathrm{j2}<br><br>\text{\hspace{1em}<br><br>}<br><br>\mathrm{\pi <br><br>}<br><br>\text{\hspace{1em}<br><br>}<br><br>{\mathrm{kf}}_0\right)<br><br>{\stackrel{\to <br><br>}{g}}_a^e<br><br>\left(k\right)<br><br>{\mathrm{<br><br>}}^{\mathrm{j2}<br><br>\text{\hspace{1em}<br><br>}<br><br>\mathrm{\pi <br><br>}<br><br>\text{\hspace{1em}<br><br>}<br><br>{\mathrm{kf}}_0<br><br>t}$

where k is the index of a Fourier series.

15. The method as defined by claim 12 wherein in step c) the motion is modeled as being periodic with period 1/fo, and wherein step d) includes processing displacement data through displacement field modeling as follows:

$\mathrm{\Delta <br><br>}<br><br>\text{\hspace{1em}<br><br>}<br><br>{\stackrel{\to <br><br>}{x}}^e<br><br>\left(\stackrel{\to <br><br>}{\mathrm{\xi <br><br>}},t\right)=\underset{a=1}{\overset{n_{\mathrm{en}}}{\sum <br><br>}}<br><br>N_a<br><br>\left(\stackrel{\to <br><br>}{\mathrm{\xi <br><br>}}\right)<br><br>\underset{k=-M}{\overset{M}{\sum <br><br>}}<br><br>{\stackrel{\to <br><br>}{g}}_a^e<br><br>\left(k\right)<br><br>\left({\mathrm{<br><br>}}^{\mathrm{j2}<br><br>\text{\hspace{1em}<br><br>}<br><br>\mathrm{\pi <br><br>}<br><br>\text{\hspace{1em}<br><br>}<br><br>{\mathrm{kf}}_0<br><br>t}-1\right)$

where Δ

{right arrow over (x)}^e is the displacement of a specific point.

16. The method as defined by claim 12 wherein step d) includes calculation of a deformation gradient tensor F in each element as follows:

${\hat{F}}_{\mathrm{iJ}}<br><br>\left(\stackrel{\to <br><br>}{X},t\right)=\frac{\partial <br><br>{\hat{x}}_i}{\partial <br><br>X_J}=\underset{a=1}{\overset{n_{\mathrm{en}}}{\sum <br><br>}}<br><br>{\hat{x}}_i^e<br><br>\left({\stackrel{\to <br><br>}{\mathrm{\xi <br><br>}}}_a,t\right)<br><br>\left(\underset{j=1}{\overset{n_{\mathrm{sd}}}{\sum <br><br>}}<br><br>\frac{\partial <br><br>N_a}{\partial <br><br>{\mathrm{\xi <br><br>}}_j}<br><br>\frac{\partial <br><br>{\mathrm{\xi <br><br>}}_j}{\partial <br><br>X_J}\right)$

wherein n_sd denotes the number of spatial dimensions.

17. The method as defined by claim 16 wherein step d) includes calculation of a Lagrangian finite strain tensor E as a function of deformation gradient as follows:

18. The method as defined by claim 12 wherein step d) includes calculation of a displacement field as follows:

$\widehat{\stackrel{\to <br><br>}{x}}<br><br>\left(\stackrel{\to <br><br>}{X},t\right)={\widehat{\stackrel{\to <br><br>}{x}}}^e<br><br>\left(\stackrel{\to <br><br>}{\mathrm{\xi <br><br>}},t\right)=\underset{a=1}{\overset{n_{\mathrm{en}}}{\sum <br><br>}}<br><br>N_a<br><br>\left(\stackrel{\to <br><br>}{\mathrm{\xi <br><br>}}\right)<br><br>{\widehat{\stackrel{\to <br><br>}{x}}}^e<br><br>\left({\stackrel{\to <br><br>}{\mathrm{\xi <br><br>}}}_a,t\right).$

19. The method of claim 1 wherein the motion data are displacement data.

20. The method of claim 19 wherein the displacement data are measured using magnetic resonance tagging.

21. The method as defined by claim 1 wherein step d) uncertainty in motion data is used to estimate uncertainty in measured kinematics of the region.