Defogging Kinect: Simultaneous Estimation of Object Region and Depth in Foggy Scenes

TL;DR

This paper presents a novel method using a ToF camera to simultaneously estimate object regions and depth in foggy scenes, overcoming challenges posed by participating media for 3D reconstruction.

Contribution

The method uniquely leverages scattering saturation and light attenuation phenomena to estimate scattering and depth simultaneously in foggy environments.

Findings

Effective depth and object region estimation in fog using Kinect v2

Robust IRLS optimization improves estimation accuracy

Validated with real fog scenes and synthesized data

Abstract

Three-dimensional (3D) reconstruction and scene depth estimation from 2-dimensional (2D) images are major tasks in computer vision. However, using conventional 3D reconstruction techniques gets challenging in participating media such as murky water, fog, or smoke. We have developed a method that uses a time-of-flight (ToF) camera to estimate an object region and depth in participating media simultaneously. The scattering component is saturated, so it does not depend on the scene depth, and received signals bouncing off distant points are negligible due to light attenuation in the participating media, so the observation of such a point contains only a scattering component. These phenomena enable us to estimate the scattering component in an object region from a background that only contains the scattering component. The problem is formulated as robust estimation where the object region is regarded as outliers, and it enables the simultaneous estimation of an object region and depth on the basis of an iteratively reweighted least squares (IRLS) optimization scheme. We demonstrate the effectiveness of the proposed method using captured images from a Kinect v2 in real foggy scenes and evaluate the applicability with synthesized data.

Tables

Table 1. Table 1: Mean depth error on each object under different density conditions. w/o method denotes error before applying the proposed method.

		Plane	Chair	Desk	Hand	Duck
Figure 9(a)	w/o method	117.26 mm	198.79 mm	459.73 mm	425.21 mm	574.28 mm
(thin)	proposed	17.97 mm	65.01 mm	83.40 mm	32.04 mm	45.08 mm
Figure 8	w/o method	253.38 mm	372.02 mm	656.11 mm	669.64 mm	798.17 mm
(medium)	proposed	20.67 mm	60.27 mm	106.68 mm	50.63 mm	118.21 mm
Figure 9(b)	w/o method	421.31 mm	531.48 mm	798.85 mm	844.32 mm	953.56 mm
(thick)	proposed	19.95 mm	70.63 mm	202.36 mm	75.79 mm	143.20 mm

Equations

$$\alpha e^{j\phi }\in \mathbb{C}.$$

$$z=\frac{c\phi }{4\pi f},$$

$$\tilde{\alpha }(u,v;\mathbf{p})e^{j\tilde{\phi }(u,v;\mathbf{p})}={\alpha }_d(u,v;\mathbf{p})e^{j{\phi }_d(u,v;\mathbf{p})}$$

$$+{\int }_{\Vert \mathbf{x}\Vert =\Vert {\mathbf{x}}_0(u,v)\Vert }^{\Vert \mathbf{p}\Vert }\alpha (u,v;\mathbf{x})e^{j\phi (u,v;\mathbf{x})}d\Vert \mathbf{x}\Vert ,$$

$$\Vert \mathbf{x}\Vert \ge \Vert {\mathbf{x}}_{saturate}\Vert \Rightarrow \alpha (u,v;\mathbf{x})=0.$$

$$\tilde{\alpha }(u,v;\mathbf{p})e^{j\tilde{\phi }(u,v;\mathbf{p})}={\alpha }_d(u,v;\mathbf{p})e^{j{\phi }_d(u,v;\mathbf{p})}$$

$$+\underset{={\alpha }_s(u,v)e^{j{\phi }_s(u,v)}}{\underbrace{{\int }_{\Vert \mathbf{x}\Vert =\Vert {\mathbf{x}}_0(u,v)\Vert }^{\Vert {\mathbf{x}}_{saturate}\Vert }\alpha (u,v;\mathbf{x})e^{j\phi (u,v;\mathbf{x})}d\Vert \mathbf{x}\Vert }},$$

$${\alpha }_d(u,v;{\mathbf{p}}_{far})=0.$$

$$\tilde{\alpha }(u,v;{\mathbf{p}}_{far})e^{j\tilde{\phi }(u,v;{\mathbf{p}}_{far})}={\alpha }_s(u,v)e^{j{\phi }_s(u,v)}.$$

$${\alpha }_d=\sqrt{(\tilde{\alpha }\cos \tilde{\phi }-{\alpha }_s\cos {\phi }_s{)}^2+(\tilde{\alpha }\sin \tilde{\phi }-{\alpha }_s\sin {\phi }_s{)}^2},$$

$${\phi }_d=\arg \left((\tilde{\alpha }\cos \tilde{\phi }-{\alpha }_s\cos {\phi }_s)+j(\tilde{\alpha }\sin \tilde{\phi }-{\alpha }_s\sin {\phi }_s)\right),$$

$$x_k(u,v)$$

$$\underset{\mathbf{x},{\mathbf{a}}_1,\cdots ,{\mathbf{a}}_K}{\min }$$

$$+{\gamma }_1\sum _{k=1}^K\Vert \mathbf{U}{\mathbf{a}}_k-{\mathbf{x}}_k{\Vert }^2$$

$$+{\gamma }_2\Vert \mathbf{Fx}-\mathbf{x}{\Vert }^2$$

$$+{\gamma }_3\Vert \nabla \mathbf{x}{\Vert }^2.$$

$$\underset{\mathbf{x},{\mathbf{a}}_1,\cdots ,{\mathbf{a}}_K}{\min }$$

$$+{\gamma }_1^{\prime }\sum _{k=1}^K\Vert \mathbf{U}{\mathbf{a}}_k-{\mathbf{x}}_k{\Vert }^2$$

$$+{\gamma }_2^{\prime }\Vert \mathbf{Fx}-\mathbf{x}{\Vert }^2$$

$$+{\gamma }_3^{\prime }\Vert \nabla \mathbf{x}{\Vert }^2,$$

$$w_i=\frac{{\rho }^{\prime }\left((x_i-{\tilde{x}}_i)/{\sigma }_1\right)}{(x_i-{\tilde{x}}_i)/{\sigma }_1}.$$

$$\rho(x)=\left\{\begin{array}[]{ll}\frac{c^{2}}{6}\left[1-\left\{1-\left(\frac{x}{c}\right)^{2}\right\}^{3}\right]&if\;\;|x|\leq c\\ \frac{c^{2}}{6}&otherwise.\end{array}\right.$$

$$w_{i}=\left\{\begin{array}[]{ll}\left\{1-\left(\frac{r_{i}}{c}\right)^{2}\right\}^{2}&if\;\;|r_{i}|\leq c\\ 0&otherwise,\end{array}\right.$$

$$\underset{\mathbf{x},{\mathbf{a}}_1,\cdots ,{\mathbf{a}}_K}{\min }$$

$$+{\gamma }_1\sum _{k=1}^K\Vert \mathbf{U}{\mathbf{a}}_k-{\mathbf{x}}_k{\Vert }^2$$

$$+{\gamma }_2\Vert \mathbf{Fx}-\mathbf{x}{\Vert }^2$$

$$+{\gamma }_3\Vert \nabla \mathbf{x}{\Vert }^2.$$

$$\frac{\partial }{\partial \mathbf{x}}\sum _{k=1}^K\rho \left(\frac{\Vert {\mathbf{x}}_k-{\tilde{\mathbf{x}}}_k\Vert }{{\sigma }_2}\right)$$

$$=\sum _{k=1}^K{\mathcal{F}}_k\left(\frac{1}{{\sigma }_2^2}{w}_k({\mathbf{x}}_k-{\tilde{\mathbf{x}}}_k)\right)$$

$$=\frac{1}{{\sigma }_2^2}\mathbf{W}(\mathbf{x}-\tilde{\mathbf{x}}),$$

$$w_k=\frac{{\rho }^{\prime }(\Vert {\mathbf{x}}_k-{\tilde{\mathbf{x}}}_k\Vert /{\sigma }_2)}{\Vert {\mathbf{x}}_k-{\tilde{\mathbf{x}}}_k\Vert /{\sigma }_2}.$$

$${\sigma }_1$$

$${\sigma }_2$$

$${\alpha }_d(u,v)=e^{-2\beta d(u,v)}\hat{\alpha }(u,v),$$

$$\beta =\frac{1}{|\Omega |}\sum _{(u,v)\in \Omega }\frac{\log \hat{\alpha }(u,v)-\log {\alpha }_d(u,v)}{2d(u,v)},$$

$$P(\theta )=\frac{1}{4\pi }\frac{1-g^2}{(1+g^2-2g\cos \theta {)}^{3/2}},$$

$${\alpha }_s(z)e^{j{\phi }_s(z)}={\int }_{z_0}^z\frac{1}{z^2}\beta P(\pi )e^{-2\beta z}{e}^{j\frac{4\pi f}{c}z}dz.$$

$${\alpha }_d(z)e^{j{\phi }_d(z)}=\frac{I}{z^2}{e}^{-2\beta z}{e}^{j\frac{4\pi f}{c}z},$$

Full text

1 Introduction

Three-dimensional (3D) reconstruction and scene depth estimation from 2-dimensional (2D) images are important tasks in computer vision. These techniques can be utilized for a variety of applications such as robot vision and self-driving vehicles. However, using 3D reconstruction techniques gets challenging in participating media such as murky water, fog, or smoke. As shown at the bottom left of Fig. 1, the contrast of a captured image in participating media is reduced because light propagating through the media gets attenuated and the received signal contains not only reflected light from the object surface but also scattered light due to suspended particles. These effects make it difficult to use conventional 3D reconstruction techniques such as structure from motion or shape from shading.

In response to the issue above, we have developed a method that enables us to acquire 3D geometry in participating media. Specifically, we use a time-of-flight (ToF) camera to detect an object and estimate the distance in participating media (see Fig. 2). The proposed method enables simultaneous automatic obstacle detection and depth reconstruction in challenging environments such as disaster sites where light is scattered.

There are several architectures for ToF cameras. We use Microsoft Kinect v2, a continuous-wave ToF camera that emits a modulated sinusoid signal into a scene and then measures the amplitude of light that bounces off an object surface and the phase shift between the illumination and received signal. These observations are represented as an amplitude image and a phase image as shown in Fig. 2. Since the phase shift depends on an optical path, we can reconstruct the depth from the phase shift. In this study, we denote the observation of an object surface by a direct component.

This architecture assumes that each camera pixel observes a single point in a scene. As mentioned previously, however, the observed signal in participating media includes a scattering component due to light scattering as well as a direct component. The amplitude and phase shift suffer from the scattering effect, and this causes error of depth measurement (see Fig. 1).

We aim to recover the direct component, and this is an ill-posed problem because the separation of two components is required. To deal with this problem, we leverage the saturation of a scattering component and light attenuation in participating media. Given a near light source in participating media, a scattering component is saturated close to the camera (Treibitz and Schechner, 2009; Tsiotsios et al., 2014). This means the scattering component does not depend on an object in the scene if the object is located at a certain distance from the camera. Moreover, the intensity of light propagating through participating media is attenuated exponentially relative to distance. Thus, reflected light from a distant point is negligible and the observation of such a point includes only a scattering component. In this paper, we assume that a target scene consists of an object region and a background that only contains a scattering component. This scattering component can then be estimated simply by observing the background. In addition to the above assumption, we introduce two priors to estimate the scattering component: first, a scattering component can be approximated using a quadratic function in a local image patch (local quadratic prior), and second, a scattering component has a symmetrical characteristic in an overall image (global symmetrical prior).

The estimation of a scattering component is formulated as robust estimation, where the object region is regarded as outliers because it contains a direct component. We propose an optimization scheme based on iteratively reweighted least squares (IRLS) (Holland and Welsch, 1977; Fox and Weisberg, 2002; Chartrand and Yin, 2008; Wipf and Nagarajan, 2010), which minimizes weighted least squares iteratively. The object region can then be extracted via the IRLS weight as outliers.

In section 2 of this paper, we briefly overview previous studies on computer vision application in participating media. In section 3, an image formation model in participating media for ToF measurement is introduced. In section 4, we describe the proposed method for simultaneous estimation of the object region and depth, and in section 5, we present the experimental results. We conclude the paper in section 6 with a breif summary and mention of future work.

2 Related work

2.1 Scattering removal

In participating media, the contrast of a captured image is reduced due to scattered light and light attenuation. In computer vision and image processing, several methods have been proposed to recover an original image from a degraded one (He et al., 2011; Nishino et al., 2012; Fattal, 2014; Berman et al., 2016). These methods are referred to as defogging or dehazing, and the problem is known to be ill-posed. Many priors have been introduced to address the ill-posed nature of the problem. For example, He et al. (2011) and Berman et al. (2016) respectively proposed a dark channel prior and a haze-line prior. Recently, many learning-based methods using neural networks have also been proposed (Cai et al., 2016; Ren et al., 2016; Li et al., 2017; Zhang and Patel, 2018; Yang and Sun, 2018).

2.2 3D reconstruction in participating media

Our goal is to reconstruct a 3D scene in participating media. There have been several works on the design of 3D reconstruction techniques for participating media. For example, Li et al. (2015) proposed simultaneous dehazing and 3D reconstruction using multi-view stereo. Photometric stereo in participating media has been proposed by Tsiotsios et al. (2014); Murez et al. (2017); Fujimura et al. (2018), in which light transport models were built in participating media for photometric stereo application. Other studies have been based on structured light (Narasimhan et al., 2005; Gu et al., 2013), light field (Tian et al., 2017), and the absorption of infrared light in underwater scenes (Asano et al., 2016).

2.3 Multipath interference of ToF

A ToF camera assumes that each camera pixel observes a single point in a scene. In participating media, however, the measurement also includes scattered light. This problem is known as multipath interference (MPI). MPI is caused not just by light scattering in participating media but also by subsurface scattering or interreflection in common scenes. Thus, many previous studies have tackled MPI compensation (Fuchs, 2010; Freedman et al., 2014; Naik et al., 2015; Kasambi et al., 2016; Guo et al., 2018).

In this paper, we limit our focus to MPI caused by light scattering in participating media. ToF measurement in participating media has been proposed by Heide et al. (2014); Satat et al. (2018). Heide et al. (2014) developed a scattering model based on exponentially modified Gaussians for transient imaging using a photonic mixer device (PMD) (Heide et al., 2013). Satat et al. (2018) demonstrated that scattered photons observed with a single photon avalanche diode (SPAD) have gamma distribution and leveraged this observation to separate received photons into a directly reflected component and a scattering component. Our method differs from these approaches in that we just use an off-the-shelf ToF camera (Kinect v2) with no special hardware modification.

3 ToF observation in participating media

In this section, we describe our image formation model for a continuous-wave ToF camera in participating media. As with many previous studies (Narasimhan et al., 2005, 2006; Treibitz and Schechner, 2009; Tsiotsios et al., 2014), we assume here that forward scattering and multiple scattering are negligible and that the density of a participating medium in a scene is homogeneous.

A continuous-wave ToF camera illuminates a scene with amplitude-modulated light and then measures the amplitude of received signal $\alpha$ and phase shift $\varphi$ between the illumination and received signal. This observation can be described using a phasor (Gupta et al., 2015), as

[TABLE]

Since the phase shift is proportional to the depth of an object, we can compute the depth as

[TABLE]

where $z$ is depth, $c$ is the speed of light, and $f$ is the modulation frequency of the camera.

In participating media, the observation contains scattered light. As shown in Fig. 3, light interacts with the medium on the line of sight and then arrives at the camera pixel. Thus, the observed scattering component is the sum of scattered light on the line of sight. Now, we consider the 3D coordinate, the origin of which is the camera center. When the camera observes a surface point $\mathbf{p}\in\mathbb{R}^{3}$ at a camera pixel $(u,v)$, the total observation $\tilde{\alpha}(u,v;\mathbf{p})e^{j\tilde{\varphi}(u,v;\mathbf{p})}$ can be written as

[TABLE]

where $\alpha_{d}(u,v;\mathbf{p})$ and $\varphi_{d}(u,v;\mathbf{p})$ are the direct components. $\alpha_{d}(u,v;\mathbf{p})$ depends on the surface albedo, shading, and attenuation, which is caused by the medium as well as the inverse square law. $\alpha(u,v;\mathbf{x})e^{j\varphi(u,v;\mathbf{x})}$ is the observation of scattered light at a position $\mathbf{x}$. Note that although the scattering component can be written using an integral, the domain of the integral (red line in Fig. 3) depends on the relative position between the light source and camera pixel. This is because an ideal point light source irradiates a scene with isotropic intensity, while a practical illumination such as a spotlight has a limited beam angle (Tsiotsios et al., 2014).

Assuming a near light source in participating media, an observed scattering component is saturated close to the camera (Treibitz and Schechner, 2009; Tsiotsios et al., 2014). That is, there exists $\mathbf{x}_{saturate}$ for which

[TABLE]

Therefore, we can rewrite Eq. (3) as

[TABLE]

where $\alpha_{s}(u,v)$ and $\varphi_{s}(u,v)$ are the scattering components, which depend on only the camera pixel $(u,v)$ rather than the object depth.

Although the observation consists of the direct component $\alpha_{d}(u,v;\mathbf{p})e^{j\varphi_{d}(u,v;\mathbf{p})}$ and the scattering component $\alpha_{s}(u,v)e^{j\varphi_{s}(u,v)}$, the attenuation due to the medium reduces the direct component dramatically. Thus, if the camera observes a distant point $\mathbf{p}_{far}$, the amplitude of the reflected light fades away, that is,

[TABLE]

Therefore, the observation of the distant point includes only a scattering component:

[TABLE]

Figure 4 shows amplitude and phase images when the camera observes a black surface in a foggy scene. The intensity of reflected light from the black surface is very small, so this approximates a distant observation where only a scattering component can be observed. As discussed above, in both the amplitude and phase images, the scattering component is inhomogeneous because the illumination has a limited beam angle.

The measurement range of our method is between a saturation point and a background point that has no direct component. More details about the saturation of the scattering component and the measurement range are provided in section 5.4.

4 Simultaneous estimation of object region and depth

As explained in the previous section, a scattering component depends on the position of a camera pixel rather than a target object. In addition, only the scattering component is observed in the background where an object is farther away. Thus, our goal is to estimate the scattering component in an object region from the observation of the background. After estimating scattering components $\alpha_{s}(u,v)$ and $\varphi_{s}(u,v)$ at each pixel, we compute the amplitude and phase shift of a direct component from Eq. (5):

[TABLE]

where an operator ${\rm arg}$ returns the argument of a complex number. Then, depth is recovered substituting the phase into Eq. (2).

In this section, we describe how our method divides camera pixels into an object region and a background, and simultaneously estimates the scattering component in the object region. First, we introduce two priors to estimate the scattering component, and then the problem is formulated as robust estimation, which allows us to extract the object region as outliers. In the following, with a slight alteration of notation, we refer to both an amplitude image and a phase image as an image, since we process both images in the same manner.

4.1 Prior of scattering component

We can estimate the scattering component of an object region from a background because the component does not depend on the object. Tsiotsios et al. (2014) approximated backscatter as a quadratic function in a captured image. Similarly to their work, we also introduce priors, local quadratic prior and global symmetrical prior, that allow us to estimate the scattering component.

Local quadratic prior.

In our ToF setting, we found that a scattering component cannot be approximated globally with a simple function. Thus, as shown in Fig. 5, we assume that a scattering component can be represented with a quadratic function in a local image patch, that is,

[TABLE]

where $x_{k}(u,v)$ is the value at a pixel $(u,v)$ in a local image patch $k$. $\mathbf{u}=[u^{2}\;\;uv\;\;v^{2}\;\;u\;\;v\;\;1]^{\top}$ is a 6-dimensional vector and $\mathbf{a}_{k}=[a_{1}^{k}\;\;a_{2}^{k}\;\;a_{3}^{k}\;\;a_{4}^{k}\;\;a_{5}^{k}\;\;a_{6}^{k}]^{\top}$ denotes the coefficients of the quadratic function in patch $k$.

Global symmetrical prior.

However, this local prior is not useful when there exists a large object region and a quadratic function is also fitted into the values in that region. To address this problem, we introduce a global prior to the scattering component.

As discussed in section 3, a scattering component depends on the relative position between a camera pixel and a light source. This is because the individual starting points of the integral in Eq. (3) differ from each other. Meanwhile, as shown in Fig. 6, we assume that the camera and light source are collocated on the line that is parallel to the horizontal axis of the image. Kinect v2 has this setting, and other devices can easily be built on the basis of this setting. In this case, the integral domain of a pixel is consistent with that of the symmetrical pixel with respect to the central axis of the image. Thus, the observed scattering component also has symmetry, and we leverage this symmetry as a global prior.

4.2 Formulation as robust estimation

We formulate the scattering component estimation problem as robust estimation. Specifically, we solve the following optimization problem:

[TABLE]

The first term of Eq. (4.2) is a data term where $\tilde{\mathbf{x}}=[\tilde{x}_{i}\;\cdots\;\tilde{x}_{N}]^{\top}$ and $\mathbf{x}=[x_{i}\;\cdots\;x_{N}]^{\top}$ are a captured image and a scattering component, respectively. $N$ is the number of camera pixels, and $\sigma_{1}$ is a scale parameter. We use a nonlinear differentiable function $\rho(x)$ rather than square error $x^{2}$, which allows us to make the estimation robust against outliers. In this study, we simply use the residual of the observation and the scattering component as the data term, i.e., pixels that contain a direct component are regarded as outliers.

We use three additional regularization terms. The second term represents the local prior. $K$ is the number of patches for local quadratic function fitting. $\mathbf{U}$ is an $N_{k}\times 6$ matrix where $N_{k}$ is the number of pixels in patch $k$ and each row of $\mathbf{U}$ is a vector $\mathbf{u}$ that corresponds to each pixel coordinate. In this study, these patches do not overlap each other. The third term represents the global prior where $\mathbf{F}\in\mathbb{R}^{N\times N}$ is a matrix that flips an image vertically. The last term is a smoothing term where $\nabla$ denotes a gradient operator. This smoothing accelerates the optimization. Hyperparameters $\gamma_{1},\gamma_{2},\gamma_{3}$ control the contribution of each term.

4.3 IRLS and object region estimation

We minimize Eq. (4.2) with respect to a scattering component $\mathbf{x}$ and the coefficients of quadratic functions $\mathbf{a}_{1},\cdots,\mathbf{a}_{K}$. However, the nonlinearity of $\rho(x)$ makes it difficult to obtain a closed-form solution. For efficient computation, the IRLS optimization was developed in the literature (Holland and Welsch, 1977; Fox and Weisberg, 2002). IRLS minimizes weighted least squares iteratively and the weight is updated using the current estimate in each iteration. The objective function in Eq. (4.2) is transformed into weighted least squares as follows:

[TABLE]

where $\mathbf{W}={\rm diag}(\mathbf{w})$ is an $N\times N$ matrix and $\mathbf{w}=[w_{1},\cdots,w_{N}]^{\top}$ is the weight for each error $x_{i}-\tilde{x}_{i}$. Hyperparameters are given as $\gamma^{\prime}_{*}=2\sigma_{1}^{2}\gamma_{*}$. Equation (4.3) is quadratic with respect to the scattering component $\mathbf{x}$, and thus is easy to optimize. In each iteration, we solve Eq. (4.3) for $\mathbf{x}$ and $\mathbf{a}_{1},\cdots,\mathbf{a}_{K}$, and the weight can be updated using the current estimate as

[TABLE]

The specific update rule of the weight depends on the nonlinear function $\rho(x)$. In this study, we use the following function as $\rho(x)$:

[TABLE]

This function yields the following update:

[TABLE]

where $r_{i}=(x_{i}-\tilde{x}_{i})/\sigma_{1}$, and $c$ is a tuning parameter. This update is referred to as Tukey’s biweight (Beaton and Tukey, 1974; Fox and Weisberg, 2002), where $0\leq w_{i}\leq 1$.

The weight controls the robust estimation, that is, a large error term reduces the corresponding weight. In this study, we consider the object region as outliers, and thus the weight in the object region should be small. Therefore, we can leverage the IRLS weight to extract the object region from the image.

4.4 Coarse-to-fine optimization

The accurate object region extraction is critical for the effectiveness of the scattering component estimation. In section 4.1, we introduced the local and global priors of the scattering component to deal with a large object region. To make the region extraction more robust, we developped a coarse-to-fine optimization scheme. Before solving Eq. (4.2), we optimize the following objective function:

[TABLE]

This is similar to the patch-based robust regression proposed by Chaudhury and Singer (2013). The difference from Eq. (4.2) is that the data term consists of patch-wise errors. Equation (4.4) can be transformed into IRLS as well as Eq. (4.2), and the IRLS weight is updated patch-wise rather than pixel-wise. Differentiating the first term of Eq. (4.4) with respect to $\mathbf{x}$, we can obtain

[TABLE]

where an operator $\mathcal{F}_{k}:\mathbb{R}^{N_{k}}\to\mathbb{R}^{N}$ embeds an input patch into an overall image with zero padding and returns its vectorized form. The weight matrix is $\mathbf{W}={\rm diags}(\mathbf{w})={\rm diags}([w_{1}\mathbbm{1}^{\top}_{1}\cdots w_{K}\mathbbm{1}^{\top}_{K}]^{\top})$ where $\mathbbm{1}_{k}\in\mathbb{R}^{N_{k}}$. The weight $w_{k}$ for each patch $k$ is given as

[TABLE]

Therefore, we can obtain the same objective function as Eq. (4.3) where $\gamma^{\prime}_{*}=2\sigma_{2}^{2}\gamma_{2}$.

Algorithm 1 shows the overall procedure of the simultaneous estimation of a scattering component and an object region. We first solve Eq. (4.4) for the weight in a patch level and then solve (4.2) in a pixel level. Each scale parameter is computed only at the first iteration and is fixed during subsequent iterations. We compute the scale parameters using a median absolute deviation, which is the robust measure of a deviation, as follows (Fox and Weisberg, 2002):

[TABLE]

At the end of the algorithm, we binarize the IRLS weight to generate an object mask. This procedure is applied to an amplitude and a phase image in the same manner, and thus we can obtain the object mask in each domain. In this study, we determine a final object mask as their intersection.

5 Experiments

5.1 Controlled scene

Experimental setup.

We first evaluated the effectiveness of the proposed method in a controlled environment. The experimental environment is shown in Fig. 7. We set up a fog generator and a Kinect v2 in a closed space sized 186 $\times$ 161 cm with black walls and floor. The observation of the wall includes only a scattering component because incident light into the wall is absorbed. The Kinect v2 has three modulation frequencies: 120, 80, and 16 MHz. We used images obtained with 16 MHz. To acquire an amplitude and phase image, we used the source code given by Tanaka et al. (2017). Their code provides the average image of several frames, and we modified the code so that only a single frame was input. To compensate for high frequency noise, we used a bilateral filter as preprocessing.

The spatial resolution of an image captured by Kinect v2 is 424 $\times$ 512 pixels. We divided a captured image into 4 $\times$ 4 patches ($K=16$) for local quadratic prior. In section 4.1, we assumed that the camera and the light source are collocated on a line that runs parallel to the horizontal axis of the image. In practice, the camera and light source are slightly out of alignment. Although this violates the symmetry of the scattering component, we found that error due to this misalignment is negligibly small. In our implementation, we defined $\mathbf{F}$ as a matrix that flips an image with respect to the 200th row of the image. In addition, we did not use the 24 rows of the lower part of the image for the third term of Eq. (4.2), as these pixels have no information of global symmetrical prior. For amplitude images, we set the hyperparameters of the objective function as $[\gamma^{\prime}_{1},\gamma^{\prime}_{2},\gamma^{\prime}_{3}]=[0.1,0.1,10]$, and the tuning parameter of the funcion $\rho(x)$ is set as $c=4,7$ in the coarse and fine level optimization, respectively. For phase images, we set $[\gamma^{\prime}_{1},\gamma^{\prime}_{2},\gamma^{\prime}_{3}]=[0.01,0.1,50]$ and $c=2,3$.

Experimental result.

The results are shown in Fig. 8. Figure 8(a) shows the foggy scene, which has five target objects: “plane,” “chair,” “desk,” “hand,” and “duck.” Figure 8(b)(c) show the estimation of a scattering component and a object region for the amplitude and phase image, respectively. The object region depicted here is the IRLS weight before binarization. As shown, the proposed method effectively estimated the object region via the weight. Of particular note is that thin regions such as the legs of the chair could be extracted. The estimation of the scattering component in the object region was also successful. Figure 8(d) shows the results of the depth reconstruction, where we show the depth measurement without and with fog, and the reconstructed depth. The depth measurement in the foggy scene had large error here due to fog. In contrast, the proposed method reconstructed the object depth correctly.

We tested the proposed method under different density conditions. Figure 9 shows the results under thin fog and thick fog. In a highly foggy scene like the one in Fig. 9(b), the accuracy of depth reconstruction was reduced where a scattering component had a large effect, but the proposed method could estimate an object region and improve the depth measurement regardless of medium density. The mean depth error on each object with and without our method in scenes under different density conditions is listed in Table 1; here, we define the ground truth as the measured depth without fog. As shown, the proposed method could reduce the error significantly regardless of fog density.

5.2 Realistic scene

Next, we tested the proposed method in a more realistic scene. Figure 10(a) shows the target scene. We artificially generated fog in the same manner as the controlled experimental setting, although the scene in Fig. 10 has neither dark walls nor floor. Note that the scene has materials with various types of reflectance, including a lamp made from paper, a glossy vase, and a wooden ornament. The results of the scattering component and object region for the amplitude and phase image are shown in Fig. 10(b)(c), respectively, and the estimation of the depth reconstruction is shown in 10(d). The results of this experiment showed that the proposed method could also extract the object region and improve the depth measurement in a scene that has a general background.

5.3 Experiments with synthesized data

Synthesized data.

To investigate the effectiveness in more varied scenes, we evaluated the proposed method with synthesized data. The procedure of generating the synthesized data is shown in Fig. 11. We assume that a scattering component does not depend on object depth, and thus we observed a direct component and a scattering component separately and then combined them into a synthesized image.

To synthesize an image, we have to know the scattering coefficient in a scene. First, we captured a foggy scene that includes calibration objects, and the region of the calibration objects was masked manually. After that, we compensated for the defective region by solving Eq. (4.3) to estimate the scattering component. The weight $\mathbf{w}$ in Eq. (4.3) corresponds to the mask, and we only used the regularization of the global symmetrical prior and smoothing prior. Using the estimated scattering component, we can compute the direct component of the amplitude using Eq. (8). Now, the relationship between the amplitude without fog and the attenuated direct component in a foggy scene is given as

[TABLE]

where $\hat{\alpha}(u,v)$ is the amplitude at a pixel $(u,v)$ without fog, $d(u,v)$ is the distance between the camera and calibration object, and $\beta$ is a scattering coefficient. Therefore, we computed the scattering coefficient as

[TABLE]

where $\Omega$ denotes the set of pixels in the mask and $|\Omega|$ is the number of pixels in $\Omega$. In the scene in Fig. 11, a scattering coefficient was computed as $\beta=3.5\times 10^{-4}$ /mm.

We also observed a scene without fog, which was used for the direct component after being attenuated by the scattering coefficient. We combined the attenuated signal and the scattering component to synthesize amplitude and phase images.

Experimental result.

The results are shown in Figs. 13, 13, and 15. In each figure, (b)(c) show the estimated scattering component and the IRLS weight for the amplitude and phase image. We show the output of both the coarse and fine level optimization. (d) shows the result of the depth reconstruction. In each scene, the proposed method effectively extracted the object region and estimated the scattering component. (e) shows the result of just applying the fine optimization. In the scenes in Figs. 13 and 15, performing only the fine-level optimization failed to detect the accurate object region, and this deteriorated the scattering component estimation. In contrast, the coarse-to-fine approach made the region extraction more robust.

We show a failure case in Fig. 15. In a scene that has a large object region, our method was less effective because a quadratic function also fits to values in the object region. In Fig. 15, a large textureless object region exists on the left side. In addition, the global symmetrical prior did not work in this region because the object occupied the pixels from top to bottom in the image.

5.4 Discussion of measurement range

We assume that a scattering component is saturated close to a camera and there exists a background that has only a scattering component. In this section, we discuss the measurement capability of our method in this context.

As shown in Fig. 16, a scene has a saturation point and a background point denoted as $\mathbf{x}_{saturate}$ and $\mathbf{x}_{background}$. For simplicity, the camera and light source are assumed to be collocated in the same place. Behind $\mathbf{x}_{saturate}$, a scattering component is constant due to its saturation, while a direct component fades away behind $\mathbf{x}_{background}$. Therefore, the measurement range of our method is between $\mathbf{x}_{saturate}$ and $\mathbf{x}_{background}$.

We simulated the measurement range to evaluate the capability. Similarly to the process of synthesizing images, a scattering coefficient was computed for the scene in Fig. 8 to use for the simulation ($\beta=3.2\times 10^{-4}$ /mm). We use the following Henyey-Greenstein phase function for scattering property:

[TABLE]

where $\theta$ is a scattering angle. The parameter $g$ was set as $0.9$ for fog (Narasimhan and Nayar, 2003). A scattering component from a camera to depth $z$ is given as

[TABLE]

We set the starting point of the integral as $z_{0}=10$ mm. A direct component from depth $z$ is computed as

[TABLE]

where $I$ consists of a surface albedo and shading, and we set $I=1$ in this simulation. The total observation $\tilde{\alpha}(z)e^{\tilde{\varphi}(z)}$ is the sum of these components as with Eq. (3).

Figure 17 shows the simulation results. In (a) and (b), the horizontal axis denotes depth $z$ and the vertical axis denotes a scattering component for amplitude and phase. These figures validate the saturation characteristic. Meanwhile, in 17(c) and (d), the vertical axes denote the residual of the observation and scattering component. $\Delta\varphi$ is given by the residual angle of $\tilde{\alpha}(z)e^{j\varphi(z)}$ and $\alpha_{s}(z)e^{j\varphi_{s}(z)}$ on the complex plane. These values represent the remaining direct components from depth $z$.

Now, we can set $\|\mathbf{x}_{background}\|=2500$ mm and $5000$ mm for amplitude and phase from Fig. 17(c) and (d) because the direct component is close to zero. In contrast, in Fig. 17(a) and (b), if we set $\|\mathbf{x}_{saturate}\|=1000$ mm, the estimation error of the scattering component for amplitude due to the saturation assumption can be considered almost zero because $1-\alpha_{s}(1000)/\alpha_{s}(8000)\approx 0$, and for phase, the error is $1-\varphi_{s}(1000)/\varphi_{s}(8000)\approx 6.0\%$.

In the experiments with real data, all of the target objects were located between $1000$ mm and $2000$ mm. If we assume the measurement range between $\|\mathbf{x}_{saturate}\|=1000$ mm and $\|\mathbf{x}_{background}\|=2500$ mm, the target objects are located in that range, and from above discussion, we have just $6.0\%$ error of the scattering component estimation for phase due to the saturation assumption. As shown in the experiments, we can effectively reconstruct the object depth regardless of the error. The measurement range depends on the density of a participating medium, and it will get larger under thinner fog.

6 Conclusion

In this paper, we proposed a method that simultaneously estimates an object region and depth by using a continuous-wave ToF camera. We leveraged the saturation of a scattering component and the attenuation of a direct component from a distant point in a scene. The formulation with a robust estimator and the IRLS optimization scheme allow us to estimate the scattering component and object region simultaneously.

The limitation of the proposed method is that we assume the density of the participating medium in the scene to be homogeneous; thus, it cannot be applied to an inhomogeneous medium or a dynamic floating medium. In these environments, global symmetrical prior does not hold. In addition, we assume that a scene has a background region, which makes it difficult to apply the method to scenes filled with objects. In future work, we will address these problems in order to further enhance the real-world applicability of the proposed method.

Acknowledgements.

This work was supported by JSPS KAKENHI Grant Number 18H03263.