The 3D object (
) used in the simulations is a cylinder 25 mm in diameter and 16 mm in height and represents a finite element mesh of an idealized section of a mouse abdomen where adipose, bowel, kidney, and bone tissues are identified and assigned realistic wavelength-dependent optical properties. It is assumed that different regions are distinguished through CT x-ray imaging or other imaging technique and all elements within each region have the same optical properties. For the BLT, 64 uniformly distributed detectors are used around the circumference and 15 layers of these 64 detectors separated by 1 mm are used in the

*z*-direction as shown in . The total number of detector readings at each wavelength is

.

shows three cross sections of the object that illustrate the positions of the bioluminescence sources to be reconstructed by the proposed algorithm. The circular transverse cross section is taken at the middle of the object (*z* = 0) parallel to the upper and lower surfaces of the object. is a sagittal cross section that passes through the bowel and bone while shows a coronal cross section passing through the two kidneys. shows 3 spherical bioluminescence sources of 5 mm and 3 mm diameter centered in the bone at (0, −7, 0) and two kidneys (−4.5, −2.5, 0) and (4.5, −2.5, 0), respectively.

In the first step, the optical properties of all tissues at the wavelengths 500, 550, 600, 650, and 700 nm were reconstructed using the DOT algorithm described in [

17]. 16 external sources uniformly distributed around the middle of object have been used for the illumination. Associated with each source, 15 layers of detectors were uniformly distributed across the object height with a separation between layers of 1 mm. Each layer contains 49 detectors forming a field of view of 270 degree. The light fluence rate at the “detectors” could be derived from the corresponding pixels of the image captured by a CCD camera. It is assumed that an accurate segmented 3D image set is provided by x-ray CT, MR or by a deformable atlas of mouse anatomy. All elements within each region in the object are assumed to have the same optical properties. The light fluence rate at the boundary is written as a Taylor series expansion around an initial guess corresponding to a homogenous object with the same optical properties in all regions. The fluence rate is approximated by the first three terms, including the first and second order derivatives which are calculated by the direct method and give the change of light fluence rate at the boundary due to small changes in the tissue optical properties in each region [

17]. The first and second order derivatives are used in an iterative algorithm to reconstruct the tissue optical properties. For this simulation, forward finite element calculations of the light fluence rate were performed using the true optical properties and Gaussian noise, which is an approximation to the Poisson measurement error distribution, of 2% of the signal was added to obtain the data for reconstruction,

*i.e*.

, where randn is a Matlab function that generates random numbers drawn from a normal distribution with mean zero and standard deviation one. In
, the true and reconstructed scattering coefficients,

and

respectively, are presented while
shows the true and reconstructed absorption coefficients,

and

respectively. The true scattering and absorption coefficients were calculated using the data reported in [

30,

31]. The maximum relative error obtained for scattering coefficients is 22% at 600 nm while the maximum relative error obtained for the absorption coefficients is 21% at 700 nm.

| **Table 2**True [30,31] and reconstructed scattering coefficients (mm^{−1})for each region in the heterogeneous object at five different wavelengths |

| **Table 3**True [30,31] and reconstructed absorption coefficients (mm^{−1}) for each region in the heterogeneous object at five different wavelengths |

The true optical properties at each wavelength and the true source distribution were used by the forward model to calculate the bioluminescence fluence rate, and 2% Gaussian noise was added to simulate realistic data:

in

Eq. (6). The Green’s function at each wavelength

was calculated using the reconstructed optical properties in and obtained by the DOT algorithm. The Green’s functions and light fluence rates were normalized as shown in

Eq. (6). The permissible region

*R* was initially chosen to be the whole object except for nodes within one transport length

of the surface. At each iteration, the permissible region is slightly shrunk by an arbitrary factor β. For a permissible region initially containing 1050 mesh nodes (as in the current problem) that is finally shrunk to 10 nodes in 60 iterations, β = (1050/10) ^ (1/60) = 1.081. The 60 iterations took approximately 7 minutes on a pc with Intel processor with a clock of 2.66 GHz and 2GB RAM. The execution time is less than that required by our previous algorithm in [

18].

shows the actual and reconstructed bioluminescence sources. The magnitudes of the three sources are 1 nW/mm^{3}, and the total powers of the three sources are 5.65, 5.65, and 40.71 nW, respectively. The transverse, sagittal, and coronal cross sections of the reconstructed sources in show good agreement with the actual sources in terms of source magnitude, total power and size. The total power of the reconstructed sources are 5.31, 5.49, and 41.4 nW. The source magnitude varies between 0.8 to 1.2 nW/mm^{3} which is within 20% of the true value.

The evolution of the objective function with the number of source nodes in the permissible region is illustrated in
. shows the objective function

as a function of the number of source nodes in the permissible source region. The permissible region has been shrunk in 60 iterations by removing nodes corresponding to low source magnitude, assuming that source nodes with low magnitude are less likely to be in the optimized source region. and show the reconstructed source distribution in nW/mm

^{3} after the first iteration which includes a permissible region of 1050 nodes and at the optimized permissible region which includes 126 nodes and minimizes the objective function. After the first iteration, the reconstructed source powers in the two kidneys and bone are 3.37, 3.8, and 51.93 nW, respectively which are within 40% of their actual values. The source magnitudes and sizes cannot be reconstructed perfectly after the first iteration due to the degeneracy of the solution. When the size of the source permissible region becomes comparable to the size of the actual source, the objective function has its lowest value and good accuracy in reconstructing the source power, magnitude, and size can be obtained, as shown in . When the permissible region is further reduced, the size of the source permissible region becomes smaller than the size of the actual source. Hence, the objective function increases again and the reconstruction of the source power, magnitude and size becomes less accurate.

shows the BLT reconstruction for multiple sources with different magnitudes and sizes. The sources in the left and right kidneys are 2 and 3 nW/mm

^{3} respectively, and the source in the bone is 1 nW/mm

^{3}. The three sources have total powers of 11.31, 16.96, and 40.71 nW. The reconstructed sources in have powers of 11.56, 18.81, and 38.92 nW respectively, which are within 11% of the true value. The source magnitudes are reconstructed within 20% of the true values. As a comparison with the result obtained using our previous algorithm in [

18], the maximum errors in the reconstructed source power and magnitude are 16% and 60% respectively.

shows the BLT reconstruction for multiple sources with different magnitudes and sizes. In this case, one source is greater in both size and magnitude and hence has much larger total power than other sources and dominates the measured light fluence rate. The sources in the left and right kidneys are 1 nW/mm^{3} and 2 nW/mm^{3} respectively and the source in the bone is 3 nW/mm^{3}. The three sources have total powers of 5.65, 11.31, and 122.14 nW respectively. The reconstructed sources are shown in . The reconstructed sources powers are 8.7, 10.94, and 120.14 nW respectively. The size and magnitude of the smallest source in the left kidney could not be reconstructed in this case. When there are large differences between source powers or there is a large inhomogeneous source distribution, the shrinking of the permissible region cannot be done correctly. The assumption that nodes with low source magnitude are least likely to be in the optimized source region is no longer true.

shows the BLT reconstruction for a large non-uniform source that fills the entire left kidney. There is a 3mm diameter “hot spot” of 2 nW/mm

^{3} in the center of the organ and a uniform “background” of 1 nW/mm

^{3} elsewhere. The total powers of the actual and reconstructed sources integrated over the kidney volume are 53.67 and 53.71 nW. The reconstructed source in and is in good agreement with the actual source in terms of size and total power. The inhomogeneity of the source distribution could be reconstructed and the high spot in the center of the kidney could be distinguished from the background. The magnitude of the reconstructed source at the center of the hot spot is within 20% of its actual value. As a comparison with the result obtained using our previous algorithm in [

18], the inhomogeneity of the reconstructed source could not be distinguished. An average value for the source magnitude was obtained rather than the hot spot shown in .

shows the BLT reconstruction for a non-uniform source distribution where the inhomogeneity of the source is high. In this case, the hot spot has a magnitude of 5 nW/mm^{3} which is 5 times that of the surrounding background. The total powers of the actual and reconstructed sources are 70.63 and 70.8 nW. However, the inhomogeneity of the source could not be reconstructed correctly in this case. The uniform background shows a smaller size and a higher magnitude, which varies from 1.5 to 2.5 nW/mm^{3}, than the actual source. The objective function curve with respect to the number of source nodes for the reconstructed source shown in does not show a deep valley with a clear global minimum as shown in . Instead, the curve is slowly increasing after the global minimum indicating that different source distributions can generate the same error and hence are indistinguishable. Therefore, shrinking the permissible region reduces the possible solutions, and if the permissible region equals the actual source size, a unique solution can be obtained for uniform and low inhomogeneity source distributions as shown in , , and . However, for high inhomogeneity source distributions as in and , the solution is still not unique and different size and magnitude sources that have the same total power can generate the same light fluence rate at the object boundary. One way to improve the uniqueness of the solution is to increase the available data points by doing measurements at more wavelengths.