1.

Introduction

To tackle the principal problem of optical tomographic imaging of human tissue, which has the dominance of scattering, it is essential to have an accurate model for light propagation. Grunbaum’s model was one of the earliest, wherein photon propagation is considered as a one- or two-step Markov process.¹ This comes under one of the two main types of forward models, namely stochastic models, under which come the other well-known models, such as Monte-Carlo² and random walk³ as well. The second main type is the analytic models, based on partial differential equations, under which comes the Boltzmann transport model⁴ and its approximation using the diffusion equation (DE).^{5 6} Diffusion approximation of the light transport has been investigated thoroughly,^{6 7 8} and it is found to be inaccurate to describe photon propagation near the boundary within distances that are comparable to transport mean free path, and in regions where the reduced scattering coefficient (μ_s ^′) is less than or comparable to the absorption coefficient (μ_a). ⁹

On the other hand, the Boltzmann transport equation holds good even in the previously described situations where the diffusion equation fails. The transport equation (TE) distinguishes between the scattering and absorbing inhomogeneities inside the medium through use of the directional information of photons. Both the forward and inversion methods based on the TE model of light propagation were tried by many researchers.^{4 10 11 12 13} Dorn⁴ developed recently a TE-based method for optical tomographic inversion, wherein the needed derivatives were calculated using the adjoint form of the TE. Optical property estimation, based on solving the forward problem and inversion using gradients calculated through adjoint differentiation, both based on TE, was attempted by Hielscher et al.;^{10 11 12 13} Both the time-independent and time-dependent cases were solved for isotropic and anisotropic scattering background and inhomogeneities. The results from both the synthetic and experimental data are very promising, which provide a means of tomographic imaging of objects where the DE fails to hold. One drawback that runs through all these methods is the relatively large computational time for inversion, which needs repeated solution of the forward problem, on account of the need to take care of the new set of variables in the TE, namely the angles of scattering.

Inversion strategies based on stochastic models of light transport have also been attempted with varying degrees of success in the past. Methods based on Monte-Carlo (MC) simulation,^{9 14} which also takes into account the directional information of photons, have yielded reconstruction results that should pave the way for developing robust algorithms for reconstructing absorption and scattering coefficients. Therefore, these methods should have larger applicability like the ones using the TE model. Here again we have the drawback of a large computational time involved to take many millions of photons through the tissue to complete an MC simulation. An inversion method that makes use of repeated application of the forward model and updating of the derivative needed for guiding the solution to convergence would become too expensive to be a useful option. This problem has been recently addressed⁹ by using the perturbation Monte-Carlo (PMC) method to extract Frechet derivative information, which is the derivative of the forward model with respect to the optical absorption (μ_a) and scattering coefficient (μ_s). (In a general case, the Frechet derivative is the derivative of the forward model used with respect to the optical properties of the medium. After discretizing the forward model and the properties in a finite dimensional space, the derivative can be represented by a Jacobian matrix, representing the rate of change of data with respect to the optical properties in the pixels.) In Ref. 9, an analytical expression for the perturbation in the detected photon weight consequent to changes in optical properties in a subregion, borrowed from the neutron-transport theory, was used to update both Jacobian and the computed forward data. After discretization, such a procedure will require only a single Monte-Carlo (MC) simulation with the derivative as well as forward projection data extraction handled by the PMC, which is rapidly done. In Ref. 9, this was used to solve a simple two region inverse problem of photon migration in a heterogeneous slab of thickness comparable to the transport mean free path, an object that falls under the transport regime, where DE fails to hold.

In this work, using the basic MC simulation as a forward model, we extend the PMC approach to construct a Jacobian matrix for use in a perturbation-based method to reconstruct transport-regime low-scattering objects with discretized absorption and scattering coefficient distributions. The object was discretized to 67×67 pixels (or grid points), and for a single view with one source and 11 detectors, we construct a Jacobian matrix of dimension either 11×4489 or 11×8978, depending on whether we reconstruct one of the optical properties or both. Similarly, for 11 such views, to reconstruct from the full dataset, we construct Jacobian matrices of dimension 121×4489 or 121×8978. The details of arriving at the Jacobian matrix are given in Sec. 3. For this reconstruction problem, we did not tailor the perturbation equation [Eq. (1) of Ref. 9] to extract changes in the computed forward projection data to update it without recourse to a full MC. The forward projection data modified this way were found to be inaccurate (compared with data obtained by running a full MC simulation), and application of such computed projection data in the iterative reconstruction algorithm led to solutions that were erroneous. Therefore, we chose to repeat the MC with a relatively smaller number of photons in each of the iterations (compared to the original MC done on the homogeneous phantom to extract the information needed to construct and update the Jacobian). We show that the inversion algorithm using the Jacobian matrix obtained from this MC-based method is able to give superior reconstructions for low-scattering objects (considering both the value of the inhomogeneity, its location, and shape) compared to methods that use the diffusion equation as the propagation model.

The work presented here is as follows. In Sec. 2, the basic iterative reconstruction scheme is introduced that uses the Jacobian to calculate the update vectors. The Jacobian matrix is calculated using the PMC in Sec. 3. Then, in Sec. 4, the details of the numerical experiments conducted are described. The test objects and the geometry of data collection are also given in Sec. 4. We consider objects with both absorption and scattering inhomogeneities, and two slightly different strategies for reconstruction. 1. The propagation-backpropagation (PBP) method¹⁵ is considered when data from a single source (and many detectors) are reconstructed to arrive at an updated object, which is used as the initial estimate for reconstructing data from the next source. 2. Data from all the sources are assembled and solved as one problem. Since the number of equations in the first case is smaller, the reconstruction problem is too ill posed to converge to a solution when both scattering and absorption inhomogeneities are considered simultaneously without a priori knowledge on location. When they are considered separately, the PBP method converges to a solution faster than when all the data are considered together. The convergence to a solution was obtained for both inhomogeneities simultaneously in the second case, when the data from all the sources are considered in a single problem. (Details are given in Sec. 4.) Further discussion of the results and conclusions are given in Sec. 5.

2.

Iterative Reconstruction Algorithm

Given that the forward model F is acceptable as accurate, the methods in optical tomographic reconstruction seek for a solution set (μ_a,μ_s) that minimizes an objective functional ϕ(μ_a,μ_s). The objective functional is usually the mean-squared norm of the difference between experimental data W^e and W^c, the computed data obtained by applying the forward model F on the current estimate (μ_a ⁱ,μ_s ⁱ). In this work, data (W) used in the reconstruction of μ_a and/or μ_s is integrated intensity. A Newton type of algorithm to solve this minimization problem leads to an iterative step,¹⁶ given by

The steps involved in the iterative reconstruction procedure are shown in the flow chart of Fig. (1). It has two iterations, one main part (the outer iteration) and another subsidiary part (the inner iteration).¹⁶ The i’th solution vector (μ_a ⁱ,μ_s ⁱ) given to the forward model (the Monte-Carlo procedure is repeated with 1 million photons to model the forward propagation) predicts the computed data (W^c). The experimental data (W^e) when plugged in helps us find ΔW(=W^e−W^c), which is used to set up the update equation [Eq. (1)]. The inner iteration outputs the update vector (Δμ_a,Δμ_s), which is used to arrive at the new solution μ_a ⁱ⁺¹=μ_a ⁱ+Δμ_a and μ_s ⁱ⁺¹=μ_s ⁱ+Δμ_s. In the next section, we elaborate on the construction of the Jacobian matrix used in Eq. (1).

Figure 1

The flowchart used in the iterative reconstruction. The inner loop uses a gradient search algorithm to output update vectors for the optical properties.

3.

Jacobian Calculation Through Perturbation Monte-Carlo Method

Light transport through tissue can be studied through stochastically following the absorption and scattering of photons in turbid media. A large number of photons are launched into the tissue and traced through the probabilistically defined absorption and scattering until they either exit the tissue or are absorbed. A large number, usually 6 to 10 million, of photons are required to extract statistically meaningful average quantities such as photon arrival times and arrival histograms, pathlength distribution, and a temporal point spread function (TPSF). The distance traveled between two successive scattering and/or absorbing events s, the angle of scattering, one azimuthal (ψ) and one polar (θ) are assigned values following three uniformly distributed random numbers x₁, x₂, and x₃ (for example). The following equations are used:

The basic MC simulation has the drawback of a very large computation time to arrive at time histograms of photons. To overcome this limitation, a hybrid Monte-Carlo procedure was developed that makes use of only one full MC simulation.¹⁴ This was done for a homogeneous background with μ_a ^b and μ_s ^b as the absorption and scattering coefficients, respectively, and the coordinates of the interaction points of all the detected photons are stored. Thereafter, perturbation in μ_a ^b and μ_s ^b are introduced in specified locations, and the new photon weight due to these perturbations is evaluated using the already stored coordinates and path lengths, and two scaling laws connecting the modified photon weights to optical property changes in specified regions.¹⁴ Analytic expressions for changes introduced in detected photon weights owing to optical property perturbation in a subregion in the object are used to calculate both the rate of change of photon weights with respect to optical properties and the new photon weights.^{9 20 21} Specifically, if μ_a and μ_s are perturbed in certain locations and become μ¯_a=μ_a+δμ_a and μ¯_s=μ_s+δμ_s, then the detected photon weight w changes to w¯ as

We have also tried to use the PMC to compute the perturbed photon weight (W^c) owing to changes in optical properties (brought about by updates during the iterations) using the scaling laws of Eq. (5). The procedure used is given next. Subregions of pixels along typical detected photon trajectories (selected at random) are identified, and the corresponding n and l values are used in Eq. (5) to find the modified photon weight at the detector caused by changes in μ_a and μ_s values in all the pixels along the trajectories. Since the subregions of adjacent pixels overlap, and also because the area covered by the subregion is more than the area covered by a pixel, the change in photon weight calculated W is scaled down by a factor equal to the ratio between the area of a pixel to the subregion. The modified W is used to compute changes in photon weights at other detectors by scaling this W proportional to the previous photon weights at these detectors. The updated weights at the detectors had an overall error, when compared to the modified weights obtained by running a full MC simulation with a smaller number of photons, equal to the detected photons. This error is unacceptably large, for the iterative reconstruction algorithm did not converge, in this case, to optical property distributions close enough to the original distributions.

4.

Numerical Simulations

We have considered a 2-D cross section of a cylindrical tissue-mimicking phantom of diameter 66 mm. The center of the cylinder is chosen to be at (34,34). The background absorption and scattering coefficients are represented by μ_a ^b and μ_s ^b (in mm⁻¹), respectively. The circular inclusions of diameter 12 mm with optical properties denoted by μ_a ⁱⁿ and μ_s ⁱⁿ (in mm⁻¹) took different values for the different objects considered. The background as well as the inclusions has low-scattering coefficients, which make them transport-regime objects where the DE fails. The object was divided into 67×67 pixels of uniform size and the optical properties were considered constant in each of these pixels (grid elements). The anisotropy factor g is kept at 0.90 and the refractive index of the tissue at 1.33, for all the cases discussed here. The number of detectors we have considered for each of the sources are 11, and the detector diameter is 4 mm. The source-detector assembly is shown in Fig. (2). For collecting data for different views, this source-detector assembly is rotated in unison by 32.73 deg, such that with 11 source positions, the entire object is covered all around. The data for reconstruction (W^e) is obtained by running MC simulations through the inhomogeneous object distribution we want to reconstruct. One million photons are launched from the source, and the arrival times at the parallel detectors are found from the distances traversed by the photons. The temporal point spread functions (TPSFs) are obtained in this way, which are integrated to give us the measurement, the integrated intensity. This procedure is repeated for all the source locations. To the overall set of measurements, a 2 Gaussian noise was added, which gave us the experimental projection dataset (W^e). In all the numerical simulations, we have considered two types of objects: 1. objects with only one inhomogeneity either in μ_a or μ_s, and 2. objects with two inhomogeneities, either two absorption changes at two spatial locations or an absorption and a scattering inhomogeneity at different locations. For all these cases, we used the homogeneous background optical properties (μ_a ^b and μ_s ^b) as the initial estimate to start the iteration.

Figure 2

The data gathering geometry: the source S and detectors D1 through D11 are rotated in unison for data gathering at different views.

We first run an initial Monte Carlo with 9.6 million photons with extended boundary and homogeneous background values of μ_a ^b and μ_s ^b and store n and l, the average number of collisions and length of traverse of photons, in the 12-pixel subregions surrounding each of the pixels. This dataset is needed for calculation of the Jacobian, as given in Sec. 3 (see Fig. 1). For forward projection data, 1 million photons are taken across the object. For objects that have either only one type of inhomogeneity (either μ_a or μ_s) without a priori information on location, or that have both inhomogeneities at a priori known locations, the PBP strategy gives good reconstructions. Here data from one view is considered sufficient to solve the iteration that resulted in an updated object for use with data from the next view. All the views are considered cyclically. With μ_a - and μ_s inhomogeneities simultaneously present, without assumption of locations, the PBP-based problem was too ill posed to result in a solution. For dealing with such an object, we combined the data from all the views, which are solved as one problem. For such objects, the iteration with full data converged is discussed later in this section. The Jacobian matrix needed in the iterations, for the PBP case or for the full-data case, is computed using the method in Sec. 3. The convergence criterion used was that the norm of the difference between computed data and the experimental data should be below a small preset value. Since we used Monte-Carlo simulation, a stochastic forward model, the reconstructed images required postprocessing. We used a local median filter of dimension 5×5 pixels to smooth the images.

To show the effectiveness of this algorithm, we have considered six different objects, which are discussed separately in the following sections. All the simulations are carried out on a Pentium IV 2.40 GHz processor.

4.1.

Case 1

We assume that a priori information about the spatial location of the inhomogeneity is available. The unknowns are the value of μ_a and μ_s at these spatial locations. What is assumed a priori is that the inhomogeneity is somewhere inside a region of 30×30 pixels centered at (52,34). Here, we have used the PBP approach for reconstruction with 900×2 unknowns, and the Jacobian for this problem has a dimension of 11×1800. The original object shown in Fig. 3(a) has a μ_a inhomogeneity. The reconstructed μ_a distribution is shown in Fig. 3(b), and horizontal cross sections at y=34 mm of Figs. 3(a) and 3(b) are shown in Fig. 3(c). The solution converged in eleven iterations when we finished a complete dataset from all the sources. Each iteration took 17 min. The original object μ_a -inhomogeneity value is 0.14 mm⁻¹. The scattering coefficient in the homogeneous region was μ_s ^b, and the background value was 0.50 mm⁻¹. Reconstructing μ_s is also kept as an unknown in the known locations of inhomogeneity, and the reconstructed μ_a ⁱⁿ and μ_s ⁱⁿ at the center are 0.143 and 0.483 mm⁻¹, respectively.

Figure 3

Reconstruction obtained from the PBP approach with Jacobian calculated through PMC with a priori information on the spatial location of inhomogeneity. (a) Original object distribution: background μ_a ^b and μ_s ^b are 0.04 and 0.5 mm⁻¹, respectively, and the inclusion has μ_a ⁱⁿ=0.14 mm ⁻¹ and μ_s ⁱⁿ=0.5 mm ⁻¹. (b) Reconstruction of (a). The reconstructed inhomogeneity optical properties at its center are μ_a ⁱⁿ=0.143 mm ⁻¹ and μ_s ⁱⁿ=0.483 mm ⁻¹. (c) Horizontal cross sections at y=34 mm through the center of images (a) and (b). The solid line shows the original, and the dashed line shows the reconstructed object profile.

4.2.

Case 2

Here we have both μ_a - and μ_s inhomogeneities located at different a priori known locations (as before, the locations are specified to be within 30×30 pixel regions). Here μ_a ^b is assumed to be 0.04 mm⁻¹ and μ_s ^b=0.30 mm ⁻¹. We have used the PBP approach for reconstruction. There are 900×2 unknowns, and the Jacobian size is 11×1800. The original μ_a ⁱⁿ(=0.14 mm ⁻¹) and μ_s ⁱⁿ(=0.70 mm ⁻¹) images are shown in Figs. 4(a) and 4(b), respectively. The reconstructions are shown in Figs. 4(c) and 4(d), with the central value of μ_a ⁱⁿ - and μ_s ⁱⁿ reconstructions 0.132 and 0.72 mm⁻¹, respectively. The solution converged in 11 iterations, each iteration taking 17 min. Figures 4(e) and 4(f) are horizontal cross sections at y=44 mm and y=24 mm through the reconstruction shown in Figs. 4(c) and 4(d), respectively.

Figure 4

Simultaneous reconstruction of μ_a and μ_s inclusions using the PBP approach with a priori information on the spatial location of inhomogeneities. (a) Original object showing μ_a inhomogeneity with background μ_a ^b=0.04 mm ⁻¹, μ_s ^b=0.30 mm ⁻¹, and inclusion has μ_a ⁱⁿ=0.14 mm ⁻¹. (b) Original object showing μ_s distribution: background is same as (a) and the inclusion μ_s ⁱⁿ=0.70 mm ⁻¹. Reconstruction of (c) μ_a inhomogeneity and (d) μ_s inhomogeneity. The reconstructed optical properties at the centers of the inhomogeneities are μ_a ⁱⁿ=0.132 mm ⁻¹ and μ_s ⁱⁿ=0.72 mm ⁻¹. (e) Horizontal cross sections at y=44 mm through the center of images (a) and (c). The solid line shows the original, and the dashed line shows the reconstructed object profile. (f) Horizontal cross sections at y=24 mm through the center of images (b) and (d). The solid line shows the original, and the dashed line shows the reconstructed object profile.

4.3.

Case 3

This is a case where we do not assume any a priori information about the spatial location of the inhomogeneity. We assume that the inhomogeneity is in μ_a only. Using the PBP strategy with 67×67 unknown μ_a values, the Jacobian dimension is 11×4489. For forward calculations, μ_a ^b and μ_s ^b are kept at 0.04 and 0.50 mm⁻¹, respectively. The original μ_a distribution is shown in Fig. 5(a), and the reconstructed μ_a distribution is in Fig. 5(b). The solution converged in 15 iterations, each iteration taking 20 min. The original object μ_a inhomogeneity value is 0.14 mm⁻¹ and the reconstructed object μ_a inhomogeneity value at its center is 0.148 mm⁻¹. Cross sectional plots at y=34 mm across the original and the reconstruction are shown in Fig. 5(c).

Figure 5

Reconstruction without a priori knowledge of inhomogeneity location using the PBP approach. (a) Original object distribution: background μ_a ^b and μ_s ^b are 0.04 and 0.5 mm⁻¹, respectively, and the inclusion is an absorbing inhomogeneity with μ_a ⁱⁿ=0.14 mm ⁻¹. (b) Reconstruction of (a). The reconstructed inhomogeneity optical property at its center is μ_a ⁱⁿ=0.148 mm ⁻¹. (c) Horizontal cross sections at y=34 mm through the center of images (a) and (b). The solid line shows the original, and the dashed line the reconstructed object profile.

4.4.

Case 4

Here we consider another case where the inhomogeneity location is not known. We have considered both the PBP approach using data from one view at a time, and then reconstruction using the full dataset. The μ_a ^b and μ_s ^b for this case are 0.04 and 0.50 mm⁻¹, respectively. The μ_a inhomogeneity is as shown in Fig. 6(a). The reconstructions are as shown in Fig. 6(b) (using the full dataset) and Fig. 6(c) (using the PBP approach). The values of μ_a at the center of reconstructed inhomogeneity in the prior two cases (0.153 and 0.156 mm⁻¹) are higher compared to the original object inhomogeneity (0.14 mm⁻¹), but are not quite different from one another. However, the time taken by the approach using one view at a time (3.67 h for 11 iterations) is much smaller-compared to that for the full dataset (15.6 h for four iterations). Cross sectional plots at y=24 mm across the reconstructions in Figs. 6(b) and 6(c) are shown in Fig. 6(d).

Figure 6

Comparison of the PBP approach with the procedure using full data. In both the approaches, location information was not assumed. (a) Original object distribution: background μ_a ^b and μ_s ^b are 0.04 and 0.5 mm⁻¹, respectively, and the inclusion was an absorption inhomogeneity of μ_a ⁱⁿ=0.14 mm ⁻¹. (b) Reconstruction of (a) from full data considered together. The reconstructed inhomogeneity optical properties at its center, μ_a ⁱⁿ=0.153 mm ⁻¹. (c) Reconstruction of (a) with the PBP approach. The reconstructed inhomogeneity optical properties at its center, μ_a ⁱⁿ=0.156 mm ⁻¹. (d) Horizontal cross sections at y=24 mm through the center of images (a), (b), and (c). The solid line shows the original, the dotted line shows the reconstructed object profile as shown in (b), and the dashed line shows the cross section across (c).

Figure 7 shows the results for a similar effort done on data from an object with only μ_s inhomogeneity with location unknown. Figures 7(b) and 7(c) give the results from the full dataset and the PBP-based algorithm, respectively. The μ_a ^b and μ_s ^b are kept at 0.04 and 0.30 mm⁻¹, respectively. An original μ_s inhomogeneity of 0.70 mm⁻¹ [Fig. 7(a)] was reconstructed as 0.74 and 0.72 mm⁻¹ at the center of the inhomogeneity for the two cases [Figs. 7(b) and 7(c)]. We repeated the PBP-based algorithm after re-estimating the set of n and l values at the end of each iteration by quick MC simulations using 0.57 million photons to update the Jacobian. The computational times involved are for Figs. 7(b), 7(c), and 7(d): 21 h (converged in five iterations), 6.97 h (converged in 19 iterations), and 7.5 h (converged in 15 iterations), respectively. The value of μ_s inhomogeneity in Fig. 7(d) at its center is 0.73 mm⁻¹. Cross sectional plots at y=34 mm across the reconstructed images in Figs. 7(b), 7(c), and 7(d) are shown in Fig. 7(e).

Figure 7

Same as the result of Fig. 6. Comparison of the PBP approach with the procedure using full data with no a priori assumption on location of inhomogeneity. The result is for an object with only μ_s inhomogeneity. (a) Original object with μ_s ⁱⁿ=0.70 mm ⁻¹ and background μ_a ^b and μ_s ^b are 0.04 and 0.30 mm⁻¹, respectively. (b) Reconstruction of (a) from full data considered together. The reconstructed inhomogeneity optical property at its center is μ_s ⁱⁿ=0.74 mm ⁻¹. (c) Reconstruction of (a) with the PBP approach. The reconstructed inhomogeneity optical property at its center is μ_s ⁱⁿ=0.72 mm ⁻¹. (d) Reconstruction of (a) with the PBP approach with re-estimation of n and l after each view. The reconstructed inhomogeneity optical property at its center is μ_s ⁱⁿ=0.73 mm ⁻¹. (e) Horizontal cross sections at y=34 mm through the center of images (a), (b), (c), and (d). The solid line shows the original, the dotted line shows the reconstructed object profile as shown in (b), the dashed line is for that shown in (c), and the dash-dot line is for that shown in (d).

4.5.

Case 5

Here we consider the example of an object with both μ_a and μ_s inhomogeneities, where the location of these inhomogeneities are left as unknowns. The number of unknowns is 2×4489, and the PBP-based algorithm is too ill posed to converge and is abandoned. The full-data problem, however, did converge in 13 iterations, each iteration taking 6.7 h. The background values μ_a ^b and μ_s ^b in this case are 0.04 and 0.30 mm⁻¹, respectively. The results are shown in Fig. 8(c), reconstructing μ_a inhomogeneity at the center as 0.138 mm⁻¹ (original value 0.14 mm⁻¹), and Fig. 8(d) reconstructing μ_s inhomogeneity at the center as 0.58 mm⁻¹ (original value 0.60 mm⁻¹). Figure 8(e) gives the cross sectional plots at y=34 mm through images in Figs. 8(a) and 8(c), and Fig. 8(f) gives the same at y=34 mm across the images in Figs. 8(b) and 8(d).

Figure 8

Reconstruction of an object with both μ_a and μ_s inhomogeneities with no a priori knowledge on location of the inhomogeneities using the data from all views simultaneously. Original object background μ_a ^b=0.04 mm ⁻¹ and μ_s ^b=0.30 mm ⁻¹ with inclusion: (a) μ_a ⁱⁿ=0.14 mm ⁻¹ and (b) μ_s ⁱⁿ=0.60 mm ⁻¹. (c) Reconstruction of (a). The reconstructed inhomogeneity optical properties at its center are μ_a ⁱⁿ=0.138 mm ⁻¹. (d) Reconstruction of (b). The reconstructed inhomogeneity optical property at its center is μ_s ⁱⁿ=0.58 mm ⁻¹. (e) Horizontal cross sections at y=34 mm through the center of images (a) and (c). The solid line shows the original, and the dashed line the reconstructed object profile. (f) Horizontal cross sections at y=34 mm through the center of images (b) and (d). The solid line shows the original, and the dashed line the reconstructed object profile.

4.6.

Case 6

Here the aim is to demonstrate the capability of these reconstruction algorithms to handle multiple inhomogeneities of the same type (either μ_a or μ_s) . The test object has two μ_a inhomogeneities of 0.45 mm⁻¹ in the left and 0.30 mm⁻¹ in the right, with μ_a ^b and μ_s ^b as 0.14 and 0.70 mm⁻¹, respectively [Fig. 9(a)]. The location is left as unknown. The number of unknowns to be reconstructed are 4489 μ_a values. We used the PBP-based algorithm, which converged in 23 iterations, each iteration taking 20 min. The reconstructed object is shown in Fig. 9(b) with the value of inhomogeneities at the center being 0.41 and 0.29 mm⁻¹, respectively. Cross sectional plots at y=34 mm across Figs. 9(a) and 9(b) are shown in Fig. 9(c).

Figure 9

Reconstruction of an object with two μ_a inhomogeneities done with the PBP approach. The location of the inhomogeneities is not assumed. (a) Original object distribution: background μ_a ^b and μ_s ^b are 0.14 and 0.5 mm⁻¹, respectively, and the inclusions have absorption coefficients of 0.45 and 0.30 mm⁻¹. (b) Reconstruction of (a). The reconstructed inclusions have values μ_a ⁱⁿ=0.41 mm ⁻¹ and μ_a ⁱⁿ=0.29 mm ⁻¹ at the center. (c) Horizontal cross sections at y=34 through the center of images (a) and (b). The solid line shows the original, and the dashed line the reconstructed object profile.

Finally, to show that a diffusion equation-based algorithm does not converge to meaningful results for these objects that fall under the transport regime, we have reconstructed the data for simulation used in case 3 [reconstruction using the MC-based method is shown in Fig. 5(b)] using the time-received optical absorption and scattering tomography (TOAST) algorithm.²² The reconstruction obtained is shown in Fig. 10, which is unable to distinguish the inhomogeneity.

Figure 10

Reconstruction of the object of Fig. 5(a) using a DE-based reconstruction algorithm, which failed to locate the inhomogeneity.

5.

Discussion and Conclusions

We confirm by numerical simulation the usefulness of a Monte-Carlo-based method to reconstruct low-scattering objects, where the diffusion equation fails to hold. The method is tested with objects with up to two inhomogeneities, either in μ_a and/or μ_s, and the reconstructions obtained are reasonable. Data are handled either one view at a time or in the usual way with all views put together in one problem. The one view at a time method has the advantage of reduced computation time, but cannot handle, because of ill posedness, two inhomogeneities (μ_a and μ_s) in unknown locations, for in this case, all the μ_a and μ_s values at all pixels should be considered as unknowns. Reconstruction in the previously mentioned case is obtained by considering all the views simultaneously.

When handling the PBP-based reconstruction, the number of iterations for convergence is dependent on which view we started the iterations from. For example, if we started with a view for which the sum of integrated intensities in all the detectors is a minimum, the convergence is the fastest. The variation in the sum of the integrated intensities as the source-detector combination is moved around the object because of the inhomogeneity location in the object. If, for a view, either the source or some of the diametrically opposite detectors is close to the inhomogeneity (or one of the inhomogeneities), the sum of the integrated intensities becomes a minimum. Starting off with such a view will be favorable to the PBP algorithm for convergence. For example, for the object of case 3 in Sec. 4.3, the algorithm converged in 15 iterations when we started with this favorable view. On the other hand, when we started from a view at right angles to the earlier mentioned one, the convergence was obtained in 19 iterations. An updated optical property vector, which resulted from solving the data from a particular view, always reduced the error norm |ΔW| if we started with the correct view. If we did not, |ΔW| did not always show a decrease, and at times the decrease was very sluggish. For an object with a centrally, symmetrically placed inhomogeneity, the error |ΔW| always decreased irrespective of which view we started our iterations with.

In the present work, we use only one of the two advantages offered by the PMC approach in the iterative reconstruction procedure. That is constructing and then updating the Jacobian right through the iterations. It should be possible to use the other advantage, namely getting the updated integrated intensity without running an MC procedure each time. But our experience so far has been that the perturbation approach modified to this reconstruction problem results in photon-weight updates that are not accurate enough (as evidenced by comparison with updates obtained through direct MC simulation) for use in the iterative reconstruction procedure. In this work, we use the MC simulation for updating the photon weights with a reduced number of photons. This is similar to the hybrid Monte-Carlo procedure without the need for storing the detected photon paths. Meaningful average quantities, such as expected photon weights in our case, can still be ascertained with reasonable accuracy, because we are dealing with objects with highly forward scattering anisotropy and low-scattering coefficients.