Time-correlated single photon counting for evaluation of cavity dimensions with multibounce photons

Abstract. Cavity exploration is an important task for optical remote sensing that finds application in many fields, such as the assessment of disaster areas or the evaluation of subsurface cavities in terrestrial and extraterrestrial scenarios. In our work, we use single photon-counting avalanche diode cameras to record temporal signatures in the picosecond time domain and to exploit the impulse response of cavities. When illuminated with a short optical pulse, the photons are reflected several times inside the cavity before exiting it and being detected. In our datasets, we observe several intensity modulations correlated to specific optical paths with different numbers of reflections inside the cavity. We developed an analysis model based on simulation of the radiant transfer that provides a quick rough estimate of the cavity dimensions without costly measurements or computational reconstruction. Specifically, we developed a spherical model and a hemispherical model that can be used to estimate the cavity dimensions (diameter or height). Finally, we present experimental studies of eight different cavities with different dimensions, shapes, and surface materials. Our evaluations show that we can apply our analytical models to signatures of different cavities, which gives us an initial estimate of the cavity sizes before a detailed mapping is carried out.


Introduction
Exploration of natural and artificial cavities (for instance, caves and tunnels) is a common task in many scientific and technical fields, such as archaeology, 1,2 civil engineering, environmental monitoring, 3 geophysics, 4 mining, and mineral deposit exploration. 5 In addition, cavity detection and classification is a crucial ability in security and search and rescue scenarios, 6 such as disaster management, environmental monitoring, and law enforcement and military [7][8][9][10][11] operations.
In these terrestrial applications, cavity detection can be accomplished by microgravity and magnetic field measurements, 1,7,10 ground penetrating radar, [11][12][13][14] ultrasound, 15 acoustic 9 and seismic measurements, 4,8 or earth electrical resistivity profiling. 16,17 Cavity mapping, in turn, is often carried out by optical means, which can be used inside the cavities because of their compact dimensions and are preferred over other methods because of their high spatial resolution. These optical mapping methods include the application of laser scanners, 18 photogrammetry, 2 and stereo imaging. 3 Further, often these methods are discussed as part of autonomous robots 19 and sensing platforms. 6 Recently, there has been discussion of exploring extraterrestrial cavities for future lunar and other space missions to protect against threats, such as cosmic rays and micrometeorites. [20][21][22][23][24][25][26] Further, a satellite-based laser imaging, detection, and ranging (LIDAR) system has been proposed to remotely detect temporal signatures of lunar caves and reconstruct their shape by computational imaging. 27 Our work is related to these approaches as we also use a LIDAR-like acquisition device to record transient signatures, but in contrast, we focus on remote classification of cavities only, as illustrated in Fig. 1(a). We estimate the dimensions, i.e., the cavity diameter or height, and the overall geometry (spherical or hemispherical). We do not *Address all correspondence to Martin Laurenzis, martin.laurenzis@isl.eu intend to reconstruct or map the cavity shape. Some preliminary results were published in a conference paper. 30 We assume an application that focuses on specific features of cavities, as it is the case when exploring lava tube on the moon or searching for buried people in collapsed buildings. Therefore, our approach represents an intermediate step that occurs after the discovery of cavity openings and before insertion of mapping instruments into the cavity or the use of detailed computational imaging measurements. A fast initial classification of the cavity dimensions may reduce the number of candidates for further investigation and thus reduce the overall data acquisition and processing load. Figure 1(a) shows two different application scenarios for classification of natural cavities: (b) a lava tube on the Moon 28 and (c) a terrestrial cave in collapsed Karstic rock. 29 2 Temporal Response of Cavities

Temporal Response of Integrating Spheres
We assume a cavity as a hollow structure in which a collimated laser pulse can enter through an opening. This light is reflected multiple times (n-bounce path) with an efficiency depending on the surface reflection coefficient and eventually leaves the cavity through an opening in the direction of the sensing system. Therefore, a cavity acts similar to an integrating sphere.
Due to the multiple reflection processes, the exiting light intensity distribution is spatially and angular homogenized and temporally stretched. As given in Eq. (1), the temporal response I exit of an integrating sphere equals a convolution of the input light pulse I input and its temporal transfer function T. In Eq. (2), T is defined as an exponential decay with a time constant τ given by the diameter ∅ of the integrating sphere, the speed of light c, and the average surface reflectivity ρ ∈ ½0;1: 31 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 2 0 4 I exit ðtÞ ¼ I input ðtÞ Ã TðtÞ; (1) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 1 6 0 Due to the multireflections, a short laser pulse is stretched in a relatively long pulse. This stretching effect is increased with larger diameter ∅ (longer traveling paths) and higher reflectivity ρ (more bounces). We show hereafter that this empiric equation is valid only for pulse length equivalent to or greater than the sphere diameter. At shorter pulse lengths, the signals from different multiple reflection paths mix less and can be distinguished by high-resolution measurements as additional modulations of the temporal signal amplitude. Fig. 1 Cavity exploration (a) follows the steps of detection, classification, and mapping. 25,26 Natural cavities can be found, for example, in (b) lava tubes such as the lunar pit in the Mare Tranquillitatis crater 28 or (c) caves formed by collapsed structures in Karstic rock such as the Brillenhöhle cave (Germany). 29

Simulation of the Temporal Response by Transient Rendering
The temporal response of cavities can be simulated, e.g., by Monte Carlo simulations 32 or transient rendering, 33 taking into account the light paths and the interaction of light with different surface facets. Light is partially reflected, expressed by the surface reflectivity ρ, and can bounce off the cavity multiple times before it reaches the receiver area, i.e., the opening port of the cavity. Multibounce light arrives temporally later at the receiver due to its longer traveling path. Both effects, multibounce and traveling path lengths, lead to a certain temporal stretching of the impulse excitation.
In our first simulation, we use a transient rendering software 33 that is capable of calculating the temporal response of a scene due to certain illumination and sensing conditions. The scene has to be described as a set of mesh structures that can have any geometrical shape. Each mesh is defined by its vertex points (defining the area), surface orientation, and optical properties. The rendering software calculates the photon propagation within the scene even with multiple reflections. At each surface, the reflection is calculated taking into account a selected surface model (e.g., Lambertian surface). Further, for each photon, the path length is calculated, and the transient count rate and the number of bounces are recorded.
In Fig. 2, we show some results for three different cavities: (a) sphere, (b) hemisphere, and (c) model of a real cave. 29 In each scenario, S indicates the spot of laser irradiation, and P is the sensing area where radiation is reflected toward the sensor. We assume illumination with a Gaussian laser pulse. In the diagrams, we show the temporal return of photons with different numbers of bounces n (colored lines) and the total response (dashed) of the structures. The dotted line shows the total response when the structure is illuminated by a longer laser pulse.
In our simulation, we observe photons with n-bounces, n ∈ ½2;7. A single bounce (n ¼ 1) return is not visible as we assume a bistatic configuration and the impact point is outside the sensor field of view. Further, we observe an oscillation at the onset of the returning signal due to the return of photons with a low number of bounces n (n ∈ ½2;3; 4). At larger times, the signatures start to mix and form an exponential decay as expected in the empirical model; see Eq. (1).
For a spherical cavity [ Fig. 2(a)], we distinguish between photons with n ¼ 2, 3, and 4 bounces. For the hemisphere [ Fig. 2(b)], we observe a similar behavior except that we do not observe a signature for n ¼ 2 due to the ideal flat surface neglecting quasi-evanescent radiation transfer. In contrast, in later experiments, we do observe this type of signature in most cases. Finally, for the cave model, we observe a behavior close to the hemisphere model except the fact that we observe an additional 2-bounce signature due to the rough surface structure. Fig. 2 Numerical simulation of temporal response of (a) a sphere, (b) a hemisphere, and (c) a cave model. 29

Simulation of the Radiant Transfer by Random Sampling
In the second simulation, we adopted a Monte Carlo approach 32 to evaluate the radiant transfer between two surface elements. In our simulation, we consider only spherical and hemispherical cavities due to their known and simple geometries to estimate the signal peak position for different multibounce paths. We assume a cavitiy in the shape of a unit sphere with radius r ¼ 1 and the center in the origin of the coordinate system. Then, we randomly select two points on the cavity surface using a random sampling process (s← S ½−1;1) to determine for each sampling point a set of azimuth and elevation angles with ϕ ¼ 2πs and θ ¼ arccosðsÞ.
Further, in the hemisphere model, we distinguish between the hemisphere dome and the base. In both cases, the azimuth angle is determined as aforementioned, but we limit the sampling of second parameter to the upper half of the sphere. For the hemisphere base, we set the heigth z ¼ 0 and randomly sample the radius weighted by r ¼ ffiffi ffi s p . We randomly sample two points on the cavity surface and evaluate the path length (i.e., analogous to the response time) and the probability of the exchange of optical energy due to the radiant transfer. 34 Here, we evaluate the orientation of the emitting and receiving surface elements to calculate their effective size and estimate the amount of radiation transferred. Our findings are summarized in Fig. 3.
As depicted in Figs. 3(a) and 3(b), light enters the structure through an opening and is reflected multiple times before it exits toward the sensor. For each n-bounce path, the returning peak position is denoted by t n [with n indicating the number of internal reflections (bounces)]. Eventually, we observe light that is reflected at the opening rim to return a signal at t 0 without reflections inside the cavity.
Before the photon reaches sensing area P, we identify some significant photon paths. In the spherical model, the radiant energy is transferred between two points on the spherical hull (sphere-sphere), whereas in the hemispherical model, the energy is transferred between two points on the hemispheric dome (hemisphere-hemisphere) or between the dome and the base plate (hemisphere-base). Initially, we neglect radiant exchange within the base plate (base-base) because no energy can be transferred within an ideal plane ðlim θ i →π∕2 cosðθ i Þ → 0Þ.
Our random sampling approach evaluates the flight path (i.e., the L 2 -norm distance) and radiation transfer 34 (i.e., the amount of transferred energy by randomly selecting a pair of two points or surface elements, δA i ). Figure 3(d) shows the probability of a photon traveling a certain distance within the cavity hull. For each range, the radiant transfer weighs the probability of the exchange process, as illustrated in Fig. 3(c). The radiant transfer takes into account the face orientation (θ i ) and, thus, the effective size of the surface elements.
For the sphere-sphere exchange (blue), it can be clearly stated that the most probable distance is the sphere diameter, D sph−sph ¼ ∅ sph . In the hemisphere case, we obtain a most probable distances of D hemi−hemi ≈ 1.7r for the hemisphere-hemisphere exchange (orange) and D hemi−base ≈ 1r for the exchange with the base plate (green), with r ≡ h hemi being the hemisphere height.

Path Length Model for Spherical and Hemispherical Cavities
Let us evaluate the path length in the cavity models. In the sphere model, as illustrated in Fig. 3(a), a small portion of light is scattered or diffracted while propagating toward the cavity. Therefore, we expect to observe a first signature at t 0 due to a reflection of light at opening O. Within the sphere, first, the light passes the cavity by distance D OS to illuminate point S. From this point, the light can take various paths to reach P, the sensing area, and leaves the cavity opening toward the detector.
Due to the aforementioned analysis, most probable paths form the oscillations. A 2-bounce path is identified as the direct interaction of S and P, distance D SP to form a signature at t 2 . In a 3-bounce path, the light most likely interacts with face S 0 opposing S, D SS 0 ¼ ∅ sph , and reaches P via D S 0 P . Then, (4-bounce path) the light passes two times the diameter and so on. Each path with n bounces lead to a most probable time of arrival t n as summarized in Eqs.
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 3 8 9 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 3 0 7 Although we do not know the exact values of D OS , D SP , and D S 0 P , we determine the sphere diameter ∅ sph from the time of arrival difference between t 2 and t 4 ; see Eq. (6).
Further, we developed a similar model for hemisphere cavities, as illustrated in Fig. 3(b). Again, we name the distances from the opening to the illumination point as D OS and the distance to the sensing area as D SP . Similar to above, we distinguish the radiant transfer between the base and hemisphere dome (D hemi−base ≈ h hemi ) and interhemisphere dome (D hemi−hemi ≈ 1.7h hemi ).
The hemisphere model is summarized in Eqs. (7)- (10). In contrast to the sphere model, we find two 4-bounce paths, t 4 and t 4 0 , defining two different path lengths. We assume the second to be longer so that we distinguish both paths.
Similar to the sphere model, we determine the hemisphere height h hemi from the time of arrival difference between t 2 and t 4 [Eq. (11)].
Path length model for hemispherical cavities: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 1 4 3 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 9 9 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 7 3 5 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 6 8 2 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 6 ; 6 4 1 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 6 ; 6 1 1 Both models, sphere [Eq. (6)] and hemisphere [Eq. (11)], give the same value. Later in the analysis, we therefore discuss only the target size and use sphere diameter ∅ sph and hemisphere height h hemi equally, depending on the shape of the target. Further, an alternative model to determine the hemisphere dimensions is given in Eq. (12), using the multiphoton path that includes a radiant transfer within the hemisphere dome [see Eq. (10)].

Experimental Setup
Experimental investigations were carried out with a time-correlated single photon-counting setup consisting of a pulsed laser source (Picoquant LDH-P-C-640B, 80 ps, 640 nm, 2.5 MHz) and a photon-counting camera with a 32 × 32 single photon-counting avalanche diode (SPAD) array (Photon Force PF32), as shown in Fig. 4. The camera has a sampling resolution of 56 ps per time bin with a 10 bit sampling depth (1024 time bins). Thus, the camera can sample a time window of 57.3 ns. Overall, the setup has a temporal resolution of about 150 ps (i.e., a path length of 4.5 cm) due to the instrument response function (IRF). Further, ambient light is filtered by a single bandpass filter (SBF) with a spectral width of about 20 nm.
The camera and laser were used in a bistatic configuration and were separated by a base line d LC ¼ 30 cm. The cavity investigated was placed at a distance of D range ¼ 200 cm. The collimated laser illuminates a point within the cavity through an opening outside the sensor's field of view. The light propagates within the cavity, exits at the opening, and is eventually projected by a lens onto the SPAD sensor array. For every photon event, the camera records the event time and position on the array. The field of view of the camera was chosen to match the opening port of the investigated cavities, and for the analysis, the signal was averaged over the whole image. Illumination and sensing paths cross and use the same opening port.
Each sensor element records the time for every detection event. In our experiences, we use 500.000 measurement cycles to enable statistical analysis. The datasets are converted to Fig. 4 Our setup consists of a SPAD camera (C) and a collimated pulsed laser source (L). In a bistatic configuration (baseline d LC ), the laser illuminates a point inside the cavity at a distance D range . The camera captures reflected light and filters ambient light by a SBF.
histograms to obtain a transient signal, i.e., the count rate for every time bin within the sampling window. Thus, for every measurement, we obtain a multidimensional dataset giving the count rate for each sensor element (position u, v) and time bin (t j ). The data structure is illustrated in Video 1 in the Supplementary Material; see Fig. 5. The still image [ Fig. 5(a)] shows a frame of the count rate at a certain instant (u; v-cross-section at t j ) and [ Fig. 5(b)] the mean transient signal indicating main signal features. Within the first frames (t j < 7 ns), we observe the laser pulse propagating toward the cavity (light in flight) and the reflection at the cavity rim (t 0 ). Later, the signatures of multibounce photons (t 2 , t 3 , and t 4 ) are visible before the exponential decay of the signal amplitude (t j > 16 ns).

Experimental Investigation of Cavities
We studied a set of eight cavities, as depicted in Fig. 6, with different shapes (sphere, hemisphere, and polygonal), sizes (∅ ∈ ½15 cm; 3 m) and surface materials [spectralon, expanded polystyrene (EPS), white paper, 3D printed nylon, and fabric]. The three sphere targets have different sizes (∅ sph ∈ ½15;46; 50 cm) and materials (spectralon and EPS), and the EPS hemisphere has a height of h hemi ¼ 23 cm. The polygonal cavities are a cube and a dodecahedron, both made of white paper.
Finally, we investigate two use cases. First, a 3D printed cave model made of selective laser sintered (SLS) nylon. This cavity has an irregular surface (nonconvex roughness of about 1 cm) that can be roughly described as a hemisphere. The base resembles an oval with a depth and width of 30 cm × 40 cm, and the height of the hemisphere is about 15 AE 0.5 cm. The printed material has a smooth, diffusely reflective surface with high reflectivity (ρ ≈ 0.9). In contrast, the surfaces of natural caves have a more complex texture with a lower reflectivity. Both surface properties could change the signal amplitude but should not affect the round-trip time used in our analysis.
As a second use case, we consider a room size cavity with a volume of 3 m × 3 m × 2 m to demonstrate the up-scaling capability of our method. The room target was built of large white paper and aluminum compound plates coated with a diffuse reflective paint. In addition, we covered the floor with a white curtain fabric. All cavity parameters are summarized in Table 1. Figure 7 shows the transient signatures recorded at the different targets. In all cases, we observe oscillations due to multibounce photon paths and at least determine the position of the first peaks (t 0 , t 2 , t 3 , and t 4 ) using a multipeak detection algorithm. These values were used in our path length model for spherical and hemispherical cavities [Eqs. (6) and (11)] to estimate the cavity dimensions, sphere diameter ∅ sph , or hemisphere height h hemi , respectively.

Discussion
For the discussion of experimental results, first, we define the relative deviation error σ of the results from the actual dimensions of the target; see Eq. (13). Good results are obtained when this deviation error is close to zero: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 1 1 6 ; 1 5 9 The derived peak positions, the results from the analysis, and the deviation between the results and real dimensions of the cavity are summarized in Table 2. For all targets, we obtain a good estimation of the target dimensions, which is shown in the deviation errors with values below 20%. Only for the hemisphere target, we observe an error of 37%. Further, some peaks   In contrast, at larger cavities, the IRF influence on the determination of the peak positions is reduced as the peaks are well separated. We find good values for the sphere targets independent from the size and material. For the polygonal targets (cube and dodecahedron), again, we observe a good agreement between the model and target dimensions. Only the dimensions of the dedecahedron are estimated too low, resulting in a deviation of 12%. Overall, we can see that the polygonal structure does not affect the analysis results.
In the use cases, cave and room, we estimate the depth nearly perfectly (σ ≈ 3%). This demonstrates the usefulness of the hemisphere and sphere models [Eqs. (6) and (11)] for estimating the size of a cavity even in the presence of structured or textured surfaces.
Although the cave has a nearly hemispheric shape, the analysis of the ideal hemisphere target shows a significant deviation from the actual dimensions. This observation feeds the assumption that, in the case of the ideal hemisphere, further processes may play a role. For instance, due to the large transmission angles [θ i , see Fig. 3(c)], radiant transfer could be reduced within the hemisphere base. Then, we may not observe an oscillation at t 2 or t 4 . In that case, these peak positions (t 2 and t 3 in Table 2) would become t 3 and t 4 0 , and the hemisphere model, Eq. (12), gives a height h hemi 0 ¼ 24.7 cm with σ ¼ 7.4%. This value is much closer to the real hemisphere shape.

Conclusion
In this paper, we have established a method to characterize the dimensions of a cavity by measuring its temporal response to a short optical pulse. Using time correlated single photon counting imaging, we observe an oscillation of the reflected light due to different multibounce paths. Further, we have established two empirical models, the sphere and the hemisphere model, to determine the cavity's diameter or height, respectively.
In our experimental investigation, we have studied the temporal response of eight different cavities, which represent a variation of shape, size, surface structure, and material. Until now, we have only investigated targets with high diffuse reflective surfaces.
We showed that we can reliably estimate the dimensions of the sphere, polygon, and other cavities. Our method was not affected by the material (spectralon, paper, EPS, and fabrics) or target size (dimensions between 15 cm and 3 m). Although we established a hemisphere model, the results of the hemisphere target had the largest error of 37%. This significant deviation may Table 2 Position of peaks in the signature of different targets and estimation of the target dimensions using the sphere and hemisphere models. be due to the hemisphere having an almost perfectly flat base in the experiments. Here, we showed that assuming a reduced transfer of radiant energy within the base plane and the use of alternative n-bounce paths led to a better estimation of the dimension. In contrast, the cave cavity showed good results with the sphere model even if it had a nearly hemisphere shape. Our experiments with the room target showed that the method was also applicable for larger cavities and led to a very accurate estimate of the cavity size. Unfortunately, the SPAD camera used was limited to an observation window of 57 ns, which was sufficient for a investigation of target determination of 3 m. To sense longer cavities, we have to use a photon-counting device with a larger sampling depth or lower sampling resolution, which would effect the overall depth resolution.
Martin Laurenzis received his MS degree in physics from the University of Dortmund, Germany, in 1999 and his PhD in Electrical Engineering from the University of Aachen, Germany, in 2005. Since 2004, he has been at the French-German Research Institute Saint-Louis as a senior researcher in the Advanced Visionics and Processing group. He is a member of SPIE.
Frank Christnacher received his MS degree in physics in 1988 and his PhD in optical data processing and pattern recognition in 1992 from the University of Haute-Alsace, France. He joined French-German Research Institute of Saint-Louis in 1988 and currently leads the Advanced Visionics and Processing group. He is a member of SPIE.