Open Access
Add to BibSonomy
Abstract
Time-domain near-infrared optical tomography (TD NIROT) techniques based on diffuse light were gaining performance over the last years. They are capable of imaging tissue at several centimeters depth and reveal clinically relevant information, such as tissue oxygen saturation. In this work, we present the very first in vivo results of our SPAD camera-based TD NIROT reflectance system with a temporal resolution of ∼116 ps. It provides 2800 time of flight source-detector pairs in a compact probe of only 6 cm in diameter. Additionally, we describe a 3-step reconstruction procedure that enables accurate recovery of structural information and of the optical properties. We demonstrate the system’s performance firstly in reconstructing the 3D-structure of a heterogeneous tissue phantom with tissue-like scattering and absorption properties within a volume of 9 cm diameter and 5 cm thickness. Furthermore, we performed in vivo tomography of an index finger located within a homogeneous scattering medium. We employed a fast sampling rate of 2.5 Hz to detect changes in tissue oxygenation. Tomographic reconstructions were performed in true 3D, and without prior structural information, demonstrating the powerful capabilities of the system. This shows its potential for clinical applications.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Near-infrared optical tomography (NIROT) is a powerful tool that enables imaging biomarkers in tissue, such as the tissue oxygen saturation. This parameter is clinically highly relevant, because it reflects the balance between oxygen consumption and oxygen supply. Additionally, NIROT is capable of determining the presence and the concentration of tissue components such as oxy-, deoxy- and total hemoglobin ([O2Hb], [HHb], [tHb] in µM), water and lipids. Being non-invasive and relatively low cost, this technology will add a new dimension to biomedical imaging in research and clinical application [1–5]. NIROT relies on near-infrared light penetrating deep into tissue and interacting with different tissue components. The tomogram of the inner structure is reconstructed by use of the light that is detected at the surface after it travelled through the tissue. Among the three NIROT modalities, namely continuous wave (CW), frequency domain (FD) and time domain (TD), TD provides the richest information about the tissue. It employs time-of-flight (ToF) detectors and benefits from the ability to acquire the temporal distribution of the photons that have arrived [5].
Common implementations of TD NIROT employ a limited number of sources and detectors (e.g. [6]), because these instruments are based on optical fibers and photomultiplier tubes [7,8] or silicon photomultipliers (SiPM) [9]. The low amount of information also limits the image quality. However, these systems are bulky, and it is not possible to scale the number of detectors up. Several research groups worked on new systems to solve this problem. For example, an imaging system based on single photon avalanche diodes (SPADs) achieved a large number of virtual sources and detectors by scanning an object’s surface and measuring at many positions [10,11]. However, this system has to be applied in free space, which hinders the use in a clinical environment, where ambient light is present. A handheld reflection TD NIROT system based on a time-gated intensified CCD camera helped diminishing the influence of ambient light [12] and showed sensitivity to inclusions placed at various depths. Since it included only one wavelength, it was unable to determine the oxygenation and an accurate 3D image reconstruction was not attempted.
The recent progress in CMOS technology has brought a new detector type into the field: the SPAD array [13,14]. Featured with low noise and high dynamic range, it provides a previously unachievable density of detectors and thus promises to drastically increase the spatial resolution [15] and depth sensitivity [16] of TD NIROT. Although, recent publications underline the progress in NIROT, the results were achieved under specific conditions, i.e. in transmission mode or employing prior knowledge of both depth of the inhomogeneity and optical properties of the background [15]. Such conditions are unachievable in real clinical scenarios where measurements in reflection mode are commonly required and a priori information is limited.
The goal of this work is to present a new TD NIROT system called Pioneer and to report its performance in various tests. The system works in reflection mode. It is based on in-house developed CMOS SPAD 32 × 32 array with high sensitivity to NIR light [14,17]. Here, we also aim at showing the benefits of the CMOS SPAD technology in NIROT, namely the compact design, the high number of detectors and the ToF information. We discuss several complications due to the tremendous amount of information that the sensor generates and provide solutions. We performed various experiments in vitro to assess all stages of tomogram generation, starting from system validation and going all the way to the image reconstruction of the inner structure of an object. We proceeded to in vivo imaging and monitoring tissue oxygenation during an occlusion test. We demonstrated the ability to reconstruct images rapidly and accurately by employing state-of-the-art hardware with GPUs and massively parallelized computation [18]. We believe this work paves the way for clinical application of NIROT in the near future.
2. Method
2.1 The novel TD NIROT system called Pioneer
The schematic of Pioneer and a photo of the system are shown in Fig. 1(a),(c). It employs a custom designed 32 × 32 pixel TD image sensor called Piccolo [17] based on SPADs with dynamic sharing of time-to-digital converters (TDCs) [19]. Piccolo is the key element of Pioneer, since it provides a large number of detectors and hence source-detector pairs. To take fully advantage of these, we implemented a novel concept of NIROT in Pioneer. In this and the following sections, we describe how we harvested the benefits.
Pioneer employs a supercontinuum laser SuperK Extreme EXR-15 with acousto-optic tunable filter (AOTF) (NKT, Denmark). The laser emits ∼5 ps pulses with a repetition rate of 80 MHz at multiple wavelengths selectable freely in the range of 650-900 nm. An optical fiber guides the light to one of 11 source positions, which is selected by a fiber optical switch (Laser Components, Germany). The use of a tunable pulsed laser enables us to achieve two goals: (i) it provides time-correlated light pulses required for TD signal acquisition, and (ii) the multispectral light enters the tissue at exactly the same location, in our case through the same optical fiber, as described later. This enables a wavelength based coupling correction as previously described in Jiang et al. [20]. Additionally, the ability to freely choose a wavelength makes the system flexible in terms of detecting different chromophores.
Optical fibers guide the light from the optical switch to the source ring. The latter is the only element of Pioneer in contact with an object or subject. It has a rigid structure enclosed in soft black biocompatible silicone and forms a ring with an inner diameter of 25 mm (Fig. 1(b)). This aperture defines the field-of-view (FoV) of Piccolo. The source fibers are equally distributed around the FoV in a circle of 45 mm diameter. They are fixed in the source ring and terminated with transparent windows to guide the light to the tissue. Soft biocompatible silicone Silpuran 2400 (Wacker, Germany) is used for both the optical windows and the black ring enclosure to achieve gentle but tight contact with the tissue. The manufacturing process of the source ring is described in [21]. We used multi-mode optical fibers with a core diameter of 62.5 µm and 0.23 NA (Thorlabs, US) with FC/APC connectors.
We employed a NIR lens SY110M with a short focal length of 1.67 mm and a relatively large aperture of F/1.8 (Theia Technologies, US) to collect light from the tissue and to project it onto the Piccolo sensor. This open space arrangement enables us to interrogate a relatively large FoV of 25 mm diameter with this tiny sensor array. Each SPAD pixel covers an area of 1.06 × 1.06 mm2 of object’s surface. Due to the short focal length and wide angle of view, the lens enables a working distance of only ∼6 cm and a distance from an object to the image plane of ∼12 cm. The space between the object and the lens is shielded by a metal tube from ambient light. It is coated inside with Vantablack (Surrey NanoSystems, UK) to eliminate internal reflections. The tube also acts as a mechanical interface between the camera and the source ring. The case is in-house designed and 3D-printed. It contains the SPAD array and all accompanying electronics. A laptop controls this hardware by an in-house developed C++ / Qt software. Data post processing and image reconstruction is performed on a separate workstation (128 GB RAM, GPU Nvidia Gtx Titan Xp).
Pioneer is applicable inside an incubator for preterm infants. We achieved this by compact probe design keeping all bulky components next to the cot on a movable rack. The case of 21 × 12.5 × 6 cm3 contains the Piccolo camera and its electronics. The optical components of the probe, i.e., the lens, the Vantablack tube and the optical source ring form a cylinder with a diameter of 6.5cm and 11.5 cm length. The probe is connected by wires (SS-USB and power supply) and optical fibers, which are flexible enough to pass through the openings of the cot. The whole probe weighs ∼0.85 kg, which is light enough to be handled manually during data acquisition. The probe can be placed on a pillow and held to the baby’s head. It does not require a fixation. Thus, Pioneer is able to perform tomography of the neonatal brain at the bedside.
2.2 Data analysis and image reconstruction
The Pioneer system generates a huge amount of data, i.e., 32 × 32 pixels of 256 time gates per wavelength per source location. For each pixel the ToFs of the arriving photons are assembled to a histogram inside the Piccolo camera. The histograms are transferred to a powerful workstation (Intel Xeon CPU E5-1620 v4, 128 GB RAM, Titan Xp GPU, Ubuntu 16.04) for image reconstruction. Dark photon counts are eliminated by subtracting the median value of the photon counts distributed over all time gates. The instrumental response function is removed by wavelength normalization (WN) [20] or calibration on a homogeneous phantom with known optical properties. The corrected TD data are used for the image reconstruction. Here, we apply the physical model based image reconstruction. It requires an accurate description of light transport in tissue, which is called the forward problem. Here, we employ the diffusion equation (DE), which is valid in biological tissue due to its high reduced scattering coefficient ${\mu _s}^{\prime}$ and low absorption coefficient ${\mu _a}$ in the NIR. The time-dependent DE is expressed as follows:
To recover accurate optical properties we apply a region-prior, i.e., prior information about the regions in the medium. Structural information drastically reduces the ill-posedness of the inverse problem, since only ${\mu _s}^{\prime}$ and ${\mu _a}$ of a few regions instead of thousands of nodes are recovered [28–30]. This information is retrieved from the volumetric reconstruction. The inverse problem is solved with the MATLAB function lsqcurvefit [28] and the forward problem by NIRFAST. Given that Pioneer features a large number of source-detector pairs (∼2800 in our experiments) and temporal data of a reasonable signal-to-noise ratio (SNR), the inverse problem is less ill-posed, and this improves the accuracy of the reconstructed optical properties.
2.3 Validation of Pioneer by phantom experiments
To validate the performance of Pioneer, we performed measurements on a silicone phantom of 120 mm diameter by 60 mm thickness with a cylindrical inclusion of 10 mm diameter and 30 mm length at a depth of 15 mm. We used Silpuran 2420 two-component silicone (Wacker, Germany) to produce this phantom. We added a defined quantity of silicone-compatible white ink ELASTOSIL RAL 9010 (Wacker, Germany) to obtain the specified $\mu _s^{\prime}$. ${\mu _a}$ was adjusted by carbon black powder. Since carbon black has little wavelength dependence in the NIR range, we added another ink, such as ELASTOSIL Heliogruen RAL 6004 (Wacker, Germany), to obtain a wavelength dependent ${\mu _a}$ of the inclusion. The manufacturing of heterogeneous silicone phantoms was described previously [20,31–33].
We applied the Pioneer probe to the surface of the phantom. One of the 11 light sources was switched on at a time and a temporal histogram was acquired by each detector simultaneously. It took ∼25 s to obtain 2 × 108 photons in the FoV. The whole measurement at 11 sources and 3 wavelengths (689 nm, 802 nm and 840 nm) took ∼14 min. Prior to the measurement, the system was calibrated using a reference homogeneous phantom with known optical properties. The optical properties of the phantom were measured with a commercial device Imagent (ISS Inc., Champaign, IL, USA), as describe e.g. here [33]. The $\mu _s^{\prime}$ and ${\mu _a}$ of the phantom and of the inclusion are shown in Table 1. The inclusion had a higher ${\mu _a}$ compared to the bulk, whereas the $\mu _s^{\prime}$ was similar. The term ‘bulk’ refers to the volume of the phantom without the inclusion, i.e., its background.
Table 1. Optical properties of the phantom
The cylindrical phantom with a dimension of 90 mm diameter and depth of 50 mm was modeled by a fine mesh of with an average node distance of 0.9 mm (231’977 nodes) for obtaining the forward results, i.e., to model the light propagation in the object. This mesh represented only the central superficial part and not the whole silicone phantom (Fig. 2). Limiting the simulation process to this region of interest (ROI) facilitates the computation. In a first analysis to obtain regional information, we applied a coarser mesh of 3’784 nodes to perform a no-prior reconstruction of the spatial distribution of the ${\mu _a}$ in the phantom. The optical properties of the bulk were measured and set as initial guess for the whole volume. Since there is little variation of ${\mu _s}^{\prime}$ in tissue we reconstructed the ${\mu _a}$ only. The ${\mu _s}^{\prime}$ was fixed to the bulk ${\mu _s}^{\prime}$ in the whole volume. We set the maximum number of iterations to 15 and the initial regularization parameter was 10.
Fig. 2. Schematic of the mesh used for modelling. Black contour line represents the geometry and relative position of the actual phantom.
The reconstructed 3D image was created and segmented. The root mean squared error (RMSE) and Dice similarity were calculated to quantify the 3D image quality. We denote the ground-truth image as the array ${\boldsymbol{\mu}}_{\boldsymbol{a}}$, where ${\mu _{a}}_{i},\; i{\; } \in \; \textrm{N} = {\; }\{{1,{\; }2,{\; }\ldots ,\; \textrm{n}} \}$ is the value in the i-th mesh node, and the reconstructed image as $\boldsymbol{\hat{{\mu }}}_{\boldsymbol{a}}$ and ${\hat{\mu}}_{{a}_{i}},\; i{\; } \in \; \textrm{N}$. The RMSE is a measure of the error of the reconstruction:
The smaller the RMSE value is, the smaller is the deviation of the reconstructed image from the ground-truth. The Dice similarity quantifies the similarity of the structure of the reconstructed 3D image to the ground-truth [34].
In the second analysis, the aim is to quantitatively obtain the ${\mu _a}$ and ${\mu _s}^{\prime}$ for all the three wavelengths: 689, 802, and 840 nm in the inclusion and bulk. This second analysis is based on the regional information known from the first analysis, i.e., a shape for the inclusion. The phantom provided us with a wavelength dependent ${\mu _a}$ contrast of the inclusion, which varied from 6.7 ${\times} $ compared to the bulk at 689 nm to 2.5 ${\times} $ at 840 nm (Table 1). The initial guess values of ${\mu _a}$ and $\mu _s^{\prime}\; $ for this region-prior reconstruction were ${\mu _a}\; $ = 0.01 mm−1 and $\mu _s^{\prime}\; $ = 0.5 mm−1 for both regions at all wavelengths. ${\mu _a}$ for both the bulk and the inclusion was optimized in the range of [0.001, 0.04] mm−1, and $\mu _s^{\prime}\; $ in the range of [0.3, 1.2] mm−1 [29].
2.4 Imaging hemodynamics of a finger in a diffusive medium
We designed an experiment to prove that Pioneer is able to image internal structures and show hemodynamic information about a diffusive medium. We manufactured a silicone phantom with a cavity with optical properties similar to tissue. The outer shape of this phantom is similar to the phantom described above. The cavity has the shape of an index finger of the left-hand. This finger of a healthy volunteer (male adult, 35-year-old, 183 cm, 85 kg) was placed in the cavity and was tightly surrounded by the soft silicone. We also inserted a small amount of water-Intralipid emulsion with Indian ink into the cavity to exclude air gaps between the silicone and the skin. The optical properties of this liquid matched the ones of the silicone. In this experiment, we used the same probe as described above (Fig. 2), i.e., 11 light-sources in a circle with a diameter of 45 mm arranged around the FoV with a diameter of 25 mm, which is covered by 257 detectors. The Pioneer probe was placed in the center of the flat surface of the cylindrical phantom, with the finger located ∼ 11 mm below the surface (Fig. 3). In addition, we repeated the experiment and measured with the commercial oximeter OxyPrem 1.4 (OxyPrem AG, Switzerland). This study was cleared by the Ethical Committee of the Kanton Zurich, Switzerland.
Fig. 3. Illustration of the finger measurement with Pioneer (a). The Pioneer probe was placed on the top surface of the phantom. The view of the surface shows the locations of sources and the FoV of the detectors. A conventional meter measured the pressure in the cuff on the arm. A similar measurement was performed with the oximeter OxyPrem 1.4 (b). OxyPrem 1.4 features two detectors and two sources with several wavelengths.
2.3.1 Measurement protocol
To elicit controlled hemodynamic changes, we induced venous and arterial occlusions by inflating a cuff that was wrapped around the upper arm of the subject, who was sitting.
Venous occlusion was achieved by rapidly inflating the cuff to 60 mmHg. We performed three repetitions of 2 min of venous occlusion and 2 min of rest. We also measured the baseline conditions for 20 s before and after the series of occlusions. At 60 mmHg, the venous outflow from the tissue is blocked, while the inflow of highly oxygenated arterial blood is unimpeded. Therefore, the [tHb] increases, which reflects the blood flow. Due to the O2 consumed by the tissue, also the [HHb] increases. The ratio of the increase in [O2Hb] divided by the increase in [tHb] reflects the venous O2 saturation [35].
Arterial occlusion: A higher pressure of 250 mmHg was applied to stop the arterial flow. Since an arterial occlusion causes discomfort, the measurements were limited to two repetitions of occlusion for 5 min and recovery for 3 min. Again, the baseline measurements for 20 s were added before and after the occlusion measurements.
For Pioneer, a short exposure time of 400 ms was used per source position and wavelength (725 nm and 802 nm) to capture these rapid changes. These wavelengths were chosen for their sensitivity to hemoglobin and high SNR. The total time for a complete measurement (2 wavelengths, 11 source positions) was 10 s, which includes data acquisition and transfer to the PC. The OxyPrem instrument has a time resolution of 1s.
2.3.2 Image reconstruction for in vivo measurements by Pioneer
Here we employed a novel three-step 3D image reconstruction procedure. It does not require calibration on a homogeneous phantom and enables imaging without prior information about the structure of the object. As a first step, we measured the average optical properties with a long exposure time. As a second step, we performed a no-prior reconstruction of the internal structure of the object. As a third step, we segmented the volume based on this reconstruction and performed a reconstruction of the optical properties of the regions, using only the measured data acquired at short exposure time to obtain a high sample rate.
The first step was to determine the average optical properties. TD data were transformed into FD by a fast Fourier transform. The FD data $\mathrm{\Phi \;\ consist\;\ }$ of the amplitude, i.e. $\ln ({I \cdot {r^2}} )$, where I is the amplitude of the light intensity at 100 MHz and r the scalar source detector distance, and the phase $\theta \; $[20]. The average ${\mu _a}$ and ${\mu _s}\mathrm{^{\prime}}$ were calculated by the multi-distance approach as described in detail in [36]. In short, it is based on the solution to the DE for the semi-infinite boundary condition under the assumption that the coupling from the tissue to the active area of the camera is the same for all pixels [37]. This theory shows that $\textrm{ln}({I\cdot \; {r^2}} )\propto r$ and $\theta \propto r{\; }$ are both linear functions and that the ${\mu _a}$ and ${\mu _s}^{\prime}$ can be calculated from the slopes of these lines. We determined these slopes in our data by a linear fit (MATLAB function fit). We calculated the average slope for 11 sources. This enabled us to determine the average ${\mu _a}$ and ${\mu _s}^{\prime}$ constituting our initial guess for the 3D image reconstruction.
As second step a volumetric image was reconstructed. We accumulated photon counts for an exposure time of 25 s for each source and wavelength. We employed these data for reconstructing the internal structure. For this reconstruction, two cylindrical meshes with dimensions of 90 mm diameter and 40 mm thickness were created in NIRFAST. The fine mesh consisting of 185608 nodes was used for the forward calculation, and the coarser mesh consisting of 2274 nodes for the inverse problem. We applied the previously proposed wavelength-normalization (WN) approach to avoid a tedious calibration for the system response functions [20]. We reconstructed the image from the ratio of the FD data at two different wavelengths. The Tikhonov regularization with an initial parameter of $\beta = 10$ was applied for iterative optimization. The ${\mu _s}\mathrm{^{\prime}}$ was assumed constant in the whole volume and only the distribution of ${\mu _a}$ was recovered to obtain the shape of the finger in the phantom. The reconstructed volume was segmented by setting a threshold for ${\mu _a}$ of twice the median ${\mu _a}$. The segmented mesh was utilized as region-prior for reconstructing ${\mu _a}$ and ${\mu _s}\mathrm{^{\prime}}$ of the finger, and to access the hemodynamic in vivo.
As third step, the optical properties ${\boldsymbol{\mu}}_{\boldsymbol{a}}$ and ${{\boldsymbol{\mu}}_{\boldsymbol{s}}}^{\prime}$ in the finger region were recovered using TD data with a much shorter acquisition time of only 400 ms. The image reconstruction was based on region-priors, i.e., the anatomic information recovered in the second step. This enables the detection of rapid changes in the oxygenation of the finger.
2.3.3 Evaluation of the results of image reconstruction
For the evaluation of the 3D image reconstruction we only evaluated the structural information with Dice similarity, since the ground truth optical properties of the finger are unknown. We converted changes in the ${\mu _a}$ into changes in [O2Hb], [HHb] and [tHb] in the finger region using the absorption spectrum O2Hb and HHb. We calculated the correlation coefficient between the changes in [O2Hb] and [HHb], and the inflation status of the cuff (1 for inflated and 0 for rest) for both the venous and arterial occlusion tests.
Additionally, we obtained the changes in [O2Hb], [HHb] and [tHb] from the output of the OxyPrem 1.4 oximeter. We calculated the correlation coefficients of the time traces between Pioneer and OxyPrem 1.4.
3. Results
3.1 Results of in-phantom validation
The Pioneer system measures the ToF of photons detected in all pixels of the Piccolo camera, i.e. all pixels are active at the same time. The onboard FPGA assembles this timing information from each pixel in histograms, the distribution of ToFs, which shows the number of detected photons per time bin. An example of such a response is shown in Fig. 4(a).
Fig. 4. ToF histogram representing the distribution of photons as function of time at the central pixel of the Piccolo array. It was measured on a homogeneous silicone phantom and is compared to the distribution expected from simulations (main plot a). Instrumental response function (inset b), normalized intensity (c) and phase delay (d) over the field of view, measured on the same phantom. Each pixel corresponds to an area of 1.06 × 1.06 mm2 on the surface of the phantom.
The typical distributions of the intensity and phase delay (i.e. ToF) over the circular FoV are shown in Fig. 4(c),(d). As expected, the intensity exponentially decreases with distance between the source and detector, whereas the phase delay increases linearly. Note that only the central part of the sensor was illuminated, thus some pixels appear outside of the FoV. There are a few ‘hot’ pixels in Piccolo, which are masked and therefore visible as pixels with missing data in the 3D surface plots (Fig. 4(c),(d)). 257 detectors provide valid data based on their SNR. This results in a total of ∼2’800 source-detector pairs.
The inner structure of the heterogeneous phantom was reconstructed from the acquired data. The optimization ended after 13 iterations when it reached the criterion of no further improvement. The whole process took <10 minutes. The location and shape of the inclusion were successfully recovered (Fig. 5(a)). The RMSE was 0.0023 mm−1 (Eq. (4)). After the segmentation, we obtained a high Dice similarity of 0.65 between the reconstructed image and the ground truth (Fig. 6(a)). As second analysis, we used the segmented volume as region-prior for the reconstruction of the optical properties of the two regions. The reconstructed and ground-truth optical properties are shown in Fig. 5(b) in comparison with the initial guesses. In addition, it is seen that a better agreement appears for ${\mu _a}$ than for ${\mu _s}^{\prime}$ because the contrast between the regions is higher in ${\mu _a}$.
Fig. 5. Distribution of reconstructed ${\mu _a}$ (a), with insets showing the ground truth (blue) and reconstructed (red) values along the dashed black lines drawn through the center of the inclusion. The dotted black frames show projections of the inclusion (a). Results of the region-based reconstruction for both ${\mu _s}\mathrm{^{\prime}}$ and ${\mu _a}$ for the two regions (bulk and inclusion) at three different wavelengths (b): comparison of the initial guess, ground truth and reconstructed values.
Fig. 6. Reconstructed finger region (white cloud) and the ground-truth finger shape (blue surface) (a) with the inset showing the top view (b). Note that the modeled volume is smaller than the real phantom. Hemodynamic responses during the arterial occlusion (c) and venous occlusion tests (d). The values were set to zero at the first data point. Shaded areas indicate an inflated cuff.
3.2 Results of the in vivo occlusion tests
The 3D image reconstruction took 74 s for one iteration and the optimal solution was reached after the 24th iteration. The segmented region shows good agreement with the true finger shape (Fig. 6(a),(b)). To calculate the Dice similarity, we converted the reconstructed distribution of the ${\mu _a}$ into a region of voxels and compared it to the reference model. The reference model was constructed from a 3D-scan of the mold that we employed to create the phantom. The Dice similarity between reconstructed 3D image and the reference model was 0.83.
The hemodynamic in the finger was recorded at a sample rate of 2.5 Hz. The acquisition of one complete set of data took 10 s: 0.4s x 11 sources x 2 wavelengths. The optical properties of the finger region were reconstructed for every measurement. [O2Hb] and [HHb] were calculated from the absorption coefficients of hemoglobin at two wavelengths [18]. The concentrations were set to zero at the first data point. The time series of [O2Hb] and [HHb] were smoothed with a 12th order Butterworth IIR filter (cut-off frequency 0.4 Hz) using the MATLAB zero-phase filtering function filtfilt. The changes in [O2Hb] and [HHb] during the arterial and venous occlusion are shown in in Fig. 6(c),(d) for both Pioneer and OxyPrem.
The correlation coefficients between the inflation status (0 for rest and 1 for inflated) of the cuff and Δ[O2Hb] were −0.68 (p<1%) and 0.68 (p<1%) for the arterial and venous occlusion tests, respectively. The correlation coefficients between the inflation status and Δ[HHb] were 0.66 (p<1%) and 0.50 (p<1%) for the arterial and venous occlusion tests, respectively. These high correlations that the expected hemodynamic changes were consistently detected.
Additionally, we compared the results of Pioneer to the CW oximeter OxyPrem 1.4. The correlation coefficients between them for Δ[O2Hb] were 0.8 (p<1%) and 0.63 (p<1%) for the arterial and venous occlusion tests, respectively. The correlation coefficients for Δ[HHb] were 0.61 (p<1%) and 0.39 (p<1%).
4. Discussion
In this paper, we described the TD NIROT system Pioneer, which is based on Piccolo camera. It has been shown that an increase in the number of source-detector pairs enables a higher spatial resolution of NIROT [15]. Thus, we designed the Pioneer system to take advantage of 1024 pixels of the ToF camera Piccolo, of 11 light sources and of a free choice of the wavelength in the NIR wavelength range. In this work, we used ∼2.8k source-detector pairs, only a fraction of the total number of available source-detector pairs. While Lyons et.al. demonstrated the potential of such a system in a specific 2D imaging task with a priori knowledge of the depth of a highly absorbing inclusion in transmission mode [15], we solved the more general problem of showing the 3D inner structure of an object in reflection mode and without priors.
The SPAD camera Piccolo enabled a novel concept of a TD NIROT system, which resulted in a very compact instrument. This differs fundamentally from previously developed NIROT systems, which used separate detectors, such as photomultiplier tubes, single SPADs or silicon photomultipliers, and optical fibers for emission and detection [6,8,9]. The novel Pioneer concept implements a scalable number of detectors, i.e., a larger SPAD camera with 36’288 pixels can be used [38]. In addition, the number of light sources may be expanded (up to 24 in our case). Notably, a similar arrangement of a ToF camera and light sources was employed by Sawosz et al. [12], but the image reconstruction was not performed. Our results demonstrate that it is possible to handle the abundance of data generated by the Pioneer system and to reconstruct reasonable images. Even dynamic measurements are feasible, which has been demonstrated by Jiang et al. in phantoms [39]. In this paper, we demonstrated the dynamic in vivo measurement with Pioneer for the first time. We employed repetitive measurements with a short exposure time of 400 ms to observe hemodynamic changes in an index finger during arterial and venous occlusions. The finger was located inside a silicone phantom with the optical properties similar to biological tissue. With the proposed 3D image reconstruction method, a volumetric image was obtained, and the finger region was correctly determined.
Furthermore, we demonstrated the capability of Pioneer to 3D image hemodynamic changes. With the reconstructed regional information, the optical properties were recovered at two wavelengths and converted to [O2Hb] and [HHb]. We observed reasonable changes in [O2Hb] and [HHb] between the occlusion and rest states. The non-imaging conventional oximeter reported data similar to those obtained by Pioneer. Why are these changes quantitatively similar? The conventional CW NIRS device OxyPrem does not provide images and quantifies the [O2Hb] and [HHb] assuming that the change is evenly distributed (homogeneity) in the whole object (the phantom containing the finger). Since changes in [O2Hb] and [HHb] actually only occur in the finger, we expected the OxyPrem to underestimate these changes. Since the Pioneer system provides images of the finger, we expected it to be more sensitive to the finger region and hence to yield larger concentration changes than OxyPrem. Possible explanations could be that i) the multi-distance approach of OxyPrem leads to a higher sensitivity to the tissue of the finger than being assumed by homogeneity, ii) OxyPrem assumes the optical properties of the neonate’s head, while the finger has different optical properties and iii) the two instruments did not measure simultaneously. Still the similar shape of the time series obtained by the two instruments confirms that the reconstructed images provide reasonable physiological information.
The no-prior image reconstruction is known to be inaccurate in the reconstructed values due to its ill-posed and ill-conditioned nature [18]. The tomograms tend to be blurry if no prior information is provided. As mentioned above, we have applied a strong regularization on the scattering. However, no further prior information was supplied to the reconstruction algorithm. Due to the tremendous amount of the acquired data, we were able to achieve high localization accuracy. This shows again that an increase in the amount of the obtained data enables higher image quality.
The inverse problem becomes well-posed and well-conditioned for the region-prior reconstruction [28–30]. Here we demonstrated that an image reconstructed without priors can be obtained first and then be employed as anatomical prior for a quantitatively correct image reconstruction based on this region prior. The challenge of this approach is to achieve an accurate segmentation result. We used a threshold of twice the median ${\mu _a}$ without priors for the segmentation and this yielded a reasonable structural 3D image. In a clinical environment, the region-prior can also be obtained from other imaging techniques, such as ultrasound, X-ray CT or MRI. The structure of an object can also be resolved in the future by machine learning methods [40,41]. The region-prior approach enabled us to reconstruct the optical properties of the two regions in the phantom at all tree wavelengths with a high precision.
We would like to emphasize, that we intended to challenge the Pioneer system and the reconstruction algorithms, and therefore intentionally set the initial guess values far from the correct ones. In the future, a more careful choice of an initial guess will facilitate the reconstruction process and is expected to result in an even more accurate solution.
5. Conclusion
To conclude, we reported in this paper the first in vivo results of our TD NIROT system Pioneer, which employs several cutting-edge technologies. Its Piccolo camera simultaneously counts photons in a large number of ToF detectors, which have a temporal resolution in the order of ps and a high sensitivity to near-infrared light. This enabled the collection of a massive amount of information from the object under investigation. The state-of-the-art reconstruction software enabled accurate and fast no-prior reconstruction of the structure of the object. Furthermore, we presented here the 3-step reconstruction approach: 1) The average optical properties of the object were determined and used as an initial guess. 2) The structure of the object was recovered without any priors. 3) This structural information was used to segment the object and to accurately determine the ${\mu _a}$ and ${\mu _s}^{\prime}$ in each region at all 3 wavelengths.
We employed this approach to image dynamic changes of tissue oxygenation in vivo. Pioneer successfully reconstructed the shape and the location of the finger positioned in a diffusive medium. With the exposure time of only 0.4 s (2.5 Hz), Pioneer detected the changes in the hemoglobin concentrations in the finger during arterial and venous occlusions. These results show a great potential of Pioneer to be applied in clinical environments where the detection of rapid changes in hemodynamics is important.
Funding
Swiss Cancer League (KFS-3732-08-2015); Swiss National Science Foundation (SNF) (159490 and 197079); Innovationspool of the University Hospital Zurich; National Competence Center in Biomedical Imaging, Swizerland (NCCBI); National Council of Science and Technology, Mexico (CONACyT) (CVU-627802).
Acknowledgment
We thank Dr. Juan Mata Pavia and Dr. Salvador Sánchez Majos for their contribution to the work at its early stages.
Disclosures
Martin Wolf: OxyPrem AG (I, P, S), Alexander Kalyanov: OxyPrem AG (E).
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. A. T. Eggebrecht, S. L. Ferradal, A. Robichaux-Viehoever, M. S. Hassanpour, H. Dehghani, A. Z. Snyder, T. Hershey, and J. P. Culver, “Mapping distributed brain function and networks with diffuse optical tomography,” Nat Photonics 8(6), 448–454 (2014). [CrossRef]
2. P. Antonio, C. Davide, M. Alberto Dalla, F. Andrea, S. Lorenzo, and T. Alessandro, “New frontiers in time-domain diffuse optics, a review,” J. Biomed. Opt. 21, 1–17 (2016). [CrossRef]
3. C. W. Lee, R. J. Cooper, and T. Austin, “Diffuse optical tomography to investigate the newborn brain,” Pediatr. Res. 82(3), 376–386 (2017). [CrossRef]
4. H. Zhao and R. J. Cooper, “Review of recent progress toward a fiberless, whole-scalp diffuse optical tomography system,” Neurophotonics 5, 011012 (2018). [CrossRef]
5. Y. Yamada, H. Suzuki, and Y. Yamashita, “Time-domain near-infrared spectroscopy and imaging: a review,” Applied Sciences-Basel 9(6), 1127 (2019). [CrossRef]
6. R. J. Cooper, E. Magee, N. Everdell, S. Magazov, M. Varela, D. Airantzis, A. P. Gibson, and J. C. Hebden, “MONSTIR II: A 32-channel, multispectral, time-resolved optical tomography system for neonatal brain imaging,” Rev. Sci. Instrum. 85(5), 053105 (2014). [CrossRef]
7. Y. Ueda, K. Yoshimoto, E. Ohmae, T. Suzuki, T. Yamanaka, D. Yamashita, H. Ogura, C. Teruya, H. Nasu, E. Ima, H. Sakahara, M. Oda, and Y. Yamashita, “Time-resolved optical mammography and its preliminary clinical results,” Technol. Cancer Res. Treat 10(5), 393–401 (2011). [CrossRef]
8. P. Taroni, A. M. Paganoni, F. Ieva, A. Pifferi, G. Quarto, F. Abbate, E. Cassano, and R. Cubeddu, “Non-invasive optical estimate of tissue composition to differentiate malignant from benign breast lesions: A pilot study,” Sci. Rep. 7(1), 40683 (2017). [CrossRef]
9. R. Re, E. Martinenghi, A. D. Mora, D. Contini, A. Pifferi, and A. Torricelli, “Probe-hosted silicon photomultipliers for time-domain functional near-infrared spectroscopy: phantom and in vivo tests,” Neurophotonics 3(4), 045004 (2016). [CrossRef]
10. A. D. Mora, D. Contini, S. Arridge, F. Martelli, A. Tosi, G. Boso, A. Farina, T. Durduran, E. Martinenghi, A. Torricelli, and A. Pifferi, “Towards next-generation time-domain diffuse optics for extreme depth penetration and sensitivity,” Biomed. Opt. Express 6(5), 1749–1760 (2015). [CrossRef]
11. Y. Mu and M. Niedre, “A fast SPAD-based small animal imager for early-photon diffuse optical tomography,” Conf. Proc. IEEE Eng. Med. Biol. Soc. 2014, 2833–2836 (2014). [CrossRef]
12. P. Sawosz, M. Kacprzak, N. Zolek, W. Weigl, S. Wojtkiewicz, R. Maniewski, and A. Liebert, “Optical system based on time-gated, intensified charge-coupled device camera for brain imaging studies,” J Biomed Opt 15(6), 066025 (2010). [CrossRef]
13. A. Tosi, A. Dalla Mora, F. Zappa, A. Gulinatti, D. Contini, A. Pifferi, L. Spinelli, A. Torricelli, and R. Cubeddu, “Fast-gated single-photon counting technique widens dynamic range and speeds up acquisition time in time-resolved measurements,” Opt. Express 19(11), 10735–10746 (2011). [CrossRef]
14. S. Lindner, C. Zhang, I. M. Antolovic, A. Kalyanov, J. Jiang, L. Ahnen, A. di Costanzo, J. M. Pavia, S. S. Majos, and E. Charbon, “A novel 32× 32, 224 Mevents/s time resolved SPAD image sensor for near-infrared optical tomography,” in Optics and the Brain, (Optical Society of America, 2018), JTh5A. 6.
15. A. Lyons, F. Tonolini, A. Boccolini, A. Repetti, R. Henderson, Y. Wiaux, and D. Faccio, “Computational time-of-flight diffuse optical tomography,” Nat. Photonics 13(8), 575–579 (2019). [CrossRef]
16. A. Puszka, L. Di Sieno, A. Dalla Mora, A. Pifferi, D. Contini, G. Boso, A. Tosi, L. Herve, A. Planat-Chretien, A. Koenig, and J. M. Dinten, “Time-resolved diffuse optical tomography using fast-gated single-photon avalanche diodes,” Biomed. Opt. Express 4(8), 1351–1365 (2013). [CrossRef]
17. C. Zhang, S. Lindner, I. M. Antolovic, M. Wolf, and E. Charbon, “A CMOS SPAD imager with collision detection and 128 dynamically reallocating TDCs for single-photon counting and 3d time-of-flight imaging,” Sensors-Basel 18(11), 4016 (2018). [CrossRef]
18. H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction,” Commun. Numer. Methods Eng. 25(6), 711–732 (2009). [CrossRef]
19. S. Lindner, C. Zhang, I. M. Antolovic, J. M. Pavia, M. Wolf, and E. Charbon, “Column-parallel dynamic TDC reallocation in SPAD sensor module fabricated in 180 nm CMOS for near infrared optical tomography,” in Int. Image Sensor Workshop, 2017), 86–89.
20. J. Jiang, A. Di Costanzo Mata, S. Lindner, C. Zhang, E. Charbon, M. Wolf, and A. Kalyanov, “Image reconstruction for novel time domain near infrared optical tomography: towards clinical applications,” Biomed. Opt. Express 11(8), 4723–4734 (2020). [CrossRef]
21. L. Ahnen, H. Stachel, S. Kleiser, C. Hagmann, J. Jiang, A. Kalyanov, S. Lindner, M. Wolf, and S. Sanchez, “Development and Validation of a Sensor Prototype for Near-Infrared Imaging of the Newborn Brain,” Adv. Exp. Med. Biol. 977, 163–168 (2017). [CrossRef]
22. K. Furutsu and Y. Yamada, “Diffusion-approximation for a dissipative random medium and the applications,” Phys. Rev. E 50(5), 3634–3640 (1994). [CrossRef]
23. R. Pierrat, J. J. Greffet, and R. Carminati, “Photon diffusion coefficient in scattering and absorbing media,” J. Opt. Soc. Am. A 23(5), 1106–1110 (2006). [CrossRef]
24. S. R. Arridge and M. Schweiger, “Photon-measurement density functions. Part 2: Finite-element-method calculations,” Appl. Opt. 34(34), 8026–8037 (1995). [CrossRef]
25. S. R. Arridge, “Photon-measurement density-functions Part 1. analytical forms,” Appl. Opt. 34(31), 7395–7409 (1995). [CrossRef]
26. P. K. Yalavarthy, B. W. Pogue, H. Dehghani, and K. D. Paulsen, “Weight-matrix structured regularization provides optimal generalized least-squares estimate in diffuse optical tomography,” Med Phys 34(6Part1), 2085–2098 (2007). [CrossRef]
27. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems 15(2), R41–R93 (1999). [CrossRef]
28. J. Jiang, L. Ahnen, A. Kalyanov, S. Lindner, M. Wolf, and S. S. Majos, “A new method based on graphics processing units for fast near-infrared optical tomography,” Adv. Exp. Med. Biol. 977, 191–197 (2017). [CrossRef]
29. J. Jiang, M. Wolf, and S. S. Majos, “Fast reconstruction of optical properties for complex segmentations in near infrared imaging,” J Mod Optic 64(7), 732–742 (2017). [CrossRef]
30. M. A. Mastanduno, J. Xu, F. El-Ghussein, S. Jiang, H. Yin, Y. Zhao, K. Wang, F. Ren, J. Gui, B. W. Pogue, and K. D. Paulsen, “MR-guided near-infrared spectral tomography increases diagnostic performance of breast MRI,” Clin. Cancer Res. 21(17), 3906–3912 (2015). [CrossRef]
31. J. Jiang, S. Lindner, A. Di Costanzo Mata, C. Zhang, E. Charbon, M. Wolf, and A. Kalyanov, “In phantom validation of time-domain near-infrared optical tomography pioneer for imaging brain hypoxia and hemorrhage,” Adv. Exp. Med. Biol. 1269, 341–346 (2021). [CrossRef]
32. W. Ren, J. Jiang, A. D. C. Mata, A. Kalyanov, J. Ripoll, S. Lindner, E. Charbon, C. Zhang, M. Rudin, and M. Wolf, “Multimodal imaging combining time-domain near-infrared optical tomography and continuous-wave fluorescence molecular tomography,” Opt. Express 28(7), 9860–9874 (2020). [CrossRef]
33. J. Jiang, A. Di Costanzo Mata, S. Lindner, E. Charbon, M. Wolf, and A. Kalyanov, “Dynamic time domain near-infrared optical tomography based on a SPAD camera,” Biomed. Opt. Express 11(10), 5470–5477 (2020). [CrossRef]
34. J. Jiang, A. Di Costanzo Mata, S. Lindner, M. Wolf, and A. Kalyanov, “Localization of deep ischemia and hemorrhage in preterm infants’ head with near-infrared optical tomography: a numerical case study,” Adv. Exp. Med. Biol. 1269, 131–136 (2021). [CrossRef]
35. U. Wolf, M. Wolf, J. H. Choi, L. A. Paunescu, A. Michalos, and E. Gratton, “Regional differences of hemodynamics and oxygenation in the human calf muscle detected with near-infrared spectrophotometry,” J. Vasc. Interv. Radiol. 18(9), 1094–1101 (2007). [CrossRef]
36. D. M. Hueber, S. Fantini, A. E. Cerussi, and B. B. Barbieri, “New optical probe designs for absolute (self-calibrating) NIR tissue hemoglobin measurements,” in BiOS ‘99 International Biomedical Optics SymposiumI li (SPIE, 1999), 14.
37. S. Fantini, M. A. Franceschini, and E. Gratton, “Semi-infinite-geometry boundary problem for light migration in highly scattering media: a frequency-domain study in the diffusion approximation,” J. Opt. Soc. Am. B 11(10), 2128–2138 (1994). [CrossRef]
38. C. Zhang, S. Lindner, I. M. Antolovic, J. M. Pavia, M. Wolf, and E. Charbon, “A 30-frames/s, 252 × 144 SPAD flash LiDAR with 1728 dual-clock 48.8-ps TDCs, and pixel-wise integrated histogramming,” IEEE J. Solid-State Circuits 54(4), 1137–1151 (2019). [CrossRef]
39. J. Jiang, A. di Costanzo, S. Lindner, M. Wolf, and A. Kalyanov, “Tracking objects in a diffusive medium with time domain near infrared optical tomography,” in Biophotonics Congress: Biomedical Optics 2020 (Translational, Microscopy, OCT, OTS, BRAIN), OSA Technical Digest (Optical Society of America, 2020), JTu3A.18.
40. Ö. Çiçek, A. Abdulkadir, S. S. Lienkamp, T. Brox, and O. Ronneberger, “3D U-Net: learning dense volumetric segmentation from sparse annotation,” in Medical Image Computing and Computer-Assisted Intervention–MICCAI 2016 (Springer International Publishing, 2016), 424–432.
41. O. Ronneberger, P. Fischer, and T. Brox, “U-Net: convolutional networks for biomedical image segmentation,” Lect. Notes Comput. Sc. 9351, 234–241 (2015). [CrossRef]
