Fast reconstruction for multichannel compressed sensing using a hierarchically semiseparable solver

Abstract

Purpose: The adoption of multichannel compressed sensing (CS) for clinical magnetic resonance imaging (MRI) hinges on the ability to accurately reconstruct images from an undersampled dataset in a reasonable time frame. When CS is combined with SENSE parallel imaging, reconstruction can be computationally intensive. As an alternative to iterative methods that repetitively evaluate a forward CS+SENSE model, we introduce a technique for the fast computation of a compact inverse model solution.

Methods: A recently proposed hierarchically semiseparable (HSS) solver is used to compactly represent the inverse of the CS+SENSE encoding matrix to a high level of accuracy. To investigate the computational efficiency of the proposed HSS-Inverse method, we compare reconstruction time with the current state-of-the-art. In vivo 3T brain data at multiple image contrasts, resolutions, acceleration factors, and number of receive channels were used for this comparison.

Results: The HSS-Inverse method allows for >6× speedup when compared to current state-of-the-art reconstruction methods with the same accuracy. Efficient computational scaling is demonstrated for CS+SENSE with respect to image size. The HSS-Inverse method is also shown to have minimal dependency on the number of parallel imaging channels/acceleration factor.

Conclusions: The proposed HSS-Inverse method is highly efficient and should enable real-time CS reconstruction on standard MRI vendors' computational hardware.

Keywords: SENSE; compressed sensing; hierarchically semiseparable; parallel imaging.

Figure 1

Split Bregman CS+SENSE implementations are illustrated. The required pre-computation is shown above the corresponding flow diagrams. The Matrix Free and Matrix methods rely on iterative CG solutions, while the HSS-Inverse method gives a direct solution for each Split Bregman iteration. The CG based approaches are optimized with the diagonal Jacobi pre-conditioner.

Figure 2

Computational scaling with respect to image size for non-iterative inverse methods. The time for a single inverse evaluation is shown for the optimized sparse Cholesky decomposition CHOLMOD and the HSS solver using a 10⁻⁶ tolerance. The matrix is associated with a R = 3 acceleration and 32 channels. The images range in size from 112 × 112 to 448 × 448.

Figure 3

CS+SENSE reconstructed images and error for T2 and FLAIR imaging contrasts. The dynamic range for the error images is scaled to 1/8 of the fully sampled and sensitivity combined ground truth images. (A) shows the R = 1 sensitivity combined images for the T2 contrast at a resolution of 0.8 × 0.8 × 3mm³. The reconstructed images and error are shown below for R = 3 and 4 accelerations using either the coil compressed 8-channel under-sampled data or the full 32-channel data. Similar results are shown for the FLAIR images at a resolution of 1.0 × 1.0 × 3mm³ in (B).

Figure 4

Relative difference in CS+SENSE reconstructed images for T2 and FLAIR imaging contrasts between the Matrix Free and HSS-Inverse methods. The relative difference is shown for R = 3 and 4 accelerations using either the coil compressed 8-channel under-sampled data or the full 32-channel data.

Figure 5

Computational scaling with respect to image size for CG and HSS based reconstruction methods, see Figure 1 for algorithm flow-diagrams. R = 3 acceleration is applied to the T2 weighted images. A 10⁻⁶ tolerance is assumed for all algorithms to ensure consistent final image error. All methods include 5 iterations of Split Bregman with a TV weighting β = 3 · 10⁻³ and soft-thresholding ε = 2 · 10−¹. The Jacobi pre-conditioner is used for all CG methods. The use of Cartesian optimized coil compression from 32 to 8-channels is explored for the Matrix Free method. The smallest and largest reconstruction times for HSS-Inverse are identified with arrows.

Figure 6

Computational scaling of the HSS-Inverse method with respect to the number of parallel imaging channels and acceleration factor. A 10⁻⁶ tolerance is assumed for 5 iterations of Split Bregman with a TV weighting β = 3 · 10⁻³ and soft-thresholding ε = 2 · 10⁻¹. Cartesian optimized coil compression is used to reduce from 32 to 8-channels. R = 2, 3, and 4 under-sampling is examined.