Reference-Guided Sparsifying Transform Design for Compressive
Sensing MRIS. Derin Babacan, Xi Peng, Xian-Pei Wang, Minh N. Do, Zhi-Pei Liang
Abstract— Compressive sensing (CS) MRI aims to accurately
reconstruct images from undersampled k-space data. Most
CS methods employ analytical sparsifying transforms suchas total-variation and wavelets to model the unknown image
and constrain the solution space during reconstruction. Recently, nonparametric dictionary-based methods for CS-MRI
reconstruction have shown significant improvements over theclassical methods. These existing techniques focus on learning
the representation basis for the unknown image for a synthesisbased reconstruction. In this paper, we present a new framework for analysis-based reconstruction, where the sparsifying
transform is learnt from a reference image to capture theanatomical structure of unknown image, and is used to guide thereconstruction process. We demonstrate with experimental data
the high performance of the proposed approach over traditional
methods.

I. INTRODUCTION
The emerging theory of compressive sensing (CS) has
created significant interest in magnetic resonance imaging(MRI), as it offers a robust and systematic framework forincreasing the imaging speed which is inherently low in MRI.
This slow acquisition is prohibitive in many important clinical applications, especially in dynamic and interventional
imaging, where near real-time acquisition speeds are desired.
Despite the advances in imaging hardware, the acquisition
speed is limited by MR physics and prohibitive in many
applications.
Compressive sensing provides a systematic methodology
to recover images from a significantly smaller number ofmeasurements than that required by Nyquist sampling. The
underlying principle to make this possible is the sparsityof the images in some transform domain. An extensive
body of work investigated a number of different sparsifying
transforms that exploit the image characteristics for accurate
reconstruction. Classical approaches consider non-adaptive,
analytical transforms such as wavelets, total-variation, andcontourlets [1], [2]. More recently, adaptive transforms (ordictionaries) have become increasingly popular as they canmore accurately capture characteristics of the images ofinterest in a particular application. Early work in this area
focused on analytical derivation of parametric bases [3],whereas more recently, nonparametric dictionary learning has
been developed and applied to many problems in imaging
S. D. Babacan, M. N. Do and Z.-P. Liang are with the Beckman Institutefor Advanced Science and Technology and the Department of Electrical andComputer Engineering, University of Illinois at Urbana-Champaign, Urbana,IL, USA. e-mail: dbabacan@illinois.edu
X. Peng and X.-P. Wang are with the School of Electronic Information,
Wuhan University, Wuhan, China. X. Peng is now a visiting scholar atthe Beckman Institute for Advanced Science and Technology, University ofIllinois at Urbana-Champaign.

[4]–[7]. The work in [8] proposed reference-guided weights
to analytical transforms during reconstruction.
Most of the work in dictionary-based recovery utilizes
learning a dictionary from a collection of precollected reference images. The trained dictionary is then used as thesparse representation basis for the image, and which isthen employed in a synthesis-based reconstruction algorithm
typically minimizing the lp -norm (0 ≤ p ≤ 1) of therepresentation coefficients. The work [9] also proposed joint
dictionary learning and sparse reconstruction which showedsignificant improvements over existing methods.
The main principle behind these methods is to constrain
the characteristics of the reconstructed image to certain
anatomical structures, rather than blindly and uniformly
enforcing local smoothness as in traditional methods. The
morphology of the subject of interest is generally known toan high extent, which provides valuable information to guidethe reconstruction process. In this paper, we propose a novel
algorithm based on the same principle. However, rather thanconstraining the representation basis, an adaptive sparsifying
transform is learnt from a reference image and used in
an analysis-based reconstruction. This transform is essentially spatially-varying high-pass filters which enforce local
anatomical structure (learnt from a reference image) duringreconstruction. We demonstrate with experimental resultsthat the proposed method is very effective for compressive
sensing MRI reconstruction while maintaining robustness.
The rest of this paper is organized as follows. The proposed formulation is described in Section II. Experimental
validation of the approach is done in Section III, followedby concluding remarks in Section IV.
II. PROPOSED FORMULATION
The goal in compressive sensing MRI is to reconstruct theunknown image ρ ∈ Cp from limited k-space measurementsd ∈ Cm with m  p. The imaging system is commonlymodeled asd = Fu ρ + n,(1)mxp

where Fu ∈ C
is the undersampled Fourier encoding
matrix, and n is the observation noise. Typical formulationsof the CS reconstruction problem use either the analysisbased regularization as [10]ρ̂ = argmin k d − Fu ρ k22 + λ R (Ψρ) ,(2)ρor, synthesis-based sparse reconstruction usingα̂ = argmin k d − Fu Φα k22 + λ R (α) ,α(3)and setting ρ̂ = Φα̂. The parameter λ controls the strengthof the regularization. In both cases, the regularization functional R(·) is generally a variant of the lp -norm with 0 ≤p ≤ 1, which is known to enforce sparsity in the solutions. Although the difference between the two approaches
might appear to be subtle, they deviate significantly bothin terms of modeling and algorithm development [10]. The
analysis-based algorithms are common in many inverse problems; algorithms based on total-variation, and x-let (wavelet,curvelet, etc.) forward (or analysing) transforms for Ψ, fall
into this category. The synthesis-based approach is relatively
more recent, and it is based on the assumption that the signal
ρ can be represented as a linear combination of columns ofΦ, and this representation α is sparse.
In both approaches, the transforms Φ and Ψ are chosen
traditionally as data-independent, i.e., analytical transforms
with certain desirable properties such as fast evaluationand inversion. The development of data-dependent, nonparametric dictionaries showed considerable improvement overanalytical ones due to better and more specific modelingof image classes [3]–[7]. As the synthesis-based approachis constructive, and therefore more intuitive for modeling,
all dictionary learning approaches are developed using thisapproach. To our knowledge, dictionary learning based on theanalysis approach is not investigated, although analysis-based
optimization is easier to develop and solve. Our proposal
in this work is to develop such an approach for referenceguided reconstruction of MR images from undersampled kspace measurements.
The most important property in designing the analysing
transform Ψ is its sparsifying ability, i.e., a good transformshould produce highly sparse coefficients such that imposing
R (Ψρ) significantly constrains the solution space of ρ.
The transform Ψ should accurately capture the local spatialactivity of image ρ, and in the ideal case, its application tothe image will result in uncorrelated noise.
Analytical analysis operators such as total variation andwavelets impose sparseness on all preselected filters. For
instance, total-variation suppresses finite differences in bothhorizontal and vertical directions in the whole image domain. However, if prior knowledge on the edge structure isavailable, the suppression of image variation can employ thisinformation and follow the edge structure, which in turn will
result in higher sparsity levels.
In most MRI applications, a reference fully-sampled image
ρref can be obtained to extract the local anatomical structure
within the image. This information can be effectively used toadaptively design filters to impose the learnt local anatomical
structure during the reconstruction process.
Based on this, we define the analysis operator as spatiallyvarying high-pass filters defined at each pixel location r:1
1 X
refδ(s)f (r, s) g ρref, (4)Ψ(r)=δ(r)−r , ρsKrs∈N (r)s6=r

1 Note that the r th row of the operator Ψ defines a filter for a particularlocation r.

(a)

(b)
0.45
0.35
0.4

50

50
0.3

0.35

0.3
100

0.25
100

0.25
0.2
0.2

150

150
0.15

0.15
0.1
200

0.1

200

0.05
0.05
250

250
050

100

150

200

250

50

100(c)

150

200

250(d)

Fig. 1. (a) Original image, and (b) reference image taken from a diffusion
MR experiment. Both images are normalized to unit maximum intensity.
The magnitudes of the coefficients of the output of the filters learnt fromthe reference image are shown in (c), and the magnitudes of horizontal
gradients in the original image are shown in (d).

where δ(r) is the delta function, N (r) defines a rectangularneighborhood around pixel r, and Kr is the normalization
factor. The filter Ψ(r) consists of two parts: A spatiallyinvariant filter f (r, s), used to weigh closer pixels higher thanref,pixels far apart; and a spatially-varying filter g ρrefr , ρscomputed using the local spatial variations
inthereferrefence image. Conceptually, g ρrefmust assign largerr , ρsweights to pixels with similar intensity values, such thatthe filtering is performed along the edges, and filtering
across edges is suppressed. Based on this, we define g(·)
to be a function of the radiometricbetween pixels
 distance
ref
ref
ref,
ρ=g|ρ
r and s, such that g ρrefsr − ρs | . Anyrmonotonically decreasing function can be used for g(·); in
this work, we use the absolute differences such thatrefg ρref=r , ρs1ref|ρrefr − ρs + |,(5)where  is a small number which accounts for the observationnoise in the reference image, and for possible anatomical
discrepancies between the reference and unknown images.
As an example of how much structural information can beobtained using the proposed analysis operator, consider theimage pair in Fig. 1, extracted from a dynamic MR sequence.
We have constructed the filter Ψ using the reference shownin Fig. 1(b) and applied it to the image in Fig. 1(a). The
reference and original images have different local and globalcontrast, and they have different structure around the sculp.
The magnitudes of the coefficients at the filter output isshown in Fig. 1(c). For comparison, the absolute values of
the horizontal differences are shown in Fig. 1(d). Noticethat the proposed filter produces a much sparser output, andcaptures the characteristics of the image to a greater extent.

ρ̂ = argmin k d − Fu ρ k22 + λ k Ψρ k1 ,(6)

ρwhich can be solved using a number of existing methods.
Here we use the reweighted least squares method as a reference method to compare the proposed design with traditionalalgorithms.

30
P−OriginalP−Noisy
P−ReferenceTV
l1−Noisy

25

Relative Error (%)

In addition, the filtered image have noise-like characteristics,
and hence much of the anatomical structure is captured by
the proposed filter, whereas the horizontal gradients exhibit
significant image structure, as expected.
Using the reference-guided analysis transform describedabove, the reconstruction of the image ρ can be performed
using the classical optimization shown in (2). Employing thestandard l1 -minimization formulation, (2) becomes

l1−Reference

20151050510

III. RESULTS AND DISCUSSION
A. Simulation Results
The proposed method has been evaluated systematicallyusing different reference images with varying degrees of similarity to the original image. To demonstrate the advantage ofthe method, we include comparisons with two methods; thetraditional TV-based CS reconstruction [1], and an l1 -basedalgorithm that solvesρ̂ = argmin k d − Fu ρ k22 + λ k ρ − ρref k1 .

(7)

ρThis method uses the reference image directly in the reconstruction as an additional constraint. In all algorithms,
the algorithmic parameters are optimized in each case as torepresent their best reconstruction results.
The unknown target image is chosen from a time series dynamic MRI experiment, and is shown in Fig. 1(a). A variable
density sampling pattern is used to acquire undersampled
k-space data at various undersampling ratios. We test theproposed approach with three different reference images:
the target image, which is used to provide an “oracle”
performance bound; the target image corrupted with additive
Gaussian noise of variance 0.01; and another frame fromthe same dynamic MR sequence, shown in Fig. 1(b). The
noisy target image as the reference is used to assess therobustness of the proposed method to noise in the reference.
As mentioned above, the reference image in Fig. 1(b) contains the general anatomical structure, but it has different
local intensity variations and observation noise. Therefore, it
represents a realistic reference image that can be acquired inmany applications.
As the quantitative quality measure we use the relative
2error, defined by kρ̂−ρkkρk2 , where ρ is the original image, andρ̂ is the reconstructed image, respectively. Figure 2 shows
the performance of the algorithms with different reference
images and different undersampling ratios. The proposed
algorithm is denoted with preceding P-, and the formulationin (7) is denoted with preceding l1 −, with different reference
images. Reconstructed images at undersampling factors of
15% and 25% are shown in Fig. 3 and Fig. 4, respectively.
Several important remarks can be made from the results.
First, the proposed approach provides better reconstruction

152025Sampling Percentage

3035Fig. 2.
Relative reconstruction error vs. sampling percentage for all
algorithms with different reference images.

performance at all undersampling factors and with all reference images. Second, while the performance of the approach
degrades when reference images with lower similarity areused, this degradation is not severe (around 2-3%), and themethod provides performance close to the oracle bound.
It can be stated that if the reference image contains thegeneral anatomical structure (albeit degraded), the proposed
approach is able to exploit this information to an highextent during reconstruction. On the contrary, the l1 -based
reconstruction using the reference image is highly prone
to the noise in the reference image, and the relative reconstruction error is mostly determined by the noise in thereference image. In addition, when an image with different
contrast characteristics is used, its performance decreases
significantly and becomes close to the TV reconstruction.
Finally, the TV-reconstruction is not able to provide accuratereconstructions; the reconstructed images either exhibit large
aliasing artifacts or oversmoothing.
B. Discussion
The proposed method can be employed for a varietyof imaging applications where reference images can becollected that contain the anatomical structure. Potential suchapplications include interventional MRI, diffusion-weighted
MRI, and MR spectroscopy imaging. To more effectively
utilize the proposed approach, a variety of experiments (e.g.,with different contrast weightings or diffusion tensor imaging) can be performed to acquire additional information forthe design of the analysis transform (see a related discussionin [8]).
In cases where the expected similarity between the target
and reference images can be estimated, this knowledge canbe incorporated in the transform design using the filter
parameters (such as  in (5)). The parameter can be adjusted
to obtain more robustness to discrepancies between the target
and reference images. In addition, similar to the method in[9], the filters learnt from the reference image can be adaptedusing the reconstructed images during the iterative procedure.

(a)

(b)

(a)

(b)

(c)(d)(c)

(d)

(e)

(e)

Fig. 3. Reconstructed images at 15% sampling ratio. Top row: Proposed
method, middle row: l1 -based method using (7), bottom row: TV reconstruction. (a,c) is obtained with the noisy target image as the reference, and(b,d) is obtained with the image in Fig. 1(b) as the reference.

Fig. 4. Reconstructed images at 25% sampling ratio. Top row: Proposed
method, middle row: l1 -based method using (7), bottom row: TV reconstruction. (a,c) is obtained with the noisy target image as the reference, and(b,d) is obtained with the image in Fig. 1(b) as the reference.

These approaches will be considered in future work. Finally,
although the focus of this work is compressive sensing
reconstruction, the proposed construction of the analysis
operator can also be used for other recovery applications,such as denoising and motion artifact removal.

[2] Y. Kim, M. I. Altbach, T. P. Trouard, and A. Bilgin, “Compressedsensing using dual-tree complex wavelet transforms,” in Proc. ISMRM
2009, Honolulu, HI, April 2009.
[3] Z.-P. Liang and P. Lauterbur, “An efficient method for dynamic
magnetic resonance imaging,” IEEE Trans. Med. Imaging, vol. 13,
pp. 677–686, 1994.
[4] M. Aharon, M. Elad, and A. Bruckstein, “K-SVD: An Algorithm
for Designing Overcomplete Dictionaries for Sparse Representation,”
IEEE Trans. on Signal Processing, vol. 54, no. 11, pp. 4311–4322,
2006.
[5] A. Bilgin, Y. Kim, F. Liu, and M. S. Nadar, “Dictionary design forcompressed sensing MRI,” in Proc. of ISMRM, Stockholm, Sweden,May 2010.
[6] R. Otazo and D. K. Sodickson, “Adaptive compressed sensing MRI,”in Proc. of ISMRM 2010, 2010, p. 4867.
[7] Y. Chen and X. Ye, “A novel method and fast algorithm for MR
image reconstruction with significantly under-sampled data,” Inverse
Problems and Imaging, vol. 4, no. 2, pp. 223 – 240, May 2010.
[8] J. P. Haldar, D. Hernando, S.-K. Song, and Z.-P. Liang, “Anatomicallyconstrained reconstruction from noisy data,” Magn. Reson. Med.,
vol. 59, no. 4, pp. 810–818, 2008.
[9] S. Ravishankar and Y. Bresler, “MR image reconstruction from highly
undersampled k-space data by dictionary learning,” IEEE Trans. on
Med. Imaging, vol. PP, no. 99, p. 1, 2010.
[10] M. Elad, P. Milanfar, and R. Rubinstein, “Analysis versus synthesis
in signal priors,” Inverse Problems, vol. 23, no. 3, pp. 947–968, June
2007.

IV. CONCLUSIONSWe presented a novel method for compressive sensing
MRI based on the design of reference-guided analysis transforms. We have demonstrated that this formulation achieves
a high degree of sparsity and is able to capture and enforce
the salient anatomical structure during the reconstruction
process. The proposed method has shown to provide higher
reconstruction performance compared to methods using analytical transforms and to methods using the reference images
directly as constraints.
R EFERENCES[1] M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The applicationof compressed sensing for rapid mr imaging,” Magn. Reson. Med.,
vol. 58, no. 6, pp. 1182–1195, Dec 2007.