Compressed Sensing of Spatial Electron Paramagnetic Resonance Imaging

David H. Johnson, Rizwan Ahmad, Guanglong He, Alexandre Samouilov, Jay L. Zweier

This is an author's electronic pre-print of the following publication. MRM allows pre-prints to be shared through individual author's websites; for terms, please see this website's license. Please cite the published work as:

Compressed Sensing of Spatial Electron Paramagnetic Resonance Imaging. Johnson DH, Ahmad R, He G, Samouilov A, Zweier JL. Magn Reson Med 2013. Pubmed: not assigned yet. Journal's version.

Address for correspondence: Jay L. Zweier, MD. Jay.Zweier@osumc.edu.

ABSTRACT

Purpose: To improve image quality and reduce data requirements for spatial electron paramagnetic resonance imaging (EPRI) by developing a novel reconstruction approach using compressed sensing (CS).

Methods: EPRI is posed as an optimization problem, which is solved using regularized least-squares with sparsity promoting penalty terms, consisting of the l1 norms of the image itself and the Total Variation of the image. Pseudo-random sampling was employed to facilitate recovery of the sparse signal. The reconstruction was compared to the traditional Filtered Back-Projection reconstruction for simulations, phantoms, isolated rat hearts, and mouse GI tracts labeled with paramagnetic probes.

Results: A combination of pseudo-random sampling and CS was able to generate high-fidelity EPR images at high acceleration rates. For 3D phantom imaging, CS-based EPRI showed little visual degradation at 9-fold acceleration. In rat heart datasets, CS-based EPRI produced high quality images with 8-fold acceleration. A high resolution mouse GI tract reconstruction demonstrated a visual improvement in spatial resolution and a doubling in SNR.

Conclusion: A novel 3D EPRI reconstruction utilizing compressed sensing was developed and offers superior SNR and reduced artifacts from highly undersampled data.

INTRODUCTION

Continuous wave spatial electron paramagnetic resonance (EPR) imaging is a form of tomography where a 3D object is recovered from projections blurred by a known point spread function (i.e., the spatially invariant electron spectrum) (1,2). EPR has the potential to measure in vivo free radicals, which provide important physiological information (e.g., oxygen, pH, and cellular redox status) (3-6). Spatial EPRI can be used for studies of probe uptake, clearance, metabolism, and biodistribution. Purely spatial EPR imaging can determine the rate of clearance of a probe, such as a nitroxide, which is useful for characterizing tumors redox state. Dynamic changes in signal intensity are important for the study of cancer metabolism or myocardial ischemia and reperfusion (7-13). The absence of EPR probe signal or absence of metabolism may itself be informative; for example, in myocardial ischemia spatial EPRI allows visualization of the risk region and region of necrosis (8). This purely spatial EPRI experiment is sufficient to determine where the loss of viable myocardium has occurred.

Historically, applications of spatial EPRI have been limited by a lack of projections, low spatial resolution, and limited SNR. Recent advances have improved the quality of EPR images by expediting the acquisition using rapid scan or spinning gradient, choosing more effective sets of projections to acquire, and developing narrow linewidth probes with high spin densities (14-16). Yet, the image reconstruction for in vivo experiments is still based on Filtered Back-Projection (FBP) or Algebraic Reconstruction Techniques (ART), both of which require a relatively large number of projections. Reconstructions based on the maximum entropy method and Tikhonov regularization recently demonstrated improvements to spectral-spatial EPR reconstructions in phantom experiments (17-19). Another recent approach is to restrict the EPR probe to multiple discrete sites with non-overlapping spatial locations and perform a lineshape-fitting based reconstruction (20,21). There is a potential to further improve spatial in vivo EPR imaging with diffuse probes by exploiting the inherent sparsity in the images, and this is the main focus of this study.

Compressed Sensing (CS) has been used to reconstruct images from only a small amount of acquired data, yielding a large reduction in acquisition time (22-24). For example, MRI angiography has been accelerated by a factor of 10, which enables imaging of the entire heart in a single breath-hold acquisition (25). The requirements for CS are that the image is sparse in some domain (e.g., finite differences between adjacent pixels, wavelets, etc.), the noise in the image is incoherent, and a nonlinear reconstruction is used to recover the image. We propose that a CS reconstruction can produce high quality EPR images from a fraction of the number of projections which would normally be required. CS has been applied to spectral-spatial EPRI before, but only in the context of multisite oximetry with known lineshape (21,26,27). The CS methods presented in this work are less restrictive as they do not strictly rely on the spatially sparse distributions encountered in multisite oximetry. Another limitation of the multisite CS approach is the requirement to uniquely identify the parameters of the lineshape model, which can be difficult in mixtures of Lorentzian and Gaussian shapes with both dispersion and absorption components.

The goal of this study was to reduce the number of projections required for 3D spatial EPRI by the development and implementation of a reconstruction which uses l1 penalty terms in the optimization and also explicitly incorporates the measured EPR spectrum into the signal model. An exact adjoint operator is defined and used to recover deblurred images via iterative estimation. Sparsifying transformations for EPRI are proposed, and the effects of incoherent sampling are investigated by comparing different gradient angle distributions. The CS reconstruction is tested on a variety of typical EPR imaging experiments using common probes and compared to the FBP reconstruction.

THEORY

Problem Formulation

For spatial EPRI, the forward model is given by Eq (1):

y = C R x

(1)

where the 3D image volume (x) is projected by the Radon transform (R) and those projections are blurred by the convolution operator (C) using the known spectrum to produce the observed projections (y). In the classic FBP reconstruction used in EPRI, image is recovered by applying multiple disjoint steps

x = R^{-1} C^{-1} y \approx R^{*} \frac{d^{2}}{{ds}^{2}} W C^{-1} y

(2)

where the second derivative with respect to projection coordinate (s) and projection weighting (W) operators are needed along with the deconvolution operator (C-1) and the backprojection operator, R*, which is adjoint of R. The weighting operator is necessary when the distribution of projections is not uniform (14). Deconvolution, C-1y, can be performed in Fourier domain using point-by-point division of the Fourier transform of y with the Fourier transform of the spectrum. This implementation, however, amplifies noise and may lead to division-by-zero problems (1). Therefore, a lowpass filter generally precedes the point-by-point division in the Fourier domain, yielding

C^{-1} y \approx F^{-1} (\frac{V (F y)}{F z})

(3)

where F denotes the discrete Fourier transform, z indicates the EPR spectrum, and V represents a lowpass filter in the Fourier domain.

The FBP reconstruction is non-iterative and computationally fast. However, the limitations of a finite number of projections and involvement of lowpass filter V often result in a loss of spatial resolution, streaking artifacts, and SNR degradation. Nevertheless, it remains the mainstay of EPR imaging experiments because it is simple and fast.

We propose a novel EPRI reconstruction based on searching for the solution (xCS), which is not only consistent with the acquired data but only prefers piece-wise constancy and spatial sparsity. This can be achieved by introducing l1 regularization terms to the data fidelity term of the cost function. This approach, which exploits compressibility of EPR images, allows robust reconstruction from a small number of projections. Deconvolution is entirely avoided by incorporating the convolution into the signal model. Instead, the reconstruction directly uses the measured first derivative EPR projections without needing integration, thresholding, or other denoising procedures. Therefore, we solve the optimization problem given by Eq (4):

x_{CS} = \arg \min_{x} {∥C R x - y∥}_{2} + λ_{1} {∥x∥}_{1} + λ_{2} {∥TV (x)∥}_{1}

(4)

where TV stands for the Total Variation. The TV norm favors piece-wise constant images, and it is frequently applied to image reconstructions of MRI and x-ray computed tomography data (25). The norms given above are defined as:

{∥x∥}_{p} = {(\sum_{i = 1}^{N} {|x_{i}|}^{p})}^{\frac{1}{p}}

(5)

where p is the order of the norm (e.g., p=1 for l1) and N is the number of elements in the vector x. The first term with the l2 norm ensures that the solution is consistent with the acquired projections. The weights λ₁ and λ₂ reflect the tradeoffs between favoring solutions which are mostly zeros and more strongly piece-wise constant, respectively. For the degenerate case of λ1 = λ2 = 0, the reconstruction is equivalent to the Algebraic Reconstruction Technique with lineshape information (28).

CS-based EPRI Reconstruction

To solve the optimization problem posed by Eq (4), we introduce a temporary image (u) and use an alternating minimization algorithm to solve the two sub-problems:

\arg \min_{x} {∥C R x - y∥}_{2} + λ_{1} {∥x∥}_{1}

(6)

\arg \min_{u} {∥x - u∥}_{2} + λ_{2} {∥TV (u)∥}_{1}

(7)

The solution to Eq (6) fits the observed projections and also favors images which are mostly zeros. Solving Eq (7) performs TV denoising of the current estimate of the image. Both problems are solved using the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) (29). The intermediate result (z) of TV denoising is used to re-initialize the solution of Eq (6), which reinforces the data consistency and l1 norm. The reconstruction alternates between solving Eq (6) and (7) until the solution converges, which occurs when the l1 norm of the current estimate changes less than a fixed threshold in two successive iterations (e.g., 0.1%).

FISTA is used to solve Eq (6) using a Landweber update to the current estimate of the image (z) followed by the shrinkage operator (S):

x_{n} = S (z_{n} - μ (R^{*} C^{*} C R z_{n} - R^{*} C^{*} y), \frac{λ_{1} μ}{2})

(8)

S (x, λ) = sign (x) \max (|x| - λ, 0)

(9)

where S performs soft-thresholding on each voxel in the image, sending it to zero if it is smaller than the weight λ. Since the spectrum is purely real, C* corresponds to convolution with reversed spectrum, which can be readily implemented using point-by-point multiplication in the Fourier domain. Likewise, C*C is also block circulant and can be implemented by a multiplication in the Fourier domain. It should be noted that the adjoint reconstruction of the projections (R*C*y) does not change, so it is computed only once for speed. The parameter μ is related to the maximum eigenvalue of the forward operator (CR), but its value is assigned manually since it depends on the gradient angles and it is not generally known in advance.

After updating the image with soft-thresholding, a non-negativity constraint is imposed because the spin density cannot assume a negative value. The reconstruction continues with the TV denoising algorithm developed by Taboole et al. (30). We used an isotropic 3D TV computed along the X, Y, and Z axes. An overview of the entire reconstruction is given in Listing 1.

$y_{1} = R^{} C^{} y$ $x_{0} = z_{0} = {[0]}_{N \times N \times N}$ $n = 0$ $a_{0} = 1$ Repeat
	$x_{n} = S (z_{n} - μ (R^{} C^{} C R z_{n} - R^{} C^{} y), \frac{λ_{1} μ}{2})$ $x_{n} = \max (x_{n}, 0)$ $a_{n + 1} = \frac{1 + \sqrt{1 + a_{n}^{2}}}{2}$ $u_{n} = x_{n} + \frac{a_{n} - 1}{a_{n + 1}} (x_{n} - x_{n - 1})$ $z_{n + 1} = \arg \min_{z} {∥z - u_{n}∥}_{2} + λ_{2} {∥TV (z)∥}_{1}$ $n = n + 1$

Until z_n satisfies a stopping criterion.

x_{CS} = z_{n}

Listing 1. Pseudo-code Overview of the CS EPRI Reconstruction

Implementation of a 3D Radon operator and its Exact Adjoint

The 3D Radon operator and its adjoint were constructed using implicit operators rather than explicit matrix representations. The size of the matrices required to explicitly represent these operators as matrices would have far exceeded computational memory for even modestly-sized 64x64x64 image volumes. To handle this limitation, we adopted a memory efficient implementation of R and R* as described by De Man and Basu (31,32). We extended their algorithm to project 3D image volumes to a line using the distance-driven projection and backprojection operators. The distance-driven projection and backprojection operators have the advantages of being fast, accurate, and exact adjoints. Entire subvolumes of the 3D image volume can be projected or backprojected in parallel, and the cache coherence of the algorithm scales well on modern multi-core CPUs. The distance-driven kernel operator was implemented in C for speed, and the remaining functions in the reconstruction were implemented in Matlab for simplicity. The L2 of the projections were normalized to unity before reconstruction to narrow the range of λ1 and λ2 between datasets which had diverse spectral amplitudes.

METHODS

Simulations

Simulations were performed to measure sampling incoherence and sparsity in EPRI. While EPRI cannot achieve truly incoherent sampling because data are collected in the form of projections, there may still be differences in incoherence from choosing different gradient angles. The maximum amplitude of the sidelobe-to-peak ratio in the point spread function (PSF) was used as a measurement of the sampling incoherence (22). A noiseless simulation of the PSF of R was computed by projecting and backprojecting an impulse function in the image domain (Figure 1). The simulation was performed for 100 projections from both the Equal Linear Angle (ELA) sampling strategy and the 3D Golden Means (GM) sampling strategy (33).

Figure 1. Point spread functions (PSF) along X, Y, and Z axes for different sampling strategies (a-c). The GM PSF approaches zero rapidly along every axis, whereas the sidelobes of ELA are bigger and wider. The maximum of the sidelobe-to-peak ratio for GM and ELA are 0.218 and 0.241, respectively, indicating that GM sampling has more incoherence than ELA sampling. Individual 2D slices of the 3D PSFs (d) show ELA's PSF in the Z plane has strong sidelobes.

The sparsity of EPR images under Eq (4) was examined using a synthetic phantom similar to physical phantoms used in previous EPRI publications (1). A noiseless 3D simulation of ellipsoids in the shape of the letters "EPR" was fed through a simulated acquisition using the forward model (RC), and then the CS reconstruction was performed to examine the effects of ELA and GM sampling (Figure 2). The simulation was performed using a first derivative Lorenztian lineshape generated by the Robinson-Mailer model (34) with a HWHM linewidth of 0.01 mT, 0.04 mT modulation amplitude, 2 mT sweep width, 500 projections from either ELA or GM, 25 mT/m gradient amplitude, and 30 mm field of view. The CS reconstruction was performed with λ1 = 0.01, λ2 = 0.005, and μ = 0.4. The algorithm was terminated when there was at most a 0.1% change in the l1-norm of the image for two successive iterations. The normalized mean-squared error (NMSE) was computed relative to the original unblurred phantom.

Figure 2. Compressed Sensing (CS) reconstructions of a 3D phantom with different sampling strategies. Projections (N=625) were simulated with EPR first-derivative Lorentzian lineshape blurring to investigate sampling incoherence and sparsity in EPRI. The normalized mean-squared error (NMSE) is indicated in the above each slice. The spatial resolution and reconstruction accuracy is improved in GM as compared to ELA sampling due to the greater sampling incoherence in GM sampling.

Phantom Experiments

To empirically determine how much acceleration can be applied to EPR acquisitions with real noise, a phantom consisting of 19 tubes (43 mm long, 4.0 mm outer diameter, 3.81 mm inner diameter) of 1 mM tetrathiatriarylmethyl radical (TAM) were imaged on a 1.2 GHz EPR imaging system. The approximate linewidth was 0.01 mT, and the acquisition parameters were 0.012 mT modulation amplitude, 100 kHz modulation frequency, 3 mT sweep width, 43.3 mT center field, 5.24 s per projection, 1,764 projections, 26.7 mT/m gradient amplitude, as previously reported (35). The dataset was retrospectively downsampled by reducing the number of projections while maintaining an approximately uniform distribution. Images were reconstructed with both CS and FBP for comparison of the two reconstruction methods.

To quantify the accuracy of the reconstructions, the sizes and positions of the tubes were measured in the full dataset by least-squares fitting of the center image slice using 2D Gaussian functions. The average of the Full-Width, Half-Maximum (FWHM) of each Gaussian was compared between the CS and FBP reconstructions. The physical dimensions of the tubes were measured using a dial caliper with 0.02 mm accuracy. The physical inner diameter of the tube was compared to the FWHM of the Gaussians, and the physical outer diameter of the tube was compared to the spacing in-between the centers of the Gaussians.

Isolated Rat Heart Imaging

Ex vivo rat heart experiments were performed on the same imaging system to determine the effectiveness of the CS reconstruction with biologically limited SNR. Lithium octa-n-butoxynaphthalocyanine (LiNc-BuO) suspension was infused into an isolated rat heart, and then the heart was subjected to global ischemia. The linewidth of the probe narrowed to 0.01 mT before the image acquisition began (36). Projections (N=256) from ELA were acquired in 697 s with 0.02 mT modulation amplitude, 1.0 mT sweep width, 40 mT/m gradient amplitude, and 25 mm field of view. To examine the effects of different sparsifying transforms, the CS reconstruction was performed with different reconstruction weights. The dataset was reconstructed with λ1 = λ2 = 0 (equivalent to ART), and then reconstructed again with λ1 = 2 and λ2 = 0.02 (both penalties). In both cases, the images were reconstructed on an 80x80x80 volume with μ = 0.5.

A second isolated rat heart was perfused with LiNc-BuO, and 1,024 projections were acquired with GM sampling in 266 s. A linewidth of 0.01 mT was measured. The imaging acquisition consisted of 0.04 mT modulation amplitude, 1.4 mT sweep width, 47 mT/m gradient amplitude, and 30 mm field of view. Subsets of the projections were reconstructed to determine the minimum amount of data needed to reconstruct a high quality image. The acceleration factor was defined as the fraction of projections retained from the original 1,024 projections (e.g., an acceleration factor of 4 corresponded to using only the first 256 projections to perform the reconstruction). FBP reconstructions were performed on the same subsets. The peak SNR was calculated in each reconstructed image by dividing the maximum amplitude of the image by the standard deviation of the non-signal regions of the image. The non-signal regions were determined by thresholding the FBP reconstruction from the full dataset, and this mask was applied to all subsequent calculations. The CS reconstruction was performed with λ1 = 0.5, λ2 = 0.005, and μ = 1.

High Resolution In Vivo Mouse GI Tract Imaging

To examine the effects of the CS reconstruction on spatial resolution, a high resolution gastro-intestinal imaging experiment was performed by feeding mice a mixture of activated charcoal and glucose, as previously described (37). Projections (N=1,024) from ELA were acquired in 91 min with a linewidth of 0.1 mT, 0.14 mT modulation amplitude, 6.0 mT sweep width, 150 mT/m gradient amplitude, and 40 mm field of view. The CS reconstruction was performed with λ1 = 0.5, λ2 = 0.02, and μ = 1.

Image Sharpness Index

To quantify the sharpness of edges in the GI tract and heart images, the image sharpness index method was employed (38). A sigmoid function (L) given by

L (x) = I_{0} + \frac{I_{1}}{1 + 10^{- w (x - x_{0})}}

(10)

was fitted along the normalized position (x) in 1-dimensional line profiles orthogonal to image edges using a least-squares fit. The fitted parameters of the model consisted of an offset intensity (I0), an intensity amplitude (I1), a center offset (x0), and a sigmoid width (w). Larger values of w corresponded to sharper edges. The mean and standard deviation of the w in all edges in each image were computed, reported as image sharpness index (ISI), and compared using a Student t-test.

RESULTS

Simulations

Incoherent sampling in EPRI was examined by comparing the point spread function (PSF) resulting from the ELA and GM sampling strategies (Figure 1). Individual line profiles along the X, Y, and Z axes show that the GM sampling is more incoherent than ELA sampling due to the smaller sidelobes. The maximum of the sidelobe-to-peak ratio for GM and ELA are 0.218 and 0.241, respectively, indicating that GM sampling has more incoherence than ELA sampling. The PSF of GM sampling also goes to zero much more rapidly than ELA sampling. Both distributions have the same center peak width, which is dominated by the linewidth of the Lorentzian distribution used in the simulation.

The effects of incoherent sampling were demonstrated in a noiseless 3D simulation of ellipsoids in the shape of the letters "EPR" (Figure 2). As the phantom was entirely synthetic and noiseless, the artifacts observed in the reconstructions are entirely due to differences in the ELA and GM sampling strategies which serves as a proxy for incoherence. The same CS reconstruction parameters were used for both simulations. Slices from the ELA reconstruction showed strong blurring and streaking along the X and Y axes. Artifacts from the GM reconstruction were smaller and isotropic. The NMSE was approximately 25% lower in the GM reconstruction as compared to the ELA reconstruction.

Phantom Experiments

FBP and CS reconstructions were compared in the 19 tube TAM phantom (Figure 3). The original 1,764 projections were retrospectively accelerated by the acceleration factors of 2, 4, 9, and 18 corresponding to 882, 441, 196, and 98 projections, respectively. The CS images show little visual degradation at 9-fold acceleration, whereas the FBP images have decreasing SNR and increasing streaking artifacts due to angular undersampling.

Figure 3. FBP and CS reconstructions of the 19 tube TAM phantom. The original 1,764 projections were retrospectively accelerated by the factors indicated above to 882, 441, 196, and 98 projections, respectively. The CS images show little degradation at 9-fold acceleration, whereas the FBP images have rapidly decreasing SNR and increasing streaking artifacts due to angular undersampling.

After fitting the center image slice with 2D Gaussian functions, the FWHM of the CS images was 3.89 ± 0.06 mm (mean ± standard deviation) with a spacing of 5.98 ± 0.09 mm. The FWHM of the FBP images was 4.20 ± 0.07 mm with a spacing of 5.95 ± 0.10 mm. The caliper measurements indicated a 3.81 mm inner diameter and a 5.97 mm outer diameter. The spacings between Gaussians of both CS and FBP reconstructions matched the physical outer diameter within one standard deviation, but the FWHM of the FBP Gaussian was bigger than the physical inner diameter by 0.39 mm (10%). In contrast, the CS FWHM was within 0.09 mm of the true inner diameter (2%).

Isolated Rat Heart Imaging

The effects of regularization on the CS reconstruction were demonstrated in a LiPc rat heart dataset (Figure 4). In the absence of any regularization, some streaking artifacts and high frequency noise appeared in the reconstructed images (4-a). When using both the l1 and TV penalty terms, the image SNR increased from 11.6 ± 3.8 to 21.8 ± 7.7 and background artifacts were decreased (4-b).

Figure 4. Effects of regularization using a LiNc-BuO rat heart dataset. In the absence of any regularization, the reconstruction is equivalent to ART with lineshape information (a). When using both the l 1 and TV penalty terms, the image has increased SNR and decreased artifacts (b).

Acceleration of the acquisition was compared by retrospectively reconstructing only a small fraction of the projections originally acquired in a LiNc-BuO rat heart dataset (Figure 5). The FBP reconstruction (top row) degraded as the acceleration factor increased from 2 (i.e., using 512 out of the 1,024 original projections) to 32 (using only 32 projections). In contrast, the CS reconstruction (bottom row) also degraded, but the 8-fold acceleration was visually similar to the fully sampled image. The ISI for the FBP reconstruction from the full 1,024 projections was 0.56±0.02 (mean ± standard deviation), which was not different (P>0.25) than any of the other FBP reconstructions up to an acceleration factor of 16X. In contrast, the ISI for the CS reconstruction with the full 1,024 projections was 0.88±0.25, which was significantly different (P<0.01) than the FBP reconstructions. There was no significant difference in ISI between CS reconstructions with an acceleration factor of 16X or less. We repeated the CS reconstructions after dividing both reconstruction weights λ1 and λ2 by 10, and the conclusions did not change.

Figure 5. Comparison of Compressed Sensing (CS) and Filtered Backprojection (FBP) reconstructions in a LiNc-BuO heart dataset. As the acceleration factor increases, the visual quality of both algorithms degrade. Both algorithms fail at 32-fold acceleration.

The effects of acceleration on SNR were examined in the LiNc-BuO rat heart dataset on a log-log plot (Figure 6). The predicted square root of time reduction in SNR was similar in both the FBP and CS reconstructions, though the CS reconstruction started with higher SNR and ISI. When the CS reconstruction was used, the SNR dropped from 70.9 (fully sampled) to 11.6 (8-fold acceleration). In comparison, the FBP reconstruction dropped from 12.9 to 2.8 (fully sampled vs. 16-fold).

Figure 6. SNR in 3D reconstructed heart images vs. acceleration factor on a log-log scale. Increasing the acceleration factor degrades both reconstructions. Least squares fits of SNR vs. acceleration factor indicate a slightly faster degradation of the CS reconstruction, but the CS reconstruction starts at a higher SNR.

In Vivo Mouse GI Tract Imaging

The effects of deconvolution were examined by comparing the FBP and CS reconstructions of activated charcoal in the GI tract of a mouse (Figure 7). The spatial resolution of the FBP reconstruction is limited by explicit deconvolution, which results in blurring and negative values in the image. The resolution in the CS reconstruction is limited only by the original projections and the linewidth of the char. The SNR improved from 10.6 ± 3.1 to 20.7 ± 6.1 in the FBP and CS reconstructions, respectively. The ISI for the FBP reconstruction was 0.56±0.05, which was significantly different (P<0.001) than the 1.13±0.66 ISI for the CS reconstruction.

Figure 7. EPR imaging of GI tract from a mouse fed charcoal EPR probe. CS improves the spatial resolution because convolution with the EPR spectrum is included in the forward model, whereas FBP must use explicit deconvolution. The boundaries of the GI tract are sharpened by the l 1 prior. The SNR was doubled. Arrows indicate “spike” artifacts removed by the CS reconstruction.

DISCUSSION

A novel CS reconstruction algorithm for 3D spatial EPR imaging was developed and evaluated in simulations, phantoms, and murine imaging models of cardiac and intestinal anatomy. Both penalty terms were necessary to sparsify EPR images and produce robust, high SNR images given the constraints of typical EPRI linewidths and undersampling. The incoherence of sampling strategies was characterized by the PSF and by simulations. The FBP reconstruction was inferior to the CS reconstruction in terms of SNR. An 8-fold acceleration was possible in the cardiac images with no visual loss of quality, and SNR was doubled in gastro-intestinal tract images. Background signals were very effectively suppressed in the animal experiments. Image sharpness index was approximately doubled in both cardiac and gastro-intestinal EPR images by the CS reconstruction as compared to the FBP reconstruction.

Compressed sensing has not been applied to EPRI before in the form proposed here. An earlier study applied parametric modeling to reconstruct 1D spatial, 1D spectral profiles, but 3D reconstructions were deemed too large to solve (17). The use of the TV operator is an advantage unexplored by previous applications of CS to EPRI. The TV operator removes oscillations and other background signals from the reconstructed images, which makes the CS reconstruction less sensitive to noise. This effect is visually apparent in the GI tract images where “spike” artifacts in the FBP reconstructions are removed in the CS reconstructions (Fig 7, arrows). Previous reconstruction algorithms have tried to improve the spatial resolution of EPRI by improving the deconvolution algorithm or by incorporating the lineshape information into the reconstruction with varying degrees of success (39). In contrast, this method avoids explicit deconvolution entirely. Instead, the lineshape information is seamlessly incorporated directly into the forward model, which avoids choosing a cut-off frequency that reduces SNR for improved spatial resolution. Signals that extend beyond the field segment collected (i.e., due to hyperfine splitting, or intentionally truncated to reduce acquisition time) can also be incorporated into the forward model. Specifically, the operators C and C* are implemented by the same functions when hyperfine splitting occurs in the EPR spectrum. The phantom experiments indicated that the spatial resolution is better in the CS reconstruction than the FBP reconstruction when the physical dimensions of the tubes are assumed to be known.

There are several complimentary approaches available to expediting the EPR acquisition. If only a small region of the anatomy contains a EPR signal of interest, then the magnetic field sweep width can be tailored to cover only that region, resulting in a corresponding reduction in acquisition time (40,41). If a dynamic EPR signal is being investigated for probe kinetics, principal components analysis can be applied to the temporal signal evolution to suppress noise (42). These approaches can be combined with the CS method developed here to exploit the signal sparsity and further accelerate the EPR acquisition. Time domain EPR is already taking advantage of CS reconstructions to great effect (43,44).

The choice of gradient angles is important (45). For static EPR imaging with the CS reconstruction, GM sampling is recommended over other sampling strategies. Traditional sampling strategies like ELA can be used with the CS reconstruction, but a loss of efficiency occurs, and spatial resolution is compromised even in the absence of noise. The GM datasets in this work were much easier to accelerate (i.e., retrospectively downsample) because subsets of projections maintain the same approximate uniformity due to the recursive nature of the golden ratio.

There are some limitations of the new reconstruction. Choosing the reconstruction weights λ1 and λ2 is difficult to do systematically. L curve methods can be used for picking a single reconstruction weight (46), but their generalization to multiple weights remains an open research topic (47). A complete and rigorous evaluation of different choices of these values is outside the scope of this preliminary study. Normalizing the l2 of the projections to unity made the values of λ1 and λ2 more consistent between datasets. Another consideration is the additional computational time required for the reconstruction, which was approximately six minutes for the CS images shown in Figure 7 (bottom row) but less than one minute for FBP (top row). Computational power continues to improve, and the distance-driven Radon transform can be implemented using graphics card computation for additional speed (31). A further reduction in reconstruction time can be achieved using the Fourier transform to compute the Gram operators (R*C*CR) in Eq (8) with Toeplitz embedding (48).

Probes with multiple peaks in the EPR spectrum (e.g., nitroxides) can be used in the same forward model developed here. Some spatial resolution will be lost to the broader linewidths, but it may be possible to reconstruct images from even fewer projections. Accelerating the EPR acquisition further is possible, especially for EPRI experiments that repeatedly acquire the same image to capture the kinetics of redox metabolism. The method proposed here can also be extended to spectral-spatial EPRI by introducing a sparsifying transform for the EPR spectrum, enabling rapid measurements of pH, oxygen, etc. Instead of acquiring individual projections with fast-scan, the spinning gradient acquisition can provide nearly limitless projections very quickly. Recent work has shown how to make the sampling strategy in the spinning gradient acquisition both uniform and efficient (36). Future work could combine a CS reconstruction with the uniform spinning gradient acquisition to produce fast, high quality spectral-spatial images.

CONCLUSION

Compressed sensing has been applied to accelerate and improve spatial EPRI with diffuse probes for the first time. This approach exploits the structure present in EPR images via l1 regularized least-squares. Leveraging efficient implementation of forward operators and their exact adjoints as well as fast optimization routines allows reconstructing large EPR images in a reasonable time. Both simulation and in vivo data indicate that CS-based reconstruction allows reconstructing high fidelity EPR images from highly undersampled data. The resulting images also exhibit high spatial resolution and are devoid of artifacts.