Quantitative Susceptibility Map Reconstruction Using Annihilating Filter-based Low-Rank Hankel Matrix Approach

Quantitative Susceptibility Map Reconstruction Using Annihilating Filter-based Low-Rank Hankel Matrix Approach

Hyun-Seo Ahn Department of Bio and Brain Engineering, Korea Advanced Institute of Science and Technology, Daejeon, South Korea Sung-Hong Park sunghongpark@kaist.ac.kr Jong Chul Ye jong.ye@kaist.ac.kr

Quantitative susceptibility mapping (QSM) inevitably suffers from streaking artifacts caused by zeros on the conical surface of the dipole kernel in k-space. This work proposes a novel and accurate QSM reconstruction method based on a direct k-space interpolation approach, avoiding problems of over smoothing and streaking artifacts. Inspired by the recent theory of annihilating filter-based low-rank Hankel matrix approach (ALOHA), QSM reconstruction problem is formulated as deconvolutional problem under low-rank Hankel matrix constraint in the k-space. To reduce the computational complexity and the memory requirement, the problem is formulated as successive reconstruction of 2-D planes along three independent axes of the 3-D phase image in Fourier domain. Extensive experiments were performed to verify and compare the proposed method with existing QSM reconstruction methods. The proposed ALOHA-QSM effectively reduced streaking artifacts and accurately estimated susceptibility values in deep gray matter structures, compared to the existing QSM methods. Our suggested ALOHA-QSM algorithm successfully solves the three-dimensional QSM dipole inversion problem without additional anatomical information or prior assumption and provides good image quality and quantitative accuracy.

Quantitative susceptibility mapping, Dipole inversion, Low-rank Hankel matrix completion
journal: NeuroImage

1 Introduction

Magnetic susceptibility provides physiological property of tissue different from that of conventional MRI contrasts such as T1w, T2w, and proton density. Magnetic susceptibility is a potential biomarker providing clinically relevant information including tissue iron concentration (Bilgic et al., 2012) and venous oxygen saturation level (Haacke et al., 2010). Also, distribution of susceptibility is related to tissue abnormality, which can be used as diagnosis tool for neuro-degenerative disease (Acosta-Cabronero et al., 2013) and neuro-surgical planning (Eskreis-Winkler et al., 2017). Phase imaging and susceptibility weighted imaging (SWI) have been used to investigate susceptibility of tissues (Haacke et al., 2004); however, susceptibility signals in these methods are not quantitative and are slightly different from the original susceptibility sources. Quantitative susceptibility mapping (QSM) has been developed to map the original susceptibility sources quantitatively in their exact shapes and locations (De Rochefort et al., 2008). QSM is widely used since it can be acquired with conventional gradient-echo imaging protocols and offer quantitative values of the tissue susceptibility directly with no additional scan.

From the physical viewpoint, every susceptibility source builds dipole-shaped local phase distribution when placed in a strong magnetic field. Phase images can thus be expressed as the convolution of the susceptibility source and the dipole kernel. This complicated convolutional relationship between phase images and susceptibility images can be simplified as a simple division after the Fourier transform (Salomir et al., 2003; Marques and Bowtell, 2005). However, in Fourier domain, the dipole kernel generates signal voids near the surface forming an angle of 54.7 along the axis (magic angle), which makes the problem ill-posed. The ill-posed problem generated by the dipole kernel can be solved completely by acquiring additional data at three different angles (Liu et al., 2009; Wharton and Bowtell, 2010); however, it has practical limitations on a clinical application. In many studies, it has been tried to attain QSM images from data acquired along single head orientation by simply replacing k-space following the magic angle (the conical surface) with tuning parameters (threshold based k-space division; TKD (Shmueli et al., 2009; Schweser et al., 2013)), using intrinsic information of gradient echo phase images for homogeneity reconstruction (homogeneity-enabled incremental dipole inversion; HEIDI (Schweser et al., 2012)), exploiting sparsity of its estimated susceptibility maps (compressed sensing compensation; CSC (Wu et al., 2012)), or using structural edge information of magnitude images as an additional regularization term (morphology enabled dipole inversion; MEDI (Liu et al., 2012a, b)). However, current QSM methods still suffer from streaking artifacts, over-smoothing, and/or magnitude image induced errors (Wang and Liu, 2015).

Recently, annihilating filter-based low-rank Hankel matrix approach (ALOHA) was proposed for sparse MR image reconstruction (Jin et al., 2016; Ye et al., 2017). Conventional image reconstruction methods based on compressed sensing theory reconstruct images from down-sampled data by exploiting sparse representation of the data in a specific transform including wavelet transforms, and total variation (TV), etc (Lustig et al., 2008). ALOHA is also based on the sparsity of the data; however, instead of using them as regularization term, ALOHA translates the issue of image reconstruction into k-space interpolation using the low-rank weighted Hankel matrix (Jin et al., 2016; Ye et al., 2017). ALOHA showed superior performance over conventional CS methods, and it has many applications such as acceleration of MRI image acquisition (Jin et al., 2016), EPI ghost correction (Lee et al., 2016b), MR parameter mapping (Lee et al., 2016a), MRI artifact correction (Jin et al., 2017), and image inpainting (Jin and Ye, 2015). However, k-space interpolation algorithms including ALOHA have never been applied to solving the dipole inversion problem of QSM. Also, ALOHA requires a large amount of memory and computational resources to store and decompose Hankel matrix; therefore, it has been mainly applied for two-dimensional datasets.

In this study, we propose a new k-space based QSM image reconstruction method, referred to as ALOHA-QSM, to obtain QSM images with relatively high accuracy and suppression of streaking artifacts by applying ALOHA for 3-D phase images. ALOHA-QSM has been evaluated using a public domain numerical phantom and human brain data as well as additional human brain data acquired from five subjects. Extensive comparison results with other existing QSM methods are also provided.

2 Theory

2.1 Dipole inversion

When a substance is placed in a magnetic field, it tends to repel or assist the surrounding magnetic field, a phenomenon called magnetic susceptibility (). Specifically, paramagnetism () and diamagnetism (), respectively, tend to attract and repel surrounding magnetic fields. With a gradient echo sequence, MR signal phases induced by this field distribution evolve during the echo time, resulting in a dipole-shaped phase distribution around the susceptibility source.

Hence, the relationship between the magnetic susceptibility source and its resultant phase images can be expressed as a convolution as follows:


where is the phase signal generated by the magnetic susceptibility source, is the gyromagnetic ratio, is strength of the applied magnetic field, TE is an echo time of the gradient echo, is a convolution operator, and is the dipole kernel which represents local magnetic field shift induced by the susceptibility source, given by


where is the angle between and the direction of main magnetic field.

Based on Eq. (1), magnetic susceptibilities can be estimated through the deconvolution of the phase image with the dipole kernel. However, as all of these calculations are based on 3-D, it is computationally too demanding to take deconvolution on Eq. (1) directly. Instead, by taking the Fourier transform, we can change this problem into the simple division problem on the Fourier domain. For simplicity, we will denote Fourier domain of phase images and susceptibility images as ”k-space” and their axes as in k-space and in the image domain. Hence, the dipole kernel in Eq. (2) can be transformed into the simple k-space filter as follow:


where is a k-space vector , and . Also, by denoting , the relationship between the phase and the magnetic susceptibility in Eq. (1) can be simplified into:


where and are the k-space representation of the susceptibility image and phase image, respectively. However, becomes zero on the conical surface: or , on which the direct inversion using Eq. (4) is not valid. Many studies have conducted to solve this problem. Calculation of susceptibility through multiple orientation sampling (COSMOS) is the most representative and obvious method to solve the problem (Liu et al., 2009). It oversamples data using multiple head orientations and eliminate all the null points at the conical surface. COSMOS can serve as the ideal solution for the dipole inversion problem. However, it requires patients to be scanned at least three times with different head orientations inside the magnet bore. From this reason, it is not easy to implement COSMOS clinically.

There are many approaches to address this problem using the data acquired at a single head orientation. Truncated k-space division (TKD, (Shmueli et al., 2009; Schweser et al., 2013)) is a simple and intuitive method. Some portions of the dipole kernel filter near the conical surface is substituted into the simple threshold value () as follows:


However, it causes discontinuity and vanishing surfaces on the resultant k-space, which leads to streaking artifacts in the corresponding image domain. Also, susceptibility values are often underestimated when a high threshold value is used (Schweser et al., 2013).

Alternatively, additional information from the data can be exploited to handle this uncertainty near the conical surface. One of the popular and powerful way is to use magnitude images. Morphology enabled dipole inversion (MEDI) assumes that structures of susceptibility sources are well reflected on the magnitude information (Liu et al., 2011, 2012a). Specifically, the edges of the susceptibility sources are assumed to coincide with those on the magnitude images:


where is the weighting mask derived from the magnitude image, is the inverse Fourier-transform of the matrix, is the regularization parameter, is a three-dimensional gradient operator, and is a binary mask of the magnitude edge information. Since MEDI uses fully-sampled magnitude images, most streaking artifacts are removed well. Streaking artifacts in susceptibility images can be also minimized by subtracting streaking artifacts iteratively (iLSQR) (Li et al., 2011, 2015). This method assumes that the streaking artifacts are originated from inaccuracies of inversion at the truncated k-space regions. By iteratively subtracting this artifact-originated image, artifact-reduced susceptibility image is obtained.

These aforementioned algorithms are reported to estimate susceptibility values well (Wang and Liu, 2015). However, they still have many limitations. K-space based methods such as TKD or iLSQR show remaining streaking artifacts in images. For MEDI, magnitude images may induce some bias originated from the difference between the magnitude and susceptibility images.

2.2 Deconvolution for QSM using ALOHA

Recently, compressed sensing (CS) based image reconstruction algorithms are widely used to solve the dipole inversion problem in QSM (Wu et al., 2012), and they have shown good performance (Langkammer et al., 2017). CS is one of the most powerful signal reconstruction methods that can recover accurate MR images from down-sampled k-space measurements (Lustig et al., 2008). It exploits the sparsity of the data in some transform domain. In the case of MR image reconstruction, the total variation (TV) or the wavelet transform domain is mainly used for acceleration of data acquisition typically by a factor of 2-5.

An important requirement for applying CS algorithm on MR image reconstruction is that the artifacts originated from the k-space downsampling are distributed incoherently. From this reason, Gaussian random or Poisson disc shaped sampling patterns are often used in CS algorithms. However, the dipole inversion in QSM is far from this requirement: the under-determined regions of the dipole inversion filter in Fourier domain () is clearly form a conical shape. Fortunately, its under-sampling ratio is relatively low (about 10%) compared to conventional CS acceleration problems (over 50%), and a previous study showed CS could successfully reconstruct the susceptibility image from the phase image (Wu et al., 2012).

In contrast to the conventional CS approach, k-space direct interpolation methods determine missing k-space data using intrinsic structural information of the k-space. K-space interpolation methods including simultaneous auto-calibration and k-space estimation (SAKE) (Shin et al., 2014), or Low-rank modeling of local k-space neighborhoods (LORAKS) (Haldar, 2014) outperformed conventional CS algorithms. However, these k-space interpolation methods require either parallel acquisition using multiple coil or the image support to be finite.

This limitations have been overcome by the recently proposed algorithm, called annihilating filter-based low-rank Hankel matrix (ALOHA) (Jin et al., 2016; Ye et al., 2017). ALOHA exploits the fact that the redundancy in the spatial domain results in the low-rank Hankel structured matrix in the weighted k-space domain. This implies that an artifact-free k-space data can be recovered from the incomplete or distorted k-space measurements by constructing Hankel structured matrix and exploiting its low-rankness (Jin et al., 2016). It clearly suggests an answer toward the current limitation in conventional k-space interpolation approaches by uncovering the link between CS and k-space interpolation approaches. ALOHA was shown to be superior to the conventional CS and parallel-MRI (p-MRI) (Jin et al., 2016), and it has been actively applied to several MR-related problems (Jin et al., 2017; Lee et al., 2016a, b). Accordingly, we adopted this idea to solving the QSM dipole inversion problem.

Figure 1: Reconstruction flow of ALOHA-QSM for dipole inversion. The 3-D phase image is first Fourier transformed into 3D k-space. Then, the Haar wavelet spectrum is used for weighting each slice of k-space, which is used as the measurement term in Eq. 7. Eq. 8. was solved at each slice along the selected direction, which is repeated along the other two remaining directions. For the first loop, truncated dipole kernel was used for data initialization, and for the next loops, results from the previous loop was used. Finally, a correction factor was multiplied for compensation of the systematic underestimation due to the k-space null on the conical surface at the original k-space.

In ALOHA, the artifact-free weighted k-space data can be estimated for susceptibility images by formulating an equation as follows:


where is the weighted k-space of the measured phase image, is the weighted k-space of the estimated susceptibility image, is the regularization parameter, is the element-wise multiplication operation, () denotes a Hankel structured matrix of the weighted susceptibility k-space, and RANK() refers to the matrix rank. The first term of Eq. (7) is the data fidelity term, and the second term is the additional regularization term based on the low-rankness of the Hankel structured matrix. Note that the problem is formulated entirely in the Fourier domain, and the weighting comes from Haar wavelet spectrum.

However, it is difficult to estimate the low-rank matrix by solving Eq. (7) directly, since the rank term is not convex. Hence, as suggested in the original ALOHA paper (Jin et al., 2016; Ye et al., 2017), the matrix factorization approach is used for . Although still non-convex, it can be relatively easily solved as shown below. Specifically, Eq. (7) can be re-written as:


where () is the Hermitian transpose of a matrix. There are three variables in this equation: , , and . Hence, we employed the alternating direction method of multiplier (ADMM) algorithm with the following Lagrangian:


where is the hyper-parameter of ADMM. ADMM update of Eq. (9) is given by


Since the dipole inversion requires a calculation of 3-D k-space data, it is difficult to apply ALOHA to the 3-D data directly. More specifically, it is difficult to store three dimensional Hankel matrix in memory and its matrix factorization is also computationally demanding. While circulant approximation of Hankel matrix is recently proposed to handle this issue (Ongie and Jacob, 2017), we overcame this limitation by applying the 2-D ALOHA algorithm to the 3-D data successively along three different directions () in this study.

Specifically, when we take 2-D inverse Fourier transform on a plane of the 3-D k-space along direction, the central and peripheral slices contain low-frequency and high-frequency information, respectively, along the axial direction, which is in domain. This image can be also sparsified using Haar wavelet transform, indicating the existence of the corresponding annihilating filter. It implies that when we take the Hankel matrix on the weighted k-space data along the direction, it also shows low-rank property. Hence, the 3-D k-space can be split up into multiple 2-D planes and each 2-D plane is processed using the 2-D version of ALOHA-QSM algorithm. The same trick can be applied along the other two directions as well, i.e., or direction. To minimize discrepancy according to the choice of the processing direction, this 2-D processing was applied along each of the three directions successively in this study.

2.3 Overview of the ALOHA-QSM method

Specifically, Fig.1 illustrates the reconstruction flow of the quantitative susceptibility maps using the proposed ALOHA-QSM method. Note that 3-D k-space of phase image is used as input for the algorithm. The proposed algorithm consists of three steps. In the first step, one of the k-space directions () was selected, and each 2-D plane of the 3-D k-space along the selected direction was processed by 2-D ALOHA. Here, the 3-D k-space was initialized with the truncated k-space division (TKD) method, where the threshold value was set to 0.1. To exploit the sparsity of the selected plane, we applied Haar wavelet weighting to the plane (Jin et al., 2016). Then an artifact-suppressed k-space plane was reconstructed from the weighted plane by updating Eq. (10) iteratively, and then, the weighting was removed. A new artifact-suppressed 3-D k-space of susceptibility image was formulated by successively repeating the above procedure for all the planes of the original 3-D k-space along the selected direction. This procedure of formulating the artifact-suppressed 3D k-space was applied along each of the other two directions successively as well.

Due to the null points at the conical surface of the dipole kernel’s k-space, the estimated susceptibility values are often underestimated. To address this problem, we also estimated the correction factor by calculating the slope of a least-squares regression between the original phase image and the phase image obtained by applying the forward model (Marques and Bowtell, 2005) on the reconstructed ALOHA-QSM image. The final QSM image was generated after compensation of the intensity scale with the correction factor.

3 Methods

3.1 Public data from repository

To compare the performance of ALOHA with other methods, the numerical phantom and the in-vivo human brain were downloaded from the repository (http://weill.cornell.edu/mri/pages/qsmreview.html) (Wang and Liu, 2015). For the phantom and in vivo data, the resolution was 0.9375 0.9375 1.5 mm and 0.9375 0.9375 1 mm, respectively, and the matrix size was 256 256 98 and 256 256 146, respectively. For the in vivo human brain, datasets acquired with COSMOS along five different head orientations with respect to the direction of main magnetic field were additionally available and downloaded.

3.2 Experimental data for In-vivo human brain

In order to conduct further detailed quantitative comparison between ALOHA-QSM and other QSM methods, in-vivo human GRE images were obtained on 12-channel head coil on a 3T MRI scanner (MAGNETOM Verio, Siemens Healthcare, Erlangen, Germany). Five healthy volunteers (5 males, 21-25 years old) were scanned using conventional multi-echo GRE sequences at a single head orientation. Sequence parameters were: bandwidth Hz/px for all echoes, TR ms, TE ms, TE ms, number of echoes = 3, matrix size , resolution , partial Fourier along phase and slice encoding directions with factor , flow compensation along all the three directions, slice oversampling , and flip angle = . Total acquisition time was min sec for each subject.

3.3 Data Processing

For human brain data acquired from multi-coils, coil combination was executed first. All the data were processed as follows. Brain mask was extracted from magnitude images, using brain extraction tool (BET) of FSL (Smith, 2002), the phase images from multi-coils were combined by Hermitian inner product (Bernstein et al., 1994), multi-echo phase images were corrected with non-linear frequency map estimation (Liu et al., 2013), Laplacian-based phase unwrapping was employed to avoid abrupt phase jumps, and varied version of Sophisticated Harmonic Artifact Reduction on Phase data (V-SHARP) method was applied for eliminating background phase signals (Li et al., 2014, 2011; Wu et al., 2012).

GeForce GTX 1080 graphic card was used for graphic processor unit (GPU), and i7-6700k was used for central processing unit (CPU). All the codes were written in MATLAB 2014b (Mathwork, Natick). Computed Unified Device Architecture (CUDA) was employed for GPU to accelerate the algorithm.

Optimal hyper-parameters for ALOHA-QSM reconstruction were determined by evaluating root-mean squared error (RMSE) referenced to the true susceptibility map for the numerical phantom and the COSMOS image for the in-vivo human brain. The performance of the ALOHA-QSM method was evaluated in comparison with three representative QSM reconstruction methods: TKD (Shmueli et al., 2009), iLSQR (Li et al., 2011, 2015), and MEDI (Liu et al., 2011, 2012a). Reconstruction parameters were as follows: for TKD, for MEDI, and 30 iteration steps for iLSQR, which were confirmed from previous studies (Li et al., 2015; Wang and Liu, 2015).

Mean susceptibility values of several brain anatomical structures were acquired from the regions of deep brain gray matter structures: substantia nigra, red nucleus, globus pallidus, putamen, and head of caudate nucleus. ROIs of gray matter structures were manually segmented from a slice of brain susceptibility images. The mean susceptibility values of TKD, iLSQR, MEDI, and ALOHA in the structures were compared each other for the five in-vivo human brains. Also, the susceptibility values from ALOHA-QSM was compared with those from literatures (Bilgic et al., 2012; Lim et al., 2013; Liu et al., 2011; Wharton and Bowtell, 2010). In addition, for the data from repository linear regression was used to indicate how well the QSM methods estimate susceptibility as reference to the ground truth provided in the repository, and the root-mean squared error was calculated for the quantitative evaluation.

4 Results

Figure 2: Impact of the reconstruction parameters in Eq. 9. (a,c) From left to right columns: root mean squared error of ALOHA reconstructed susceptibility map referenced to COSMOS image, L2-norm of the data fidelity term (cost1), and the matrix factorization of the rank estimation (cost2) of the reconstructed k-space with ALOHA were evaluated. (b,d) Sagittal slice of the susceptibility map.
Figure 3: Susceptibility maps of the numerical 3-D brain model were reconstructed with various methods. From left to right: Truth, TKD, iLSQR, MEDI, and ALOHA was used for reconstruction. Images are scaled from -0.10 to 0.20 ppm.
Figure 4: (a) Comparison of the human brain susceptibility map determined with COSMOS, TKD, iLSQR, MEDI, and ALOHA. (b) Zoomed image of the axial slice, and (c) susceptibility profile along the red line of (b).

4.1 Parameter determination

From Eq. (9), we can find that the controls weighting of the matrix factorization of the rank estimation and the determines iterative updates on the ADMM algorithm. Fig.2 represents how RMSE changed as a function of hyper-parameters and in various ranges. A data fidelity term for the phase-susceptibility relationship (cost1; ) and the regularization term from the matrix factorization of the rank estimation (cost2; ) were also analyzed for comparison. When was below , RMSE was low and the susceptibility image was reconstructed well (Fig.2(a,b)(\romannum1,\romannum2)). However, when the exceeds , RMSE increased rapidly (Fig.2(a)(\romannum3)), and the image showed severe streaking artifacts (Fig.2(b)(\romannum3)).

Effects of are displayed in Fig.2(c)-(d). Good susceptibility images could be acquired when was adjusted to minimize RMSE () (Fig.2(d)(\romannum2,\romannum3)) . As became smaller than the optimal point, streaking artifacts started to appear in the images due to weighting of the data fidelity term (Fig.2(d)(\romannum1)). As increased above the optimal point, the images were over-smoothed due to weighting of the regularization term from the matrix factorization (Fig.2(d)(\romannum3)).

From these observations, RMSE graph were analyzed with L-curve fitting method, and then we found and as optimal values. Similar procedures were applied for the numerical phantom from the repository, and the in-vivo human brain acquired from five normal volunteers (phantom: , ; human: , ). Note that the parameter values satisfying L-curve fitting criteria of RMSE graph are close to the maximum points in the cost2 graph.

Figure 5: (a) Comparison of the human brain susceptibility map determined with TKD, iLSQR, MEDI, and ALOHA. (b) Estimated susceptibility value of deep gray matter structures from five healthy volunteers were compared for various QSM reconstruction methods.

4.2 Comparison of QSM reconstruction methods

For the numerical phantom, all the methods successfully reconstructed susceptibility images without severe streaking artifacts. Among them, MEDI showed least RMSE values and images with least fluctuation (Fig.3, Table.1).

However, the results were quite different for the in vivo data as shown in Fig.4 and Table.1, where TKD, iLSQR, MEDI, and ALOHA were compared each other and COSMOS was considered as reference. TKD showed severe streaking artifacts. iLSQR reconstructed susceptibility images well, but showed minor artifacts. MEDI eliminated streaking artifacts perfectly, while it showed awkward delineation on the boundaries of brain structures along with artificial over smoothing effects. In addition, structures of gray matter and white matter in cortical areas were hard to distinguish. Overall ALOHA removed streaking artifacts effectively with no over-smoothing effect. ALOHA reduced artifacts better than TKD or iLSQR, although not as perfect as MEDI. However, ALOHA provided visually more natural images, where structural details were delineated well. Fig.4(c) represents the susceptibility profile along the red line in Fig.4(b). Line profiles of the ALOHA and MEDI matched those of the reference (COSMOS) well, while the profiles of TKD and iLSQR were quite different from those of the reference.

Figure 6: Impact of range of hyper parameters for ALOHA reconstruction: (a) and (b) . X-axis is log-scaled value of each parameter. Data attained from repository (top) and gradient echo with single-head orientation (bottom) was used for processing. From left to right columns: error from the data fidelity term and the matrix factorization form of the rank estimation were evaluated with range of parameters. Red point indicates optimal point for ALOHA reconstruction (, ).
Figure 7: Sagittal view of the reconstructed k-space with each QSM reconstruction methods. From the image of TKD, iLSQR, and ALOHA, abnormal k-space energy distribution along magic angle are shown.

For the data acquired at a single head orientation on the five human subjects (Fig.5(a)), the results were overall similar to those for the human data from the repository shown in Fig.4: TKD showed severe streaking artifacts, whereas iLSQR and ALOHA reconstructed the QSM images well. However, some streaking artifacts were observed in MEDI (Fig.5(a)), in contrast to the repository data (Fig.4). For deep gray matter structures, ALOHA showed susceptibility values similar to those of the other QSM methods (Fig.5(b)), while iLSQR slightly underestimated the results. Even there were variations in the susceptibility values for each of QSM reconstruction methods among literature, ALOHA-QSM showed a reasonable range of quantitative estimates (Table.2).

5 Discussion

In this study, we proposed a novel QSM image reconstruction algorithm using the ALOHA. The proposed algorithm successfully reconstructed susceptibility images for the numerical phantom and human brains acquired at a single head orientation.

Parameters for ALOHA-QSM were determined by analyzing RMSE with L-curve fitting criteria. For this procedure, the reference images acquired with COSMOS are necessary to calculate the RMSE of the reconstructed image. Since COSMOS needs to be scanned at multiple head orientations, however, it cannot be applied on a dataset acquired at a single head orientation. In Fig.2, we found that the regularization term from the matrix factorization of the rank estimation in the objective function (; ) was maximized at the optimal value of parameters, which satisfies L-curve fitting criteria of the RMSE graph. From this observation, we used the as criteria of the parameter optimization, rather than RMSE. Based on this criteria, we tried to analyze the optimal parameters for human brain data acquired at a single head orientation (Fig.6). All graphs from the five subjects showed similar patterns with difference in scale. As we expected, parameter values that maximized were most appropriate for the image reconstruction (data not shown). Hence, we selected and for optimization, and we confirmed that the ALOHA-QSM could be utilized for any phase measurements acquired at a single head orientation.

Figure 8: To overcome underestimating issue in ALOHA reconstructed QSM image, different schemes were evaluated: optimal regularization parameter with correction factor (; ALOHA(scaled)), without correction factor (L1.5), and smaller parameter (; L1.4). (a) Regularization term from the matrix factorization. (b) Susceptibility values in gray matter structures. (c) (Top) Regression plot, (middle) sagittal slice of reconstructed k-space image, (bottom) coronal slice of reconstructed image.
Simulated Brain In Vivo Brain
slope RMSE(%) slope RMSE(%)
TKD 0.89 0.90 1.62 0.34 0.88 4.80
iLSQR 0.94 0.90 1.18 0.55 0.86 3.09
MEDI 0.99 0.99 0.44 0.58 0.93 3.07
(without correction factor)
0.89 0.91 1.63 0.51 0.63 2.79
(with correction factor)
0.89 0.92 1.63 0.51 0.92 3.51
Table 1: Regression coefficients (, slope) and root mean squared error (RMSE(%)) for each reconstruction methods in simulated brain and in-vivo human brain.

width=1 Method N Substantia nigra Red nucleus Globus pallidus Putamen Caudate nucleus ALOHA (current study) 5 0.110.04 0.090.01 0.120.02 0.050.01 0.060.02 COSMOS (Lim et al., 2013) 5 0.090.01 0.060.00 0.110.01 0.030.01 0.020.00 COSMOS (Liu et al., 2011) 9 0.130.03 0.090.02 0.190.02 0.090.04 0.080.02 COSMOS (Wharton and Bowtell, 2010) 5 0.170.02 0.140.02 0.190.02 0.100.01 0.090.02 L2 RSO (Bilgic et al., 2012) 11 0.080.04 0.070.03 0.120.02 0.060.02 0.070.02 L2 RSO (Wharton and Bowtell, 2010) 5 0.160.03 0.130.02 0.190.02 0.090.01 0.090.01 MEDI (Bilgic et al., 2012) 11 0.100.04 0.090.04 0.140.02 0.070.02 0.080.02 MEDI (Liu et al., 2011) 9 0.120.03 0.080.05 0.190.02 0.080.02 0.090.02

Table 2: Comparison of the estimated mean susceptibility (ppm) and standard deviation of deep gray matter structures with literatures.
Figure 9: Comparison of in-vivo human brain susceptibility images reconstructed with different processing order of ALOHA. (a)-(c): reconstruction was implemented along single direction; (a) sagittal direction; (b) coronal direction; (c) axial direction. (d)-(f): ALOHA reconstruction was implemented successively along (d) coronal-axial-sagittal direction; (e) axial-sagittal-coronal direction; (f) sagittal-coronal-axial direction. (g) average of susceptibility images reconstructed with single direction. (h) average of successively reconstructed results. Note that the successive ALOHA reconstruction is effective in reducing image artifacts, including streaking artifact in sagittal view, and is stable to the reconstruction direction compared to the single directional ALOHA reconstruction.

Figure 7 shows a sagittal plane of a reconstructed k-space for each method. For ALOHA-QSM, signals were slightly suppressed near the region of conical surface. The image of TKD and iLSQR also showed similar problems near the conical surface. These unnatural energy level distributions in k-space are caused by insufficient information on the conical surface, leading to systematic underestimation of the susceptibility values. It commonly happens for k-space based QSM reconstruction methods. Many studies were conducted to compensate for the underestimation. In TKD, the intensity of point spreading function constructed from the difference between the original dipole kernel and the corresponding truncated kernel was used as correction factor (Schweser et al., 2013). For post processing, the method enforcing structural consistency on the conical surface with structural prior, referred to as consistency on the cone data (CCD), was suggested (Wen et al., 2016). Alternatively, the underestimation issue could be overcome by using smaller , the optimization parameter for the regularization term from the matrix factorization. In deep gray matters, the estimated susceptibility values with a smaller value of were higher than those estimated with the original value of (Fig.8(b,c)). Also, we can confine the missing k-space information to the values estimated with a smaller value. However, this approach still underestimated the susceptibility values (Fig.8(b)) and the image showed more artifacts (Fig.8(c)). To address this issue, we introduced the correction factor to compensate for the underestimation problem. As shown in Fig.8(b)(c), this provides better correlation with COSMOS data.

Previously, ALOHA has been restricted to application to 2-D datasets due to the memory requirement to store Hankel matrix. In this study, we resolved this issue by applying the ALOHA algorithm three times successively along three axes of the 3-D k-space. Successive implementation contributed to enhance processing speed and quality of the resulting image. When the algorithm was implemented along a single axis, images showed a large variation along the processing direction (Fig.9(a)-(c)). This variation decreased with successive reconstruction along the three directions (Fig.9(d)-(f)), and also image quality was improved with successive reconstruction. Furthermore, although it takes longer processing time, taking average on the result from each directions can further enhance the data consistency and quality of the image (Fig.9(g,h)).

Existing QSM reconstruction algorithms can be classified into two mainstreams: k-space based and magnitude image guided methods. K-space based methods have a fast processing time and provide fine vascular structures, but are prone to streaking artifacts as shown for TKD in Fig.4. These streaking artifacts disturb evaluation of susceptibility images. The methods using magnitude images as prior information for regularization such as MEDI can generate susceptibility maps free of streaking artifacts, but the resulting maps suffer from over-smoothing which can obscure the detailed structures. In addition, magnitude images are often quite different from susceptibility images, which may distort the structural information in susceptibility maps from MEDI. Our suggested algorithm, ALOHA-QSM, reconstructed susceptibility maps with suppression of streaking artifacts while maintaining detailed structures and free of over-smoothing. This method can be classified into the k-space based method; however, ALOHA-QSM is a kind of CS-theory based algorithms. While most CS algorithms reconstruct images in the image domain, ALOHA directly interpolates k-space domain by exploiting the low-rankness of weighted Hankel matrix. It is, therefore, apparently different from the conventional CS approach, and more closely related to the basic principle of susceptibility estimation. This new approach may open up new direction on QSM reconstruction.

Many studies were conducted to measure venous-blood oxygen saturation levels from the magnetic susceptibility property of the hemoglobin in blood using phase imaging or T2*-weighted imaging. However, locations and sizes of susceptibility sources in phase images are different from their actual values because of effects of convolution with the dipole kernel. T2*-weighted imaging is not solely dependent on a tissue iron property. Quantitative susceptibility mapping (QSM) has a potential to resolve this problem. Specifically, ALOHA-QSM can be a good candidate for mapping venous oxygen saturation levels, because both streaking artifacts and over-smoothing effects are minimal. Further research should be conducted to investigate the possibility of venous-blood oxygen level mapping from the susceptibility image generated by ALOHA-QSM.

6 Conclusion

We demonstrated that the annihilating filter-based low-rank Hankel matrix approach (ALOHA) can be successfully used to solve the three-dimensional QSM dipole inversion problem. This new k-space based method reconstructed accurate susceptibility maps with suppression of streaking and smoothing artifacts while requiring no anatomical information. Therefore, it has high potential for QSM image reconstruction.

7 Acknowledgement

This work was supported by the National Research Foundation of Korea (NRF-2016R1A2B3008104, NRF-2017R1A2B2006526) and the Korea Health Technology R&D Project through the Korea Health Industry Development Institute (KHIDI), funded by the Ministry of Health & Welfare of South Korea (HI16C1111).


  • Acosta-Cabronero et al. (2013) Acosta-Cabronero, J., Williams, G. B., Cardenas-Blanco, A., Arnold, R. J., Lupson, V., Nestor, P. J., 2013. In vivo quantitative susceptibility mapping (QSM) in Alzheimer’s disease. PloS one 8 (11), e81093.
  • Bernstein et al. (1994) Bernstein, M. A., Grgic, M., Brosnan, T. J., Pelc, N. J., 1994. Reconstructions of phase contrast, phased array multicoil data. Magnetic resonance in medicine 32 (3), 330–334.
  • Bilgic et al. (2012) Bilgic, B., Pfefferbaum, A., Rohlfing, T., Sullivan, E. V., Adalsteinsson, E., 2012. MRI estimates of brain iron concentration in normal aging using quantitative susceptibility mapping. Neuroimage 59 (3), 2625–2635.
  • De Rochefort et al. (2008) De Rochefort, L., Brown, R., Prince, M. R., Wang, Y., 2008. Quantitative MR susceptibility mapping using piece-wise constant regularized inversion of the magnetic field. Magnetic Resonance in Medicine 60 (4), 1003–1009.
  • Eskreis-Winkler et al. (2017) Eskreis-Winkler, S., Zhang, Y., Zhang, J., Liu, Z., Dimov, A., Gupta, A., Wang, Y., 2017. The clinical utility of QSM: disease diagnosis, medical management, and surgical planning. NMR in Biomedicine 30 (4).
  • Haacke et al. (2010) Haacke, E., Tang, J., Neelavalli, J., Cheng, Y., 2010. Susceptibility mapping as a means to visualize veins and quantify oxygen saturation. Journal of Magnetic Resonance Imaging 32 (3), 663–676.
  • Haacke et al. (2004) Haacke, E. M., Xu, Y., Cheng, Y. N., Reichenbach, J. R., 2004. Susceptibility weighted imaging (SWI). Magnetic resonance in medicine 52 (3), 612–618.
  • Haldar (2014) Haldar, J. P., 2014. Low-rank modeling of local k-space neighborhoods: from phase and support constraints to structured sparsity. In: SPIE Optical Engineering+ Applications. International Society for Optics and Photonics, pp. 959710–959710–12.
  • Jin et al. (2016) Jin, K. H., Lee, D., Ye, J. C., 2016. A general framework for compressed sensing and parallel MRI using annihilating filter based low-rank Hankel matrix. IEEE Transactions on Computational Imaging 2 (4), 480–495.
  • Jin et al. (2017) Jin, K. H., Um, J.-Y., Lee, D., Lee, J., Park, S.-H., Ye, J. C., 2017. MRI artifact correction using sparse+ low-rank decomposition of annihilating filter-based Hankel matrix. Magnetic resonance in medicine 78 (1), 327–340.
  • Jin and Ye (2015) Jin, K. H., Ye, J. C., 2015. Annihilating filter-based low-rank hankel matrix approach for image inpainting. IEEE Transactions on Image Processing 24 (11), 3498–3511.
  • Langkammer et al. (2017) Langkammer, C., Schweser, F., Shmueli, K., Kames, C., Li, X., Guo, L., Milovic, C., Kim, J., Wei, H., Bredies, K., 2017. Quantitative susceptibility mapping: Report from the 2016 reconstruction challenge. Magnetic Resonance in Medicine.
  • Lee et al. (2016a) Lee, D., Jin, K. H., Kim, E. Y., Park, S., Ye, J. C., 2016a. Acceleration of MR parameter mapping using annihilating filter‐based low rank Hankel matrix (ALOHA). Magnetic resonance in medicine 76 (6), 1848–1864.
  • Lee et al. (2016b) Lee, J., Jin, K. H., Ye, J. C., 2016b. Reference-free single-pass EPI Nyquist ghost correction using annihilating filter-based low rank Hankel matrix (ALOHA). Magnetic resonance in medicine 76 (6), 1775–1789.
  • Li et al. (2014) Li, W., Avram, A. V., Wu, B., Xiao, X., Liu, C., 2014. Integrated Laplacian-based phase unwrapping and background phase removal for quantitative susceptibility mapping. NMR in Biomedicine 27 (2), 219–227.
  • Li et al. (2015) Li, W., Wang, N., Yu, F., Han, H., Cao, W., Romero, R., Tantiwongkosi, B., Duong, T. Q., Liu, C., 2015. A method for estimating and removing streaking artifacts in quantitative susceptibility mapping. Neuroimage 108, 111–122.
  • Li et al. (2011) Li, W., Wu, B., Liu, C., 2011. Quantitative susceptibility mapping of human brain reflects spatial variation in tissue composition. Neuroimage 55 (4), 1645–1656.
  • Lim et al. (2013) Lim, I. A. L., Faria, A. V., Li, X., Hsu, J. T., Airan, R. D., Mori, S., van Zijl, P. C., 2013. Human brain atlas for automated region of interest selection in quantitative susceptibility mapping: application to determine iron content in deep gray matter structures. Neuroimage 82, 449–469.
  • Liu et al. (2012a) Liu, J., Liu, T., de Rochefort, L., Ledoux, J., Khalidov, I., Chen, W., Tsiouris, A. J., Wisnieff, C., Spincemaille, P., Prince, M. R., 2012a. Morphology enabled dipole inversion for quantitative susceptibility mapping using structural consistency between the magnitude image and the susceptibility map. Neuroimage 59 (3), 2560–2568.
  • Liu et al. (2011) Liu, T., Liu, J., De Rochefort, L., Spincemaille, P., Khalidov, I., Ledoux, J. R., Wang, Y., 2011. Morphology enabled dipole inversion (MEDI) from a single-angle acquisition: comparison with cosmos in human brain imaging. Magnetic resonance in medicine 66 (3), 777–783.
  • Liu et al. (2009) Liu, T., Spincemaille, P., De Rochefort, L., Kressler, B., Wang, Y., 2009. Calculation of susceptibility through multiple orientation sampling (COSMOS): a method for conditioning the inverse problem from measured magnetic field map to susceptibility source image in MRI. Magnetic Resonance in Medicine 61 (1), 196–204.
  • Liu et al. (2013) Liu, T., Wisnieff, C., Lou, M., Chen, W., Spincemaille, P., Wang, Y., 2013. Nonlinear formulation of the magnetic field to source relationship for robust quantitative susceptibility mapping. Magnetic resonance in medicine 69 (2), 467–476.
  • Liu et al. (2012b) Liu, T., Xu, W., Spincemaille, P., Avestimehr, A. S., Wang, Y., 2012b. Accuracy of the morphology enabled dipole inversion (MEDI) algorithm for quantitative susceptibility mapping in MRI. IEEE transactions on medical imaging 31 (3), 816–824.
  • Lustig et al. (2008) Lustig, M., Donoho, D. L., Santos, J. M., Pauly, J. M., 2008. Compressed sensing MRI. IEEE signal processing magazine 25 (2), 72–82.
  • Marques and Bowtell (2005) Marques, J., Bowtell, R., 2005. Application of a fourier-based method for rapid calculation of field inhomogeneity due to spatial variation of magnetic susceptibility. Concepts in Magnetic Resonance Part B: Magnetic Resonance Engineering 25 (1), 65–78.
  • Ongie and Jacob (2017) Ongie, G., Jacob, M., 2017. A fast algorithm for convolutional structured low-rank matrix recovery. IEEE Transactions on Computational Imaging 3 (4), 535–550.
  • Salomir et al. (2003) Salomir, R., de Senneville, B. D., Moonen, C. T., 2003. A fast calculation method for magnetic field inhomogeneity due to an arbitrary distribution of bulk susceptibility. Concepts in Magnetic Resonance Part B: Magnetic Resonance Engineering 19 (1), 26–34.
  • Schweser et al. (2013) Schweser, F., Deistung, A., Sommer, K., Reichenbach, J. R., 2013. Toward online reconstruction of quantitative susceptibility maps: superfast dipole inversion. Magnetic resonance in medicine 69 (6), 1581–1593.
  • Schweser et al. (2012) Schweser, F., Sommer, K., Deistung, A., Reichenbach, J. R., 2012. Quantitative susceptibility mapping for investigating subtle susceptibility variations in the human brain. Neuroimage 62 (3), 2083–2100.
  • Shin et al. (2014) Shin, P. J., Larson, P. E., Ohliger, M. A., Elad, M., Pauly, J. M., Vigneron, D. B., Lustig, M., 2014. Calibrationless parallel imaging reconstruction based on structured low-rank matrix completion. Magnetic resonance in medicine 72 (4), 959–970.
  • Shmueli et al. (2009) Shmueli, K., de Zwart, J. A., van Gelderen, P., Li, T., Dodd, S. J., Duyn, J. H., 2009. Magnetic susceptibility mapping of brain tissue in vivo using MRI phase data. Magnetic resonance in medicine 62 (6), 1510–1522.
  • Smith (2002) Smith, S. M., 2002. Fast robust automated brain extraction. Human brain mapping 17 (3), 143–155.
  • Wang and Liu (2015) Wang, Y., Liu, T., 2015. Quantitative susceptibility mapping (QSM): decoding MRI data for a tissue magnetic biomarker. Magnetic resonance in medicine 73 (1), 82–101.
  • Wen et al. (2016) Wen, Y., Wang, Y., Liu, T., 2016. Enhancing k-space quantitative susceptibility mapping by enforcing consistency on the cone data (CCD) with structural priors. Magnetic resonance in medicine 75 (2), 823–830.
  • Wharton and Bowtell (2010) Wharton, S., Bowtell, R., 2010. Whole-brain susceptibility mapping at high field: a comparison of multiple-and single-orientation methods. Neuroimage 53 (2), 515–525.
  • Wu et al. (2012) Wu, B., Li, W., Guidon, A., Liu, C., 2012. Whole brain susceptibility mapping using compressed sensing. Magnetic resonance in medicine 67 (1), 137–147.
  • Ye et al. (2017) Ye, J. C., Kim, J. M., Jin, K. H., Lee, K., 2017. Compressive sampling using annihilating filter-based low-rank interpolation. IEEE Transactions on Information Theory 63 (2), 777–801.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description