Home

Introduction
  About ODT
  About TOAST
  Matlab toolbox

Download
  Toast sources
  Matlab toolbox
  Installation

Documentation
  Getting Started
  Matlab demos
  FAQ
  Change log
  References
  License

Support
  Links
  Message board
  Contacts
  Vacancies

Toast toolbox tutorial: Building a DOT forward solver

This example shows how to write a simple DOT steady-state forward solver in Matlab using the Toast toolbox. To run this code yourself, you need Matlab, and you have to install the Toast toolbox on your computer.

To save you typing, the full Matlab script for this example can be downloaded here.

Step 1: Create a mesh

The forward solver uses the finite element method to solve the diffusion equation in a domain. An unstructured mesh is required to define the domain and distribution of parameters. The Toast toolbox contains a few helper functions to create simple meshes. More complex meshes can be generated with external mesh generators and loaded into Toast.

For this example, we create a simple two-dimensional circular mesh of radius 25 mm, using the following commands:

where nsect, nring, and nbnd define the element structure of the mesh.

You can retrieve the geometry parameters from an existing mesh:

where 'vtx' is the list of node coordinates, 'idx' is the list of node indices for each element, 'eltp' is the list of element type identifiers for each element, 'nnode' is the number of nodes, and 'nel' is the number of elements in the mesh.

To visualise the new mesh, use the 'Display' method of the toastMesh class:

which will open the following figure:

You can change the mesh resolution by varying the parameters used for mkcircle. Here is a mesh created with nsect=8, nring=16 and nbnd=4:

Step 2: The optical coefficients

The forward solver requires the distribution of absorption and scattering coefficient, as well as the refractive index.

The mesh defines a basis u(r) with which functions f(r) defined in the domain can be approximated with a finite dimensional array of coefficients:

In this example, we assume a piecewise linear basis, where coefficients f_i are defined on the vertices of each mesh element. (Toast also allows to define a piecewise constant basis, as well has higher-order polynomial bases. We'll come to those later.)

Toast constructs the basis u_i(r) automatically. All we have to provide are the index arrays for absorption μ_a, scattering μ_s and refractive index n:

mua, mus and ref now homogeneous parameter distributions. We will look at creating inhomogeneous distributions later.

Source and detector locations

We need to define at least one source distribution and one detector profile for the forward solver. For practical applications, multiple source and detector locations are usually employed.

For this example, we assumed that sources and detectors are arranged with equidistant spacing around the circumference of the mesh:

The SetQM method attaches the source and detector locations to the mesh. You can display the source and detector location superimposed on the mesh display from step 1:

which will show the source locations as red circles, and the detector locations as blue squares:

From the source and detector locations, you create the source and boundary projection vectors using the Qvec and Mvec functions:

where Neumann defines the source as an incoming flux (you could instead using Isotropic for internal source distributions). Gaussian defines the boundary profiles for the source and detector response functions, and '2' is their width.

Running the forward solver

We now have everything in place to set up and run the forward solver to simulate the measuremenets for our problem. The FEM formulation leads to a linear system of the form

where K is a system matrix depending on the parameters x (absorption, scattering, refractive index), Q is a matrix of column vectors, where each column consists of one source distribution, M is a matrix of column vectors for the detector response distributions, Φ is the matrix of photon density distributions for each source, and Y contains the measurements for each source and detector combination.

With the Toast toolbox, this system is represented by the following code

where dotSysmat builds DOT system matrix K from the mesh geometry and parameter coefficients (the '0' at the end refers to the modulation frequency, which is zero for a steady-state problem). The linear system is solved for Φ with the backslash operator. For larger problems, you may have to switch to an iterative scheme such as bicgstab or gmres. The measurements are obtained from the photon density distributions by projecting onto the measurement profiles M.

You can display the measurements as a sinogram:

or as a boundary profile as a function of source detector separation:

Note that all 16 measurement profiles coincide, because of the symmetry of the problem. For inhomogeneous problems or more complex shapes, this will no longer be the case.

We can now write out the generated data as a linear vector with source-major ordering for later re-use: