Iterative Reconstruction from Undersampled Radial Data Using a Total-Variation Constraint, MATLAB Code

A nonlinear-optimization-based image reconstruction that allows for use of prior knowledge via penalty functions with data from multiple coils.

We are making available an algorithm for iterative reconstruction from undersampled radial MRI with a total-variation (TV) constraint.

Get the Code

The software available on this page is provided free of charge and comes without any warranty. CAI²R and NYU Grossman School of Medicine do not take any liability for problems or damage of any kind resulting from the use of the files provided. Operation of the software is solely at the user’s own risk. The software developments provided are not medical products and must not be used for making diagnostic decisions.

The software is provided for non-commercial, academic use only. Usage or distribution of the software for commercial purpose is prohibited. All rights belong to the author (Florian Knoll) and NYU Grossman School of Medicine. If you use the software for academic work, please give credit to the author in publications and cite the related publications.

Please spell out your affiliation (e.g. “New York University” rather than “NYU”).

Notes on Use

The package contains two different MATLAB programs:

  • iterativeradial_phantom.m demonstrates the reconstruction of a numerical Shepp-Logan phantom using iterative reconstruction with a TV constraint on the real-part of the image (corresponding to Figure 3 in the related publication).
  • iterativeradial_multicoil.m demonstrates the full algorithm with initial coil-estimation step and final joint-coil image calculation (see Figures 6-8).

Real MRI datasets are provided from a phantom and a brain scan, both acquired using the golden-angle ordering scheme, which allows retrospective undersampling by truncating spokes at the beginning of the data (note that it is better to truncate at the beginning due to steady-state effects). The optimization is done using the L-BFGS algorithm. The cost functions are implemented in separate files (costfunction_phanom.m for the Shepp-Logan case; costfunction_coils.m and costfunction_image.m for the full algorithm). The estimated vector x and the to-be-calculated gradient vector g are passed as 1D vectors and are reshaped using the helper functions img_to_vec and vec_to_img. In all cost functions, preconditioning with the DCF is used to accelerate convergence. The penalty terms are implemented in external helper functions (TV.m for the total variation of the real part, L2D.m for the L2-based smoothness penalty for the coil profiles). These methods have been implemented for simplicity and are not optimized at all for performance.

Reconstruction parameters for each MATLAB program are defined at the beginning of the program using the structure param. Retrospective undersampling can be defined using the variable param.undersampleSpokesTo. To increase calculation speed, the radial data is downsampled from 2x oversampling prior to the reconstruction if param.readoutDownsampling is set to 1. At the end of the calculation, the iterative reconstruction will be shown in comparison to a gridding solution. The estimated coil profiles will also be shown. The Shepp-Logan phantom program additionally shows the Fourier transform of both reconstructions, which allows verifying that applying the TV constraint in image space leads to compensation of data gaps in k-space.

If the source code is used for reconstruction of other datasets than the provided files (the MATLAB program allows reading files from Siemens MR systems with software version VBxx / VDxx), the penalty weights (lambdaxxx) and stopping criteria for the optimizer (stopTolxxx) must be adjusted accordingly.

Additional comments on the implementation can be found in the source code and in the related publication.


The source code uses the following external packages:

Version 12.10.15.


Questions about this resource may be directed to Kai Tobias Block, PhD.