WO2016014354A1 - Iterative analysis-based non-convex prior for enhanced sparse recovery - Google Patents

Abstract

A system for performing non-convex signal recovery using magnetic resonance imaging data includes an imaging device and a central computer unit. The imaging device comprises a plurality of coils configured to acquire a k-space dataset acquired using a magnetic resonance imaging device. The central computer unit is configured to apply an iterative reconstruction process to generate an image based on the k-space dataset, wherein the iterative reconstruction process minimizes a cost function that uses a non-convex pseudo-norm applied to a tight frame combined with a data fitting term.

Description

Iterative Analysis-Based Non-Convex Prior for Enhanced Sparse Recovery

CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] This application claims priority to U.S. provisional application Serial No.

62/028,846 filed July 25, 2014, which is incorporated herein by reference in its entirety. TECHNOLOGY FIELD

[0002] The present invention relates generally to methods, systems, and apparatuses for performing Magnetic Resonance Imaging (MRI) reconstruction with an iterative analysis-based non-convex prior. The disclosed techniques may be applied to enhance data recovery in sparse sampling applications.

BACKGROUND

[0003] The injection of convex priors for linear inverse problems marked a paradigm shift in signal processing reconstruction. The ability to inject sparsity using the Lj norm as a proxy for LQ norm marked the rapid development of convex optimization tools for sparse signal recovery. While the Li norm serves as an adequate proxy for the truly sparse L0 regularizer, the use of the iteratively reweighted Lj norm brings the Lj norm closer to its sparse L0 counterpart. Reweighted Li has proven to be quite successful in a host of applications from MRI reconstruction to Radio Interferometric Imaging.

[0004] Recent work has demonstrated the ability of non-convex prior (L_p) based recovery to surpass standard convex prior (Li) based recovery. Non-convex problems of this form have been shown to be especially beneficial in medical imaging. In enforcing these sparse priors, a subtlety revolves in the formulation of the problem. The problem can be formulated to recover either the signal itself assuming that its analysis coefficients in a given transform are sparse (analysis), or a sparse subset of transform coefficients from which the signal can be reconstructed (synthesis). Analysis-based signal priors outperform synthesis based priors in many image processing applications. Accordingly, it is desired to develop an efficient analysis-based technique for implementing non-convex prior based recovery.

SUMMARY

[0005] Embodiments of the present invention address and overcome one or more of the above shortcomings and drawbacks, by providing methods, systems, and apparatuses related to performing signal processing reconstruction with an iterative analysis-based non-convex prior. The disclosed techniques may be applied, for example, to enhance data recovery in Magnetic Resonance Imaging (MRI) sparse sampling applications.

[0006] According to some embodiments, a system for performing non-convex signal recovery using magnetic resonance imaging data includes an imaging device and a central computer unit. The imaging device comprises coils configured to acquire a k-space dataset using MRI techniques. The central computer unit is configured to apply an iterative reconstruction process to generate an image based on the k-space dataset, wherein the iterative reconstruction process minimizes a cost function that uses a non-convex pseudo-norm applied to a tight frame combined with a data fitting term. In some embodiments, the non-convex pseudo-norm is a p- norm with a p value less than 1. The iterative reconstruction process may minimize the cost function using any technique known in the art including, without limitation, a forward-backward splitting algorithm or an implementation of Dykstra's projection algorithm. [0007] In some embodiments of the aforementioned system, the cost function comprises a total variation (TV) prior. In these embodiments, the iterative reconstruction process may minimize the cost function using an implementation of Dykstra's projection algorithm which enforces the non-convex pseudo-norm. In some embodiments, a gradient TV proximal function (e.g., a Chambolle Pock gradient TV projection) is applied prior to enforcing the non-convex pseudo-norm. In some embodiments, the implementation of Dykstra's projection algorithm comprises proximity operators applied in parallel in a distributed computing environment. These proximity operators may include, for example, a first proximity operator indicator function of a convex set generated by the tight frame, a second proximity operator corresponding to a sparse analysis prior, and a third proximity operator corresponding to a TV regularizer.

[0008] In some embodiments of the aforementioned system, the cost function results in super-resolution of the image and the cost function comprises a downsampling matrix, a blur matrix, and a concatenation of a plurality of wavelets. In some embodiments, the cost function comprises one or more values applicable to denoise the k-space dataset during generation of the image.

[0009] According to other embodiments of the present invention, an article of manufacture for performing non-convex signal recovery comprises a non-transitory, tangible computer- readable medium holding computer-executable instructions for performing a method which includes receiving a sparse image dataset and applying an iterative reconstruction process to generate an image based on the sparse image dataset. The iterative reconstruction process minimizes a cost function that uses a non-convex pseudo-norm applied to a tight frame combined with a data fitting term. In some embodiments, the non-convex pseudo-norm is a p-norm with a p value less than 1. The iterative reconstruction process may minimize the cost function, for example, using a forward-backward splitting algorithm or an implementation of Dykstra's projection algorithm.

[0010] In some embodiments of the aforementioned article of manufacture the cost function comprises a total variation (TV) prior. In these embodiments, the iterative reconstruction process may minimize the cost function using an implementation of Dykstra's projection algorithm which enforces the non-convex pseudo-norm. In one embodiment, the implementation of Dykstra's projection algorithm comprises a plurality of proximity operators applied in parallel in a distributed computing environment.

[0011] According to other embodiments of the present invention, a system for performing non-convex signal recovery using imaging data includes an imaging device and a computer device comprising one or more processors. The imaging device is configured to acquire a sparse image dataset. The computer device is configured to apply an iterative reconstruction process to generate an image based on the sparse image dataset. This iterative reconstruction process minimizes a cost function that uses a p-norm with a p value less than 1 applied to a tight frame combined with a data fitting term.

[0012] Additional features and advantages of the invention will be made apparent from the following detailed description of illustrative embodiments that proceeds with reference to the accompanying drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

[0013] The foregoing and other aspects of the present invention are best understood from the following detailed description when read in connection with the accompanying drawings. For the purpose of illustrating the invention, there is shown in the drawings embodiments that are presently preferred, it being understood, however, that the invention is not limited to the specific instrumentalities disclosed. Included in the drawings are the following Figures:

[0014] FIG. 1 shows a system for ordering acquisition of frequency domain components representing magnetic resonance image data for storage in a k-space storage array, as used by some embodiments of the present invention;

[0015] FIG. 2A shows an example method which implements IAN with forward-backward splitting, according to some embodiments;

[0016] FIG. 2B shows an example method which implements IAN with Dykstra-like splitting, according to some embodiments;

[0017] FIG. 3 shows a method which implements IAN with a total variation (TV) prior, according to some embodiments;

[0018] FIG. 4 shows an alternative method where a parallel proximal Dykstra framework is used to implement IAN with a total variation (TV) prior, according to some embodiments;

[0019] FIG. 5 includes a table showing experimental results comparing IAN to other super- resolution techniques;

[0020] FIG. 6 shows MRI angiography super-resolution results using IAN and state of the art super-resolution techniques;

[0021] FIG. 7 shows a comparison of MRI tumor super-resolution results generated by IAN and conventional super-resolution techniques; [0022] FIG. 8 provides a comparison of IAN-based MRI undersampled reconstruction to conventional undersampled reconstruction techniques;

[0023] FIG. 9A shows the results of undersampled computed tomography images;

[0024] FIG. 9B shows the quantitative results for each recovery shown in FIG. 9A;

[0025] FIG. 10 shows a set of images demonstrating the use of IAN to denoise noisy X-ray fluoroscopic data, according to some embodiments;

[0026] FIG. 11 shows of set of images demonstrating the use of IAN to denoise noisy MRI angiography data, according to some embodiments; and

[0027] FIG. 12 illustrates an exemplary computing environment within which embodiments of the invention may be implemented.

DETAILED DESCRIPTION OF ILLUSTRATIVE EMBODIMENTS

[0028] The following disclosure describes the present invention according to several embodiments directed at methods, systems, and apparatuses related to solving a generic, analysis-based, non-convex prior for ill-posed inverse problems related to Magnetic Resonance Imaging (MRI) image reconstruction. The techniques described herein are generally referred to as Iterative Analysis-Based Non-Convex Prior for Enhanced Sparse Recovery, or "IAN." IAN leverages the benefits of tight frames and utilizes an efficient computational architecture for analysis-based signal recovery.

[0029] As described in greater detail below, this generic non-convex framework may be applied to three challenging image processing applications: single image super-resolution, medical image undersampled recovery, and image denoising. For single image super-resolution, IAN may be used to super-resolve pathologically relevant features in medical images to aid in clinical diagnosis. For undersampled medical image recovery, IAN outperforms competing methods while circumventing the standard artifacts that are common to subsampled imaging. In some embodiments, IAN may also be extended to incorporate a total variation regularizer to deal with intrinsic hardware noise generated in MRI, CT and Fluoroscopic medical images. IAN, with a TV regularizer, is able to denoise raw medical images and can be used without oracle knowledge of the noise variance. In order to overcome the challenge of parameter tuning for image denoising, a Monte Carlo based Stein Unbiased Risk Estimator (SURE) is used in some embodiments for automated parameter tuning. Finally, in order to tackle instances in which the forward imaging model is unknown, IAN can be used with existing blind deconvolution methods in order to estimate the parameters of the underlying model. IAN comes at a computational cost that is at least a factor of 2 less than that of competing non-convex algorithms.

[0030] FIG. 1 shows a system 100 for ordering acquisition of frequency domain components representing MRI data for storage in a k-space storage array, as used by some embodiments of the present invention. In system 100, magnetic coils 12 create a static base magnetic field in the body of patient 1 1 to be imaged and positioned on a table. Within the magnet system are gradient coils 14 for producing position dependent magnetic field gradients superimposed on the static magnetic field. Gradient coils 14, in response to gradient signals supplied thereto by a gradient and shim coil control module 16, produce position dependent and shimmed magnetic field gradients in three orthogonal directions and generates magnetic field pulse sequences. The shimmed gradients compensate for inhomogeneity and variability in an MRI device magnetic field resulting from patient anatomical variation and other sources. The magnetic field gradients include a slice-selection gradient magnetic field, a phase-encoding gradient magnetic field and a readout gradient magnetic field that are applied to patient 1 1.

[0031] Further radio frequency (RF) module 20 provides RF pulse signals to RF coil 18, which in response produces magnetic field pulses which rotate the spins of the protons in the imaged body of the patient 1 1 by ninety degrees or by one hundred and eighty degrees for so- called "spin echo" imaging, or by angles less than or equal to 90 degrees for so-called "gradient echo" imaging. Gradient and shim coil control module 16 in conjunction with RF module 20, as directed by central control unit 26, control slice-selection, phase-encoding, readout gradient magnetic fields, radio frequency transmission, and magnetic resonance signal detection, to acquire magnetic resonance signals representing planar slices of patient 1 1.

[0032] In response to applied RF pulse signals, the RF coil 18 receives magnetic resonance signals, i.e., signals from the excited protons within the body as they return to an equilibrium position established by the static and gradient magnetic fields. The magnetic resonance signals are detected and processed by a detector within RF module 20 and k-space component processor unit 34 to provide a magnetic resonance dataset to an image data processor for processing into an image. In some embodiments, the image data processor is located in central control unit 26. However, in other embodiments such as the one depicted in FIG. 1 , the image data processor is located in a separate unit 27. ECG synchronization signal generator 30 provides ECG signals used for pulse sequence and imaging synchronization. A two or three dimensional k-space storage array of individual data elements in k-space component processor unit 34 stores corresponding individual frequency components comprising a magnetic resonance dataset. The k-space array of individual data elements has a designated center and individual data elements individually have a radius to the designated center. [0033] A magnetic field generator (comprising coils 12, 14, and 18) generates a magnetic field for use in acquiring multiple individual frequency components corresponding to individual data elements in the storage array. The individual frequency components are successively acquired in an order in which the radius of respective corresponding individual data elements increases and decreases along a substantially spiral path as the multiple individual frequency components is sequentially acquired during acquisition of a magnetic resonance dataset representing a magnetic resonance image. A storage processor in the k-space component processor unit 34 stores individual frequency components acquired using the magnetic field in corresponding individual data elements in the array. The radius of respective corresponding individual data elements alternately increases and decreases as multiple sequential individual frequency components are acquired. The magnetic field acquires individual frequency components in an order corresponding to a sequence of substantially adjacent individual data elements in the array and magnetic field gradient change between successively acquired frequency components is substantially minimized.

[0034] Central control unit 26 uses information stored in an internal database to process the detected magnetic resonance signals in a coordinated manner to generate high quality images of a selected slice(s) of the body (e.g., using the image data processor) and adjusts other parameters of system 100. The stored information comprises predetermined pulse sequence and magnetic field gradient and strength data as well as data indicating timing, orientation and spatial volume of gradient magnetic fields to be applied in imaging. Generated images are presented on display 40 of the operator interface. Computer 28 of the operator interface includes a graphical user interface (GUI) enabling user interaction with central control unit 26 and enables user modification of magnetic resonance imaging signals in substantially real time. Continuing with reference to FIG. 1, display processor 37 processes the magnetic resonance signals to reconstruct one or more images for presentation on display 40, for example. Various techniques may be used for reconstruction. For example, as described in greater detail below, an optimization algorithm is applied to iteratively solve a cost function which results in the reconstructed image.

[0035] According to various embodiments described herein, the system 100 illustrated in FIG. 1 is configured to apply IAN to determine an algorithmic solution for the following analysis non-convex problem:

where 0 < p < l ; ^4 is a linear operator which maps (Cⁿ→ C^m(m « n), and y e (C^m. W is a tight frame such that W^TW = I, and y G (C^mxl. In some embodiments, IAN is combined with a Total Variation (TV) prior. This combination is referred to herein as "IAN with TV." As is generally understood in the art, a TV prior seeks to smooth the reconstructed image while preserving high frequency edges. IAN with TV seeks to solve the following optimization problem: min_xeCn ( ^ \\Ax - y\\ +λ₁ \\Wx\\ (2)

The TV operator takes into account the horizontal differences and the vertical differences in the image. The manner in which these differences are combined define isotropic and anisotropic TV. In isotropic TV, the differences are combined by taking the square root of the squared and summed vertical and horizontal differences. In anisotropic TV, the absolute values of the horizontal and vertical difference are summed in the computation of the TV operator. [0036] IAN solves Equations (1) and (2) when p < 1 by integrating a proximal splitting framework with an iteratively estimated optimal shrinkage operator for a given L_p pseudo-norm. IAN leverages the benefit of tight frames and uses a computationally efficient architecture that is intrinsically parallelizable and is significantly cheaper than conventional convex and non-convex regularized inverse problem solvers. IAN may be applied, for example, to ill-posed inverse problems such as single image super-resolution, undersampled medical image reconstruction, and image denoising. In order to handle image denoising, IAN may be extended to handle a total variation regularizer (both isotropic and anisotropic). In order to handle realistic scenarios when the variance of the noise is unknown for image denoising, IAN can be incorporated into an existing SURE for estimation of the mean square error. This allows for parameter tuning without explicit need for the ground truth. In instances when the parameters of the forward model are unknown for super-resolution (e.g., point spread function of a system), IAN can be used with standard blind deconvolution splitting approaches to recover the unknown signal and point spread function.

[0037] When solving the optimization problem presented in Equation 1 , a natural first pass algorithmic formulation would be the forward-backward splitting approach. An example method 200 which implements forward-backward splitting is shown in FIG. 2A. The proximal L_p pseudo-norm used in FIG. 2A (see line 4 of FIG. 2A) is general and any shrinkage function generally known in the art can be used.

[0038] FIG. 2B shows an alternative method 205 which implements IAN with Dykstra-like splitting to solve Equation 1 , according to some embodiments of the present invention. In this algorithm, IAN seeks to solve the non-convex optimization problem set forth in Equation 1. As with the example shown in FIG. 2A, the proximal L_p pseudo-norm used in FIG. 2B (see line 6 of FIG. 2B) is general and any shrinkage function generally known in the art can be used. The method 205 presented in FIG. 2B adapts the splitting methodology used in some conventional forward splitting algorithms with a proximal step to ensure that the enforcement of the L_p pseudo-norm is constrained to an e ball centered on the previous fidelity enforcement. This form of Dykstra-like algorithm solves the problem min_xeCn ( (x) + g(x) + - \\x— r|| ^) with the imposed feasibility condition that dom( ( )) Π dom(,g( ))≠ 0. Within the context of solving an analysis non-convex prior, f(x) takes the form of an indicator function of the convex set generated by the tight frame, g(x) takes the form of the sparse analysis prior and V is obtained from the previous splitting step which enforces data fidelity. Thus, the optimization problem can be expressed as min_xeCn (^t-c + \\Wx \\_p + - \\x— HI It should be noted that the t parameter used in the algorithm 205 is the step size which is defined to be the inverse of the Lipschitz constant. The Lipschitz constant can be exactly computed using the power iteration on the A^TA operator to extract the largest eigenvalue of the symmetric operator.

[0039] FIG. 3 shows a method 300 that may be used for solving the optimization problem set forth in Equation 2. The method 300 in FIG. 3 is similar to that presented in FIG. 2. However, in FIG 3, a gradient TV proximal initial step prior is used to enforce the proximal L_p pseudo- norm (see line 3 of FIG. 3). Any proximal gradient TV step can be used in the method including, for example, the Chambolle Pock gradient TV projection.

[0040] FIG. 4 shows an alternative method 400 where a parallel proximal Dykstra framework is used to solve the optimization problem set forth in Equation 2, according to some embodiments of the present invention. The parallel Dykstra proximal framework is designed to solve problems of the form min_x£Cn (∑f₌₁ ω^{χ) + - ||x - r|| ¾ with the imposed feasibility condition that dom( j(:t))≠ 0 and = 1. Within the context of solving an analysis non-convex prior, takes the form of an indicator function of the convex set generated by the tight frame (see line 6 of FIG. 4), f₂ (x)takes the form of the sparse analysis prior (see line 7 of FIG. 4), ₃ ( ) is the imposed TV regularizer (see line 7 of FIG. 4). The variable V used at line 4 represents the result of the previous fidelity enforcement. Let ω_έ be defined as the normalized regularization parameters. Thus, the problem then takes the form min_xeCn +

<^y ₂l| x||_p + ω₃ ||χ||_Γ + i ||x - ΗΙ^)· Iⁿ addition to the advantages previously enumerated, custom edge detection filters can be developed and utilized in some embodiments, provided the filters meet the constraint tr(HV^T) = 0 where H and V are the horizontal and vertical edge detection filters, respectively. Depending on the size of the edge detection filters used, 3 x 3 or 2 X 2 based block processing may be necessary for this approach. In order to prevent artifacts along block boundaries, this method 400 may be repeated for shifts of each block and with subsequent averaging of the results obtained from each stage.

[0041] The method 400 shown in FIG 4 offers several benefits over the method 300 in FIG. 3 including greater versatility in TV filter selection, choice of isotropic versus anisotropic TV, a comparatively better theoretical guarantee, and the ability to be implemented in parallel processing environments. Thus, in embodiments of the present invention implemented on a distributed computing architecture the method 400 shown in FIG. 4 may be used. Conversely, the method 300 shown in FIG. 3 may be used in embodiments where IAN is implemented on a single core architecture.

[0042] The follow paragraphs described three potential applications for the non-convex prior used in IAN: single image super-resolution, undersampled medical image reconstruction, and image denoising. It should be noted that, although these applications are discussed independently, the various techniques associated with each application may be combined in some embodiments (e.g., to result in an undersampled medical image reconstruction that applies both super-resolution and denoising functionality).

[0043] For super-resolution applications, the problem formulation with IAN may be described with the following equation: m!n \\y— DBx jj^ H- λ|| τ ? ||_ρ (3)

In Equation 5, y is the observed low- resolution image, D is the downsampling matrix, B is the blur matrix, and ψ is a concatenation of multiple wavelets. The problem formulation for IAN with TV may be expressed as: nam j j y - DBx \ | | + Λχ j j φχ j j _p + A₂11 TV(x) | j ₁ (4)

Equations 3 and 4 may be solved using the general IAN and IAN with TV techniques described above with respect to FIGS. 2 - 4, respectively. The applicability of IAN single image super- resolution is motivated by examining medical images for which super-resolution would greatly aid clinical diagnosis. For example, the techniques may be applied in some embodiments to super-resolve an MRI image of cerebral tumor and angiographic images with fine capillary details. Also in the medical domain, IAN may be applied in some embodiments to super-resolve guidewires in CT fluoroscopic images which is of significant clinical interest to physicians.

[0044] FIG. 5 includes a table showing experimental results comparing IAN to other super- resolution techniques. To generate a low resolution image for each of these comparisons, images were blurred with a 3x3 Gaussian kernel and then downsampled in the pixel domain by a factor of 2 for each dimension. The reconstructed results from IAN are compared to standard interpolation based approaches (i.e., bilinear, bicubic, box kernel, Lanczos 2 kernel, Lanczos 3 kernel), the method of sparse mixing estimators (SME), sparse representation based nonlocal autoregressive modelling (NARM), spatially adaptive iterative singular value thresholding (SAIST), and the fast iterative shrinkage and thresholding algorithm (FISTA). To quantify the accuracy of image recovery, two values were calculated. First, the peak signal to noise ratio was calculated according to the following equation:

PSNR (x, x) = (20 1og₁₀ max ( )) + 10 log₁₀ (MSE(x, )), where MSE is mean square error, x is the reference signal and x is the recovered signal. The second value used for comparison is the structure similarity metric (SSIM), defined using techniques generally known in the art.

[0045] FIG. 6 shows MRI Angiography super-resolution results using IAN and state of the art super-resolution techniques. Image 605 shows an observed low resolution image. Image 610 shows the results recovered using the method of sparse mixing estimators (PSNR: 34.67 dB, SSIM: 0.916. Image 615 shows the results recovered using Sparsity Averaging Reweighted Analysis (SARA) (PSNR: 37.64 dB, SSIM: 0.942), while Image 620 is the results recovered using IAN (PSNR: 38.03 dB, SSIM: 0.943).

[0046] Continuing with reference to FIG. 6, Image 625 shows the SME residual relative to ground truth. Image 630 is the SARA residual relative to ground truth, while Image 635 IAN residual relative to ground truth. (All residuals have been amplified by a factor of 2 to aid in visualization). The highlighted regions in the residual images draw attention to overall arterial structure and bifurcations that are super-resolved with highest fidelity using IAN relative to the next best reconstruction method (i.e., SARA). These regions may be used, for example, by clinicians dealing with potential embolic plaques / aneurysms.

[0047] FIG. 7 shows a comparison of MRI tumor super-resolution results generated IAN and conventional super-resolution techniques. The tumor can be observed as the small area of bright intensity in the right temporal lobe. Image 705 shows the observed low resolution image. Image 710 shows the results recovered using the method of Sparse Mixing Estimators (SME) (PSNR: 40.31 dB, SSIM: 0.987), while Image 715 was recovered using SARA (PSNR: 45.05 dB, SSIM: 0.995). Finally Image 720 was recovered using IAN (PSNR: 47.02 dB, SSIM: 0.995).

[0048] Continuing with reference to FIG. 7, Image 725 shows the SME residual relative to ground truth. Images 730 and 735 show SARA and IAN, respectively, residual relative to ground truth. All residuals have been amplified by a factor of 10 to enhance visual observation of reconstruction differences. From the residuals and numerical reconstruction metrics, it is evident that IAN recovers the high resolution image with the highest fidelity. Of greater clinical interest is IAN's recovery of tumor details and margins with the least errors of all techniques.

[0049] Undersampling refers to acquisition of data with far fewer measurements than that dictated by the Nyquist rate. In medical imaging, undersampled acquisition results in faster scan times which, in turn, mitigate motion artifacts. Motion artifacts are a significant challenge in medical imaging and are the cause of image streaks and ghosting artifacts. Patient motion (internal and external) necessitates image registration between long acquisitions which further reduces the interpretability of the data. IAN is capable of reconstructing undersampled pathological and regular MRI and CT medical images with greater fidelity than both state of the art convex and non-convex methods. [0050] FIG. 8 provides a comparison of IAN-based MRI undersampled reconstruction to conventional undersampled reconstruction techniques. Image 805 shows a ground truth image. Image 810 is reconstructed using SARA (PSNR: 34.70 dB) and image 815 is reconstructed using IAN but with a thresholding function (PSNR: 36.88 dB). Image 820 was reconstructed using IAN (PSNR: 38.73 dB). Image 825 is a SARA residual image. Image 830 shows IAN with RC shrinkage function residual image, and Image 835 is the IAN residual image. The overset zoomed in image of the two insets represent the wavelet artifacts that occur, on occasion, when using SARA that is otherwise absent in IAN. The same artifacts were seen in some of SARA's super-resolution reconstruction (results not shown). It is interesting to note that IAN reconstructs the image with the highest quantitative and qualitative fidelity. A cursory examination of the residual image shows that the thresholding function causes reconstruction errors along the arterial walls of the aorta.

[0051] FIG. 9A shows the results of undersampled (45 views) computed tomography images. The first column depicts the ground truth. The second column shows the reconstruction afforded by standard filtered back projection method. The third column contains the images reconstructed using standard least squares (gradient descent). Sparse reconstruction with a synthesis prior (using FISTA) is shown in the fourth column. Sparse non-convex reconstruction using an analysis reweighted L_\ prior is shown in the fifth column while the sixth column contains the reconstruction afforded by analysis-based non-convex L_p prior (IAN TV). The quantitative results for each recovery can be examined in the table shown in FIG. 9B. In each case, it is clear that the combination of non-convex L_p and a TV prior outperforms conventional reweighted L_\ and other sparse regularization techniques. The pure L_v (without TV) is also shown in FIG. 9B. For the images shown, it is clear that streaking artifacts that are common to undersampled CT reconstruction can be eliminated using a non-convex L_p and TV prior.

[0052] Medical images are noisy given the inability to control factors such as thermal noise, electromagnetic interference and other non-laboratory conditions. IAN can be used to denoise ground truth MRI and CT Fluoroscopic images with the explicit goal of aiding clinical diagnosis. For denoising applications, the problem formulation with IAN may be described with the following equation: miii \\y— + λ { | ψ,τ [ [_ρ (5)

In Equation 6, y is the observed noisy image and ψ is a concatenation of multiple wavelets. The problem formulation for IAN with TV may be expressed as: miii \\y ------ x\\ ¾ --^■}--^■ λ-j ||^> ^*|L -}-- A2 ||TV^r(a;) || ] ⁽⁶⁾

Equations 5 and 6 may be solved using the general IAN and IAN with TV techniques described above with respect to FIGS. 2 - 4, respectively. For image denoising applications where the noiseless ground truth is not available, SURE can be used to achieve reference less parameter tuning for IAN with TV.

[0053] FIG. 10 demonstrates the use of IAN to denoise a noisy X-ray fluoroscopic image with no a priori knowledge about noise statistics. In fluoroscopic images, contrast is of critical importance as physicians inject contrast agents continuously to determine catheter location relative to anatomical markers. Noise in such images reduces contrast and significantly increases the need for exposure to large quantities of contrast agents. FIG. 10 shows that IAN denoises the provided fluoroscopic image better than conventional methods. Image 1005 is the ground truth noisy image. Image 1010 shows the results of denoising using BM3D, a technique generally known in the art. Images 1015 and 1020 show the results of denoising with IAN for the best visual reconstruction result and the best IQM result, respectively. The insets show that IAN visually denoises relevant features better than BM3D.

[0054] FIG. 1 1 shows an example of denoising an MRI angiography image with nonuniform noise. In this specific case, the patient suffers from a cerebral aneurysm which can be seen in the right hemisphere near the corpus callosum. Clear boundaries are of vital clinical importance as it used to determine clinical prognosis and treatment paradigms. Image 1 105 shows the ground truth noisy image. Image 1 1 10 is an image denoised using BM3D. Image 1 1 15 is denoised using IAN (best visual reconstruction). A cursory examination of the image reveals the superior visual denoising by IAN relative to BM3D. For example, the arterial out pouching and its anatomical boundaries are clearest in the image denoised by IAN.

[0055] FIG. 12 illustrates an exemplary computing environment 1200 within which embodiments of the invention may be implemented. For example, this computing environment 1200 may be used to implement the methods 200, 300, 400 described in FIGS. 2, 3, and 4, respectively. In some embodiments, the computing environment 1200 may be used to implement one or more of the components illustrated in the system 100 of FIG. 1. The computing environment 1200 may include computer system 1210, which is one example of a computing system upon which embodiments of the invention may be implemented. Computers and computing environments, such as computer system 1210 and computing environment 1200, are known to those of skill in the art and thus are described briefly here. [0056] As shown in FIG. 12, the computer system 1210 may include a communication mechanism such as a bus 1221 or other communication mechanism for communicating information within the computer system 1210. The computer system 1210 further includes one or more processors 1220 coupled with the bus 1221 for processing the information. The processors 1220 may include one or more central processing units (CPUs), graphical processing units (GPUs), or any other processor known in the art.

[0057] The computer system 1210 also includes a system memory 1230 coupled to the bus 1221 for storing information and instructions to be executed by processors 1220. The system memory 1230 may include computer readable storage media in the form of volatile and/or nonvolatile memory, such as read only memory (ROM) 1231 and/or random access memory (RAM) 1232. The system memory RAM 1232 may include other dynamic storage device(s) (e.g., dynamic RAM, static RAM, and synchronous DRAM). The system memory ROM 1231 may include other static storage device(s) (e.g., programmable ROM, erasable PROM, and electrically erasable PROM). In addition, the system memory 1230 may be used for storing temporary variables or other intermediate information during the execution of instructions by the processors 1220. A basic input/output system (BIOS) 1233 containing the basic routines that help to transfer information between elements within computer system 1210, such as during start-up, may be stored in ROM 1231. RAM 1232 may contain data and/or program modules that are immediately accessible to and/or presently being operated on by the processors 1220. System memory 1230 may additionally include, for example, operating system 1234, application programs 1235, other program modules 1236 and program data 1237.

[0058] The computer system 1210 also includes a disk controller 1240 coupled to the bus 1221 to control one or more storage devices for storing information and instructions, such as a hard disk 1241 and a removable media drive 1242 (e.g., floppy disk drive, compact disc drive, tape drive, and/or solid state drive). The storage devices may be added to the computer system 1210 using an appropriate device interface (e.g., a small computer system interface (SCSI), integrated device electronics (IDE), Universal Serial Bus (USB), or FireWire).

[0059] The computer system 1210 may also include a display controller 1265 coupled to the bus 1221 to control a display 1266, such as a cathode ray tube (CRT) or liquid crystal display (LCD), for displaying information to a computer user. The computer system includes an input interface 1260 and one or more input devices, such as a keyboard 1262 and a pointing device 1261 , for interacting with a computer user and providing information to the processor 1220. The pointing device 1261 , for example, may be a mouse, a trackball, or a pointing stick for communicating direction information and command selections to the processor 1220 and for controlling cursor movement on the display 1266. The display 1266 may provide a touch screen interface which allows input to supplement or replace the communication of direction information and command selections by the pointing device 1261.

[0060] The computer system 1210 may perform a portion or all of the processing steps of embodiments of the invention in response to the processors 1220 executing one or more sequences of one or more instructions contained in a memory, such as the system memory 1230. Such instructions may be read into the system memory 1230 from another computer readable medium, such as a hard disk 1241 or a removable media drive 1242. The hard disk 1241 may contain one or more datastores and data files used by embodiments of the present invention. Datastore contents and data files may be encrypted to improve security. The processors 1220 may also be employed in a multi-processing arrangement to execute the one or more sequences of instructions contained in system memory 1230. In alternative embodiments, hard-wired circuitry may be used in place of or in combination with software instructions. Thus, embodiments are not limited to any specific combination of hardware circuitry and software.

[0061] As stated above, the computer system 1210 may include at least one computer readable medium or memory for holding instructions programmed according to embodiments of the invention and for containing data structures, tables, records, or other data described herein. The term "computer readable medium" as used herein refers to any medium that participates in providing instructions to the processor 1220 for execution. A computer readable medium may take many forms including, but not limited to, non-volatile media, volatile media, and transmission media. Non-limiting examples of non-volatile media include optical disks, solid state drives, magnetic disks, and magneto-optical disks, such as hard disk 1241 or removable media drive 1242. Non-limiting examples of volatile media include dynamic memory, such as system memory 1230. Non-limiting examples of transmission media include coaxial cables, copper wire, and fiber optics, including the wires that make up the bus 1221. Transmission media may also take the form of acoustic or light waves, such as those generated during radio wave and infrared data communications.

[0062] The computing environment 1200 may further include the computer system 1210 operating in a networked environment using logical connections to one or more remote computers, such as remote computer 1280. Remote computer 1280 may be a personal computer (laptop or desktop), a mobile device, a server, a router, a network PC, a peer device or other common network node, and typically includes many or all of the elements described above relative to computer system 1210. When used in a networking environment, computer system 1210 may include modem 1272 for establishing communications over a network 1271 , such as the Internet. Modem 1272 may be connected to bus 1221 via user network interface 1270, or via another appropriate mechanism.

[0063] Network 1271 may be any network or system generally known in the art, including the Internet, an intranet, a local area network (LAN), a wide area network (WAN), a metropolitan area network (MAN), a direct connection or series of connections, a cellular telephone network, or any other network or medium capable of facilitating communication between computer system 1210 and other computers (e.g., remote computer 1280). The network 1271 may be wired, wireless or a combination thereof. Wired connections may be implemented using Ethernet, Universal Serial Bus (USB), RJ-1 1 or any other wired connection generally known in the art. Wireless connections may be implemented using Wi-Fi, WiMAX, and Bluetooth, infrared, cellular networks, satellite or any other wireless connection methodology generally known in the art. Additionally, several networks may work alone or in communication with each other to facilitate communication in the network 1271.

[0064] The embodiments of the present disclosure may be implemented with any combination of hardware and software. In addition, the embodiments of the present disclosure may be included in an article of manufacture (e.g., one or more computer program products) having, for example, computer-readable, non-transitory media. The media has embodied therein, for instance, computer readable program code for providing and facilitating the mechanisms of the embodiments of the present disclosure. The article of manufacture can be included as part of a computer system or sold separately.

[0065] While various aspects and embodiments have been disclosed herein, other aspects and embodiments will be apparent to those skilled in the art. The various aspects and embodiments disclosed herein are for purposes of illustration and are not intended to be limiting, with the true scope and spirit being indicated by the following claims.

[0066] An executable application, as used herein, comprises code or machine readable instructions for conditioning the processor to implement predetermined functions, such as those of an operating system, a context data acquisition system or other information processing system, for example, in response to user command or input. An executable procedure is a segment of code or machine readable instruction, sub-routine, or other distinct section of code or portion of an executable application for performing one or more particular processes. These processes may include receiving input data and/or parameters, performing operations on received input data and/or performing functions in response to received input parameters, and providing resulting output data and/or parameters.

[0067] A graphical user interface (GUI), as used herein, comprises one or more display images, generated by a display processor and enabling user interaction with a processor or other device and associated data acquisition and processing functions. The GUI also includes an executable procedure or executable application. The executable procedure or executable application conditions the display processor to generate signals representing the GUI display images. These signals are supplied to a display device which displays the image for viewing by the user. The processor, under control of an executable procedure or executable application, manipulates the GUI display images in response to signals received from the input devices. In this way, the user may interact with the display image using the input devices, enabling user interaction with the processor or other device. [0068] The functions and process steps herein may be performed automatically or wholly or partially in response to user command. An activity (including a step) performed automatically is performed in response to one or more executable instructions or device operation without user direct initiation of the activity.

[0069] The system and processes of the figures are not exclusive. Other systems, processes and menus may be derived in accordance with the principles of the invention to accomplish the same objectives. Although this invention has been described with reference to particular embodiments, it is to be understood that the embodiments and variations shown and described herein are for illustration purposes only. Modifications to the current design may be implemented by those skilled in the art, without departing from the scope of the invention. As described herein, the various systems, subsystems, agents, managers and processes can be implemented using hardware components, software components, and/or combinations thereof. No claim element herein is to be construed under the provisions of 35 U.S.C. 1 12, sixth paragraph, unless the element is expressly recited using the phrase "means for."

Claims

CLAIMS We claim:

1. A system for performing non-convex signal recovery using magnetic resonance imaging data, the system comprising:

an imaging device comprising a plurality of coils configured to acquire a k-space dataset acquired using a magnetic resonance imaging device; and

a central computer unit configured to apply an iterative reconstruction process to generate an image based on the k-space dataset, wherein the iterative reconstruction process minimizes a cost function that uses a non-convex pseudo-norm applied to a tight frame combined with a data fitting term.

2. The system of claim 1 , wherein the non-convex pseudo-norm is a p-norm with a p value less than 1.

3. The system of claim 2, wherein the iterative reconstruction process minimizes the cost function using a forward-backward splitting algorithm.

4. The system of claim 2, wherein the iterative reconstruction process minimizes the cost function using an implementation of Dykstra's projection algorithm.

5. The system of claim 2, wherein the cost function comprises a total variation (TV) prior.

6. The system of claim 5, wherein the iterative reconstruction process minimizes the cost function using an implementation of Dykstra's projection algorithm which enforces the non- convex pseudo-norm.

7. The system of claim 6, wherein a gradient TV proximal function is applied prior to enforcing the non-convex pseudo-norm.

8. The system of claim 7, wherein the gradient TV proximal function comprises a

Chambolle Pock gradient TV projection.

9. The system of claim 6, wherein the implementation of Dykstra's projection algorithm comprises a plurality of proximity operators applied in parallel in a distributed computing environment.

10. The system of claim 9, wherein the plurality of proximity operators comprise:

a first proximity operator indicator function of a convex set generated by the tight frame, a second proximity operator corresponding to a sparse analysis prior, and

a third proximity operator corresponding to a TV regularizer.

11. The system of claim 1 , wherein application of the cost function results in super-resolution of the image and the cost function comprises a downsampling matrix, a blur matrix, and a concatenation of a plurality of wavelets.

12. The system of claim 1 , wherein the cost function comprises one or more values applicable to denoise the k-space dataset during generation of the image.

13. An article of manufacture for performing non-convex signal recovery, the article of manufacture comprising a non-transitory, tangible computer-readable medium holding computer-executable instructions for performing a method comprising:

receiving a sparse image dataset;

applying an iterative reconstruction process to generate an image based on the sparse image dataset, wherein the iterative reconstruction process minimizes a cost function that uses a non-convex pseudo-norm applied to a tight frame combined with a data fitting term.

14. The article of manufacture of claim 13, wherein the non-convex pseudo-norm is a p-norm with a p value less than 1.

15. The article of manufacture of claim 14, wherein the iterative reconstruction process minimizes the cost function using a forward-backward splitting algorithm.

16. The article of manufacture of claim 14, wherein the iterative reconstruction process minimizes the cost function using an implementation of Dykstra's projection algorithm.

17. The article of manufacture of claim 14, wherein the cost function comprises a total variation (TV) prior.

18. The article of manufacture of claim 17, wherein the iterative reconstruction process minimizes the cost function using an implementation of Dykstra's projection algorithm which enforces the non-convex pseudo-norm.

19. The article of manufacture of claim 18, wherein the implementation of Dykstra's projection algorithm comprises a plurality of proximity operators applied in parallel in a distributed computing environment.

20. A system for performing non-convex signal recovery using imaging data, the system comprising:

an imaging device configured to acquire a sparse image dataset; and

a computer device comprising one or more processors configured to apply an iterative reconstruction process to generate an image based on the sparse image dataset, wherein the iterative reconstruction process minimizes a cost function that uses a p-norm with a p value less than 1 applied to a tight frame combined with a data fitting term.