Systems and methods for fast magnetic resonance image reconstruction using a heirarchically semiseparable solver

Systems and methods for reconstructing images using a hierarchically semiseparable (“HSS”) solver to compactly represent the inverse encoding matrix used in the reconstruction are provided. The reconstruction method includes solving for the actual inverse of the encoding matrix using a direct (i.e., non-iterative) HSS solver. This approach is contrary to conventional reconstruction methods that repetitively evaluate forward models (e.g., compressed sensing or parallel imaging forward models).

BACKGROUND OF THE INVENTION

In clinical applications of structural magnetic resonance imaging (MRI), there exists a multi-objective trade-off between image quality, imaging time, and reconstruction time. Reducing imaging time for a given protocol is clearly beneficial from a cost prospective, and can also facilitate more detailed studies with the same patient throughput. Image quality tends to be a firm barrier placed by radiologists or researchers based upon requirements for data analysis. Finally, stringent hardware limitations exist for clinical FDA approved scanners. It is important to note that advances in MRI sequences and hardware continue to increase the computational burden for image reconstruction, e.g. large coil arrays, increased resolution, and multi-contrast studies. In this work, we investigate a highly scalable inverse algorithm intended to ameliorate the computational challenges associated with accurate compressed sensing (CS) reconstruction.

Sparse signal reconstruction has been introduced for MRI as a method to improve imaging time through random under-sampling of k-space. By assuming a sparsity inducing L1 image prior, the reconstruction problem can be formulated as an unconstrained optimization problem. This problem incorporates fidelity against the observed k-space samples with a penalty imposed on the sparsity prior. These methods have been shown to provide good image accuracy, but can significantly increase the computational burden for image reconstruction. This is especially evident with the inclusion of parallel imaging techniques (e.g., SENSE, GRAPPA, L1-Spirit). Several attempts have been made to reduce the computational requirements associated with sparse signal reconstruction. These iterative techniques rely on repetitive evaluation of forward CS parallel imaging models, however, and thus still have a computational burden.

SUMMARY OF THE INVENTION

The present invention overcomes the aforementioned drawbacks by providing a method for reconstructing an image of a subject from data acquired using a magnetic resonance imaging (MRI) system, in which data acquired from a subject using the MRI system is provided and an inverse of an encoding matrix is computed using a hierarchically semiseparable solver. An image of the subject is then reconstructed from the provided data using the computed inverse of the encoding matrix.

DETAILED DESCRIPTION OF THE INVENTION

Described here are systems and methods for reconstructing images using a hierarchically semiseparable (“HSS”) solver to compactly represent the inverse encoding matrix used in the reconstruction. The reconstruction method described here thus includes solving for the actual inverse of the encoding matrix using a direct (i.e., non-iterative) HSS solver. This approach is contrary to conventional reconstruction methods that repetitively evaluate forward models (e.g., compressed sensing or parallel imaging forward models). The methods described here are capable of achieving a linear computational scaling with respect to the size of the system being solved.

Specifically, the HSS-Inverse described here scales efficiently with the number of imaging voxels and minimizes the influence of the acceleration factor and the number of parallel imaging channels toward the reconstruction time. This technique is therefore capable of achieving upwards of a six-fold speedup over iterative methods, even when those iterative methods take advantage of state-of-the-art pre-conditioning and coil compression techniques.

The adoption of compressed sensing (“CS”) for clinical 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. Rather than solving the resulting linear system repetitively with a conjugate gradient (“CG”) type method, the HSS-Inverse approach described here solves for the actual inverse of the encoding matrix. The proposed HSS inverse model allows for a greater than six-fold speed-up when compared to current state-of-the-art reconstruction methods, as mentioned above, and should enable real-time CS reconstruction on standard MRI vendors' computational hardware.

Referring now toFIG. 1, a flowchart setting forth the steps of an example method for reconstructing images from data acquired with an MRI system using a fast compressed sensing algorithm that implements an HSS solver is illustrated. The method begins by providing data, from which an image is to be reconstructed, as indicated at step102. This data can be provided by retrieving previously acquired data from storage, or may be provided by acquiring the data from a subject using an MRI system.

As indicated at step104, the image reconstruction process includes computing, or otherwise generating, the inverse of an encoding matrix, A, by using an HSS solver to compute a compact representation for the inverse operator, A−1, which is the inverse of the encoding matrix, A. As an example, the inverse operator can be computed by performing a structured factorization of the encoding matrix, A, into lower diagonal, L, and diagonal, D, components,
A=LDLH(1).

Here, many of the terms in the lower diagonal matrix, L, can be easily inverted or can be represented using low-rank modeling. The low-rank properties of the lower diagonal matrix can allow for efficient evaluation of the inverse model.

It is worth noting that the computation of the inverse operator is independent of the acquired, or otherwise provided, data and thus the inverse operator can be pre-computed. By using a pre-computed inverse operator, the total reconstruction time can be significantly reduced. It is also worth noting that the computation of the inverse operator, A−1, has minimal dependency on the number of parallel imaging channels, or on the acceleration factor. It will be appreciated by those skilled in the art that the HSS solver used to compute the inverse operator can be readily made parallel for improved performance.

After the inverse operator has been computed, or otherwise provided, the image reconstruction proceeds, as indicated at step106. In some embodiments, the image reconstruction process is based on a compressed sensing reconstruction. As one example, the image reconstruction can include using Split Bregman (“SB”) reconstructions, which may also incorporate compressed sensing, SENSE, or other, parallel imaging techniques.

CS reconstructions generally involve solving an inverse problem in order to match an observed subset of data under an assumed sparsity prior. As an example, images obtained with magnetic resonance imaging (“MRI”) can be assumed to be sparse or compressible under a total variation (“TV”) transformation, a wavelet transformation, or both.

Generally speaking, the CS formulation for MRI is an unconstrained optimization problem that involves penalty terms based upon assumed TV and wavelet sparsity. By pre-defining penalty weights α and γ, an example CS optimization can estimate the true image, xϵN, as,

where, FΩϵM×Nis the undersampled Fourier operator that transforms the image, x, into k-space to match the observations, yϵM; ψTis the wavelet transform that is applied to the image, x; and TV (x) is the total variation norm of the image, x. In this example, the acceleration factor is R=N/M.

In Eqn. (2), the data fidelity is measured using the L2metric to represent RMSE against the observations, and the L1metric is used to promote sparsity in the wavelet transform domain. The TV operator computes a finite difference across the image, x, to promote sparsity in this spatial smoothness domain. In the examples provided below, a more general parallel imaging problem is focused on and, for ease of illustration, only TV sparsity is considered. It will be appreciated by those skilled in the art, however, that other sparsity can be considered, such as sparsity in the wavelet transform domain, as noted above.

By introducing complex coil sensitivity profiles {Ci}i=1,κ, the SENSE parallel imaging model can be incorporated into the CS formulation as follows:

In Eqn. (3), the TV operator has been re-written as a sum of horizontal and vertical finite difference operators, Ghand Gv. The SB approach relaxes the L1penalties through the iterative construction of L2targets:

By way of example, the targets, gvand gh, can be updated using a soft thresholding truncation parameter, ε. For example, the following update can be implemented:

This operation is linear-time and, thus, the computational cost is dependent on the quadratic minimization shown in Eqn. (4). The explicit solution of this minimization problem is,

where TFis the Fourier operator that has been combined with the Laplacian operator, TS. The inverse problem can thus be referred to as the solution of,

or as evaluating,
x=A−1b(8).

It is important to note that TF, TS, and the parameter, β, are constant with respect to the SB iteration and only depend on the protocol and coil sensitivity maps.

The main computational burden for SB is the repeated solution of Eqn. (6). This computation occurs during updates to the gradients targets, ghand gv, which change b. The constant forward model is represented by the encoding matrix, A. Often times, Eqn. (6) is solved using a Jacobi pre-conditioned conjugate gradient method; however, a compact representation for the SB model inverse, A−1, can be computed using an HSS solver, as described above with respect to step104, and this inverse operator, A−1, can be used to efficiently solve Eqn. (6).

In the following examples, the performance of the HSS-Inverse method was compared against several CG-based approaches to highlight the computational trade-offs for reconstruction. The computational scaling for all approaches was analyzed with respect to the image size and the number of parallel imaging channels. The image accuracy for all methods was computed as RMSE against the complex coil combined images from the fully sampled data. Therefore, the results include computational aspects for the algorithms and analysis of the methods using the acquired data. Exhaustive sweeps of TV and soft-thresholding parameters were performed for “best case” accuracy. As both the CG and HSS solvers have controllable accuracy, a typical 10−6criteria was used for all methods to ensure consistent results across the reconstructions.

In order to accurately compare the different CS approaches, multi-contrast in vivo data were acquired from a healthy volunteer subject to institutionally approved protocol consent. The data were acquired on a 3T MRI system with a 32-channel head array coil. T2-weighted and Fluid Attenuated Inversion Recovery (“FLAIR”) images were acquired with a 224×224 mm2FOV, across 35 slices with a 30% distance factor. The T2-weighted scan uses a Turbo Spin Echo (“TSE”) sequence with imaging parameters TR=6.1 s, TE=98 ms, flip angle=150 degrees, and a resolution of 0.5×0.5×3.0 mm3, with a matrix size of 448×448.

The FLAIR scan used a TSE sequence with an inversion pulse and imaging parameters TR=9.0 s, TE=90 ms, TI=2.5 s, flip angle=150 degrees, and a resolution of 0.9×0.9×3.0 mm3, with a matrix size of 256×256. The fully sampled uncombined complex k-space data were retrospectively under-sampled for all computational experiments. In order to examine the computational scaling of the CS reconstruction algorithms, data sets of consistent size were generated across the multiple imaging contrasts. Where applicable, matrix sizes of 112×112, 168×168, 224×224, and 280×280 were constructed by down-sampling the coil data. These images were utilized to represent resolutions of 0.8, 1.0, 1.33, and 2.0 mm for the same FOV.

FIG. 2illustrates an example of the computational scaling of several Split Bregman optimization techniques with respect to image size. The Matrix Free and Matrix methods rely on pre-conditioned CG to solve Eqn. (6), and the proposed HSS-Inverse method uses the HSS direct solver described above. To ensure consistent reconstruction error, all numerical approaches assumed a 10−6tolerance for the solution. The times reported inFIG. 2correspond to five iterations of Split Bregman with a TV weighting β=3·10−3and soft-thresholding ε=2·10−1. The Jacobi pre-conditioner is used for all the CG methods. The use of Cartesian optimized coil compression from 32 to 8-channels is explored for the Matrix Free method. The HSS-Inverse method had times of 1.1 s and 5.4 s for in-plane resolutions of 2×2 mm2to 0:8×0:8 mm2. With the use of 4× channel compression the Matrix Free method became the best performing alternative to HSS-Inverse.

FIG. 3illustrates the lack of dependence of the HSS-Inverse method on the number of parallel imaging channels and acceleration factor. Here, reconstruction parameters consistent with the results shown inFIG. 2were used. R=2, 3, and 4 accelerations were examined across in-plane resolutions from 2×2 mm2to 0:8×0:8 mm2. The deviation in reconstruction time for HSS-Inverse was under 0.7 s for all cases considered. This small deviation in time can be considered a constant based upon the numerical conditioning of the matrix, A.

It is important to note that the Matrix Free methods compared inFIGS. 2 and 3utilize highly optimized FFT code and by re-implementing the MATLAB code used to implement the HSS-Inverse method inFIGS. 2 and 3into a lower level programming language, it is expected that further improvements in the speed-up times can be achieved.

Thus, an efficient CS reconstruction strategy for MRI assuming SENSE parallel imaging has been described. The proposed HSS-Inverse method exploits the fact that the Split Bregman framework produces a series of least squares problems with a fixed reconstruction operator.

The proposed HSS-Inverse method demonstrates minimal computational dependency with respect to both the acceleration factor and the number of parallel imaging channels. Given the dependence of the CG-based Matrix Free method on channel count, it is contemplated that the speed-up of the HSS-Inverse method will increase when considering larger array coils as additional compressed channels will be required for similar accuracy.

As was alluded to above, the efficiency of HSS-Inverse does not substantially change as the CS acceleration factor is increased, which will allow for consistent reconstruction times regardless of the protocol. In addition, the HSS solver is non-iterative and the computational time should not be significantly affected by choice of CS penalty parameters. In the context of pre-specified MRI acquisition protocols, many factors for the HSS-Inverse method can be pre-computed and should enable clinically relevant reconstruction times (e.g., the computation of the inverse encoding matrix can be computed as part of a separate adjustment scan).

The idea of compact representations for the inverse of CS+SENSE reconstruction operators has been described above. This is accomplished through the use of a non-iterative HSS numerical technique. It is contemplated that the methods described here will be applicable to many reconstruction operators that rely on locality of interactions (e.g., wavelet transformations and GRAPPA-based parallel imaging).

Finally, the proposed HSS-Inverse method should be amenable to computationally demanding applications, such as cardiac imaging, and magnetic resonance fingerprinting (“MRF”), where the problem size can become very large due to the additional time dimension.

Referring particularly now toFIG. 4, an example of a magnetic resonance imaging (“MRI”) system400is illustrated. The MRI system400includes an operator workstation402, which will typically include a display404; one or more input devices406, such as a keyboard and mouse; and a processor408. The processor408may include a commercially available programmable machine running a commercially available operating system. The operator workstation402provides the operator interface that enables scan prescriptions to be entered into the MRI system400. In general, the operator workstation402may be coupled to four servers: a pulse sequence server410; a data acquisition server412; a data processing server414; and a data store server416. The operator workstation402and each server410,412,414, and416are connected to communicate with each other. For example, the servers410,412,414, and416may be connected via a communication system440, which may include any suitable network connection, whether wired, wireless, or a combination of both. As an example, the communication system440may include both proprietary or dedicated networks, as well as open networks, such as the internet.

The pulse sequence server410functions in response to instructions downloaded from the operator workstation402to operate a gradient system418and a radiofrequency (“RF”) system420. Gradient waveforms necessary to perform the prescribed scan are produced and applied to the gradient system418, which excites gradient coils in an assembly422to produce the magnetic field gradients Gx, Gy, and Gzused for position encoding magnetic resonance signals. The gradient coil assembly422forms part of a magnet assembly424that includes a polarizing magnet426and a whole-body RF coil428.

RF waveforms are applied by the RF system420to the RF coil428, or a separate local coil (not shown inFIG. 4), in order to perform the prescribed magnetic resonance pulse sequence. Responsive magnetic resonance signals detected by the RF coil428, or a separate local coil (not shown inFIG. 4), are received by the RF system420, where they are amplified, demodulated, filtered, and digitized under direction of commands produced by the pulse sequence server410. The RF system420includes an RF transmitter for producing a wide variety of RF pulses used in MRI pulse sequences. The RF transmitter is responsive to the scan prescription and direction from the pulse sequence server410to produce RF pulses of the desired frequency, phase, and pulse amplitude waveform. The generated RF pulses may be applied to the whole-body RF coil428or to one or more local coils or coil arrays (not shown inFIG. 4).

The pulse sequence server410also optionally receives patient data from a physiological acquisition controller430. By way of example, the physiological acquisition controller430may receive signals from a number of different sensors connected to the patient, such as electrocardiograph (“ECG”) signals from electrodes, or respiratory signals from a respiratory bellows or other respiratory monitoring device. Such signals are typically used by the pulse sequence server410to synchronize, or “gate,” the performance of the scan with the subject's heart beat or respiration.

The pulse sequence server410also connects to a scan room interface circuit432that receives signals from various sensors associated with the condition of the patient and the magnet system. It is also through the scan room interface circuit432that a patient positioning system434receives commands to move the patient to desired positions during the scan.

The digitized magnetic resonance signal samples produced by the RF system420are received by the data acquisition server412. The data acquisition server412operates in response to instructions downloaded from the operator workstation402to receive the real-time magnetic resonance data and provide buffer storage, such that no data is lost by data overrun. In some scans, the data acquisition server412does little more than pass the acquired magnetic resonance data to the data processor server414. However, in scans that require information derived from acquired magnetic resonance data to control the further performance of the scan, the data acquisition server412is programmed to produce such information and convey it to the pulse sequence server410. For example, during prescans, magnetic resonance data is acquired and used to calibrate the pulse sequence performed by the pulse sequence server410. As another example, navigator signals may be acquired and used to adjust the operating parameters of the RF system420or the gradient system418, or to control the view order in which k-space is sampled. In still another example, the data acquisition server412may also be employed to process magnetic resonance signals used to detect the arrival of a contrast agent in a magnetic resonance angiography (“MRA”) scan. By way of example, the data acquisition server412acquires magnetic resonance data and processes it in real-time to produce information that is used to control the scan.

The data processing server414receives magnetic resonance data from the data acquisition server412and processes it in accordance with instructions downloaded from the operator workstation402. Such processing may, for example, include one or more of the following: reconstructing two-dimensional or three-dimensional images by performing a Fourier transformation of raw k-space data; performing other image reconstruction algorithms, such as iterative or backprojection reconstruction algorithms; applying filters to raw k-space data or to reconstructed images; generating functional magnetic resonance images; calculating motion or flow images; and so on.

Images reconstructed by the data processing server414are conveyed back to the operator workstation402where they are stored. Real-time images are stored in a data base memory cache (not shown inFIG. 4), from which they may be output to operator display412or a display436that is located near the magnet assembly424for use by attending physicians. Batch mode images or selected real time images are stored in a host database on disc storage438. When such images have been reconstructed and transferred to storage, the data processing server414notifies the data store server416on the operator workstation402. The operator workstation402may be used by an operator to archive the images, produce films, or send the images via a network to other facilities.

The MRI system400may also include one or more networked workstations442. By way of example, a networked workstation442may include a display444; one or more input devices446, such as a keyboard and mouse; and a processor448. The networked workstation442may be located within the same facility as the operator workstation402, or in a different facility, such as a different healthcare institution or clinic.

The networked workstation442, whether within the same facility or in a different facility as the operator workstation402, may gain remote access to the data processing server414or data store server416via the communication system440. Accordingly, multiple networked workstations442may have access to the data processing server414and the data store server416. In this manner, magnetic resonance data, reconstructed images, or other data may exchanged between the data processing server414or the data store server416and the networked workstations442, such that the data or images may be remotely processed by a networked workstation442. This data may be exchanged in any suitable format, such as in accordance with the transmission control protocol (“TCP”), the internet protocol (“IP”), or other known or suitable protocols.

Referring now toFIG. 5, a block diagram of an example computer system500that can be configured to reconstruct magnetic resonance images using a hierarchically semiseparable (“HSS”) solver, as described above, is illustrated. The data from which the magnetic resonance images are reconstructed can be provided to the computer system500from the respective MRI system and received in a processing unit502.

In some embodiments, the processing unit502can include one or more processors. As an example, the processing unit502may include one or more of a digital signal processor (“DSP”)504, a microprocessor unit (“MPU”)506, and a graphics processing unit (“GPU”)508. The processing unit502can also include a data acquisition unit510that is configured to electronically receive data to be processed, which may include magnetic resonance image data. The DSP504, MPU506, GPU508, and data acquisition unit510are all coupled to a communication bus512. As an example, the communication bus512can be a group of wires, or a hardwire used for switching data between the peripherals or between any component in the processing unit502.

The DSP504can be configured to receive and processes the magnetic resonance data or reconstructed magnetic resonance images. The MPU506and GPU508can also be configured to process the magnetic resonance data or reconstructed magnetic resonance images in conjunction with the DSP504. As an example, the MPU506can be configured to control the operation of components in the processing unit502and can include instructions to perform reconstruction of the magnetic resonance image data on the DSP504. Also as an example, the GPU508can process image graphics.

In some embodiments, the DSP504can be configured to process the magnetic resonance image data received by the processing unit502in accordance with the algorithms described above. Thus, the DSP504can be configured to reconstruct magnetic resonance images using an HSS solver as described above.

The processing unit502preferably includes a communication port514in electronic communication with other devices, which may include a storage device516, a display518, and one or more input devices520. Examples of an input device520include, but are not limited to, a keyboard, a mouse, and a touch screen through which a user can provide an input.

The storage device516is configured to store images, whether provided to or processed by the processing unit502. The display518is used to display images, such as images that may be stored in the storage device516, and other information. Thus, in some embodiments, the storage device516and the display518can be used for displaying reconstructed magnetic resonance images.

The processing unit502can also be in electronic communication with a network522to transmit and receive data, including CT images, MR images, and other information. The communication port514can also be coupled to the processing unit502through a switched central resource, for example the communication bus512.

The processing unit502can also include a temporary storage524and a display controller526. As an example, the temporary storage524can store temporary information. For instance, the temporary storage524can be a random access memory.