1 Introduction

When an image is recorded by a general camera, various optical phenomena are collectively recorded. Light rays leaving from the light source change in intensity by repeating reflection and scattering, and the final intensity is recorded by the camera.

In recent years, temporal point spread function (PSF) of a scene, which is a response to the impulsive light, is attracting attention because it can be used for analysis of interreflection and scattering. By analyzing temporal PSFs, how the light transports in the scene can be analyzed in detail.

Temporal PSF of the scene can be obtained using an interferometer [1], holography [2], and femtosecond-pulsed laser [3], which requires complex optics and expensive devices. Heide et al. [4] firstly recover a temporal PSF in an inexpensive way using an amplitude-modulated continuous wave ToF camera. Using a continuous wave ToF camera, the propagation of light in units of nanoseconds can be visualized [5–7].

In this paper, we propose a method for estimating the temporal PSF using a pulsed ToF camera. Using pulsed ToF cameras, low-resolution temporal PSFs can be straightforwardly obtained, e.g., tens of nanoseconds, which is insufficient for analyzing scattering that is finished within 1 ns. To achieve sub-nanosecond resolution, we introduce an additional circuit for delaying the light emission. By delaying the light emission in sub-nanoseconds, we recover temporal PSFs via computation and achieve sub-nanosecond resolution.

Our contributions are twofold. Firstly, the impulse response of the scene can be estimated even with the pulsed ToF camera. Secondly, we achieve sub-nanosecond resolution. We show a simple PSF recovery by delaying the light emission and computation.

2 PSF estimation using ToF camera

In this section, we explain the method that estimates temporal PSF using pulsed ToF camera. Before explaining the method, we start from the basic theory of the pulsed ToF camera.

2.1 Principle of ToF camera operation

A ToF camera is originally a device for distance measurement by measuring the time between the light emission and reflecting back to the camera. The distance to the object d is obtained as

$$\begin{array}{*{20}l} d = \frac{c \tau}{2}, \end{array} $$
(1)

where τ is the measured round trip time of light and c is the speed of light. Because the observation of the ToF camera is independent for each pixel, the camera pixel p is omitted for simplicity. In this paper, we assume a pulsed ToF camera, which emits a square pulsed light for a few tens of nanoseconds as shown in Fig. 1. Synchronizing with the light emission, two images I 1 and I 2 exposed at the same time as the light emission width are obtained. From these two images, the round trip time τ can be obtained as

$$\begin{array}{*{20}l} \tau = \frac{I_{2}}{I_{1}+I_{2}}T, \end{array} $$
(2)
Fig. 1
figure 1

Principle of ToF camera operation

Full size image

where T is the width of light emission and exposure. Moreover, because the observable range of the depth is limited to the width of the square wave, the start of the exposure is shifted to change the region of interest. In such a case, Eq. (2) can be written as

$$\begin{array}{*{20}l} \tau = \frac{I_{2}}{I_{1}+I_{2}}T + t_{\text{offset}}, \end{array} $$
(3)

where t offset is the shift of the exposure.

2.2 Distortion of reflected light due to PSF

The reflected light is assumed to be an ideal square wave in principle. However, in practice, the waveform of the reflection light is distorted according to the geometrical and optical characteristics of the scene as shown in Fig. 2. For example, if the scene includes subsurface scattering or inter-reflection, the temporal PSF r(t) of the scene is spread along with the time domain; therefore, the waveform of the reflected light is no longer the square wave.

Fig. 2
figure 2

Waveform of the reflection is distorted due to subsurface scattering

Full size image

The reflection waveform reaching the ToF camera can be represented by the convolution of the emission wave l(t) and the temporal PSF r(t). Observation i of the ToF camera is then represented as

$$\begin{array}{*{20}l} i &=\int g(t) \left(l \otimes r(t) \right) dt, \end{array} $$
(4)

where ⊗ is the convolution operator and g(t) is the exposure function of the ToF camera representing the exposure sensitivity at a certain time t. We assume that g(t) is a binary function, whether the exposure is performed or not, but in general, it can take a continuous value.

2.3 PSF estimation with delayed light emission

Temporal PSF varies depending on the surface shape and translucency of the object; it has important information about the material and structure of the object. In order to analyze the properties of the object in detail, we are interested in recovering the temporal PSF with high resolution.

If the emission and exposure of the ToF camera can be freely controllable, the temporal PSF of the scene could be easily obtained. Shortening the light emission and exposure widths to the infinitesimal time, the impulsive light emission l(t) and the exposure g(t) can be regarded as the delta function δ(t). Appropriately shifting the start of the exposure, the temporal PSF can be directly obtained as

$$\begin{array}{*{20}l} i(t) &= \int \delta(s - t) ((\delta \otimes r)(s)) ds\\ &= r(t), \end{array} $$
(5)

where s is the variable of integration on the time domain and t corresponds to the shift of exposure. However, it is difficult to shorten the emission and exposure width to the infinitesimal time because of the limitation of the logical gates. Using the typical ToF cameras, only tens of nanoseconds resolution can be achieved.

Even if the light emission and exposure width can not be shortened, the PSF can be recovered from observations with shifting the exposure. Shifting the exposure finely, a blurred PSF \(\bar {r}(t)\) can be observed as

$$\begin{array}{*{20}l} \bar{r}(t) &= \int g(s - t) (l \otimes r (s)) ds \\ &= g^{*} \otimes l \otimes r(t), \end{array} $$
(6)

where g ∗ is a exposure function whose time axis is inverted. Because \(\bar {r}(t)\) is convoluted by illumination and exposure, the original PSF r(t) can be estimated by deconvolution.

To realize the fine shift of the exposure, we adopt to delay the light emission by inserting an existing delay circuit as shown in Fig. 3. Delaying the light emission is equivalent to reversely shifting the exposure; hence, the shift of the exposure is easily controlled. The circuit delays the synchronization signal in sub-nanosecond; hence, we can recover high-resolution PSF by a deconvolution approach.

Fig. 3
figure 3

Light emission is delayed by an inserted circuit between the camera and the light source. Delayed light emission is equivalent to reversely shifting the exposure

Full size image

2.4 Numerical implementation

The observation i j at j-th delay setting is given by

$$\begin{array}{*{20}l} i_{j} &=\int g_{j}(t) ((l \otimes r)(t)) dt, \end{array} $$
(7)

where g j is the j-th exposure function corresponding to the j-th delay of light emission. Descretizing Eq. (7),

$$\begin{array}{*{20}l} i_{j} = \mathbf{g}_{j}^{T} \mathbf{L} \mathbf{r}, \end{array} $$
(8)

where g j is a vector representing j-th exposure, L is a convolution matrix of emitted wave, and r is a vector of discretized temporal PSF, whose resolution is the same as the width of the delay. Superposing all the observations i

$$\begin{array}{*{20}l} \mathbf{i} = \mathbf{G} \mathbf{L} \mathbf{r}, \end{array} $$
(9)

where

$$ \left\{\begin{array}{ll} \mathbf{i} &= \left[ i_{1} \;\; i_{2} \;\; \cdots\;\; i_{m}\right]^{T}\\ \mathbf{G} &= \left[\mathbf{g}_{1}^{T} \;\; \mathbf{g}_{2}^{T} \;\; \cdots \;\; \mathbf{g}_{m}^{T} \right]^{T}, \end{array}\right. $$
(10)

and m is the number of observations.

Since both the exposure setting and the emission waveform are known, the matrix G L is a known matrix. Therefore, if the number of the observation is sufficiently large, the temporal PSF r can be estimated by the least squares method as

$$\begin{array}{*{20}l} \mathbf{r} = (\mathbf{G}\mathbf{L})^{+} \mathbf{i}, \end{array} $$
(11)

where (G L)+ is the pseudo-inverse matrix.

In practice, matrix G L is ill-conditioned hence the calculation as shown in Eq. (11) is unstable. However, using the property that the intensity of light does not have a negative value, it is expected to estimate robustly against the instability of calculation. Estimation of temporal PSF can be estimated by the least square method with the non-negative constraint as

$$ \underset{\mathbf{r}}{\text{argmin}} \; \;\left\lVert\mathbf{i} - \mathbf{G}\mathbf{L}\mathbf{r}\right\rVert_{2}^{2}\;\; \textrm{s.t.} \;\;\mathbf{r} \succeq 0. $$
(12)

Since Eq. (12) is a convex optimization problem, the optimal solution can be obtained in a polynomial time.

3 Experiments

We confirm the effectiveness of our method via real-world experiments. We build a measurement system as shown in Fig. 4. A delay circuit as shown in (c), where DS1023 IC (Maxim Integrated) and Arduino micro-controller are equipped, is inserted between a ToF camera S11963-01CR (Hamamatsu Photonics) as shown in (a), and the light source as shown in (b). The delay circuit can control the delay of the synchronization signal in steps of 0.25 nanoseconds from the computer. Using this system, we obtain 256 images by changing the delay from 0 to 64 nanoseconds by 0.25 nanosecond.

Fig. 4
figure 4

Experimental setting. a ToF camera. b Light source. c Delay circuit with controller

Full size image

We measure a scene including wood, wax, soap, and books as shown in Fig. 5 a, which cause subsurface scattering and inter-reflection. Firstly, we confirm the effectiveness of the non-negative constraint. Figure 5 b shows raw observation corresponding to the low-resolution PSF, estimated PSF by pseudo-inverse, and estimated PSF by non-negative least squares of the wood region. The observation is quite blurred due to the convolution effect of the width of light pulse and exposure time. While pseudo-inverse estimation is suffered from over-fitting, estimated result with the non-negative constraint is much stable, and the effectiveness of non-negative constraint is confirmed. We plot the PSFs of different translucent objects in the detailed scale in Fig. 5 c. Different shapes of PSFs are estimated, especially for wood and the others. The temporal spread of light inside translucent objects is finished within 1 ns; hence, our sub-nanosecond estimation has an advantage. Secondly, we recover the temporal PSFs for all pixels as shown in Fig. 5 d. This is so called “light-in-flight” images and shows how the light propagates in the scene.

Fig. 5
figure 5

Experimental result. a The scene. b Comparison with observation and recovered results with and without non-negative constraint. Non-negativity contributes to the stability. c PSFs of different translucent objects. Different shape of PSFs are recovered. d Recovered PSFs for all pixels, as known as light-in-flight images. Light propagation of the scene can be visualized

Full size image

We have also measured another scene as shown in Fig. 6 a, which consists of a screen and a mannequin. Figure 6 b shows a captured image with a long exposure, where the reflection on the screen and the mannequin are superposed. Recovered PSF is shown in Fig. 6 c, where two peaks are confirmed. Separating two peaks of the estimated PSF, reflections on the screen and mannequin can be separated as shown in Fig. 6 d, e. Two layers are clearly separated, and the sharpened image of the mannequin is obtained.

Fig. 6
figure 6

Result for a mannequin’s scene. a A mannequin is placed behind the screen. b Normal photo. Reflections on the mannequin and the screen is superposed. c Recovered PSF has two peaks. d Foreground reflection components (reflection of the screen). e Background components (Reflection on the mannequin)

Full size image

4 Conclusion

In this paper, we propose a method to estimate temporal PSF with high temporal resolution by combining a pulsed ToF camera and a simple delay circuit. We use the delay circuit to control the light emission timing in units of sub-nanoseconds and recover the temporal PSF by a least square method with the non-negative constraint. We have conducted some real-world experiments and confirm the effectiveness of our method. We are interested in improving more higher temporal resolution in the future research so that the subsurface scattering can be fully analyzed.