Toast toolbox tutorial: Forward solver for inhomogeneous parameters
This example builds on the homogeneous parameter example described in the previous example, by setting up an inhomogeneous distribution for the absorption and scattering parameters.
The Matlab script for this example can be downloaded here.
Step 1: Creating parameter maps
The simplest way to define an inhomogeneous parameter distribution for a 2D problem is by creating an image and mapping it onto the mesh basis. Below are two 128x128 bitmaps to be used as the absorption and scattering parameter distributions.
The images are stored in a common bitmap format (in this case, greyscale .png), and loaded into Matlab with the standard imread command:
bmua = imread('demo_matlab_fwd2_mua.png');
bmus = imread('demo_matlab_fwd2_mus.png');
The images are scaled to the desired parameter values:
bmua = double(bmua)./255.*0.02 + 0.01;
bmus = double(bmus)./255.*1.0 + 1.0;
axis equal tight; colorbar
axis equal tight; colorbar
Step 2: Mapping the images in the mesh basis
The images now have to be mapped into the mesh basis. To do this, we first have to create a mapping object:
grd = size(bmua);
basis = toastBasis(mesh,grd);
where 'mesh' is the mesh object for a circular mesh created as in the previous example, and grd are the dimensions of a regular grid basis (set here to coincide with the target image sizes).
toastBasis creates the mapper object for mapping between unstructured mesh and regular pixel basis.
We can now use the Map method of the basis object to map the parameter images into the mesh basis, and display them with mesh.Display:
mua = basis.Map('B->M',bmua);
mus = basis.Map('B->M',bmus);
Note that the edges of the inclusions are slightly jagged and washed out, because the outlines cannot be mapped exactly into the piecewise linear mesh basis. If required, a higher precision can be obtained with a more refined mesh or adaptive refinement along the edges.
The third parameter distribution, for refractive index n, is kept constant. So as in the previous example, we can write
nnd = mesh.NodeCount;
ref_bkg = 1.4;
ref = ones(nnd,1) * ref_bkg;
Step 3: Invoking the forward solver
The rest of the code for this example is identical to the previous one: create the source and measurement vectors, build the system matrix, and solve the linear system.
The sinogram and measurement profiles of the resulting boundary measurements look like this:
It can be seen now that there are (small) variations between the measurement profiles of individual source distributions, caused by the parameter inhomogeneities.
Step 4: Differences to homogeneous results
To see the effect more clearly, we can display the differences to the homogeneous results (this assumes that you have already run Demo 1):
data_homog = toastReadVector('demo_matlab_fwd1.dat');
logYhomog = reshape(data_homog,nq,nq)';
dlogY = log(Y)-logYhomog;
xlabel('source index q');
ylabel('detector index m');
axis equal tight;
ywrap = [dlogY(i:end,i); dlogY(1:i-1,i)];
axis([0 360 -1 0.1]);
xlabel('angular source-detector separation');
ylabel('log intensity perturbation');
It can be seen that the perturbation consist mainly of a reduction in the measured intensity, caused by the increased absorption inclusions. Only for a few measurements in the vicinity of the scattering object close to the boundary, a small increase in intensity due to backscattering is observed.
The task of the inverse solver discussed in the next tutorial will be to reconstruct the internal parameter distributions from these measurement perturbations.