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: A simple parameter reconstruction example

This example shows how to write a simple DOT inverse solver for reconstructing optical parameters from boundary data 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 1: Defining a few helper functions

Project a photon density distribution to boundary measurements. This function takes an array of column vectors of complex nodal photon densities phi and maps them to boundary measurements, using measurement operators mvec. The measurements are returned as log amplitude and phase data:

Compute the objective function for an array of boundary data proj computed from a set of trial parameters, given measurement data data and scaling parameters sd:

Compute the objective function, given vector of logarithmic parameter solutions logx. This function is called by the line search step in the iterative loop:

For assessing the reconstruction results, evaluate the image error as a norm between target and reconstructed parameter distributions:

Step 2: Generating target data

For this example, we generate the target data directly with the forward solver. We use a high-resolution mesh for this, while the following reconstruction will be performed on a lower-resolution mesh to avoid an inverse crime. In addition, the simulated target data can be contaminated with noise. For the target data we use an inhomogeneous parameter distribution. The details of the python forward solver can be found here.

Step 3: Set up inverse problem

We first load the FEM mesh to be used during the reconstruction, and construct the source and measurement operators:

The initial parameter estimates are set to homogeneous distributions:

We also create a basis mapper object that allows to map parameter and field distributions between the mesh basis of the forward solver, and the regular pixel basis of the inverse problem:

The forward operator is applied to obtain boundary data from these parameter estimates:

We use the differences between target and projection data to define a data scaling vector:

We now map the parameters into the basis of the inverse solver:

and construct the solution vector:

Calculate the initial objective value:

Step 4: The iterative solver loop

We now construct the iterative nonlinear conjugate gradient (NCG) solver to find parameter distributions that minimise the objective function. This has the form:

where itrmax is an iteration limit, itr is the iteration counter, and step is the initial step length for the line search.

The body of the NCG contains the following components:

1. Construction of the Jacobian matrix (later we will discuss methods to avoid the explicit storage of the Jacobian):

The gradient of the cost function is given by J^TΔy, where Δy is the data difference:

We implement a Polak-Ribiere conjugate gradient to obtain an update direction. The step length is obtained from a non-exact one-dimensional line search:

Direction and step size are used to update the estimate of the parameter vector:

We can now map x back to the nodal absorption and scattering parameter distributions:

We can also map the solutions to images for display and error analysis during the reconstruction: