Reconstruction

The main reconstruction tool you should use is `riesling recon-rlsq`_. This solves the regularized least-squares reconstruction problem using the Alternating-Directions Method-of-Multipliers algorithm. It supports several common regularizers including L1-wavelets, Total Variation and Total Generalized Variation. To run an unregularized reconstruction use `riesling recon-lsq`_. recon-lsq and recon-rlsq will self-calibrate sensitivity maps from the input data if no maps are supplied. These can be calculated explicitly with the sense-calib command, see there for details of the calibration.

Common Options

recon-lsq

Reconstructs the image as the least-squares solution to an inverse problem. The system matrix includes sensitivity maps and the NUFFT. The LSMR algorithm is used as this does not require forming the normal equations, keeping the condition number low. In combination with pre-conditioning in k-space this ensures convergence in a handful of iterations.

Usage

riesling recon-lsq input.h5 output.h5

Important Options

  • --max-its=N, --atol=A, --btol=B, --ctol=C

    Termination conditions. Reasoable image quality can be obtained in less than five iterations. The a and b tolerances are relative to how accurate the solution has become, c is a tolerance on the condition number of the system.

  • --lambda=L

    Apply basic Tikohonov/L2 regularization to the reconstruction.

  • --tacq=T, --Nt=N

    Off-resonance correction parameters. To use these you must also supply an off-resonance map in the input. This should be a dataset labelled f0map and must be the same matrix size as your image at the reconstruction FOV. These two parameters are the total acquisition time and the number of time segments you want for the reconstruction. The units of tacq must match f0map, i.e. if your off-resonance map is in Hz then your acquisition time is in seconds.

recon-rlsq

By default, uses the Alternating Directions Method-of-Multipliers (ADMM) to add regularizers to the least-squares reconstruction problem. This is similar to the BART pics command. See S. Boyd, ‘Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers’ doi: 10.1561/2200000016. The Primal Dual Hybrid Gradient (PDHG) algorithm is also implemented and faster than ADMM, but it requires the calculation of the maximum eigenvalue of the encoding operator (using riesling eig).

Usage

riesling recon-rlsq input.h5 output.h5 --tgv=1e-3

Regularizers

See denoise. The same regularizers are available for recon-rlsq.

Common Options

  • --scale=bart/otsu/S

    The optimal regularization strength λ depends both on the particular regularizer and the typical intensity values in the unregularized image. To make values of λ roughly comparable, it is usual to scale the data such that the intensity values are approximately 1 during the optimization (and then unscale the final image). By default riesling will perform a NUFFT and then use Otsu’s method to find the median foreground intensity as the scaling factor (specify otsu to make this explicit). The BART automatic scaling can be chosen with bart. Alternately a fixed numeric multiplicative scaling factor can be specified, which will skip the initial NUFFT. If you already know the approximate scaling of your data (from a test recon), this option will be the fastest.

ADMM Options

The Alternating-Directions Method-of-Multipliers is a very robust algorithm for solving non-smooth regularized least-squares. However, it requires solving an inverse problem on every iteration, which itself must be solved using an iterative scheme. This means it can be very slow. However, there is only parameter for the algorithm (rho) and the adaptive scheme implemented in riesling means that you should not have to adjust the default parameter. ADMM is hence currently the default choice as it is essentially guaranteed to converge to a sensible answer, given enough iterations.

  • --max-its1=N, --max-its0=N–atol=A``, --btol=B, --ctol=C

    These control the inner optimization for ADMM (the x update step), which is solved with LSMR. As this step is warm-started, it is possible to control the maximum number of iterations for the zeroth and subsequent steps independently.

  • --max-outer-its=N

    The maximum number of ADMM iterations. The default is 20 but a higher number (50 or more) may be required for optimal image quality.

  • --eps=E

    Primal and dual convergence tolerance for ADMM. Default value is 0.01.

  • --rho=P

    Coupling factor for ADMM. The default value of 1 is robust, and will be adjusted inside the algorithm according to ADMM Penalty Parameter Selection by Residual Balancing.

PDHG Options

The preconditioned Primal-Dual Hybrid Gradient is potentially must faster than ADMM as it does not require an inner solve. However, if the step lengths are incorrectly chosen it will not converge.

recon-rss

Perform a basic reconstruction using root-sum-of-squares channel combination. Very fast but worst image quality. Does not calculate or use sensitivity maps. Useful for testing.

Usage

riesling recon-rss input.h5 output.h5

sense-calib

Sensitivity maps are an integral part of any reconstruction from a multi-channel coil. Calculating high quality sensitivity maps is a difficult and open research question for two reasons. First, the multi-channel reconstruction problem is ill-posed as there is no unique solution (if the sensitivities are multiplied and the image divided by an arbitrary field the same data will result), and second because sensitivities exist in the background of the image where we cannot acquire signal.

RIESLING estimates sensitivities assuming that a fully-sampled calibration region with consistent contrast has been acquired in the data. This is true for the majority of non-cartesian sequences, see E. N. Yeh et al., ‘Inherently self-calibrating non-cartesian parallel imaging’, Magnetic Resonance in Medicine, vol. 54, no. 1, pp. 1–8, Jul. 2005,, and this step is hence incorporated into the reconstruction commands. However, there are many situations where it is beneficial to calculate the sensitivities up-front, potentially from alternate data. There is hence an explicit sense-calib command for this. All the relevant options to this command are also exposed for the reconstruction commands.

Note that RIESLING calculates and stores the sensitivity kernels in k-space, not the maps themselves. If you want to see the maps, a separate sense-maps command is provided to convert between them.

The FOV and oversampling used in the calibration must match your reconstruction.

The algorithm used by RIESLING consists of these steps:

  1. Reconstruct low-resolution images for each channel from a fully-sampled calibration region (inverse NUFFT).

  2. Obtain a reference image either from fully-sampled single channel data or by taking the root-sum-squares across the multi-channel images.

  3. Solve the inverse problem \(c = RFPs\) where \(s\) are the sensitivity kernels in k-space, \(c\) are the channel images, \(P\) is a padding operator, \(F\) is an FT, and \(R\) is an operator that multiplies each channel by the reference image.

To ensure the maps are smooth and have support in the background region, the forward model is modified to incorporate regularization with a Sobolev Norm term \(λW = (1 + |k|^2)^{l/2}\) (where \(k\) is k-space co-ordinate, i.e. penalises high frequency terms) and a mask \(M\) over the object:

\[\begin{split}c' = A s\\ c' = \begin{bmatrix} c\\ 0 \end{bmatrix}\\ A = \begin{bmatrix} M R F P\\ λW \end{bmatrix}\end{split}\]

This problem is badly conditioned, and even with a preconditioner can take approximately 100 iterations to converge. However due to the small matrix sizes this should only take a few seconds. See H. C. M. Holme, S. Rosenzweig, F. Ong, R. N. Wilke, M. Lustig, and M. Uecker, ‘ENLIVE: An Efficient Nonlinear Method for Calibrationless and Robust Parallel Imaging’, Scientific Reports, vol. 9, no. 1, Dec. 2019, for the regularizer.

Usage

riesling sense-calib input.h5 kernels.h5

Important Options

  • --ref=reference.h5

    Use the supplied data to reconstruct the reference image (i.e. from a body coil acquisition) instead of using the root-sum-squares of the channels.

  • --sense-lambda=λ

    The amount of regularization to apply to the sensitivities. Over regularization will result in the per-voxel sensitivities reducing.

  • --sense-l=L

    The L parameter to the Sobolev Norm weights. Higher numbers increase the regularization strength in a highly non-linear fashion.

  • --sense-res=R

    The resolution of the initial reconstructions for the sensitivity maps. Because sensitivities are generally agreed to be smooth, only a low resolution reconstruction is required and the default is 6mm isotropic. However, the resulting images must have a sufficiently large matrix size to extract the kernels from.

  • --sense-width=K

    The width of the sensitivity kernels in k-space on the nominal grid. The value specified here will be mulitipled by the oversampling factor to produce the final kernel size. Hence, if you override the default oversampling in the main reconstruction you must also do so here.

  • --sense-tp=T

    If the input data contains multiple timepoints, use this one to calculate the sensitivities (default is first volume).

denoise

Denoise an image using a proximal operator possibly combined with the Primal Dual Hybrid Gradient algorithm. See F. Ong, Accelerating Non-Cartesian MRI Reconstruction Convergence Using k-Space Preconditioning’ doi: 10.1109/TMI.2019.2954121

Usage

riesling denoise input.h5 output.h5 --tgv=1e-3

Important Options

  • --max-its=N, --res-tol=r, --delta-tol=d, --lambda-G=l

    These control the PDHG algorithm for regularizers that cannot be solved using only a proximal operator (i.e. those that have a non-invertible transform). The residual tolerance is calculated relative to the initial data, and the delta-tolerance will terminate the algorithm when the x update becomes small enough. --lambda-G is the maximum eigenvalue for the regularizer transform and is used to determine the step-size. The default value of 16 is sufficiently large to cover all the regularizers implemented in riesling.

  • --scale=bart/otsu/S

    The optimal regularization strength λ depends both on the particular regularizer and the typical intensity values in the unregularized image. To make values of λ roughly comparable, it is usual to scale the data such that the intensity values are approximately 1 during the optimization (and then unscale the final image). By default riesling will perform a NUFFT and then use Otsu’s method to find the median foreground intensity as the scaling factor (specify otsu to make this explicit). The BART automatic scaling can be chosen with bart. Alternately a fixed numeric multiplicative scaling factor can be specified, which will skip the initial NUFFT. If you already know the approximate scaling of your data (from a test recon), this option will be the fastest.

Regularization Options

Multiple regularizers can be specified simultaneously, each with a different regularization strength λ and options. At least one regularizer must be specified, there is no default option at present.