Multi Simulation Platform for Time Domain Diffuse Optical Tomography: An Application to a Compact Hand-Held Reflectance Probe

Featured Application

The present work represents a preliminary performance assessment of a simulated system for Time Domain Diffuse Optical Tomography, in view of the upcoming development of a real system for breast cancer imaging, which shall feature a reduced number of sources and detectors, compact silicon-based detectors, and that shall avoid exceedingly long breast examination times. These features make the proposed system suitable to be integrated into a compact hand-held probe for breast lesion characterization (malignant vs. benign). To improve the reconstructions, morphological information will be provided by ultrasound imaging, integrating an ultrasound transducer in the probe along with the TD-DOT instrumentation.

Abstract

Time Domain Diffuse Optical Tomography (TD-DOT) enables a full 3D reconstruction of the optical properties of tissue, and could be used for non-invasive and cost-effective in-depth body exploration (e.g., thyroid and breast imaging). Performance quantification is crucial for comparing results coming from different implementations of this technique. A general-purpose simulation platform for TD-DOT clinical systems was developed with a focus on performance assessment through meaningful figures of merit. The platform was employed for assessing the feasibility and characterizing a compact hand-held probe for breast imaging and characterization in reflectance geometry. Important parameters such as hardware gating of the detector, photon count rate and inclusion position were investigated. Results indicate a reduced error (<10%) on the absorption coefficient quantification of a simulated inclusion up to 2-cm depth if a photon count rate ≥ 10⁶ counts per second is used along with a good localization (error < 1 mm down to 25 mm-depth).

1. Introduction

Medical imaging often employs tomography as an effective tool to help detecting and characterizing both pathological and physiological conditions in tissues. Tomography provides 2D/3D images of the tissue under investigation and is performed exploiting approaches that typically involve ionizing radiations or complex magnetic field analysis (e.g., X-Ray CT, MRI).

A different non-invasive experimental approach involves injecting low-power laser light in the Near-Infrared (NIR) spectrum and studying photon paths through the tissue to derive its optical properties, i.e., absorption and scattering. Thus, performing optical measurements with a set of properly arranged sources and detectors enables a full 3D reconstruction of the optical properties that correlate with the nature of imaged tissue. Such an approach is known as Diffuse Optical Tomography (DOT), because light propagation follows a quasi-diffusion like model when working in the NIR on large (cm³) tissue volumes (e.g., for clinical applications). Other approaches involve the detection of the diffused fluorescence signal (Fluorescence DOT) emitted by a specific fluorescent marker after being injected into the tissue [1].

Specifically, Time Domain DOT (TD-DOT) is performed by injecting pulsed laser light and detecting the effects of propagation through the tissue (i.e., amplitude attenuation, delay and temporal broadening of the pulse) to retrieve the optical properties of the tissue [2].

TD-DOT is potentially powerful for applications needing in-depth body exploration (such as brain activation [3], thyroid and breast imaging), and can be implemented in compact hand-held probes for easier non-invasive clinical examinations [4]. The main advantage of the TD approach is that it enables one to uncouple the effects of scattering from the effects of absorption in a single measurement, with no need for multiple distance measurements. Furthermore, when reflectance geometry is used, the photon Time-of-Flight (TOF) provides information on the depth probed by the measurement. More generally, the TD approach, exploiting time as an additional dimension, reduces —at least in principle—the required number of sources and detectors with respect to continuous wave or frequency domain operation, and that enables simpler implementations for medical use [5].

In the last years an important technology push has made the TD approach more attractive: Detectors like Silicon Photomultipliers (SiPMs), adapted from different applications, proved to be cheaper, more compact and sensitive than conventional detectors (e.g., Photomultiplier Tubes) and started to be widely applied as single-photon detectors in optical measurements [4,6]. At the same time, SiPMs can be developed for hardware gating of the photon distribution, thus opening to higher penetration depths, but still keeping small dimension, which is not trivial for gated detectors [5].

Despite recent technology innovations, there is plenty of room for improvement in TD-DOT in reflectance geometry on the software side (e.g., reconstruction algorithms), as well as on the hardware side (e.g., more compact multiwavelength or tunable sources). Thus, anticipating a consistent framework for this modality, a clear and widely shared approach to performance quantification is essential for comparing results from different implementations, driving improvements in hardware and data analysis, and anticipating problems that may arise in the clinical environment.

Breast cancer applications are crucial, as it is the most common cancer amongst women, with about 1 in 8 females being diagnosed with the disease during their lifetime [7,8]. Survival rate significantly improves when an early diagnosis is reached. Nowadays screening mammograms play a key role in the prevention programs, but X-ray mammography still suffers from a low sensitivity that generally leads to unnecessary, social-care expensive, painful and invasive examinations [9,10]. Apart from Magnetic Resonance Imaging (MRI) [11,12] and Ultrasound (US) imaging [13,14], that are already applied in clinics as supplementary imaging tools, several techniques are presently being investigated for the quantitative imaging of breast cancer [15]. Among them, a hand-held optical probe for TD-DOT could prove useful to improve diagnostic sensitivity, performing non-invasive, cheap and physiological-related examinations to better characterize and classify breast lesions after a positive screening mammography.

The purpose of the paper is twofold. First, we present a simulation platform suitable for simulations of large multi-parameter scenarios in DOT, and then we exploit this platform for a systematic study of the performances and feasibility of a TD-DOT hand-held probe to work in reflectance geometry for breast cancer imaging and characterization. More in detail, the key features and the novelty of the platform are: (i) Fully automated scheme which permits to run automatically large multi-parametric scenarios in batch processing; (ii) an objective, quantitative and operator-independent assessment of the outcome of the reconstruction in view of the desired clinical application; (iii) attempt to enforce standardization and performance assessment criteria, following internationally agreed protocols (e.g., BIP [16], MEDPHOT [17], NEUROPT [18]). As compared to existing tomographic packages (e.g., TOAST [19] or NIRFAST [20]) our platform does not offer novelty in the algorithms, but rather aims at studying automatically and with an objective and quantitative outcome the effects of different parameters in a flexible and systematic way.

The second aim of the paper is to exploit the simulation platform to investigate TD-DOT in reflectance geometry with a limited set of source-detector couples towards applications for breast diagnostics. Many parameters with a complex interplay were considered simultaneously, at the level of hardware implementation (number and temporal width of the time gates, counts per second), experimental conditions (depth of the absorption inhomogeneity), and reconstruction parameters (regularization, voxel size). Also, an objective quantification of the tomography reconstruction is needed to assess the capability of the system to classify different lesions. The results will display the limits of applicability, the expected performances, and the key bottlenecks of such a compact 3D probe with respect to realistic scenarios.

The rest of this paper is organized as follows: In Section 2 we present the architecture of the developed multipurpose simulation platform, followed by a detailed discussion about the specific kernel for time domain breast tomography, with problem specifications. Finally, we describe the systematic set of simulations we performed for investigating a wide combination of variables. In Section 3 we report the main results we obtained along with the accuracy in the estimate of optical parameters.

2. Materials and Methods

2.1. Multi Simulation Platform

The (Matlab^®-based) simulation platform was developed to perform systematic simulations mimicking realistic instrumentation for different kinds of applications (simulation kernel), and a generic software architecture was exploited to build a flexible and multipurpose tool. At the same time, we focused on the quantitative assessment of the outcomes and on their comprehensive visualization to help better understand system performances. Quantification is reached by the identification either of some meaningful characteristics or figure of merit (FoM) that well describe system performances, or by comparison with what is defined in well-established standardization protocols. Finally, a visualization section provides a graphical representation of FoMs when up to 5 parameters are changed in the simulations.

Figure 1 depicts the software architecture. The scheme can be divided into 5 main parts:

• Setting of the simulation framework variables;
• Simulation parameters update and override onto the default values;
• The kernel: The single-shot simulation core code for the specific application the user wishes;
• Quantification and extraction of the FoMs;
• Visualization of the FoMs.

The first step is to set the entire variables framework the user wishes to use for the simulation batch and flag only those for which the user wants to study a specific range of values. The space of variables is described in an abstract way, so that the user can decide for each variable either to fix a constant value, or to specify a range of discrete values. For up to five independent variables, the program will study all possible combinations of values in the range.

The number of total iterations of the kernel code is given by the product of the numericity of the flagged variables. Each iteration thus corresponds to a specific combination of the values of the flagged variables. The values of all the variables are set just before the kernel runs.

At the end of each iteration, FoMs are calculated and stored for following visualization.

2.2. Tomographic Kernel

The kernel is divided into two sections: (i) Forward problem and (ii) Inverse problem.

For the forward problem, analytical solutions to the Diffusion Equation under Extrapolated Boundary Conditions for a semi-infinite medium have been computed both for the homogeneous case and for the heterogeneous case. In the latter case, a perturbative solution under the Born approximation has been used over a voxelated volume [21]. Notwithstanding the large availability of more sophisticated models, such as finite-elements (FEM) [22] or Monte-Carlo [23,24], the use of an analytical solution has the advantage to be quite faster than numerical methods, and to be independent of the temporal discretization. After the computation of the forward time-resolved curves, they are converted to photon-counts, accordingly to the desired count-rate, responsivity of the system, measurement time, source/detector areas, and source power. More in detail: the fluence is translated to actual photon counts, considering an average injection power of 1 mW at 800 nm. The acquisition time is set to 1 s, while the detectors are simulated with a 7 mm² active area and a 0.04 quantum efficiency. The simulated values represent a realistic estimate based on actual system setups. No optical losses are considered, since the detector is supposed to be in contact with the medium with full acceptance angle. These figures are used to estimate the maximum number of detected photons. Then, to account for the maximum count rate sustainable for the detection electronics, that number is clipped to an upper limit (typically, this value is 10⁶ counts per second). If desired, an experimentally-derived Instrumental Response Function (IRF) can also be convolved to all the curves for simulating a specific detector. Afterwards, multiplicative Gaussian noise (0.1%) is added to take into account amplitude fluctuations and in general other instabilities of the system. Finally, Poisson noise is added to simulate the shot noise associated with photon detection. The two types of noise simulate the modification of an ideal distribution of time of flight due to the non-ideal experimental setup (Gaussian noise) and photon detection (Poisson noise). More details can be found in [5]. In particular, measurements with a gated detector can be also simulated [25].

Concerning the reconstruction of the absorption coefficient (${\mu }_a$) in the volume, a Jacobian matrix is first computed by solving the following time convolution [21]:

$$J\left({\mathit{r}}_{\mathit{i}},t_k;{\mathit{r}}_{\mathit{s},}{\mathit{r}}_{\mathit{d}}\right)={{\displaystyle \int }}_{{\mathit{t}}_{\mathit{k}}}^{t_k+\mathsf{\Delta }t}G\left({\mathit{r}}_{\mathit{s}},{\mathit{r}}_{\mathit{i}},t\right)⨂G\left({\mathit{r}}_{\mathit{i}},{\mathit{r}}_{\mathit{d}},t\right)⨂h\left(t;{\mathit{r}}_{\mathit{s}},{\mathit{r}}_{\mathit{d}}\right)dt,$$

(1)

where r_i is the position of the i^th voxel, $t_k$ is the k^th $\mathsf{\Delta }t$ broad temporal window, ${\mathit{r}}_{\mathit{s}}$ and ${\mathit{r}}_{\mathit{d}}$ the source/detector position, respectively, $G\left({\mathit{r}}_{\mathit{l}},{\mathit{r}}_{\mathit{m}},t\right)$ is the Green’s Function for the semi-infinite medium calculated at the position ${\mathit{r}}_{\mathit{m}}$ for a source placed at ${\mathit{r}}_{\mathit{l}}$, and $h\left(t;{\mathit{r}}_{\mathit{s}},{\mathit{r}}_{\mathit{d}}\right)$ is the IRF related to the source/detector pairs placed at ${\mathit{r}}_{\mathit{s}}$ and ${\mathit{r}}_{\mathit{d}}$. The first convolution between the Green’s functions has been analytically solved to speed up the calculation [21].

The distribution of the absorption perturbation with respect to the background optical properties is retrieved by solving the following problem:

$$\mathsf{\Delta }{\mu }_a=arg\underset{\mathsf{\Delta }{\mu }_a}{\min }\left\{{{\displaystyle \sum }}_i{\left(\frac{y_i-J_i\mathsf{\Delta }{\mu }_a}{{\sigma }_i}\right)}^2+\tau \Vert \mathsf{\Delta }{\mu }_a\Vert {}^2\right\},$$

(2)

where $\mathsf{\Delta }{\mu }_a$ is the vector of the absorption perturbations in each voxel ($\mathsf{\Delta }{\mu }_a={\mu }_a^{perturbation}-{\mu }_a^{background}$), y_i is the simulated photon-count in each time window and for every source/detector pairs, J_i is the i^th row of the Jacobian, ${\sigma }_i=\sqrt{y_i}$ is the standard deviation accordingly to the Poisson statistics, and $\tau$ is a regularization parameter. Even though the chosen regularization is basic, it provides a solid basis for the future usage of more sophisticated methods, including total variation or ${\mathrm{l}}_{\mathrm{p}}$ sparsity [26,27,28,29]. Finally, it is worth noticing that the whole simulated distribution time of flight is used for the reconstruction, thus exploiting the information carried by the time domain data.

Problem (2) is solved by a direct solution of the normal equation:

$$\left({\tilde{J}}^T\tilde{J}+\tau I\right)\mathsf{\Delta }{\mu }_a={\tilde{J}}^T\tilde{y},$$

(3)

where $\tilde{J}$ and $\tilde{y}$ are the Jacobian matrix and the simulated data vector normalized by standard deviation accordingly to the following:

$$\{\begin{array}{c}\tilde{J}={\mathsf{\Sigma }}^{-1}J\\ \tilde{y}={\mathsf{\Sigma }}^{-1}y\end{array},\mathsf{\Sigma }=diag\left({\sigma }_i\right),$$

(4)

2.3. Simulation Overview

For our Time Domain Diffuse Optical Tomography (TD-DOT) simulations in reflectance geometry we used a fixed geometry of sources and detectors placed in a rectangular pattern (30 mm × 20 mm). On each long side of the rectangle 4 pairs of sources and detectors are placed, with overlapping source and detector positions. Spacing between two nearest source-detector pairs is of 10 mm. The simulated probed volume is a parallelepiped (60 mm side and 30 mm) and the rectangular pattern of sources and detectors is located on its top surface. A localized spherical heterogeneity (~5 mm radius, ${\mu }_a$ = 0.04 mm⁻¹) is moved within the simulated volume (${\mu }_a$ = 0.01 mm⁻¹), thus probing the system performances in the region underlying the sources and detectors arrangement. Figure 2 depicts the phantom and measure geometry.

Both sources and detectors are simulated using realistic signal and noise (e.g., source power, Instrument Response Function, detector area and quantum efficiency, and signal noise), mimicking the typical performances of a TD instrument. We modeled the breast as a semi-infinite medium, with a local absorbing perturbation with respect to the background.

To provide a full characterization of our system we performed a thorough series of simulations studying the main parameters that could affect system performances both on the hardware side (e.g., temporal gating) and the software side (e.g., reconstruction parameters). The set of simulations is summarized in

Table 1. Set of performed simulations.

	Parameter #
Simulation Name	1	2	3	4	5	Added Poisson Noise
A	Xp	Yp	Zp	Regularization	Voxel Size	No
B	Xp	Yp	Zp	Time Channel	Number of Temporal Steps	No
C	Xp	Zp	Number of Time Windows	Regularization	-	Yes
D	Xp	Zp	Number of Delays	Regularization	-	Yes
E	Xp	Zp	Counts per Second	Regularization	-	Yes
F	Xp	Yp	Zp	Regularization	-	Yes

Where: Xp, Yp, Zp are the coordinates of the perturbation within the volume; “Regularization” refers to the regularization parameters (named $\tau$) indicated as a fraction of the maximum singular-value of the Jacobian $\tilde{J}$ in Equation (3); “Number of Time Windows” is the number of software gating windows within the temporal scale; “Number of Delays” is the number of hardware gating delays of the detectors; “Counts per Second” is the photon count rate of the emitting source; “Time Channel” is the temporal width of the single channel of the temporal axis; “Number of Temporal Steps” is the total number of channels of the temporal axis.

.

Table 2. Ranges and default values of the parameters under investigation.

Parameter Name	Range	Default Value	Unit
Xp	−20:5:20	-	mm
Yp	−20:5:20	0 (C, D, E)	mm
Zp	5:5:25	.	mm
Counts per Second	105, 106, 107, 108	106	cps
Number of Delays	1, 3, 5, 10	1 (A, B), 3 (C, D, E, F)	-
Regularization	10−4, 10−3, 10−2, 10−1, 1	10−4 (B)	-
Voxel size (1 side of a cube)	1, 2, 5	2	mm
Time Channel	6, 12, 24, 48	12	ps
Number of Temporal Steps	200, 400, 600	600	-
Number of Time Windows	5, 10, 20, 40	20	-

lists the values of the parameters that have been studied.

Quantification of the simulations outcome is reached through the definition of FoMs:

• The error on the localization of the reconstructed inclusion along the three axes x, y and z with respect to the simulated position according to the center of mass.

  $${ϵ}_{loc}={{\displaystyle \int }}_{ROI}^{}{r}_i^{rec}-{{\displaystyle \int }}_{ROI}^{}{r}_i^{true}i=x,y,z,$$

  (5)
• The volumetric error on $\mathsf{\Delta }{\mu }_a$ is calculated as:

  $${ϵ}_{vol}=\frac{{{\displaystyle \int }}_{ROI}^{}\mathsf{\Delta }{\mu }_a^{rec}dV-{{\displaystyle \int }}_{ROI}^{}\mathsf{\Delta }{\mu }_a^{true}dV}{{{\displaystyle \int }}_{ROI}^{}\mathsf{\Delta }{\mu }_a^{true}dV},$$

  (6)

  where ROI is a spherical region of interest (ROI) of radius 2 cm cantered on the center of mass of the reconstructed inclusion and $\mathsf{\Delta }{\mu }_a^{true}$ and $\mathsf{\Delta }{\mu }_a^{rec}$ are the true and reconstructed absorption perturbation values.

Each simulation set is made of ~2000 independent reconstructions and requested an average of 12 h to complete on a dedicated work station (Dell Precision Tower Workstation 7910, mounting Dual Intel Xeon processors—10 cores, 2.35 GHz—and 64 Gb of RAM memory).

3. Results and Discussion

For the sake of clarity, in the following we will focus on the most meaningful simulations (simulations E and F) and report the main conclusions for the simulations A, B, C and D. For simulations from A to D, we provide external references to FigShare links, which give access to the relative figures that are not reported in this paper, so as not to overload the discussion [30,31,32,33,34,35].

Simulation A has provided indications as for the best voxel size to use during the reconstruction algorithm showing how the 1 mm and 2 mm voxel size are totally equivalent, while the 5 mm voxel size slightly worsens the localization errors and spatial resolution in favor, though, of a significant reduction of calculation time [30,32,35]. We used the 2 mm voxel size for the rest of simulations.

Simulation B showed no significant changes with the Time Channel duration up to 48 ps. A default 12 ps time channel duration was chosen. The number of temporal steps was fixed to be 600, thus leading to a ~7 ns long temporal scale, which fully includes the dynamic range of the simulated curves [33].

Simulation C suggested that the optimal temporal width of the software gating windows should be less than 500 ps. A wider time window introduces errors especially on the absorption coefficient estimation and reduces the penetration depth. The time window duration is given by the duration of the time scale (i.e., the number of time channels multiplied by the time channel duration) divided by the number of time windows. As we fixed the time channel duration to 12 ps and the number of temporal steps to 600, the minimum number of temporal windows is 14: A rounded value of 20 was chosen [31].

Simulation D showed how hardware gating of the detector improves the penetration depth. The relative volumetric error on the absorption coefficient, without gating, is up to 30% at 2 cm depth, while it decreases as close as 10% when several gating windows are used. We identified three as the optimal number of hardware gates: A greater number possibly reduces the photon counts in each window, which may increase error and worsen localization [34].

The previous simulations helped us to define the best variables set to provide a lower absorption quantification error and a better localization of the perturbation. Simulation E is of great interest as it studies the effects of the signal level, which is often a concern in optical breast imaging due to high absorption (e.g., at long wavelengths due to water concentration) and breast thickness variability. In the following we will discuss in detail the figures of merit.

Each panel of Figure 3 shows the error of the reconstructed x coordinate of the inclusion with respect to the simulated one, i.e., the localization on x, expressed in millimeters. The horizontal axis of each panel corresponds to the simulated x coordinate of the inclusion (from −20 mm to 20 mm at 5 mm step), while the Y axis is the depth of the inclusion, i.e., the z coordinate (from 5 mm to 25 mm in 5 mm steps). Different rows correspond to different count rates (10⁵, 10⁶, 10⁷, 10⁸ from top to bottom) while columns represent different values of the regularization parameter (10⁻⁴, 10⁻³, 10⁻², 10⁻¹, 10⁻⁰ from left to right). The error is color-coded: Red colors correspond to positive errors, while blue colors are negative errors; small errors are identified by greenish colors. The thick black line is the error-free contour line, while a label indicates the error value on the other contour lines.

Figure 3 suggests a general high localization of the x coordinate in most of the conditions as the green color (indicating low error, close to 0 mm) is spread over most of the panels except for low regularization (τ = 10⁻⁴). This is expected, as a very small value of the regularization parameter leaves excessive room to the effects of the reconstruction noise. At the same time, the localization error is spatially more uniform for higher regularization: At a fixed count rate, the error dispersion tends to reduce as the contour lines are straighter and zones with different error value are merged together.

However, values exceeding 10⁻¹ could worsen the figure of merit, indicating that a high value of regularization might induce artifacts. The count rate value does not seem to affect the x coordinate localization error for τ ≥ 10⁻³, as within each column a low error (close to 0 mm) is obtained for every experimented condition. A good localization is achieved even in depth down to 25 mm for τ ≥ 10⁻³.

The same considerations apply to the localization along the y coordinate (data not shown).

Figure 4 shows the localization on the z coordinate, i.e., the reconstructed z coordinate of the inclusion with respect to the simulated one. Similar conclusion to Figure 3 can be drawn for Figure 4, but it is worth noticing the low localization error (around 0 mm) even in depth (25 mm) for τ = 10⁻³, which seems to be the best value for the regularization parameter. This is a key aspect as the maximum probed depth is a typical issue when dealing with DO systems, and 25 mm is considered to be a significant value. Moreover, higher regularization values cause the inclusion to be reconstructed at shallower depths, as the negative (blue) error increases for a fixed z coordinate value (same row). As an example, let us consider the last row (count rate = 10⁸ cps) and Zp = 15 mm. The error increases (in absolute value) from 0 mm at τ = 10⁻³ to −2 mm at τ = 10⁻¹, and reaches −3 mm at τ = 1.

Figure 5 shows the relative volumetric error on the absorption coefficient. The figure clearly highlights that high regularization values cause a negative error—i.e., an underestimation of the absorption coefficient—rapidly increasing with the regularization parameter (i.e., left to right within each row). Again, the best value of $\tau$ seems to be 10⁻³ as the volumetric error is kept as low as 20% (in absolute value), even at 2 cm-depth. Clearly, the error on the quantitation of $\mathsf{\Delta }{\mu }_a$ can be related to the underestimation of z. The regularization tends to attract the perturbation towards the surface and consequently a lower $\mathsf{\Delta }{\mu }_a$ is enough to match the observed contrast. For $\tau$ = 10⁻³, a photon count rate ≥ 10⁶ leads to a wide low-error region (<10%) up to 2 cm depth, which is a positive result as ~10⁶ counts per second (cps) is a realistic signal level.

Figure 6 sums up the previous results and corresponds to simulation F. In this case the Y axis of each subplot is the y coordinate, while rows correspond to a different z coordinate of the inclusion. Count rate is now fixed to 10⁶ cps. A low error region is located within the rectangle identified by sources and detectors positions, and a 2-cm depth can be reached with a ~80% accuracy level in the estimate of the absorption coefficient if a $\tau$ = 10⁻³ is chosen. Higher regularization values cause underestimation of the absorption coefficient as well as worse localization.

Figure 7 is obtained in the same condition as Figure 6, but it shows the Contrast-to-Noise Ratio defined as:

$$CNR=\frac{{\overline{\mathsf{\Delta }\mu }}_{rec,roi}-{\overline{\mathsf{\Delta }\mu }}_{rec,\overline{roi}}}{\sqrt{w_{roi}{\sigma }_{roi}^2+w_{\overline{roi}}{\sigma }_{\overline{roi}}^2}}$$

where ${\overline{\mathsf{\Delta }\mu }}_{rec,roi}$ and ${\sigma }_{roi}$ are the mean value and standard deviation of the reconstructed $\mathsf{\Delta }{\mu }_a^{rec}$ calculated in a region of interest; ${\overline{\mathsf{\Delta }\mu }}_{rec,\overline{roi}}$ and ${\sigma }_{\overline{roi}}$ are calculated on the complementary region of interest ($\overline{roi}=1-roi$); $w_{roi}$ and $w_{\overline{roi}}$ are weights defined as the volume rations $w_{roi}=\frac{V_{roi}}{V_{tot}}$ and $w_{\overline{roi}}=\frac{V_{\overline{roi}}}{V_{tot}}$. The region of interest is defined as the region of space in which the reconstructed $\mathsf{\Delta }{\mu }_a^{rec}$ is greater than half of the maximum reconstructed value.

This figure of merit is an indicator of the goodness of the reconstruction for specific structures within the volume [36]. The CNR strongly decreases with the depth for any regularization value, but quicker for higher $\tau$ values. The column $\tau ={10}^{-3}$ and $\tau ={10}^{-2}$ provide the best trend for the CNR, leading to a CNR close to 10 at a depth of 2 cm, even though $\tau ={10}^{-2}$ seems to offer a more homogeneous distribution at all depths.

Overall, Figure 6 and Figure 7 suggest that the simulated system would be capable of detecting inclusions up to 2 cm depth with a significant accuracy and reasonable contrast.

The previous figures show the spatial distribution of the figures of merit in different experimental conditions (i.e., photon count rate, regularization parameter). However, the most common method of assessing the goodness of a reconstruction is through the cross sections of the overall tomographic image. Thus, in the following Figure 8 and Figure 9, we show an example of the simulated tomography problem and of the obtained reconstruction, respectively. Each subplot corresponds to a different depth (expressed in millimeters) and the X/Y axis is the x/y coordinate in space. The color bar shows the value of the absorption coefficient in mm⁻¹ units. The count rate is fixed at 10⁶ cps and three hardware gates of the detectors are simulated. The regularization parameter is set to 10⁻². We simulated the spherical inclusion (whose features are described in Section 2.3) in the center of the x-y plane and at a 10 mm-depth, thus spanning from 5 mm to 15 mm.

The reconstruction in Figure 9 shows small surface artifacts, but the localization in the x-y plane is good as the inclusion is set in the center of such plane, as expected. The z coordinate is reconstructed at a slightly shallower position, as the center of the inclusion is at 9 mm, but the overall height is properly identified. The quantification error is limited to 5%.

All of the presented simulations were carried out at 1 s acquisition time per source-detector pair over the number of hardware gates of the detector. Overall, considering the 64 combinations of the source-detector pair, the measurement time would be around 1 minute. That measurement time was estimated as an acceptable trade-off between the collected signal and measurement duration. Furthermore, different source-detector implementations may allow shorter times. For example, efficient choosing of the most meaningful source-detector pair could reduce the measurement time, at least by a factor of 2.

4. Conclusions

In our work on breast imaging and lesion characterization, we wanted to assess the performances of a Time Domain Diffuse Optical Tomography (TD-DOT) system working in reflectance geometry, formed by eight couples of sources and detectors placed in a rectangular pattern, with a particular focus on the quantification of the system figures of merit. Thus, we developed a general-purpose simulation platform to carry out a quantitative performance assessment of clinical TD-DOT systems. The platform is able to study up to five parameters at the same time and to extract the corresponding figures of merit. It also provides an immediate and concise visualization tool for an analysis and comparison of results. We used the platform to provide a thorough characterization and validation of our TD-DOT system using realistic levels of signal and noise and studying important parameters, such us inclusion position, hardware gating of the detector, regularization parameter ($\tau$) and the total photon count rate. The main figures of merit used are the localization error of the reconstructed inclusion with respect to the simulated one and the relative volumetric error on the absorption coefficient.

We found that the hardware gating of the detector significantly improves absorption quantification, enabling one to probe the tissue down to 2 cm-depth with a volumetric error as low as 10%. At the same time a software gating of the reconstructed photon distribution of time of flights improves the penetration depth if the time window is shorter than 500 ps.

Better inclusion localization in the x-y plane (error < 1 mm) was reached under most of the experimented conditions when $\tau$ ≥ 10⁻³ and down to 25 mm-depth in the source-detector rectangular region. The localization error is low along the z axis for $\tau$ = 10⁻³ (error < 1 mm down to 25 mm-depth) but increases with higher values of $\tau$, reaching a −4 mm value when $\tau$ = 1. Thus, the inclusion is reconstructed at a shallower depth and there is a strong underestimation of the absorption coefficient. Finally, once we fixed the other parameters, a photon count rate ≥ 10⁶ cps leads to a wide volume (under the source-detector pattern and down to a 2 cm-depth) in which the relative volumetric error is low (<10%), indicating a high accuracy of the proposed TD-DOT system and fostering its use for breast imaging.

To further improve these results multiple options are being studied: The use of a priori information coming from US images should further improve the localization and penetration depth, and a finite element fitting procedure should improve the quantification of the absorption coefficient as inversion algorithms are avoided for the reconstruction in favor of a more robust data fitting procedure, using an analytical solution of the problem. These approaches could also benefit from the implementation of a discretization error model, thus further improving its quantification capability [37]. The use of a multi-wavelengths approach is being considered to improve quantification. A spectral fitting procedure would provide the concentrations of the main components of the breast which results in a more direct method for breast and lesion type identification.