HARP (algorithm)

HARP processing

Demonstration of the HARP processing of a tagged cardiac MRI slice. (a) An MR image with vertical SPAMM tags. (b) shows the magnitude of its Fourier transform. By extracting the spectral peak inside the circle, a complex image is produced with a magnitude (c) and a phase (d).

A tagged MRI showing motion of a human heart is shown in the image (a). The effect of tagging can be described as a multiplication of the underlying image by a sinusoid tag pattern having a certain fundamental frequency, causing an amplitude modulation of the underlying image and replicating its Fourier transform into the pattern shown in (b).

HARP processing uses a bandpass filter to isolate one of the spectral peaks. For example, the circle drawn in (b) is the -3 dB isocontour of the bandpass filter used to process this data. Selection of the filters for optimal performance is discussed in this paper.^[14] The inverse Fourier transform of the filtered image yields a complex harmonic image ${\displaystyle I_{k}(\mathbf {y} ,t)}$ at image coordinates ${\displaystyle \mathbf {y} =[y_{1},y_{2}]^{T}}$ and time ${\displaystyle t}$:

${\displaystyle I_{k}(\mathbf {y} ,t)=D_{k}(\mathbf {y} ,t)e^{j\phi _{k}(\mathbf {y} ,t)}}$

where ${\displaystyle D_{k}}$ is called the harmonic magnitude image and ${\displaystyle \phi _{k}}$ is called the harmonic phase image. The harmonic magnitude image in (c) extracted from a using the filter in (b) shows the geometry of the heart. And the harmonic phase image in (d) contains the motion of the myocardium in horizontal direction. In practice, tagged images from two directions (both horizontal and vertical, i.e., ${\displaystyle k}$ is 1 and 2) are processed to provide a 2D motion map in the image plane. Notice that the harmonic phase images are computed by taking the inverse tangent of the imaginary part divided by the real part of ${\displaystyle I_{k}(\mathbf {y} ,t)}$, such that the range of this computation is only in ${\displaystyle [-\pi ,+\pi )}$. In other words, d is only the wrapped value of the actual phase. We denote this principle value by ${\displaystyle a_{k}(\mathbf {y} ,t)}$; it is mathematically related to the true phase by:

${\displaystyle a_{k}(\mathbf {y} ,t)=mod(\phi _{k}(\mathbf {y} ,t)+\pi ,2\pi )-\pi }$

Either ${\displaystyle \phi _{k}}$ or ${\displaystyle a_{k}}$ might be called a harmonic phase (HARP) image, but only ${\displaystyle a_{k}}$ can be directly calculated and visualized. It is the basis for HARP tracking.

HARP tracking

Result of HARP tracking of a tagged cardiac MRI slice

For a fixed material point with a HARP value, only one of the points sharing the same HARP value in a later time frame is the correct match. If the apparent motion is small from one image to the next, it is likely that the nearest of these points is the correct point. The tracking result is very accurate in this case.

Consider a material point located at ${\displaystyle \mathbf {y} _{m}}$ at time ${\displaystyle t_{m}}$. If ${\displaystyle \mathbf {y} _{m+1}}$ is the apparent position of this point at time ${\displaystyle t_{m+1}}$, we have:

${\displaystyle \phi _{k}(\mathbf {y} _{m+1},t_{m+1})=\phi _{k}(\mathbf {y} _{m},t_{m})}$

The Newton–Raphson interactive method is used to find a solution, which is:

${\displaystyle y^{(n+1)}=y^{(n)}-[\nabla \phi _{k}(\mathbf {y} ^{(n)},t_{m+1})]^{-1}[\phi _{k}(\mathbf {y} ^{(n)},t_{m+1})-\phi _{k}(\mathbf {y} _{m},t_{m})]}$

In practice, since ${\displaystyle \phi _{k}}$ is not available, ${\displaystyle a_{k}}$ is used in its place. This equation can be rewritten after a few derivations considering the "wrapping" relation between ${\displaystyle \phi _{k}}$ and ${\displaystyle a_{k}}$.

The result of HARP tracking of one frame of cardiac MRI is shown in the figure. It is obtained by calculating both motions from horizontal direction and vertical direction, resulting in a 2D vector field showing the motion of every material point on the myocardium at this time frame.

The entire HARP algorithm takes only a few minutes to perform on a normal computer and the motion tracking result is accurate (with a typical error range of ${\displaystyle \pm 1}$ pixel). As a result, it is now widely adopted by the medical image analysis community as a standard processing technique for tagged MRI.