Proc SPIE Int Soc Opt Eng 7259: 72590V

PMID: 19946575

Shortest Path Refinement For HARP Motion Tracking

1. INTRODUCTION

Magnetic resonance (MR) imaging is capable of directly imaging motion of soft tissues. Traditional MR images carry information about motion only at the boundaries of tissues because the image intensities of the interior of the tissues is pretty homogeneous. MR tagging¹2 places non-invasive and temporary markers (tags) inside the soft tissues in a pre-specified pattern, yielding images that carry information about motion within homogeneous tissues. This complements traditional anatomical images, and enables imaging the detailed information about the motion of tissues such as the heart and the tongue throughout the time. Displacement, velocity, rotation, elongation, strain, and twist are just some of the quantities that be computed from it.

As one of the methods for estimating motion from tagged MR images, the harmonic phase (HARP) method³5 has proved to be valuable for both scientific and clinical applications. HARP tracking implicitly assumes that tissue points do not move much from one time to the next. If the tissue moves too fast, or the temporal resolution is too low, or the MR tag imaging parameters are not selected correctly, this assumption is violated and HARP tracking will fail. Manual correction is required when a point is mistracked, and this can be very time-consuming. There have been some previous e orts to identify and automatically correct mistracked points. Khalifa et al.⁶ used an active contour model to correct HARP tracking for cardiac motion. The approach is limited to the circular geometry, however, and is not easily generalized for other applications for example the tongue. Tecelao et al.⁷ proposed an extended HARP tracking method to correct the mistracking caused by through-plane motion and near the cardiac boundaries, but it did not completely address the mistracking problem. Use of spatial continuity of motion, generically called refinement, was described in Osman et al.³ as a process that could be employed to alleviate the mistracking described above. At the time, however, it was thought to be overly time consuming to employ on a routine basis. It was also developed on circular geometries, and is not straightforward to extend to arbitrary tissue points. Recently, we reported a region growing refinement method,⁸ which automatically tracks the whole tissue after a seed point is manually specified, and does this extremely rapidly. It has been successfully applied to both tongue and cardiac motion tracking. Though applied in order to compute motion, the region-growing refinement process can also be thought of as an application-specific harmonic phase unwrapping process. With this interpretation in hand, it is clear that the flood-fill algorithm used to unwrap DENSE phase images (cf.⁹), reported in Spottiswoode et al.,¹⁰ is another example of a refinement algorithm for improved motion estimation.

In both the region growing HARP refinement method and the flood-fill algorithm in DENSE phase unwrapping method, tissue points are tracked or phase-unwrapped in an order based solely on the local smoothness of the phase image. Because of image noise and the existence of other tissues besides that of primary interest (e.g., the heart or tongue), the paths from the seed to any given points of interest can be quite erratic, and incorrect tracking may result. Though refinement is still a highly desirable strategy, it is clear better results should be expected if an optimal path from the seed to each point can be found.

In this paper, we propose a new refinement method that represents the image as a graph, and solves for the optimal refinement path by solving a shortest path problem. In this way, each tissue point will be accessed—i.e., tracked or phase unwrapped—via an optimal path from the seed point (which is assumed to be correctly tracked). Cost functions defined on both the edges and vertices of the (image-based) graph encourage the shortest path to stay within the same tissue region as the seed point so that incorrect estimation from long paths in adjacent regions are much less likely. The proposed method is shown to be very robust and fast. Though not demonstrated herein, the method should be applicable to DENSE phase unwrapping and motion tracking as well.

2. BACKGROUND

2.1 MRI tagging and HARP method

The spatial modulation of magnetization (SPAMM) method¹2 is a popular MR tagging protocol. As the most basic form of SPAMM, the 1–1 SPAMM generates smoothly varying sinusoidal tag pattern in the tissue using two 90° RF pulses. With a [+90°, +90°] RF pulse pair, the acquired image can be described as:

A(x, t) = M₀(1 − e^{+ cos(}ϕ(x, t))e^),

(1)

where M₀ is the intensity of underlying MR image, e^{represents the signal decay due to T1 recovery, and}ϕ(x, t) is the tagging phase. At the time t₀ immediately after the tagging is applied, the tagging phase is a linear function of the point coordinate x:

ϕ(x, t₀) = kx ⋅ n + ϕ₀,

(2)

where k is the tagging frequency, n is the normal unit vector that represents the tagging direction and ϕ₀ is the phase o set. When using a [+90°, −90°] RF pulse pair, the image is:

B(x, t) = M₀(1 − e^{− cos(}ϕ(x, t))e^).

(3)

In order to improve the tag contrast and tag persistency, these two images can be combined together to generate the CSPAMM (¹¹) or MICSR (¹²) images:

\begin{matrix} I_{CSPAMM} = A - B = 2 M_{0} e^{- t ∕ T_{1}}) \cos (ϕ (x, t)), \\ I_{MICSR} = A^{2} - B^{2} = 4 M_{0}^{2} (1 - e^{- t ∕ T_{1}}) e^{- t ∕ T_{1}} \cos (ϕ (x, t)) . \end{matrix}

(4)

An example of a SPAMM cardiac image is shown in Fig. 1. The Fourier transform of a SPAMM image has two harmonic peaks. To estimate the phase value ϕ(x, t) from the tagged images, the HARP method applies a bandpass filter in the Fourier domain to extract just one of the harmonic peaks. The resulting filtered complex image is given by³

I(x, t) = D(x, t)e^,

(5)

where D(x, t) is called the harmonic magnitude image, and ϕ(x, t) is the harmonic phase (HARP) image. The magnitude image reflects the tissue anatomy, and the HARP image contains the tissue motion information. Because the HARP phase must be computed using an arctangent operation, its value θ(x, t) is the principal value of the true phase. It is therefore restricted to take on values in the interval [−π, +π), and is related to the true phase by

θ(x, t) = W(ϕ(x, t))

(6)

where W (·) is the wrapping function defined as:

W(ϕ) = mod(ϕ + π, 2π) − π

(7)

The true phase and harmonic phase are both material properties of tissue point. Thus, the HARP values of a material point do not change as the point moves around in space. This property is called phase invariance and is the basis of HARP motion tracking.

An external file that holds a picture, illustration, etc.
Object name is nihms-158349-f0001.jpg

Figure 1

(a) Part of a SPAMM tagged MR image, (b) the magnitude of its Fourier transform, and (c) the HARP image after applying the bandpass filter as in (b).

2.2 HARP motion tracking

The HARP tracking method is described in detail in³ and we give only a brief description herein. The HARP image θ(x, t) contains information about tissue motion in the tagging direction n. To track the 2D apparent motion of a material point, one needs two HARP images which are generally acquired with orthogonal tagging directions. Let Φ = [ϕ₁, ϕ₂]^{be the tagging phases of the two orthogonally tagged images. For a material point located at}x_{t_i} at time frame t_i, we look for a point x_{t_i+1} at time t_i+1 such that

Φ(x_{t_i+1}, t_i+1) = Φ(x_{t_i}, t_i)

(8)

This is a multidimensional root finding problem, and can be solved iteratively using the Newton-Raphson technique:

x_{t_{i + 1}}^{(n + 1)} = x_{t_{i + 1}}^{(n)} - \nabla Φ {(x_{t_{i + 1}}^{(n)}, t_{i + 1})}^{- 1} (Φ (x_{t_{i + 1}}^{(n)}, t_{i + 1}) - Φ (x_{t_{i}}, t_{i}))

(9)

with $x_{t_{i + 1}}^{(0)} = x_{t_{i}}$ .

Since the corresponding HARP values ϴ = [θ₁, θ₂]^{are just the principal value of the tagging phase Φ, (}9) cannot be used directly. If one assumes, however, that any material point moves less than half of the tag separation from one time frame to the next in both tag orientations—i.e., |ϕ_k(x_{t_i+1}, t_i+1) − ϕ_k(x(t_i), t_i)) < π, k = 1, 2—then it can be shown that

Φ (x_{t_{i + 1}}^{(n)}, t_{i + 1}) - Φ (x (t_{i}), t_{i}) = W (ϴ (x_{t_{i + 1}}^{(n)}, t_{i + 1}) - ϴ (x (t_{i}), t_{i}))

(10)

Moreover, the gradient of Φ can be written as:

∇Φ = ∇^ϴ = ∇^[θ₁,θ₂]^T

(11)

where

\nabla^{*} θ_{k} = {\begin{matrix} \nabla θ_{k}, & ∣ ∣ \nabla θ_{k} ∣ ∣ \leq ∣ ∣ \nabla W (θ_{k} + π) ∣ ∣ \\ \nabla W (θ_{k} + π) ∣ ∣, & otherwise \end{matrix} k = 1, 2

(12)

Equation (9) can then be written as

x⁼x^{− ∇}^ϴ(x^,t_i+1)^W(ϴ(x^,t_i+1) − ϴ(x(t_i), t_i))

(13)

which is computable from the underlying images. HARP tracking is therefore the iteration of Eqn. (13) until ||x⁻x^{|| is below a pre-specified small number, or until a specific number of iterations is reached.}

2.1 MRI tagging and HARP method

A(x, t) = M₀(1 − e^{+ cos(}ϕ(x, t))e^),

(1)

ϕ(x, t₀) = kx ⋅ n + ϕ₀,

(2)

where k is the tagging frequency, n is the normal unit vector that represents the tagging direction and ϕ₀ is the phase o set. When using a [+90°, −90°] RF pulse pair, the image is:

B(x, t) = M₀(1 − e^{− cos(}ϕ(x, t))e^).

(3)

In order to improve the tag contrast and tag persistency, these two images can be combined together to generate the CSPAMM (¹¹) or MICSR (¹²) images:

\begin{matrix} I_{CSPAMM} = A - B = 2 M_{0} e^{- t ∕ T_{1}}) \cos (ϕ (x, t)), \\ I_{MICSR} = A^{2} - B^{2} = 4 M_{0}^{2} (1 - e^{- t ∕ T_{1}}) e^{- t ∕ T_{1}} \cos (ϕ (x, t)) . \end{matrix}

(4)

I(x, t) = D(x, t)e^,

(5)

θ(x, t) = W(ϕ(x, t))

(6)

where W (·) is the wrapping function defined as:

W(ϕ) = mod(ϕ + π, 2π) − π

(7)

Figure 1

(a) Part of a SPAMM tagged MR image, (b) the magnitude of its Fourier transform, and (c) the HARP image after applying the bandpass filter as in (b).

2.2 HARP motion tracking

Φ(x_{t_i+1}, t_i+1) = Φ(x_{t_i}, t_i)

(8)

This is a multidimensional root finding problem, and can be solved iteratively using the Newton-Raphson technique:

x_{t_{i + 1}}^{(n + 1)} = x_{t_{i + 1}}^{(n)} - \nabla Φ {(x_{t_{i + 1}}^{(n)}, t_{i + 1})}^{- 1} (Φ (x_{t_{i + 1}}^{(n)}, t_{i + 1}) - Φ (x_{t_{i}}, t_{i}))

(9)

with $x_{t_{i + 1}}^{(0)} = x_{t_{i}}$ .

Φ (x_{t_{i + 1}}^{(n)}, t_{i + 1}) - Φ (x (t_{i}), t_{i}) = W (ϴ (x_{t_{i + 1}}^{(n)}, t_{i + 1}) - ϴ (x (t_{i}), t_{i}))

(10)

Moreover, the gradient of Φ can be written as:

∇Φ = ∇^ϴ = ∇^[θ₁,θ₂]^T

(11)

where

\nabla^{*} θ_{k} = {\begin{matrix} \nabla θ_{k}, & ∣ ∣ \nabla θ_{k} ∣ ∣ \leq ∣ ∣ \nabla W (θ_{k} + π) ∣ ∣ \\ \nabla W (θ_{k} + π) ∣ ∣, & otherwise \end{matrix} k = 1, 2

(12)

Equation (9) can then be written as

x⁼x^{− ∇}^ϴ(x^,t_i+1)^W(ϴ(x^,t_i+1) − ϴ(x(t_i), t_i))

(13)

3. METHOD

Standard HARP tracking only considers the temporal dependency of a point trajectory—i.e., that a point should not move much from one time frame to the next. The idea of using spatial dependency of motion was first proposed by Osman et al.³ and was called HARP refinement. HARP refinement was first proposed to track points on concentric circles that are placed manually within the left ventricular (LV) myocardium, and was recently extended by region growing HARP refinement method (RG-HR)⁸ to track any tissue point. RG-HR tracks tissue points in an order based solely on the local spatial smoothness of the phase images. After a seed is manually picked, this method automatically tracks all pixels in an image by using a region growing algorithm. Unfortunately, the motion estimate computed at any given point may depend on the specific path over which the motion is estimated from the seed to the given point. It is possible, for example, for a point within the free wall of the left ventricle to be assigned a displacement or phase value based on a seed in the septum and a path that travels through the liver for some distance rather than entirely through the myocardium. Following such paths may yield incorrect tracking. Therefore, region growing according to local HARP phase smoothness is not enough to promise an optimal path. A very similar strategy was employed to unwrap phase images for motion and strain estimation using the DENSE imaging framework,¹⁰ and will therefore suffer from the same problem. To overcome this problem, we developed a new HARP refinement method by formulating the problem as a single source shortest path problem; we call this the shortest path HARP refinement (SP-HR) method. Rather than making local decisions on how the region should be grown relative to the region’s current boundary, SP-HR determines the optimal path from the seed to each point in the image.

In SP-HR, the image is represented as an undirected graph G = (V, E), where V is the set of vertices, and E is the set of edges, as shown in Fig. 2(a). Consider tracking points within an image at time t₁ to a later time frame t₂. Each point (pixel) x(t₁) in the image at time t₁ is represented as a vertex v in the graph, and each edge e_ij = 〈v_i, v_j〉 in E corresponds to a neighboring vertex pair v_i and v_j. Each edge has an edge cost C_E(e_ij), which is non-negative and measures the dissimilarity between the two end vertices. As well, each vertex is associated with a vertex cost C_V (v). Both types of costs are defined below. With this framework, the HARP tracking refinement problem can be formulated as a single source shortest path problem, which lends itself to an optimal solution that can be computed very efficiently.

An external file that holds a picture, illustration, etc.
Object name is nihms-158349-f0002.jpg

Figure 2

The graph representation of the refinement problem and numerical example.

3.1 Cost functions

We define the edge cost using both the harmonic magnitudes and the implied motions at the two end vertices. Since the motion within a single tissue is smooth, then the displacements of neighboring tissue points should be similar. Accordingly, if x_i(t₁) and x_j(t₁) are adjacent tissue points at time t₁, and at time t₂ they move to x_i(t₂) and x_j(t₂), respectively, then the difference between their displacements, |(x_i(t₂) − x_i(t₁)) − (x_j(t₂) − x_j(t₁))| should be small. If x_i(t₁) is correctly tracked to x_i(t₂), then

x_{j}^{'} (t_{2}) = x_{j} (t_{1}) + x_{i} (t_{2}) - x_{i} (t_{1})

(14)

should be approximately equal to x_j(t₂). In fact, if $x_{j}^{'} (t_{2}) \approx x_{j} (t_{2})$ then both of their harmonic phase values should be approximately equal and the following phase similarity function

Δ (x_{i}, x_{j}) = \sum_{k = 1}^{2} ∣ W (ϕ_{k} (x_{j} (t_{1}), t_{1}) - ϕ_{k} (x_{j}^{'} (t_{2}), t_{2})) ∣

(15)

should be small. This function defines that part of the edge cost that depends on the smoothness of harmonic phases.

Pixels that fall outside of tissue—e.g., air or bone—have very weak MR signal and therefore very unreliable harmonic phases. We use the harmonic magnitude images to determine whether adjacent points (in both time frames) are likely to be within tissue or not by defining the weighting functions

w_{1} (x_{i}, x_{j}) = (\overset{‒}{D} (x_{i} (t_{1}), t_{1}) + \overset{‒}{D} (x_{j} (t_{1}), t_{1})) ∕ 2 and

(16)

w_{2} (x_{i}, x_{j}) = (\overset{‒}{D} (x_{i} (t_{2}), t_{2}) + \overset{‒}{D} (x_{j}^{'} (t_{2}), t_{2})) ∕ 2,

(17)

where D‾(x, t) = D₁(x, t) + D₂(x, t) is normalized to the interval [0, 1].

When Δ(x_i, x_j) is small and both w₁(x_i, x_j) and w₂(x_i, x_j) are large, then the edge cost should be small. Accordingly, we define the edge cost function associated with the two neighboring points x_i and x_j as

C_{E} (e_{i j}) = \frac{Δ (x_{i}, x_{j})}{w_{1} (x_{i}, x_{j}) w_{2} (x_{i}, x_{j})} .

(18)

The edge cost function penalizes both a lack of harmonic phase smoothness and edges that cross between tissue and background.

Given a seed vertex v₀, we define the vertex cost function C_V (v_i) at any other vertex v_i as the accumulated edge cost along the shortest path from v₀ to v_i. For any path p = 〈v₀, v₁, v₂, ..., v_i〉 in G from v₀ to v_i, its accumulated cost is $C (p) = \sum_{k = 1}^{i} C_{E} (〈 v_{k - 1}, v_{k} 〉)$ . Therefore the vertex cost function at v_i is

C_V(v_i) = min{C(p):p is any path from v₀tov_i}

(19)

where C_V (v₀) is set to zero. The vertex cost function helps to establish the best path starting from the seed by which to define the initial estimated displacement at any pixel within the image. The actual estimated displacement at any pixel is found by iterating (13) starting from the initial estimate determined by the best path.

3.2 Motion tracking via shortest path following

In SP-HR, the shortest path from the single manually specified source (seed) to every point is resolved using Dijkstra’s algorithm.¹³ The overall algorithm is a region growing algorithm in the sense that boundary pixels are successively added to the growing list of points comprising a region. But in addition to keeping track of the region itself, for points on the region boundary the shortest path costs and the nearest neighbors (toward the seed)—we call them the predecessors—along the shortest paths are also computed and stored. When a new point is added on the basis that its shortest path cost is the next smallest, it is immediately tracked using traditional HARP tracking initialized using the displacement of the point’s predecessor on the shortest path.

To carry out Dijkstra’s algorithm, the vertices are classified into three disjoint sets: the boundary vertex set V_b, the tracked vertex set V_t and the unvisited vertex set V_u. The boundary vertices are maintained in a linked list structure that is sequentially sorted based on their vertex costs—i.e., the first vertex in the list has the smallest vertex cost. We denote N(v) to be the predecessor of v on its shortest path, and u(v) = x(t₂) − x(t₁) to be the displacement of the point x associated with v. Given these notations, the SP-HR algorithm is summarized as follows.

Shortest Path Harp Refinement (SP-HR)

Pick a seed vertex v₀. Set N(v₀) = v₀, and u(v₀) = 0. Initialize V_t = ∅, V_b = {v₀}, and V_u = V\{v₀}. Set C_V (v₀) = 0, and C_V (v) = ∞, for ∀_v ∈ V_u.
Remove the first vertex v_k in V_b. Set V_t = V_t ∪ {v_k}. Track the point x_k associated with v_k using the traditional HARP method, but starting from the initialization of $x_{k}^{'} = x_{k} (t_{1}) + u (N (v_{k}))$ instead of x_k(t₁).
For each v_i such that 〈v_k, v_i〉 ∈ E and v_i ∉ V_t, compute the edge cost value C_E(e_ki), and update its vertex cost value C_V (v_i):
1. Compute the new cost $C_{V}^{'} (v_{i}) = C_{V} (v_{k}) + C_{E} (e_{k i})$ .
2. If it is larger than the existing cost, i.e., $C_{V}^{'} (v_{i}) > C_{V} (v_{i})$ , do nothing.
3. Otherwise, set $C_{V} (v_{i}) = C_{V}^{'} (v_{i})$ , and N(v_i) = V_k. If v_i ∈ V_b, let V_b = V_b\{v_i}, then re-insert v_i into the sorted list V_b based on the updated C_V (v_i); if v_i ∈ V_u, let V_u = V_u\{v_i}, and insert v_i into V_b based on C_V (v_i).
If V_b ≠ ∅, go to Step 2. Otherwise, exit the algorithm.

A numerical example is given in Fig. 2(a) and 2(b). In this example, v₃₃ is selected as the seed. Fig. 2(a) shows that at some iteration, four vertices (v₂₂, v₃₂, v₄₃, and v₃₄) besides the seed have been tracked and marked in yellow. The green vertices are the current boundary vertices, and the white are the current unvisited vertices. Among the boundary vertices, v₂₄ has the lowest cost (1.6+0.5=2.1) via path 〈v₃₃, v₃₄, v₂₄〉. Therefore it is tracked, and the edge costs on 〈v₂₄, v₂₃〉 and 〈v₂₄, v₁₄〉 are computed (marked in red). The vertex costs of v₁₄ and v₂₃ are updated accordingly, and v₂₃ becomes the predecessor of both vertices. In addition, v₁₄ is changed to a boundary vertex. The edge and vertex costs of the graph at the end of this iteration is shown in Fig. 2(b).

Previous region growing motion estimation methods used an edge criterion alone to decide what point to track next. The inherent problem with this approach is that determination of what point to track next depends on what points are currently on the boundary, and these points are (potentially) far from the seed and have no tight relationships with the seed. This permits potentially unnatural effective paths to determine how a given point in the image is tracked. In contrast, the SP-HR algorithm ties the tracking of every point throughout the image directly to the seed point. Since the seed’s displacement is certified by the user to be correctly tracked, this means that a given point is more likely to be tracked correctly since it is tied by reliable estimates throughout its own optimal path—determined by motion smoothness and reluctance to cross boundaries—back to the seed. Therefore, if a point belongs to the same segmented region as the seed, it is far less likely that the point’s displacement will be determined by a path that leaves the segmented region and then returns. This property makes it far less likely that gross tracking errors will occur within the region of interest defined by the seed—e.g., the myocardium or the tongue. We will also see that this approach permits us to correctly track points that are very near to a tissue boundary. This happens because these points will almost always be tied back to the seed through the region of interest defined by the seed; continuity of motion within the region will therefore play a leading role rather over the simpler criterion of similarity of harmonic phase.

3.1 Cost functions

x_{j}^{'} (t_{2}) = x_{j} (t_{1}) + x_{i} (t_{2}) - x_{i} (t_{1})

(14)

Δ (x_{i}, x_{j}) = \sum_{k = 1}^{2} ∣ W (ϕ_{k} (x_{j} (t_{1}), t_{1}) - ϕ_{k} (x_{j}^{'} (t_{2}), t_{2})) ∣

(15)

should be small. This function defines that part of the edge cost that depends on the smoothness of harmonic phases.

w_{1} (x_{i}, x_{j}) = (\overset{‒}{D} (x_{i} (t_{1}), t_{1}) + \overset{‒}{D} (x_{j} (t_{1}), t_{1})) ∕ 2 and

(16)

w_{2} (x_{i}, x_{j}) = (\overset{‒}{D} (x_{i} (t_{2}), t_{2}) + \overset{‒}{D} (x_{j}^{'} (t_{2}), t_{2})) ∕ 2,

(17)

where D‾(x, t) = D₁(x, t) + D₂(x, t) is normalized to the interval [0, 1].

C_{E} (e_{i j}) = \frac{Δ (x_{i}, x_{j})}{w_{1} (x_{i}, x_{j}) w_{2} (x_{i}, x_{j})} .

(18)

The edge cost function penalizes both a lack of harmonic phase smoothness and edges that cross between tissue and background.

C_V(v_i) = min{C(p):p is any path from v₀tov_i}

(19)

3.2 Motion tracking via shortest path following

Shortest Path Harp Refinement (SP-HR)

Pick a seed vertex v₀. Set N(v₀) = v₀, and u(v₀) = 0. Initialize V_t = ∅, V_b = {v₀}, and V_u = V\{v₀}. Set C_V (v₀) = 0, and C_V (v) = ∞, for ∀_v ∈ V_u.
Remove the first vertex v_k in V_b. Set V_t = V_t ∪ {v_k}. Track the point x_k associated with v_k using the traditional HARP method, but starting from the initialization of $x_{k}^{'} = x_{k} (t_{1}) + u (N (v_{k}))$ instead of x_k(t₁).
For each v_i such that 〈v_k, v_i〉 ∈ E and v_i ∉ V_t, compute the edge cost value C_E(e_ki), and update its vertex cost value C_V (v_i):
1. Compute the new cost $C_{V}^{'} (v_{i}) = C_{V} (v_{k}) + C_{E} (e_{k i})$ .
2. If it is larger than the existing cost, i.e., $C_{V}^{'} (v_{i}) > C_{V} (v_{i})$ , do nothing.
3. Otherwise, set $C_{V} (v_{i}) = C_{V}^{'} (v_{i})$ , and N(v_i) = V_k. If v_i ∈ V_b, let V_b = V_b\{v_i}, then re-insert v_i into the sorted list V_b based on the updated C_V (v_i); if v_i ∈ V_u, let V_u = V_u\{v_i}, and insert v_i into V_b based on C_V (v_i).
If V_b ≠ ∅, go to Step 2. Otherwise, exit the algorithm.

Shortest Path Harp Refinement (SP-HR)

Pick a seed vertex v₀. Set N(v₀) = v₀, and u(v₀) = 0. Initialize V_t = ∅, V_b = {v₀}, and V_u = V\{v₀}. Set C_V (v₀) = 0, and C_V (v) = ∞, for ∀_v ∈ V_u.
Remove the first vertex v_k in V_b. Set V_t = V_t ∪ {v_k}. Track the point x_k associated with v_k using the traditional HARP method, but starting from the initialization of $x_{k}^{'} = x_{k} (t_{1}) + u (N (v_{k}))$ instead of x_k(t₁).
For each v_i such that 〈v_k, v_i〉 ∈ E and v_i ∉ V_t, compute the edge cost value C_E(e_ki), and update its vertex cost value C_V (v_i):
1. Compute the new cost $C_{V}^{'} (v_{i}) = C_{V} (v_{k}) + C_{E} (e_{k i})$ .
2. If it is larger than the existing cost, i.e., $C_{V}^{'} (v_{i}) > C_{V} (v_{i})$ , do nothing.
3. Otherwise, set $C_{V} (v_{i}) = C_{V}^{'} (v_{i})$ , and N(v_i) = V_k. If v_i ∈ V_b, let V_b = V_b\{v_i}, then re-insert v_i into the sorted list V_b based on the updated C_V (v_i); if v_i ∈ V_u, let V_u = V_u\{v_i}, and insert v_i into V_b based on C_V (v_i).
If V_b ≠ ∅, go to Step 2. Otherwise, exit the algorithm.

4. EXPERIMENT RESULTS

We applied SP-HR to track the motion of the left ventricle (LV) of the human heart. Experiments were carried out on 13 short axis cine image sequences acquired over time from four different subjects covering different parts of the heart. In each case the tag period was 12 mm and the slice thickness was 8 mm. The numbers of time frames varied from 15 to 20, and the temporal resolutions varied from 20 ms to 43 ms. The image sizes were all 256×256, with FOVs either 300 mm or 320 mm.

To quantitatively evaluate our method, we asked 20 volunteers to delineate the LVs from all the image sequences. We picked one time frame from each of the 13 image sequences. After some training (on other images), each volunteer manually drew the epicardial and endocardial contours of the left ventricles together with two insertion points of the right ventricle on each of the 13 selected images. For each image sequence, we determined the LV myocardium region by averaging all 20 rater’s epicardial and endocardial contours, and including all pixels within these average boundaries. A seed point was manually specified at the first time frameand tracked through all other time frames. The displacement fields from the first frame to all later time frames were then computed using both RG-HR and SP-HR, and the results were compared. Tracking was considered to be successful if the computed motion field inside LV has no visible abrupt jump. Over all the image sequences and all time frames, the success rate of RG-HR was 93.8%, while SP-HR was 99.5%. Fig. 3 shows an example in which RG-HR fails while SP-HR succeeds. Fig. 3(b) and (c) show the magnitude of the displacement field from the the image displayed in Fig. 3(a) to a later time frame using RG-HR and SP-HR methods, respectively. We observe that RG-HR result contains abrupt jumps of displacement inside the myocardium region, which indicates errors in tracking. In contrast, the SP-HR method result is smooth throughout the whole tissue.

An external file that holds a picture, illustration, etc.
Object name is nihms-158349-f0003.jpg

Figure 3

Comparison of RG-HR and the SP-HR method. (a) A grid tag image generated by multiplying the orthogonally tagged image pair; (b) magnitude of the displacement field computed using RG-HR; and (c) magnitude of the displacement field computed using SP-HR. Units are in mm.

Fig. 4 shows the intermediate results of the two methods on the same example as in Fig. 3. The evolution of the tracked regions is depicted as iteration increases from left to right. It is observed that in RG-HR, the liver is tracked before the LV is completely tracked while in SP-HR, points in the liver are tracked only after nearly all LV points are tracked. Because of the globally defined vertex cost function in SP-HR, the tracked region grows inside the LV with a similar speed on all sides of the seed, which is not true in RG-HR. Thus, points inside LV are reached via an optimal path in SP-HR and are more likely to be correctly tracked.

An external file that holds a picture, illustration, etc.
Object name is nihms-158349-f0004.jpg

Figure 4

Illustration of the tracking order of (top row) RG-HR and (bottom row) SP-HR. The red cross marks the seed.

To evaluate the effectiveness of our method when there is large motion and low temporal resolution, we used subsets of the image sequences obtained by dropping intermediate images. In particular, for a time step n, the points were tracked in the image series at time frames 1, n+1, 2n+1, ... only. The images at the last time frames were always included so the final computed displacements could be directly compared. By using time step otherthan n = 1, we mimic the situation of large motion and/or poor temporal resolution. In our experiments we computed the displacement fields from the first time frame to the last time frame using both SP-HR method and traditional HARP tracking with time step 1, 2 and 3 on all the 13 image sequences. From the tracking results, we computed the ratio of correctly tracked points inside the LV myocardium.

Fig. 5 (left to right) shows the magnitude of the displacement fields on one of the 13 test image sequences computed with SP-HR and traditional HARP tracking with time steps 1, 2 and 3. Fig. 6 compares the correctly tracked regions on one of the 13 test image sequences using SP-HR with traditional HARP with different time steps. We observe that with time step more than 1 traditional HARP mistracks a large portion of the LV region and performs much worse than SP-HR, and for time step 1, traditional HARP performs slightly worse than SP-HR (small number of yellow pixels). Fig. 7 illustrates the ratios of correctly tracked points in the LV myocardium of 13 image sequences when using SP-HR and traditional HARP tracking with time steps 1, 2, and 3. It can be seen from Fig. 7 that SP-HR worked better than traditional HARP in all image sequences. On average, 98.4% points were correctly tracked using SP-HR, 93.1% in traditional HARP tracking with time step 1, 79.3% with time step 2, and 58.7% with time step 3.

An external file that holds a picture, illustration, etc.
Object name is nihms-158349-f0005.jpg

Figure 5

Example of the magnitude of displacement field from the first to the last time frame computed using different methods. (a) SP-HR; (b-d) traditional HARP tracking with time steps (b) 1, (c) 2 and (d) 3.

An external file that holds a picture, illustration, etc.
Object name is nihms-158349-f0006.jpg

Figure 6

Comparison of the correctly tracked regions using SP-HR and traditional HARP with time steps (a) 1, (b) 2 and (c) 3. Blue: correctly tracked regions using traditional HARP with different time steps; yellow+blue: correctly tracked regions using SP-HR; Red: average LV myocardium boundaries.

An external file that holds a picture, illustration, etc.
Object name is nihms-158349-f0007.jpg

Figure 7

The percentage of correctly tracked points when using SP-HR and traditional HARP with different time steps on the 13 test image sequences, which is shown from left to right.

We implemented the SP-HR method in C and compiled and ran it in Matlab (Mathworks, Natick MA). It works as fast as the RG-HR method, because the computation caused by the shortest path calculation is negligible in comparison to the remaining computation. On a computer with Intel Core Duo 1.83 GHz processor and 1.0 GB of RAM, our implementation took less than 0.2 second to track an 128 by 128 image for one time frame.

5. CONCLUSION

We presented a new refinement method for HARP tracking by formulating a single source shortest path problem and solving it using Dijkstra’s algorithm. Experimental results show that this method can reliably track every point inside the tissue. This method is also computationally fast and makes it feasible to reliably compute other useful quantities in functional analysis of motion using MRI.

Image Analysis and Communication Laboratory, Johns Hopkins University, Baltimore, MD USA

Send correspondence to Jerry L. Prince, ude.uhj@ecnirp; Phone: +1(410)516-5192; Fax: +1(410)516-5566.

Abstract

Harmonic phase (HARP) motion analysis is widely used in the analysis of tagged magnetic resonance images of the heart. HARP motion tracking can yield gross errors, however, when there is a large amount of motion between successive time frames. Methods that use spatial continuity of motion—so-called refinement methods—have previously been reported to reduce these errors. This paper describes a new refinement method based on shortest-path computations. The method uses a graph representation of the image and seeks an optimal tracking order from a specified seed to each point in the image by solving a single source shortest path problem. This minimizes the potential for path dependent solutions which are found in other refinement methods. Experiments on cardiac motion tracking shows that the proposed method can track the whole tissue more robustly and is also computationally efficient.

Keywords: Motion tracking, HARP, refinement, single source shortest path problem, Dijkstra’s algorithm

Abstract

REFERENCES

References

1. Axel L, Dougherty LMR imaging of motion with spatial modulation of magnetization. Radiology. 1989;171:841–845.[PubMed][Google Scholar]
2. Axel L, Doughety LHeart wall motion: improved method of spatial modulation of magnetization for MR imaging. Radiology. 1989;172:349–350.[PubMed][Google Scholar]
3. Osman NF, Kerwin WS, McVeigh ER, Prince JLCardiac motion tracking using CINE harmonic phase (HARP) magnetic resonance imaging. Mag. Res. Med. 1999;42:1048–1060.[Google Scholar]
4. Osman NF, Prince JLVisualizing myocardial function using HARP MRI. Phys. Med. Biol. 2000;45:1665–1682.[PubMed][Google Scholar]
5. Osman NF, McVeigh ER, Prince JLImaging heart motion using harmonic phase MRI. IEEE Trans. Med. Imag. 2000;19:186–202.[PubMed][Google Scholar]
6. Khalifa AM, Youssef AM, Osman NF. Improved harmonic phase (HARP) method for motion tracking a tagged cardiac MR images; 2005 IEEE EMBS Ann. Conf.; 2005.pp. 4298–4301. [[PubMed]
7. Tecelao SR, Zwanenburg JJ, Kuijer JP, Marcus JJExtended harmonic phase tracking of myocardial motion: Improved coverage of myocardium and its effect on strain results. Journal of Magnetic Resonance Imaging. 2006;23(5):682–690.[PubMed][Google Scholar]
8. Liu X, Murano E, Stone M, Prince JL. HARP tracking refinement using seeded region growing; Proc. of IEEE Int. Sym. Med. Imag.; 2007.pp. 372–375. [PubMed]
9. Aletras AH, Ding S, Balaban RS, Wen HDENSE: Displacement encoding with stimulated echoes in cardiac functional MRI. J. Mag. Res. 1999;137:247–252.[Google Scholar]
10. Spottiswoode BS, Zhong X, Hess AT, Kramer CM, Meintjes EM, Mayosi BM, Epstein FHTracking myocardial motion from cine DENSE images using spatiotemporal phase unwrapping and temporal fitting. IEEE Tran. Med. Imag. 2007;26(1):15–30.[PubMed][Google Scholar]
11. Fischer SE, McKinnon GC, Maier SE, Boesiger PImproved myocardial tagging contrast. Magnetic Resonance in Medicine. 1993;30(2):191–200.[PubMed][Google Scholar]
12. NessAiver M, Prince JLMagnitude image CSPAMM reconstruction (MICSR) Mag. Res. Med. 2003;50(2):331–342.[PubMed][Google Scholar]
13. Dijkstra EWA note on two problems in connexion with graphs. Numerische Machematik. 1959;1:269–271.[PubMed][Google Scholar]