3.1.

Forward Modeling

The forward model in DOT describes the propagation of near-infrared light through biological tissues from a distribution of light sources illuminating the surface, to a collection of detectors measuring the light exiting the tissue. The sources can either consist of optical fibers delivering light to a small skin area, or of devices that illuminate a larger area of the surface possibly using structured light patterns.¹⁰ Data acquisition systems may employ continuous sources (continuous wave imaging), periodic at MHz to GHz frequencies (frequency-domain imaging) or pulsed (time-domain imaging). At the wavelength used in DOT, tissue is highly scattering, and photon transport can often be described as a diffusion process. For the frequency domain case, the forward model is then given by

Analytic solutions of Eqs. (1) and (2) only exist for simple geometries. For practical applications, a numerical model must be employed that can incorporate inhomogeneous parameter distributions and complex boundaries. The Toast++ package uses a Galerkin FEM to implement the DOT forward problem.

3.2.

Time-Domain Modeling

As an alternative to measuring a modulated signal, time-domain data acquisition systems use short pulses of light $Q(\mathbf{r},t)$ in the picosecond range for sources, and measure the temporal dispersion of the transilluminated signal at the detector sites. The light transport model in this case is the time-dependent diffusion equation,

Numerically, the time derivative term can be approximated by a finite-difference scheme. The Toast++ software allows the simulation of the time-dependent measurement signals, using either implicit or explicit schemes for propagation in time. In addition, Toast++ allows the direct computation of the temporal moments, as well as the Laplace transform and its moments (the Mellin–Laplace transform) of the measurement signal, without the need for a full computation of the temporal profile.^12,13

3.3.

Inverse Problem Solver

The Toast++ suite contains the building blocks to construct an iterative inverse solver for reconstructing the distribution of model parameters $\mathbf{x}=\{{\mu }_a,\kappa \}$ or $\mathbf{x}=\{{\mu }_a,{\mu }_s\}$ from boundary measurements $\mathbf{y}={\langle M,J\rangle }_{\partial \mathrm{\Omega }}$ for some measurement profile $M$. The inverse problem is defined by a regularized least-squares approach, where an objective function $\mathrm{\Psi }$ is minimized,^1,14 given by

Toast++ provides a range of regularization options, including Tikhonov and total variation (TV) methods, to enforce different smoothness conditions in the solution. In addition, spatially varying regularization parameters can be applied to take into account prior knowledge about the internal structure of the target.

Toast++ provides functions for constructing the direct and adjoint fields for a given set of source and detector locations, and it allows the solution of the adjoint problem either by explicitly constructing the Jacobian and Hessian matrices of the forward operator, or by an implicit, matrix-free approach. Using this functionality, gradient-based inverse solvers can be constructed by making use of first derivatives, such as nonlinear conjugate gradients (NCG),¹⁵ or second derivatives such as Gauss-Newton (GN) or Levenberg–Marquardt approaches.¹⁴

4.1.

Matrix Library

The Toast++ matrix library (libmath) contains a class hierarchy for sparse and dense matrix and vector types, as shown in Fig. 2. It exposes a common matrix interface via the abstract Matrix class, which hides the implementation details of specialized matrix subtypes such as dense, compressed sparse row (CSR), and coordinate-indexed matrix classes, from the application layer. The matrix library is templated and can be instantiated for real and complex-valued data types.

Fig. 2

Matrix class hierarchy in the Toast libmath library.

Public methods include element, row, and column access via iterators, matrix-vector products, as well as direct and iterative linear solvers including conjugate gradients (CG), generalized minimum residual methods, LU, and Cholesky decompositions. Many of the methods are implemented as wrapper functions to external libraries, such as the SuperLU library for LU solution or preconditioning of sparse matrix classes.

The MATHLIB library also contains a distributed sparse matrix class (TCompRowMatrixMPI) using the MPI protocol, which allows the matrix data to be distributed over multiple processor nodes, while maintaining the common matrix interface, and thus being transparent to the application programmer. This matrix class can be used to represent the linear system of the finite-element model in conjunction with a partitioned distributed mesh class, as discussed in the next section.

4.2.

Finite-Element Library

Using an $N$-dimensional polynomial basis expansion with local support $u_i(\mathbf{r})$, $i=1\dots N$ over domain $\mathrm{\Omega }$, a piecewise polynomial approximation ${\varphi }^h$ of field $\varphi$ can be expressed as

The realization of the integrals in Eq. (8) depends on the type of element and the shape functions used. In some cases, such as triangles or tetrahedra with polynomial basis functions, analytic integration formulae may be available. Otherwise, numerical quadrature rules must be employed.

The Toast++ finite-element subsystem (libfe) contains the mesh and element classes for defining a numerical representation of the forward problem. The Mesh class represents the discretization of the computation domain and it is a container for a list of node coordinates and a list of elements. It provides functions for returning the sparsity pattern of the system matrices for the linear FEM system, and for populating the matrices by performing the appropriate integrals over elements and mapping the results from the local element matrices to the global system matrices. Toast also contains methods for node reordering, e.g., for bandwidth minimization or minimizing the fill-in of the system matrix decomposition for direct solvers and preconditioners.

Toast++ provides support for different element types and shape functions in two and three dimensions. In 2-D, 3-, 6- and 10-noded subparametric and isoparametric triangles are available, providing linear, quadratic and cubic polynomial shape functions, respectively. In addition, Toast provides 4-noded rectangular elements for representing structured meshes. In 3-D, Toast offers 4- and 10-noded sub- and isoparametric tetrahedra, as well as wedge and regular voxel elements. The element class hierarchy is shown in Fig. 3. The element types are represented in a hierarchical class structure, and expose a common interface via the abstract Element base class. This includes methods for mapping from local to global coordinates, Jacobian computation, as well as integrals of shape functions, shape function derivatives and combinations thereof, over the element volume.

Fig. 3

Toast element class hierarchy. Abstract classes are shown in gray, nonabstract classes in yellow together with graphical shape representations.

Elements use analytic integration rules where available, including subparametric triangle and tetrahedron elements, and resort to numerical quadrature rules for isoparametric elements with curved surfaces. The implementation details of the element integrals and other methods are hidden from the application layer by use of polymorphism. It is, therefore, possible to write applications independent of element type, and to define meshes with mixed elements.

Toast++ also contains a distributed mesh class (MeshMPI) which can be used to automatically partition a mesh and distribute it over multiple processor nodes. Toast provides the possibility to link to the Zoltan partitioning library¹⁷ for partitioning and node migration. In combination with Toast’s distributed sparse matrix support, the linear FEM system can be assembled and solved transparently to the application layer. This allows the deployment of Toast for large-scale problems on computer clusters without additional coding effort from the user.

Different options are available for the choice of source distributions $Q$. Toast supports point sources and extended source profiles with a Gaussian or cosine shape to simulate optical fiber-based light delivery systems.¹⁸ In addition, Toast also allows modeling of illumination of large areas of the surface with shaped light patterns.¹⁰ Source patterns are projected onto the surface of a mesh of arbitrary shape by using an instance of one of Toast’s projector classes, which support either pinhole or orthographic projection. The projectors are implemented with the use of OpenGL routines via the Mesa-3-D library which can also make use of graphics hardware acceleration. MATLAB bindings for the projector functions are supported in the Toast toolbox.

The same projectors can be used to map a distribution of transilluminated light from the surface to an imaging system such as a charge coupled device (CCD) camera. Figure 4 shows an example of a source pattern projected onto a cylindrical surface, and the resulting surface exitance on the opposite side of the cylinder mantle projected back onto a CCD camera. The use of shaped light illumination and CCD detectors can reduce measurement times significantly.

Fig. 4

Example of a source pattern (top left) projected on a cylinder with Toast’s pinhole projector operator (top right), and resulting surface exitance (bottom right) projected back into a CCD detector array positioned opposite the source projector (bottom left).

4.3.

Alternative Numerical Modeling Schemes

In addition to the FEM scheme described above, Toast also supports alternative numerical subsystems for simulating light propagation in scattering media. The software package includes a discontinuous Galerkin (DG) discretization scheme,¹⁹ which can be used to accurately simulate discontinuities in the internal optical parameter distributions, including the refractive index. In addition, meshes for DG computation can be refined more easily than those used in standard FEM methods. DG, therefore, supports efficient adaptive refinement methods, e.g., as a function of the gradient of the evolving parameter distribution during a reconstruction.

Toast also contains a modeling module using the boundary element method (BEM),²⁰ which allows efficient simulation of piecewise constant parameter distributions. This scheme has been applied to shape reconstruction of internal boundaries of regions with piecewise constant optical parameter values. BEM can also be used to reduce the complexity of reconstruction by restricting the solution on interface layers, such as the cortical surface.²¹ Toast has been used to combine FEM and BEM models for DOT reconstruction in layered media, where superficial tissue layers are modeled with BEM, while the region of interest is reconstructed with a spatially resolved FEM approach.²²

4.4.

DOT Inverse Solver Library

Low-level support for building an inverse solver for parameter estimation in DOT is provided by the libstoast library. The library contains a class hierarchy of basis mapping operators that allow a field or parameter distribution on $\mathrm{\Omega }$ to be mapped from the mesh basis to an independent basis used for the reconstruction. Given the mesh basis $\{u_i(\mathbf{r})\}$, $i=1\dots N$, and inverse solution basis $\{v_j(\mathbf{r})\}$, $j=1\dots M$, the mapping from nodal coefficients ${\mathrm{\Phi }}_i^{(u)}$ to inverse basis coefficients ${\mathrm{\Phi }}_j^{(v)}$ is given by²³

Toast supports a range of inverse basis functions $v$, including linear and cubic voxel bases, as well as blob bases with radially symmetric basis functions of different shape, such as truncated Gaussian and Kaiser–Bessel window functions.²³ Figure 5 shows examples of the radial profiles of the supported blob basis functions.

Fig. 5

Radial profiles of blob basis types supported by Toast.

The separation of forward and inverse basis allows the reconstruction to be tailored to the required performance and convergence criteria, while conserving the fidelity of the forward solver. For example, it is then possible to adaptively refine the mesh without also affecting the reconstruction component.

libstoast contains support for defining the parameter estimation problem by minimizing the regularized least-squares cost functional Eq. (5). Methods for computing direct and adjoint fields are provided by linking to the forward solver module. Computation of the Jacobian matrix of the forward operator for gradient-based methods is supported as well as matrix-free methods. MATLAB and PYTHON bindings for the inverse solver utilities are provided, allowing the iterative optimization loop to be defined in application space directly by the user.

4.5.

Support for Parallel Computation

Toast++ can be used on different parallel architectures for performance improvement. It contains suppport for shared memory systems by using threaded solutions for the available Krylov solvers, using a master-worker structure to distribute multiple right-hand sides of the linear FEM problem over multiple threads.⁷ This high-level parallelization strategy is easy to implement and has little computational overhead, resulting in good scalability, in particular where the FEM problem is to be solved iteratively for a large number of right-hand sides, as is typically the case in DOT. When compiled with thread support, parallel versions of time-critical routines, such as forward solution and gradient computation, are also available from the MATLAB and PYTHON interfaces. Figure 6 shows timings for FEM forward solution with a BiCGSTAB iterative linear solver (tolerance ${10}^{-14}$) on cylindrical meshes of different resolution, using a threaded implementation with 1, 2, 4, 8, and 16 threads on a Xeon E5-2665 workstation with 16 2.4 GHz processors. When using the MATLAB and PYTHON interfaces, no additional setup code is required to make use of the multithreaded Toast function support. The number of processors can be adjusted manually, e.g., in MATLAB with

Fig. 6

Timing comparsions for FEM forward computations of 64 source positions and three different mesh resolutions using a BiCGSTAB iterative solver with thread support.

A graphics processor hardware accelerated version of the FEM forward solver for Toast has been developed using the CUDA framework for NVidia graphics processors, and the CUSP library²⁴ for sparse linear system computation on the graphics processor. The Toast CUDA implementation moves the solution of the linear system to the graphics hardware after the sparse system matrices and right-hand sides have been assembled. This approach is particularly efficient for the iterative solution of the time-domain problem, because the stiffness and mass matrices in Eq. (8) remain unchanged and can be reused for each iteration after having been copied once to the device memory. We have reported speed improvements of a factor of up to 8 for a frequency-domain solution, and up to 13 for a time-domain solution on an NVidia GTX 285 GPU, compared to a single-threaded CPU solution on an Intel Xeon processor clocked at 2.0 GHz with 4 MB cache and 12 GB main memory.⁸

4.6.

MATLAB Bindings

The functionality of the Toast core libraries is exposed to the MATLAB scripting environment as a toolbox containing a collection of compiled mex files and utility scripts that are available to MATLAB as callable functions. The MATLAB interface to Toast makes use of an object-oriented approach, by representing meshes, elements, solvers, etc., as MATLAB objects. Each object has properties and methods that closely resemble the class definitions of the underlying C++ code. As an example, a subset of methods of the toastMesh class is shown in Fig. 7. The interface consists of a single mex file that handles the communication between the MATLAB Toast toolbox and the C++ libraries. For example, the ToastMesh class contains a handle to a C++ mesh object, and a series of methods for computing mesh geometry and transformations, matrix assembly, extracting element data, surface integrals, and so on. Each of these methods is implemented by calling the mex file to pass a request to the C++ Toast libraries. The libraries process the request, and pass the results back to the MATLAB method. In this way, the efficiency of the Toast library is maintained, while allowing the user the convenience of the MATLAB script interface.

Fig. 7

Methods in the ToastMesh class of the Toast MATLAB toolbox.

4.7.

PYTHON Bindings

Similar to the MATLAB toolbox, the PYTHON script interface is provided by compiled PYTHON modules that can be imported by a PYTHON script. The functionality of the MATLAB and PYTHON interfaces is similar. Sharing of the underlying core libraries minimizes code duplication and improves stability and simplifies the process of code generation. Toast functionality is available in PYTHON via a plugin module that links to the Toast core libraries. Toast makes use of the numpy and scipy modules for the sparse matrix data structures and linear solvers. The interface is similar to the MATLAB Toast toolbox, with methods provided for constructing or reading meshes, defining the FEM problem, assembling, and solving the linear system. An additional advantage of the MATLAB and PYTHON interfaces is the availability of additional functionalities, such as sparse matrix solvers and visualization tools.

The Toast toolbox functions in both interfaces can be used as building blocks for creating scripts for solving both the forward and inverse problems. Section 5 contains script examples for different types of problems.

4.8.

Data Structures and Efficiency

Both MATLAB and PYTHON support their own data structures to represent matrix objects which are then exposed to external modules via a programming interface. To avoid replication and copying of data between MATLAB or PYTHON and the Toast core libraries, Toast allows the construction of matrix and vector objects that link to external data buffers. This allows the efficient direct manipulation of MATLAB and PYTHON objects while maintaining a consistent interface in the Toast matrix and vector class structure. The following code example from a PYTHON module shows the construction of a Toast vector that directly links to a PYTHON vector:

Subsequently, vec behaves like an ordinary Toast vector, but shares its data with the py_vec PYTHON object. Dense matrices are wrapped in a similar way. For sparse matrices, Toast uses a CSR format which is also supported by PYTHON with the scipy module. Scipy does not have a C interface; therefore, the data and index arrays are extracted from the PYTHON object at script level, and passed separately to the interface module

Where functions 1 and 2 are the corresponding function in the toast PYTHON module. The CSR matrix construction in their implementation would then be

Where py_mat_data, py_mat_rp and py_mat_ci are the representations of the CSR matrix arrays passed to the function. Again, this construct avoids the duplication and copying of data between matrix representations.

MATLAB only supports the compressed sparse column matrix format. While this could be accommodated by a transposition flag in the Toast sparse matrix class, MATLAB in addition uses a storage format for complex data types that is incompatible with the Toast format, necessitating the conversion of data structures at the Toast-MATLAB interface, which makes the Toast-MATLAB interface slightly less efficient than PYTHON.

4.9.

Visualization Tools

Both MATLAB and PYTHON provide visualization tools that can be utilized in connection with Toast functionality to display mesh geometry and results. MATLAB provides a toastShowMesh function to display 2-D and 3-D meshes and nodal funtions interpolated over the mesh area or surface. In addition, Toast can write out mesh geometries and nodal solutions in visualization toolkit (VTK) format, which can be visualized in a variety of standard tools. Figure 8 shows the surface of a head mesh exported from Toast and displayed in the MayaVi2 Data Visualizer.

Fig. 8

Head mesh exported to VTK format and displayed in MayaVi2 Data Visualizer.

5.1.

DOT Forward Solver

We consider the problem of computing the complex photon density field $\varphi (\mathbf{r},\omega )$ in Eq. (1) in a complex head mesh using the Toast-MATLAB toolbox. A schematic of the data components comprising the forward problem is shown in Fig. 9, including the definition of the mesh geometry from node coordinates and element connectivity graph, the definition of source and boundary projection operators and the nodal parameter distribution. The forward solver uses this information to assemble the system matrices of the FEM problem, the right-hand sides and boundary operators. Using an appropriate preconditioned linear solver, the resulting fields and boundary data can be computed.

Fig. 9

Layout of the Toast forward solver. The mesh object, together with parameter distributions provided as nodal coefficients, source and boundary conditions is used to construct the FEM linear system.

The problem geometry is represented in a Mesh class which acts as a container for element and node coordinate lists, it maintains the mapping from individual element to global DOF, and provides methods for querying the sparsity structure of the global system matrices and their assembly. Toast only contains limited support for mesh generation, but can import meshes generated from external meshing tools, which can then be saved in Toast format. Access to the mesh geometry and integration methods in the MATLAB-Toast framework is provided via the toastMesh class. A mesh object may be loaded from a file, or constructed on the fly, and then can be queried with the Data method to extract node coordinates and element connectivity graph:

The Data method retrieves the node (nd) and element arrays (elidx), together with an element type list (eltp). The Display command shows the mesh geometry in a 3-D surface plot. An example mesh view for a head mesh with internal structure is shown in Fig. 10.

Fig. 10

MATLAB mesh viewer: Visualization of region interfaces for skin, skull, gray and white matter of a volume mesh of the head.

To perform the system matrix assembly for the DOT forward problem in Eq. (7), we now define an empty sparse matrix to hold the DOT stiffness matrix, and loop over all mesh elements to compute the element integrals in Eq. (8), before mapping them to the global DOF:

Where mua, kap, and zeta are the nodal coefficient vectors of the parameter distributions of the DOT forward problem, and omega is the modulation frequency. Instead of computing the integrals per element and assembling them manually into the global stiffness matrix, Toast also allows to assemble each global matrix component in a single step with the SysmatComponent method, performing the mapping from local to global DOF in the library core:

As an additional level of abstraction, Toast also provides a convenience function, dotSysmat that simplifies the matrix assembly specifically for the DOT problem to a single function call:

where mua, mus, and ref are the nodal absoption, scattering, and refractive index coefficient vectors.

The right-hand sides are assembled from a predefined list of source positions read from a data file:

where qvec is a sparse matrix representing multiple boundary source distributions as column vectors. The complex fields for all sources can now be computed with a MATLAB linear solver, such as GMRES (generalized minimum residuals):

As an example, the logarithmic amplitude on the interfaces of a four-layer head model, arising from a modulated point source on the skin surface, is shown in Fig. 11. The fields can now be projected to measurement data of boundary exitance by applying the measurement operators, which can then be displayed as source-detector maps, as shown in Fig. 12:

Fig. 11

Visualization of logarithmic amplitude on the skin surface and the interfaces between skin/skull, CSF/gray matter and gray matter/white matter arising from a point source at the outer surface (indicated by a blue dot).

Fig. 12

Top: sinograms of relative logarithmic amplitude (left) and phase boundary data differences (right) between background parameters and an inclusion of increased absorption (“haemorrhage”) in the cortical layer. The images in the bottom row show the projections onto the surface of the field differences from a single source opposite the inclusion. (An animation showing the projected inclusion from all sources is available online.)

To demonstrate the use of the Toast-PYTHON module, the equivalent PYTHON script is shown in the following listing. Both script interfaces link to the same Toast core libraries, and therefore produce equivalent results.

For ease of use, the MATLAB toolbox also provides a graphical interface for providing the input parameters to the DOT forward solver, as shown in Fig. 13.

Fig. 13

Example for a MATLAB graphical interface to the DOT forward solver.

5.2.

Iterative Parameter Reconstruction

The inverse problem for the DOT problem consists of recovering the spatial distributions of the optical parameters ${\mu }_a(\mathbf{r})$ and $\kappa (\mathbf{r})$ [or equivalently, ${\mu }_a(\mathbf{r})$ and ${\mu }_s(\mathbf{r})$] of the model Eq. (1) in domain $\mathrm{\Omega }$ from boundary measurements of light transmission at the surface $\partial \mathrm{\Omega }$. The iterative approach provided by the Toast suite solves this nonlinear problem by finding the minimizer $\{{\mu }_a,\kappa \}=\widehat{\mathbf{x}}=\arg {\min }_{\mathbf{x}}\mathrm{\Psi }(\mathbf{x})$ to the objective function defined in Eq. (5).

Toast provides the building blocks to minimize $\mathrm{\Psi }$ with a choice of iterative optimization schemes, such as NCGs, or a GN approach.¹⁴ This includes a method for the computation of the direct and adjoint fields (toastFields), the gradient of the objective function $\mathrm{\Psi }$ with respect to the model parameters $x$ (toastGradient), the Jacobian of the forward model (toastJacobian), and an inexact line search for obtaining a step length in a given search direction (toastLineSearch).

Further, Toast provides a selection of regularization schemes, such as Tikhonov and TV, for implementing the regularization term $\mathcal{R}$ in Eq. (5). The Toast-MATLAB and PYTHON bindings provide a convenient way to implement the inverse solver, although it is also possible to write a solver that directly links to the C++ interface. A simple example for a nonlinear Polak-Ribiére conjugate gradient scheme with line search follows:

As a further example, the following MATLAB code fragment implements a damped GN scheme with line search, given by the update rule

Where toastJacobian returns the Jacobian matrix $\mathsf{J}$ for the DOT problem for a given set of absorption and diffusion parameter coefficients, toastKrylov solves $({\mathsf{J}}^T\mathsf{J}+\tau {\mathcal{R}}^{\prime \prime })\mathbf{x}=\mathbf{g}$, toastLineSearch implements an inexact 1-D line search for step length $\alpha$, and hreg is a regularization object that supports the evaluation of value, gradient ${\mathcal{R}}^{\prime }$ and Hessian ${\mathcal{R}}^{\prime \prime }$ of the regularization term $\mathcal{R}$. In this example, the parameter vector $\mathbf{x}$ is transformed to logarithmic values, ensuring a restriction to positive values.

Where the computation of the Jacobian matrix is impractical due to its size, the solver can be modified to use a matrix-free approach, where the Frechet derivative $\mathsf{J}\mathbf{x}$ and adjoint Frechet derivative ${\mathsf{J}}^T\mathbf{y}$ are evaluated implicitly by function calls rather than using an explicit matrix representation. This allows the application of the GN solver to large-scale problems with a combination of high-grid resolutions and a large number of measurements.

A graphical user interface (GUI) MATLAB application for the DOT inverse problem included in the Toast++ package is shown in Fig. 14.

Fig. 14

GUI sample application for the DOT inverse solver included in the Toast-MATLAB toolbox.

An example of a reconstruction of the spatial distribution of the two coefficients in a cylindrical domain is shown in Fig. 15. In this example, frequency-domain forward data were generated on a mesh with 83,142 nodes and 444,278 tetrahedral elements, for all combinations of 80 sources and 80 detectors, arranged in five rings around the cylinder mantle. The forward mesh geometry, as well as source (red) and detector positions (blue) are shown in Fig. 15(a).

Fig. 15

Reconstruction of absorption and scattering features in a cylindrical object. (a) Left: target with absorption (red) and scattering features (blue). Source and detector positions on the surface are marked with white and cyan points, respectively. Middle: isosurfaces of features reconstructed with a generic TV regularizer. Right: reconstruction results with structural MRF prior. (b) Cross sections of absorption distributions. Rows from top to bottom: target, generic TV reconstruction, structural MRF prior reconstruction. (c) Cross sections of scattering distributions. Rows from top to bottom: target, generic TV reconstruction, and structural MRF prior reconstruction.

The data were then contaminated with 1.0% additive random Gaussian noise and used as an input for a GN reconstruction with TV regularization into a regular $48\times 48\times 48$ grid of trilinear basis functions. The mesh used by the inverse solver was coarser with 27,084 nodes and 141,702 tetrahedra elements. Figures 15(b) and 15(c) show cross sections of the target and reconstructed parameter distributions after 10 GN iterations. Figure 15(a) shows isosurfaces of the target and recovered inclusions at value 0.012 for ${\mu }_a$ and 1.2 for ${\mu }_s$.

5.3.

Regularization and Structural Priors

The Toast library contains regularization classes for augmenting the objective function of the inverse problem with a regularization term. Supported methods include Tikhonov, TV, and Markov random field (MRF) regularizers. In addition, structural information, such as internal boundaries with step changes in parameters can be incorporated by spatially varying regularization. Boundary information can be provided as a raster image, which can for example be obtained from a segmented magnetic resonance imaging or computed tomography image.

Using the MATLAB bindings, a TV regularization object can be created with

Where ‘TV’ denotes the regularization method (here a TV scheme), hbasis is the handle of a basis mapper object, x is the initial parameter vector, and tau is the regularization hyperparameter. In addition to these fundamental parameters required by all regularization constructors, individual regularization schemes may require addtional arguments. In this case, the TV object is supplied with a threshold parameter beta. After construction, the regularization object can be queried to return a value, gradient, second derivative, or its diagonal, given a parameter vector $\mathbf{x}$:

In addition to the TV prior, Toast supports a range of other regularization strategies, including Tikhonov, Huber, or Perona-Malik methods. With each regularization method, Toast also allows the inclusion of structural prior information about internal boundaries. This information is provided to the regularization constructor as an image of a diffusivity field. This allows application of the regularization term in a spatially varying manner. For the definition of a TV prior, the constructor would then be specified as

The right image of Fig. 15(a) and the bottom rows of Figs. 15(b) and 15(c) show the results of a GN reconstruction of the absorption and scattering distributions in a cylinder, using a MRF regularization scheme that is based on the neighborhood graph of the reconstruction basis. Structural prior information was included, where the shapes and locations of the internal features were assumed known. Consequently, the results are significantly better than the generic TV reconstruction shown in the same image. All objects are recovered with correct shape and locations. Recovery of feature contrasts is also very good, except, in the case with no prior, for the most challenging features consisting of the small high-contrast absorption feature, and the absorption feature embedded in a highly scattering sphere.