Biomedical Image Reconstruction: From the Foundations to Deep Neural Networks

TL;DR

This tutorial reviews the evolution of biomedical image reconstruction methods from traditional techniques to modern deep learning approaches, highlighting common principles and the integration of different generations.

Contribution

It unifies diverse research on biomedical image reconstruction, clarifies connections between classical, sparsity-based, and learning-based methods, and provides a comprehensive framework for understanding their development.

Findings

Classical methods serve as modules in modern algorithms

Deep learning approaches outperform traditional techniques

Unified framework aids in designing advanced reconstruction algorithms

Abstract

This tutorial covers biomedical image reconstruction, from the foundational concepts of system modeling and direct reconstruction to modern sparsity and learning-based approaches. Imaging is a critical tool in biological research and medicine, and most imaging systems necessarily use an image-reconstruction algorithm to create an image; the design of these algorithms has been a topic of research since at least the 1960's. In the last few years, machine learning-based approaches have shown impressive performance on image reconstruction problems, triggering a wave of enthusiasm and creativity around the paradigm of learning. Our goal is to unify this body of research, identifying common principles and reusable building blocks across decades and among diverse imaging modalities. We first describe system modeling, emphasizing how a few building blocks can be used to describe a broad range of imaging modalities. We then discuss reconstruction algorithms, grouping them into three broad generations. The first are the classical direct methods, including Tikhonov regularization; the second are the variational methods based on sparsity and the theory of compressive sensing; and the third are the learning-based (also called data-driven) methods, especially those using deep convolutional neural networks. There are strong links between these generations: classical (first-generation) methods appear as modules inside the latter two, and the former two are used to inspire new designs for learning-based (third-generation) methods. As a result, a solid understanding of all of three generations is necessary for the design of state-of-the-art algorithms.

Figures

Tables

Table 1. Table 4.1 : Infinitely divisible distributions and their corresponding potential functions. B 𝐵 B is the beta function.

	$p_U\left(u\right)$	${\Phi }_U\left(u\right)$
Gaussian	$\frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{\frac{-u^2}{2{\sigma }^2}}$	$\frac{1}{2{\sigma }^2}{u}^2+C$
Laplace	$\frac{\lambda }{2}{e}^{-\lambda \left|u\right|}$	$\lambda \left|u\right|+C$
Student	$\frac{1}{B(r,\frac{1}{2})}{\left(\frac{1}{u^2+1}\right)}^{r+\frac{1}{2}}$	$\left(r+\frac{1}{2}\right)\log \left(1+u^2\right)+C$

Equations

$$\mathbf{g}=\mathcal{H}\{f\}\in {\mathbb{R}}^M.$$

$$\mathcal{H}\{f\}\approx H\{f\}+n.$$

$$H\{af+bg\}=aH\{f\}+bH\{g\},$$

$$[H\{f\}{]}_m=[\mathbf{g}{]}_m=⟨{\eta }_m,f⟩={\int }_{{\mathbb{R}}^d}{\eta }_m(\mathbf{x})f(\mathbf{x})d\mathbf{x},$$

$$[S_X\{f\}{]}_m=f({\mathbf{x}}_m).$$

$$\mathcal{F}\{f\}(\mathbf{\omega })=\hat{f}(\mathbf{\omega })={\int }_{{\mathbb{R}}^d}f(\mathbf{x}){\text{e}}^{-\text{j}⟨\mathbf{\omega },\mathbf{x}⟩}d\mathbf{x}.$$

$$\mathbf{g}=(S_X\circ \mathcal{F})\{f\},$$

$$f(\mathbf{x})={\mathcal{F}}^{-1}\{\hat{f}\}(\mathbf{x})=\frac{1}{(2\pi {)}^d}{\int }_{{\mathbb{R}}^d}\hat{f}(\mathbf{\omega }){\text{e}}^{\text{j}⟨\mathbf{\omega },\mathbf{x}⟩}d\mathbf{\omega }.$$

$$M_h\{f\}(\mathbf{x})=f(\mathbf{x})h(\mathbf{x}),$$

$$\mathrm{rect}(\mathbf{x})=\{\begin{array}{ll}1& \mathbf{x}\in [-\frac{1}{2},\frac{1}{2}{]}^d;\\ 0& \text{otherwise}.\end{array}$$

$$h(\mathbf{x})={\text{e}}^{\text{j}⟨{\mathbf{\omega }}_0,\mathbf{x}⟩}.$$

$$(\mathcal{F}\circ M_{{\text{e}}^{\text{j}⟨{\mathbf{\omega }}_0,\mathbf{x}⟩}})\{f\}(\mathbf{\omega })=\mathcal{F}\{f\}(\mathbf{\omega }-{\mathbf{\omega }}_0).$$

$$\mathbf{g}=(S_X\circ \mathcal{F}\circ M_h)\{f\},$$

$$C_h\{f\}(\mathbf{x})={\int }_{{\mathbb{R}}^d}h(\mathbf{x}-{\mathbf{x}}^{\prime })f({\mathbf{x}}^{\prime })d{\mathbf{x}}^{\prime }.$$

$$C_h\{f\}(\mathbf{x})=({\mathcal{F}}^{-1}\circ M_{\hat{h}}\circ \mathcal{F})\{f\}(\mathbf{x}),$$

$$g=C_h\{f\}.$$

$$h(\mathbf{x})={\left(\frac{2J_1(r)}{r}\right)}^2,r=\frac{2\pi a\Vert \mathbf{x}{\Vert }_2}{\lambda R},$$

$$\hat{h}(\mathbf{w})=\{\begin{array}{ll}C,& \Vert \mathbf{\omega }{\Vert }_2\le w_{\text{max}};\\ 0,& \text{otherwise}.\end{array}$$

$$\mathcal{F}g$$

$$g=(C_h\circ M_p)\{f\},$$

$$\begin{gathered} \mathcal{F}\{g\}(\mathbf{\omega })=(M_{\hat{h}}\circ \mathcal{F}\circ M_p)\{f\}(\mathbf{\omega }) \\ =\hat{h}(\mathbf{\omega })\mathcal{F}\{f\}(\mathbf{\omega }-{\mathbf{\omega }}_0). \end{gathered}$$

$$f(\mathbf{x}-{\mathbf{x}}_0)={\Phi }_{{\phi }_{{\mathbf{x}}_0}}\{f\}(\mathbf{x}),\text{with}{\phi }_{{\mathbf{x}}_0}(\mathbf{x})=\mathbf{x}-{\mathbf{x}}_0.$$

$$X_{\theta }\{f\}(y)={\int }_{\mathbb{R}}f(y\mathbf{\theta }+t{\mathbf{\theta }}^{\perp })dt,$$

$$({\mathcal{F}}_{\text{1D}}\circ X_{\theta })\{f\}(\omega )=({\Phi }_{{\phi }_{\theta }}\circ {\mathcal{F}}_{\text{2D}})\{f\}(\omega ),$$

$$({\Phi }_{{\phi }_{\theta }}\circ {\mathcal{F}}_{\text{2D}})\{f\}(\omega )$$

$$={\int }_{{\mathbb{R}}^2}f(\mathbf{A}\mathbf{a}){\text{e}}^{-\text{j}⟨\omega \theta ,\mathbf{A}\mathbf{a}⟩}d\mathbf{a}$$

$$={\int }_{{\mathbb{R}}^2}f(\mathbf{A}\mathbf{a}){\text{e}}^{-\text{j}\omega [\mathbf{a}{]}_1}d\mathbf{a}$$

$$={\int }_{\mathbb{R}}{\int }_{\mathbb{R}}f(\mathbf{A}\mathbf{a})d[\mathbf{a}{]}_2{\text{e}}^{-\text{j}\omega [\mathbf{a}{]}_1}d[\mathbf{a}{]}_1$$

$$={\mathcal{F}}_{\text{1D}}\circ X_{\theta }f(\omega ),$$

$$f(\mathbf{x})=T_{\{{\beta }_n{\}}_{n=1}^N}\{\mathbf{f}\}(\mathbf{x})=\sum _{n=1}^N[\mathbf{f}{]}_n{\beta }_n(\mathbf{x}),$$

$$[\mathbf{H}{]}_{m,n}=⟨{\eta }_m,{\beta }_n⟩={\int }_{{\mathbb{R}}^d}{\eta }_m(\mathbf{x}){\beta }_n(\mathbf{x})d\mathbf{x}.$$

$$\beta (\mathbf{x})=\beta (x_1)\beta (x_2)\cdots \beta (x_d),$$

Full text

\newunicodechar

ℓ$\ell$ \newunicodechar₁1 \newunicodechar₂2

Biomedical Image Reconstruction: From the Foundations to Deep Neural Networks

Michael T. McCann

  

Michael Unser

(École polytechnique fédérale de Lausanne)

Abstract

This tutorial covers biomedical image reconstruction, from the foundational concepts of system modeling and direct reconstruction to modern sparsity and learning-based approaches.

Imaging is a critical tool in biological research and medicine, and most imaging systems necessarily use an image-reconstruction algorithm to create an image; the design of these algorithms has been a topic of research since at least the 1960’s. In the last few years, machine learning-based approaches have shown impressive performance on image reconstruction problems, triggering a wave of enthusiasm and creativity around the paradigm of learning. Our goal is to unify this body of research, identifying common principles and reusable building blocks across decades and among diverse imaging modalities.

We first describe system modeling, emphasizing how a few building blocks can be used to describe a broad range of imaging modalities. We then discuss reconstruction algorithms, grouping them into three broad generations. The first are the classical direct methods, including Tikhonov regularization; the second are the variational methods based on sparsity and the theory of compressive sensing; and the third are the learning-based (also called data-driven) methods, especially those using deep convolutional neural networks. There are strong links between these generations: classical (first-generation) methods appear as modules inside the latter two, and the former two are used to inspire new designs for learning-based (third-generation) methods. As a result, a solid understanding of all three generations is necessary for the design of state-of-the-art algorithms.

List of Abbreviations

ADMM alternating direction method of multipliers CCD charge-coupled device CG conjugate gradient CNN convolutional neural network CT computed tomography DCT discrete cosine transform ET electron tomography FBP filtered back projection FFT fast Fourier transform GPU graphics processing unit i.i.d. independent and identically distributed ISTA iterative shrinkage and thresholding MAP maximum a posteriori MMSE minimum mean square error MRI magnetic resonance imaging MSE mean squared error PDF probability distribution function PET positron emission tomography PSF point spread function RKHS reproducing kernel Hilbert spaces SGD stochastic gradient descent SIM structured-illumination microscopy SNR signal-to-noise ratio SPECT single-photon emission computed tomography SSIM structural similarity index TCIA The Cancer Imaging Archive TV total variation USC-SIPI University of Southern California Signal and Image Processing Institute

Chapter 1 Introduction

Biomedical imaging is a vast and diverse field: there are a plethora of imaging devices using, e.g., light, X-rays, sound waves, magnetic fields, electrons, or protons, to measure structures ranging from nano- to macroscale. In many cases, computer software is needed to turn the signals collected by the hardware into a meaningful image. These computer algorithms are similarly diverse and numerous. For example, Google Scholar lists around 1,600 papers with the words “MRI reconstruction” in the title.

In this tutorial, we aim to present a wide swath of biomedical image reconstruction algorithms under a single framework—using a unified notation across domains, modeling many modalities with a few basic operators, and emphasizing commonalities among reconstruction algorithms. Our goal is not to review the totality of the biomedical image reconstruction literature, nor even to provide a comprehensive tutorial (though we have included references throughout that we hope will be useful pointers for a reader wishing to learn more). Instead, we have focused on creating a brief and coherent narrative that traverses some six decades of research. The result is naturally a little stylized, and reflects the biases of the authors. For another recent perspective on these topics, see [1].

Focusing on the similarities across modalities is not only a didactic choice, it is also a practical one. While developing high-quality reconstruction algorithms requires a deep understanding of both the physics of the imaging system and the biomedical questions at hand, we have found great benefit in our own work from adopting a unified point of view. We model imaging modalities as combinations of a small set of building blocks, which allows us to see connections between modalities and quickly port ideas and computer code from one to the next. Reconstruction algorithms can treat the imaging model as a black box, meaning that one algorithm can work for many modalities. Such modular solutions may not be completely optimal, but they are often sufficient for a given application. When they are not, they still serve as a valuable starting point.

The structure of the tutorial is as follows. We begin by introducing the concept of forward models and a set of building blocks for creating them in Section 2. In the next three sections we introduce reconstruction algorithms in a roughly chronological way, beginning with direct reconstruction and $\ell_{2}$ regularization (Section 3), followed by modern, sparsity-based techniques (Section 4), and ending with learning-based methods (Section 5). We compare these approaches and give our outlook on the next steps for the field in Section 6.

Chapter 2 Forward Models

The goal of any imaging system is recover an image, which is a map of a meaningful physical quantity that varies over space. In biomedical imaging, this might be the concentration of fluorophores in a cell or the density of the tissues in a person’s body. Mathematically, we describe an image as a function, $f:\mathbb{R}^{d}\to\mathbb{R}$, that maps a point in space ${\boldsymbol{x}}\in\mathbb{R}^{d}$ to a real number, $y=f({\boldsymbol{x}})\in\mathbb{R}$. We leave the dimension, $d$, of ${\boldsymbol{x}}$ unspecified; it is usually between one and four: up to three spatial dimensions possibly plus time.

In any practical setting, the image of interest is not directly observed, but must be reconstructed from measurements. We view $f$ through a physical imaging system, $\mathcal{H}$, that maps $f$ to a vector of measurements,

[TABLE]

The measurement process encompasses all of the physics of the imaging system, including the sources of noise. It is nondeterministic, nonlinear, and impossible to perfectly characterize or simulate. The measurements, on the other hand, are necessarily finite, quantized, and discrete because they are stored digitally.111Throughout the text, we use bold font for these finite-length vectors and matrices to emphasize that they are concrete objects that we can store and use in computations. Biomedical imaging includes a broad diversity of imaging modalities, e.g., systems that rely on electromagnetic radiation (visible light, X-rays, gamma rays, etc.), electrons, magnetic fields, or sound waves; see Figure 2.1 for examples.

Most222but not all, as we will discuss in Section 5 methods to recover an image from its measurements rely on modeling the imaging system, $\mathcal{H}$, with a deterministic operator $H$ and a random noise term, $n\in\mathbb{R}^{M}$, so that

[TABLE]

Let us first focus on the deterministic part, $H\{f\}$. The key challenge is to develop an operator $H$ that is simultaneously mathematical convenient, computationally tractable, and physically accurate is a key challenge; fortunately, it turns out that many biomedical imaging modalities can be approximated using combinations of just a few basic linear operators. In this section, we will introduce these operators along with examples.

2.1 Vector Spaces

Before we can describe the action of an operator, $H$, we need to specify its domain and range, i.e., where it maps from and to. Collectively, these are known as vector spaces, which are the central topic of linear algebra. For our purposes, it suffices to put the mathematical details aside and define two special vector spaces that we will use throughout.

The first is $\mathbb{C}^{N}$, which is the space of finite-length lists of $N$ complex numbers.333In an engineering context, the term vector usually refers to an element of $\mathbb{C}^{N}$ or $\mathbb{R}^{N}$; in the following, we have always tried to specify vector spaces explicitly to avoid ambiguity. Generally, we use elements of this space to represent measurements or coefficients. The use of complex-valued vectors is mathematically convenient, but can be a source of frustration during implementation—on one hand, some programming languages do not natively handle complex values, one the other, some papers specify algorithms only for real values. When trouble arises, it is useful to remember that a vector of $N$ complex numbers can be represented by a vector of $2N$ real numbers, provided all arithmetic performed on the vector respects the convention used.

The second is the space of square-integrable functions $L_{2}(\mathbb{R}^{d})=\{f:\int_{\mathbb{R}^{d}}|f({\boldsymbol{x}})|^{2}d{\boldsymbol{x}}<\infty\}$. It is also called the space of finite-energy functions and is closely related to the concept of energy in physics. We use elements of this space to represent continuous objects, including the image we are trying to reconstruct. While the choice of $L_{2}(\mathbb{R}^{d})$ cannot always be physically justified, it is mathematically convenient and hardly restrictive on the resulting reconstructions. Other restrictions, (e.g., on the size of the non-zero support of a function, its bandwidth, or its smoothness) may also be important depending on the specifics of the imaging system; we will mention these as they arise. We can generalize the concept of $L_{2}(\mathbb{R}^{d})$ to $L_{p}(\mathbb{R}^{d})$ in a straightforward way: $L_{p}(\mathbb{R}^{d})=\{f:\int_{\mathbb{R}^{d}}|f({\boldsymbol{x}})|^{p}d{\boldsymbol{x}}<\infty\}$, with $L_{1}(\mathbb{R}^{d})$ being the most common.

2.2 Linear Operators

Our focus here is on operators that are linear.

Definition 1**.**

Let $\mathcal{X}$ and $\mathcal{Y}$ be two vector spaces. An operator, $H:\mathcal{X}\rightarrow\mathcal{Y}$, is a linear operator from $\mathcal{X}$ to $\mathcal{Y}$ when it satisfies

[TABLE]

for all scalars $a$ and $b\in\mathbb{R}$ and images $f$ and $g\in\mathcal{X}$.

Linear operators provide a good approximation for many real imaging systems. At the same time, they are mathematically convenient, e.g., they allow the reconstruction problem to be solved efficiently using convex optimization (Section 3.5). Another implication of selecting $H$ to be linear is that, when $\mathcal{X}$ is $L_{2}(\mathbb{R}^{d})$ and $\mathcal{Y}$ is $\mathbb{C}^{M}$, each measurement can be written as the inner product between the image and a detector function

[TABLE]

where $[{\boldsymbol{\mathbf{g}}}]_{m}$ designates the $m$th element of a vector ${\boldsymbol{\mathbf{g}}}$.444This result holds for certain vector spaces other than $L_{2}(\mathbb{R}^{d})$ as well; we give a special case here for the sake of simplicity. The specific form of each $\eta_{m}$ is determined by the system model, $H$. This result follows from a well-known theorem in analysis called the Riesz representation theorem [5, 6].

2.3 Building Blocks

The composition of linear operators is also linear, which means that it is useful to build up a toolbox of simple, well-understood operators from which to build more complicated ones. We begin with sampling.

Definition 2**.**

Let $\mathcal{X}$ be a vector space of continuous functions. The sampling operator, $S_{X}:\mathcal{X}\rightarrow\mathbb{R}^{M}$, with $X={\{{\boldsymbol{x}}_{m}\}_{m=1}^{M}}$, returns a vector of the values of $f$ at $M$ known locations,

[TABLE]

See Figure 2.2 for an example.

For the sampling operator to be well-defined on $\mathcal{X}$, the elements of the space need to be continuous functions, This is not the case for $L_{2}(\mathbb{R}^{d})$, however, many subspaces of $L_{2}(\mathbb{R}^{d})$ work, e.g., the space of bandlimited functions. A Hilbert space over which sampling can be defined is called @iaciRKHS reproducing kernel Hilbert spaces (RKHS), see [7] for more details.

Sampling is not a particularly useful imaging model on its own; it is useful because it allows us to express imaging as a sampled version of a well-understood continuous operation, such as the Fourier transform.

Definition 3**.**

The Fourier transform, $\mathcal{F}:L_{2}(\mathbb{R}^{d})\cap L_{1}(\mathbb{R}^{d})\rightarrow L_{2}(\mathbb{R}^{d})$, expresses $f$ in terms of its frequency components,

[TABLE]

Here, we have introduced the common shorthand that $\hat{f}$ stands for the Fourier transform of $f$. See Figure 2.3 for an example.

The restriction to $f\in L_{2}(\mathbb{R}^{d})\cap L_{1}(\mathbb{R}^{d})$ has two advantages: First, it ensures that the integral (2.6) is well defined. Second, it guarantees that $\hat{f}$ is continuous, which then makes sampling feasible.

Example 1** (Magnetic resonance imaging).**

The Fourier transform and sampling operators are sufficient to model MRI; this is known as the $k$-space formalism [9]. In proton density-weighted MRI, the physical quantity of interest is the proton density. The measured signal is given by

[TABLE]

where $\circ$ indicates operator composition,i.e., $(A\circ B)\{f\}=A\{B\{f\}\}$. The locations of the Fourier samples depend on the hardware setup. Typical choices include Cartesian sampling, random sampling, and spiral sampling.

Definition 4**.**

The inverse Fourier transform, $\mathcal{F}^{-1}:L_{2}(\mathbb{R}^{d})\rightarrow L_{2}(\mathbb{R}^{d})$, recovers $f$ from its Fourier transform,

[TABLE]

The fact that the Fourier transform is invertible suggests a method for MRI reconstruction: approximate the integral (2.8) with a sum over the known $[{\boldsymbol{\mathbf{g}}}]_{m}=\hat{f}({\boldsymbol{\omega}}_{m})$. As the measurement density increases, this method becomes more and more accurate. We will explore this approach in more detail in Section 3.1.

Definition 5**.**

The multiplication operator, $M_{h}:L_{2}(\mathbb{R}^{d})\rightarrow L_{2}(\mathbb{R}^{d})$, multiplies $f$ pointwise by another image,

[TABLE]

where $h$ must be bounded: for all ${\boldsymbol{x}}$, $|h({\boldsymbol{x}})|<\infty$.

One example of multiplication is windowing (Figure 2.4), where $h$ is a function taking values between zero and one, such as a $\operatorname{rect}$.

[TABLE]

While it is mathematically convenient to work with infinite signals (e.g., note the integral over all space in the definition of the inverse Fourier transform), real imaging systems have a limited spatial extent; windowing with $\operatorname{rect}$ conveniently models this limit.

Another example of multiplication is modulation, where $h$ is a complex exponential,

[TABLE]

The effect of modulation is easier to see in the frequency domain,

[TABLE]

Thus, multiplication by a complex exponential results in a shift in the frequency content of the image.

Example 2** (MRI with coil sensitivity).**

A more accurate model of MRI accounts for spatial variations in the coil sensitivity, which can be modeled by a of $f$ prior to Fourier sampling,

[TABLE]

where $h$ is measured as part of calibration. For a detailed treatment of the MRI forward model, see [10].

Definition 6**.**

The convolution-by-$h$ operator, $C_{h}:L_{2}(\mathbb{R}^{d})\to L_{2}(\mathbb{R}^{d})$, implements the linear, shift-invariant system with impulse response $h$,

[TABLE]

The impulse response, $h$, must be chosen so that the integral always exists; e.g., it is sufficient for $h$ to be absolutely integrable ($\int_{\mathbb{R}^{d}}|h({\boldsymbol{x}})|d{\boldsymbol{x}}<\infty$). Convolution can be equivalently defined in the frequency domain,

[TABLE]

where $\hat{h}$ is the frequency response of the system, $\hat{h}=\mathcal{F}\{h\}$.

Although convolution is a composition of Fourier transforms and a multiplication, it deserves its own symbol because of its ubiquitous role in imaging. What makes it so useful is that many measurement systems are essentially shift-invariant, because the underlying physics work the same everywhere and the sensors are placed on a regular grid, e.g., pixels on a charge-coupled device (CCD) camera. We say essentially because a shift-invariant model implicitly assumes that our measurements extend infinitely in all directions—while this is not the case, it can usually be handled properly with the help of boundary conditions. Convolution is also important because discrete convolutions can be implemented using the fast Fourier transform (FFT), which can result in fast reconstruction algorithms.

Example 3** (Brightfield microscopy).**

In brightfield microscopy, what appears at the eyepiece can be modeled as a blurred version of what is on the specimen plane,

[TABLE]

In two dimensions, the blurring function is called the Airy pattern,

[TABLE]

where $J$ is the Bessel function of the first kind, $a$ is the radius of the aperture, $\lambda$ is the wavelength of the illumination, and $R$ is the focal length. The width of the central lobe of this function plays an important role in the achievable resolution—the minimum distance that two point sources must be separated for them to be distinguished. For visible light, this limit is around 250 nm.

For another perspective of the resolution limit, we note that the Fourier transform of the Airy pattern is a scaled version of the indicator function of a circle, i.e.,

[TABLE]

Functions for which the Fourier transform is zero outside of some radius, $w_{\text{max}}$, are called low-pass functions. Writing the measurement process in the Fourier domain we have

[TABLE]

and therefore the high-frequency content of $f$ has no effect on $g$ because everything outside of $w_{\text{max}}$ is multiplied by zero during the measurement. For a textbook treatment of the pattern disk, see, e.g., [11]; for excellent introductory information on microscopy, see [12].

Example 4** (Structured-illumination microscopy).**

The idea behind structured-illumination microscopy (SIM) is to use a special illumination pattern to modulate the image $f$, rearranging its frequency content so that its high frequencies are no longer destroyed by the measurement process. The continuous forward model is

[TABLE]

where $p$ is a complex exponential, $p({\boldsymbol{x}})=\textrm{e}^{\textrm{j}\left\langle{\boldsymbol{\omega}}_{0},{\boldsymbol{x}}\right\rangle}$. Viewed in the Fourier domain,

[TABLE]

To recover both the high-frequency and low-frequency content of $f$, the imaging process can be repeated with different illumination patterns, resulting in different values for the frequency shift, ${\boldsymbol{\omega}}_{0}$. For more information on SIM, see [13]

Definition 7**.**

A change of variables operator, $\Phi_{\varphi}:L_{2}(\mathcal{X})\rightarrow L_{2}(\mathcal{Y})$, is a linear operator such that $\Phi_{\varphi}f=f(\varphi(\cdot))$, where $\varphi:\mathcal{Y}\rightarrow\mathcal{X}$ is a diffeomorphism (a smooth function with a smooth inverse).

We have already seen one change of variables: the frequency shift caused by multiplication by a complex exponential. Using a change of variables operator, we can express a shift as

[TABLE]

See Figure 2.5 for another example.

Example 5** (X-ray computed tomography).**

The idealized forward model for a variety of CT imagining modalities, including X-ray CT, PET, SPECT, and electron tomography (ET) is called the X-ray transform. It measures the line integrals of an image and is usually given in 2D by

[TABLE]

where $\theta$ is an angle, ${\boldsymbol{\theta}}=(\cos\theta,\sin\theta)$ is a unit vector, and ${\boldsymbol{\theta}}^{\perp}=(-\sin\theta,\cos\theta)$ is its perpendicular. It turns out that the X-ray transform can also be written in terms of the building blocks we have already defined as

[TABLE]

where $\varphi(\omega)=\omega{\boldsymbol{\theta}}$ and the subscripts on the Fourier transforms remind us that that the image is two dimensional while its X-ray transform is one dimensional. We verify this by writing

[TABLE]

where the trick is that we can always use the orthogonal matrix ${\boldsymbol{\mathbf{A}}}=\begin{bmatrix}{\boldsymbol{\mathbf{\theta}}}&{\boldsymbol{\mathbf{\theta}}}^{\perp}\end{bmatrix}$ to express ${\boldsymbol{x}}$ as ${\boldsymbol{\mathbf{A}}}{\boldsymbol{\mathbf{a}}}$ for some ${\boldsymbol{\mathbf{a}}}$. The result is called the Fourier slice theorem or central slice theorem, because it relates the 1D Fourier transform of the X-ray transform of an image to a slice through its 2D Fourier transform. For more information see [14] for X-ray CT, with additional details on PET in [15] and ET in [16].

2.4 Discretization

We have, so far, presented a toolbox of operators that can be combined to model many biomedical imaging modalities. Generally speaking, we will need to implement these operators in software to solve image reconstruction problems. But, how can we do this when the definitions are in the continuous domain? That is, how can we represent $f\in L_{2}(\mathbb{R}^{n})$ in a computer? This is the problem of discretization.

One approach to the problem is to discretize at the end: formulate the reconstruction algorithm in the continuous domain and develop discrete versions of each of the necessary operations. Anywhere a function is needed, it is represented by a finite vector of its samples; the evaluation of the function at other points is achieved via interpolation, and integrals are replaced with sums by the trapezoidal rule, etc. We will discuss this approach further in Section 3.1.

We focus here on discretizing at the beginning via the synthesis operator.

Definition 8**.**

The synthesis operator $T_{\{\beta_{n}\}_{n=1}^{N}}:\mathbb{R}^{N}\rightarrow L_{2}(\mathbb{R}^{d})$, expresses a function as a sum of a finite number of basis functions,

[TABLE]

where $\beta_{n}\in L_{2}(\mathbb{R}^{d})$, $f$ denotes the original function and ${\boldsymbol{\mathbf{f}}}$ is the vector of its coefficients.

The synthesis operator allows for a simple discretization of any continuous forward model, $H:L_{2}(\mathbb{R}^{d})\rightarrow\mathbb{R}^{M}$, according to ${\boldsymbol{\mathbf{H}}}=H\circ T$, where ${\boldsymbol{\mathbf{H}}}$ is an $M\times N$ matrix, called the system matrix. The discretized imaging model is then ${\boldsymbol{\mathbf{g}}}={\boldsymbol{\mathbf{H}}}{\boldsymbol{\mathbf{f}}}$. This discretization approach is also called finite series expansion; see [17] for an early discussion. Note the dimensions of the system matrix: $M$ is the number of measurements and $N$ is the number unknowns to reconstruct; these can both easily be on the order of millions. So, while ${\boldsymbol{\mathbf{H}}}$ is a matrix, it is usually too large to store. Instead, its elements are computed whenever they are needed.

Recalling that we can view the linear measurement process as taking inner products with detector functions (2.4), we can express the elements of the system matrix as inner products between the basis functions and these detector functions,

[TABLE]

For simplicity, the basis functions are often chosen to be shift-invariant, $\beta_{n}({\boldsymbol{x}})=\beta({\boldsymbol{x}}-{\boldsymbol{x}}_{n})$, where the ${\boldsymbol{x}}_{n}$ form a regular grid, e.g., ${\boldsymbol{x}}_{n}\in\{{\boldsymbol{x}}=T{\boldsymbol{k}}+{\boldsymbol{x}}_{0}:{\boldsymbol{k}}\in[0,1,\dots,N_{1}-1]\times[0,1,\dots,N_{2}-1]\}$. They may also be multiplicatively separable,

[TABLE]

where (in a slight abuse of notation) $\beta({\boldsymbol{x}})$ is the $d$-dimensional separable basis function and $\beta(x)$ is the 1D function from which it is built.

Example 6** (Pixel model).**

Because the ultimate result of image reconstruction is a digital image, it is natural to use the square pixel as a basis function, i.e., $\beta({\boldsymbol{x}})=\operatorname{rect}({\boldsymbol{x}})$, see Figure 2.6 for an example. The downside of the pixel model is, generally, its lack of smoothness: It is not differentiable, so it is not appropriate for continuous forward models involving differentiation. It is also not bandlimited: its Fourier transform is

[TABLE]

and therefore may not interact well with forward models that amplify high frequencies. Despite this, the pixel model is ubiquitous in practice.

Although they are beyond the scope of the current discussion, there are many more tools that can be used for discretization, e.g., wavelets [18] and dictionary methods [19].

2.5 Summary

We have presented a set of building blocks—linear operators that can be composed to model a variety of physical measurement processes. In the following, we give a nonexhaustive list of these imaging modalities. It is important to note that these examples are simplified: state-of-the-art models typically account for additional physical effects. One must also include the effect of the measurement noise, which we discuss in the following sections.

2D or 3D tomography

uses X-rays. We can model the measurements as $\{(S_{X_{i}}\circ X_{\theta_{i}})\{f\}\}$, which is a set of sampled X-ray projections taken from different directions, where $i$ indexes the projection direction. Variations include parallel (where the sampling locations do not depend on the projection directions, $X_{i}=X$), cone beam, and spiral sampling patterns.

3D microscopy

uses fluorescence. We can model the measurements as $(S_{\mathbb{Z}^{3}}\circ C_{h})\{f\}$, which is a blurred, sampled version of $f$. Variations include brightfield, confocal, and light sheet microscopy.

SIM

uses fluorescence. We can model the measurements as $\{(S_{\mathbb{Z}^{2}}\circ C_{h}\circ M_{w_{i}})\{f\}\}_{i}$, which is a set of modulated, blurred, and sampled versions $f$, where $i$ indexes the modulation pattern. Variations include performing 3D reconstructions or using non-sinusoidal patterns.

PET

uses gamma rays. We can model the measurements as $\{(S_{X_{i}}\circ X_{\theta_{i}})\{f\}\}$, which, like in tomography, are sampled X-ray projections. We use the subscript $X_{i}$ to indicate that the sampling positions are dependent on the projection direction, and note that there may be as few as one sample per direction. Variations include time-of-flight PET.

MRI

uses magnetic fields and measures radio waves. We can model the measurements as $\{(S_{X}\circ\mathcal{F}\circ M_{w_{i}})\{f\}\}$, which are a set of samples of the Fourier transform of $f$ after weighting by several different coil sensitivities, indexed by $i$. There are several variations in the sampling pattern, including Cartesian and radial sampling, and the method can also be used for dynamic processes, e.g., cardiac MRI.

Optical diffraction tomography

uses coherent light. Although the problem is, in general, nonlinear, using the first Born approximation [20], we can model the measurements as $(S_{X}\circ\mathcal{F}\circ M_{h})\{f\}$, where $h$ is related to the incident wave. Variations include holography and grating-based interferometry.

In each case above, stated a measurement model that goes from a continuously defined $f:\mathbb{R}^{d}\to\mathbb{R}$ to a discrete ${\boldsymbol{\mathbf{g}}}\in\mathbb{C}^{M}$. By additionally making use of a synthesis operator, each of these can be converted to a system matrix that relates the discrete image we want to recover to the discrete measurements we have access to. This matrix plays a key role in reversing the measurement process to arrive at an estimate of the image of interest, which is the topic of the following sections.

2.5.1 Further Reading

Many of these building blocks are fundamental concepts in signal processing; for a textbook treatment, see the classic [21] or [22]. For more information on the modalities we have discussed, see the references listed at the end of each example. For a MATLAB software library based on this methodology, see [23].

Chapter 3 Classical Image Reconstruction

Assuming that we have an accurate, implementable forward model, $H$, and measurements ${\boldsymbol{\mathbf{g}}}$, our task now is to reconstruct an image. Our focus is on discretizing immediately; i.e., using a synthesis operator to create a system matrix, ${\boldsymbol{\mathbf{H}}}$, because doing so gives access to simple, implementable reconstruction algorithms. But, before dispensing with the operator $H$ in favor of the matrix ${\boldsymbol{\mathbf{H}}}$, we need to discuss direct inversion methods, which are derived in the continuum.

3.1 Direct Inversion

Many of the imaging models we have discussed are invertible under the assumption that there is no noise and a continuum of measurements: there is a formula that relates the measurements to the image we are trying to recover. As a concrete example, we recall that the Fourier transform is invertible. Thus, since MRI is based on the Fourier transform, it seems that the inverse Fourier transform should provide an MRI reconstruction technique. To develop a direct inversion algorithm, we formulate a forward model in the continuum (i.e., ignoring sampling), $H:L_{2}(\mathbb{R}^{d})\rightarrow L_{2}(\mathbb{R}^{d})$ and derive (analytically) its inverse, $H^{-1}$. Such an inverse will usually involve integrals of a continuous version of the measurements. To actually implement $H^{-1}$, we replace these integrals with sums over the known measurements and use interpolation whenever we need to evaluate a quantity in between its samples. See [24] for a early discussion of this approach.

Example 7** (Filtered back projection).**

The canonical example of direct inversion is the filtered back projection (FBP) algorithm for X-ray CT, which dates back to 1971 [25]. Starting with the X-ray transform as a forward model, we use the Fourier slice theorem to derive the following inverse

[TABLE]

with ${\boldsymbol{x}}\in\mathbb{R}^{2}$ and ${\boldsymbol{\theta}}^{\perp}=(-\sin\theta,\cos\theta)$. This equation tells us to take each projection (from 0 to $\pi$), convolve by a filter with impulse response $|\omega|$, and, finally, back project it (smear it across the image). The name filtered back projection is therefore a bit misleading because the filtering occurs before the back projection.

The formula (3.1) looks simple, but several approximations are needed to implement it from discrete measurements. First, the integral over $\theta$ must be replaced with a sum. Second, the filter must be implemented digitally, with careful handling of the boundary, high frequencies (because they are amplified), and zero frequency (because it is removed). Finally, the filtered projections must be interpolated to back project them. Despite these approximations, the FBP is a robust and high-quality method, especially when the number of projections is large.

3.2 Variational Methods

To avoid the challenges associated with discretizing a direct inversion method, we prefer to work with a linear forward model that takes the form of a matrix, ${\boldsymbol{\mathbf{H}}}$. In this setting, the simplest approach to image reconstruction would be to solve

[TABLE]

for ${\boldsymbol{\mathbf{f}}}$.111While the components of ${\boldsymbol{\mathbf{f}}}$ are actually expansion coefficients in some underlying basis (see Section 2.4), we shall, with a slight abuse of language, refer to them as pixels. The first problem with this formulation is that ${\boldsymbol{\mathbf{H}}}$ is rarely invertible—in most scenarios, we want to get as many image pixels, $N$, out of as few measurements, $M$, as possible, leading to a short, wide system matrix with no left inverse. Practically, this means that there will be many ${\boldsymbol{\mathbf{f}}}$’s that satisfy (3.2). If any of the measurement vectors are linearly dependent (which can happen in any setting, but is guaranteed when $M>N$), any amount of measurement noise will be enough to move ${\boldsymbol{\mathbf{g}}}$ out of the span of ${\boldsymbol{\mathbf{H}}}$, leading to (3.2) having no solution. Finally, even when a unique solution does exist ($M=N$), the solution can be very sensitive to noise, to the point were direct inversion of ${\boldsymbol{\mathbf{H}}}$ is useless (see Figure 3.1 for an example). Thus, (3.2) is an example of an ill-posed problem.

Definition 9**.**

A well-posed problem satisfies:

A solution exists. 2. 2.

The solution is unique. 3. 3.

The solution is a continuous function of the measurements.

A problem that is not well-posed is call ill-posed.

A more promising approach, which is the focus of this section and the next, is to formulate reconstruction as a minimization problem. These approaches are sometimes called variational methods; while this term has precise meanings in fields such as differential equations, quantum mechanics, and Bayesian inference, in our context it only means we aim to minimize something. As a first example, we can seek a reconstruction that minimizes the disagreement with the measurements according to

[TABLE]

From a geometric perspective, we see that this problem must have at least one solution: we can project ${\boldsymbol{\mathbf{g}}}$ onto the span of ${\boldsymbol{\mathbf{H}}}$, resulting in a denoised measurement, ${\boldsymbol{\mathbf{g}}}^{\prime}$ and solve ${\boldsymbol{\mathbf{H}}}{{\boldsymbol{\mathbf{f}}}}={\boldsymbol{\mathbf{g}}}^{\prime}$. However, the problem may still be ill-posed. Equating the gradient (i.e., the vector of partial derivatives with respect to the unknown ${\boldsymbol{\mathbf{f}}}=(f_{1},\dots,f_{N})$) of the objective function in (3.3) with zero results in the so-called normal equation,

[TABLE]

which ${\boldsymbol{\mathbf{f}}}$ must satisfy to be a solution. Here, ${\boldsymbol{\mathbf{H}}}^{*}$ denotes the conjugate transpose. When the number of linearly-independent measurements (rank of ${\boldsymbol{\mathbf{H}}}$) is at least the number of unknowns (length of ${\boldsymbol{\mathbf{f}}}$), then there is a unique solution $({\boldsymbol{\mathbf{H}}}^{*}{\boldsymbol{\mathbf{H}}})^{-1}{\boldsymbol{\mathbf{H}}}^{*}{\boldsymbol{\mathbf{g}}}={\boldsymbol{\mathbf{f}}}$. Otherwise, there are infinitely-many solutions. Specifically, the solutions form an affine subspace of $\mathbb{R}^{N}$, $\mathcal{N}({\boldsymbol{\mathbf{H}}})+\tilde{{\boldsymbol{\mathbf{f}}}}$, where $\mathcal{N}({\boldsymbol{\mathbf{H}}})=\{{\boldsymbol{\mathbf{f}}}:{\boldsymbol{\mathbf{H}}}{\boldsymbol{\mathbf{f}}}={\boldsymbol{\mathbf{0}}}\}$ is called the null space of ${\boldsymbol{\mathbf{H}}}$ and $\tilde{{\boldsymbol{\mathbf{f}}}}$ is one solution of (3.3).

Remark 1** (Does the square matter?).**

One challenge in understanding the biomedical reconstruction literature is the plethora of similar-looking problem statements that may or may not be equivalent. The first such example is

[TABLE]

These expressions are equivalent because the function $|\cdot|^{2}:x\mapsto|x|^{2}$ is monotonically increasing over the nonnegative real numbers. We will see a less-trivial generalization of this idea in Remark 2.

The version without the square is more natural: it gives a reconstruction that agrees with the measurements in the sense of Euclidean distance. The benefit of the square is that it cancels the square root in the norm, which leads to a simpler expression for the gradient.

Unfortunately, the infinitely-many solutions case is more the norm than the exception. This is, again, because we want the highest possible resolution, $N$, from a fixed number of measurements, $M$, determined by the hardware. As a result, the number of unknowns is often larger than the number of measurements. Thus, (3.3) is generally an ill-posed problem.

Example 8** (Null space).**

Consider a toy problem where

[TABLE]

The true measurement process is ${\boldsymbol{\mathbf{g}}}={\boldsymbol{\mathbf{H}}}{\boldsymbol{\mathbf{f}}}_{\text{true}}+{\boldsymbol{\mathbf{n}}}=\begin{bmatrix}3&-1&2.1\end{bmatrix}^{T}$ and we want to reconstruct by solving (3.3). If we use a standard solver (MATLAB’s pcg function) on the normal equation (3.4), the result is $\tilde{{\boldsymbol{\mathbf{f}}}}\approx\begin{bmatrix}1.70&0.37&1.33\end{bmatrix}^{T}$, giving a sum of squared errors of 0.0033. But, by changing the initialization of the solver, infinitely many equally good solutions can be obtained:

[TABLE]

This effect occurs because ${\boldsymbol{\mathbf{H}}}$ has a non-zero null space, specifically ${\boldsymbol{\mathbf{H}}}\begin{bmatrix}1&-1&-1\end{bmatrix}^{T}={\boldsymbol{\mathbf{0}}}$; note that each solution can be obtained from another by adding a scaled version of this null space vector.

3.3 Tikhonov regularization

When a formulation is ill-posed in the sense of having no solution, it is a dead end. But, when a formulation merely has infinitely-many solutions, all we need is a mechanism to pick one of them to obtain a well-posed problem; Tikhonov regularization provides a classical mechanism to do this.

The idea is to introduce a regularization functional, $\|{\boldsymbol{\mathbf{L}}}{\boldsymbol{\mathbf{f}}}\|^{2}_{2}$, where ${\boldsymbol{\mathbf{L}}}$ is a linear operator that measures the “roughness” of ${\boldsymbol{\mathbf{f}}}$ (or some other undesirable property). The reconstruction can then be formulated as

[TABLE]

where $\sigma$ is a scalar constant selected based on the expected level of noise. An equivalent formulation is

[TABLE]

where $\lambda$ is a scalar constant controlling the strength of the regularization (see Figure 3.2). See Remark 3 for a discussion of the equivalence of these two problems.

Remark 2** (Does the square matter now?).**

We saw in Remark 1 that the square in the least squares formulation is only a mathematical convenience. Now, with the addition of the regularization term, we can compare a more general case,

[TABLE]

*where $J_{1}$ and $J_{2}$ are positive, differentiable functionals, $\lambda_{\text{L}}$ and $\lambda_{\text{R}}$ are positive numbers, and $\Psi_{1}$ and $\Psi_{2}$ are monotonically increasing, differentiable functions. Intuitively, these are different problems when $\lambda_{\text{L}}=\lambda_{\text{R}}$ because the *$\Psi$s will change the relative importance of $J_{1}$ and $J_{2}$ (which we can think of as representing the data and regularization terms, respectively). The problems are equivalent in the sense that for any value of $\lambda_{\text{L}}$, we can always find a value of $\lambda_{\text{R}}$ such that the problems are the same. To see this, take a solution of the left problem, $\tilde{{\boldsymbol{\mathbf{f}}}}$. At this point, the gradients of $J_{1}$ and $\lambda_{\text{L}}J_{2}$ are equal and opposite. In the right problem, at $\tilde{{\boldsymbol{\mathbf{f}}}}$, the chain rule tells us that the gradient of the $J_{1}$ term has been multiplied by $\Psi_{1}^{\prime}(J_{1}(\tilde{{\boldsymbol{\mathbf{f}}}}))$ and the gradient of the $J_{2}$ term has been multiplied by $\Psi_{2}^{\prime}(J_{2}(\tilde{{\boldsymbol{\mathbf{f}}}}))$. Because both of these multipliers are positive, we can set

[TABLE]

which makes $\tilde{{\boldsymbol{\mathbf{f}}}}$ a critical point for the right problem in (3.9). We can make the same argument in the reverse direction, and a similar argument can be made with nondifferentiable functions by using subgradients.

In practice, we set $\lambda$ by hand to give the best-looking reconstruction; this argument tells us that if we can find a good $\lambda$ for the left formulation, we will be able to find a good $\lambda$ in the right formulation. As a result, we have a lot of flexibility to add or remove, e.g., squares, from data fidelity and regularization terms without affecting the set of solutions we get.

In the Tikhonov formulation, the normal equation becomes

[TABLE]

This expression reveals the main advantage of the formulation: the problem will have a unique solution when the intersection of the null spaces of ${\boldsymbol{\mathbf{H}}}$ and ${\boldsymbol{\mathbf{L}}}$ is zero. One easy way to achieve this is to let ${\boldsymbol{\mathbf{L}}}$ be full rank, e.g., when ${\boldsymbol{\mathbf{L}}}$ is the identity matrix, the solution is unique and the regularization penalizes large values in ${\boldsymbol{\mathbf{f}}}$.

Remark 3** (Constraints versus regularization).**

We have now seen (and will continue to see) data and regularization terms on ${\boldsymbol{\mathbf{f}}}$ expressed as both functionals (also called penalties), e.g., $\|{\boldsymbol{\mathbf{L}}}{\boldsymbol{\mathbf{f}}}\|_{2}$, and as constraints, e.g., $\|{\boldsymbol{\mathbf{g}}}-{\boldsymbol{\mathbf{H}}}{\boldsymbol{\mathbf{f}}}\|_{2}\leq K$. The natural question is to compare

[TABLE]

where $J_{1}$ and $J_{2}$ are positive functionals and $\lambda$ and $\sigma$ are positive numbers. We will call these the penalized and constrained formulations; the first is also called Tikhonov regularization and the second Ivanov or Morozov regularization, depending on whether $J_{1}$ is the data term (Ivanov) or the regularization term (Morozov). We can afford to be loose in interchanging these forms because, roughly, solutions of one are also solutions of the other (when the $\lambda$ and $\sigma$ are chosen correctly). We make this notion precise in the next two theorems; for a discussion of the engineering implications, skip ahead.

Theorem 1**.**

For every choice of $\lambda$, for each solution of the penalized problem, there is a choice of $\sigma$ that makes it a solution of the constrained problem.

Proof.

Fix $\lambda$ and assume ${\boldsymbol{\mathbf{f}}}^{P}$ is a solution to the penalized problem. Set $\sigma=J_{2}({\boldsymbol{\mathbf{f}}}^{P})$. Now for sake of contradiction, assume there exists ${\boldsymbol{\mathbf{f}}}^{C}$ that is a better solution of the constrained problem than ${\boldsymbol{\mathbf{f}}}^{P}$, i.e., $J_{1}({\boldsymbol{\mathbf{f}}}^{C})<J_{1}({\boldsymbol{\mathbf{f}}}^{P})$ and $J_{2}({\boldsymbol{\mathbf{f}}}^{C})\leq\sigma$. Then we have $J_{1}({\boldsymbol{\mathbf{f}}}^{R})+\lambda J_{2}({\boldsymbol{\mathbf{f}}}^{R})<J_{1}({\boldsymbol{\mathbf{f}}}^{L})+\lambda J_{2}({\boldsymbol{\mathbf{f}}}^{L})$, which contradicts ${\boldsymbol{\mathbf{f}}}^{P}$ being a solution to the penalized problem. Thus, there is no such ${\boldsymbol{\mathbf{f}}}^{C}$, which means ${\boldsymbol{\mathbf{f}}}^{P}$ is a solution to the constrained problem. ∎

Theorem 2**.**

Assume that the penalized problem has a solution (not necessarily unique) for every $\lambda\geq 0$. Then, for every choice of $\sigma$ and any each solution of the constrained problem, there is a choice of $\lambda$ that makes it a solution of the penalized problem.

Proof.

Fix $\sigma$ and assume ${\boldsymbol{\mathbf{f}}}^{C}$ is a solution to the constrained problem. Consider the function $F(\lambda)=\min_{{\boldsymbol{\mathbf{f}}}}J_{1}({\boldsymbol{\mathbf{f}}})+\lambda J_{2}({\boldsymbol{\mathbf{f}}})$. We are searching for a value, $\lambda^{*}$, such that $F(\lambda^{*})=J_{1}({\boldsymbol{\mathbf{f}}}^{C})+\lambda^{*}J_{2}({\boldsymbol{\mathbf{f}}}^{C})$, because this will make ${\boldsymbol{\mathbf{f}}}^{C}$ a solution to the penalized problem. We know that, for any $\lambda^{*}$, $F(0)=\min_{{\boldsymbol{\mathbf{f}}}}J_{1}({\boldsymbol{\mathbf{f}}})\leq J_{1}({\boldsymbol{\mathbf{f}}}^{C})+\lambda^{*}J_{2}({\boldsymbol{\mathbf{f}}}^{C})$ and that we can pick $\lambda^{\wedge}$ large enough such that $F(\lambda^{\wedge})\geq J_{1}({\boldsymbol{\mathbf{f}}}^{C})+\lambda^{*}J_{2}({\boldsymbol{\mathbf{f}}}^{C})$. We can show that $F$ is continuous; thus, the intermediate value theorem gives us that there exists $\lambda^{*}$ such that $F(\lambda^{*})=J_{1}({\boldsymbol{\mathbf{f}}}^{C})+\lambda^{*}J_{2}({\boldsymbol{\mathbf{f}}}^{C})$, which makes ${\boldsymbol{\mathbf{f}}}^{C}$ a solution of the penalized problem. ∎

Practically, these problems are equivalent; theoretically, the relationship is more subtle (e.g., we have not shown that the solution sets can be made the same, only that the problems can be made to share a solution). When more is assumed about the functionals, we can say more. Generally, which form is chosen is a matter of personal preference. For a much deeper look at this issue, see **[26]**.

3.4 Bayesian Formulation

We now consider a statistical measurement model

[TABLE]

where ${\boldsymbol{\mathbf{H}}}$ remains a deterministic measurement matrix, but where ${\boldsymbol{\mathbf{f}}}$ and ${\boldsymbol{\mathbf{n}}}$ (and therefore ${\boldsymbol{\mathbf{g}}}$) are random variables. If we know (or can model or estimate) the distribution of ${\boldsymbol{\mathbf{n}}}$, then the model (3.13) lets us determine the probability distribution function (PDF) of the measurements given a fixed image, $p({\boldsymbol{\mathbf{g}}}\mid{\boldsymbol{\mathbf{f}}})$.

To do reconstruction, we are more interested in the distribution of ${\boldsymbol{\mathbf{f}}}$ conditional on the measurements, which we can obtain using Bayes’ rule as

[TABLE]

So, for a fixed set of measurements ${\boldsymbol{\mathbf{g}}}$, (3.14) is a function that states how likely a given ${\boldsymbol{\mathbf{f}}}$ is to have given rise to ${\boldsymbol{\mathbf{g}}}$. It requires a prior—i.e., a model for $p({\boldsymbol{\mathbf{f}}})$—which is analogous to the choice of a regularizer in the variational formulations we already presented. Unfortunately, it also requires computation of $p({\boldsymbol{\mathbf{g}}})$, which is typically intractable (because it requires integrating over all ${\boldsymbol{\mathbf{f}}}$’s that can create ${\boldsymbol{\mathbf{g}}}$). So, while (3.14) in principle gives all the information we might want in a reconstruction, we usually can not use it directly.

Using the same statistical model, we can build algorithms that seek a single image satisfying some measure of optimality. In one such approach, we seek the minimum mean square error (MMSE) solution, which is given by

[TABLE]

where here $\tilde{{\boldsymbol{\mathbf{f}}}}({\boldsymbol{\mathbf{g}}})$ is the reconstruction (which depends on the measurements, ${\boldsymbol{\mathbf{g}}}$) and ${\boldsymbol{\mathbf{f}}}$ is the ground truth. One way to think of the MMSE solution is that if it is used to perform many reconstructions from many different measurements, it will have the lowest average error among all algorithms (assuming that all the model assumptions are correct). By equating the gradient of the functional in (3.15) with zero, we find that the MMSE is given by $\tilde{{\boldsymbol{\mathbf{f}}}}=\mathbb{E}({\boldsymbol{\mathbf{f}}}\mid{\boldsymbol{\mathbf{g}}})$, which is the expectation of the conditional PDF given by (3.14). Unfortunately, it is usually difficult to build algorithms to find the MMSE solution, unless both ${\boldsymbol{\mathbf{f}}}$ and ${\boldsymbol{\mathbf{n}}}$ are multivariate Gaussians.

A different approach is to find the maximum a posteriori (MAP) solution, which is the mode of (3.14),

[TABLE]

Put another way, this is the ${\boldsymbol{\mathbf{f}}}$ that is most likely to have generated the measured ${\boldsymbol{\mathbf{g}}}$. As a word of caution: when the prior on ${\boldsymbol{\mathbf{f}}}$ and noise model are correct, MAP solutions do not, in general, follow them; see [27] for more discussion. However, the MAP approach is useful because it provides a statistically coherent way to account for noise and priors on the image, while also leading to optimization problems that can be efficiently solved.

As a specific example of MMSE and MAP solutions, we return to the measurement model (3.13) and consider the case where the elements of the noise ${\boldsymbol{\mathbf{n}}}$ are independent and identically distributed (i.i.d.) Gaussian with zero mean and variance $\sigma^{2}$, and ${\boldsymbol{\mathbf{f}}}$ is a zero-mean Gaussian process with invertible covariance matrix $\mathbb{E}({\boldsymbol{\mathbf{f}}}{\boldsymbol{\mathbf{f}}}^{*})={\boldsymbol{\mathbf{C}}}$. In this case, we can transform the MAP objective (3.16) using the negative log and arrive at

[TABLE]

which has a closed-form solution

[TABLE]

One can also derive the corresponding MMSE solution by finding the expectation of (3.14), which essentially requires some lengthy algebra involving the underlying multivariate Gaussian distributions. The result is

[TABLE]

which is also called the Wiener filter. It turns out that (3.18) and (3.19) are equal. (One way to verify this is to check that both satisfy the normal equation of (3.17), ${\boldsymbol{\mathbf{H}}}^{*}{\boldsymbol{\mathbf{g}}}=({\boldsymbol{\mathbf{H}}}^{*}{\boldsymbol{\mathbf{H}}}+\sigma^{2}{\boldsymbol{\mathbf{C}}}^{-1}){\boldsymbol{\mathbf{f}}}$, which is meaningful because the solution to (3.17) is unique.) And, further, note that (3.17) is equivalent to the Tikhonov formulation (3.8) with ${\boldsymbol{\mathbf{L}}}={\boldsymbol{\mathbf{C}}}^{-\frac{1}{2}}$ and $\lambda=\sigma^{2}$; thus, the optimal choice of ${\boldsymbol{\mathbf{L}}}$ is the whitening operator matched with ${\boldsymbol{\mathbf{f}}}$. So, for the special case of Gaussian denoising of a Gaussian process, the MAP and MMSE solutions are the same, and they can also be obtained from a specific variational formulation. But, for many other image reconstruction problems, this equivalence does not hold. For a more in-depth discussion of the relationship of MAP and MMSE, see [28].

3.5 Iterative Reconstruction

We have now seen several reconstruction problems that have a closed-form solution. Unfortunately, these usually involve inverting the normal matrix, ${\boldsymbol{\mathbf{H}}}^{*}{\boldsymbol{\mathbf{H}}}$, which is large (its size is number of pixels squared). So, for practical implementations, we turn to iterative algorithms. Iterative algorithms for solving convex problems of this form are, by now, a textbook topic [29]; we present one such algorithm here to give a sense of their form.

Given a Tikhonov reconstruction problem,

[TABLE]

The gradient descent algorithm involves iterates of the form

[TABLE]

where the gradient, $\nabla J({\boldsymbol{\mathbf{f}}})$, is given by

[TABLE]

Setting $\gamma$ to a small constant in (3.21) gives a workable algorithm that is also rather flexible. For instance, convex constraints, such as the positivity of ${\boldsymbol{\mathbf{f}}}$, can be enforced by projecting onto the feasible region after each iteration, an algorithm known as projected gradient descent. Gradient descent is good starting point for building iterative algorithms, but there are many alternative. The choice of the best algorithm will generally depend on the specifics of problem at hand. For example, the conjugate gradient (CG) algorithm solves the same problem more efficiently, at the cost of increased memory usage [30].

Remark 4** (Convolutional normal matrix).**

It turns out that the normal matrix, ${\boldsymbol{\mathbf{H}}}^{*}{\boldsymbol{\mathbf{H}}}$, appears in many iterative reconstruction schemes and is often the computational bottleneck. The case of gradient descent (3.22) is typical: reconstruction requires a single back projection, ${\boldsymbol{\mathbf{H}}}^{*}{\boldsymbol{\mathbf{g}}}$, followed by one normal operation $({\boldsymbol{\mathbf{H}}}^{*}{\boldsymbol{\mathbf{H}}}+\lambda{\boldsymbol{\mathbf{L}}}^{*}{\boldsymbol{\mathbf{L}}}){\boldsymbol{\mathbf{f}}}$ per iteration. Often, we select ${\boldsymbol{\mathbf{L}}}$ to be a convolution (e.g., a high-pass filter), which means that ${\boldsymbol{\mathbf{L}}}^{*}{\boldsymbol{\mathbf{L}}}$ is a convolution as well. If ${\boldsymbol{\mathbf{H}}}^{*}{\boldsymbol{\mathbf{H}}}$ is a convolution, then the entire operation is a convolution. This is important because a discrete, linear convolution can be computed using the FFT, resulting in a fast algorithm for the normal operation. Due to boundary conditions, the matrices may be Toeplitz (diagonal-constant) rather than circulant, but FFT-based convolution can still be used for fast algorithms provided that proper zero-padding is applied; see **[21]** Section 8.7.

For example, with a specific choice of discretization, it is possible to model parallel-ray X-ray CT with a convolutional normal matrix **[31]**. This approach has been used in cryo-EM **[32]** and synchrotron microtomography **[33]**.

3.6 Summary

Classical image reconstruction methods fall into either direct methods or variational (regularized, iterative) methods. Direct methods are derived in the continuous domain; they are fast to apply and give good results when the number of measurements is high. Variational methods involve minimizing an objective function that usually consists of a data term and one or more regularization terms. These methods tend to be more robust to noise than direct methods, at the cost of increased computation.

Chapter 4 Sparsity-Based Image Reconstruction

Over the past two decades, numerous new reconstruction methods have emerged, many of which continue to use the variational framework, i.e., reconstruction is performed by minimizing a cost functional. Of these, one dominant paradigm has been that of sparsity—the idea that a high dimensional image, ${\boldsymbol{\mathbf{f}}}$, can be represented by only a small number of nonzero coefficients. The concept of sparsity is included in reconstruction problems by replacing the quadratic regularization functional, $\|{\boldsymbol{\mathbf{L}}}{\boldsymbol{\mathbf{f}}}\|^{2}_{2}$, of Tikhonov regularization to one that promotes sparsity. Especially when the number or quality of the measurements is low, such a change can drastically improve the reconstruction quality. In this section, we will introduce these sparsity-based image reconstruction techniques.

4.1 Sparsity and Compressive Sensing

As we saw in the previous section, using classical reconstruction techniques, recovering an image of $N$ pixels requires on the order of $N$ measurements. The motivating idea behind compressive sensing is that many images can be compressed by applying a suitable linear transform and throwing away small values. One way to model this compressibility is to say that, when expressed in the correct basis, many images have few nonzero coefficients. Mathematically, one can write that

[TABLE]

It turns out that, if ${\boldsymbol{\mathbf{f}}}$ fits this model, it can be recovered from a number of measurements much smaller than its number of pixels. Popular choices for the sparsifying transform include finite differences, the Fourier transform, the discrete cosine transform, and the wavelet transform; we compare a few of these in Figure 4.1.

We begin our discussion by considering ${\boldsymbol{\mathbf{L}}}={\boldsymbol{\mathbf{I}}}$, so that $\|{\boldsymbol{\mathbf{f}}}\|_{0}=K$. Given measurements, ${\boldsymbol{\mathbf{g}}}={\boldsymbol{\mathbf{H}}}{\boldsymbol{\mathbf{f}}}$, the compressive sensing problem is

[TABLE]

where the $\ell_{0}$ norm, $\|{\boldsymbol{\mathbf{f}}}\|_{0}$, counts the number of nonzero entries of ${\boldsymbol{\mathbf{f}}}$. In general, (4.2) is very challenging to solve, however under certain conditions on ${\boldsymbol{\mathbf{H}}}$ [34, 35], the solution is unique and the problem is equivalent to

[TABLE]

where $\|{\boldsymbol{\mathbf{f}}}\|_{1}=\sum_{n}|[{\boldsymbol{\mathbf{f}}}]_{n}|$, for a particular choice of $\lambda$. (Note the change from the $\ell_{0}$ to $\ell_{1}$ norm; this is what makes the problem tractable.)

A more general formulation of the sparse recovery problem includes a sparsifying transform, ${\boldsymbol{\mathbf{L}}}$, and is formulated as

[TABLE]

which, under certain conditions [36], we can equivalently reformulate into an $\ell_{1}$ problem,

[TABLE]

While the formulation (4.5) is quite similar to the Tikhonov formulation (3.8), changing to the $\ell_{1}$ norm has important consequences; we contrast these two formulations in Figures 4.2 and 4.3.

The formulation (4.5) is called the analysis form of the regularization, because the matrix ${\boldsymbol{\mathbf{L}}}$ retrieves the sparse coefficients from the original signal. This is in contrast to the synthesis form,

[TABLE]

where the matrix $\tilde{{\boldsymbol{\mathbf{L}}}}$ now acts to construct the signal from its sparse coefficients. If the synthesis transform has a left inverse, this is merely a change in notation; otherwise, the two forms are meaningfully different. For more discussion of these two forms, see [37].

4.2 Representer Theorems for $\ell_{2}$ and $\ell_{1}$ Problems

Another perspective on sparsity-promoting regularization is given by representer theorems that specify the form of solutions to certain minimization problems. For example, for Tikhonov regularization with ${\boldsymbol{\mathbf{L}}}$ being the identity matrix, we can state the following representer theorem, which is a simplified special case of a result from [38].

Theorem 3** (Convex Problem with $\ell_{2}$ Minimization).**

The problem

[TABLE]

has a unique solution of the form

[TABLE]

for a suitable set of coefficients, ${\boldsymbol{\mathbf{a}}}\in\mathbb{R}^{M}$.

The useful insight here is that the solution to an $\ell_{2}$-regularized problem always has the form of a weighted sum of the original measurement vectors, i.e., the columns of ${\boldsymbol{\mathbf{H}}}$. Moreover, the number of elements in the sum is equal to the number of measurements, $M$, so there is no reason to expect ${\boldsymbol{\mathbf{f}}}$ to be sparse, unless these measurement vectors are themselves sparse in some transform domain. If we choose instead to minimize a function including an invertible regularization operator, $\|{\boldsymbol{\mathbf{L}}}{\boldsymbol{\mathbf{f}}}\|^{2}_{2}$, we can state a similar theorem with ${\boldsymbol{\mathbf{H}}}^{*}$ replaced by $({\boldsymbol{\mathbf{L}}}^{*}{\boldsymbol{\mathbf{L}}})^{-1}{\boldsymbol{\mathbf{H}}}^{*}$.

Contrast Theorem 3 with the $\ell_{1}$ representer theorem [38],

Theorem 4** (Convex Problem with $\ell_{1}$ Minimization).**

The set

[TABLE]

is convex, compact, and has extreme points of the form

[TABLE]

where $\{{\boldsymbol{\mathbf{e}}}_{n}\}_{n=1}^{N}$ is the standard basis for $\mathbb{R}^{N}$ (i.e., unit vectors pointing along each of the axes) and for a suitable set of coefficients ${\boldsymbol{\mathbf{a}}}\in\mathbb{R}^{K}$ and locations ${\boldsymbol{\mathbf{n}}}\in\{1,2,\dots,N\}^{K}$ with $K\leq M$.

Looking at the form of $\tilde{{\boldsymbol{\mathbf{f}}}}$, we see that it is sparse: it has fewer nonzero terms than the number of measurements. The amplitudes of its nonzero terms are given by ${\boldsymbol{\mathbf{a}}}$ and their locations are given by ${\boldsymbol{\mathbf{n}}}$. One complication is that the solution to the $\ell_{1}$ problem is not, in general, unique (though with additional conditions it is). Thus, the theorem is stated in terms of the extreme points of the solution set, i.e., those solutions that cannot be expressed as linear combinations of other solutions. When ${\boldsymbol{\mathbf{L}}}$ is not the identity, a similar theorem can be stated where the extreme points of the solution set are built out of a sparse linear combination of dictionary vectors that depend on ${\boldsymbol{\mathbf{L}}}$ [38].

4.3 Bayesian View

We can turn to a statistical perspective on sparsity-promoting regularization, basing the regularization term on a statistical model of the signal. When we expect the signal model to be accurate, this approach allows us to design reconstruction algorithms that are optimal with regard to chosen statistical criteria. Even when the signal model is only an educated guess, the statistical formulation can be helpful in designing new regularization terms with good empirical performance. Here, we follow the approach of [39], which places both $\ell_{1}$ and $\ell_{2}$ regularization along a spectrum of sparse regularizers.

Working for the moment in 1D, we model the signal as a stochastic process specified by the innovation model

[TABLE]

where the discrete innovation, ${\boldsymbol{\mathbf{u}}}$, has i.i.d. elements with an infinitely divisible PDF, $p_{U}$, and ${\boldsymbol{\mathbf{L}}}$ is a whitening operator. The restriction to this family is motivated by the property that this is the only configuration where (4.7) has a continuous-domain counterpart as a stochastic differential equation, $w=L\{f\}$, that specifies the family of sparse stochastic processes [39]. Yet, the family of infinitely divisible PDFs is still rather larger: it includes the Gaussian distribution, as well as several sparser (heavy-tailed or with more mass at the origin) distributions, e.g., the Laplace, compound-Poisson, and Cauchy.

Combining this signal model with a forward operator ${\boldsymbol{\mathbf{H}}}$ and a Gaussian noise model, we can express the posterior

[TABLE]

Taking the negative log and minimizing (i.e., performing MAP estimation), we have

[TABLE]

where $\Phi_{U}(u)=-\log p_{U}(u)$ is called the potential function.

We give a few examples of infinite divisible distributions and the corresponding potential functions in Table 4.1. In the first line, note that the Gaussian potential corresponds to the squared Euclidean norm (with scaling and plus a constant); thus, $\ell_{2}$ regularization corresponds to MAP estimation of a stochastic process with a Gaussian innovation (${\boldsymbol{\mathbf{u}}}$ in (4.7)). This is a generalization of the example from Section 3.4 because, here, our model if ${\boldsymbol{\mathbf{f}}}$ involves the operator ${\boldsymbol{\mathbf{L}}}$. In the second line, we see that the Laplace potential corresponds to the $\ell_{1}$ norm, which means that $\ell_{1}$ regularization corresponds to MAP estimation of a stochastic process with a Laplace-distributed innovation. This links with the concept of sparsity because the Laplace distribution is sparser than the Gaussian. The Student distribution (third line of Table 4.1) is the sparsest of the three, and, again, its potential function is closely related to sparsity-promoting regularization. Specifically, it is an upper bound on the $\log$ prior,

[TABLE]

which can be used as a relaxation of the $\ell_{0}$ norm in certain compressive sensing problems without changing the minima [40].

4.4 Algorithms

A wide range of algorithms have been developed for solving problem like (4.5). One key concept for understanding them is that of the proximal operator.

Definition 10**.**

The proximal operator, $\operatorname{prox}_{\Phi}:\mathbb{R}^{N}\to\mathbb{R}^{N}$, of a convex function, $\Phi:\mathbb{R}^{N}\to\mathbb{R}$, is defined by

[TABLE]

Intuitively, it finds an ${\boldsymbol{\mathbf{f}}}$ with a low $\Phi({\boldsymbol{\mathbf{f}}})$, while being close (proximal) to ${\boldsymbol{\mathbf{u}}}$.

If $\Phi$ acts pixel by pixel, then the optimization problem that the $\operatorname{prox}$ solves can be decomposed into $N$ scalar problems, one for each pixel. And, in many cases, these scalar problems can be solved efficiently and in parallel.

The $\operatorname{prox}$ is used in a class of algorithms called forward-backward splitting [41] or iterative shrinkage and thresholding (ISTA) [42] that solve problems of the form

[TABLE]

with $J_{1}$ differentiable and $J_{2}$ having a fast $\operatorname{prox}$. For example, these approaches can solve the sparse reconstruction problem (4.5) when ${\boldsymbol{\mathbf{L}}}$ is the identity. They achieve this by alternating a step down the gradient of $J_{1}$ with an application of the $\operatorname{prox}$ of $J_{2}$, which can be shown to converge. In fact, this method is closely related to projected gradient descent (discussed in Section 3.5) because when a constraint is expressed as a regularization functional (which takes the value infinity whenever the constraint is violated) the associated $\operatorname{prox}$ is just projection onto the constraint region.

When the operator ${\boldsymbol{\mathbf{L}}}$ is involved, the sparse reconstruction problem is typically much more difficult to solve. A useful paradigm is to split the data term and regularization terms by introducing an auxiliary variable. Using the sparse MAP formulation (4.11) as an example, splitting results in

[TABLE]

This problem is equivalent to the original one, but written in the form required by the alternating direction method of multipliers (ADMM) [43]. The ADMM then specifies an iterative algorithm,

[TABLE]

It turns out that these three problems are significantly simpler to solve than the original. The update (4.16) is a quadratic problem that can be solved using the techniques from classical image reconstruction (Section 3.5), e.g., the CG algorithm. For the update (4.17) we can use the $\operatorname{prox}$ of $\Phi_{U}$, which is known in closed-form for many potential functions. Finally, the update (4.18) only requires application of ${\boldsymbol{\mathbf{L}}}$.

4.5 Summary

Modern reconstruction techniques focus on sparsity-promoting regularization, with the goal of reconstructing images from very few measurements. This is possible only when the image to be recovered can be sparsely represented in some transform domain. See Figure 4.5 for a comparison between classical variational reconstructions and modern sparsity-based ones. The algorithms to solve sparse reconstruction problems often involve splitting the objective function; each iteration of these algorithms requires solving a classical reconstruction problem, which makes them computationally heavy.

4.5.1 Further Reading

Our discussion of regularization is a simplified version of a much more complex story; see [45] for a textbook treatment and [46] for more historical details and mathematical depth. See [47] for an introduction to compressive sensing, and [48] for an overview of the fields where it is applied. For examples of the statistical view presented here applied to real problems, see [39].

Chapter 5 The Learning (R)Evolution

The field of machine learning, which deals with creating computer programs that can improve with training [49], has existed for decades and has even been used for solving image reconstruction problems as early as in the 1980s [50]. Beginning with the emergence of one particular machine learning algorithm, the convolutional neural network (CNN), as a powerful and practical tool for object recognition in 2012 [51], there has been a surge in interest in applying learning-based methods to a wide variety of problems in image processing and computer vision, including image reconstruction. Our goal in this section is to give a broad overview of these approaches.

So far, we have considered a setting where we have measurements, ${\boldsymbol{\mathbf{g}}}$, and an understanding of the measurement device from which we build a forward model, ${\boldsymbol{\mathbf{H}}}$. Additionally, we may have some knowledge of the image we want to reconstruct, allowing us to design a variational regularization term. In the learning formulation, we assume that we also have access to training data: $T$ pairs of measurements and their corresponding ground truth reconstructions, $\{{\boldsymbol{\mathbf{g}}}_{t},{\boldsymbol{\mathbf{f}}}_{t}\}_{t=1}^{T}$. The inclusion of training data changes the task of the engineer from designing a reconstruction procedure that maps measurements to images, $\mathcal{R}:{\boldsymbol{\mathbf{g}}}\mapsto{\boldsymbol{\mathbf{f}}}$, to designing a learning procedure that maps training sets to reconstruction procedures, $\mathcal{L}:\{{\boldsymbol{\mathbf{g}}}_{t},{\boldsymbol{\mathbf{f}}}_{t}\}_{t=1}^{T}\mapsto\mathcal{R}$.

We have organized this section along a spectrum of how much is learned. On one end of the spectrum are methods where the training data is used to improve a part of a variational method, usually the forward model or the regularization term. On the other end are pipeline or end-to-end methods, where most or all of the reconstruction procedure is a generic regression function with its parameters learned during training. We also discuss a few methods that do not fit into our learning formulation as well as the problem of how to acquire training data.

5.1 Learning the Forward Model

Instead of relying on a physical model of the imaging setup alone to design the forward model, ${\boldsymbol{\mathbf{H}}}$, we can use training data to estimate the forward model. This approach is advantageous because the accuracy of the forward model has a large impact on the quality of the reconstruction. Learning the forward model is a daunting task for large problems because the number of elements of the matrix is the number of image pixels times the number of measurements—roughly the number of image pixels squared, easily terabytes of data. Often, knowledge of the structure of the forward model can greatly reduce the number of measurements needed to estimate it.

One example [52] comes from deconvolution of microscopy images. Using our knowledge of the underlying physics, we can assume that the system matrix is a convolution—a Toeplitz (diagonal-constant) matrix. If we further assume that the impulse response (known in microscopy as the point spread function (PSF)) is sufficiently limited in space, then the system matrix is completely specified by its middle rows. By imaging small beads, we can create a training set with images of the form ${\boldsymbol{\mathbf{f}}}_{t}=\begin{bmatrix}0&\dots&0&1&0&\dots&0\end{bmatrix}^{T}$. The resulting measurements, $\{{\boldsymbol{\mathbf{g}}}^{t}\}$, each provide an independent estimate of the PSF that can be averaged to reduce the effect of noise.

Another example [53], comes from PET where measurements of a carefully positioned point source were used to estimate the rows of the forward model. While estimating the entire forward model would require approximately four million measurements, the authors were able to use the geometrical symmetries of the system plus a physical model to interpolate the forward model from just 1,599 measurements.

A final example comes from MRI, where coil sensitivity maps are a key piece of the forward model that must be estimated from data. The standard approach is to collect data from a body coil, which has uniform sensitivity, and use it to estimate the sensitivity maps. The downside of this approach is that the SNR of the body coil data is low, making the sensitivity estimation noisy. The problem remains an area of active research [10].

5.2 Learning the Regularization Term

There are a variety of ways to use training data to improve the regularization term. One straightforward approach is to use training to optimally adjust the regularization strength, i.e., tune the parameter $\lambda$. In the approach of [54], the training takes the form of a bilevel optimization,

[TABLE]

with

[TABLE]

where $p$ is one or two. The solution of this problem is the value of $\lambda$ that gives optimal performance on the training set. Reconstructions of new images should then use this optimal $\lambda$. The same method can be extended multiple ${\boldsymbol{\mathbf{L}}}$ operators, each with their own weight; e.g., [54] uses the set of $5\times 5$ discrete cosine transform (DCT) filters.

Going a step further, we can try to learn the form of the potential function, i.e., the $\Phi$ in a regularization term $\Phi({\boldsymbol{\mathbf{L}}}{\boldsymbol{\mathbf{f}}})$. Recalling the ADMM algorithm from Section 4.4, we know that solving a regularized inverse problem involves applying the proximal operator associated with the potential function. We can implicitly specify the potential function by learning the proximal operator, e.g., by parameterizing it using 1D B-splines [55, 56]. We can see $\operatorname{prox}$ learning as a generalization of tuning the regularization weight: instead of scaling a known function, we deform a parametric function. If the learned functions satisfy mild constraints, the resulting algorithm is guaranteed to converge.

Or, as a different approach, we can learn the filters in the regularization term leading to a learning procedure [57],

[TABLE]

where the matrix ${\boldsymbol{\mathbf{L}}}$ must be constrained to lie in a set $\mathcal{C}$ to avoid the trivial solution ${\boldsymbol{\mathbf{L}}}={\boldsymbol{\mathbf{0}}}$. A variation of this approach is to allow for noisy training data, leading to

[TABLE]

Recalling the connection between the $\ell_{1}$ norm and sparsity, we can view this learning procedure as seeking a sparsifying transform for the images in the training set. Intuitively, this learned transform should allow for sparser representations (and therefore better regularization) than a fixed one.

We can also learn a sparsifying transform in the synthesis formation. Doing this is known as sparse dictionary learning (or sparse coding) [58]. To learn the ${\boldsymbol{\mathbf{L}}}$ in a synthesis formulation (4.6), we solve

[TABLE]

where $p$ is zero or one. Although for simplicity we have written the data fidelity term at the whole-image level, it is usually evaluated at the patch level, so that ${\boldsymbol{\mathbf{L}}}$ is a dictionary that can sparsely represent, e.g., $15\times 15$ images patches. These patches must then be combined, usually by averaging, to create an image. Convolutional sparse coding, which more-closely matches the formulation (5.5) has also been explored [59], as has a patch-based method with additional image-level regularization to smooth patch boundaries [60]. See [61] for an example of this approach used in X-ray CT.

Finally, we can learn the entire regularization term. Again, this can be done using a bilevel optimization [62]

[TABLE]

with

[TABLE]

where ${\boldsymbol{\mathbf{\theta}}}$ is a vector of parameters that controls both the penalty function $\phi$ and the analysis filters, ${\boldsymbol{\mathbf{L}}}$. Regularizers of other forms can also be learned, e.g., [63] use a pixel-wise autoregressive model as a regularizer.

As we saw in Section 4.4, certain algorithms for solving regularized reconstruction problems like (5.7) make use of the proximal operator of the regularization term. One very flexible method for creating new regularizers is to replace the proximal operator in these algorithms with, e.g., a CNN; this implicitly defines a new regularization term. The approach is sometimes called plug-and-play regularization. For examples of work along these lines, see [64, 65, 66, 67]. A further generalization of this idea is consensus equilibrium [68], which seeks reconstructions that balance data fidelity and regularization without necessarily minimizing an optimization functional.

Another approach [69, 70, 71], which is closer in spirit to a constrained minimization problem, is to learn a image-to-image regression function, $\mathcal{R}_{{\boldsymbol{\mathbf{\theta}}}}$, typically a CNN, that maps each ${\boldsymbol{\mathbf{f}}}$ to the nearest member of the set of plausible images, $\mathcal{S}$. This function allows us to solve

[TABLE]

which is conceptually attractive because it returns the image that best fits the measurements among all plausible images. The challenge here is that plausible is defined by the training data, which may not be sufficiently representative to capture the notion. Yet another approach [72] is to combine the ideas of sparsity and projection by a learning parametric image (or patch) encoder and decoder pair according to

[TABLE]

Reconstruction can then be achieved via

[TABLE]

5.3 Going Outside the Variational Framework

One alternative to the variational formulation that has been intensely explored in the last few years is to perform a linear reconstruction followed by a learned image-to-image regressor, typically a CNN. The linear reconstructions often result in heavy artifacts when the number of measurements is low (streaks in CT or aliasing in MRI); the job of the regressor is to remove these artifacts; therefore the method is sometimes referred to as learned denoising or artifact removal. Learning takes the form

[TABLE]

where ${\boldsymbol{\mathbf{H}}}^{\dagger}$ is a matrix that is an approximate inverse of the forward model, e.g., back projection, ${\boldsymbol{\mathbf{H}}}^{*}$, or filtered back projection, ${\boldsymbol{\mathbf{W}}}{\boldsymbol{\mathbf{H}}}^{*}$ or ${\boldsymbol{\mathbf{H}}}^{*}{\boldsymbol{\mathbf{W}}}$, for some appropriate matrix ${\boldsymbol{\mathbf{W}}}$. The same idea can be used with mixtures of linear (or nonlinear) reconstructions, with the idea that the regressor will learn to spatially vary the regularization strength based on the image content [73]. The main advantage of these methods, as compared to fully end-to-end methods, is that the design of the regressor is simplified: because it always acts on images (rather than vectors of measurements), it can can employ well-studied images processing tools such as multiscale processing and convolutions. This approach has been implemented with a variety of network architectures for X-ray CT [74, 75, 76, 77] and MRI [78] reconstruction.

Example 9** (FBPConvNet).**

As an example of some of the practical considerations that go into designing and training a learning-based method, we present details of one such method, FBPConvNet [74], which performs X-ray CT reconstruction. The FBPConvNet reconstruction method consists of FBP followed by a CNN. The FBP part is a standard algorithm which can be performed with MATLAB’s iradon function. On the other hand, special software [79] is used to define the CNN’s structure and train it efficiently. The structure of the CNN comes from [80], with modifications to improve the training stability and performance. Designing the architecture of the CNN is, at this point, an art and requires some trial and error.

In a simulation experiment, the training set consisted of 500 pairs of low-dose measurements and the corresponding ground truth reconstruction. These data were augmented by flipping each ground truth image in the horizontal and vertical directions, resulting in 2,000 training pairs. The training itself is a challenging optimization problem that is solved using stochastic gradient descent (SGD)—a version of gradient descent that approximates the gradient with a few training samples at a time. The SGD algorithm requires tuning of several parameters for good convergence. Among these are the learning rate, which is the step size in gradient descent; the learning rate decay, which is a process for the learning rate to decrease after some iterations to improve convergence; and the batch size, which is the number of training points to use when computing a stochastic gradient. In the case of FBPConvNet, 101 iterations through the training data were sufficient for good performance and could be completed in about 15 hours using a graphics processing unit (GPU). After training, the algorithm was run on 25 unseen images to assess its performance.

We show the results of this experiment in Figure 5.1, with a comparison to FBP and a total variation (TV)-based approach (with regularization parameters chosen to maximize performance on the testing set). The results are typical in that the learning-based method outperforms the variational one, both quantitative and qualitatively, giving slightly sharper-looking reconstructions.

A common extension to the learned denoising approach is to also learn a data fidelity term. This can be done via perceptual loss [81], where images are used as inputs to a pretrained CNN and the distance between them is computed as the Euclidean distance between their intermediate feature vectors; the intuition being that a network trained, e.g., for object recognition, should building perceptually meaningful intermediate representations of the images it operates on. Alternatively, adversarial loss can be employed. In this case, an image-to-scalar regressor, called a discriminator, is trained to give low values to real (in the training set) images and high values to fake (reconstructed images). The discriminator is trained alongside the reconstruction algorithm, pushing it to create ever more plausible reconstructions. See [82, 83, 84] for examples of these approaches for biomedical image reconstruction, as well as [85] for a comparison of different loss functions for image restoration.

Finally, we can consider learning the entire reconstruction pipeline end-to-end,

[TABLE]

Here, all of the work goes into designing the structure of $\mathcal{R}_{{\boldsymbol{\mathbf{\theta}}}}$ and finding sufficient training data. One approach is to be problem agnostic, e.g., using a fully-connected neural network [86], or a combination of a fully-connected neural network and a CNN [87]. Another is to include some fixed linear layers that are problem-dependent [74, 88].

An alternative design strategy, called unrolling, is to build the network structure from a preexisting reconstruction algorithm by turning its iterations into layers in the network. In effect, the network is a traditional algorithm for solving an inverse problem, except that some or all of the linear operations are learned (rather than being specified by ${\boldsymbol{\mathbf{H}}}$) and the number of iterations/layers is usually low, less than ten, to prevent the number of parameters in the network from becoming too large. For example, [74] notes that when ${\boldsymbol{\mathbf{H}}}^{*}{\boldsymbol{\mathbf{H}}}$ and ${\boldsymbol{\mathbf{L}}}$ are both convolutions, the typical iterate solutions to the $\ell_{1}$-regularized variational problem have the form of a CNN. This equivalence allows some learned denoising algorithms to be interpreted as unrolled iterative methods. Reference [89] provides another perspective, proposing to also learn the nonlinear activation functions inside the neural network as a way to learn the regularization functionals, and [90] explores the idea of parameterizing the ADMM algorithm. An intermediate approach is to allow ${\boldsymbol{\mathbf{H}}}$ to remain fixed, but to learn, e.g., the gradient step in unrolled gradient descent [91] or the proximal operators in an unrolled primal-dual algorithm [92]. Reference [93] explores several unrolled architectures that vary in terms of how many of the parameters are learned, concluding that the best performance comes from a balance approach (i.e., learning some, but not all, of the possible parameters).

The main advantage of learning most or all of the reconstruction pipeline is that the results can be excellent, outperforming a well-tuned $\ell_{1}$ result both quantitatively and qualitatively. Depending on their specific architecture, the algorithms can be orders of magnitude faster than iterative methods. The drawback is that training the algorithms is laborious: training can take days and is not guaranteed to converge, meaning that it must be repeated to find suitable parameters (including training parameters as well as the network architecture itself). These problems are compounded by the fact that the network may need to be trained again if aspects of the imaging scenario, e.g., noise level or body part, change [94]. Another drawback is that as more of the reconstruction pipeline is learned, the data term seems to disappear more and more, meaning that the returned solution may not be consistent with the measurements. On one hand, this mismatch may be advantageous: if the forward model is inaccurate, enforcing data consistency may make solutions worse (farther from the ground truth). On the other hand, reconstructions that are plausible-looking without explaining the measurements are hardly useful images. One way of enforcing consistency is to mix the learned reconstruction with a conventional one [95]

[TABLE]

5.4 Other Designs

There are many more designs that do not fit well into the previous sections. We share a few here to give a sense of the diversity of the field.

In many problems, we might be able to access high-quality measurements, ${\boldsymbol{\mathbf{g}}}_{t}^{\text{HQ}}$ at the cost of a longer or more-costly acquisition, e.g., by taking longer X-ray exposures or sampling more of the $k$-space in MRI. The goal is to use these measurements to learn to reconstruct from low-quality (low dose, fast) measurements. This situation is different from the standard formulation because now the training measurements are the same as the measurements we aim to reconstruct from. One solution is simply to generate a training pair $({\boldsymbol{\mathbf{f}}}_{t},{\boldsymbol{\mathbf{g}}}_{t})$ by performing a reconstruction (usually linear) on ${\boldsymbol{\mathbf{g}}}_{t}^{\text{HQ}}$ and by downsampling ${\boldsymbol{\mathbf{g}}}_{t}^{\text{HQ}}$. But, ${\boldsymbol{\mathbf{g}}}_{t}^{\text{HQ}}$ can also be used explicitly to learn how to inpaint missing measurements [96]

[TABLE]

where $\mathcal{D}$ is a model of measurement degradation (e.g., downsampling) and $\mathcal{R}_{\theta}$ is a parametric inpainting function. Then, to reconstruct from low-quality measurements, we apply $\mathcal{R}_{{\boldsymbol{\mathbf{\theta}}}}$ to generate high-quality measurements and reconstruct using a conventional method. This is sometimes called data domain learning and has been used for metal artifact removal in X-ray CT [97]. Another, similar, variation on the setting is learning to regress from a reconstruction from one measurement type to a reconstruction from another, for example, to infer one MRI scan type from another [98, 99].

Another closely-linked topic that does not exactly fit the paradigm is blind or semiblind image reconstruction. These involve forward model learning without training data. Although the problem is superficially impossible—recover ${\boldsymbol{\mathbf{f}}}$ from ${\boldsymbol{\mathbf{g}}}={\boldsymbol{\mathbf{H}}}{\boldsymbol{\mathbf{f}}}+{\boldsymbol{\mathbf{n}}}$ with unknown ${\boldsymbol{\mathbf{H}}}$—it can be solved provided enough constraints on ${\boldsymbol{\mathbf{H}}}$ and prior knowledge about ${\boldsymbol{\mathbf{f}}}$. The archetypal example is blind deconvolution [100], where, as with PSF estimation from training, the fact that ${\boldsymbol{\mathbf{H}}}$ is a convolution is key. Many learning algorithms can be adapted to work without training in a similar way, by alternating a reconstruction step with a learning step that uses the reconstructed images themselves as training; the challenge is that such a procedure can easily diverge and give nonsensical results.

5.5 Where to Get the Training Data

A key ingredient in any learning method is the training data. But, how can we acquire training data for inverse problems without already having a working reconstruction algorithm to make images from measurements? If we have such a system, why do we need to design another one? The easiest route out of this chicken-and-egg problem, as we already mentioned in the previous section, is to work in the low-quality (few, noisy measurements, sometimes called compressed-sensing) regime, and to exploit the fact that we can access high-quality measurements (and reconstructions) at the cost of longer scans. We can then simulate low-quality measurements by degrading the high-quality ones. Thus we train a low-quality reconstruction system from a high-quality one.

But, this is not the only way. In [101], the authors acquire training data by using both a low and high numerical aperture microscope, thus removing the need to simulate low-quality measurements. The cost of this more-realistic training data is that the images from both microscopes needs to be correctly registered before learning. Training data can also come from a highly-accurate physical simulation [98] that may be impossible to use directly as a forward model because of its computational cost or nonlinearity.

Another way is to build algorithms that do not require a training set, or, at least, not a paired training set (i.e., measurements and their corresponding reconstructions). One way to do this is by enforcing cycle consistency [102], where a pair of algorithms are built, one for reconstruction and one for simulating data. These algorithms are trained together from unpaired data such that they are inverses. This approach is similar to generating training from simulation, except that, here, the simulator itself is learned from data and, therefore, may be more realistic.

5.6 Summary

Learning-based reconstruction involves using a training set to tune a parametric reconstruction algorithm. There is a huge variety in the specific architecture of these algorithms, with some being closely related to direct or variational methods and others not. In any case, creating a suitable training set is a key challenge.

5.6.1 Further Reading

The recent special issue [103] focuses on original research on learning-based tomographic reconstruction. There are several recent reviews around the topic of learning for image reconstruction: [104] focuses on biomedical imaging, [105] on general image processing, [106] on MRI reconstruction, and [107] on computer vision and medical imaging. For a long-term perspective on learning in imaging, see [108], and for a look at how these learning approaches might affect the practice of radiology, see [109]. For another clinical perspective, [110] provides a double-blind study in which radiologists compared learning-based methods to commercial iterative reconstruction methods; the learning-based methods were shown to be at least as good as the commercial ones.

Chapter 6 Conclusion

In this tutorial, we have given a roughly chronological overview of biomedical image reconstruction. We began with a toolbox of mathematical operators that can be used to build models of many physical imaging systems. We then showed how to use these forward models to solve reconstruction problems, either via direct inversion or variational formulations, with a focus on the paradigm of sparsity-promoting regularization. We ended by covering some of the many ways that training data can be used to develop new reconstruction methods via machine learning. Here, we provide a brief comparison of these methods, commenting on the strength of their theoretical underpinnings, their complexity, speed, performance, and robustness. We conclude with comments on future directions for the field.

6.1 Comparisons

Our theoretical understanding is strongest for classical reconstruction methods and weakest for recent learning-based methods. The theory behind both classical reconstruction algorithms and many sparsity-based ones is, for the most part, settled. Thus, we can often make definitive statements about these algorithms, e.g., we know that certain iterative procedures converge to global optimums. We can sometimes even state that a method is optimal for a certain class of problems. That said, there is plenty of work to be done, especially in pushing our understanding to ever more complex and realistic data models. There has been significant progress in explaining the strong performance of learning-based methods on many imaging processing tasks, including work from the perspective of approximation theory [111], unrolling [112], and invariants [113]. Despite these efforts, it seems that theory will continue to lag behind practice for years to come. We should accept the possibility that learning approaches may never be understood as fully as we might like, and, further, that understanding may become increasingly secondary to performance. The fact that the best-understood methods are not the top performers is not a new state of affairs in imaging, but it is understandably disturbing to see the gap widen; see [114] for a clear expression of the sentiment.

The trend in complexity among the various reconstruction methods is the same: classical methods are simple to implement and tune, while learning-based methods are notoriously hard to train properly. For classical methods, the number of parameters is small and can sometimes be calibrated, e.g., by estimating the measurement noise. Superficially, learning can remove the problem of parameter selection by fixing parameters during training. However, training requires its own set of parameters: the initialization, gradient step size, the network architecture, etc. Because the training problem is usually nonlinear and nonconvex, all of these parameters affect the result. However, in a practical settings, such as industrial or medical imaging, it can be worth investing the engineering effort to train the system.

In terms of speed, conventionally we think of trading speed for reconstruction quality: direct reconstruction is fast but gives the worst results, increasingly sophisticated regularization improves these results at the cost of increased runtime. Learning-based methods upend this order. Most of the runtime is spent during training, which happens only once, after which reconstructions are much quicker. In fact, this happens because only network architectures that allow very fast evaluations can be trained at all. It is this mode of thinking—fix the computational budget, then optimize the algorithm within it—that underlies the unrolling concept [112], and it is a valuable one.

In terms of reconstruction performance, direct methods are excellent when the number of measurements is high and the noise is low; regularization and learning only make sense in the low-dose (high-noise) regime. (In the full-dose regime, it is not clear how end-to-end learning can be used at all because the training data itself would need to come from some reconstruction algorithm.) In the low-dose regime, however, regularization or learning is key, and learning-based methods show the best performance on almost every reconstruction benchmark. One perspective is that (when sufficient training is available) learning-based methods provide an upper bound on performance, which variational methods might achieve if their regularization and data terms are perfectly suited to the task at hand. One example of this is in [74], where the variational approach beats the learned one when the data fits the regularization model (TV-based reconstruction of piecewise constant images).

In terms of robustness, learning-based algorithms are currently the weakest. As we discussed in the last section, slight changes in the imaging parameters (noise level, body part, etc.) can severely reduce the performance of the network, or, worse, lead to the appearance of realistic-looking structures in the image that are not supported by the measurements at all. Also, in the context of classification, there is the concept of adversarial examples [115], which are images that are designed to fool a neural network. While there is debate about why adversarial examples can be created and what their existence means, they might indicate that a learning-based image reconstruction algorithm will occasionally (and inexplicably), give a spurious result. Recent work has shown examples of this effect in image reconstruction [116]. For direct and variational methods, on the other hand, the reconstruction quality is consistent across datasets and (at least in the case of direct and convex variational methods) degrades gracefully as the measurements degrade. In particular, the degradation has a predictable form: noise or specific artifacts (e.g., aliasing in MRI, streaks in CT).

6.2 Future Directions

The future of the field remains to produce the highest-quality image using a fixed imaging budget, e.g., scanning time or radiation dose. Doing this will continue to require realistic physical models, though these models may play a new role in reconstruction. For example, they might be used to create training data, or they themselves might be learned from training data. The focus on signal modeling will remain, and should be pursued both with sophisticated mathematical modeling and with data-driven approaches. Finally, we will always seek practical algorithms with few hand-tuned parameters.

Specifically in the context of biomedical imaging, we advocate for a shift towards task-oriented evaluations. The images that we reconstruct in the biomedical domain are usually used for a specific, predefined purpose, e.g., to grade a tumor or to measure the quantify the effect of a genetic manipulation on a model organism. In this setting, the traditional measures of reconstruction performance (SNR, mean squared error (MSE), structural similarity index (SSIM)) are insufficient and should be replaced with task-oriented evaluations. Such measures may be expensive (if they require human experts) and they may reveal difficult truths (e.g., that even a few dB improvement in SNR is meaningless), but they will be critical in creating algorithms that truly advance the state of biomedical imaging in practice.

Ultimately, future reconstruction algorithms should employ educated learning: a careful fusion between model-based approaches and data-driven ones. We believe that this approach will provide both excellent performance and the robustness and performance guarantees needed for biomedical applications. In biology, this will enable high-quality imaging of a wider range of biological structures and processes; in medicine, this will allow diagnostically useful images to be created more quickly, cheaply, and with less harm to the patient.

Acknowledgements

This work was made possible thanks to the CIBM Center for Biomedical Imaging, founded and supported by Vaud University Hospital Centre (CHUV), University of Lausanne (UNIL), Swiss Federal Institute of Technology Lausanne (EPFL), University of Geneva (UNIGE), University Hospitals of Geneva (HUG) and the Leenaards and the Louis-Jeantet Foundations and thanks to the European Research Council through the European Union’s Horizon 2020 Research and Innovation Program under Grant 692726, GlobalBioIm. The authors thank Philippe Thévenaz, Julien Fageot, Avrajit Ghosh, and Pol del Aguila Pla.