Quantitative Susceptibility Map Reconstruction Using Annihilating Filterbased LowRank Hankel Matrix Approach
Abstract
Quantitative susceptibility mapping (QSM) inevitably suffers from streaking artifacts caused by zeros on the conical surface of the dipole kernel in kspace. This work proposes a novel and accurate QSM reconstruction method based on a direct kspace interpolation approach, avoiding problems of over smoothing and streaking artifacts. Inspired by the recent theory of annihilating filterbased lowrank Hankel matrix approach (ALOHA), QSM reconstruction problem is formulated as deconvolutional problem under lowrank Hankel matrix constraint in the kspace. To reduce the computational complexity and the memory requirement, the problem is formulated as successive reconstruction of 2D planes along three independent axes of the 3D phase image in Fourier domain. Extensive experiments were performed to verify and compare the proposed method with existing QSM reconstruction methods. The proposed ALOHAQSM effectively reduced streaking artifacts and accurately estimated susceptibility values in deep gray matter structures, compared to the existing QSM methods. Our suggested ALOHAQSM algorithm successfully solves the threedimensional QSM dipole inversion problem without additional anatomical information or prior assumption and provides good image quality and quantitative accuracy.
keywords:
Quantitative susceptibility mapping, Dipole inversion, Lowrank Hankel matrix completion1 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 neurodegenerative disease (AcostaCabronero et al., 2013) and neurosurgical planning (EskreisWinkler 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 gradientecho imaging protocols and offer quantitative values of the tissue susceptibility directly with no additional scan.
From the physical viewpoint, every susceptibility source builds dipoleshaped 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 illposed. The illposed 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 kspace following the magic angle (the conical surface) with tuning parameters (threshold based kspace division; TKD (Shmueli et al., 2009; Schweser et al., 2013)), using intrinsic information of gradient echo phase images for homogeneity reconstruction (homogeneityenabled 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, oversmoothing, and/or magnitude image induced errors (Wang and Liu, 2015).
Recently, annihilating filterbased lowrank 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 downsampled 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 kspace interpolation using the lowrank 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, kspace 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 twodimensional datasets.
In this study, we propose a new kspace based QSM image reconstruction method, referred to as ALOHAQSM, to obtain QSM images with relatively high accuracy and suppression of streaking artifacts by applying ALOHA for 3D phase images. ALOHAQSM 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 dipoleshaped 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:
(1) 
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
(2) 
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 3D, 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 ”kspace” and their axes as in kspace and in the image domain. Hence, the dipole kernel in Eq. (2) can be transformed into the simple kspace filter as follow:
(3) 
where is a kspace vector , and . Also, by denoting , the relationship between the phase and the magnetic susceptibility in Eq. (1) can be simplified into:
(4) 
where and are the kspace 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 kspace 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:
(5) 
However, it causes discontinuity and vanishing surfaces on the resultant kspace, 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:
(6) 
where is the weighting mask derived from the magnitude image, is the inverse Fouriertransform of the matrix, is the regularization parameter, is a threedimensional gradient operator, and is a binary mask of the magnitude edge information. Since MEDI uses fullysampled 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 kspace regions. By iteratively subtracting this artifactoriginated image, artifactreduced susceptibility image is obtained.
These aforementioned algorithms are reported to estimate susceptibility values well (Wang and Liu, 2015). However, they still have many limitations. Kspace 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 downsampled kspace 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 25.
An important requirement for applying CS algorithm on MR image reconstruction is that the artifacts originated from the kspace 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 underdetermined regions of the dipole inversion filter in Fourier domain () is clearly form a conical shape. Fortunately, its undersampling 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, kspace direct interpolation methods determine missing kspace data using intrinsic structural information of the kspace. Kspace interpolation methods including simultaneous autocalibration and kspace estimation (SAKE) (Shin et al., 2014), or Lowrank modeling of local kspace neighborhoods (LORAKS) (Haldar, 2014) outperformed conventional CS algorithms. However, these kspace 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 filterbased lowrank 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 lowrank Hankel structured matrix in the weighted kspace domain. This implies that an artifactfree kspace data can be recovered from the incomplete or distorted kspace measurements by constructing Hankel structured matrix and exploiting its lowrankness (Jin et al., 2016). It clearly suggests an answer toward the current limitation in conventional kspace interpolation approaches by uncovering the link between CS and kspace interpolation approaches. ALOHA was shown to be superior to the conventional CS and parallelMRI (pMRI) (Jin et al., 2016), and it has been actively applied to several MRrelated problems (Jin et al., 2017; Lee et al., 2016a, b). Accordingly, we adopted this idea to solving the QSM dipole inversion problem.
In ALOHA, the artifactfree weighted kspace data can be estimated for susceptibility images by formulating an equation as follows:
(7) 
where is the weighted kspace of the measured phase image, is the weighted kspace of the estimated susceptibility image, is the regularization parameter, is the elementwise multiplication operation, () denotes a Hankel structured matrix of the weighted susceptibility kspace, 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 lowrankness 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 lowrank 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 nonconvex, it can be relatively easily solved as shown below. Specifically, Eq. (7) can be rewritten as:
(8) 
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:
(9) 
where is the hyperparameter of ADMM. ADMM update of Eq. (9) is given by
(10) 
Since the dipole inversion requires a calculation of 3D kspace data, it is difficult to apply ALOHA to the 3D 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 2D ALOHA algorithm to the 3D data successively along three different directions () in this study.
Specifically, when we take 2D inverse Fourier transform on a plane of the 3D kspace along direction, the central and peripheral slices contain lowfrequency and highfrequency 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 kspace data along the direction, it also shows lowrank property. Hence, the 3D kspace can be split up into multiple 2D planes and each 2D plane is processed using the 2D version of ALOHAQSM 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 2D processing was applied along each of the three directions successively in this study.
2.3 Overview of the ALOHAQSM method
Specifically, Fig.1 illustrates the reconstruction flow of the quantitative susceptibility maps using the proposed ALOHAQSM method. Note that 3D kspace of phase image is used as input for the algorithm. The proposed algorithm consists of three steps. In the first step, one of the kspace directions () was selected, and each 2D plane of the 3D kspace along the selected direction was processed by 2D ALOHA. Here, the 3D kspace was initialized with the truncated kspace 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 artifactsuppressed kspace plane was reconstructed from the weighted plane by updating Eq. (10) iteratively, and then, the weighting was removed. A new artifactsuppressed 3D kspace of susceptibility image was formulated by successively repeating the above procedure for all the planes of the original 3D kspace along the selected direction. This procedure of formulating the artifactsuppressed 3D kspace 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 kspace, the estimated susceptibility values are often underestimated. To address this problem, we also estimated the correction factor by calculating the slope of a leastsquares regression between the original phase image and the phase image obtained by applying the forward model (Marques and Bowtell, 2005) on the reconstructed ALOHAQSM 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 invivo 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 Invivo human brain
In order to conduct further detailed quantitative comparison between ALOHAQSM and other QSM methods, invivo human GRE images were obtained on 12channel head coil on a 3T MRI scanner (MAGNETOM Verio, Siemens Healthcare, Erlangen, Germany). Five healthy volunteers (5 males, 2125 years old) were scanned using conventional multiecho 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 multicoils, 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 multicoils were combined by Hermitian inner product (Bernstein et al., 1994), multiecho phase images were corrected with nonlinear frequency map estimation (Liu et al., 2013), Laplacianbased phase unwrapping was employed to avoid abrupt phase jumps, and varied version of Sophisticated Harmonic Artifact Reduction on Phase data (VSHARP) 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 i76700k 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 hyperparameters for ALOHAQSM reconstruction were determined by evaluating rootmean squared error (RMSE) referenced to the true susceptibility map for the numerical phantom and the COSMOS image for the invivo human brain. The performance of the ALOHAQSM 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 invivo human brains. Also, the susceptibility values from ALOHAQSM 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 rootmean squared error was calculated for the quantitative evaluation.
4 Results
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 hyperparameters and in various ranges. A data fidelity term for the phasesusceptibility 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 oversmoothed due to weighting of the regularization term from the matrix factorization (Fig.2(d)(\romannum3)).
From these observations, RMSE graph were analyzed with Lcurve fitting method, and then we found and as optimal values. Similar procedures were applied for the numerical phantom from the repository, and the invivo human brain acquired from five normal volunteers (phantom: , ; human: , ). Note that the parameter values satisfying Lcurve fitting criteria of RMSE graph are close to the maximum points in the cost2 graph.
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 oversmoothing 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.
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, ALOHAQSM 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 ALOHAQSM were determined by analyzing RMSE with Lcurve 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 Lcurve 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 ALOHAQSM could be utilized for any phase measurements acquired at a single head orientation.
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  

0.89  0.91  1.63  0.51  0.63  2.79  

0.89  0.92  1.63  0.51  0.92  3.51 
Figure 7 shows a sagittal plane of a reconstructed kspace for each method. For ALOHAQSM, 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 kspace are caused by insufficient information on the conical surface, leading to systematic underestimation of the susceptibility values. It commonly happens for kspace 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 kspace 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 2D 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 3D kspace. 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: kspace based and magnitude image guided methods. Kspace 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 oversmoothing 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, ALOHAQSM, reconstructed susceptibility maps with suppression of streaking artifacts while maintaining detailed structures and free of oversmoothing. This method can be classified into the kspace based method; however, ALOHAQSM is a kind of CStheory based algorithms. While most CS algorithms reconstruct images in the image domain, ALOHA directly interpolates kspace domain by exploiting the lowrankness 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 venousblood 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, ALOHAQSM can be a good candidate for mapping venous oxygen saturation levels, because both streaking artifacts and oversmoothing effects are minimal. Further research should be conducted to investigate the possibility of venousblood oxygen level mapping from the susceptibility image generated by ALOHAQSM.
6 Conclusion
We demonstrated that the annihilating filterbased lowrank Hankel matrix approach (ALOHA) can be successfully used to solve the threedimensional QSM dipole inversion problem. This new kspace 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 (NRF2016R1A2B3008104, NRF2017R1A2B2006526) 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).
References
 AcostaCabronero et al. (2013) AcostaCabronero, J., Williams, G. B., CardenasBlanco, 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 piecewise constant regularized inversion of the magnetic field. Magnetic Resonance in Medicine 60 (4), 1003–1009.
 EskreisWinkler et al. (2017) EskreisWinkler, 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. Lowrank modeling of local kspace 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 lowrank 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+ lowrank decomposition of annihilating filterbased Hankel matrix. Magnetic resonance in medicine 78 (1), 327–340.
 Jin and Ye (2015) Jin, K. H., Ye, J. C., 2015. Annihilating filterbased lowrank 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. Referencefree singlepass EPI Nyquist ghost correction using annihilating filterbased 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 Laplacianbased 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 singleangle 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 fourierbased 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 lowrank 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 lowrank 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 kspace 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. Wholebrain susceptibility mapping at high field: a comparison of multipleand singleorientation 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 filterbased lowrank interpolation. IEEE Transactions on Information Theory 63 (2), 777–801.