Image reconstruction for novel time domain near infrared optical tomography: towards clinical applications

: Near infrared optical tomography (NIROT) is an emerging modality that enables imaging the oxygenation of tissue, which is a biomarker of tremendous clinical relevance. Measuring in reﬂectance is usually required when NIROT is applied in clinical scenarios. Single photon avalanche diode (SPAD) array technology provides a compact solution for time domain (TD) NIROT to gain huge temporal and spatial information. This makes it possible to image complex structures in tissue. The main aim of this paper is to validate the wavelength normalization method for our new TD NIROT experimentally by exposing it to a particularly diﬃcult challenge: the recovery of two inclusions at diﬀerent depths. The proposed reconstruction algorithm aims to tackle systematic errors and other artifacts with known wavelength-dependent relation. We validated the device and reconstruction method experimentally on a silicone phantom with two inclusions: one at depth of 10 mm and the other at 15 mm. Despite this tough challenge for reﬂectance NIROT, the system was able to localize both inclusions accurately.


Introduction
Imaging the body with harmless near infrared optical tomography (NIROT, also called diffuse optical tomography) is appealing to medical professionals and patients.Due to its intrinsic sensitivity to tissue oxygen saturation, it has potential applications for diagnosis of ischemic strokes, hypoxic tumors and early detection of neonatal brain injury [1][2][3].A typical NIROT device consists of sources emitting light and detectors capturing the backscattered photons at some distances.Temporal and spatial information of light is later utilized for image reconstruction [4].In general, achieving both high resolution and accuracy in image reconstruction is challenging.Researchers have been investigating methods to improve image quality by increasing the variety of measurements, for example, with denser arrangement of sources and detectors [5].However, such systems so far have been usually bulky and require many optical fibers which limits their usability in a clinical environment.One way is to enrich the information is by employing the temporal information of light [6].Such systems are called time domain (TD) NIROT.It has been proven that TD data with proper signal to noise ratio improves the depth sensitivity and resolution [7] [8].The recently developed TD NIROT Pioneer system utilizes a 32×32 array of single-photon avalanche diodes (SPADs) measuring time of flight (TOF) for 12.5 ns with time resolution 116 ps [9].It is equipped with 11 sources of pico-second ultra-fast super-continuum laser radiation.Pioneer provides 2.88 × 10 6 measurement points for a single wavelength [10].Such an approach enables an unprecedented amount of data, which will lead to more accurate images [11].Despite the advantages of the TD NIROT, the calibration process is non-trivial for a SPAD camera based imaging system.NIROT devices are usually calibrated on a homogeneous phantom [12].However, measurement artifacts occur in TD data even for a system of which the instrumental response function (IRF) is well calibrated on phantoms [13].For example, when the imaging device is applied on a surface with different curvatures, the varying distances cause different time offsets in the TD data.The time offsets alter the previously calibrated IRFs and the use of inaccurate IRFs will deteriorate the reconstructed image quality.These artifacts must be eliminated in a real case scenario and therefore a robust calibration algorithm has to be developed.Previously, relative images with different source positions were used as the calibrated data for a fiber and charge-coupled device (CCD) camera based system to eliminate the artifacts from the imperfect calibration [14].Nevertheless, it demands to know the accurate relative power of each fiber optic source and relative offsets in time delay caused by, e.g.different fiber lengths, for a SPAD system.
We have previously proposed an auto-calibration method called wavelength normalization (WN) for region-based image reconstruction [13].It utilized ratios of Fourier transformed TD data at different wavelengths without tedious calibration in advance.This type of calibration also has been studied for the recovery of the optical properties in near infrared spectroscopy (NIRS) [15].This way of calibration naturally compresses the signal from original data volume for two wavelengths to one.Fortunately, the data volume for solving the inverse problem for both a region-based NIROT and NIRS is not demanding.Because there are only two unknown (absorption µ a and scattering property µ s ) for NIRS and a few values dependent on the number of segments for the region-based NIROT.While reconstruction methods have been established for TD NIRS [15] and structure-guided TD NIROT [13], 3D image reconstruction for NIROT without any prior has not been reported.This is a non-trivial problem because the number of unknowns is equivalent to the number of quantization element, which can be thousannds and taking the ratios lead to a more ill-posed problem.In this paper, we designed an image reconstruction method based on WN for the TD NIROT Pioneer.The aim was to validate the method experimentally by exposing it to a particularly difficult challenge: recovery of two inclusions at different depths from measured reflectance [16].

System description
We have recently developed a TD NIROT system named Pioneer and the schematic of the system is displayed in Fig. 1(a).In this Pioneer system, pico-second light pulses of a broad spectral range are generated by a super-continuum laser SuperK Extreme EXR-15 (NKT, Denmark) at the repetition rate of 80 MHz.Wavelengths in the near infrared range are selected with an acousto-optic tunable filter (AOTF) which enables multispectral measurements.A splitter guides 1% of the output light to a power meter for laser stability control and the remaining 99% of the output light enters an optical switch to enable multiple source locations.To make the system friendly for clinical applications, we designed a tube-like probe of inner diameter ∼2.5 cm (Fig. 1(b)).The interface of the tube to tissue was a ring-shaped structure (schematics illustrated in Fig. 1(a)).11 metallic ferrules for holding 11 multi-mode 62.5 µm core fibers (Thorlabs, US) were fixed in the rigid ring evenly at a distance of 1 cm around the tube.To ensure a comfortable touch, the rigid structure was coated with soft black biocompatible silicone.Transparent silicone windows were inserted at each fiber location to direct the light to the tissue.The detailed manufacturing process of the ring is described in [17].A 32×32 SPAD camera is placed on the other side of the tube.It features a low dark count rate (DCR) of 113 counts/s and a wide spectral range and high photon detection probability (PDP) of 12% at 800 nm.The sensor includes integrated timing electronics in the form of 128 50ps time-to-digital converters (TDCs) [9], such that the arrival time of each photon at any pixel is measured with high time resolution.The maximum number of photon throughout is 465 million counts per second for the whole sensor and the dead time between events is 1.2 ns [18].An IR-corrected lens SY110M 1.68mm F1.8 (Theia Technologies, US) projects the circular field of view (FOV) defined by the inner boundary of the tube on the SPAD array.Eventually, the Pioneer system will provide > 10k source-detector pairs for each measurement.We developed a software, based on Qt with Cpp, to control the hardware and the data acquisition.The data processing was performed in Matlab.More technical details of the system are described in [10].

Image reconstruction
The TD data or the enormous set of histograms obtained with the Pioneer device for the sourcedetector pairs are utilized for image reconstruction, i.e., the inverse problem.It relies on an accurate description of light propagation in tissue, which is called the forward problem.The solution to the forward problem provides the survival photon density and timing at each detector.The survival photon density distribution φ(r, t) at position r and time point t is discribed by the time-dependent diffusion equation (DE) as below [19]: with a boundary condition of: where q is a source term and in TD NIROT it is a pulsed laser source.The model parameters include the absorption coefficient µ a , diffusion coefficient κ(r) = [3(µ a + µ s )] −1 , with reduced scattering coefficient µ s and speed of light c 0 (r) in medium.The boundary condition contains the refractive index mismatch ζ and the outward boundary normal ∂ν.For complex geometries as in most clinical applications, the DEs are solved numerically with finite-element-methods (FEM) which is based on mesh generation techniques for dividing a complex tissue volume, for example brain or breast, into small elements.An FEM based simulator NIRFAST is therefore utilized to model the light transport [20].TD data in the form of histograms was transformed into frequency domain (FD) with MATLAB s fast Fourier transform fft function and sampled at a number of selected frequencies for image reconstruction.Frequencies up to Nyquist = f s /2 are theoretically suitable for reconstruction.However, the very high frequency region tend to be noisy and therefore are excluded.The forward problem in FD is described as follows [20], where Φ is the fluence and Q is the source term in FD.
The image reconstruction or so-called inverse problem can be written as follows and for simplicity we omit the parameters r and ω in the forward result Φ(r, ω).
where µ is the three dimensional distribution of optical properties and the goal of the image reconstruction is to find the optimal µ * where the calculated forward result Φ S (µ * ) matches the measured data Φ M .Γ(µ) is the regularization term and we use the Tikhonov regularization in this work.The optimization algorithm used iteratively updates µ to minimize the difference between Φ S and Φ M (µ) and it terminates once the alteration of µ generates larger deviation than the one at the previous iteration.To resolve small details of tissue, a high signal-to-noise ratio (SNR) of the measured TD data is necessary.The image quality is impaired by systematic errors, e.g. the instrumental response function, inaccurate calibration of source detector positions, and superficial obstructions such as skin features (moles, birthmarks, hair).The contaminated signal from the true signal f (t) can be expressed as f = αf (t) * g(t), where α stands for the damping factor and g(t) means the instrumental response function.It is non-trivial to calibrate the system to extract the true signal from the contaminated measured data, which are encountered in clinical scenarios.
We have previously studied an auto-calibration method called wavelength normalization (WN) for NIROT reconstruction with region-based prior [13].In this section, we applied the method to tackle a complete 3D reconstruction problem without the need of calibration in advance for a 3D image reconstruction without any shape prior as in [13] or depth prior in [21].We use a wavelength normalized (WN) data -the ratio of forward results at different wavelengths for the inverse problem where Φ λ i is the reflectance for all source-detector pairs at wavelengths λ i , i = 1, 2, . . ., n in the multi-spectral measurement.When we treat the TD data in Fourier domain, the convolution becomes simple multiplication.The WN data is therefore as follows, To eliminate the influence from IRF g(t) or G(ω) in frequency domain, the system requires characteristics of wavelength-invariance in the time response.Our Pioneer system features wavelength-invariance in the IRF as illustrated in Fig. 2 where the laser-SPAD responses are almost identical at different wavelengths.The IRFs are therefore effortlessly cancelled out by taking the ratio.The damping factor α is usually wavelength invariant and it also can be eliminated simultaneously.Consequently we obtain a true signal R i from Ri automatically.Thereafter, we reformulate the inverse problem as the optimization process to find the optimal set of µ distribution for both wavelengths µ λ i ,λ 0 * by minimizing the difference between reference measured ratio data R M i and the ratio of the forward results R S i (µ λ i ,λ 0 ) with estimated variables µ λ i ,λ 0 , being the updated optical properties for both wavelengths λ i and λ 0 .
For simplicity, we will assume i = 1, but the solution can be easily extended to multiplewavelength measurement.The target function for FEM-based image reconstruction algorithms therefore can be formulated as follows.
Tikhonov regularization is applied with parameter β.The goal of the inverse problem is to find the optimal set of optical properties µ in all N n nodes for both wavelengths λ 1 and λ 0 in order to match the simulated ratio data R S 1 to the measured ratio R M 1 for N m measurements, i.e. source-detector pairs.(µ λ 0 ,λ 1 ) 0 are the initial guess of optical properties for both wavelengths.To solve the above problem, a Newton like iterative process is taken to update the optical properties by quantities determined by where I is the identity matrix and J R is the Jacobian for the ratio data with respect to optical properties.The Jacobian J R consists of derivatives of logarithmic amplitude ln |R 1 | and phase θ(R 1 ) with respect to the optical properties µ λ 0 , µ λ 1 at two wavelengths.ln denotes the natural (base e) logarithm.In this paper, the image reconstruction was performed with the assumption of the same scattering coefficient in the whole volume and µ therefore is written as µ a .We omit the explicit dependency of frequency in the notation for convenience.These quantities can be calculated from Jacobians for both wavelengths.We write the derivatives to form the Jacobian matrix with respect to µ a for ratio data of log amplitude (lnI λ 0 , lnI λ 0 )) and phase (θ λ 0 , θ λ 1 ) as follows, The previously developed reconstruction algorithm with region priors includes only ∼10 unknowns [13].The proposed new method for a 3D volume without any shape prior needs to solve thousands of unknowns, which is obviously more ill-posed.In addition, the TD data volume is compressed by taking the ratios of measured data at two wavelengths to obtain clean signals, which further makes the inverse problem more difficult to solve.However, the high number of channels and excellent sensitivity of the imager as in the Pioneer system provide rich enough information and therefore the accurate 3D reconstruction with the auto-calibrated method can be achieved.

Phantom experiment
To validate the NIROT system, reflected light was measured on a phantom of two small inclusions at two different depths (Fig. 3(a)).This poses a high difficulty for NIROT in reflection mode, because the sensitivity decreases with increasing depth for NIROT in reflection mode [16].
As illustrated in Fig. 3(b), the silicone phantom contains two spheres of radius 5 mm embedded at 10 mm and 15 mm respectively.To produce the phantom, we firstly 3D printed molds for both two spherical inclusions and the cylinder bulk.We then prepared two mixtures of silicone material which have different absorption coefficients and similar scattering properties.The absorption of Fig. 2. Laser-SPAD response functions at four sample wavelengths measured in reflection mode where the laser light was pointed to a white paper and the reflected light was captured by the SPAD camera; Note that the wavelength dependent dispersion in time, and additional time delays caused by different fiber lengths were eliminated by numerically applying pre-measured fixed time shifts to the TD data.The functions are so similar that with wavelength ratios they can be calibrated out easily. the material was tuned by adding black ink and RAL 6004 Blue Green (Wacker, Germany) in the silicone mixture and scattering was tuned by adding RAL 9010 White (Wacker, Germany).The silicone mixture of higher absorption (mixture S2) was cast into the spherical mold to produce two inclusions shown on the right of Fig. 3(a).These two spheres were then fixed to the cylindrical mold with fishing lines at the desired locations.Eventually a big cylindrical phantom (Fig. 3(a) Left) was produced by casting the less absorbing silicone (mixture S1) into the cylindrical mold with the two already embedded inclusions.To verify the target optical properties, a commercial device Imagent (ISS Inc., Champaign, IL, USA) was utilized.The device measured optical properties of a large homogeneous phantom.Therefore we cast two homogeneous phantoms from mixtures S1 and S2 in the shape of the cylindrical phantom (Fig. 3(a) Left).Two wavelengths λ 1 = 689 nm and λ 2 = 725 nm were chosen for the measurement to validate the WN and the measured optical properties with the Imagent are listed in Table 1.
In the phantom experiment, the detection probe was placed on the surface of the phantom as depicted in Fig. 1(b) and data acquisition was performed at 689 nm and 725 nm.The SPAD camera collected 2 × 10 8 photons for each source which took ∼ 1 minute.During the process, the powermeter showed stable laser output.The variation of the average intensity of the laser did not exceed 2%.The variation of time delay is lower than the temporal resolution of the sensor and did not affect experimental results, i.e. it was below ∼100 ps.Eventually, 2431 histograms of good signal to noise ratio from 221 detectors within the circular FoV and 11 sources were selected.Dark photon counts were removed from each measured histogram by substracting its median value.Then the histograms which still contain errors from surface features such as dust and are distorted by the IRF were transformed into frequency domain in the image reconstruction.In this work, we used the frequencies 100 MHz, 200 MHz and 500 MHz.The initial regularization parameter was 10.We also simulated the phantom experiment in order to validate the feasibility of the autocalibrated image reconstruction with the same provided data volume.The current WN based reconstruction aims to recover a 3D distribution of optical properties at two wavelengths, although employing the data volume equivalent to a single wavelength reconstruction.This makes the reconstruction more difficult to achieve because the doubled number of unknowns increases the dimension of the optimization problem.To help the optimization algorithm to approach the correct solution easier, the absolute forward results at both wavelengths Φ i with lower weights η i (0<η i <1, i = 0, 1) can be combined with the ratio data to enrich the target data as Therefore, the final merit function to minimize is The values of weights depend on the signal quality of Φ i .The higher the η i is used, the more the inverse problem leans toward the original one described by Eq. ( 4).In this work, η i = 0.1 was used.

Image evaluation
The quality of reconstructed images were assessed with two widely applied metrics: mean squared error (MSE) and contrast-to-noise ratio (CNR).The metrics MSE and CNR were applied on the reconstructed 3D map of μa and ground truth µ a defined as follows: MSE measures the difference between the actual values and the estimated values predicted by the reconstructions.CNR indicates the ability to visualize different tissues or inclusions from the background.Therefore, the image quality is in general better with lower MSE and higher CNR.

Numerical validation
The phantom measurement was simulated in NIRFAST to validate the proposed imaging method.A cylindrical mesh of Ø60 × 35 mm 3 was generated to represent the phantom in the measurement and two spherical regions at two distinctive depths 10 mm and 15 mm with lateral separation of 10 mm were defined.Both spheres had the same radius of 5 mm.Note that the size of the simulated cylinder was smaller than the actual silicone phantom in the experiment to facilitate the simulation process.The optical properties were assigned to the bulk and to the 2 spheres according to Table 1 for the bulk and two spheres.As in the experiment, 11 light sources and 221 selected detectors within a circle of diameter = ∼17.5 mm were placed on the top surface of the mesh (Fig. 3(c)).A fine mesh of 100210 nodes was used to simulate the forward process.To decrease the ill-posedness of the inverse problem, a second coarser mesh of 942 nodes was generated and employed for the reconstruction.Technically, the use of a coarser mesh is computationally easier in terms of memory usage and speed.
The image reconstruction was completed with an initial guess of a homogeneous volume of the same µ s and µ a as the bulk.The whole process underwent 20 iterations.To later compare with the experimental results, we denote the reconstruction results from simulated data as "Sim".Cross-sectional distributions of reconstructed values were plotted in Fig. 4(b)(f)(j)(n) for the two wavelengths.The shapes and positions of the inclusions were correctly recovered compared to the ground truth values (GT) (Fig. 4(a)(e)(i)(m).The successful 3D reconstruction by the simulation proved the theoretical validity of the proposed imaging methods.

Experimental results
In the experiment, ToFs were acquired for 11 sources at two wavelengths.We show an example of ToFs from 3 detectors (pixels) located at different distances from the active light source (Fig. 5).Two types of reconstruction were performed: non-calibrated and auto-calibrated methods.For the non-calibrated case, TD data measured at 689 nm and 725 nm were employed directly without decoupling the IRFs or removal of artifacts such as abnormal weak signals in pixels due to dirt or dust on the phantom surface.Non-calibrated histograms were transformed in frequency domain as the target forward results.Images were reconstructed separately for the measurement at 689 nm and 725 nm.The homogenous case was the initial guesses where we assigned the same µ B a and µ B s for all elements in the mesh.The same number of iterations of 20 and regularization parameter of 10 were used for both non-calibrated and auto-calibrated data.For the auto-calibrated case the ratio of the Fourier data at two wavelengths combined with weighted forward results was used as the target forward results.The inverse problem solver was the same for both the non-calibrated and the simulated case.For simplicity, we denote results related to the non-calibrated case as "Meas" and auto-calibrated case as "MeasC".
We plotted the cross-section of the reconstructed µ a at depth 10 mm and 15 mm for both wavelengths in Fig. 4. The inclusions were clearly visible in both Sim (Fig. 4(b)(g)(j)(n)) and MeasC (Fig. 4(d)(h)(l)(p)) for both depths.However, the non-calibrated reconstruction failed to generate accurate shape of the inclusion at 10 mm depth (Fig. 4(c)(k)) and completely failed to generate the inclusion at the depth of 15 mm (Fig. 4(g)(o)).The one-dimensional distributions of µ a were also depicted in Fig. 6.The values estimated from the auto-calibrated ratio data were closer to the ground truth than the non-calibrated ones at both depths and both wavelengths.The MSE and CNR of the reconstructed results are displayed in Fig. 7. MeasC outperformed Meas in both evaluations and is close to the Sim.

Discussion
This study presented an WN based image reconstruction method for the novel TD NIROT system Pioneer based on a SPAD camera.We validated the system in reflection-mode measurements on a silicone phantom with two small absorbers embedded at depths of 10 mm and 15 mm.The sensitivity, known as the amount of changes in the measured signals caused by the changes in optical properties, is lower in deeper tissue [16].This makes the optimization algorithm in the image reconstruction easily approach a solution, which only recovers the shallower inclusion and omits the deeper inclusion.However, in our results both inclusions were evidently visible and their positions were accurately recovered.The spherical shape of the inclusion at 10 mm was perfectly reconstructed whereas the deeper one suffered some extent of deformation, which again was caused by the sensivity difference.One can further improve the depth sensivity by applying different optimization strategies.For example, the regularization method in this work was the common Tikhonov regularization.It is possible that the use of WN for image reconstructions with alternative regularization methods, e.g.spatial variant regularization [22] or region-based methods [23], may yield higher image accuracy.
Benefiting from the WN, Pioneer is able to be used without tedious prior calibration, which makes it more robust in a clinical environment.WN eliminates the IRFs automatically by taking the ratios of the Fourier transformed TD data obtained at multiple wavelengths.As our results show, the 3D image reconstruction based on WN clearly outperforms the non-calibrated one based on single wavelength.Despite the fact that in reality the IRF is never perfectly spectrally nor temporally invariant, the experimental results showed that the image reconstruction was much improved by WN without compensating for this variation.The image reconstruction might be even better if this effect was compensated better.Even though the inverse problem is more ill-posed due to the reduction in the data volume by taking the ratio, the rich information provided with TD NIROT enables an accurate solution of the inverse problem.

Conclusion
We propose a robust reconstruction method for the novel SPAD camera based TD NIROT system Pioneer.The large volume of data obtained with the SPAD camera enables recovering depth information in reflectance measurement.Reflectance is usually more accessible in clinical applications.The image reconstruction based on wavelength normalization approach avoids a tedious calibration in advance.It handles different sources of artifacts such as instrumental response functions and, we speculate, other linear artifacts, e.g.surface features.We validated the method on a phantom with two inclusions embedded at different depths and generated images of high quality.This brings the Pioneer system closer to the real clinical applications.

Fig. 1 .
Fig. 1.(a) Schematics of Pioneer system; and (b) picture of the Pioneer probe placed on a cylindrical silicone phantom.

Fig. 3 .
Fig. 3. (a) The cylindrical silicone phantom used in the experiment has a diameter of 114 mm and a height of 55 mm and both two small spheres have a radius of 5 mm.(b) The two inclusions were embedded at 10 mm and 15 mm, respectively and the lateral distance between them was 10 mm.(c) Mesh created in simulation for image reconstruction.Note that we made it smaller than the actual volume in (a): cylinder of diameter 60 mm and height 35 mm.221 detectors were selected within a circle of diameter = ∼17.5 mm and all 11 sources were utilized.

Fig. 5 .
Fig. 5. ToF histograms from three pixels (locations shown in the inset) for the phantom measurement at 689 nm.The time bin size is 48.8 ps.