ABSTRACT

Phase Unwrapping (PU) is an ill-posed problem in Magnetic Resonance Elastography (MRE). The phase information is not usable until the phases are retrieved by using PU algorithms. In this present study, we attempt to determine the ideal PU method for MRE using both phantom and volunteer psoas major (PM) muscle images. All the MRE experiments were carried out in Philips MRI (Achieva 3.0 T, Best, The Netherlands). A multi-echo gradient-echo MRE pulse sequence was employed and the four PU methods were considered based on their easy user platform. They are namely, Minimum Discontinuity (MD), Laplacian-Based Estimate (LBE), Region Growing (RG) and Dilate-Erode (DE) Propagate. Phantom images were successfully unwrapped by all four methods, whereas MD and LBE could only unwrap PM muscle images properly. RG and DE failed to unwrap the PM muscle images.

Keywords:

Phase Unwrapping, Magnetic Resonance Elastography, Minimum Discontinuity, Laplacian-Based Estimate, Region Growing, Dilate-Erode Propagate

1. Introduction

Magnetic Resonance Elastography (MRE) [1] is a dedicated imaging modality that uses phase contrast imaging to determine the viscoelastic properties of soft tissues. The stiffness measurement of the tissue is a fundamental indicator for disease diagnosis [2]. In the past, manual palpation had been considered a gold standard technique to obtain such information. With an inception of MRE, palpation has become possible virtually and is even extendable to internal organs and tissues [1]. MRE primarily involves wave image processing and phase unwrapping (PU) is an initial process in it [3]. The phase carries essential information regarding the motion of a test subject being imaged [4]. However, the phase is wrapped within the range [−π, π] in the image and the wrapped phase is not usable until the 2π phase discontinuities are not fixed by using phase unwrapping methods [5]. It requires dedicated methods to extract this delicate information. The true phase retrieving process is not straightforward and this problem has been stated radically intractable [6].

There are plenty of literatures published in phase unwrapping methods. Yet, the meta-analysis is scarce. Most of the researchers and laboratories still rely on algorithms outlined by Ghiglia and Pritt 20 years ago [7]. During this time, PU has gained huge attention due to its wide applications in Synthetic Aperture Radar (SAR), Magnetic Resonance Imaging (MRI), and optical imaging [8]. From the time of advent of MRE, phase wraps have become a troublesome. However, many researchers have attempted to solve this problem. Still, the PU remained an ill-posed problem. In the present context where MRE has been proved as a substantial imaging method for staging of liver fibrosis [9] , diagnosis of certain brain pathologies [10] , detection of breast lesion [11] and many others, it would be better to establish a robust and computationally efficient method of PU for all types of organ imaging. The main purpose of this study was to determine the best PU method for MRE imaging using phantom and volunteer psoas major (PM) muscle.

2. Theory

2.1. Phase Unwrapping (PU)

PU is a method of determining the total accumulated phase [7]. Mathematically, PU can be operated by obtaining an unwrapped estimate $\varphi \left(t\right)$ from the wrapped function $\psi \left(t\right)$ . That is,

$\psi \left(t\right)=\varphi \left(t\right)+2\text{π}k\left(t\right)$ (1)

$\varphi \left(t\right)$ is some unknown phase function before unwrapping and $k\left(t\right)$ is an integer function that forces $-\text{π}<\psi \le \text{π}$ .

The phase of a signal measurement can vary from 0 to 2π or −π to π. “Phase wrap” means that the signal keeps on rotating for more than one complete rotation. Figure 1 illustrates the wrapped and unwrapped phases in a one dimensional signal. Whenever the phase value is greater than π, the value is wrapped back into the negative axis and vice-versa. The values outside the range of ? π to π are wrapped.

2.2. Types of PU Methods

Four types of phase unwrapping methods are used in this study. They are listed as follows: Minimum Discontinuity (MD), Laplacian-Based Estimate (LBE), Region Growing (RG) and Dilate Erode Propagate (DE).

2.2.1. Minimum Discontinuity (MD)

MD, coined by Flynn in 1997 is popularly known as Flynn’s algorithm [12] and is incorporated in MREWave [13] , a free wave inversion software for MRE with phase unwrapping and image filtering steps which can be freely available from [https://www.mayo.edu/research/labs/advanced-medical-imaging-technology/focus-areas/mre-wave].

The basis mechanism of MD is that it adds an integer multiple of 2π to each phase wrapped pixel to minimize the phase difference between adjacent pixels less than π. A discontinuity is a term used when the phase difference between a pair of adjacent pixels exceed π. The jump count of the pixels to be the nearest integer is a integer multiple of 2π of the phase difference. Thus, we can write,

horizontal jump count, $v_{m,n}=\frac{Int\left({\varphi }_{m,n}-{\varphi }_{m+1,n}\right)}{\text{2π}}$ (2)

vertical jump count, $z_{m,n}=\frac{Int\left({\varphi }_{m,n}-{\varphi }_{m,n-1}\right)}{\text{2π}}$ (3)

where, Int(.) denotes the rounding off to the nearest integer. A nonzero jump count indicates that the phase difference between adjacent pixels exceeds π.

A grid plane is introduced in the image that has nodes, where four voxels touch, and is offset by 1/2 with respect to image voxel in both directions. An edge is occurred whenever the jump count between adjacent voxels is non-zero. The edge can have four possible directions (up, down, right and left). Thus, a closed path is searched by following the edges. A closed path is like a closed fringe line. The strategy is to add an appropriate multiple of 2π to all the voxels. The process of adding 2π is iterative and is looped within the total image. The search for discontinuity usually terminates at the boundary of the image or is balanced by an analogous pixel within the image [7] [12].

2.2.2. Laplacian-Based Estimate (LBE)

LBE exploits the fact that $\exp \left(i{\phi }^{\text{wrap}}\right)=\exp \left(i{\phi }^{\text{unwrap}}\right)$ [6]. That is,

${\nabla }^2\varphi =\cos {\varphi }_w{\nabla }^2\left(\sin {\varphi }_w\right)-\sin {\varphi }_w{\nabla }^2\left(\cos {\varphi }_w\right)$ (4)

With ∇² the Laplacian operator. This allows for an estimate of the true phase ${\varphi }^{\prime }$ , to be derived by application of an inverse Laplacian operator ∇⁻² to the right-hand side, Fourier-domain forward and inverse Laplace operators enable a direct solution for ${\varphi }^{\prime }$ in 4D as:

$\begin{array}{l}{\varphi }^{\prime }\left(x,y,z,t\right)\\ =FC{T}^{-1}\left(\frac{FCT\left(\cos {\varphi }_w\right)\left(FC{T}^{-1}\left[\left(p^2+q^2+r^2+s^2\right)FCT\left(\sin {\varphi }_w\right)\right]\right)}{p^2+q^2+r^2+s^2}\right)\\ -FC{T}^{-1}\left(\frac{FCT\left(\sin {\varphi }_w\right)\left(FC{T}^{-1}\left[\left(p^2+q^2+r^2+s^2\right)FCT\left(\cos {\varphi }_w\right)\right]\right)}{p^2+q^2+r^2+s^2}\right)\end{array}$ (5)

with FCT and FCT⁻¹ the forward and inverse fast cosine transforms (x, y, z, t) the time-domain coordinates and (p, q, r, s) the Fourier domain coordinates.

LBE calculates the unwrapped phase as a global approximation of true phase. The exact phase is not determined by Laplacian phase unwrapping. Though LBE calculates the approximate phase values, it preserves much of the phase variation information in the image. LBE does not depend on local phase information and is considered “noise-immune”, in the sense that unwrap failures do not proliferate spatially during the phase unwrapping process. The above equation shows three types of operations: forward and inverse cosine transforms, trigonometric operations, and the masking expression ( $p^2+q^2+r^2+s^2$ ). As the masking terms are in both numerator and denominator, their impact on the spatial domain value range cancels out. Readers who desire to understand comprehensive information about LBE can learn from references [6] [14].

2.2.3. Region Growing (RG)

This algorithm put down by Herráez et al. [5] and Abdul-Rahman et al. [15] depends on phase quality. A simple two-dimensional method of RG is explained below.

Herráez et al. [5] used second difference method for two-dimensional image. The second difference D of an (i, j) pixel can be calculated by the equation:

$D\left(i,j\right)={\left[H^2\left(i,j\right)+V^2\left(i,j\right)+D_1^2\left(i,j\right)+D_2^2\left(i,j\right)\right]}^{\frac{1}{2}}$ (6)

where,

$H\left(i,j\right)=\gamma \left[\phi \left(i-1,j\right)-\phi \left(i,j\right)\right]-\gamma \left[\phi \left(i,j\right)-\phi \left(i+1,j\right)\right]$ (7)

$V\left(i,j\right)=\gamma \left[\phi \left(i,j-1\right)-\phi \left(i,j\right)\right]-\gamma \left[\phi \left(i,j\right)-\phi \left(i,j+1\right)\right]$ (8)

$D_1\left(i,j\right)=\gamma \left[\phi \left(i-1,j-1\right)-\phi \left(i,j\right)\right]-\gamma \left[\phi \left(i,j\right)-\phi \left(i+1,j+1\right)\right]$ (9)

$D_2\left(i,j\right)=\gamma \left[\phi \left(i-1,j+1\right)-\phi \left(i,j\right)\right]-\gamma \left[\phi \left(i,j\right)-\phi \left(i+1,j-1\right)\right]$ (10)

where, $\gamma (\cdot )$ is an unwrapping operation to remove any 2π steps between two consecutive pixels. An edge exists at the intersection of two pixels that are connected in both directions. The unwrapping of edges leads to unwrapping of the image. The highest-quality voxels with the highest-reliability edges are unwrapped first and the lowest-quality voxels with the lowest-reliability edges are unwrapped last to prevent error propagation. It follows a non-continuous path during unwrapping process.

The algorithm fails when the sum of mean phase gradient and the standard deviation of the gaussian noise exceeds π.

$\text{mean phase gradient}+\text{noise}>\text{π}$ (11)

2.2.4. Dilate-Erode Propagate (DE)

DE dilates the most reliable voxels by syncing their neighbours to them and erodes the least reliable voxels by syncing their neighbours to them [6]. Contrary to RG, DE does not form groups. DE requires domain knowledge of choosing a slice by hand to propagate the unwrapping process [16]. In case of light wraps, DE can be a stand-alone unwrapping method. For the propagation, all the wrapped regions on the selected slice must be smaller than the chosen dilate-erode window. A slice can be used as a preconditioner for complex images.

The algorithm is failed when

$\frac{\text{mean phase gradient}+\text{noise}}{{\left(r\right)}^2}>\text{π}$ (12)

The “r” is the radius of the DE window which is self-adjustable. For this study, the optimal value of “r” could be selected, that could efficiently remove the phase wraps.

LBE, RG and DE are integrated in ImageJ plugin named “Phase Tools”, developed by Eric Barnhill [6]. ImageJ [17] is a public image processing and analysis program written in JAVA by National Institute of Health, Bethesda, MD, USA). It can be freely accessible via the website [http://rsbweb.nih.gov/ij/]. Phase Tools can be obtained for free from [https://github.com/ericbarnhill/phaseTools].

3. Methods

3.1. MRE Experiments

All MRE experiments were carried out on a clinical MR scanner (Achieva 3.0 T; Philips Healthcare, Best, The Netherlands). A self-designed waveform generation system (LabView, USB-6221; National Instruments, TX, USA) was used to generate the vibration waveform. This system could generate sinusoidal waveforms with arbitrary frequencies and phases. Power amplifier (XTi 1000; Crown, IN, USA) and Pneumatic pressure generator (Subwoofer TIT320C-4 12’’; Dayton Audio, OH, USA) units were used to supply vibrations to a vibration pad. The vibration pad was designed using a 3D printer. The configuration of MRE devices are shown in Figure 2.

3.2. MEG-Less Multi-Echo MRE Sequence

The multi-echo gradient-echo type MRE pulse sequence [18] was applied. The parameters of the sequence are listed in Table 1. This sequence can produce multiple echoes, however lacks motion-encoding gradient (MEG). The read-out gradient (RO) itself acts as MEG gradient and higher the TE time, greater the MEG-like effect. The MEG-like effect of 2^nd echo is greater than 1^st echo, 3^rd echo > 2^nd echo and so on. The 2^nd echo phase image is selected because transverse relaxation affects higher echo images and 1^st echo has inadequate MEG-like effect.

3.3. Data Acquisition

Phantom and ten healthy young volunteer (10 men, age range: 20 - 25 years) PM muscle images were acquired from the MRI suite. A homogeneous polyacrylamide phantom with concentration (3% w) was taken. The phantom MRE experiments were carried out at 50 Hz vibration using flex coil. The phantom positioning technique is shown in Figure 3(a) & Figure 3(b). The phantom was placed on to the vibration pad and the coils are adjusted in RL direction. Figures 4(a)-(e) represents the PM muscle positioning technique. The silicon sheet was placed between the vibration pad and the skin of the volunteer. The silicon sheet

FA: Flip angle, FOV: Field of view, SENSE: Sensitivity encoding.

was used to prevent air from leaking through the interstices between the vibration pad and the skin on the PM muscle. Psoas Major (PM) muscle [19] is a pair of long fusiform muscle located on either side of lumbar vertebral column. It originates from the transverse processes of T12-L5 and the lateral aspects of intervertebral discs (IVDs) between them and inserts in the lesser trochanter of the femur. Figure 5(a) represents the PM muscle and the slice location at the level of L3-L4 IVD whereas Figure 5(b) is an axial magnitude image that shows the vibration of lumbar vertebrae and the shear wave propagation to the bilateral PM muscle. PM muscle was considered because the magnetic susceptibility of PM muscle is greatly variable than its surrounding tissues. PM muscle MRE experiments were performed at 50 Hz with breath-holding manoeuvre during each MRE acquisition. Vibration must be adequate to reach to PM muscle. Torso Phased-array coil was used for the volunteer studies. An ethical approval of informed consent was obtained from the institutional ethical review board and consent to participate was obtained from volunteers before conducting MRE experiments.

3.4. Phase Unwrapping Methods

Considering as a need to decide the best PU method for MRE, this study has included four PU methods that are easily available and straight forward user platform. They are MD, LBE, RG and DE. The details are described in theory section above. We intend to unwrap the phantom and PM muscle images by these four methods and attempt to determine the best PU method.

4. Results

The results are classified based on phase unwrapping methods that are mentioned

below.

4.1. MD

MD could properly unwrap both the phantom and PM muscle images. There was clear wave visualization in both the phantom and the PM muscle as shown in Figure 6, row (b) and Figure 7, row (b) respectively. The phase disruptions at the location of vortices were anticipated in PM muscle image. The artificial phase offset at the locations of phase disruptions and at the boundary interface of the image was greater than π.

4.2. LBE

LBE could decently unwrap both the phantom and PM muscle images. The wave

propagation was distinct in both the phantom and the PM muscle images as illustrated in Figure 6, row (c) and Figure 7, row(c) respectively. The phase disruptions could not be identified, and the artificial phase offset was not generated.

4.3. RG

This method could successfully unwrap the wrapped phantom data, shown in Figure 6, row (d). Phase residues (discontinuities) in the form of small black dots were noticed in the unwrapped images of the phantom. However, the RG was unable to unwrap the 0, π and 3π/2 phase offset PM muscle images properly. The areas in those three images where RG failed to unwrap are shown by white arrows in Figure 7, row (d).

4.4. DE

RG was used as a preconditioner for DE method. DE was able to unwrap the wrapped phantom data as illustrated in Figure 6, row (e). Phase residues (small black dots) were noticed in the unwrapped images of the phantom. Again, DE could not unwrap 0, π and 3π/2 phase offset PM muscle images appropriately. The areas in those three images where DE failed to unwrap are shown inside white circles in Figure 7, row (e).

5. Discussion

All four methods performed well to retrieve the wrapped phases in the phantom data, whereas MD and LBE could only unwrap the volunteer PM muscle images properly. In other words, RG and DE failed to unwrap the PM muscle images. All the four methods have concomitant groups of strengths and weaknesses.

5.1. Remarks for the Usage of MD, LBE, RG and DE

5.1.1. MD

MD has been used extensively in MRE due to its implementation in MREWave software. However, it possesses some drawbacks. It requires a huge amount of computations and has low efficiency in searching discontinuity area for creating loop within an image. Nevertheless, the operation of adding 2π and creating loop that followed discontinuity curves proved successful during unwrapping. In this study, both phantom and human PM muscle images were successfully unwrapped by MD. In PM muscle, due to proximity of intestine, abdominal aorta (AA) and inferior venacava (IVC), the phase disruptions in Flynn’s algorithm are unavoidable in PMMRE. The phase disruptions are the consequences of an incorrect solution and are generally classified as phase errors. Due to large and rapidly varying phase difference between adjoining pixels, MD could generate the phase errors. That is, MD could create an arbitrary global phase offset of an integer multiple of π. It is expected that the arbitrary phase offset occurred at the regions of phase disruptions. Moreover, those arbitrary phase offset was also noticed at the boundary in both phantom and PM muscle images.

5.1.2. LBE

LBE is a newly developed Laplacian type phase unwrapping method in which fast cosine transform (FCT) is employed to calculate the final approximate solution. It successfully unwrapped both phantom and PM muscle images since it is a robust method. Robust in a way that LBE is entirely a spectral method and is not influenced by local phase disruptions. It is immune to the phase disruptions due to the bowel gas and the blood vessels present at juxtaposition of PM. Another benefit of LBE is that the unwrapped solution is not much greater than [−π, π]. Moreover, LBE possessed lowest variance results.

5.1.3. RG

Though RG could completely unwrap the phantom images, the high image gradient forbade the region-forming basis in PM muscle images. In other words, good phase quality map could not be generated during unwrapping process. Thus, the phase discontinuities could not be resolved, and the areas shown by white arrows in Figure 7, row (d) in the unwrapped images appeared like the corresponding areas in the wrapped images. Again, the phase residues were noticed in the unwrapped images of the phantom due to unresolved phase discontinuities.

5.1.4. DE

In case of light wraps, DE could be used as an independent unwrapping method [6]. However, in this study, RG was required as a preconditioner for DE method. We expect the high image gradient prevented DE to unwrap properly in PM muscle images. Furthermore, the phase residues were noticed in the unwrapped images of the phantom due to the unsolved phase discontinuities.

5.2. Phase Unwrapping in Phantom

Phantom is a homogeneous structure, made up of single material, namely polyacrylamide, because of which, the magnetic environment inside the phantom became constant. The presence of light wraps and the absence of highly varying image gradient made it easy to unwrap. Hence, all four phase unwrapping methods could properly unwrap the phantom images.

5.3. Phase Unwrapping in Psoas Major Muscle

In PM muscle, the viscera and the blood vessels present anterior to it, behave differently in magnetic environment than the back lying PM muscles. The viscera mostly comprising gases and the main blood vessels are the abdominal aorta (AA) and the inferior venacava (IVC) that have drastically varied magnetic susceptibility than the PM muscles. They unambiguously play an unfavourable role. The gases and the blood have propensity to yield higher gradient in the phase image. At the same time, it is inevitable that gut is irregular in shape and size and due to its serpentine and void nature, the PU became difficult in PM muscle. Because of this, PU is considered a very important method in PMMRE. Our primary aim is to remove phase wraps from PM muscle.

The phase disruptions anterior to the PM muscle cannot be a ruling factor to persuade us that the performance of MD is facile. Both MD and LBE performed well in PU of PM muscle. On the contrary, the phase values by LBE are only the approximation values. In another word, LBE can never output exact values, but it did not influence the stiffness measurement from wave-length estimation method, that is local frequency estimate (LFE) as it preserved the phase variation. Beside this, LBE has already been in practice in muscle elastography software for the Siemens IDEA platform [20].

It is observed that the rigidity of PM muscle could be a strong sign of low back pain due to the compression of lumbar IVDs [21]. Recently, PMMRE has captured a huge attention to investigate the non-specific back pain. Numano et al. [22] proposed that PMMRE could prove to be a new bioimaging marker for the back pain. PMMRE, being an emerging technique needs robust PU method. Future MRE research using MD and LBE undertaking wave inversion methods are recommended to decide the better one among these two.

6. Conclusion

Thus, at the meantime, the selection of phase unwrapping method fully depends upon the requirement of the user. If the operator is content with true phase values with phase disruptions anterior to the PM muscle, and the artificial phase offset at the boundary, MD can be used. However, if the user needs robust, with no phase disruptions and do not need exact phase values, LBE is recommended. The structure of PM muscle, its surrounding tissue, the excitation frequency of vibration and the insertion location of vibration pad could possess direct implication to the phase wrap in PMMRE. In general, considering the multi-factorial dependency of PU, it is affected by the composition of the test object, MRE sequence parameters, positioning technique, wave image processing parameters and the choice of best PU method.

Acknowledgements

The authors are thankful to Tokyo Metropolitan Government for research grant.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

Maharjan, S., Numano, T., Habe, T., Ito, D., Ueki, T., Igarashi, K. and Maeno, T. (2018) Phase Unwrapping in Magnetic Resonance Elastography. Open Journal of Medical Imaging, 8, 111-125. https://doi.org/10.4236/ojmi.2018.84014

References

Abbreviations

AA Abdominal Aorta

DE Dilate-Erode

FCT Fast Cosine Transform

IVC Inferior Venacava

IVD Inter Vertebral Disc

LBE Laplacian-Based Estimate

LFE Local Frequency Estimate

MD Minimum Discontinuity

MEG Motion Encoding Gradient

MRI Magnetic Resonance Imaging

MRE Magnetic Resonance Elastography

PM Psoas Major

PMMRE Psoas Major Magnetic Resonance Elastography

PU Phase Unwrapping

RG Region Growing

RO Read Out

SAR Synthetic Aperture Radar