ECMath CH12: Advanced Magnetic Resonance Imaging: Fingerprinting and Geometric Quantification

Project head Prof. Dr. Michael Hintermüller (1,2)
Staff Dr. Guozhi Dong (1)
Project associate Dr. Kostas Papafitsoros (2)
Project period 1 June 2017 - 31 December 2018
Affiliations (1) Humboldt-Universität zu Berlin
(2) Weierstrass Institute Berlin



x x

Background of the project

Magnetic resonance imaging

Magnetic resonance imaging (MRI) is a powerful imaging tool for a detailed visualization of soft tissues, leading to accurate diagnoses and thus to more effective therapies. It is based on the interaction of an externally controlled magnetic field B and the magnetization field $m$ of the tissue, both related by the Bloch equations [Wri97]
\begin{equation} \frac{\partial m}{\partial t}(t)=m(t)\times \gamma B(t)-\left (\frac{m_{x}(t)}{T_{2}}, \frac{m_{y}(t)}{T_{2}}, \frac{m_{z}-m_{eq}}{T_{1}} \right)^{\top}.\qquad\qquad \qquad (1) \end{equation} Here, $\gamma$ and $m_{eq}$ are known parameters while $T_{1}$, $T_{2}$ are the longitudinal and transverse relaxations times respectively which depend on the tissue type. Thus, precise estimations of these quantities are necessary for its accurate description. Focusing on a specific tissue slice, the transverse magnetization component ($m_{x}$, $m_{y}$) is excited with radio frequency pulses, and on its way to equilibrium it precedes around the $z$-axis, emitting an electromagnetic signal which is then detected. The obtained data $D$, consist of an incomplete collection (by applying a projection $P$) of the coefficients of the Fourier transform $\mathcal{F}$ of the net magnetization $M=\rho m$ where $\rho$ is the matter density, i.e., $D=P\mathcal{F} M$, $M=\rho m$ where $m$ solves (1). The reconstruction of $M$ requires "inverting" the Fourier transform (in the spirit of an inverse problem). Due to undersampling, compressed sensing techniques have been recently employed for this final reconstruction step. In classical MRI, multiple excitation pulses are used to obtain a sufficiently large sample of coefficients for a reliable reconstruction. This typically results in long scan times as, after each excitation, sufficient time must be allowed for the magnetization to achieve equilibrium again.

Magnetic resonance fingerprinting

Magnetic resonance fingerprinting (MRF) has been recently introduced as a highly promising MRI acquisition scheme which allows for the simultaneous quantification of the tissue parameters (e.g. $T_{1}$, $T_{2}$ and others) using a single acquisition process [MGS+13]. In MRF, the tissue of interest is excited through a random sequence of pulses without needing to wait for the system to return to equilibrium between pulses. After each pulse, a subset of the signal's Fourier coefficients is collected, as in classical MRI, and a reconstruction of the net magnetization image is performed. These reconstructions suffer from the presence of artifacts since the Fourier coefficients are not fully sampled. The formed sequence of image elements is then compared to a family of predicted sequences (dictionary of fingerprints) each of which corresponds to a specific combination of values of the tissue parameters. This dictionary, otherwise called Bloch manifold, is computed beforehand by solving the Bloch equations. The idea is that, provided the dictionary is rich enough, every material element (voxel) can be then mapped to its parameter values. While first, very promising results have been obtained in biomedical engineering, many aspects of MRF remain widely open and require a proper mathematization for optimizing and robustifying the procedure.

Aim of the project

Physically based quantitative mathematical model for MRF

MRF research is very recent in biomedical engineering. Like for many other imaging modalities, the further advance of the technology will strongly depend on its successful mathematical treatment. In the context of MRF this starts already with setting up a proper mathematical model as no such model is known, except for the one in [DPVW14] (BLIP method) conceived in a specific discrete context only. Hence, our aim is to establish a novel physically oriented framework in function space, which directly incorporates the Bloch equations into the reconstruction process. We are also interested in incorporating regularization in the overall procedure, with the first choices being the (weighted) Total variation (TV) as well as Total Generalised Variation (TGV), a higher-order, artifact-free, generalization of TV.

Variational regularization: Fast solvers and multimodal medical imaging

Any type of regularization procedure must be performed efficiently due to the large amount of involved data. Currently, primal-dual first-order and variable splitting schemes are frequently employed for TV and TGV. We refer for instance to [CP11] for TV and analogous methods for TGV. It is known that these solvers typically admit no function space analysis, hence, resulting in resolution (mesh) dependent convergence behavior; see [HRH14] for an investigation. Moreover, proper discretization must be employ to keep numerical artifacts (such as numerical diffusion near contours) minimal. In this vein, this project will develop quasi-second order methods by pursuing a Fenchel-dual approach along with an advancing semismooth Newton methods for the resulting dual variational problems; see [HK04] for TV regularization with scalar regularization parameter.

With respect to variational regularization methods, we are also interested in multimodality in medical imaging, i.e., where, for instance, information from one modality, e.g., an MRI image, can be exploited in the reconstruction process of another modality, e.g., Positron Emission Tomography (PET). One of the targets of the project will be a to develop a rigorous mathematical framework for TV-based regularization functionals that are involved in these kind of tasks.

Geometric quantification in segmentation

Finally, we aim to investigate the uncertainty of MR-data based segmentation (not necessarily only related to MRF), given data along a known distribution of measurement errors. Concerning the latter, we assume that it can be represented as a finite sum of i.i.d. normally distributed random variables. We are in particular interested in a proper visualisation of the uncertainty in the segmentation in order to facilitate a decision support devise for medical practitioners. The underlying mathematical segmentation paradigm is the Mumford-Shah functional which contains the unknown boundaries of image segments as geometric quantities [HR04]. For the quantification of the geometric uncertainty, i.e., the uncertainty in segmentation, using the renowned Ambrosio-Tortorelli phase field approximation, the original sharp interface (contour) model will transformed into a purely functional context. The quantification will then be done with respect the posteriori related to the phase field variable using a variance reduction approach based on anthitetic sampling within a Monte Carlo framework.

Progress of the project

A novel physically oriented method for quantitative MRI

In our work [1], we performed an analysis of the MRF and BLIP algorithms in an inverse problems setting. Motivated by the fact that the Bloch manifold is non-convex, and the accuracy of the MRF type algorithms is limited by the discretization size of the dictionary, we propose here a novel physically oriented method for qMR, where a single step, dictionary-free model for estimating the values of the tissue parameters was proposed. It relies on a non- linear operator equation instead of the conventional two step models (comprised of (i) reconstruction of the magnetization and (ii) matching to dictionary). In a compact form, this operator equation reads
\begin{equation} P\mathcal{F}(\rho T_{x,y}M(\theta))=D, \qquad\qquad \qquad (2) \end{equation} where $D$ is the detected MRI signal, $\theta=T_{1},T_{2}$, $T_{x,y}M(\theta)$ is the solution of the Bloch equations corresponding to $\theta$, projected in the transverse domain, and $P\mathcal{F}$ is a subsampling of Fourier coefficients. After showing some differentiability properties for the operator in (2), we proposed a Levenberg-Marquardt method for the solution of the non-linear equation. Stability under noise and subsampling were shown and verified via numerical examples. In contrast to the MRF and BLIP algorithms, the performance of our proposed method is not restricted by the fineness of a dictionary and it outperforms them in terms of accuracy, memory and computational time.

                       Pointwise error of BLIP method [DPVW14] for the estimation of tissue parameters $T_{1}$, $T_{2}$ and density $\rho$.

                       Pointwise error of our method for the estimation of tissue parameters $T_{1}$, $T_{2}$ and density $\rho$.
                       Note the different scale between the two rows.

Geometric quantification in segmentation

In [4], we applied a Bayesian technique to image segmentation in order to obtain a better understanding of the behaviour of edges with respect to certain error types and to provide mathematical meaning to the informal notion of "random edges". Our approach is based on the Ambrosio-Tortorelli approximation of the Mumford-Shah segmentation model and allows to detect areas where it is likely to encounter an edge under the presence of uncertainty. The method was successfully applied to the segmentation of MRI images corrupted by motion blur.

             Left: Original image; Middle; Segmentation result for the noise model; Right: result for the motion blur model

Multimodal medical imaging: A function space framework for structural total variation regularization

In [5], we introduced and analyzed a function space framework for a large class of structural TV functionals that are typically used in multimodal medical imaging. We showed that a function space formulation of associated functionals can be well-defined through a relaxation process and by using a rather general duality result. One of our key outcomes guarantees that the Tikhonov regularization problem in function space can still be understood via its equivalence to a saddle-point formulation, where no knowledge of the precise formulation of this relaxation is needed. As a result, our work allows the function space formulation of a wide class of multimodal medical imaging problems, as MR-guided PET reconstruction shown below. This is important as it paves the path for the development of corresponding function space solution algorithms, where resolution independent behaviour can be guaranteed.

             MR-guided PET reconstruction using a structural TV functional. From left to right: (i) MRI prior, (ii) PET ground truth,
            (iii) Standard TV, (iv) Structural TV

Bilevel optimization for optimal selection of regularization parameters

In [7] and [8], bilevel optimization for the optimal selection of spatially dependent total variation regularization parameters was finalized. It was successfully applied to the MRI reconstruction from undersampled data with emphasis on the recovery of fine scale details. In an additional work [6], surrogate iterative methods were employed to handle given subsampled data in transformed domains, such as Fourier or wavelet data, leading to improved MRI reconstructions. This work is currently being extended to optimal selection of spatially varying parameters in TGV, and is further supported by an upcoming publication on fast dual-based TGV denoising solver. The dualization framework was established via a series of density of convex intersection results in Banach spaces in [2] and [9].

Further activities: June 2017-today


Communication of the above results to the community as well as further build-up of synergies was done in international conferences and workshops. From August 29 to December 20, 2017, we co-organized the four-month research program "Variational Methods and Effective Algorithms for Imaging and Vision" at the Isaac Newton Institute for Mathematical Sciences, University of Cambridge, U.K., where also an invited talk was given. Furthermore, supported by the ECMath and the DFG, we organized the latest version of the "Mathematics and Image Analysis" (MIA 2018) conference, January 15-17, 2018, in Berlin, where we also presented three posters. We actively participated in the last SIAM conference on Imaging Science, June 5-8, 2018, Bologna, Italy, with an organization of a minisymposium on "Learning and Adaptive Approaches in Image Processing" and three invited talks.

Internal MATHEON cooperation

An internal MATHEON cooperation with Jörg Polzehl and Karsten Tabelow (head of ECMath project OT7, was initiated in [3], where a new image denoising method, called patch-wise adaptive smoothing (PAWS), was developed. There, local patches of image intensities are compared in order to define local adaptive weighting schemes, that preserve sharp edges in the images while reducing the reconstruction bias. Connections and analogies to the relationship between TV and TGV were drawn and the method was successfully tested for the denoising of 3D MRI image volumes.


[1] G. Dong, M. Hintermüller, K. Papafitsoros, A physically oriented method for quantitative magnetic resonance imaging, WIAS preprint 2528, (2018).
[2] M. Hintermüller, K. Papafitsoros, C.N. Rautenberg, Singular mollifiers and applications, HU preprint, (2018).
[3] J. Polzehl, K. Papafitsoros, K. Tabelow, Patch-wise adaptive weights smoothing, WIAS preprint 2520, (2018).


[4] S.-M. Stengl, Bildsegmentierung im Mumford-Shah-Modell: Analysis und Numerik mit anschließender Quantifizierung von Unsicherheiten für fehlerhafte Bilder, Master thesis, Humboldt-Universität zu Berlin, (2017).

Accepted Publications

[5] M. Hintermüller, M. Holler, K. Papafitsoros, A function space framework for structural total variation regularization with applications in inverse problems, Inverse Problems, vol. 34, no. 6, 064002, (2018).
[6] M. Hintermüller, A. Langer, C.N. Rautenberg, T. Wu, Adaptive regularization for image reconstruction from subsampled data, to appear in, X.-C. Tai, E. Bae and M. Lysaker, editors, Imaging, Vision and Learning Based on Optimization and PDEs. IVLOPDE, Bergen, Norway, August 29 - September 2, 2016, Springer International Publishing, (2018).
[7] M. Hintermüller, C.N. Rautenberg, Optimal selection of the regularization function in a weighted total variation model. Part I: Modelling and Theory. Journal of Mathematical Imaging and Vision 59, no. 3, 498-514, (2017).
[8] M. Hintermüller, C.N. Rautenberg, T. Wu, A. Langer, Optimal selection of the regularization function in a weighted total variation model. Part II: Algorithm, its analysis and numerical tests. Journal of Mathematical Imaging and Vision 59, no. 3, 515-533, (2017).
[9] M. Hintermüller, C.N. Rautenberg, S. Rösel, Density of convex intersections and applications. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 473, no. 2205, (2017).


[MGS+13] D. Ma, V. Gulani, N. Seiberlich, K. Liu, J.L. Sunshine, J.L. Duerk, M.A. Griswold, Magnetic resonance fingerprinting. Nature, 495(7440):187-192, 2013.
[Wri97] G.A. Wright, Magnetic resonance imaging. IEEE Signal Processing Magazine, 14(1):56-66, 1997.
[DPVW14] M. Davies, G. Puy, P. Vandergheynst, Y. Wiaux, A compressed sensing framework for magnetic resonance fingerprinting. SIAM Journal on Imaging Sciences, 7(4):2623-2656, 2014.
[CP11] A. Chambolle, T. Pock, A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision, 40(1):120-145, 2011.
[HRH14] M. Hintermüller, C.N. Rautenberg, J. Hahn, Functional-analytic and numerical issues in splitting methods for total variation-based image reconstruction. Inverse Problems, 30(5):055014, 2014.
[HK04] M. Hintermüller, K. Kunisch, Total bounded variation regularization as a bilaterally constrained optimization problem. SIAM Journal of Applied Mathematics, 64(4):1311-1333, 2004.
[HK15] M. Hintermüller, C.N. Rautenberg, On the density of classes of closed convex sets with pointwise constraints in Sobolev spaces. Journal of Mathematical Analysis and Applications, 426(1):585-593, 2015.
[HR04] M. Hintermüller, W. Ring, An inexact Newton-CG-type active contour approach for the minimization of the Mumford-Shah functional. Journal of Mathematical Imaging and Vision, 20:19-42, 2004.