1. Introduction

Diffuse Optical Tomography (DOT) is an emerging noninvasive imaging technology that can reconstruct either 2D or 3D images (optical properties distribution) of highly diffusive (turbid) media such as bio-tissues [1–3]. DOT exploits the particle nature of light (photons) in the biological window (600 nm - 1000 nm where light scattering dominates absorption) to explore the chemical composition of tissues of up to a few cm in depth [1,4]. Conversely, other optical tomography modalities such as optical coherence tomography (OCT) are utilized to explore shallow regions (1 - 3 mm depth) of the turbid media by exploiting the wave nature of the light [1,5]. Therefore, in DOT, light loses its wave like nature and photons diffuse within turbid media thicker than 1 cm according to the Radiative Transfer Equation (RTE) [1,2,6]. Because of the complexity of the RTE, its simplified approximation, the diffusion equation (DE), is more commonly used to simulate light propagation (forward problem solver) and to reconstruct tomographic images by recovering the optical properties (OP) distributions (inverse problem solver) of turbid media [7–9]. In turbid media, the OP can be represented by several parameters, with the main ones being: the reduced scattering coefficient (${\text{μ}}_{\text{s}}^{\text{'}}$), the absorption coefficient (µ_a) and refractive index (n) [5,10].

However, this inverse problem is ill-conditioned, ill-posed and suffers from nonlinearity [10–13]. Prior knowledge of the structure (particularly hard prior) plays a significant role in improving the reconstructed images of the turbid media because it limits the number of unknowns when estimating the OP. Therefore, using prior knowledge leads to a reduction in the under-determination of the inverse problem [11,12]. Prior knowledge is classified to two types (based on the method of the minimization process) such as hard prior and soft prior [1,11,12]. In the hard prior approach, the turbid media can be segmented (using meshes or voxels) into several homogeneous regions where all nodes within a region have the same OP [1,3,12]. Therefore, the number of unknowns is reduced drastically to be the number of regions in the turbid medium instead of the total number of nodes. On the other hand, prior knowledge is called soft prior when the number of unknown nodes are preserved, whereas the structural information is incorporated into the regularization term [1,11,12]. Moreover, the anatomical information of any organ can be extracted from another imaging modality such as Magnetic Resonance Imaging (MRI) or Computed Tomography (CT) [12,14]. In DOT, the Time-Resolved (TR) approach is known for its superiority over other approaches of DOT such as Continuous Wave (CW) and Frequency Domain (FD). This is because the TR approach can be used to accurately quantify the values of the OP (${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a) separately, and to reconstruct the highest quality images of the OP distribution within the turbid medium [1,15]. However, the inverse problem solvers suffer from costly computations of the large amount of raw data that is generated by TR-DOT prototypes [1]. Therefore, during the last few decades the issue of having a huge amount of data when using Time-Resolved Diffuse Optical Tomography (TR-DOT) systems was one of the principal challenges that has hindered its expansion and usage. Recently, several methods and algorithms have been developed by researchers in the field to speed-up the data processing. One of these approaches is to reduce the amount of data that is processed by the inverse problem solver and improve the image quality by focusing on the most useful points (seven well-defined points) of the Distribution Time-of-Flight (DToF) of the reemitted photons [16]. Thereafter, we developed an inverse problem solver (utilizing seven well-defined points of DToF histograms) to reconstruct high quality images for several shapes of turbid media (simulation based), which increases the algorithm’s computational efficiency [17].

In this article, a computationally effective TR-DOT prototype was built using the accelerated inverse problem solver to reconstruct tomographic images of solid cylindrical phantoms [17]. The performance of this integrated TR-DOT prototype is investigated using the accelerated inverse problem solver to analyze the experimental raw data (DToF histograms) generated by the Time-Resolved Diffuse Optical Spectroscopy (TR-DOS) setup. The performance of the prototype is evaluated based the computational load and time (speed of image reconstruction process) along with several criteria which are important for functional imaging such as: accuracy of the OP quantification of inclusions and phantoms, inclusions detection, and localization. The paper is organized as follow. In section two, the components of the TR-DOT prototype such as light sources, photon timing subsystem, and the inverse problem solver are discussed. In section three, preparation of phantom and inclusion materials, measurement of the total Instrument Response Function (IRF_Total), and procedures for recording the DToF histograms are described. In section four, the results of the system (3D images and recovered OP of the media) are discussed and the performance of this prototype is evaluated according to different conditions for the experiments. The main contribution of this paper is to validate a computationally efficient TR-DOT which ensures fast reconstruction of high quality images using a relatively small amount of measured data from a reasonable number of source-detector pairs. This work can lead to the development of real-time TR-DOT systems in the near future.

2. Prototype architecture

The TR-DOT prototype has two main parts: TR-DOS setup and software for image reconstruction process.

2.1 TR-DOS setup

A pulsed diode laser (685 nm) is used for illumination with a maximum average power of ~10 mW [18]. Light is transmitted from the laser source to a phantom through a multi-mode optical fiber [19]. The reemitted photons from the phantom are collected by an optical fiber (Edmund Optics Inc., numerical aperture = 0.22 and diameter = 1 mm) and counted by a state of art single-photon avalanche diode (SPAD) detector (PDM series, MPD Italy) with a pixel that has a diameter of 50 µm [20,21]. The laser driver and SPAD are connected to a Time-Correlated Single Photon Counting (TCSPC) module (PicoHarp 300, Picoquant) to measure the IRF_Total of the TR-DOS setup, and to construct histograms of DToF for all source-detector pairs [22]. Figure 1 shows the main components of the TR-DOT prototype.

 

Fig. 1 Diagram of the Time-Resolved Diffuse Optical Tomography (TR-DOT) prototype.

Download Full Size | PDF

In the TR-DOS setup, a wheel is installed horizontally to adjust the positions of the detectors, while source positions are adjusted using a bipolar stepper motor to rotate the phantom [23]. The motor makes precise movements and has 200 steps per rotation (1.8 degrees per step). In order to control this motor effectively, it is paired with a motor driver and an Arduino Uno microcontroller. With this setup, it is possible to perform precise measurements to emulate a multichannel system and accurately replicate experiments.

2.2 Image reconstruction process

Light propagation in a turbid medium can be described using the time-dependent diffusion equation (TD-DE) with a boundary condition such as the Robin boundary condition at a particular wavelength [17,24]. The TD-DE and the Robin boundary condition are given by Eq. (1) and Eq. (2) respectively.

In this article, DOT images are reconstructed using an algorithm that was described and published in [17,25]. The algorithm discretizes Eq. (1) for each temporal bin in the histograms of DToF. The finite element method is then used to solve the discretized Eq. (1), for each mesh node within the turbid media. This algorithm analyzes only three points from the rising edge and four points from the falling slope of the DToF histogram for each source-detector pair, and these seven points are stored in the index matrix. These points are specified as a fraction of the maximum number of counted photons of each DToF histogram, which are 20%, 50%, and 80% on the rising edge and 80%, 50%, 20%, and 10% on the falling slope [16,17]. The Jacobians are calculated by a recursive direct method from these seven points to recover the OP of each node inside the turbid medium and to detect and localize any heterogeneity [17]. This algorithm assumes equal ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a for all nodes in the phantom in the first iteration. In the following iterations, the values of ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a are only updated for nodes within the region of interest (ROI) while the OP remain constant in the rest of the medium. During each iteration, the Jacobian matrix is calculated for all seven selected points (using a direct matrix multiplication) for all source-detector pairs which significantly reduces the computation time. The algorithm reconstructs images by recovering the OP of each mesh node within the ROI that represents a specific part of the cylindrical phantom, and source-detector plane is located in the center around ROI. This iterative inverse problem algorithm ends when either the results of the forward model (DToFs) using the assumed OP match the measured results to a certain degree (by solving the minimization problem) or when the specified maximum number of iterations is reached Fig. 4.

The minimization problem is ill-posed and requires a regularization scheme to get a convergent and stable solution. The regularization scheme used, as described in [17], is based on solving the minimization problem in a normalized space wherein the sensitivity of the Jacobian matrix is made equal for all nodes in the ROI. Therefore, at each iteration, the minimization function uses different dynamic ranges of the lower and upper limits of the optical properties at each node based on the corresponding magnitudes of the Jacobian matrix. So the minimization function, to fit the light fluence data, gives priority for changing the optical properties for deep regions where the Jacobian magnitudes are low, more than that given to the boundary regions where the Jacobian magnitudes are high. This regularization scheme does not require adding an additional penalty term to the objective function for convergence [17]. The reconstructed OP in the original space are obtained by scaling back the reconstructed OP in the normalized space as shown in [17]. In this work, the iterative solution of the minimization problem continues until the change in the maximum error of the simulated fluence relative to the measured fluence goes below 5% or the predefined maximum number of iterations are exceeded. More details about the image reconstruction algorithm that is used in this work can be found in [17,25].

3. Measurements

Measurements were taken using an identical arrangement for the source-detector pairs. Each source was monitored with seven detectors located at 90°, 120°, 150°, 180°, 210°, 240°, and 270° degrees (with respect to the source) around a cylindrical phantom. We investigated the required number of source-detector pairs (in simulation) to determine the best possible arrangements to reduce the time in the image reconstruction process and maintain high quality of the images. It was found that increasing the number of sources to more than 12 did not have a significant effect on the quality of images, probably due to the small size of the phantoms. This configuration was used at multiple positions around the phantom (scanning) to emulate a multichannel prototype. For precise positioning, a computer controlled step motor was used to rotate the phantom to adjust the positions of the injected laser pulse according to the number of views. This configuration (hard prior, 15 mm ROI, 20 is the maximum number of iterations, 12 sources with 84 detectors at the same height around the phantoms) was used for all measurements (called standard configuration) except when the impact of fewer sources or detectors, the impact of prior knowledge, and a larger ROI are investigated.

3.1 Phantoms and inclusions

Phantoms were prepared using epoxy resin with hardener, titanium dioxide (TiO₂), and India ink for the phantom matrix media, scattering agent and absorbing agent, respectively. A ratio of weight (10/3) is used for mixing Araldite® Epoxy Resin with Aradur® hardener produced by Huntsman Advanced Materials Inc., USA [26]. To adjust ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a, TiO₂ from Sigma-Aldrich Co. LLC, USA and India ink were used, respectively [27]. On the other hand, inclusion materials were prepared using mixtures of matrix media (gelatin), scattering agent (TiO₂), and absorption agent (India ink). The OP of each mixture of these materials were determined using a steady-state spatially resolved diffuse reflectance system Fig. 2. In this system, a Quartz Tungsten Halogen (QTH) lamp was used to illuminate the phantom and inclusion materials. The reemitted light was detected at different distances from the source (see the probe in Fig. 2). Then, an inverse problem solver was used to recover the OP of the phantom and inclusion materials over a broad range of wavelengths.

 

Fig. 2 Diagram of steady-state spatially resolved diffuse reflectance system used to accurately determine the optical properties (OP) of phantom and inclusion materials.

Download Full Size | PDF

In this paper, two solid cylindrical phantoms (25 mm diameter each) are used, one with one hole, and the other with two holes. The first phantom has one hole (15 mm diameter) and the second phantom has two holes (5 mm diameter each), and each hole was filled with a mixture of inclusion.

3.2 Instrument response function

For TR-DOT experiments, the IRF_Total of the system should be estimated and deconvolved from the measured DToF histograms to calibrate the raw data. The IRF_Total represents the root of sum of squares of the IRF for each equipment in the setup such as the laser source, the detector, optical fibers, and the TCSPC module using [28–31]:

 

Fig. 3 The total IRF of our TR-DOS setup (~0.3 ns) and the recorded DToF at 180°.

Download Full Size | PDF

It is worth noting that deconvolution is time consuming and an ill-posed problem [34]. Therefore, we convolved the measured IRF (IRF_Total) with the simulated DToF in the inverse problem solver to eliminate the impact of IRF_Total of the setup [36]:

3.3 Data acquisition

Light injected into the phantoms had 0.5 mW optical power measured at the tip of the fiber, and a repetition- rate of 10 MHz. The reemitted photons are counted by a SPAD detector in free running mode [37]. All DToF histograms are recorded using the same period of measurement (30 s to count the reemitted photons for each source-detector pair) for all phantoms with a temporal resolution of 16 ps for the TCSPC module. Also, the laser power is fixed for illuminating two phantoms, and all measurements are taken with no “pile-up” effect because it may significantly degrade the accuracy of the measured data [31,38]. Figure 4 shows a flowchart that describes the process, starting from data acquisition, followed by the iterative inverse problem solver, until reconstruction of tomographic images of the phantom.

 

Fig. 4 Flowchart of the integrated TR-DOT prototype.

Download Full Size | PDF

In this work, we used hard prior to segment the phantoms in all experiments. Therefore, in the first step of image reconstruction process, tetrahedral meshes are generated by Iso2mesh (open-source software) to construct the phantom, inclusions, sources, and detectors [39]. Then, the iterative inverse problem solver starts the comparison with the calibrated raw data (DToF for each source-detector pair) as described earlier (subsection 2.2). This comparison is achieved by solving the minimization problem using the function ‘fmincon’ of Matlab. In all experiments, the initial guess of ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a were 1.8 mm⁻¹ and 0.015 mm⁻¹, respectively. Also, the lower and upper limits used (in the minimization problem) are 0.25 and 7 mm⁻¹ of ${\text{μ}}_{\text{s}}^{\text{'}}$, and 0.0025 and 0.08 mm⁻¹ of µ_a, respectively. On the other hand, we have compared the results of using different number of points on the recorded DToF histograms. It was determined that using only seven points is enough to maintain the quality of the reconstructed images and speed up of the computational process.

3.4 Quantification accuracy of OP

To reconstruct the OP accurately, the linearity and the accuracy of the prototype should be evaluated and considered [40,41]. For this purpose, we have simulated the impact of changing ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a for inclusion separately, and present the recovered OP for inclusion and background in two different ways. First, we display the linearity of the recovered ${\text{μ}}_{\text{s}}^{\text{'}}$ for inclusion and background versus several values of the actual ${\text{μ}}_{\text{s}}^{\text{'}}$ of the inclusion, and the recovered µ_a for inclusion and background versus variable values of the actual µ_a of the inclusion Fig. 5 (a,b). Second, we show the crosstalk (coupling) between ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a. The crosstalk can be observed from the recovered ${\text{μ}}_{\text{s}}^{\text{'}}$ for inclusion and background versus several values of the actual µ_a of the inclusion, and the recovered µ_a for inclusion and background versus variable values of the actual ${\text{μ}}_{\text{s}}^{\text{'}}$ of the inclusion Fig. 5 (c,d). From these results we notices that the prototype maintain accurate linearity for the recovered OP of the inclusion in Fig. 5 (a,b), and the recovered OP for the background are slightly fluctuated. However, the crosstalk between ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a is negligible with small values of the actual OP, but it increases noticeably when larger values of ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a are used.

 

Fig. 5 Cross correlation of the actual OP of the inclusion versus the recovered OP for inclusion and background of the first phantom (one inclusion). (a) Actual ${\text{μ}}_{\text{s}}^{\text{'}}$ against recovered ${\text{μ}}_{\text{s}}^{\text{'}}$; (b) Actual µ_a against recovered µ_a; (c) Actual ${\text{μ}}_{\text{s}}^{\text{'}}$ against recovered µ_a; (d) Actual µ_a against recovered ${\text{μ}}_{\text{s}}^{\text{'}}$.

Download Full Size | PDF

4. Results and discussion

The TR-DOT prototype is evaluated separately for each change in the configurations of the source-detector pairs, sizes and the OP of inclusions, as well as settings chosen for the inverse problem. First, the accuracy of the OP quantification and the time of image reconstruction for two phantoms are estimated using the standard configuration (described in section 3). The 1^st phantom has one inclusion and the 2^nd phantom has two inclusions. Second, the effect on the image quality and the time of image reconstruction using a fewer number of sources (6) attached with seven detectors for each source is shown for the 2^nd phantom. Third, 12 sources are used, but each source is attached with a fewer number of detectors (5 and 4) to investigate the expected decrease in the image quality and the time of image reconstruction of the 2^nd phantom. Fourth, the impact of prior knowledge of the structure of the 1^st phantom is investigated by comparing reconstructed images with and without hard prior. Fifth, the influence of prior knowledge on the image quality and the time of image reconstruction are studied for the 2^nd phantom using different sizes for the ROI (1.5 and 1.8 cm). Figure 6 displays images of the actual OP recovered by steady-state spatially resolved diffuse reflectance system (Fig. 2) for the phantoms. For the 1^st phantom, the actual values of ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a are 1.6 mm⁻¹ and 0.011 mm⁻¹ for the background, and 4 mm⁻¹ and 0.05 mm⁻¹ for the inclusion, respectively. For the 2^nd phantom, the actual values of ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a are 1.6 mm⁻¹ and 0.011 mm⁻¹ for the background, 6 mm⁻¹ and 0.075 mm⁻¹ for the right inclusion, and 3 mm⁻¹ and 0.038 mm⁻¹ for the left inclusion, respectively.

 

Fig. 6 Cross sectional and sagittal views of images for actual ${\text{μ}}_{\text{s}}^{\text{'}}$ (left) and µ_a (right) coefficients for two phantoms: (row 1) the 1^st phantom, and (row 2) the 2^nd phantom.

Download Full Size | PDF

4.1 Image reconstruction on phantoms

The first step in evaluating the performance of the prototype is to compare the quality of the reconstructed images for both phantoms that have known OP Fig. 6. Accurately determining the OP is important; however high precision is mainly required for µ_a because small variations of µ_a is the main factor associated with hemodynamic response during muscles or brain activities in functional imaging [42–44]. 3D images were reconstructed, and Fig. 7 shows and cross sectional and sagittal views for all phantoms using the standard configuration (described in section 3). The total time for image reconstruction (shown in Fig. 7) was ~480 seconds to complete image reconstruction process of each phantom.

 

Fig. 7 Left column: cross sectional and sagittal views of images. Right column: the recovered ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a for: the 1^st phantom (rows 1 and 2), and the 2^nd phantom (rows 3 and 4).

Download Full Size | PDF

The image reconstruction using this configuration yields good results for both phantoms. The good results can be attributed to having enough information for the image reconstruction process (due to small size of the phantoms and using hard prior). For instance, hard prior is used to segment the 1^st and the 2^nd phantom into two and three homogeneous regions, respectively. Overall, the accuracy of the reconstructed ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a is very good, particulary for µ_a. Note that the recovered ${\text{μ}}_{\text{s}}^{\text{'}}$ of the inclusion material is overestimated by an acceptable range of error (< 8%) for the inclusion in the 1^st phantom, and underestimated for one inclusion for the 2^nd phantom. The overestimation of the recovered ${\text{μ}}_{\text{s}}^{\text{'}}$ of the inclusion (in the 1^st phantom) is probably caused due to the crosstalk between ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a. From Fig. 5 (a,b) the recovered ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a of the inclusion has accurate linearity particulary for the actual OP used for the inclusion in this work (${\text{μ}}_{\text{s}}^{\text{'}}$ = 4 mm⁻¹ and µ_a = 0.05 mm⁻¹). However, Fig. 5 (d) shows that the recovered ${\text{μ}}_{\text{s}}^{\text{'}}$ of the inclusion is overestimated (almost by 13%) for this the actual µ_a value which can be the reason for the overestimation in this expirement. On the other hand, the underestimation of the recovered ${\text{μ}}_{\text{s}}^{\text{'}}$ for the left inclusion (in the 2^nd phantom) is caused by the variation of ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a in the right inclsuion. This anaccurate estimation were noticed also in the simulated results when there is large variation between the OP used for the two inclusions. The experimental results in this sub-section are for the standard configuration for the prototype and will be referenced throughout this work. The results of varying the experimental conditions such as the number of sources and detectors, prior knowledge and ROI size are examined in the following subsections.

4.2 Impact of number of sources

Increasing the number of source-detector pairs plays a significant role in improving the quality of the images because a larger data set of the DToF histograms (denser source-detector pairs collected photons at more positions of the phantom surface) are processed by the inverse problem solver. In this subsection, we investigated the impact of using fewer source-detector pairs by reducing the number of sources, but keeping the number and positions of detectors for each remaining source the same as in the standard configuration (7 detectors per source). Therefore, the variation in image quality and the time of image reconstruction are compared using two measurement configurations with 12 and 6 sources attached to 84 and 42 detectors respectively. The tradeoff are investigated between reducing time of data acquisition and image reconstruction process without losing the acceptable quality of images, and the results are shown in Fig. 8.

 

Fig. 8 Left column: Cross sectional and sagittal views of images. Right column: the recovered ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a for the 2^nd phantom: 12 sources and 84 detectors (rows 1 and 2), and 6 sources and 42 detectors (rows 3 and 4).

Download Full Size | PDF

In Fig. 8, it is shown that reducing the number of sources to 6 yields nearly the same accuracy of the recovered OP of the background (phantom) as the standard configuration in section 4.1 (12 sources and 84 detectors). Also, the accuracy of the recovered ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a of the right inclusions is slightly different versus the actual values as a result of decreasing the number of measured DToF histograms that are used in the inverse problem. These underestimation of the recovered ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a for the right inclusion lead to decrease the variation between the actual and the recovered ${\text{μ}}_{\text{s}}^{\text{'}}$ of the left inclusion due to the crosstalk between the two inclusions. All 3D images were reconstructed, and total time for image reconstruction for the experiments using 12, and 6 sources required 480 s, and 300 s, respectively. Overall, for a phantom of this size, using half the number of sources (compared to the standard configuration) results in accurate images requiring a shorter time for data acquisition and analysis. This is beneficial (for this size of phantoms or tissue) as the data acquisition time and the time of image reconstruction will be reduced and the inclusions are still detected but with slightly lower accuracy of quantification of the OP.

4.3 Impact of number of detectors

Reducing the number of source detector pairs by using fewer detectors per source is expected to reduce the quality of the images and the time of image reconstruction. Here, we investigate these effects by comparing results from the standard configuration (7 detectors per source) with experiments in which 5 and 4 detectors per source are used. The first experiment has the same detector configuration as the standard configuration, except that two detectors are removed (90° and 270°) when 5 detectors are used. In the second experiment, four detectors are used (90°, 150°, 210° and 270°) to investigate the effect of reducing the number of detectors without reducing the span.

From the results in Fig. 9, it can be noticed that decreasing the total number of detectors from 84 to 60 resulted in almost the same accuracy in the OP estimation of the phantom and the inclusions. In the experiment with 5 detectors per source, ${\text{μ}}_{\text{s}}^{\text{'}}$ of the inclusions are slightly overestimated (for left inclusion) while the estimated µ_a remain almost as accurate as in the standard configuration. In contrast to the previous experiment, the accuracy of the reconstructed OP is degraded significantly when only four detectors per source (12 sources – 48 detectors) are used for the span of 90° to 270°. Therefore, using four detectors for this span automatically increases the gap between each two detectors to (60°), which reduced the accuracy of the recovered µ_a noticeably (in this experiment) for background and two inclusions. Consequently, this large overestimation of the estimated µ_a (almost 35%) for background degrades the accuracy of the recovered OP in this experiment. The total time for image reconstruction for the experiments using 84, 60 and 48 detectors was 480 s, 370 s and 310 s, respectively. Overall, by reducing the number of detectors reasonably (without increasing the separation 30° between them), the data acquisition and the time of image reconstruction can be shortened while obtaining accurate results.

 

Fig. 9 Left column: Cross sectional and sagittal views of images. Right column: the recovered ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a for the 2^nd phantom: 12 sources and 84 detectors (rows 1 and 2), 12 sources and 60 detectors (rows 3 and 4), and 12 sources and 48 detectors (rows 5 and 6).

Download Full Size | PDF

4.4 Influence of prior knowledge

Having prior knowledge of a phantom’s structure is required to generate high quality images using any inverse problem solver. Without prior knowledge, the DToF data alone is not enough to estimate the size and localize inclusions and OP accurately [42,45,46]. In this subsection we investigate the impact of the incorporated prior knowledge on the image reconstruction with our prototype. In Fig. 9, the results obtained for the 1^st phantom are compared with hard prior and without the use of prior knowledge.

The results obtained from the reconstruction process without prior knowledge demonstrate poor localization and poor size estimation of the inclusion. They also show that the recovered OP of the phantom and inclusion are not very accurate, and variations between the actual OP and the recovered OP will be higher at other positions in the cross-section. Also, we note that this prototype can detect the existence of one large inclusion in the phantom (Fig. 10, rows 3 and 4) which is very beneficial for situations where prior anatomical knowledge is not available. The poor localization and size estimation results of the inclusion agree with known limitations of the full image reconstruction problem in DOT which leads to blurry images when no prior knowledge is incorporated into the ill-posed inverse problem [47]. On the other hand, absence of prior knowledge increases the number of unknowns (mesh nodes) which dramatically increases the computation time for the image reconstruction [11]. To reconstruct images without prior structural information, the computation time increased to 3300 s, whereas it required only 480 s when hard prior knowledge was incorporated.

 

Fig. 10 Left column: cross sectional and sagittal views of images. Right column: the recovered ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a for the 1^st phantom: with hard prior (rows 1 and 2), and without prior knowledge (rows 3 and 4).

Download Full Size | PDF

4.5 Influence of ROI size

Increasing the size of the ROI is expected to increase the under-determination of the inverse problem and consequently reduce the accuracy of the estimated OP and increase the time it takes to reconstruct images from the experimental measurements. In addition, absence of prior knowledge along with enlarging ROI are expected to decrease the accuracy of the reconstructed images because of the larger amount of unknown nodes that require OP estimation. In this subsection, we investigate the ability of the prototype to detect two small inclusions in the 2^nd phantom using two different sizes of ROI (1.5 and 1.8 cm), and without incorporating hard prior. Figure 11 shows the results using 12 sources and 84 detectors.

 

Fig. 11 Left column: cross sectional and sagittal views of images. Right column: the recovered ${\text{μ}}_{\text{s}}^{\text{'}}$ and µ_a for the 2^nd phantom without prior knowledge: 1.5 cm ROI (rows 1 and 2), and 1.8 cm ROI (rows 3 and 4).

Download Full Size | PDF

These results suffer considerably from poor localization which degrades the quantification of the recovered OP in the cross-section position at the real positions of inclusions. In addition, existence of inclusions can be observed (particularly from the cross sectional images) in the upper half of the cylinder. The right inclusion can be detected without difficulty due to its large OP values. However, the absence of prior knowledge make inclusion’s detection, OP quantification, and localization less accurate, particularly for such small inclusions. Also, several source-detector planes at different heights around the phantom are required if the ROI is enlarged and prior knowledge is not incorporated. Therefore, we increased the ROI reasonably by 20% to 1.8 cm, and the results are similar of the original ROI (1.5 cm), whereas the computation time for image reconstruction increased noticeably. The total reconstruction time was 3200 s, and 3700 s for 1.5 and 1.8 cm ROI sizes, respectively.

5. Conclusion

In this work, a major step towards building real-time TR-DOT prototypes to reconstruct high quality 3D images was demonstrated through experiments. The reported TR-DOT prototype significantly reduces the total amount of time required to generate high quality images of objects using a TR approach, by selecting 7 points on each recorded DToF histogram. This reduces the computational complexity of the inverse problem solver which accelerates the image reconstruction without degrading the accuracy of the OP quantification for phantoms and inclusions.

To evaluate this prototype for diffuse optical tomography applications, several experimental variations were used with the configuration, phantoms, and inverse problem solver settings. We used two phantoms that were cylindrical in shape and had OP similar to typical human tissue. The 1^st phantom has a large inclusion (15 mm diameter) and 2^nd phantom has two small inclusions (5 mm diameter each). In the first experiment, we reconstructed images on both phantoms using the standard configuration (hard prior, 15 mm ROI, 20 is the maximum number of iterations, 12 sources with 84 detectors at the same height around the phantoms), and high quality images were reconstructed in 8 minutes for each phantom. Also, the accuracy of distinguishing the variation of the OP of two small inclusions in the same phantom was tested successfully. In the second and third experiments, we discussed the results when using fewer numbers of sources and detectors on the 2^nd phantom. We compared the accuracy of the recovered OP and the time of image reconstruction process in the 2^nd and 3^rd experiments with the 1^st experiment. Reducing the total number of sources and detectors by 50% (in the 2^nd experiment) from the standard configuration maintained the quality and accuracy of the reconstructed images and the estimated OP, respectively. Also, the time of image reconstruction process is reduced from 8 minutes in the 1^st experiment to 5 minutes in the 2^nd experiment (using 6 sources and 42 detectors). In the 3^rd experiment, reducing the number of detectors per source (when 5 detectors per source are used and the small gap 30° between detectors is maintained) preserved the quality of the reconstructed images and accuracy of the recovered OP particularly of µ_a, when compared to the results of the 1^st experiment, whereas the time of image reconstruction process is decreased from 8 minutes to almost 6 minutes. However, the accuracy of the recovered OP were degraded considerably when 4 detectors per source (larger gap 60° between detectors) are used in the 3^rd experiment, and the time of image reconstruction process is decreased from 8 minutes to almost 5 minutes. In general, the 2^nd and the 3^rd experiments proved that the data acquisition time of the measurements and the time required for image reconstruction can be reduced reasonably if enough source-detector pairs (> 40) are used and the small gap 30° between detectors is maintained without degrading the quality of the reconstructed images. In the 4th experiment, when the hard prior knowledge was not used, image reconstruction process required much longer time (~55 minutes), and poor size estimation and localization of the inclusion were observed. This requires investigation into what alternative methods can be incorporated into our inverse problem solver, to eliminate the poor results that are obtained when prior knowledge is not used. In the 5th experiment, larger ROI (1.8 cm) maintains same level of accuracy in comparison with small ROI (1.5 cm), but the time of image reconstruction process increased from 53 minutes to almost 62 minutes. The 5th experiment demonstrated that source-detector pairs in several heights are required to be utilized with larger turbid media in transmittance geometry when prior knowledge is not used, but the time of image reconstruction will be very long (> 1 hour). Overall, the variation of these conditions also helps to establish the best conditions to make image reconstruction as fast and accurate as possible to use this TR-DOT prototype for diffuse optical tomography applications.