High Resolution, Deep Imaging Using Confocal Time-of-flight Diffuse Optical Tomography
Yongyi Zhao, Ankit Raghuram, Hyun K. Kim, Andreas H. Hielscher, Jacob T. Robinson, Ashok Veeraraghavan
11 High Resolution, Deep Imaging Using ConfocalTime-of-flight Diffuse Optical Tomography
Yongyi Zhao, Ankit Raghuram, Hyun K. Kim, Andreas H. Hielscher,Jacob T. Robinson, and Ashok Veeraraghavan
Abstract —Light scattering by tissue severely limits both how deep beneath the surface one can image, and at what spatial resolutionone can obtain from these images. Diffuse optical tomography (DOT) has emerged as one of the most powerful techniques for imagingdeep within tissue – well beyond the conventional ∼ > meanfree paths). In addition, we demonstrate that two additional innovations: focusing on confocal measurements, and multiplexing theillumination sources allow us to significantly reduce the scan time to acquire measurements. Finally, we also rely on a novelconvolutional approximation that allows us to develop a fast reconstruction algorithm achieving a 100 × speedup in reconstruction timecompared to traditional DOT reconstruction techniques. Together, we believe that these technical advances, serve as the first steptowards real-time, millimeter resolution, deep tissue imaging using diffuse optical tomography. Index Terms —Time-of-Flight Imaging, Diffuse Optical Tomography, Confocal, Time Binning, Fluorescence Imaging (cid:70)
NTRODUCTION L IGHT scattering by tissue is the primary challenge limit-ing our ability to exploit non-ionizing, optical radiationin the 400-1000 nm wavelength range, to perform high-resolution structural or functional imaging, deep inside thehuman body. Most existing techniques, including confocalmicroscopy, two-photon (2P) microscopy and optical coher-ence tomography (OCT), exploit only the ballistic (or single-scattered) photons and can only be used to image withinthe ballistic regime ( < mean scattering lengths deep)[1], [2]. This limits imaging to approximately the top 1-2millimeters of tissue surface (as mean scattering lengths intissue is ≈ − µ m range [1], [3]) as seen in Fig. 1a. Manyapplications (both clinical and scientific) require imaging atmuch higher depths of penetration than can be achieved byremaining within the ballistic regime.Diffuse optical tomography (DOT) [4] has emerged asthe one of the most promising techniques (another beingphoto-acoustic tomography [5]) for high-resolution imagingdeep within tissue, in the diffusion regime (i.e., > meanscattering lengths). The idea in DOT is that even in thediffusive regime, where light-paths are highly random, thereare statistically predictable structures in its distribution inspace, and this regularity can be exploited if sufficientdiversity of measurements are obtained. DOT uses an arrayof sources and detectors placed over the imaging volume –and the light transport data acquired between each source-detector pair provides the required measurement diversity. • Y. Zhao, A. Raghuram, J. T. Robinson, and A. Veeraraghavan are withthe Department of Electrical and Computer Engineering, Rice University,Houston, TX 77005.E-mail: { yongyi, ar89, jtrobinson, vashok } @rice.edu • H. K. Kim and A. H. Hielscher are with the Department of BiomedicalEngineering, New York University, New York City, NY 11201E-mail: { hk3363, ahh4614 } @nyu.edu Challenges.
In spite of its promise, DOT systems today re-main severely limited. Firstly, existing DOT systems providelow spatial resolution. Most are limited to cm -scale spatialresolutions because of a combination of factors includinglack of sufficient measurement diversity, modeling inaccu-racies, and low SNR measurements (Fig. 1b). Second, thesequential nature of DOT measurement process introducesa trade-off between SNR and capture time, further limitingresolution (and quality) when it comes to imaging dynam-ics. Third, DOT reconstruction algorithms have to contendwith solving large-scale optimization problems with poten-tially millions of variables and therefore tend to be quiteslow, precluding real-time performance. Our goal, in thispaper, is to directly address these limitations. Key Ideas.
Our approach leverages three key ideas.
Key Idea 1 - Increased measurement diversity provided bytransients.
The primary cause of reduced spatial resolutionis understood to be the limited measurement diversity.Increasing the number of source-detector pairs improvesspatial resolution but this tends to saturate beyond a point.It becomes essential to enhance diversity of measurementsby adding additional dimensions. Time of travel betweensource and detector may be a promising additional dimen-sion that is significantly beneficial since many of the surfacescattered background photons tend to have a significantlyshorter travel time than most of the deep penetrating signalphotons that interact with the tissue of interest [6]. Wedemonstrate that exploiting this additional transient di-mension (by capturing transient histograms between everysource-detector pair), provides sufficient increase in mea-surement diversity to obtain mm spatial resolution even inthe diffusive regime. Key Idea 2 - Reduced capture time through multiplexedmeasurements.
DOT measurements are typically acquired a r X i v : . [ ee ss . I V ] J a n sequentially and this establishes a trade-off between capturetime and SNR. We propose that multiplexed acquisition,wherein multiple light sources are ’on’ simultaneously, im-proves measurement SNR. With a reconstruction algorithmthat can de-multiplex these measurements, we show thatsource multiplexing can provide a × - × reduction incapture time compared to traditional sequential DOT. Key Idea 3 - Real-time reconstruction using a novel con-volutional approximation.
Traditional DOT reconstruction al-gorithms are already computationally intensive — andwith the ∼ × increase in measurement dimensional-ity imposed by capturing transient information, this bur-den is severely exacerbated precluding any hope for nearreal-time reconstruction performance. We propose a novelconvolutional approximation for multiplexed (and non-multiplexed), confocal time-of-flight diffuse optical tomog-raphy and utilize this approximation to develop a fast,real-time reconstruction algorithm (which is a × - × speedup). Outcomes and Potential Impacts.
The primary outcomethat we are able to demonstrate is that we show millimeterspatial resolution in the diffusive regime ( > mean scat-tering lengths). This, in itself, opens up a variety of newclinical and scientific imaging applications. In particular,we believe that non-invasive brain imaging (both structuraland functional) is a critical application domain. As skullseverely attenuates acoustic waves making through-skullphoto-acoustic tomography difficult [7], DOT already is thepredominant technology for this application. Improving theachievable spatial resolution will provide us better speci-ficity potentially allowing us to image columnar fields inthe brain for the first-time. The secondary outcome is thefirst demonstration of a real-time reconstruction algorithmfor time-of-flight DOT. In addition, we also show that mul-tiplexing can significantly reduce capture-time in DOT. Fi-nally, we develop two different versions of the algorithm forboth fluorescence and absorption imaging, and demonstratereal results for both modes - expanding the potential scopeof applications. Limitations.
All our current demonstrations are in tissuesamples and phantoms (both fluorescence and absorption).We are actively working towards demonstrating the feasibil-ity in real biological tissue, as we realize that there are addi-tional challenges such as reduced fluorescence/absorptioncontrast, increased biological noise, and motion (especiallywhen imaging in vivo ) that we might need to address beforethe technology can reach its promised potential.Our current prototype is sub-optimal in many respects.While traditional DOT systems have a wearable form-factor, our laboratory prototype uses benchtop optics, witha scanned laser head and a single detector being scannedto mimic a detector array. While compact systems with asimilar wearable form-factor to existing DOT systems isindeed possible, this requires fabrication of an array ofSPAD detectors and corresponding laser diodes, with acommon shared clock – something that is beyond the scopeof this paper. Since our benchtop prototype scans a singlepixel to emulate a detector array, the total scan time of oursystem is increased by a factor that is proportional to thetotal number of detectors being emulated (typically × - × in our results). This along with scanning inefficienciesmean that the scan time in all our results are in the several Fig. 1:
Imaging depth and spatial resolution of DOT tech-niques. (a) Approximate imaging depth of optical imagingtechniques. Ballistic imaging techniques such as OCT, confocalmicroscopy, and 2P microscopy cannot image past ∼
15 meanfree paths (MFPs) . DOT approaches can achieve 10s-100sof MFPs (b) Approximate spatial resolution of different DOTtechniques. Our technique is the only method to demonstrate 1 mm spatial resolution. seconds to minutes range, precluding any ability to showreal-world dynamics in our real results. We are workingtowards realizing a compact, fabricated, prototype for awearable, brain imaging system and are hopeful that wecan demonstrate that system in action in about a year. ELATED W ORK
Imaging within the ballistic regime.
The fraction of pho-tons that enter a tissue and remains ballistic decreasesexponentially with the thickness of the tissue being imaged.Even after just 3 mean scattering lengths, the fraction ofphotons that are ballistic become 1 in 20 according to theBeer-Lambert Law [8]. As a consequence, even at thesedepths, techniques such as direct imaging, brightfield imag-ing, or fluorescence imaging that do not actively filter outthe scattered photons get overwhelmed by the backgroundfrom these multiply-scattered photons reducing the imagingcontrast to below the sensor sensitivity thresholds [1].Beyond this depth active means of rejecting the multiply-scattered photons are needed. Confocal microscopy uses aset of matched pinholes to reject a large fraction of the scat-tered light, and typically extends imaging to about 6 meanscattering lengths (1 in 400 photons are ballistic) [9]. Multi-photon microscopy techniques including 2P microscopy,rely on the non-linear excitation process to confine fluores-cent emission, and these techniques may allow imaging tobe performed as deep as 16 mean scattering lengths (1 in ∼ ∼
480 million photons remain ballistic. Going beyond asyou encroach into the diffusive regime ( ∼
50 mean scatteringlengths and beyond), techniques that rely exclusively on bal-listic (or single-scattered) photons are completely infeasibleas less than 1 in . × photons are ballistic. Beyond the ballistic regime.
As you move beyond theballistic regime, the fraction of ballistic photons is so smallthat relying on them exclusively is insufficient. Therefore, itbecomes imperative, to find ways to model the localization(even if it is only partial) of the scattered photons and exploitthese scattered photons as well.
Diffuse Optical Tomography (DOT).
DOT originated in the1990s as a way of detecting absorption changes in medical imaging applications [6]. Traditional DOT systems utilize anarray of near-infrared, continuous-wave (CW) light sourcesilluminating the tissue, resulting in multiply scattered pho-tons that arrive at an array of detectors [4], [11]. Modelsof photon propagation physics could then infer local ab-sorption and scattering properties within the tissue from themeasurements captured by the detectors. Early applicationsof DOT included imaging tumors for breast cancer andmonitoring brain bleeds for infants [4], [6]. Transmittancemeasurements of these geometries provided absorption in-formation on the whole volume of interest. However, theadult brain and internal organs must be imaged in reflectionmode due to the strong scattering and absorption properties,or limited access to the tissue of interest [6]. In the rest ofthe paper, we will refer to continuous wave DOT as DOTfor simplicity.Recent advances in DOT have been focused on algorith-mic improvements resulting in higher spatial resolution [11]and development of wearable devices [12], [13], [14]. Themost significant drawback of DOT is depth sensitivity. Fordeeper penetration in reflection mode, source and detectorseparations must be farther apart, reducing the SNR of themeasurements [11]. Frequency- and time-domain DOT havebeen developed to counteract these problems. While bothfrequency- and time-domain can capture the same informa-tion, time-domain DOT (TD-DOT) can make measurementsfaster, albeit with more expensive hardware [15].
ToF-DOT.
ToF-DOT (or TD-DOT) uses a high-power, nar-row pulse-width laser and a fast-gated detector to capturetransient light transport data [16], [17]. These transientscontain photon arrival time information for each source-detector pair, providing an additional dimension of infor-mation to improve depth sensitivity [16]. The emergenceof single-photon avalanche diodes (SPADs) in recent yearscoupled with on-chip time-correlated single photon count-ing (TCSPC) electronics has allowed for fast-gated, largedynamic range, ps resolution transient measurements inreasonable acquisition times, making ToF-DOT a promisingtechnology to explore [6]. In addition, hardware improve-ments are making wearable ToF-DOT systems feasible, andthere is some early work towards that direction [18], [19].While the predominant application of DOT and ToF-DOT has been deep tissue imaging (especially breast can-cer and through-skull imaging), the technology could po-tentially be used for other applications including imagingthrough thick scatterers, or imaging around a corner. Re-cently, [20] demonstrated a that a 3D image can be acquiredthrough thick scattering media using an imaging systemvery similar to ToF-DOT. The main difference is in this par-ticular novel application of ToF-DOT, the imaging subjectis physically separated from the scattering media, and theircomputational model accounts for this geometric change.The principal limitation of DOT and ToF-DOT remainsthe limited spatial resolution provided by this approach.Existing DOT and ToF-DOT systems [4], [21], have only beenable to demonstrate cm -scale spatial resolution. Reconstruction Algorithms.
DOT reconstruction ap-proaches have traditionally focused on iteratively solvingapproximations of the Radiative Transfer Equation. Analyti-cal solutions to the radiative transfer equation only exist forthe most simple examples, and for any scenario approachingreal-world complexity, numerical techniques are the only
Fig. 2:
ToF-DOT concept. (a) Photon trajectories for 2 source-detector pairs. A and B are sources, α and β are detectors, andP, Q, R, and S are voxels of interest. Source-detector pair A- α ismore sensitive to P and Q and source-detector pair B- β is moresensitive to R and S. (b) Photon arrival times passing throughspecific voxels associated with source-detector pair A- α and B- β . Photons passing through voxels closer to the surface (P andR) tend to arrive earlier than photons passing through voxelsdeeper inside (Q and S). alternative. This numerical process for solving the radiativetransfer equation is computationally challenging, resultingin reconstruction algorithms that take hours to converge[22]. Fortunately, photon propagation can be reformulatedas a linear system using the Born approximation. Thensolving for the optical properties can just be a linear inverseproblem, thereby speeding up reconstruction algorithms[23], [24]. However, these algorithms still require storageof an extremely large sensitivity matrix and therefore sufferfrom increases in dimensionality of the measurements. Insummary, even the fast reconstruction techniques such as[25], [26] typically end up taking several 10s of seconds tominutes per iteration.This computational challenge is further exacerbated inToF-DOT where the inclusion of transient information addsan additional dimension to the problem. As a result, naiveattempts at high-resolution reconstruction for such ToF-DOT systems, by directly incorporating time of travel infor-mation within the existing DOT algorithms, can lead to fargreater reconstruction times as a result of dimensionality.As a result, there exist no real-time (or near real-time)reconstruction algorithms for ToF-DOT systems. O F-DOT
A traditional DOT system consists of an array of lightsources and an array of detectors placed on top of theimaging volume. Shown in Fig. 2(a) is a statistical distribu-tion of light transport paths between two different source-detector pairs. Intuitively, each of these intensity light trans-port measurements can be thought of containing weightedinformation about the attenuation (absorption) or emission(fluorescence) from the voxels in the imaging volume. Theweights themselves can be intuitively thought of as beingapproximately proportional to the likelihood that light pathsfor that source-detector pair traverse through that particularvoxel. In the example shown in Fig. 2(a), intensity lighttransport measurement between source A and detector α contains more information about voxels P and Q , whileintensity light transport measurement between source B and detector β contains more information about voxels R and S . Fig. 3:
Overview of DOT forward model.
In the linear forward model, the target scene ( µ ) is mapped to a set of measurements( m ) by the Jacobian matrix ( J ) In ToF-DOT, the light sources are typically ultra-shortpulsed sources, and the detectors measure transient (or timeof travel) information in addition to the intensity. Thus, foreach source-detector pair, the transient light transport infor-mation is recorded. Shown in Fig. 2(b), are a statistical dis-tribution of light transport paths between a source-detectorpair, where the time of travel of these paths are also color-coded. Clearly, the original intuition behind DOT holds true.But in addition to that, we notice that photons with differenttravel times pass through very different locations within theimaging volume, providing us an additional informationabout spatial localization. In the example shown in Fig. 2(b),transient light transport measurement with a travel timeof ≈ τ , contains more information about voxel R , whiletransient light transport measurement with a travel time of ≈ τ , contains more information about voxel S .The primary advantage of ToF-DOT is that this addi-tional transient information has the potential to significantlyimprove spatial resolution in the reconstructions. The propagation of light though a scattering media is well-modeled using the radiative transfer equation (RTE) [8]: ∂L ( (cid:126)r, ˆ s, t ) /c∂t = − ˆ s · ∇ L ( (cid:126)r, ˆ s, t ) − µ t L ( (cid:126)r, ˆ s, t )+ µ s (cid:90) π L ( (cid:126)r, ˆ s (cid:48) , t ) P (ˆ s (cid:48) · ˆ s ) d Ω (cid:48) + S ( (cid:126)r, ˆ s, t ) (1)Where L ( (cid:126)r, ˆ s, t ) is the radiance at a particular position (cid:126)r ,solid angle ˆ s , and time t ; P ( · ) is the phase function, whichdescribes the scattering angle; S ( · ) is the source term; and µ t is the extinction coefficient. As shown in [27], using the BornApproximation, the RTE can be reformulated as a linearequation by considering the differential measurements: m = Jµ, (2)where, µ represents the spatially varying material proper-ties within the imaging volume, m is the transient lighttransport measurements acquired by ToF-DOT, and J is theJacobian, or sensitivity matrix. Going back to our intuition,the Jacobian, J , nominally represents the weights of eachvoxel in the volume to each measurement. An overview ofthis linear model is shown in Fig. 3.Human tissue and most other biological tissues (includ-ing skull for example) are predominately scattering andhave little absorption. So it is reasonable to assume thatnative tissue absorption can be ignored. In addition, tissueoptical properties are fairly uniform, with some significantheterogenities that correspond to physiologically important variations. So, these properties are modeled as the sum-mation of a spatially homogeneous background materialcoefficient ( µ ) and a foreground, spatially varying materialcoefficient that is typically the imaging property of interest( µ ). In Equation (2), µ represents the spatial distribution ofthese heterogeneities in the scattering media. These hetero-geneities can be fluorophores (emission signal) or opticalabsorbers (absorption). In a biological context, they canrepresent features of interest such as tumors, vasculatureor regions of biological activity.If the imaging volume is discretized into N voxels = L × W × H voxels, then µ is a vector of length N voxels , thatrepresents the tissue heterogenities. Let us assume a ToF-DOT system consists of N s sources of light, N d detectors,wherein each detector measures a transient that is thenbinned into one of N t time bins. In this case, the set of allmeasurements can be represented as a vector m of length N meas = N s × N d × N t . The two quantities µ and m arerelated by the Jacobian, J , which is a matrix of dimension N meas × N voxels , where J pq = ∂m p ∂µ q . Each entry of theJacobian, J pq defines the sensitivity of measurement m p isto a corresponding heterogeneity µ q . In order to leverage the linear approximation in Equation(2), one needs to first obtain an accurate estimate of thesensitivity matrix J . In practice there are two potential waysto estimate the sensitivity matrix: (a) fast analytical approx-imation, or (b) accurate but slow Monte-Carlo simulation.Note that in either case, the computation of the sensitivitymatrix is a one time process for any application and neednot be real-time. Analytical approximation.
Using the diffusion approxima-tion, we can derive a closed form approximation to the RTE[8], [11], [28]. According to [11], we can derive this equationusing the Born Approximation: m ( (cid:126)r d , (cid:126)r s ) = (cid:90) v (cid:18) Φ ( (cid:126)r v − (cid:126)r s ) R ( (cid:126)r d − (cid:126)r v ) (cid:19) µ ( (cid:126)r v ) d(cid:126)r v (3)Where (cid:126)r s , (cid:126)r d , (cid:126)r v are the positions of source s , detector d ,and voxel v respectively; m ( (cid:126)r d , (cid:126)r s ) is the measurement asa function of source-detector position; µ ( (cid:126)r v ) is the spatialdistribution of optical properties, i.e. the image of interest; Φ ( (cid:126)r v − (cid:126)r s ) and R ( (cid:126)r d − (cid:126)r v ) are the fluence rate and diffusereflectance terms. This product is the Jacobian: J ( (cid:126)r s , (cid:126)r d , (cid:126)r v ) = Φ ( (cid:126)r v − (cid:126)r s ) R ( (cid:126)r d − (cid:126)r v ) (4)Equation (4) can be adapted to time of flight measurementsby calculating the time-domain convolution of the Green’s function and reflectance rather than the direct product asshown by Hyde et al. [28]: J ( (cid:126)r s , (cid:126)r d , (cid:126)r v , t ) = Φ ( (cid:126)r v − (cid:126)r s , t ) (cid:126) t R ( (cid:126)r d − (cid:126)r v , t ) (5)In the time-domain, Equation (3) becomes: m ( (cid:126)r d , (cid:126)r s , t ) = (cid:90) v J ( (cid:126)r s , (cid:126)r d , (cid:126)r v , t ) µ ( (cid:126)r v ) d(cid:126)r v (6)Because this expression only requires a 1D convolution,it can be used to quickly calculate the Jacobian matrix.However, this expression can only be applied to simplescene geometries, such as a single homogeneous slab, andassumes that the scene is highly scattering [8]. As a con-sequence the approximation is not appropriate in manysituations such as (a) near surface, where we are not yet inthe diffuse regime, (b) inhomogeneous tissue or (c) layeredtissue that contain low-scattering regions (such as skull,cerebrospinal fluid and brain) [29]. Monte Carlo simulations of the forward model.
Whilethe closed form expressions can be calculated efficiently,they can be limited by the prior assumptions of a highlyscattering media, and slab geometry. Since our long-termgoal is to tackle brain imaging, we primarily use MonteCarlo simulations for determining the sensitivity matrix.This technique is widely regarded as the ”gold-standard”for modeling photon propagation [27]. In Monte Carlo,simulated photons are propagated through the imagingvolume. Each photon follows a random walk, which issampled from a distribution that is parameterized by theoptical parameters of the scene [30]. Finally, the aggregateinformation from many photon samples can be used to esti-mate the sensitivity matrix. More details on this procedurecan be found in Yao et al. [27].
The goal of DOT imaging systems is to produce an imagereconstruction of the spatial distribution of optical param-eters, typically the absorption coefficient, represented by µ a . This image reconstruction is done using the followingoptimization setup: min µ (cid:107) m − f ( µ ) (cid:107) + Λ( µ ) (7)Where, µ is the spatial distribution of optical parameters; m is the set of collected measurements, which describesthe intensity of light incident on the detectors; f ( · ) is theforward model, which calculates the measured intensity asa function of the optical parameters of the scene and Λ( µ ) isan appropriately chosen regularization term.Using the linear model, the image reconstruction prob-lem can be formulated as a linear inverse problem: ˆ µ = min µ (cid:107) m − Jµ (cid:107) + (cid:107) µ (cid:107) , (8)where (cid:107) m − Jµ (cid:107) is the data fidelity term, and (cid:107) µ (cid:107) is a reg-ularization term that enforces sparsity in the heterogeneityof the optical properties within the imaging volume. Thisoptimization problem is known to be convex and there are ahost of well-understood algorithms that can be used to solveit. We use the fast iterative shrinkage thresholding algorithm(FISTA) [31] to solve this optimization since it is fast, hasreasonable memory complexity and has been shown to beaccurate (and reaches the global optimal solution). Even with the use of a fast, iterative algorithm and animplementation on a multi-CPU, multi-core computationalsystem, the algorithm remains too slow to enable real-timeapplications. As an example, if we consider reconstructionof a mm × mm × mm volume at 1 mm voxel size,using a ToF-DOT system that consists of 100 sources and 100detectors and each transient being binned into 50 differenttime bins, then the corresponding sensitivity matrix J is ofsize k × k and each FISTA iteration on a Intel Xeonmachine, with 6 cores takes about 6.3 seconds. Accuratereconstruction may require hundreds of iterations for con-vergence, meaning that total reconstruction time could beon the order of an hour. ONFOCALITY AND M ULTIPLEXING IN T O F-DOT
The computational complexity of current generation ToF-DOT reconstruction algorithms preclude near real-time op-eration. A careful study of the computational complexityprovides two potential avenues that might facilitate signifi-cant improvements in computational speed.
Measurement selection.
The computational complexity ofsolving large scale linear inverse problems scales betweenquadratic and cubic in the problem size, based on the kind ofalgorithms used. This means that, in practice, while the 10 × -100 × increased measurements afforded by ToF-DOT signif-icantly improves spatial resolution of the reconstruction, italso slows down the reconstruction time by the several or-ders of magnitude compared to traditional DOT algorithms.One way to combat this is measurement selection, whereinonly a select subset of measurements are used in the recon-struction. To maintain the resolution advantages providedby ToF-DOT, one has to carefully select the measurementsso as to ensure that the maximally useful (high SNR, highinformation gain) measurements are retained. Faster forward models.
The key computational step inalmost all iterative algorithms (including FISTA) that areintended to solve the optimization problem in Eqn. (7) isthe repeated application of the forward operator (or its con-jugate or transpose). In the case of ToF-DOT, this amountsto a matrix multiplication with the corresponding sensitivitymatrix (or its transpose) and this matrix multiplication haslinear complexity in the number of elements in the matrix(or quadratic in the number of rows/columns). One key ideathat has been in many other applications is if under somerestricted regimes of operation, the general linear model canbe reduced to a convolutional form, then one could leveragefast implementations of convolutions (that rely on FFTs) tosignificantly reduce the computational burden.Here we argue that focusing on confocal ToF-DOT mea-surements allows us to leverage both these advantagessimultaneously, allowing us to achieve, real-time ToF-DOTreconstruction performance. This would correspond to re-taining all the measurements wherein the source and thedetector location are the same (or close enough to be mod-eled as confocal in a real system).
Related work.
The scanning-time and reconstruction-timechallenge in ToF-DOT is not unique to DOT, but rather com-mon across a variety of emerging applications that attemptto utilize the extra temporal dimension offered by transientdetectors such as SPADs. These applications include imag-ing around corners [32], non-line-of-sight imaging [33], and
Fig. 4:
Validity of convolutional approximation.
Visualizationof (a) rows of Jacobian for different absorber locations, and theircorresponding (b) 1D profiles along X and Y directions (coloredlines). Note: 1D profiles have been aligned for visualization. imaging through thick diffusers [20] and in all of theseexamples, the imaging geometry is somewhat similar toToF-DOT. There is an array of sources and detectors that arescanned and transient light-transport measurements are ob-tained. The principal difference between these applicationsand ToF-DOT is that in these applications, the scatteringsurface or the thick diffuser acts as an obscurant and thegoal is to image objects beyond that obscurant. In contrast,in ToF-DOT, the goal is to obtain a volumetric image of theoptical properties (scattering, absorption or fluorescence) ofthe tissue itself. Thus the computational model for lightpropagation in these different applications are quite differ-ent. That said, the symmetry in imaging geometry, betweenthese applications and ToF-DOT is quite striking.In all of these applications that use transients, recon-struction algorithms tended to be slow precluding any real-time operation. Over the last few years, confocality hasemerged as a key idea enabling real-time reconstructionin these applications. First, within non-line-of-sight recon-struction, it was shown that restricting the measurementsto confocal measurements allows both a reduction in thenumber of measurements and also enabled a convolutionalapproximation to the forward model resulting in real-timereconstruction algorithms [34], [35]. More recently, similarinsight was used to demonstrate near-real-time reconstruc-tion performance for imaging through thick obscurants [20].We are motivated by the success of these techniques andshow that this idea, when translated to ToF-DOT, allows usto obtain real time ToF-DOT reconstructions for estimating2/3D optical properties of thick tissues.
It is known that collocated source-detector pair contains themost information as it pertains to deep features [16]. Ona high-level, this claim is based on the intuition that thecollocated source-detector pair will possess greater sensi-tivity to deeper features than when a larger source-detectorseparation is used. This idea is later reinforced by our resultsin Fig. 6. We see that selecting only the measurements fromthe collocated source-detector pair leads to a more well-conditioned Jacobian matrix than selecting measurementsfrom source-detector pairs of arbitrary separation distance.
Convolutional approximation.
From equation (2) we seethat the forward model for ToF-DOT can be modeled as a linear system. Fortunately, when restricting our atten-tion to confocal measurements, the linear operator is shift-invariant. This shift-invariance allows us to develop a con-volutional approximation for the confocal ToF-DOT system.To empirically demonstrate the shift invariance of thesensitivity matrix (i.e., the matrix is doubly circulant), weuse the Monte Carlo simulator to generate different rowsof the sensitivity matrix that correspond to point targets atdifferent locations within the volume. In this simulation,we assume a confocal geometry with features fixed to aspecific depth. Fig. 4(a) shows a visualization of rows ofthe Jacobian. Additionally, from Fig. 4(b) we see that eachblur kernel has the same profile. As the feature locationis shifted, there is a corresponding shift in the measuredoutput. This indicates that the Jacobian is a doubly circulantmatrix. Therefore, when performing image reconstructionusing confocal measurements, we can apply the forwardmodel using a convolutional approximation rather than amatrix-vector product. Equation (2) can be substituted with: m ( x, y, t ) = ρ ( x, y, t ) (cid:126) µ ( x, y ) . (9)Here m ( x, y, t ) is the measurement, which is now a functionof the collocated source-detector ( x, y ) and time t ; µ ( x, y ) is the spatially varying material properties, which is nowa function of just the lateral positions ( x, y ) ; and finally ρ ( x, y, t ) is the blur kernel. The blur kernel ρ ( x, y, t ) can bedetermined using either the Monte Carlo simulator or theanalytical expressions by calculating (or estimating in thecase of Monte Carlo) the measurement for a single feature,i.e. a spatial delta function. Computational Complexity Analysis.
Using the standardforward model, the main bottleneck in solving the inverseproblem is the matrix-vector product: Jµ . This operationscales linearly with the number of sources ( N s ), numberof detectors ( N d ), time bins( N t ), and number of voxels( N voxels ). The runtime complexity is O ( N s N d N t N voxels ) .The bottleneck for memory usage is the storage of theJacobian, which is of complexity equal to the matrix size.This complexity can be significantly reduced in the con-focal mode, using a convolutional model. Convolution witha size K × K blur kernel can be efficiently implementedusing the fast fourier transform (FFT). In this case, thecomputational complexity is O ( N s N t K log( K )) . The firstimprovement in computational complexity is the reductionin the number of measurements from N s N d N t to N s N t (in confocal measurements N s = N d and only 1 transientmeasurement is obtained per source location). The secondimprovement arises because of the convolutional approxi-mation. In addition, the size of the convolutional kernal K istypically much smaller than the field of view of the volumebeing images as well resulting in additional efficiencies.Figure 7 shows the significant reduction in computa-tional complexity that is achieved due to the convolutionalmodel imposed on the confocal ToF-DOT measurements.There is two orders of magnitude speed-up in runtimecompared to existing ToF-DOT algorithms [28]. Even moreremarkable is the resultant confocal ToF-DOT algorithm isover an order of magnitude more efficient than even conven-tional DOT algorithms [8] that do not utilize any transientinformation at all (and result in worse spatial resolution). Fig. 5:
Experimental setup to test CToF-DOT. (a) Future goal to develop a wearable array of sources and detectors. This array isemulated in our testbed (b) by raster scanning a laser beam and single pixel detector. (c) shows an image of the physical setup,with the SPAD (white), galvo mirrors (red), E-ink display (orange) and tissue phantom (blue).
Traditional DOT systems use point-scanning to capturemeasurements, which can result in a long measurementcapture durations precluding the capture of dynamics. Thischallenge is compounded by the fact that DOT systems oftenrequire a long exposure duration (even for a single sourcelocation), due to the fact that only a miniscule fraction of theincident photons are sensed at the detector – meaning thatthe detectors are operating at extremely low photons levels.We demonstrate that source multiplexing can be used topotentially address both these challenges simultaneously.
Multiplexing sources far away.
Typical DOT and ToF-DOTsystems have a field of view of the order of − cmon a side to image through skull. Detectors and sourcesare typically placed on an array (anywhere from × to × arrays) with a spacing of a few mm to a cm between array elements. When a source is ’on’, all detectorsare measuring the corresponding light transport transients,but the detectors that are far away typically (i.e., with safeillumination power and within reasonable exposure dura-tions) get little to no photons making their measurementsuseless. In practice, the photon signal dies exponentiallywith separation distance and after about a − cm sep-aration there are typically very few photons measured.This means that one can safely assume that there is nocross-talk between measurements even if multiple sourcesare kept ’on’ simultaneously, as long as we can ensuresufficient separation between the sources. For each detectormeasurement, we can allocate the entire transient measure-ment to the closest source (note that this is only possiblewhen we can ensure that sources that are simultaneously’on’ are sufficiently far away). In our prototype system witha FOV of about cm, this means that we can multiplex upto4 sources simultaneously without any cross-talk. This allowsus to get a 4 × improvement in total capture time, while itdoes not affect the SNR of the individual measurements. Multiplexing sources with cross-talk.
Even in the presenceof measurement cross-talk, one can obtain significant multi-plexing gain [36], [37], [38], [39]. It is well-known that thismeasurement gain is somewhat limited at the high signallevel regimes but becomes significant in photon starved en-vironments and applications such as ToF-DOT. With practi-cal constraints on illumination intensity (set by safety limits)and detector exposure duration, we typically measure a fewthousand photons per entire transient — resulting in tens to hundreds of photons per time bin. At such extremely lowsignal levels, it is expected that source multiplexing (withappropriate post-capture de-multiplexing), will result in asignificant gain. Shown in Figure 10, is the SNR improve-ment that is obtained due to multiplexing sources.
Composite reconstruction algorithm.
In the presence ofsource multiplexing the new measurements acquires y be-come multiplexed versions of the old measurements m –wherein y and m are related by the multiplexing matrix S as y = Sm = SJµ . The combined optimization problem tobe solved becomes µ = min µ (cid:107) y − SJµ (cid:107) + (cid:107) µ (cid:107) , (10)where Jµ can be further efficiently implemented withineach iteration using the convolutional approximation. Asbefore, we use the fast iterative shrinkage thresholdingalgorithm (FISTA) [31] to solve this optimization problem. ATERIALS AND M ETHODS
Simulation setup.
We use an in-house Monte Carlo simula-tor for generating measurements and the Jacobian matricesneeded for both simulated and experimental results. Ourimplementation is based on the standard Monte Carlo forscattering samples and closely follows the details in [27],[30]. The Monte Carlo simulations are run on GPUs (NvidiaRTX 2080 Ti). Obtaining the Jacobian through analyticalexpressions was performed on CPU (Intel Xeon 3.30 GHz).Finally, our simulator can operate in both fluorescenceand absorption imaging mode. The details of extendingabsorption-based Monte Carlo to fluorescence mode aredescribed by Liebert et al. and Chen et al. [40], [41] and wefollow these to adapt our implementations as well.
Experimental setup.
To perform real-world data collection,we constructed an experimental prototype as shown in Fig.5. Two galvo mirrors raster scan the source and detectorseparately, emulating measurements that could be obtainedwith an array of light sources and detectors. A NKT Pho-tonics SuperK EXTREME supercontinuum laser produceseither 680 nm or 480 nm , 80 MHz light pulses for absorptionand fluorescence experiments, respectively. Photon arrivaltimes are detected using a MPD FastGatedSPAD single pixeldetector with a temporal jitter of < ps connected to aPicoQuant HydraHarp 400. A MPD Picosecond Delayerprovides a delay to the synchronization signal from the Fig. 6:
Jacobian matrix conditioning.
The singular values ofthe Jacobian matrix are plotted to determine the matrix con-ditioning. We compare traditional DOT (blue), ToF-DOT (red),and our CToF-DOT (yellow). We see that the introduction oftime binning (ToF-DOT) and confocal geometry (CToF-DOT)provides improvements to our matrix conditioning. laser to ensure the SPAD’s 5 ns gate encompasses the entiretransient from the scene. Scattering tissue phantoms.
We use a 3D printer (FormlabsForm 3) to synthesize the optical tissue phantoms used inour experiments. Our goal is to emulate a skull phantomand we closely mimic the known properties of the humanskull including its thickness and mean scattering length.The scattering slab is 50 mm × mm × mm with ascattering coefficient µ s = 9 mm − , corresponding to ∼ mm (correspond-ing to 45 MFPs) is used. Both the thickness and scatteringcoefficient of this skull phantom were set to be within theaccepted range for human skull [3], [42]. We adapt theprocedure used by Dempsey et al. and synthesize our ownresin for optical phantoms [43]. The scattering parameters ofthe phantom are set by controlling the volume ratio of the’white’ and ’clear’ Form resins. The scattering coefficient ofthe phantom can be determined by measuring the temporalbroadening of the transients [44]. In Fig. 9, we see thatour experimentally measured transients matches with theoutput of Monte Carlo simulations. The surface curvature ofthe human skull is something our skull phantom does notemulate, but we do not believe this has a significant effect onthe resolution or performance characteristics predicted byour phantoms. The results shown in Figures 9, 10, 11, and12 use this skull phantom as the scattering layer betweenthe target and the imaging system. Absorptive and fluorescent targets.
To emulate an absorp-tive target such as a tumor, we use a E-ink display behindthe scattering tissue sample. An E-ink display allows us toprogrammatically control the spatially varying absorption ata fine spatial resolution. The results shown in Figures 9, 10,and 12 use the E-ink display-based target behind the skullphantom. In order to emulate fluorescent targets, we embedfluorescent beads (Fluoresbrite YG 1 µ m beads) in PDMS.The spatial patterning of the fluorescent target is achievedusing a 3D printed mold. The results shown in Fig. 11 usethe fluorescent target behind the skull phantom. Fig. 7:
Algorithm runtime characterization.
The algorithm run-time was characterized as a function of source-detector arraysize (a), and the voxel grid size (b). We see almost two orders ofmagnitude decrease in runtime using our methods as comparedto traditional DOT [8] and ToF-DOT [28].Fig. 8:
Simulated spatial resolution of CTOF-DoT.
Our tech-nique is able to resolve two 0.5 mm thick lines separated by 0.5 mm . Traditional DOT fails at this task. Simulated measurementswere generated in Monte Carlo. ESULTS
We perform an extensive array of experiments, and perfor-mance characterizations both in simulation and experimen-tally using a benchtop prototype ToF-DOT system.
Inverting the Jacobian matrix is critical to our image recon-struction procedure. A well-conditioned Jacobian will allowus to improve our image reconstruction quality. As shown inFig. 6 we demonstrate that the additional information pro-vided by time-binning results in a more well-conditionedmatrix. Each Jacobian was obtained through Monte Carlosimulations. We compare three cases: 1)
Traditional DOT in which all measurements are a scalar intensity value;2)
ToF-DOT , which uses all time bins; and 3)
CToF-DOT .625 total scan points and 65 time bins were used for eachJacobian. The simulated scene was a × × grid of 1 mm cubes. All singular value plots were normalized to 1.Below a threshold − , the singular values are consideredto be below the noise floor. We see that the introduction oftime domain information improves the matrix conditioning,increasing the minimum singular value from 67 to 82. Inaddition, for the confocal geometry, because all 625 mea-surements were collocated, the minimum singular valuewas further increased to 1276. This provides additional Fig. 9:
Resolution test with experimental data.
Panels from left to right show scattering tissue phantom, imaging target, and our1D reconstructions of the imaging target (either 0.5, 1, or 2 mm spacing and linewidth). The inset image on the far left showsthat the experimentally captured TPSF matches the results predicted by Monte Carlo, thus verifying the scattering coefficient of µ s = 9mm − of the 6.5 mm skull phantom ( ∼
60 MFPs). The 1D reconstruction demonstrates that CToF-DOT is capable of mm resolution.Fig. 10: Simulated and Experimental Multiplexing Results.
Multiplexing allows comparable performance with reducedintegration time compared with single point scanning for CToF-DOT imaging. (Left) Plots shows PSNR versus integration timefor simulated and experimental results. The images correspond to the image reconstructions performed at different integrationtimes with/without multiplexing. With multiplexing, the image reconstruction is more robust to noise at lower integration times.Measurements were captured through a 5 mm phantom with µ s = 9mm − ( ∼
45 MFPs). support that collocated source-detector pairs provide moreinformation than an arbitrary set of source-detector pairs.
Additionally we test the algorithm runtime. These exper-iments were conducted on an Intel Xeon 3.30 GHz CPU.We test how the image reconstruction speed is affected by2 system parameters: the voxel size (for a fixed total area),and the number of sources and detectors. We compared tothe algorithms for traditional DOT and ToF-DOT, whichwere constructed in-house (described by Wang et al. andHyde et al. , respectively [8], [28]). In Fig. 7, we see that theconfocal geometry achieves almost 2 orders of magnitudeimprovements in speed primarily attributed to a reductionin Jacobian matrix size.
Simulation resolution test.
We test the spatial resolutionthat can be achieved by traditional DOT and our method(Fig. 8). The simulated scene consists of two fluorescentlines, with 0.5 mm line width and separation. The featureswere 6.5 mm deep and the background scattering coefficientwas set to µ s = 9 . − . With a × confocal scan, wesee that we are able to clearly resolve the two lines, whichindicates our system can resolve mm -scale features. Experimental resolution test.
In Fig. 9, we performed a 1-dimensional resolution test through a 6.5 mm thick skullphantom ( µ s = 9mm − ) by scanning 32 points in a confocalgeometry. To obtain the Jacobian experimentally, a black linewas projected on the E-ink display at 32 locations with 0.5,1, and 2 mm separation. For each line position, a 32-pointscan was captured, which becomes a column of the Jacobianmatrix. After obtaining the Jacobian, we projected the targetimage onto the E-ink display: two lines of thickness andseparation distance between 0.5-2 mm. Though there is aslight offset due to calibration, we are able to resolve thetwo lines and demonstrate mm -scale spatial resolution (Fig.9). Simulations on multiplexing.
We tested source multiplex-ing with an × array of sources, which leads to multiplex-ing with a × Hadamard matrix. The simulated mea-surements and Jacobian were generated using the analyticalexpressions. Poisson noise was applied, assuming a countrate of 5 million counts per second, the approximate inten-sity level before our SPAD experiences the pile-up effect. Inaddition, dark count noise was added to our measurements,with a rate parameter of 200 counts/sec corresponding to Fig. 11:
Fluorescence imaging with CToF-DOT. (a) Imagereconstruction of fluorescence targets (2 4 mm lines separatedby 4 mm and oval-line scene) using CToF-DOT. (b) For the oval-line image reconstruction, we show that the image cannot berecovered without time binning. the dark count rate of the FastGatedSPAD. Image recon-struction of the letter ’R’ is performed for a range of inte-gration times. In Fig. 10, we see that with multiplexing, theimage reconstruction still maintains a reasonable PSNR atshort exposure durations. This demonstrates that the imagereconstruction with multiplexing is more robust to increasednoise at lower integration times. The panels on the rightof Fig. 10 show the reconstruction results, again showingincreased robustness to noise with multiplexing. Multiplexing with experimental data.
In addition to sim-ulations, we also captured experimental data to test thebenefits of source multiplexing. The plots on Fig. 10 showthe PSNR as a function of the integration time. We show thatthe multiplexed measurements are more robust to highernoise levels at lower integration times. Therefore, multi-plexing can improve the temporal resolution by reducingthe integration time needed to maintain a threshold imagereconstruction quality. From the experimental results, wesee an order of magnitude improvement since the imagereconstruction quality at 10 ms is approximately comparableto 100 ms integration time without multiplexing (Fig. 10).We must also account for a factor of since multiplexingdoubles the total number of measurements, resulting in anoverall ∼ In the image reconstruction experiments with real-worlddata, the scattering slab is placed on top of either the E-ink display or the fluorescent target to emulate scatteringby biological tissue. In Fig. 11(a). we see that the signalfrom the fluorescent targets is significantly enhanced by re-moving early arriving photons using time-gating, albeit stillblurred due to the effects of scattering. The target imagescan be recovered by our image reconstruction algorithm.Fig. 11(b) shows measurements and a reconstruction of theoval-line scene without time-binning. The measurements aredominated by background noise from the excitation lightleading to an image reconstruction which contains virtuallyno information about the underlying scene. The oval linescene experiment was conducted with the standard 6.5 mmphantom with µ s = 9mm − , without an emission filter.However, for imaging the two lines image reconstruction,a filter was used as well as an un-calibrated skull phantomwith scattering coefficient in the range − − .Additionally, we demonstrate the benefits of using a con-focal geometry in Fig. 12(a). Here, we image absorber targets Fig. 12:
Experimental image reconstruction of absorptive tar-gets. (a) Comparison of ToF-DOT and CToF-DOT reconstruc-tion of 2D absorption targets through a 6.5 mm phantom with µ s = 9mm − ( ∼
60 MFPs). We maintain comparable imagereconstruction quality while reducing the computation time forthe inverse solver by approximately two orders of magnitude,and reducing the number of scan points by almost an order ofmagnitude. (b) Depth reconstruction of a 3D absorptive target.The triangle is located at 1 mm ( ∼
10 MFPs) depth and the circleis located 3 mm ( ∼
30 MFPs) depth with µ s = 9mm − . Ourmethod is able to accurately determine the depth and shapes ofthe absorbing targets. displayed with the E-ink display. In the standard TD-DOTset up, we capture a measurement for all pairs of sourcesand detectors. For this experiment, we use a × arrayof sources and detectors, which correspond to 10,000 totalscan points. Additionally, we capture measurements for thesame scene using a × array in the confocal geometry,corresponding to 1024 scan points. As shown in Fig. 12(a),even though the number of scan points is reduced by almostan order of magnitude, we are able to maintain comparableimage reconstruction quality using confocal ToF-DOT com-pared to ToF-DOT. Additionally, with the confocal geometry,the algorithm runtime is reduced by approximately twoorders of magnitude from ∼ ∼ ∼ ONCLUSIONS
We demonstrate that confocal and multiplexed versions ofToF-DOT have the potential to achieve millimeter resolu-tion, real-time imaging through thick scattering tissue. Withfuture developments in terms of on-chip SPAD hardwareand integrated source-detector arrays, these results can leadto wearable imaging devices paving the way for high reso-lution structural and functional imaging of the brain.
CKNOWLEDGMENTS
The authors would like to thank Prof. Aswin Sankara-narayanan for his suggestions and edits for our paper. Theauthors would like to acknowledge Biorender in providingcertain images to supplement figures. Yongyi Zhao sup-ported by a training fellowship from the NLM TrainingProgram (T15LM007093). R EFERENCES [1] A. K. Pediredla, S. Zhang, B. Avants, F. Ye, S. Nagayama, Z. Chen,C. Kemere, J. T. Robinson, and A. Veeraraghavan, “Deep imaging in scattering media with selective plane illumination microscopy,”Journal of Biomedical Optics, vol. 21, no. 12, pp. 1–14, 2016.[2] B.-H. Oh, K. H. Kim, and K.-Y. Chung, “Skin imaging usingultrasound imaging, optical coherence tomography, confocal mi-croscopy, and two-photon microscopy in cutaneous oncology,”Frontiers in Medicine, vol. 6, p. 274, 2019.[3] F. Bevilacqua, D. Piguet, P. Marquet, J. D. Gross, B. J. Tromberg,and C. Depeursinge, “In vivo local determination of tissue opticalproperties: applications to human brain,” Applied Optics, vol. 38,no. 22, pp. 4939–4950, 1999.[4] D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kilmer,R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuseoptical tomography,” IEEE signal processing magazine, vol. 18,no. 6, pp. 57–75, 2001.[5] J. Xia, J. Yao, and L. V. Wang, “Photoacoustic tomography: prin-ciples and advances,” Electromagnetic waves (Cambridge, Mass.),vol. 147, pp. 1–22, 2014.[6] A. Puszka, L. Di Sieno, A. Dalla Mora, A. Pifferi, D. Con-tini, G. Boso, A. Tosi, L. Herv´e, A. Planat-Chr´etien, A. Koeniget al., “Time-resolved diffuse optical tomography using fast-gated single-photon avalanche diodes,” Biomedical optics express,vol. 4, no. 8, pp. 1351–1365, 2013.[7] L. Nie, X. Cai, K. Maslov, A. Garcia-Uribe, M. A. Anastasio, andL. V. Wang, “Photoacoustic tomography through a whole adulthuman skull with a photon recycler,” Journal of biomedical optics,vol. 17, no. 11, p. 110506, nov 2012.[8] L. Wang and H.-I. Wu, Biomedical Optics: Principles and Imaging.Hoboken, NJ: John Wiley and Sons, 2007.[9] M. Kempe, W. Rudolph, and E. Welsch, “Comparative study ofconfocal and heterodyne microscopy for imaging through scatter-ing media,” JOSA A, vol. 13, no. 1, pp. 46–52, 1996.[10] E. A. Sergeeva, “Scattering effect on the imaging depth lim-itin two-photon fluorescence microscopy,” Quantum Electronics,vol. 40, no. 5, p. 411, 2010.[11] C. Liu, A. Maity, A. W. Dubrawski, A. Sabharwal, and S. G.Narasimhan, “High resolution diffuse optical tomography usingshort range indirect subsurface imaging,” in IEEE InternationalConference on Computational Photography. IEEE, May 2020.[12] H. Zhao and R. J. Cooper, “Review of recent progress towarda fiberless, whole-scalp diffuse optical tomography system,”Neurophotonics, vol. 5, no. 1, p. 011012, 2017.[13] E. M. Frijia, A. Billing, S. Lloyd-Fox, E. V. Rosas, L. Collins-Jones, M. M. Crespo-Llado, M. P. Amad´o, T. Austin, A. Edwards,L. Dunne et al., “Functional imaging of the developing brain withwearable high-density diffuse optical tomography: a new bench-mark for infant neuroimaging outside the scanner environment,”NeuroImage, p. 117490, 2020.[14] M. D. Wheelock, J. P. Culver, and A. T. Eggebrecht, “High-densitydiffuse optical tomography for imaging human brain function,”Review of Scientific Instruments, vol. 90, no. 5, p. 051101, 2019.[15] A. Gibson, J. Hebden, and S. R. Arridge, “Recent advances indiffuse optical imaging,” Physics in Medicine & Biology, vol. 50,no. 4, p. R1, 2005.[16] A. Pifferi, D. Contini, A. Dalla Mora, A. Farina, L. Spinelli, andA. Torricelli, “New frontiers in time-domain diffuse optics, areview,” Journal of biomedical optics, vol. 21, no. 9, p. 091310,2016.[17] A. Lyons, F. Tonolini, A. Boccolini, A. Repetti, R. Henderson,Y. Wiaux, and D. Faccio, “Computational time-of-flight diffuseoptical tomography,” Nature Photonics, vol. 13, no. 8, pp. 575–579,2019.[18] A. Farina, S. Tagliabue, L. Di Sieno, E. Martinenghi, T. Durduran,S. Arridge, F. Martelli, A. Torricelli, A. Pifferi, and A. Dalla Mora,“Time-domain functional diffuse optical tomography systembased on fiber-free silicon photomultipliers,” Applied Sciences,vol. 7, no. 12, p. 1235, 2017.[19] L. Di Sieno, J. Nissinen, L. Hallman, E. Martinenghi, D. Contini,A. Pifferi, J. Kostamovaara, and A. D. Mora, “Miniaturized pulsedlaser source for time-domain diffuse optics routes to wearabledevices,” Journal of biomedical optics, vol. 22, no. 8, p. 085004,2017.[20] D. B. Lindell and G. Wetzstein, “Three-dimensional imagingthrough scattering media based on confocal diffuse tomography,”Nature Communications, vol. 11, no. 1, p. 4517, 2020.[21] A. Puszka, L. Di Sieno, A. Dalla Mora, A. Pifferi, D. Contini,A. Planat-Chr´etien, A. Koenig, G. Boso, A. Tosi, L. Herv´e et al.,“Spatial resolution in depth for time-resolved diffuse optical to-mography using short source-detector separations,” Biomedicaloptics express, vol. 6, no. 1, pp. 1–10, 2015. [22] H. K. Kim, M. Flexman, D. J. Yamashiro, J. J. Kandel, and A. H.Hielscher, “Pde-constrained multispectral imaging of tissue chro-mophores with the equation of radiative transfer,” Biomedicaloptics express, vol. 1, no. 3, pp. 812–824, 2010.[23] D. A. Boas, T. Gaudette, and S. R. Arridge, “Simultaneous imagingand optode calibration with diffuse optical tomography,” Opticsexpress, vol. 8, no. 5, pp. 263–270, 2001.[24] T. Tarvainen, V. Kolehmainen, J. P. Kaipio, and S. R. Arridge,“Corrections to linear methods for diffuse optical tomography us-ing approximation error modelling,” Biomedical Optics Express,vol. 1, no. 1, pp. 209–222, 2010.[25] M. A. Naser and M. J. Deen, “Time-domain diffuse optical to-mography using recursive direct method of calculating jacobianat selected temporal points,” Biomedical Physics & EngineeringExpress, vol. 1, no. 4, p. 045207, 2015.[26] M. Mozumder and T. Tarvainen, “Time-domain diffuse opticaltomography utilizing truncated fourier series approximation,”JOSA A, vol. 37, no. 2, pp. 182–191, 2020.[27] R. Yao, X. Intes, and Q. Fang, “Direct approach to compute Ja-cobians for diffuse optical tomography using perturbation MonteCarlo-based photon “replay”,” Biomedical Optics Express, vol. 9,no. 10, pp. 4588–4603, 2018.[28] D. Hyde, “Improving forward matrix generation and utilizationfor time domain difiuse optical tomography,” Ph.D. dissertation,Worcester Polytechnic Institute, 2002.[29] H. K. Kim, L. D. Montejo, J. Jia, and A. H.Hielscher, “Frequency-domain optical tomographicimage reconstruction algorithm with the simplifiedspherical harmonics (SP(3)) light propagation model.”International journal of thermal sciences = Revue generale de thermique,vol. 116, pp. 265–277, jun 2017.[30] L. Wang, L., Jacques, S. L., & Zheng, “MCML - Monte Carlomodeling of light transport in multi-layered tissues,” Computermethods and programs in biomedicine, vol. 47, no. 2, pp. 131–146,1995.[31] A. Beck and M. Teboulle, “A Fast Iterative Shrinkage-ThresholdingAlgorithm,” Society for Industrial and Applied MathematicsJournal on Imaging Sciences, vol. 2, no. 1, pp. 183–202, 2009.[32] A. Velten, T. Willwacher, O. Gupta, A. Veeraraghavan, M. G.Bawendi, and R. Raskar, “Recovering three-dimensional shapearound a corner using ultrafast time-of-flight imaging,” NatureCommunications, vol. 3, no. 1, p. 745, 2012.[33] A. Pediredla, A. Dave, and A. Veeraraghavan, “Snlos: Non-line-of-sight scanning through temporal focusing,” in 2019 IEEEInternational Conference on Computational Photography (ICCP),2019, pp. 1–13.[34] B. Ahn, A. Dave, A. Veeraraghavan, I. Gkioulekas, and A. Sankara-narayanan, “Convolutional Approximations to the General Non-Line-of-Sight Imaging Operator,” in 2019 IEEE/CVF InternationalConference on Computer Vision (ICCV), 2019, pp. 7888–7898.[35] M. O’Toole, D. B. Lindell, and G. Wetzstein, “Confocal non-line-of-sight imaging based on the light-cone transform,” Nature, vol.555, no. 7696, pp. 338–341, 2018.[36] O. Cossairt, M. Gupta, and S. K. Nayar, “When does compu-tational imaging improve performance?” IEEE Transactions onImage Processing, vol. 22, no. 2, pp. 447–458, 2013.[37] A. C. Sankaranarayanan, “Hadamard multiplexing and when it isuseful,” 2018.[38] Y. Y. Schechner, S. K. Nayar, and P. N. Belhumeur, “Multiplexingfor Optimal Lighting,” IEEE Transactions on Pattern Analysis andMachine Intelligence, vol. 29, no. 8, pp. 1339–1354, aug 2007.[39] K. Mitra, O. S. Cossairt, and A. Veeraraghavan, “A framework foranalysis of computational imaging systems: Role of signal prior,sensor noise and multiplexing,” IEEE Transactions on PatternAnalysis and Machine Intelligence, vol. 36, no. 10, pp. 1909–1921,2014.[40] A. Liebert, H. Wabnitz, N. ˙Zołek, and R. Macdonald, “Monte Carloalgorithm for efficient simulation of time-resolved fluorescence inlayered turbid media,” Optics Express, vol. 16, no. 17, pp. 13 188–13 202, 2008.[41] J. Chen, V. Venugopal, and X. Intes, “Monte Carlo based methodfor fluorescence tomographic imaging with lifetime multiplexingusing time gates,” Biomedical Optics Express, vol. 2, no. 4, pp.871–886, 2011.[42] H. Li, J. Ruan, Z. Xie, H. Wang, and W. Liu, “Investigationof the critical geometric characteristics of living human skullsutilising medical image analysis techniques,” International Journalof Vehicle Safety, vol. 2, jan 2007.2