1. INTRODUCTION

Optical tomography utilizes near-infrared light to image the optical properties of biological systems [1,2] in a manner that is directly related to tissue function and structure. The use of fluorescent molecular probes has the potential to increase the sensitivity and specificity of the method [3–5]. Fluorescence optical tomography aims to reconstruct the concentration and lifetime of the fluorescent contrast agents. The conventional approaches to fluorescent tomography require the knowledge of the intrinsic optical parameters at both the excitation and emission wavelengths. The latter can be determined separately, but this typically requires some assumptions about the spectral dependence of the absorption and scattering coefficients. As a result, inaccuracies arise in the reconstructed images due to the use of simplistic models of the tissue as well as due to errors in aligning the images provided by different techniques into a single integrated image.

We have recently shown theoretically that fluorescence imaging can be achieved in weakly scattering systems whose size is of the order of a few transport mean free paths [6]. Systems satisfying this condition include epithelial tissue layers, engineered tissue, and some model organisms [7]. Several techniques have been developed for tomographic imaging in this regime [7–20]. However, the wavelength was assumed to be constant in these studies. Inelastic scattering was addressed in [21], where the inversion formula relied on the known angle-dependence energy in Compton scattering and it was assumed that the attenuation coefficient of the medium depends on the energy linearly at each point. In [6], we introduced a tomographic imaging technique, referred to as non-reciprocal broken ray tomography (NRBRT), which enables simultaneous imaging of both the intrinsic and extrinsic optical parameters of a weakly scattering medium [6]. NRBRT accounts for the change in wavelength upon scattering and does not require any assumptions about the spectral properties of the tissue. In the presence of a fluorescent contrast agent, incident light is absorbed by the contrast agent and re-emitted in a different direction and at a different wavelength. In this setting, NRBRT makes use of angularly selective measurements of the fluorescence light exiting the medium at a fixed angle relative to the incident direction. Here the path of light as it travels from the source to the detector is assumed to be a broken ray where the wavelength is changed at its vertex. Image reconstruction in NRBRT involves inversion of the broken-ray transform (BRT)—a generalization of the Radon transform to the case of broken rays, relating measurements to line integrals along broken rays. As a result of the change in the wavelength that occurs at the vertex of the broken ray, different measurements are obtained when the source and detector positions are interchanged. In other words, the resulting BRT is non-reciprocal, meaning that the measurements depend on the direction in which the broken ray is traversed. We have shown that, by exploiting the non-reciprocity of the BRT and by using a measurement geometry with three coplanar broken rays with a common vertex, it is possible to derive inversion formulas for simultaneous reconstruction of the intrinsic attenuation coefficient of the medium at both the excitation and the fluorescence emission wavelengths, as well as the concentration of the fluorophore.

The inverse problem of NRBRT is two-dimensional, since any broken ray can be confined to a plane. Thus, three-dimensional image reconstruction can be performed slice by slice. Moreover, inversion formulas for the BRT can be applied to both backscattering and transmission measurement geometries. Importantly, we note that NRBRT does not require rotating the imaging device around the sample and does not require any assumptions about the spectral dependence of the attenuation coefficient.

In this paper, we demonstrate image reconstruction in NRBRT by performing numerical experiments. For this purpose, forward data are generated by Monte Carlo simulations in inhomogeneous samples characterized by varying scattering strengths and optical parameters for both transmission and backscattering measurement geometries. Although the single-scattering approximation is employed in deriving the reconstruction formula of the NRBRT, in the simulations reported below we do not assume that only singly scattered light is measured by the detectors. Therefore, our simulations are designed to be representative of a potential practical implementation of NRBRT in which all fluorescence light within the geometrical and collimation constraints of a detector is measured.

The rest of this paper is organized as follows. In Sections 2 and 3, we review the principles of NRBRT and present the derivation of the forward data function for the case of finite-size detectors, as well as the image reconstruction formulas used in this study. In Section 4, we describe the simulation of intensity measurements using the Monte Carlo technique. Image reconstructions are presented in Section 5, and Section 6 contains a discussion and conclusions.

2. PRINCIPLES OF NRBRT

In what follows, we review the principles of NRBRT, which have been introduced in [6]. We consider the problem of fluorescence tomography of samples characterized by spectrally and spatially dependent intrinsic absorption and scattering coefficients, ${\mu _{\texttt a}}$ and ${\mu _{\texttt s}}$, and the spatially dependent fluorophore concentration $n$. The presence of the fluorophore results in additional absorption, so that the total attenuation coefficient of the sample is $\mu= {\mu _{\texttt a}} + {\mu _{\texttt s}} + n\sigma$, where $\sigma$ is the spectrally dependent fluorophore absorption cross section. NRBRT reconstructs slice by slice the total attenuation coefficient at the excitation and fluorescence wavelengths, ${\mu _{\texttt e}}$ and ${\mu _{\texttt f}}$, and the fluorophore concentration $n$. It is based on angularly selective measurements of the fluorescence light intensity on the surface of the sample involving sets of three coplanar broken rays for a given slice of the medium. To implement this, three distinct, fixed directions are used for illumination and detection. As illustrated in Fig. 1, the slice under consideration is illuminated in one of the three directions by collimated, monochromatic sources at a wavelength within the fluorophore resonant absorption band, and the fluorescence photons emerging from the slice in the other two directions are registered. The source and detector positions are chosen so that, for each source, there are sets of three broken rays, corresponding to the source-detector and detector-detector pairs, which lie in the selected slice and have a common vertex. An example of such set of three broken rays is illustrated in Fig. 2. Moreover, for a given region of interest, the number of sources and detectors are such that the region of interest is fully scanned, i.e., all pixels of the discretized grid of the region of interest contain the vertices of three broken rays corresponding to a source and detectors at the two detection angles.

 

Fig. 1. Numerical sample and scanning of a region of interest in a slice with multiple sources and detectors in the transmission measurement geometry. The sphere represents the inhomogeneity and the small cuboid the region of interest (slice) in which the reconstruction is performed. The arrows pointing into the sample represent the incident photons, and those pointing out of the sample the detected fluorescence photons.

Download Full Size | PDF

 

Fig. 2. Two measurement geometries for a set of one source and two detectors. The three broken rays corresponding to each of the source-detector and detector-detector pairs are presented by shaded lines. The numbers 1, 2, 3 label the three fixed directions used in this study.

Download Full Size | PDF

Two measurement geometries can be implemented as is illustrated in Fig. 2. To generate a complete data set, scans are performed for a given configuration of source and detector positions and directions, and then the sources and detectors are interchanged, and the sample is scanned again, in each case using the same positions and directions. This results in three scanning configurations for each measurement geometry as is shown in Fig. 3. For each measurement geometry, all the scanning configurations where the sources and detectors are interchanged, and the positions of the sources and detectors varied, represent a complete scan. We note that one of the complete scans consists of both transmission and backscattering measurements. We refer to this as the transmission scan. The other complete scan consists of only backscattering measurements, and we refer to this as the backscattering scan.

 

Fig. 3. Three measurement configurations for the transmission scan for a set of one source and two detectors. In this study, directions 2 and 3 are mutually orthogonal and are at 135º relative to 1.

Download Full Size | PDF

3. THEORY

A. Forward Data Function

In [6], the forward data function was derived for the case of infinitely narrow incident beams and infinitely small detectors. In this study, we consider finite-size and finite-collection angle detectors. The signal measured by an angularity selective detector at ${{\textbf{r}}_{\texttt{d}}}$ collimated in the direction ${{\hat{\textbf s}}_k}$ and having an acceptance angle $\beta$, for a source at ${{\textbf{r}}_{\texttt{s}}}$ and in the direction ${{\hat{\textbf s}}_l}$, can be expressed as

Here $W$ is the source power per unit area; ${\textbf{R}} = {\textbf{r}} - [|{\textbf{r}} - {{\textbf{r}}_{\texttt s}}|\sin {\theta _l}/\sin ({\theta _l} + \theta)]{\hat{\textbf s}}$ is the position of the vertex of the broken ray defined by the source and detector pair; $\tilde n({\textbf{R}}) = {\sigma ^{3/2}}n({\textbf{R}})$; ${\theta _l}$, ${\varphi _l}$ and $\theta$, $\varphi$ are the polar angles of the unit vectors ${{\hat{\textbf s}}_l}$ and ${\hat{\textbf s}}$, respectively, defined in a reference frame whose $z$ axis coincides with the line of sight ${\textbf{r}} - {{\textbf{r}}_{\texttt s}}$, and $\delta (\cdot)$ is the Dirac delta function. ${\mu _{\texttt e}}$ and ${\mu _{\texttt f}}$ are the total attenuation coefficients at the excitation and fluorescence wavelengths, and ${L_1} = |{\textbf{R}} - {{\textbf{r}}_{\texttt s}}|$ and ${L_2} = |{\textbf{r}} - {\textbf{R}}|$ represent the lengths of the two rays making up the broken ray defined by the source-detector pair. Note that Eq. (2) was derived under the assumption that ${{\hat{\textbf s}}_l}$ and ${\textbf{r}} - {{\textbf{r}}_{\texttt s}}$ lie in the same plane [6]. The presence of the delta function and step function ensure that ${\hat{\textbf s}}$ lies in the same plane as well, and that the rays in the directions ${\hat{\textbf s}}$ and ${{\hat{\textbf s}}_l}$ intersect in that plane. Using Eqs. (1) and (2), we obtain

 

Fig. 4. Illustrating the geometry relevant to Eq. (3).

Download Full Size | PDF

For small detector acceptance angles, we can approximate ${{\textbf{R}}^\prime}({\textbf{r}},\theta) \simeq {{\textbf{R}}^\prime}({\textbf{r}},{\theta _0})$, where ${\theta _0}$ is the angle between ${{\hat{\textbf s}}_k}$ and ${\textbf{r}} - {{\textbf{r}}_{\texttt{s}}}$. This approximation assumes that the sample parameters are constant in the domain within the sample probed by the signal measured at ${\textbf{r}}$ (the intersection between the cone at ${\textbf{r}}$ and the source ray). This is a good approximation apart from areas of large variations of the optical parameters such as an inhomogeneity boundary. By using this approximation, we obtain

 

Fig. 5. Illustrating detection along the direction ${{\hat{\textbf s}}_k}$ with an acceptance angle $\beta$ of photons emitted by fluorescence sources along the line defined by the incident source direction ${{\hat{\textbf s}}_l}$. For a small distance between the fluorescence source and the detector, each fluorescence source contributes to the detected signal with photons detected over areas smaller than the detector area, and not all the extent of the detector area in the $x$ direction is crossed by photons.

Download Full Size | PDF

 

Fig. 6. Geometrical factor ${f_{31}}$ as a function of the position $(x,y)$ on the detector for the transmission measurement. The source line passes through the center of the sample, and the source-detector separation on $y$ direction is, from left to right, 4.85 cm, 2.55 mm, and 0.25 mm. The values have been normalized in each case, and the other parameters are ${L_z} = 5.1\;{\rm{mm}}$, $\beta = 0.5\deg$, ${\Delta _{\texttt{d}}} = 0.07\;{\rm{mm}}$.

Download Full Size | PDF

Further, ${{\textbf{R}}^\prime}$ is approximated to be the same as the vertex position ${\textbf{R}}$ defined by the source position and direction ${{\textbf{r}}_{\texttt s}}$, ${{\hat{\textbf s}}_l}$ and the center and collimation direction of the detector ${{\textbf{r}}_{\texttt d}}$, ${{\hat{\textbf s}}_k}$. Also, we approximate the path of the fluorescence photons from the vertex position to any point ${\textbf{r}}$ on the detector surface with the path from the vertex position to the center of the detector. This is a good approximation apart from areas of large variation in the sample parameters and situations where the vertex position is close to the surface of the sample where the detector is positioned. With these approximations, we obtain

Here ${\bar L_{1,2}}$ are the lengths of the two rays making up the broken ray defined by ${{\textbf{r}}_{\texttt{s}}}$, ${{\hat{\textbf s}}_l}$ and ${{\textbf{r}}_{\texttt{d}}}$, ${{\hat{\textbf s}}_k}$, and

Based on Eq. (5), we define the projection data function corresponding to a source at position ${{\textbf{r}}_{\texttt{s}}}$ and in direction ${{\hat{\textbf s}}_l}$ and a finite-size detector at position ${{\textbf{r}}_{\texttt{d}}}$ and in direction ${{\hat{\textbf s}}_k}$ as

We note that, with the data function defined in this way, the forward model has the same form as the forward model for the case of point detectors [6]. Therefore, the image reconstruction formalism derived in [6] can be directly applied here by using the data function defined by Eq. (7).

B. Image Reconstruction

The sample attenuation coefficient at the excitation and fluorescence wavelengths ${\mu _{\texttt{e}}}$ and ${\mu _{\texttt{f}}}$ are reconstructed using the image reconstruction formulas derived in [6] for the spectrally averaged attenuation coefficient ${\mu ^{(+)}} = ({\mu _{\texttt{e}}} + {\mu _{\texttt{f}}})/2$, the spectral difference of the attenuation coefficients ${\mu ^{(-)}} = ({\mu _{\texttt{e}}} - {\mu _{\texttt{f}}})$, and $\tilde n$. The reconstruction formulas for ${\mu ^{(\pm)}}$ are

Here ${\sigma _{1,2,3}}$ are constants determined based on the three source and detector directions, and ${\Phi ^{(\pm)}}$ are linear combinations of the symmetric and anti-symmetric data functions $\phi _{\textit{lk}}^{(\pm)}$, with the coefficients of the linear combinations being two-dimensional vectors in the slice plane and depending on the source and detection directions. The details for computing the functions ${\Phi ^{(\pm)}}$ and the constants ${\sigma _{1,2,3}}$ are presented in [6]. The symmetric and anti-symmetric data functions $\phi _{\textit{lk}}^{(\pm)}$ are defined as the average and difference of the two datasets measured for a source-detector configuration and for the configuration where the sources and detectors are interchanged, $\phi _{\textit{lk}}^{(+)} = ({\phi _{\textit{lk}}} + {\phi _{\textit{kl}}})/2$, $\phi _{\textit{lk}}^{(-)} = {\phi _{\textit{lk}}} - {\phi _{\textit{kl}}}$. The data functions ${\phi _{\textit{lk}}}$ are given by Eq. (7) with the measured signal $\delta {W_{\textit{lk}}}$ computed according to Eq. (14). Once ${\mu ^{(\pm)}}$ are determined, the attenuation coefficient at the excitation and fluorescence wavelengths are obtained as ${\mu _{\texttt{e,f}}}({\textbf{R}}) = {\mu ^{(+)}}({\textbf{R}}) \pm (1/2){\mu ^{(-)}}({\textbf{R}})$.

We note that the reconstructions formulas in Eq. (8) are local, which enables reconstructing the attenuation coefficients in the domain where fluorescence contrast agent is present [6]. Also, due the differentiation in Eq. (8), the image reconstruction of ${\mu ^{(\pm)}}$ and ${\mu _{\texttt{e,f}}}$ does not depend on the constants in Eq. (7) for the data function and in the expression of measured signal (presented in Section 4), and therefore the image reconstruction is quantitative.

The reconstruction formula for the concentration of the fluorescent contrast agent reads [6]

Here the index 1 corresponds to the direction perpendicular to the sample surface, and 2, 3 to the other two directions as illustrated in Fig. 2. Thus, $I_1^{(+)}({\textbf{R}})$ in Eq. (11) is calculated along a line that is perpendicular to the sample surface and is always contained in the region where ${\mu ^{(+)}}$ is known. In this study, Eq. (11) was used to reconstruct the fluorophore concentration. Note that this formula involves (through the symmetric data functions $\phi _{\textit{lk}}^{(+)}$) pairs of measurements, with sources and detectors in directions $l$ and $k$, as well as with the sources and detectors interchanged and having the directions $k$ and $l$, respectively. The reconstruction formulas derived in [6] can also be applied to the case of conventional fluorescence optical tomography. In this case, the attenuation coefficients ${\mu _{\texttt{e,\!f}}}$ of the medium are known, and the problem is to reconstruct just the fluorophore concentration, which can be accomplished by using the formula

Note that in this case only one source-detector measurement configuration $lk$ is needed, with the sources in direction $l$ and detectors in direction $k$, and without the need to interchange the sources and detectors. Therefore, image reconstructions can be obtained for any of the six $lk$ pairs, and in practice only one such experimental setup would be needed. If measurements are available for a given source-detector configuration and for the configuration where the source and detector are interchanged, Eq. (9) can also be applied to the case of conventional fluorescence optical tomography. This could potentially lead to better reconstructions, as it can be seen in the results section. $I_{l,k}^{(+)}({\textbf{R}})$ in Eq. (9) are calculated in this case as the line integrals of the known average attenuation coefficient computed according to Eq. (10).

Table 1. Sizes of the Samples Used in This Study^a

View Table | View all tables in this article

Note that the reconstructed values of $\tilde n$ depend on the constants in the expression in Eq. (7) of the data function through a global multiplication constant, and quantitative reconstruction is not possible without calibration. The images of the reconstructed fluorophore concentration $\tilde n$ presented in the results section correspond to the values obtained based on Eqs. (11), (12), or (9) and normalized to the maximum value of the particular image.

4. SIMULATIONS

Projection data used for image reconstruction were generated using the Monte Carlo software package Geant4-based Architecture for Medicine-Oriented Simulations (GAMOS) [22,23] and the Tissue Optics plugin [24]. Monte Carlo simulations account in great accuracy for all relevant processes. Specifically, this software package was used to simulate light transport in inhomogeneous optical samples containing fluorescence contrast agents (fluorophore) and the detection of fluorescence light intensity emerging from the sample.

This study utilized three-dimensional samples with the overall shape of cuboids consisting of a homogeneous background and spherical inclusions representing inhomogeneities, as is illustrated in Fig. 1. The sample scattering strength can be quantified by the optical depth, ${L_z}{\bar\mu_{\texttt{s}}}$, where ${\bar\mu_{\texttt{s}}}$ is the intrinsic scattering coefficient of the background and ${L_z}$ is the sample length in the $z$ direction. Further, we considered isotropic light scattering, and, to improve computational efficiency, the fluorescent molecules were assumed to have the quantum efficiency of one, and fluorophore absorption and emission spectra were modelled as delta functions. The latter is equivalent to monochromatic excitation and narrowly spectrally filtered detection. Also, the absorption of the fluorescence light by the fluorophore and refraction of rays at the sample boundary have been neglected. Numerical experiments have been performed with a variety of samples characterized by various scattering strengths and contrasts between the inhomogeneity and the background. The sizes of the samples used in this study are presented in Table 1, and the optical characteristics in Table 2.

Each sample was discretized, first, by representing it as a stack of slices drawn parallel to the ${\textit{yz}}$ plane and of the thickness 0.1 mm. In the numerical examples below, we only show results for the slice containing the center of the inhomogeneity, but reconstructions in other slices can be similarly obtained. Further, each slice was discretised into voxels of the size $0.1 \times 0.1 \times 0.1\;{\rm{mm}}^3$. The sources and detectors used for scanning the slice have been placed on a line parallel to the $y$ axis on the surface of the sample. The source and detector positions were evenly separated for all configurations, with the space between these positions of 0.1 mm. One of the three scanning directions was normal to the surface of the sample. The other two directions were at angles $\alpha = \pm {45^ \circ}$ relative to the normal to the surface, for the transmission scan, and at 45º and ${-}{70^ \circ}$, for the backscattering scan. The scan type and number of source and detector positions for each sample and direction are listed in Table 1. For each source, ${10^{11}}$ photons were incident onto the medium. This was implemented by performing 1000 simulations with ${10^8}$ incident photons for each source. Each simulation used a different seed for the pseudorandom number stream and a different position in the stream for the value used to initialize the simulation. The output data were calculated as the sum of the outputs of all simulations.

Table 2. Characteristics of the Samples Used in the Simulations^a

View Table | View all tables in this article

To simulate light-intensity detection by angularly selective detectors represented whose surface is the region $\delta \Sigma$ and collimated in the direction defined by the angle $\alpha$, we have scored the fluorescence photons on the surface of the sample in the region $\delta {\Sigma ^\prime}$ representing the projection of the detector onto the sample surface as shown in Fig. 7. Note that there are photons scored close to the edge $\delta {\Sigma ^\prime}$ that do not cross $\delta \Sigma$ and also photons that cross $\delta \Sigma$ that do not cross $\delta {\Sigma ^\prime}$. However, the fraction of $\delta {\Sigma ^\prime}$ that scores photons not crossing $\delta \Sigma$, and the fraction of the area $\delta \Sigma$ crossed by photons not recorded by $\delta {\Sigma ^\prime}$ are of the order of $(1/2)\tan \alpha \sin \beta \simeq {10^{- 2}}$ for the typical values of $\beta$. For all samples except Samples 2, 3, 4, and 12, we have considered measurements with square detectors of fixed area $\Delta _{\texttt d}^2$. Correspondingly, photons were scored on the sample surface on rectangles (for detection angle $\alpha \ne 0$) of sides ${\Delta _{\texttt{d}}}$ and ${\Delta _{\texttt{d}}}/\cos \alpha$, in the $x$ and $y$ direction, respectively, or squares (for $\alpha = 0$). The variable ${\Delta _{\texttt{d}}}$ in Section 5 is the side of the squares $\delta \Sigma$. This corresponds to using the same detector for all directions. In contrast, for the backscattering scans of Samples 2, 3, 4, and 12, for all directions, photons were scored on the sample surface on squares of fixed area $\Delta _{\texttt{d}}^2$, and ${\Delta _{\texttt d}}$ referred to in Section 5 is the side of these squares. This corresponds to using different detector sizes for different directions. Image reconstruction was better for these samples with such measurements. This is so because, as discussed in Section 6, the data generated by Monte Carlo codes are noisier for backscattering measurements, which affects image reconstruction. The detector configuration that was used for samples 2, 3, 4, and 12 increases the signal-to-noise ratio while keeping the computational time within a reasonable limit.

 

Fig. 7. Photon scoring. $\delta \Sigma$ represents the collimated detector and $\delta {\Sigma ^\prime}$ the projection of the detector onto the sample surface where the photons are scored.

Download Full Size | PDF

For a source at position ${{\textbf{r}}_{\texttt{s}}}$ and illuminating in the direction ${{\hat{\textbf s}}_l}$, we recorded the number $N({{\textbf{r}}_i},{{\hat{\textbf s}}_j};{{\textbf{r}}_{\texttt{s}}},{{\hat{\textbf s}}_l};\beta)$ of photons (hits) at the fluorescence wavelength emerging from the slice at points ${{\textbf{r}}_i}$ on the surface and propagating in the direction ${{\hat{\textbf s}}_j}$ within a cone of half-angle $\beta$ and axis in the direction ${{\hat{\textbf s}}_k}$ and with the apex at ${{\textbf{r}}_i}$, where ${{\hat{\textbf s}}_{l,k}}$ are two of the three distinct directions considered in this study. The mesh for the points ${{\textbf{r}}_i}$ is fixed in GAMOS and is ${10^{- 7}}\;{\rm{mm}}$. Using Eq. (1), the signal measured by an angularly selective detector positioned at ${{\textbf{r}}_{\texttt{d}}}$ of area $\delta \Sigma$ and acceptance angle $\beta$ and collimated in the direction ${{\hat{\textbf s}}_k}$ was computed as

Table 3. RMSE and SSIM between the Reconstructed Images and the Corresponding Model for Different Numbers of Incident Photons per Source Position N for Sample 1, Using the Transmission Scan with $\beta ={ 0.5^ \circ}$, ${\Delta _{\texttt d}} = 0.070 \;{\rm{mm}}$ and $\lambda = 5$ for $N={ 10^9}$ and $N={ 10^{10}}$, and $\lambda = 20$ for $N = 5 \times {10^{10}}$

View Table | View all tables in this article

Implementing the reconstruction formulas in Eq. (8) requires numerical differentiation of the functions ${\Phi ^{(\pm)}}$. Projection data generated using Monte Carlo simulations are inherently noisy. In addition, due to the presence of sharp-edge inhomogeneities, the data functions are not smooth. Numerical differentiation methods (for example, finite-difference approximations) are known to amplify the noise in the data [21,25,26]. Denoising the data prior to differentiation was found not to always solve the problem [25]. Therefore, a technique based on the regularisation of the derivative [26] was used in this study. This method uses total-variation regularization and is suitable for noisy and discontinuous functions. We used this method following the application of a median filter to ${\Phi ^{(\pm)}}$. The total-variation regularization method is based on using periodic boundary conditions for Fourier-domain processing, which resulted in artefacts in the reconstruction at the edges of the region of interest [26]. To reduce these artefacts, after applying the median filter, the functions ${\Phi ^{(\pm)}}$ were extrapolated to relocate the area most affected by the inherent artefacts of the derivative method outside the region of interest. Specifically, 10 pixels at the edge of the region of interest on all sides were duplicated outside the region of interest. Furthermore, to improve the accuracy of the derivative computation (affected by large local variations of the data functions), the regularization method was applied using as base derivative method both the forward and backward finite difference approximations, and the derivatives in Eq. (8) have been calculated using the second-order finite difference approximation. The value of the regularization parameter $\lambda$ in the formalism can be varied to obtain an optimal reconstructed image quality.

In implementing the reconstruction formulas from Eqs. (11), (12), and (9), a median filter was applied to the all the data functions ${\phi _{\textit{lk}}}$ used and to the reconstructed ${\mu ^{(+)}}$, before the line integrals and $\tilde n$ have been computed.

5. RESULTS

In this section, we present image reconstructions for the various samples described above. For all the figures, the horizontal and vertical axes of the images correspond to the $y$ and $z$ axes in Fig. 3, respectively, and the same color scale is used for the reconstructed quantities and the corresponding models. The color scale for ${\mu ^{(\pm)}}$, ${\mu _{\texttt{e}}}$, and ${\mu _{\texttt{f}}}$ has been clipped in the range of 0 to 1.5 times the maximum value in the model. The parameters $\beta$, ${\Delta _d}$, and $\lambda$ in each case have been chosen to obtain optimal quality of the reconstructed images, based on visual inspection. Further, the quality of these reconstructed images is assessed using the root mean square error (RMSE) and the structural similarity index measure (SSIM) [27], measuring the quantitative and perceptual difference, respectively, between the reconstructed image and the model. Small (close to zero) values of RMSE and values of SSIM close to one indicate good agreement. Table 3 presents the reconstruction errors as a function of the incident number of photons per source, for Sample 1. Based on this, we obtain that acceptable reconstruction is possible for ${10^{10}}$ or more incident photons. The reconstructions presented below correspond to ${10^{11}}$ photons per source, and the error metrics for these reconstructions are shown in Table 4.

Table 4. RMSE and SSIM between the Reconstructed Image and the Corresponding Model for the Reconstructions Presented

View Table | View all tables in this article

The reconstructed optical parameters of Sample 1, characterized by an optical depth of less than 1, based on the transmission scan, are presented in Fig. 8 for varying detector acceptance angles and detector areas. We obtain a good agreement between the reconstructed quantities and the model. The reconstructed image quality is better for larger detector acceptance angles and sizes. Figure 9 presents the reconstructed optical characteristics of Samples 2, 3, and 4, based on the backscattering scan. These samples are all characterized by the same optical depth of 0.7, but they have different sizes and contrasts in optical characteristics between the inhomogeneity and the background as presented in Table 2. We obtain a reasonable agreement between the reconstructed quantities and the model, but the image quality in this case is lower than it is for transmission scans, particularly for regions deeper into the sample and behind the inhomogeneity.

 

Fig. 8. Reconstructed optical parameters of Sample 1 using the transmission scan with ${\Delta _{\texttt d}} = 0.070 \;{\rm{mm}}$ and $\beta ={ 0.75^ \circ}{,0.25^ \circ}$; $\beta ={ 0.5^ \circ}$ and ${\Delta _{\texttt d}} = 0.070 \;{\rm{mm}}$, $0.028 \;{\rm{mm}}$. For all cases, $\lambda = 50$.

Download Full Size | PDF

 

Fig. 9. Image reconstruction for Samples 2 (left), 3 (middle), and 4 (right) based on the backscattering scan with $\beta ={ 0.5^ \circ}$ and ${\Delta _{\texttt d}} = 0.1 \;{\rm{mm}}$, and for $\lambda = 30$, 100, 200, respectively.

Download Full Size | PDF

Next, we consider the image reconstruction of stronger scattering samples. Figure 10 presents the reconstructed optical parameters of Samples 5, 6, 7, and 8. All these samples are characterized by the same optical depth of 1.5, but have different sizes, contrast agent concentration in the background, and contrast in optical parameters between the inhomogeneity(ies) and the background. Specifically, Sample 6 is larger and has a larger inhomogeneity than Sample 5, and a lower concentration of the contrast agent in the background; Sample 7 has the same characteristics as Sample 5, apart from a 10 times less contrast agent concentration in the background and a higher contrast between the inhomogeneity and the background for the fluorescence contrast agent; Sample 8 has an even lower concentration of the contrast agent in the background. We obtain good quality images for transmission scans and poorer quality for backscattering scans. Also, the image quality degrades for lower fluorescence contrast agent concentration in the background. However, the inhomogeneity is still visible in all reconstructions.

 

Fig. 10. Reconstructed images for (from left to right) Sample 5 using the transmission scan with $\beta ={ 0.25^ \circ}$, ${\Delta _{\texttt d}} = 0.070\;{\rm{mm}}$, and $\lambda = 250$; Sample 6 using the transmission scan with $\beta ={ 0.5^ \circ}$, ${\Delta _{\texttt d}} = 0.070\;{\rm{mm}}$, and $\lambda = 50$; Sample 7 using the transmission scan with $\beta ={ 0.25^ \circ}$, ${\Delta _{\texttt d}} = 0.070\;{\rm{mm}}$, and $\lambda = 100$; Sample 7 using the backscattering scan with $\beta ={ 0.5^ \circ}$, ${\Delta _{\texttt d}} = 0.034\;{\rm{mm}}$ and $\lambda = 25$; and Sample 8 using the transmission scan with $\beta ={ 0.25^ \circ}$, ${\Delta _{\texttt d}} = 0.070\;{\rm{mm}}$ and $\lambda = 50$. T and B indicate the transmission and backscattering scan, respectively.

Download Full Size | PDF

The scattering strength was further increased to an optical depth of $3$, for which the single-scattering approximation used in deriving the image reconstruction algorithm is expected to start breaking down. Figure 11 presents the reconstructed optical parameters based on the transmission scan of Samples 9 and 10, differing by their size and contrast agent concentration in the background and by the optical depth of the inhomogeneity. While the quality of the reconstructed images has deteriorated compared to weaker-scattering samples for the same type of scan, the inhomogeneity can still be distinguished with good accuracy, particularly if ${\mu ^{(+)}}$, ${\mu _{\texttt{e}}}$, and $n$ are used as markers.

 

Fig. 11. Image reconstruction based on the transmission scan with ${\Delta _{\texttt d}} = 0.070 \;{\rm{mm}}$ for Sample 9 (left) for $\beta ={ 0.25^ \circ}$ and $\lambda = 50$ and Sample 10 (right) for $\beta ={ 0.5^ \circ}$ and $\lambda = 25$.

Download Full Size | PDF

Consider now the case where the contrast agent is exclusively accumulated in the inhomogeneity. In this case, as discussed above, using the reconstruction formalism and the measurements considered in this study, it is possible to reconstruct only the coefficients ${\mu ^{(\pm)}}$ and ${\mu _{\texttt{e,\!f}}}$ within the domain where the contrast agent is accumulated (the inhomogeneity), but the contrast agent concentration cannot be reconstructed anywhere. RMSE and SSIM for these reconstructed images were computed only for the domain where reconstruction was possible.

Figure 12 presents the reconstructed images for Samples 11 and 12, which have the optical depth of 3, and Sample 13, with an optical depth of 4. The transmission scan was used for Samples 11 and 13 and the backscattering scan for Sample 12. We obtain that the quality of the reconstructed images, and the agreement between the reconstructed images and the model are reasonably good even for these stronger scattering samples.

 

Fig. 12. Image reconstruction for Sample 11 (left) using the transmission scan with $\beta ={ 0.25^ \circ}$, ${\Delta _{\texttt d}} = 0.070\;{\rm{mm}}$, and $\lambda = 50$; Sample 12 (middle) using the backscattering scan with $\beta ={ 0.5^ \circ}$, ${\Delta _{\texttt d}} = 0.1\;{\rm{mm}}$, and $\lambda = 100$; and Sample 13 (right) using the transmission scan with $\beta ={ 0.25^ \circ}$, ${\Delta _{\texttt d}} = 0.070\;{\rm{mm}}$, and $\lambda = 50$.

Download Full Size | PDF

Finally, we consider the case of conventional fluorescence optical tomography. Figure 13 presents the reconstructed images for Samples 7, 9, 11, and 13 based on Eq. (12) and each of the source-detector measurement configurations that were illustrated in Fig. 3. For samples with optical depth of 1.5 (Sample 7), very good image reconstruction is obtained for all source-detector configurations, with the inhomogeneity distinguished in all reconstructed images. For stronger-scattering samples, good image reconstruction is obtained for transmission measurements. For samples of the optical depth of 3, relatively good reconstruction is possible also for backscattering measurements if the contrast agent is exclusively accumulated in the inhomogeneity. This is also indicated by the SSIM value, which increases by nearly 5 times when comparing Sample 11 to Sample 9. Figure 14 presents the reconstructions for the same samples, but based on Eq. (9) and pairs of measurements, corresponding to a given source-detector configuration and to the configuration where the source and the detector are interchanged. This shows that additional measurements enable reducing the artefacts in the background and better identifying the inhomogeneity.

 

Fig. 13. Reconstruction of the fluorescence contrast agent distribution of, from left to right, Samples 7, 9, 11, amd 13, based on the source-detector configurations in Fig. 3, with $\beta ={ 0.25^ \circ}$ and ${\Delta _{\texttt d}} = 0.070 \;{\rm{mm}}$, respectively. The labels $lk$ with $l,k = 1,2,3$ on the left indicate the source-detector measurement setup with the source in direction $l$ and detector in direction $k$.

Download Full Size | PDF

 

Fig. 14. Reconstruction of the fluorescence contrast agent distribution of, from left to right, Samples 7, 9, 11, and 13, based on the pairs of source-detector configurations in Fig. 3 with $\beta ={ 0.25^ \circ}$ and ${\Delta _{\texttt d}} = 0.070 \;{\rm{mm}}$, respectively. The labels $lk$ & $kl$ with $l,k = 1,2,3$ on the left indicate that the projection data was the average data corresponding to the source-detector measurement configurations $lk$ and $kl$.

Download Full Size | PDF

6. DISCUSSION AND CONCLUSIONS

This computational study demonstrated simultaneous reconstruction based on NRBRT of three different quantities characterizing an optical sample, namely, the attenuation coefficient at the excitation and fluorescence wavelengths, ${\mu _{\texttt{e,\!f}}}$ (retrieved via the reconstruction of the average and difference attenuation coefficients ${\mu ^{(\pm)}}$), and the fluorescence contrast agent distribution $n$. The image reconstruction for ${\mu ^{(+)}}$ and ${\mu _{\texttt{e}}}$ was generally better than for ${\mu ^{(-)}}$ and ${\mu _{\texttt{f}}}$. This is so because it involves addition of data functions, while the reconstruction of ${\mu ^{(-)}}$ and ${\mu _{\texttt{f}}}$ involves subtraction of (positive) noisy functions, resulting in a reduced signal-to-noise ratio. Similarly, the reconstruction of $n$ was relatively stable. This is so because this reconstruction utilizes the symmetric data function (involving summation) and the relatively stable reconstruction of ${\mu _{\texttt{e}}}$. This study suggests that, for situations wherein the reconstruction of ${\mu _{\texttt{f}}}$ is of poor quality, ${\mu ^{(+)}}$ and ${\mu _{\texttt{e}}}$ and $n$ can be used as three independent markers, and therefore comprehensive assessment is still possible. We emphasize that, while the image reconstruction algorithm is based on the single-scattering approximation, the reconstructed images presented above were obtained from simulated data accounting for all orders of scattering.

 

Fig. 15. Measured signal, computed based on Monte Carlo simulations and the forward model, and the noise $\varepsilon$ in the simulated values as functions of the source-detector separation on the $y$ direction for Sample 1, for (a) and (c) transmission measurements and for (b) and (d) backscattering measurements. The noise was calculated as the percentage relative difference between the simulated signal value and the model signal value. The detector characteristics are $\beta ={ 0.5^ \circ}$, ${\Delta _{\texttt d}} = 0.070 \;{\rm{mm}}$.

Download Full Size | PDF

Image reconstruction based on the backscattering scan was not as stable as that for the transmission scan. This is so due to the higher levels of noise in the projection data for backscattering measurements. To illustrate this, we considered sets of transmission and backscattering measurements for Sample 1, corresponding to the source-detector pairs 31 and 32 presented in Fig. 3, respectively. The points in the sample probed by each set of measurements (the vertex positions) were the same (and positioned on a line parallel to the $z$ axis and passing through the center of the inhomogeneity). Figure 15 presents the measured signal $\delta {\bar W_{\textit{kl}}} = \delta {W_{\textit{kl}}}/\eta {E_{\texttt f}}$, computed using Monte Carlo simulations and Eq. (14), as a function of the source-detector separation. In this figure, we also present the (scaled) theoretical prediction for the measured signal, calculated using the forward model given in Eq. (5) and the values of ${\mu _{\texttt{e,\!f}}}$ and $n$ of the model. The theoretical signal is determined based on the assumption that no scattering event takes place apart from the inelastic scattering of photons interacting with the contrast agent. This is a good approximation for the weakly scattering Sample 1, and therefore no significant systematic difference between the theoretical and simulated signals is expected. Indeed, it can be seen from the data of Fig. 15 that the theoretical and simulated signals have the same behavior. However, the Monte Carlo signal is noisy, with the noise relatively high for the backscattering measurements, particularly, when probing positions deeper into the sample (which correspond to large source-detector separations). This strongly affects the image reconstruction based on a complete scan consisting only of backscattering measurements, such as the backscattering scan. There is also an effect even for the transmission scan, due to the backscattering measurement components (the set of measurements for the source-detector pairs 23 and 32 presented in Fig. 3), and this explains the artefacts above the inhomogeneity in the reconstructed images presented in this study. While numerical simulations and the levels of noise generated using the available Monte Carlo code and computational resources are not necessarily representative of the experimental situation, indeed, in practice, a lower signal-to-noise ratio (due to both lower signal and higher noise) is expected for the backscattering measurements than for the transmission measurements proposed here. In particular, the lower signal is due in turn to the longer photon paths in the medium, and there could be situations where no signal is measured. In this study, the incomplete dataset was corrected to some extent by using the median filter. Further improvements can be achieved by interpolation or extrapolation, or by appropriately choosing the angles of illumination and detection that determine the photon path length. Also, using sensitive enough detectors and optimal denoising techniques could, in principle, reduce the noise in the measured signal. All these in turn could enable better image reconstructions than in these numerical experiments. Image reconstruction for the transmission scan could also be improved by considering a trade-off for the photons path lengths between transmission and backscattering measurement components, by using different directions of illumination and detection than presented here.

We also note that the measurement geometry with two of the relevant directions at angles ${\pm}{45^ \circ}$ and the symmetric samples considered in this study for the transmission scan is the most stringent test for image reconstruction since, in this case, the anti-symmetric data functions used in reconstructions are characterized by a relatively small signal-to-noise ratio due to subtraction of two similar noisy quantities. Further improvement in image reconstruction can be achieved by suitably choosing the directions involved to avoid symmetry in the star formed by the three broken rays defined by the illumination and detection directions. In this study, image reconstruction for the backscattering scan was indeed improved by choosing the two directions to be at angles 45º and ${-}{70^ \circ}$. In practice, insight for the optimization of the directions for a particular sample can be obtained from a preliminary scan and reconstruction based on a chosen set of directions.

Some of the artefacts in the reconstructed images, particularly the streak artefacts, are caused by the sharp discontinuity in the optical characteristics, resulting in sharp variations in the detected signal and data functions as can also be seen in Fig. 15. Although the method used to compute the derivatives involved in the reconstruction formulas is suitable for discontinuous functions [26], its performance has been compromised by the presence of high levels of noise. However, in practice, there is usually a smoother variation in the optical characteristics than the step-like discontinuity considered here. This, together with possibly lower noise levels in the measured signal, would ameliorate the artefacts.

The detector size and acceptance angle both influence the quality of image reconstruction, and a reduction of noise in the measured signal and potentially better reconstruction can be obtained as both these detector characteristics are increased. However, increasing the detector acceptance angle increases the size of the domain within the sample probed by the signal measured at a particular point on the detector (represented by $|{{\textbf{P}}_1} - {{\textbf{P}}_2}|$ in Fig. 4), and there could be situations where the optical properties of the sample in that domain are not constant or even exhibit large variations. This possibility was ignored in the derivation of the reconstruction algorithm as discussed in Sec. 3. Therefore, as the detector acceptance angle is increased, the image reconstruction could deteriorate, particularly in areas of variations in the sample parameters. Thus, a trade-off between the measurement noise and the applicability of the forward model needs to be considered when deciding the detector parameters, to obtain the best reconstructions. For this reason, the best reconstructions for some of the samples studied here were obtained for smaller acceptance angles.

Image reconstruction for weakly scattering samples was relatively good. The reconstruction quality deteriorated as the sample optical depth of the background was increased above 4. The breakdown of the reconstruction formalism based on single-scattering approximation is indeed expected as the scattering strength is increased. This is exacerbated by the fact that the detector has a nonzero acceptance angle, and, as a result, the contribution to the detected signal of multiple scattered light could be significant, particularly for stronger scattering samples. In this numerical study, the detector acceptance angle needed to be above a certain value, to maintain a reasonable signal-to-noise ratio for the detected signal while keeping the number of simulations and the computational time within a reasonable limit. However, further reducing the acceptance angle could be possible in practice, potentially extending the applicability of NRBRT to even stronger scattering samples.

In samples with lower fluorescence contrast agent uptake in the background, the inhomogeneity was reconstructed to a good quality, whereas the reconstruction of the background presented artefacts. This was due to low signal-to-noise ratio for measurements probing the background of the sample, resulting from low emitted signal. The ability to reconstruct the inhomogeneity with good accuracy highlights the advantage of an imaging technique with a local reconstruction formula such as the one considered here.

In conclusion, using numerical experiments, this study has validated the principles and image reconstruction formalism of NRBRT applied to fluorescence optical tomography. Detailed guidance was also provided for experimental implementation, as well as possibly a motivation for advancements in areas such as sensitive detection or data denoising while maintaining features and discontinuities and accurate derivative calculation for such data. Future work will include further development and validation of NRBRT to the case of contrast agent exclusively accumulated in the inhomogeneity.