COMPUTER PROCESSING AND MODELING Full Papers

Magnetic Resonance in Medicine 73:2321–2331 (2015)

Real-Time 4D Phase Unwrapping Applied To Magnetic Resonance Elastography Eric Barnhill,1* Paul Kennedy,1 Curtis L. Johnson,2 Marius Mada,3 and Neil Roberts1 Purpose: Phase amplitude is a source of signal in magnetic resonance elastography (MRE) experiments but its exploitation in experimental design has been limited due to the challenges of phase wrap. This study addressed this aspect of MRE through new developments in algorithms, heuristic strategy, and user interface. Methods: A test dataset with systematic variation of three parameters—nested wrap, gradient, and noise level—was developed to choose phase-unwrapping algorithms and to analyze their performance. A new application, PhaseTools, was developed that implemented three phase-unwrapping algorithms that adhere to a “real-time” criterion of less than 3 min for a four-dimensional MRE acquisition. Two of the algorithms extend previously published algorithms and one was newly developed. The algorithms were then applied to five datasets from MRE, two typical cases and three edge cases that were particularly challenging in one of the three parameters. Results: The performance of the PhaseTools algorithms on the test dataset was comparable to two widely cited algorithms that take hours or days to complete. Guidelines for the optimal use of each algorithm are established. Conclusion: PhaseTools enables the substantial increase of signal-to-noise in MRE experiments at negligible additional computational cost. PhaseTools is freely released with this study, making robust real-time phase unwrapping available to any group using phase-based imaging. Magn Reson Med C 2014 Wiley Periodicals, Inc. 73:2321–2331, 2015. V Key words: magnetic resonance elastography; phase unwrapping; brain; muscle; software

INTRODUCTION Phase and Signal in MR Elastography Magnetic resonance elastography (MRE) (1) is a method of phase contrast imaging that enables the determination of

1 Clinical Research Imaging Centre, School of Clinical Sciences and Community Health, College of Medicine and Veterinary Medicine, The University of Edinburgh, Edinburgh, UK. 2 Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA. 3 Wolfson Brain Imaging Centre, Department of Clinical Neurociences, The University of Cambridge, Cambridge, UK Additional Supporting Information may be found in the online version of this article. *Correspondence to: Eric Barnhill, Clinical Research Imaging Centre, The University of Edinburgh, 47 Little France Crescent, Edinburgh EH16 4TJ, Scotland, United Kingdom. E-mail: [email protected] Received 10 February 2014; revised 22 April 2014; accepted 22 April 2014 DOI 10.1002/mrm.25332 Published online 18 June 2014 in Wiley Online Library (wileyonlinelibrary. com). C 2014 Wiley Periodicals, Inc. V

soft tissue viscoelastic properties. Manduca et al. (2) first observed that greater phase amplitude yields greater signal in MRE. However, the “phase unwrapping” problem has been proven fundamentally intractable (3) and reliable phase unwrapping can require extended processing times even in two dimensions (4). A combination of higher signal levels and real-time processing would advance MRE imaging in several ways. Brain scans at higher signal levels contain detail otherwise lost to noise, enabling the detection of subtler neuroplastic changes and a richer viscoelastic brain atlas (5,6), while real-time muscle scans potentially enable an interactive evaluation of functional muscle performance within a single scan session, as well as more accurate determination of intramuscle viscoelastic changes (7), and liver MRE could be acquired at higher signal-to-noise ratios (SNR), without additional masking or other manual intervention. Accurate determination of tissue viscoelastic parameter values with MRE requires the full 3D derivative of the wave motion (8). This necessitates acquisitions of a 3D volume over a multistep time series, resulting in a 4D phase image, usually containing upward of 10 million voxels to be unwrapped. MRE phase images pose several additional challenges at once, including high and varied spatial frequencies, coarse temporal sampling, irregular boundaries, anatomical discontinuities, and signal inhomogeneities. As MRE is intended for clinical use, the increase of phase amplitude loses its clinical utility if phase unwrapping adds hours or days of processing time for a single acquisition. This study developed “real-time” solutions for the 4D phase unwrapping problem to better exploit phase amplitude as a source of signal in MRE experiments. “Realtime” here means less than the scan time for a 4D MRE acquisition. As a high resoution MRE acquisition of the full brain currently takes 8–12 min (6,9), an upper limit was set at 3 min for phase unwrapping processing time of a similar volume. While clinical MRE can be feasible with longer unwrapping times, this “real-time” requirement is an effective benchmark for clinical convenience and additionally identifies which unwrapping algorithms can function effectively in an inline processing pipeline. This study evaluated phase-unwrapping algorithms for their ability to unwrap 4D MRE data robustly in real time. To do this, they must first have good dimensional scaling properties, that is, they must be able to unwrap large numbers of 4D voxels within the processing time limits. Second, algorithm performance must be robust to high levels of three key parameters in MRE phase unwrapping—phase gradient, noise level, and number of

2321

2322

Barnhill et al.

nested wraps. A novel evaluation methodology was developed to test both dimensional scaling and robustness to these three parameters. Algorithms that performed successfully were integrated into a new software application called PhaseTools that is freely released with the study. As other domains of phase imaging face many of the same challenges, PhaseTools makes robust realtime phase unwrapping available to any group using phase-based imaging. METHODS Analytic Test Data A test dataset was developed to choose algorithms for implementation in PhaseTools and to evaluate the chosen algorithms’ performance. The set systematically varies three parameters that cause difficulty in phase unwrapping— phase gradient, noise level, and number of nested wraps— and has similar voxel number to an MRE acquisition. Because so few algorithms have been published in 4D, the test dataset was developed in 3D to facilitate comparisons between existing algorithms. Algorithms of sufficient robustness were additionally subject to a “scaling test” in which the dataset was run as a single, combined 4D volume to determine if the algorithm could perform feasibly in 4D. The test set was modeled as shear wave propagation in an equilibrium Hookean solid, which has as its equation 2 of motion F ¼ r ddtu2 , where r is density and u the displacement field. This equation takes F ¼ sin(mh) for real m as a solution (10). The test dataset was, then, generated on a uniform 64  64  10 grid using:     ðx þ y þ zÞ2p fw ðx; y; zÞ ¼ mod asin þ gðsÞ [1] l where (x, y, z) are image coordinates, fw is the wrapped phase, mod is the modulo operator wrapping the data within the range of 0  fw  2p, a is the amplitude, k is the wavelength, and g(r) is an additive Gaussian noise function with standard deviation r. a, k, and g(r) were all progressively modulated by permuting the following value sets: an ¼ 2pn; 1  n  8 ln ¼ 32n; 1  n  8 sn ¼ ½:0001; :0005; :001; :005; :01 for integer n. This resulted in 320 images of varying amplitude (with higher amplitude creating more nested phase wraps), wavelength (with shorter wavelength creating higher gradients), and Gaussian noise. The 320 64  64  10 images contained a total of 13,107,200 voxels. The dataset was then wrapped modulo 2p and unwrapped by each algorithm and compared to the true data. The dataset was designed to probe the limits of the algorithms, and as the three parameters increase, the wrap of the images gradually transitions from highly challenging to impossible. Evaluation of Analytic Data Four-dimensional scaling test. To evaluate the feasibility of the algorithms’ performance on higher-dimensional datasets, processing times were compared for the unwrap-

ping of the test dataset in 2D (as 3200 64  64 images), in 3D (as 320 64  64  10 images), and in 4D (as one 320  64  64  10 image). Algorithms that could not unwrap the dataset in 4D in real-time were excluded. As few phase-unwrapping algorithms exist in 4D, 3D algorithms were first evaluated by unwrapping the dataset as a single 3D (64  64  3200) image. If this performance was feasible in real-time, a fully 4D version of the algorithm was developed and tested for real-time performance. Robustness test. If an algorithm performed in real-time on the 4D scaling test, then for the second test the unwrapped 3D dataset was evaluated against the true 3D dataset. This evalution was kept to 3D so that the robustness of the PhaseTools algorithms could be compared with the most robust algorithms in the field, even if they could not scale to 4D. The evaluation was performed by measuring two quantities. The first was the variance between the values of the true and the unwrapped volumes after matching the means of the volumes. This is the critical measure when the absolute distribution of the phase values are relevant to the measurement. The second measure was the variance of the phase gradient as computed using a 3D second-order leapfrog finite differences scheme for the gradient (4,11), scaled to amplitude. This measure is the most indicative of success when the unwrapped image is used to study pattern or contour (such as wavelength) and the absolute value distribution of the measurements is not needed. A score of zero indicated a perfect, lossless unwrap, and only occurred in both quantities simultaneously. Mean variance per voxel above 2% was considered an unwrapping failure. If one algorithm contained more error, on all 320 test images, than another algorithm, the former was excluded from further consideration. Comparison with reference algorithms. The algorithms’ variance results were compared with the results of wellestablished optimization-based algorithms—Jenkinson’s PRELUDE algorithm (12) and Constantini’s CUNWRAP algorithm (13) in a 3D implementation—to ensure that the overall robustness was competitive with the best algorithms in the field. The 3D and 4D implementations of CUNWRAP were created for this study, and the source code is available at http://cognitive-eurhythmics.com/ code.html. Only unwrapping algorithms using  Oðn log nÞ methods, such as greedy or fast spectral techniques, were expected to perform in 4D in real time and this is verified by the processing times of PRELUDE and CUNWRAP in the Results section. Consequently other robust optimization-based approaches, e.g., (3,12–17) were excluded from consideration. As the only other article addressing MRE phase unwrapping in particular used an iterative method (18), it too was excluded. Two published algorithms passed the exclusion criteria and were consequently implemented in PhaseTools and used to unwrap the MRE datasets. The first was a Laplacian-based algorithm published in Schofield and Zhu (19) and Volkov and Zhu (20). The second was a region-growing algorithm published in 2D by Herraez et al. (21) and 3D by Abdul-Rahman et al. (22). A third algorithm, called here Dilate-Erode-Propagate, was

Real-Time 4D Phase Unwrapping

2323

developed for the project. The algorithm runs in linear time and outperformed all other real-time algorithms in delivering exact solutions in high noise with low gradients. The three algorithms are described briefly below. Descriptions of Chosen Algorithms Laplacian-Based Estimate (LBE) The Laplacian-based algorithm (19,20) relates wrapped to true phase not in terms of gradients, as is found in, for example, least-squares techniques (4), but in terms of Laplacians. For a wrapped phase field fw , the Laplacian of the true phase, f, is formulated in terms of the wrapped phase: r2 f ¼ cos fw r2 ðsin fw Þ  sin fw r2 ðcos fw Þ

[2]

with r2 the Laplace operator. This allows for an estimate of the true phase, f0 , to be derived by application of an inverse Laplace operator r2 to the right-hand side. Fourier-domain forward and inverse Laplace operators enable a direct solution for f0 in 4D as: f0 ðx; y; z; tÞ ¼ FCT 1   FCT ½cos fw ðFCT 1 ½ðp2 þ q2 þ r 2 þ s2 ÞFCT ðsin fw ÞÞ  p 2 þ q 2 þ r 2 þ s2 FCT 1

  FCT½sin fw ðFCT 1 ½ðp2 þ q2 þ r 2 þ s2 ÞFCTðcos fw ÞÞ p2 þ q2 þ r 2 þ s2

[3] with FCT and FCT1 the forward and inverse fast cosine transforms, (x, y, z, t) the time-domain coordinates and (p, q, r, s) the Fourier domain coordinates. Schofield and Zhu (19) proceed to apply an integer rounding operator to f0 and iterate until a desired level of convergence. However, a proof of the failure of these steps of the algorithm is presented as Supporting Information accompanying the online article, and will not be considered further in the body of the article. Unlike phase unwrapping techniques that rely on local information, the Laplacian-based estimate consists entirely of spectral techniques and is, therefore, considered “noise-immune,” in the sense that disruptions to the unwrap do not propagate spatially in the phase unwrapping process. However, the algorithm is not lossless: inspection of Eq. [3] shows three types of operations: forward and inverse cosine transforms, trigonometric operations, and the masking expression ðp2 þ q2 þ r 2 þ s2 Þ. As the masking terms are in both numerator and denominator, their impact on the spatialdomain value range cancels out. Consequently the unwrapped phase solution is projected onto the original phase range which was restricted by phase wrap to ½p; pÞ. The first effect of this range projection is the inability to recover true phase value as the phase will always be mapped on the range ½p; pÞ. The second effect is that, although the noise does not propagate, it is amplified: additive noise that would have a minor impact if the image were unwrapped to a larger phase range has a much greater impact on the compressed phase range.

RG Algorithm (RG) The second algorithm, published in Herraez et al. (21) and Abdul-Rahman et al. (22), relies on local quantities of voxel edge reliability (sometimes referred to as “pixel quality” or “phase quality”) rather than spectral values. Abdul-Rahman et al. (22) calculate reliability using an L2 norm second differences reliability criterion, which extended to 4D becomes: rðpÞ ¼

40  X 1 ðp  nþ Þ2 þ ðp  n Þ2 2

[4]

n¼1

with r(p) the reliability of the pixel p, nþ a neighbor on the 4D hypercube of neighbors to p, and n the opposite neighbor. In four dimensions each voxel has 80 neighbors, resulting in a reliability criterion that is the sum of 40 measurements. The algorithm grows regions through the progressive unwrapping of edges. The voxels on both sides of the edge are synchronized and their respective regions merged through a union-find routine. Edges are unwrapped with a greedy algorithm but the algorithm proceeds discontinuously in quality order. This limits the impact of noise and discontinuity that often propagate through path-guided algorithms (4) while maintaining near-linear time for the operation. The RG algorithm handles noise and discontinuity by unwrapping the least reliable voxels last; however, the algorithm’s ability to contain disruption depends on the correct formation of groups of pixels within the same phase range. To take a linear-time (i.e., greedy) approach, the algorithm uses a hard-cutoff difference value of p to decide whether a neighboring voxel is in the same group or not. Consequently, the algorithm fails over any given discrete subdomain in which the combination of contiguous phase movement and noise exceed p, or more succinctly, when the subdomain has the property:  þ gðsÞ > p D

[5]

 the mean phase gradient and g(r) the standard with D deviation of the Gaussian noise. If a dimension of the image, such as the volumetric or time-series dimensions, has a particularly high gradient, such as from a coarse time-sample, then an unwrap including that dimension is particularly prone to error: consequently, a sound strategy for use of the algorithm is to unwrap in four, then three, and then if necessary two dimensions to resolve remaining errors and unwrap the image. Dilate-Erode-Propagate (DE) Algorithm A third algorithm that performed robustly in 4D realtime was developed for this study. The DE algorithm propagates the phase of a single slice through the 4D image. This requires a “domain knowledge” step (23), of manually selecting a slice to propagate [another domain knowledge phase-unwrapping algorithm is DENSE (24)]. For the propagation all wrapped regions on the selected slice must be smaller than the chosen dilate-erode window. In images of light wrap, such a slice may already exist and the DE alone can unwrap the image. In more

2324

Barnhill et al.

FIG. 1. Windows from the PhaseTools workflow for the “typical brain” dataset in the study: (b) magnitude and (c) phase DICOM images are (e) masked and (f) normalized using a (d) mean image as reference. The dataset is then unwrapped using the (g) LBE, (h) RG, and (i) DE algorithms. Each successfully unwraps the dataset. All operations are buttons on the (a) PhaseTools GUI.

complex images, one slice can be preconditioned until the wrap satisfies this criterion. For the robustness test, one slice in each image was manually selected, preconditioned if necessary, and chosen to propagate through the image. Unlike the RG which needs to form groups, the DE dilates the most reliable voxels by synchronizing neighbors to them, then erodes the least reliable by synchronizing them to their neighbors. This allows for a resolution of less reliable voxels and clusters that is entirely data-driven, eliminating the need for any criterion for region coherence. Consequently, the algorithm is more robust in high noise that might disrupt regionforming. Pseudocode for the DE is as follows: 1. Make 2n – 1 passes evaluating the voxels with increasing windows of radius 1 to n, then decreasing windows of n – 1 back to 1 (with n useradjustable through the interface).

a. Dilate i. Evaluate the reliability of each voxel within the current window. Store voxels by reliability in a heap. ii. Until the heap is empty, extract the next voxel of highest reliability. Synchronize the phase of its eight neighbors by applying a modulo operator. Delete the neighbor voxels from the heap. (Eliminating the neighbors from further consideration ensures that moderately reliable voxels only dilate in the absence of adjacent, more reliable neighbors.) b. Erode i. Re-evaluate voxel quality, storing in a second heap. ii. Until the heap is empty, choose the next voxel of lowest reliability. If four edge neighbors around a center voxel agree with each other, but disagree with the center, synchronize the center to the neighbors.

Real-Time 4D Phase Unwrapping

2325

2. If the slice is not a boundary, propagate: a. Each voxel pn on the slice synchronizes the corresponding voxel pn þ 1 on the next slice with a modulo operator. b. pn þ 1 is evaluated with a second-differences quality measure (except the initial propagation, when a first-differences measure is used). If the value is above a user-adjustable cutoff, pn þ 1 is used to synchronize pn þ 2. If it does not make the cutoff, the new voxel is skipped in the following propagation and pn is used to synchronize pn þ 2. The error criterion is similar to the RG but limited by the distribution of the noise over the dilation/erosion window rather than by the grouping threshold:  þ gðsÞ D >p r2

[6]

with r the radius of the dilation/erosion window, which is set by user adjustable cutoff. For this study, this cutoff was set to 5, which was found to effectively remove small patches of broken wrap with minimal sacrifice of algorithm performance, and this window is used in the performance evaluation times. As the “domain knowledge” aspect of the algorithm will depend on user input, it was not controlled for experimentally in the study. Instead, the software package that is released in conjunction with the study contains a tutorial and screencast in which the DE is applied, in conjuction with manual slice editing, to unwrap the phase of a challenging phantom acquisition. The manual editing in the tutorial takes less than a minute and in the empirical experience of the authors, the DE consistently outperforms other algorithms on high noise wrap with less than 5 min editing.

Table 1 Performance Comparison for the Various Algorithms on a Dell E6430 Laptop Running 8 GB Memory Average time Algorithm

MRE Data Evaluation Five 4D MRE datasets from three different laboratories were used to evaluate the real-world performance of the PhaseTools algorithms:  Typical brain Axial brain acquisition, vibration frequency 62.5 Hz, 128  128  5 with 16 time steps and three displacement vectors, acquired with an EPI sequence at the University of Cambridge. This acquistion was low in wrap, gradient, and noise.  Edge case brain A second axial acquisition with the same specifications except vibrated at 25 Hz. This acquisition had high gradient, up to four nested

3D

4D

34 s 42 s 1 min 30 s

40 s 2 min 42 s 1 min 15 s

1 h 43 min 1 h 46 min

N/Ab N/Ac

Processing times for the  13 M pixel test dataset processed as 2D, 3D, and 4D volumes a DE maximum window radius was set to 5 pixels. b PRELUDE algorithm failed to terminate after 7 days running the dataset as a single 3D image. c CUNWRAP 4D algorithm ran out of memory even when increased to 24 GB.







PhaseTools Implementation PhaseTools was written in Java as a plug-in for the ImageJ platform and can also be run within MATLAB through the use of MIJ (http://bigwww.epfl.ch/sage/soft/mij/). The PhaseTools graphical user interface is shown in Figure 1. In addition to an implementation of the phaseunwrapping algorithms, a suite of support tools includes manual shifting of selections by 62p, masking, normalizing, and conversion between Zero and NaN. PhaseTools is released with this study and can be downloaded from http://cognitive-eurhythmics.com/code.html.

2D

PhaseTools algorithms LBE 32 s RG 49 s N/A DEa Reference algorithms PRELUDE 8 min 58 s CUNWRAP 8 min 30 s



wraps, high noise, and phase disruptions at the ventricles. Typical muscle Transverse thigh acquisition, 50 Hz, three displacement vectors, 128  128  5 with eight time steps and three displacement vectors, acquired with a 3D EPI sequence at the University of Edinburgh (7). Acquisition had moderate noise and up to one nested wrap. Edge case muscle Similar acquisition to typical muscle except in the sagittal plane. Acquisition had high noise, very low gradients, and four nested wraps. Liver An axial liver acquisition acquisition, 62.5 Hz, 128  128  5 with 16 time steps and one displacement vector, acquired with an EPI sequence at the University of Cambridge. This scan incorporates challenges unique to abdominal imaging, including phase disruption on much of the 3D boundary of the organ. Edge case phantom An agarose phantom acquisition, 100 Hz, 128  128  10 with 16 time steps and three displacement vectors, acquired with a spiral sequence at the University of Illinois (9). In contrast to the previous image, this edge case had very high gradients with low noise and up to three nested wraps.

The masked wrapped phase images from the five datasets are shown in column 1 of Figure 4. RESULTS Performance Comparison of the Algorithms The processing times of the five studied algorithms, running on a Dell Latitude E6430 notebook with 8 GB of RAM, is shown in Table 1. Test Data Evaluation A scatterplot of the performance of the LBE, RG, and DE algorithms, with axes for wavelength, amplitude, and noise level, and points colored by variance, is shown in Figure 2. The color table is a rainbow with red

2326

Barnhill et al.

FIG. 2. Three-dimensional variance scatterplot of unwrapping algorithms under systematic alteration of spatial frequency, amplitude, and Gaussian noise level. PhaseTools Algorithms: (a) LBE, (b) RG, (c) DE. Reference algorithms: (d) PRELUDE, (e) CUNWRAP. Red is lowest variance, with diamonds indicating a perfect unwrap, while blue represents most variance, with dark blue a cutoff of 2% mean variance per image per voxel.

representing normalized variance close to zero, and blue representing a cutoff of mean variance per voxel of >2%. The zero-variance perfect unwraps are additionally represented as diamonds, while imperfect unwraps are represented by circles. Plots were generated in R (25) using the scatterplot3d (26) library. A summary of the error of each algorithm in categories of severity is shown in Table 2. The PhaseTools algorithms each show different strengths and weaknesses. The LBE shows the lowest accuracy of

values but shows the highest accuracy of high-noise gradients. The RG and DE have roughly the same number of exact solutions, with the DE additionally having more accurate gradients in higher noise. The reference algorithms contain more exact solutions than any of the PhaseTools algorithms but are roughly as robust to gradient as the DE and less robust to gradient in high noise than the LBE. Sample slices for three test images, each a result in which a different algorithm performed best, are shown in Figure 3.

Real-Time 4D Phase Unwrapping

2327

Table 2 Summary Statistics for Exact Test Data

Value variance LBE RG DE PRELUDE CUNWRAP Gradient variance LBE RG DE PRELUDE CUNWRAP

Perfect (0%)

Minor error (< 0.1%)

Some error (0.1 – 2%)

Failure (>2%)

0% 34.1 35.9 46.9 40.0

5% 12.5 20.6 15.6 17.5

62.5% 27.2 17.5 17.5 21.3

32.5% 32.2 25.9 20.0 21.2

0% 31.6 33.1 44.7 36.6

13.1% 13.1 17.2 13.1 17.8

79.7% 38.1 39.7 32.2 29.7

7.2% 17.2 10.0 10.0 15.9

wavelength. In this image, the LBE image is dominated by noise, and the noise disrupts the RG from forming regions. The DE performs best as the propagation of a manually validated slice was robust to noise. MRE Data Results

Slices were chosen with a random number generator. Each image also contains a yellow plot line, and the results from the line plots are included as Supporting Information accompanying the online article for readers who wish to see the algorithm performance in more detail. Test image 45 (a ¼ 1, l ¼ 2, s ¼ 0.0001). Image 45 has light noise, moderate wavelength, and one phase wrap in the amplitude. In this image, both the RG and DE perform best delivering exact solutions. While the LBE, RG, and DE all recover phase contour accurately, the LBE fails to recover values due to the phase range projection detailed in methods. Test image 122 (a ¼ 4, l ¼ 1, s ¼ 0.0005). Image 122 has high amplitude, heavy noise, and high gradient. In this image, the LBE performs best, maintaining the contour accurately. However, the values are projected to the range ½p; þpÞ, and the noise in the image has grown relative to the signal, in comparison with the true values. Test image 256 (a ¼ 7, l ¼ 4, s ¼ 0.001). Image 256 has very high noise, very high amplitude, and broad

Single slices for the MRE results are shown in Figure 4. For typical cases, slices were selected by a random number generator as all images unwrapped perfectly. For edge cases, slices were chosen to illustrate particular aspects of the algorithms’ performance. Running times for the MRE data using the same memory specifications described above are shown in Table 3. a. Typical Brain As columns 3 and 4 of row (a) show, both RG and DE unwrapped the acquisition without error, while the LBE successfully recovered contour as seen in column 2. However, the range projection of the LBE has diminished the impact of lower spatial frequencies. b. Edge Case Brain The LBE resolved the image effectively with some increase in noise as shown in row (b), column 2. On some slices the ventricle phase disrupted the RG’s region-forming criterion as shown in column 3 of row (b). Using the RG as preconditioner, the DE resolved these slices by propagating unwrapped slices through slices containing wrap as seen in column 4. c. Typical Muscle Columns 2, 3, and 4 of row (c) show that all three algorithms unwrap the muscle successfully. As the phase did not extend much beyond 62p, the LBE range projection only slightly compressed the range of the result. d. Edge Case Muscle This image had results similar to test image 256. The combination of low gradient and high noise caused the LBE to produce an image dominated by noise as shown in column 2 of row (d). The RG had occasional disruptions to its regions from the noise as shown in column 3. With

FIG. 3. Sample slices of (a) test image 45, (b) test image 122, (c) test image 256. Row (a) shows strong results from all three algorithms (LBE, RG, and DE) while in (b) the LBE performs best and in (c) the DE performs best.

2328

Barnhill et al.

FIG. 4. Slices from MRE test datasets: (a) typical brain (b) edge case brain (c) typical muscle (d) edge case muscle (e) liver (f) edge case phantom. In typical cases (a) and (c) all three algorithms (LBE, RG, and DE) perform strongly. In (b), both LBE and DE better handle phase disruptions from the ventricles. In (d), the DE unwraps most successfully with broad waves and high noise. In (e), all algorithms unwrap the liver, with the RG occasionally disrupted by residues outside the liver, which can be resolved through use of the DE. In (f), only the LBE unwraps the high gradients successfully.

Table 3 Performance Comparison for the Various Algorithms on a Dell E6430 Laptop with 8GB RAM

Size(px) LBE RG DE

Typical Brain

Edge case Brain

Typical Muscle

Edge case Muscle

Edge case Phantom

Liver

3,923,160 16.3 s 7.2 s 10.1 s

3,923,160 16.2 s 51.9 s 98.9 sa

1,966,080 7.3 s 6.0 s 24.6 s

1,966,080 7.5 s 20.1 s 40.4 sa

3,923,160 14.5 s 21.5 s 76.4 sa

1,310,720 4.7 s 13.5 s 47.0 s

Processing times for MRE datasets. Times in seconds. DE was used with RG as preconditioner and value is total processing time for both.

a

Real-Time 4D Phase Unwrapping

2329

FIG. 5. PhaseTools algorithms’ performance as a group versus reference algorithms’ performance as a group. Each data point is colored by the group membership of the best performing algorithm. If the best unwrap came from the PhaseTools group, the point is violet, and if it came from the Reference group, it is green.

the RG as a preconditioner, the DE algorithm successfully propagated an unwrapped slice through the image as shown in column 4. e. Liver Due to proximity of the lungs, aorta, and other organs, regions on the liver boundary create phase disruptions that can interfere with the unwrap across the liver. All algorithms effectively unwrapped the liver tissue, but several slices on the RG were disrupted by phase residues outside the liver. Column 3 of row (e) shows an example of a phase residue created by signal incoherence outside the liver that the RG resolves by leaving phase wrap across the liver. In row (f) column 4, unwrap across the liver has been resolved by propagating a coherent unwrap of the liver through the volume using the DE. f. Edge Case Phantom This image had results similar to test image 122. The high gradients prevented the region-forming criteria from forming properly, and the RG and DE solutions contain as much wrap as the original image as shown in columns 3 and 4 of row (f). In contrast the LBE performed very accurately, as shown in column 2, with the range projection impacting the image little due to low noise and moderate overall phase range.

DISCUSSION Guidelines for Use of the LBE, RG, and DE Algorithms In both test and real-world data, the three algorithms implemented in PhaseTools have complementary sets of strengths and weaknesses: RG Algorithm The RG is a fully 4D implementation of an algorithm that slightly underperforms two of the most widely-cited reference algorithms, terminates in seconds rather than hours or days, and makes much lower memory demands. This makes the PhaseTools’ RG implementation a good “first stop” for MRE phase unwrapping, as is confirmed in its perfect unwrap for the typical cases, Figure 4a,b, column 3. Rapid changes in gradient can cause its

region-forming criterion to err, especially along more coarsely sampled dimensions; therefore, a best practice is to first run the algorithm in 4D, then with the majority of the phase aligned, run it a second time in 3D, and then a third time in 2D if necessary. As Table 3 shows, on a real-world image even all three runs can be expected to take less than a minute. Laplacian-Based Estimate If rapid gradient and excessive discontinuity prevent the RG from functioning well, the LBE can be expected to produce a strong solution. As the LBE projects results into a single phase range ½p; þpÞ, it is never perfect even in trivial cases. However, when the overall phase range is not much greater than ½p; þpÞ, the solution is highly accurate, produces no residues, and is immune to disruptions from noise, discontinuity, residues, and aliasing. Solutions to difficult edge cases such as those seen in test image 122 and Figure 4b, e, column 2, demonstrate why this algorithm outperforms even the reference algorithms in the most challenging data. Furthermore, its projection effects can be a useful preconditioner when low frequencies need to be removed in later stages of image processing. DE Algorithm If in contrast the image has moderate or low gradients but is disrupted by high noise, the DE is an optimal solution as it does not rely on region-forming criteria that can be disrupted by noise, and will not amplify the noise in the manner of the LBE. This is shown in test image 256 as well as Figure 4d, f, column 4. In these cases, the DE requires a second algorithm for preconditioning. In any image where at least one slice has sufficiently light wrap such that all wrapped clusters are smaller than the dilate-erode window, the DE is a robust standalone unwrapper, as shown in the MRE typical cases. In summary, each algorithm can be viewed as possessing two of three key traits: delivery of exact solutions (RG and DE), robustness in high noise (DE and LBE), or fully automated performance (RG and LBE). Optimal choice of algorithm for an MRE experiment would

2330

depend on which of the three algorithms is best equipped to address the properties of the images. Comparison of PhaseTools Algorithms and Reference Algorithms Dimensional Scaling Using techniques with running time of  Oðn log nÞ, the LBE and DE algorithms scale into higher dimensions effectively without performance loss. The RG contains one polynomial-time operation that caused an increase of performance time in higher dimensions, but remaining well within the real-time criterion of 3 min. In the realworld MRE acquisitions, all PhaseTools algorithms executed 4D unwraps in under 20 s. By contrast, the optimization-based algorithms demonstrated their difficulties in scaling beyond two dimensions, with a 4D unwrap of the test dataset proving infeasible for PRELUDE or CUNWRAP. This result corroborated previous reports in the literature of PRELUDE failing to terminate on sufficiently large or complex MRI images such as Witoszynskyj et al. (16). The use of optimization-based phase-unwrapping algorithms for MRE, then, requires forcing of phase unwrapping along a lower-dimensional subspace of the set. This results in a lack of alignment along the remaining dimensions, and in the case of high phase wrap, subsequent alignment along the remaining dimensions is nontrivial. Robustness For an overall robustness comparison between PhaseTools and reference algorithms, the lowest-variance result for either the PhaseTools group of algorithms, or the reference group, for each data point, was plotted in Figure 5. Both groups tied on the simpler wraps, those most likely to be encountered in MRE analysis. The reference algorithms outperformed the PhaseTools algorithms, in 3D without reference to running time, in some conditions of higher gradient and amplitude. However, the PhaseTools algorithms frequently outperformed the reference algorithms, especially in conditions of high noise.

Barnhill et al.

waves, wave scattering, and noise; depending on the inversion approach, the LBE’s suppression of low frequencies may aid inversion accuracy over the use of lossless unwraps, or its noise enhancement may require alteration of other noise suppression parameters. The real-time performance of PhaseTools enables more feasible implementation of scanner-based inline phase unwrapping, and the PhaseTools LBE has already been implemented inline in muscle elastography software for the Siemens IDEA platform (27). Given the advantages of GPUs for similar computations on large data blocks, it would be expected that GPU implementation would further enhance performance, and early investigations running ImageJ via Bumblebee Optirun (http://bumblebeeproject.org/) show additional performance benefits. Future work will develop releases of PhaseTools for a variety of scanner-based MRI processing applications, and incorporate the use of graphics processing units in offline processing. Despite the wide range of published phase-unwrapping algorithms, available code and metastudies are scarce, and research groups still frequently rely on the algorithms summarized 15 years ago in Ghiglia and Pritt (4). In that time, the need for high-volume phase unwrapping, often in 4D, has grown substantially. PhaseTools’ fully 4D capabilities will be of benefit to 4D flow imaging techniques such as 4D-PC-MRI (28) and 4D cardiac imaging such as 3D cine DENSE (29). Phase unwrapping in B0 field mapping has been described as artifact prone [e.g., Robinson and Jovicich (30)]; the noise-immune LBE, or the incorporation of domain knowledge with the DE, may provide more rapid resolution of phase wrap difficulties. Finally, the PhaseTools algorithms may prove of benefit to imaging techniques such as susceptibility weighted imaging (SWI), in which rapid inline phase image processing is required, enabling the elimination of lossy techniques such as high pass filtering of the original complex data (31). Overall, new clinical applications for the phase image continue to be developed rapidly; with PhaseTools, groups needing an efficient high-dimensional solution to the phase unwrapping problem will have a new freely available tool with which to aid their research.

CONCLUSIONS MRE acquisitions similar to the typical cases found in Figure 4a, c are the current common practice, with greater wrap avoided due to the risk that the scan would have to be abandoned due to intractable wrap. With PhaseTools MRE, scientists can now have confidence that wraps approaching the edge cases in Figure 4b, d, and e will deliver valid and useful data. These higher-wrap edge cases can be further expected to contain more signal and, therefore, more clinically useful detail. Future work will investigate the impact of increased phase amplitude and choice of unwrapping algorithm on MRE elasticity reconstruction maps. The projection effect of the LBE does not impact phase contour or wavelength; therefore, it is not expected to impact results from inversion methods that use wavelength estimation such as local frequency estimation or direct inversion (2). However, MRE inversion approaches vary in terms of other preprocessing steps including the suppression of bulk

ACKNOWLEDGMENT The authors wish to acknowledge discussions with Dr. Penny Davies, University of Strathclyde and Dr. Miguel Herraez, University of Valencia. The imaging research project was carried out at the Clinical Research Imaging Centre, Edinburgh (www.cric.ed.ac.uk), University of Edinburgh, which is part of the SINAPSE (Scottish Imaging Network - A Platform for Scientific Excellence) collaboration (www.sinapse.ac.uk) funded by the Scottish Funding Council and the Chief Scientist Office. Support from NHS Lothian R&D, the WTCRF are gratefully acknowledged. Eric Barnhill holds a Scottish Universities Physics Alliance (SUPA) Industrial Prize Studentship and is grateful for financial support from The Mentholatum Company. Paul Kennedy holds a SINAPSE PhD studentship. Neil Roberts is a member of The University of Edinburgh Centre for Cognitive Ageing and Cognitive Epidemiology, part of the cross council

Real-Time 4D Phase Unwrapping

Lifelong Health and Wellbeing Initiative (MR/K026992/ 1). Funding from the BBSRC and MRC is gratefully acknowledged. REFERENCES 1. Muthupillai R, Lomas DJ, Rossman PJ, Greenleaf JF, Manduca A, Ehman RL. Magnetic resonance elastography by direct visualization of propagating acoustic strain waves. Science 1995;269:1854–1857. 2. Manduca A, Oliphant T, Dresner M, Mahowald J, Kruse S, Amromin E, Felmlee J, Greenleaf J, Ehman R. Magnetic resonance elastography: noninvasive mapping of tissue elasticity. Med Image Anal 2001;5:237–254. 3. Chen C, Zebker H. Network approaches to two-dimensional phase unwrapping: intractability and two new algorithms. J Opt Soc Am A 2000;17:401–414. 4. Ghiglia DC, Pritt MD. Two-dimensional phase unwrapping: theory, algorithms, and software. New York: Wiley; 1998. 5. Johnson CL, McGarry MD, Gharibans AA, Weaver JB, Paulsen KD, Wang H, Olivero WC, Sutton BP, Georgiadis JG. Local mechanical properties of white matter structures in the human brain. NeuroImage 2013;79:145–152. 6. Guo J, Hirsch S, Fehlner A, Papazoglou S, Scheel M, Braun J, Sack I. Towards an elastographic atlas of brain anatomy. PloS One 2013;8: e71807. 7. Barnhill E, Kennedy P, Hammer S, van Beek EJ, Brown C, Roberts N. Statistical mapping of the effect of knee extension on thigh muscle viscoelastic properties using magnetic resonance elastography. Physiol Meas 2013;34:1675–1698. 8. Sinkus R, Lorenzen J, Schrader D, Lorenzen M, Dargatz M, Holz D. High-resolution tensor MR elastography for breast tumour detection. Phys Med Biol 2000;45:1649–1664. 9. Johnson CL, McGarry MDJ, Van Houten EEW, Weaver JB, Paulsen KD, Sutton BP, Georgiadis JG. Magnetic resonance elastography of the brain using multishot spiral readouts with self-navigated motion correction. Magn Reson Med 2013;70:404–412. 10. Auld BA. Acoustic fields and waves in solids, Vol. 1. New York: Wiley; 1973. 11. Grossmann C, Roos H-G, Stynes M. Numerical treatment of partial differential equations. Heidelberg, Germany: Springer; 2007. 12. Jenkinson M. Fast, automated, N-dimensional phase-unwrapping algorithm. Magn Reson Med 2003;49:193–197. 13. Costantini M. A novel phase unwrapping method based on network programming. IEEE Trans Geosci Remote Sensing 1998;36:813–821. 14. Cusack R, Papadakis N. New robust 3-D phase unwrapping algorithms: application to magnetic field mapping and undistorting echoplanar images. Neuroimage 2002;16:754–764.

2331 15. Bioucas-Dias JM, Valad~ ao G. Phase unwrapping via graph cuts. IEEE Trans Image Process 2007;16:698–709. 16. Witoszynskyj S, Rauscher A, Reichenbach JR, Barth M. Phase unwrapping of mr images using f un–a fast and robust regiongrowing algorithm. Med Image Anal 2009;13:257–268. 17. Langley J, Zhao Q. A model-based 3d phase unwrapping algorithm using gegenbauer polynomials. Phys Med Biol 2009;54:5237–5252. 18. Wang H, Weaver JB, Perreard II, Doyley MM, Paulsen KD. A threedimensional quality-guided phase unwrapping method for MR elastography. Phys Med Biol 2011;56:3935–3952. 19. Schofield MA, Zhu Y. Fast phase unwrappng algorithm for interferometric applications. Opt Lett 2003;28:1194–1196. 20. Volkov VV, Zhu Y. Deterministic phase unwrapping in the presence of noise. Opt Lett 2003;28:2156–2158. 21. Herr aez MA, Burton DR, Lalor MJ, Gdeisat MA. Fast two-dimensional phase-unwrapping algorithm based on sorting by reliability following a noncontinuous path. Appl Opt 2002;41:7437–7444. 22. Abdul-Rahman H, Gdeisat M, Burton D, Lalor M. Fast threedimensional phase-unwrapping algorithm based on sorting by reliability following a non-continuous path. Proc SPIE 2005;5856:33. 23. Cook DJ, Holder LB, Djoko S. Scalable discovery of informative structural concepts using domain knowledge. IEEE Expert 1996;11: 59–68. 24. Spottiswoode BS, Zhong X, Hess A, Kramer C, Meintjes EM, Mayosi BM, Epstein FH. Tracking myocardial motion from cine dense images using spatiotemporal phase unwrapping and temporal fitting. IEEE Transa Med Imaging 2007;26:15–30. 25. R Core Team. R: a language and environment for statistical computing. R foundation for statistical computing, Vienna, Austria, 2013. 26. Ligges U, M€ achler M. Scatterplot3d-an R package for visualizing multivariate data. J Stat Softw 8(11);1–20. Univ., SFB 475, 2002. 27. Barnhill E, Kennedy P, Brown C, van Beek E, Roberts N. Real-time MR elastography of contracted skeletal muscle. In Proceedings from the 25th Annual Meeting of the European Congress of Radiology, Vienna, Austria, 2014, Abstract Number C-1442. 28. Markl M, Frydrychowicz A, Kozerke S, Hope M, Wieben O. 4D flow MRI. J Magn Reson Imaging 2012;36:1015–1036. 29. Auger DA, Zhong X, Epstein FH, Meintjes EM, Spottiswoode BS. Semi-automated left ventricular segmentation based on a guide point model approach for 3D cine dense cardiovascular magnetic resonance. J Cardiovasc Magn Reson 2014;16:1–12. 30. Robinson S, Jovicich J. B0 mapping with multi-channel RF coils at high field. Magn Reson Med 2011;66:976–988. 31. Haacke EM, Mittal S, Wu Z, Neelavalli J, Cheng Y-C. Susceptibilityweighted imaging: technical aspects and clinical applications, Part 1. Am J Neuroradiol 2009;30:19–30.