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 python module tutorial: Building a DOT forward solver

This example shows how to write a simple DOT steady-state forward solver in python using the Toast module. To run this code yourself, you need python with the numpy, scipy and matplotlib modules, as well as the Toast python module.

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

Step 1: Preliminaries

We need to import a few standard modules:

Step 2: Load 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. For this example, we load a simple 2-D circular mesh stored in a file:

The function should respond with the message

The geometrical mesh data (node positions and element connectivity) can be extracted from the mesh handle with:

where nlist is an array of node coordinates, elist is the element connectivity list, and eltp is a list with element type identifiers. nlen is the number of nodes, used later in the script.

Step 3: 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 read a list of source and detector locations predefined in a file, and associate it with the mesh:

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

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, %Phi; is the matrix of photon density distributions for each source, and Y contains the measurements for each source and detector combination.

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

where mesh.Sysmat builds 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 sparse linear system is solved for Φ using the iterative BiCGSTAB solver provided with scipy. Other solvers such as GMRES are also available.

We now need to map the fields to boundary measurements by projecting onto the measurement profiles, and convert to logarithmic data:

You can display the measurements as a sinogram, using matplotlib plotting functions: