## Abstract

We analyze multi-bounce propagation of light in an unknown hidden volume and demonstrate that the reflected light contains sufficient information to recover the 3D structure of the hidden scene. We formulate the forward and inverse theory of secondary scattering using ideas from energy front propagation and tomography. We show that using Fresnel approximation greatly simplifies this problem and the inversion can be achieved via a backpropagation process. We study the invertibility, uniqueness and choices of space-time-angle dimensions using synthetic examples. We show that a 2D streak camera can be used to discover and reconstruct hidden geometry. Using a 1D high speed time of flight camera, we show that our method can be used recover 3D shapes of objects “around the corner”.

© 2012 Optical Society of America

## 1. Introduction

Estimating the 3D geometry of scenes is a central problem in imaging and traditional time-of flight techniques are one of the most popular methods to achieve this. The basic idea is that when a short duration pulse of light is transmitted in a particular direction, the pulse reflects isotropically from the first obstacle in that direction and returns back to the source location. A collocated receiver computes the time of flight from which the distance to the obstacle can be readily calculated. While this simple technique is extremely effective for unoccluded scenes, the task of estimating 3D geometry of hidden objects (objects with no line of sight from both source and sensor) is significantly more involved. Such a technique, if developed could have several potential applications such as surveillance, search and rescue, car navigation and even medical applications such as endoscopy. Inspite of the potential for widespread impact, there are few techniques that directly address this problem. In this paper, we extend time of flight techniques to enable accurate and robust estimation of 3D geometry of hidden objects. This is achieved by explicitly modelling the secondary and tertiary bounces from the unknown hidden objects and developing a mathematical framework based on sparse approximations and back projection to invert the effect of these secondary bounces. The result is an ability to extract 3D information for hidden objects which has hitherto been impractical.

#### Related work

Recent work in Light Transport (LT) in optics, computer graphics and computer vision has shown ability to recover surprising geometric and photometric information about the scene from diffuse scattering. The pioneering work in dual photography [1] shows one can exploit the second bounce to recover images of hidden objects. The theory of inverse light transport [2] can be used to eliminate inter-reflections from real scenes. The frequency domain properties of direct and global components of scattered light can be exploited to recover images of objects behind a shower curtain [3]. Three bounce analysis of a time-of-flight camera can recover hidden 1-0-1 barcodes [4] and single object position [5] while form factors can be recovered from just two bounces [6]. In [7], Raskar and Davis proposed inverse analysis using a 5D time-light transport matrix to recover geometric and photometric scene parameters. The recent paper [8] presents an approach to recover 3D shape from tertiary diffuse scattering.

Similar to these and other inverse light transport approaches [2], we illuminate one spot at a time with our pulsed laser projector and record the reflected light after its interaction with the scene. We record an extra time dimension of light transport with an imaging system that uses a short duration pulse and a time-of-flight camera. A schematic of our setup can be found in Fig. 1. We show that this extra temporal dimension of the observations makes the 3D structure of a hidden scene observable. The lack of correspondence between scene points and their contributions to the captured streak images after multiple bounces is the main challenge and we show computational methods to invert the process.

#### Contributions

The recent publication [8] introduces the concept of time of flight imaging and provides a proof of principle experiment using a highly robust but imprecise modified filtered backprojection algorithm. Here we go beyond this innitial experiment by evaluating several different reconstruction techniques as well as different experimental setups. In particular we

- apply a SPGL1 compressive reconstruction to new experimental data.
- apply a linear inversion algorithm to real data and show that linear inversion is possible if the calibration bias of the streak camera is taken into account.
- provide results of a superior reconstruction with Cosamp from simulated data.
- provide a mathematical model that describes our reconstruction problem as a computed tomography.
- provide additional data and results for the filtered backprojection algorithm.

#### Limitations and scope

Our multi-scattering based method is inherently limited by the signal to noise ratio due to serious loss in scattering. We require fine time-resolution and good sensitivity. Scenes with sufficient complexity in geometry (volumetric scattering, large distances) or reflectance (dark surfaces, mirrors) can pose problems. We limited our scope to scenes with approximately Lambertian reflectance (no directionally varying BRDF) and with few or no partial occlusions for simplicity. In order to achieve good spatial resolution, we require the visible area to be sufficiently large and the object to be sufficiently close to this visible area, so that we have a large ’baseline’ to indirectly view hidden objects.

## 2. Modeling propagation of a light pulse for multiple bounces

To develop a method for reconstruction of images from time-of-flight information we need to obtain an adequate, physically realistic model of light propagation in the scene. The model can then be used to understand the forward process of rendering time-of-flight streak camera images from a known geometry. Model can also be used for the inverse process of reconstructing the hidden geometry from the images. The light transport model is sketched in Fig. 1. Light emitted by a laser (B) in a collimated beam strikes a diffuse surface, the *diffuser wall* at a single point L. The laser emits pulses that are much shorter than 1 ps in time, or 0.3 mm in space. From point L light is diffusely scattered in all directions. A small fraction of the light travels a distance *r _{l}* and strikes point s. At s the light is diffusely reflected once again and some of it travels to point w covering the distance

*r*. From w some light travels back to the camera C. The camera observes points on the wall with high time resolution, such that light having traveled different total distances through the scene is detected at different times.

_{c}By scanning the position of L, and exploiting the spatial resolution of the camera we can probe different sets of light paths. 5D light transport captured in this way contains information about the hidden object and allows us to uniquely identify the geometry of hidden objects. Consider for example just one reflecting spot s in the hidden scene as indicated in Fig. 1. The space time image of this spot for a given laser position L is a hyperbola. The curvature and position of the hyperbola indicate the position of s.

In our experimental setup we make use of a streak camera with a 1 dimensional field of view. We thus only capture 4D light transport. This does, however, still allow us to reconstruct the hidden geometry. The reconstruction quality along the axis perpendicular to the camera field of view is, however, compromised.

#### 2.1. Space-time warping for bounce reduction

The light path from laser to camera can be divided into four straight line segments with three bounces in between. The first segment is the collimated laser beam travelling from the laser source to the wall. In the second segment the laser spot on the wall behaves as a point light source and light travels from there into the hidden scene. The third segment(s) involve scattering from the hidden object. For the fourth segment light travels from the wall to the camera *which is focused on the wall*. The data received by the camera *I _{C}*(

*p*,

*t*) has two degrees of freedom - space and time. Since the first and fourth segments are focused, we can apply transforms to the streak images to account for their effects (Please see Fig. 2). The effect of the first segment can be accounted for by shifting the streak images in time for each laser location on the wall. To account for the effect of the fourth segment we use the distances between the wall pixels and the camera sensors (||

*C*−

*w*||). We assume that we know the geometry of

*R*to calculate the camera-to-wall homography and the correspondence between the camera and wall pixels. The following mathematical formulation is a concise representation of this concept. Let

*H*be the projective transformation (homography) mapping coordinates on

*R*to camera coordinates. The time shift ||

*C*−

*w*|| by the distance from camera to screen varies hyperbolically with the pixel coordinate

*w*. Here we have set the speed of light to

*c*= 1 for simplicity, or in other words we measure time in distance units. Note that we don’t require to adjust for a cos(

*θ*) factor or 1/

*r*

^{2}fall off because the camera sensor integrates for more wall pixels if they are farther away. The following equation describes a coordinate transform from camera to wall where

*I*(

_{R}*w*,

*t*) is the intensity of light received by world pixels as a function of world spatial and temporal coordinates,

*I*is the streak image,

_{C}*H*is the homography.

#### 2.2. Scattering of a pulse

### 2.2.1. Generating streak photos

Consider the scenario in Fig. 1. Let us analyze the scattering of light in the second and third segments. For simplicity we model the hidden object as a collection of unfocused emitters sending impulses *I _{s}*(

*s*,

*τ*) at times

*τ*. We model the receiver as an array of unfocused receivers which capture photons at picosecond resolution. We can achieve this configuration experimentally by using a picoseond resolution streak camera focused on the diffuser wall pixels. The following equation represents the data recorded by the streak camera after making the mentioned assumptions. We also ignore the changes in normals for sender surface and receiver surface.

*w*∈

*R*,

*s*∈

*S*,

*t*,

*τ*∈ ℝ and

*r*= ||

_{c}*w*−

*s*||. and

*I*(

_{R}*w*,

*t*) is the intensity observed at

*w*∈

*R*at time

*t*. After removing the time shifts as described in the previous section and applying transforms from the calculated homography we can further simplify the equation and remove the

*δ*required to adjust for receiver camera distances. The following equation provides a mathematical summary of this analysis. Note that over all these equations we assume that the receiver and sender are perfectly Lambertian and ignore the variation in normal vectors. Equation (2) hence becomes

### 2.2.2. Hyperbolic contribution

Let us analyze the relationship between the time when a sender emits a pulse and the time and location of a receiver detecting the light. For a fixed sender the response function is a hyperboloid in space and time given by the following mathematical equation. The parameters of the hyperboloid depend on location of sender, a lateral displacement leads to shifts, while a displacement in depth corresponds to flattening. Change in sender time equates to a constant time shift for responses to any of the receivers:

where*u*,

*v*are the two coordinates of

*w*in the receiver plane. This equation describes an ellipsoid in sender location if we fix the laser and receiver location. The ellipsoid’s parameters depend on the time at which a receiver receives an impulse. The laser spot and the receiver (on the wall) constitute the two foci of this ellipsoid. The eccentricity depends on the time when the impulse is received.

## 3. Forward model: elliptical tomographic projection

In this section we rephrase the above approximation to the forward light transport using notions from tomography. In an idealized case, the inverse problem of recovering the hidden shape can be solved explicitly. We use this explicit solution to inspire our algorithm for real world scenarios.

#### 3.1. Elliptical tomography problem description

Our problem has similarities to tomographic projection. Let us rewrite Eq. (3) in the following form

*W*(

*x*) =

*I*

*δ*(

_{S}*x*) is a delta function with support on the surface

*S*to be reconstructed. Apart from the 1/

*r*

^{2}factors, which typically vary slowly over the region of interest, we hence see that individual streak image pixels measure elliptical projections of the world volume

*W*(

*x*). Due to its similarity with traditional tomography, we call this problem

*elliptical tomography*. Note however that there are also key differences to traditional tomography: (i) The recorded projections live in a higher dimensional (5D) space than the world (3D). (ii) The projections are along (2D) ellipsoids instead of (1D) lines. This makes the analysis more complicated.

#### 3.2. Challenges and missing cones

It is instructive to consider the above tomography problem in the limit when the object is small compared to the distance to the diffuser wall. In this case the elliptic tomography problem reduces to a planar tomography problem, see Fig. 3. Each pair of a camera point and a laser position on the diffuser wall (approximately) measures intersections of the target object with planes whose normals agree with the normals to the ellipsoids. By the standard Fourier slice theorem, each line of each streak image will hence measure one line of the Fourier transform of the object, in the direction of the normal. Unfortunately, in our situation these normals cover *only* a limited region of the unit sphere. Hence without additional priors it is not possible to reconstruct the Fourier transform of the target object in the missing directions. This is the missing cone problem well known from traditional tomography. Experimentally we get a very good resolution in the depth (orthogonal to the wall) direction, while transverse (parallel to the wall) high frequency spatial features tend to get lost.

## 4. Inverse algorithm: filtered back projection

In this section we give a description of our reconstruction algorithm.

#### 4.1. Overview of the algorithm

The imaging and reconstruction process consists of 3 phases:

**Phase 1: data acquisition.**We direct the laser to 60 different positions on the diffuser wall and capture the corresponding streak images. For each of the 60 positions 100 streak images are taken and overlayed to reduce noise.**Phase 2: data preprocessing.**The streak images are loaded, intensity corrected and shifted to adjust for spatiotemporal jitter.**Phase 3: 3D reconstruction.**The clean streak images are used to reconstruct the unknown shape using our backprojection-type algorithm.

The first of the three phases has been described above. Let us focus on Phases 2 and 3.

### 4.1.1. Phase 2: data preprocessing

**Timing correction.**To correct for drift in camera timing synchronization (jitter) both in space and time we direct part of the laser directly to the diffuser wall. This produces a sharp “calibration spot” in each streak image. The calibration spot is detected in each image, and the image is subsequently shifted in order to align the reference spot at the same pixel location in all streak images. The severity of the jitter is monitored in order to detect outliers or broken datasets.**Intensity correction.**To remove a common bias in the streak images we subtract a reference (background) image taken without the target object being present in the setup.**Gain correction.**We correct for non-uniform gain of the streak camera’s CCD sensor by dividing by a white light image taken beforehand.

### 4.1.2. Phase 3: 3D reconstruction

**Voxel grid setup.**We estimate an oriented bounding box for the working volume to set up a voxel grid (see below).**Downsampling (optional).**In order to improve speed the data may be downsampled by discarding a fraction of the cameras pixels for each streak image and/or entire streak images. Experiments showed that every second camera pixel may be discarded without losing much reconstruction accuracy. When discarding entire streak images, is is important that the laser positions on the diffuser wall corresponding to the remaining images still cover a large area.**Backprojection.**For each voxel in the working volume and for each streak image, we compute the streak image pixels that the voxel under consideration might have contributed. Concretely, the voxel at location*v*can contributed to a pixel corresponding to a point*w*on the wall at time*t*if Here*C*is the camera’s center of projection and*L*is the laser position as above. Let us call the streak image pixels satisfying this condition the*contributing pixels*. We compute a function on voxel space, the*heatmap H*. For the voxel*v*under consideration we assign the value Here the sum is over all contributing pixels*p*, and*I*is the intensity measured at that pixel. The prefactor corrects for the distance attenuation, with_{p}*α*being some constant. We use*α*= 1.**Filtering.**The heatmap*H*now is a function on our 3 dimensional voxel grid. We assume that the axis of that grid are ordered such that the third axis faces away from the diffuser wall. We compute the*filtered heatmap H*as the second derivative of the heatmap along that third axis of the voxel grid. This can be done by shifting and subtracting world volume along correct axes._{f}- The filtered heatmap measures the confidence we have that at a specific voxel location there is a surface patch of the hidden object.
**Thresholding.**We compute a sliding window maximum*M*of the filtered heatmap_{loc}*H*. Typically we use a 20x20x20 voxel window with resolution around 2mm. Our estimate of the 3D shape to be reconstructed consists of those voxels that satisfy the condition where_{f}*M*= max(_{glob}*H*) is the global maximum of the filtered heatmap and_{f}*λ*,_{loc}*λ*are constants. Typically we use_{glob}*λ*= 0.45,_{loc}*λ*= 0.15_{glob}**Compressive reconstruction**- We use techniques like SPGL1, and CoSAMP [9] as an alternative to back projection and filtering. We rely on the fact that the 3D voxel grid is sparsely filled, containing surfaces which can occupy only one voxel in depth. Since SPGL1 uses only matrix vector multiplications.
**Rendering.**The result of thresholding step is a 3D point cloud. We use the Chimera rendering software to visualize this point cloud.

In order to set up the voxel grid in the first step we run a low resolution version of the above algorithm with a conservative initial estimate of the region of interest. In particular, we greatly down sample the input data. By this we obtain a low resolution point cloud which is a coarse approximation to the object to be reconstructed. We compute the center of mass and principal axis of this point cloud. The voxel grid is then fit so as to align with these axis.

#### 4.2. A remark about the filtering step

To motivate our choice of filter by taking the second derivative, let us consider again the planar (“far field”) approximation to the elliptical tomography problem as discussed in section 3.2. In this setting, at least for full and uniform coverage of the sphere with plane directions, the theoretically correct filtering in the filtered backprojection algorithm is a |*k*|^{2} filter in Fourier space. In image space, this amounts to taking a second derivative along each scanline. This motivates our choice of filter above. We tested both filtering by the second derivative in image and world space and found that taking the derivative approximately along the target surfaces normal yielded the best results.

#### 4.3. Application of CoSAMP

We attempted compressive reconstructions with CoSAMP on real as well as simulated data. We found that CoSAMP performs better than backprojection on simulated data. The situation is reversed with data from the actual system. We attribute the poor performance of CoSAMP and other linear equation based methods with real data to bias in the data due to imperfect calibration in the system. A result of a simulated CoSAMP reconstruction is shown in Fig. 4. This is something that we will like to explore in the future.

## 5. Experiments

Our streak camera is a Hamamatsu C5680, with an internal time resolution of 2 picoseconds per sample and a temporal point spread function of about 15 ps. The impulse response of the camera therefore fits well to a gaussian with 15 ps full width at half maximum. The temporal point spread function can be changed at the expense collecting less light. We use a mode-locked Ti:Sapphire laser to generate pulses at 795 nm wavelength with about 50 femtosecond pulse duration at a repetition rate of 75 MHz. The laser’s average output power is about 500 mW. The streak camera has a one dimensional field of view, imaging a line in the scene. It provides a two dimensional image in which one dimension corresponds to space and the other to time [10].

We use a Faro Gauge measurement arm to calibrate the laser and camera. We treat the camera and laser as a rigid pair with known intrinsic and extrinsic parameter [11]. The visible parts of the geometry could be measured with the laser directly using time of flight with micrometer precision. Methods like LiDAR and OCT can achieve this and are well understood. In the interest of focusing on the novel components of our system we instead measure the visible parts of our scene with the Faro Gauge. We also collect ground truth data to validate our reconstructions.

#### 5.1. Results

We recorded series of streak images for several simple 3D scenes, comprised of white Lambertian objects. We used 30–60 laser positions, spread over a 20 × 40 cm wall. The reconstructed surfaces are displayed in Figs. 5, 6, 7, 8 from setup Fig. 9. To produce the 3D pictures from the volumetric data, we use the Chimera visualization software [12] for thresholding and rendering.

#### 5.2. Performance evaluation

We conducted several experiments to test the performance of our experimental setup. This includes verifying the spatial and temporal resolution of the camera and the resolution obtained in a simple hidden scene. The resolution in the hidden scene depends greatly on the position in the scene and overall scene complexity. Experiments indicate that for a simple patch a precision of about 500 *μ*m perpendicular to the wall and 1 cm parallel to the wall is achievable.

## 6. Future directions

We have shown that the goal of recovering hidden shapes is only as challenging as the current hardware. The computational approaches show great promise. But, on the hardware front, emerging integrated solid state lasers, new sensors and non-linear optics will provide practical and portable imaging devices. Our formulation is also valid for shorter wavelengths (e.g., x-rays) or for ultrasound and sonar frequencies in large scenes where diffraction can be neglected. Beyond geometry, one maybe able to recover full light transport and bidirectional reflectance distribution function (BRDF) from a single viewpoint to eliminate encircling instrumentation. Our current method assumes friendly reflectances, i.e., a non-zero diffuse component towards the wall. Non-lambertian reflectance and partial occlusions will create non-uniform angular radiance. While our method is robust towards deviations from a Lambertian reflector, reconstructions will benefit from adjusting the model to account for the particular scenario. Generally, non-Lambertian reflectors provide potential challenges by introducing the reflectance distribution as another unknown into the problem and by increasing the dynamic range of the intensities detected by the system. On the other hand, they may provide the benefit of higher reflected intensities and better defined reflection angles that could be exploited by an adequate reconstruction method. The visible wall need not be planar and one can update the *r _{l}* and

*r*distances from a known model of the visible parts. Supporting refraction involves multiple or continuous change in the path vector. Atcheson et. al. have shown a refraction tomography approach [13] which could be extended.

_{c}Initial applications may be in controlled settings like endoscopy, scientific imaging and industrial vision. This will require addressing more complex transport for volumetric scattering (e.g. for tissue) or refracting elements (e.g. for fluids). A very promising theoretical direction is in inference and inversion techniques that exploit scene priors, sparsity, rank, meaningful transforms and achieve bounded approximations. Adaptive sampling can decide the next-best laser direction based on current estimate of the carved hull. Future analysis will include coded sampling using compressive techniques and noise models for SNR and effective bandwidth.

We used the COSAMP matching pursuit algorithm which allows us to explore the sparsity of the solution. We tested both the backprojection and linear equation based methods on both artificially generated and real data. It turned out that for artificial data the COSAMP based reconstruction algorithm was generally superior the backprojection algorithm. The backprojection algorithm can recover objects front-to parallel to the wall quite well, but fails for highly sloped surfaces. However, for real data the linear equation based methods were very sensitive to calibration errors, i.e., errors in the matrix *A* above. One can obtain results for very good datasets, after changing *A* slightly to account for intricacies of our imaging system like vignetting and gain correction. On the other hand the backprojection algorithm turned out to be quite robust to calibration errors, though they can deteriorate the obtained resolution. Since we typically have to deal with some calibration error (see section 5), we used the backprojection algorithm to obtain the real data reconstruction results of this paper.

Designing a perfectly calibrated system is an interesting research challenge. Research directions here include the development of a model for lens distortions in the streak camera system and a scheme for their compensation, as well as a method of actively and autonomously recalibrating the system to account for day to day variations and drift in laser operation and laser-camera synchronization wile capturing data.

## 7. Conclusion

The ability to infer shapes of objects beyond the line of sight is an ambitious goal but it may transform recoding of visual information and will require a new set of algorithms for scene understanding, rendering and visualization. We have presented a new shape-from-x approach that goes beyond the abilities of today’s model acquisition and scanning methods. Light transport with a time component presents a unique challenge due to the lack of correspondence but also provides a new opportunity. The emphasis in this paper is to present a mathematically sound forward model and several novel inversion processes. The nonlinear component due to jitter and system point spread function makes the ultra-fast imaging equipment difficult to use. So the physical results can only be treated as a proof-of-concept.

Currently ultrafast imaging systems are mostly utilised in (bio-)chemical or physical applications. We hope that our work will spur more research into the possible applications of such imaging systems in computational photography and computer vision. To assist the community, we will make image datasets and Matlab code freely available online.

## Acknowledgments

This work was funded by the Media Lab Consortium Members, DARPA through the DARPA YFA grant, and the Institute for Soldier Nanotechnologies and U.S. Army Research Office under contract W911NF-07-D-0004. Ashok Veeraraghavan was supported by NSF Grants IIS-1116718 and CCF-1117939.

## References and links

**1. **P. Sen, B. Chen, G. Garg, S. R. Marschner, M. Horowitz, M. Levoy, and H. P. A. Lensch, “Dual photography,” ACM Trans. Graphics **24**(2), 745–755 (2005). [CrossRef]

**2. **S. M. Seitz, Y. Matsushita, and K. N. Kutulakos, “A theory of inverse light transport,” in *Proc. of IEEE ICCV* (2005), Vol. 2, pp. 1440–1447.

**3. **S. K. Nayar, G. Krishnan, M. D. Grossberg, and R. Raskar, “Fast separation of direct and global components of a scene using high frequency illumination,” ACM Trans. Graphics **25**, 935–944 (2006). [CrossRef]

**4. **A. Kirmani, T. Hutchison, J. Davis, and R. Raskar, “Looking around the corner using transient imaging,” in *Proc. of IEEE ICCV* (2009), pp. 159–166.

**5. **R. Pandharkar, A. Velten, A. Bardagjy, M. G. Bawendi, and R. Raskar, “Estimating motion and size of moving non-line-of-sight objects in cluttered environments,” in Proc. of CVPR (2011), pp. 265–272. [CrossRef]

**6. **S. Liu, T. Ng, and Y. Matsushita, “Shape from second-bounce of light transport,” in Proc. of ECCV (2010), pp. 280–293.

**7. **R. Raskar and J. Davis, “5D time-light transport matrix: What can we reason about scene properties,” MIT Technical Report (2008), http://hdl.handle.net/1721.1/67888.

**8. **A. Velten, T. Willwacher, O. Gupta, A. Veeraraghavan, M. G. Bawendi, and R. Raskar, “Recovering three-dimensional shape around a corner using ultra-fast time-of-flight imaging,” Nat. Commun. **3**, 745 (2011). [CrossRef]

**9. **D. Needell and J. Tropp, “CoSaMP: Iterative signal recovery from incomplete and inaccurate samples,” Appl. Comput. Harmon. Anal. **26**, 301–321 (2009). [CrossRef]

**10. **Hamamatsu, “Hamamatsu streak camera tutorial,” http://learn.hamamatsu.com/tutorials/java/streakcamera/.

**11. **D. Forsyth and J. Ponce, *Computer Vision, a Modern Approach* (Prentice Hall, 2002).

**12. **E. Pettersen, T. Goddard, C. Huang, G. Couch, D. Greenblatt, E. Meng, and T. Ferrin, “UCSF Chimera–a visualization system for exploratory research and analysis,” J. Comput. Chem. **25**, 1605–1612 (2004). [CrossRef] [PubMed]

**13. **B. Atcheson, I. Ihrke, W. Heidrich, A. Tevs, D. Bradley, M. Magnor, and H. Seidel, “Time-resolved 3d capture of non-stationary gas flows,” ACM Trans. Graphics **27**, 1–9 (2008). [CrossRef]