Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Single scattering models for radiative transfer of isotropic and cone-shaped light sources in fog

Open Access Open Access

Abstract

The simulation of rare edge cases such as adverse weather conditions is the enabler for the deployment of the next generation of autonomous drones and vehicles into conditions where human operation is error-prone. Therefore, such settings must be simulated as accurately as possible and be computationally efficient, so to allow the training of deep learning algorithms for scene understanding, which require large-scale datasets disallowing extensive Monte Carlo simulations. One computationally-expensive step is the simulation of light sources in scattering media, which can be tackled by the radiative transfer equation and approximated by analytical solutions in the following. Traditionally, a single scattering event is assumed for fog rendering, since it is the dominant effect for relatively low scattering media. This assumption allows us to present an improved solution to calculate the so called air-light integral that can be evaluated fast and robustly for an isotropic point source in homogeneous media. Additionally, the solution is extended for a cone-shaped source and implemented in a computer vision rendering pipeline fulfilling computational restrictions for deep learning uses. All solutions can handle arbitrary azimuthally symmetric phase functions and were tested with the Henyey-Greenstein phase function and an advection fog phase function calculated from a particle distribution using Mie’s theory. The used approximations are validated through extensive Monte Carlo simulations and the solutions are used to augment good weather images towards inclement conditions with focus on visible light sources, so to provide additional data in such hard-to-collect settings.

© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

Today’s computer vision methods enable autonomous vehicles and drones to reason about their environment. A key factor behind their success are the underlying datasets used for training, enabling to learn complex neural networks able to distill information from raw camera captures [13]. However, rare events like driving in adverse weather are underrepresented in traditional datasets [4]. Nonetheless, enabling the operation in such conditions would lead to beneficial safety gains [4]. Typical approaches to close the data gap are based on the idea of augmenting clear weather scenes to match real world inclement weather recordings [57]. Thereby, fog is rendered through empirical models applying Koschmieder’s model and extracting physical parameters from the clear reference. One of such assumptions is the constant ambient airlight, only true for daytime fog [5,6] and which does not hold anymore in the presence of active light sources. Extensions as presented in [7] introduce a handling of light sources but model the spread of the light source illumination through scattering, empirically by matching a Gaussian filter to experimental observations.

One approach to model this physically-accurate broadening is the use of Radiative Transport Equations (RTE). Such equations are also known in related theoretical fields, as neutron transport [8]. Its importance for the simulation of foggy scenarios is elaborated for example by Cerezo et al. [9]. But solving RTE equations numerically is not trivial and requires Monte Carlo simulations [10]. Such methods have an important history in many applications from medicine to astrophysics [1114], but they require a high amount of computational power and contain statistical noise in contrast to fast analytical solutions. In our example, the Monte Carlo simulation runs multiple days on graphical processing units, whereas the analytical solution takes only seconds for the single scattering solution. Analytical solutions are difficult to find but are important to improve and accelerate rendering applications containing participating media like fog [15]. Due to this challenges, the full RTE can be reduced to only consider single scattered light. It should be noted that the solution employed to compute the single-scattered radiance is typically expressed in terms of the so-called air–light integral. In fact, it has been used for fog rendering already by Nishita, Miyawaki and Nakamae [16]. Pegoraro and Parker [17] derived an air-light integral representation in isotropically and anisotropically scattering media, which they later improved together with Schott [18]. In general, the role of the RTE for light transport in fog is outlined by Bentz et al. [19]. They also developed an analytical model for weak angular dependence based on the diffusion equation, which is limited to the moderate and high scattering regime. More analytical solutions exist taking into account single scattering [20,21] as well as all scattering events [2224]. In our recent publication [25] we have derived some exact solutions for the single-scattered radiance in semi-infinite media based on radiative transport theory. Further applications of RTE solutions within the computer graphics field can be found in [2630].

In this publication, two models for the single scattering radiance due to a point source located in an infinitely extended scattering medium are presented as they offer a good approximation for perception systems located within the scattering medium. An example image with traffic and street lights in a foggy atmosphere is shown in Fig. 1. These models differ in the angular radiation characteristics of the applied source. The first model applies an isotropic point source that emits light uniformly in all directions. As stated above, the solution to this problem can be expressed in terms of the air–light integral. The numerical evaluation of this integral representation can become a challenging task. Therefore, we present a modified integral representation as an alternative to the classical air-light integral, which can be efficiently evaluated without numerical problems. The second model assumes a cone-shape source, where light is emitted in a cone with opening angle $\theta _0$. To our knowledge, this is the first solution with a cone-shaped radiation characteristic of the source. It includes the isotropic solution in the limiting cases of $\theta _0={180}^{\circ }$. All models are verified using Monte Carlo simulations and the validity of the single scattering approximation is investigated by comparing it to the solution containing all scattering events. Additionally, an application of the derived solution for fog rendering is presented. It is used to improve the rendering of light sources in adverse weather simulations.

 figure: Fig. 1.

Fig. 1. Image of active light sources in foggy atmosphere, taken in Bucharest on January 9, 2022 by fusion-of-horizons [31]. The cone shape illumination pattern of street lanterns, signal and traffic light motivate the found approximation. License: "CC BY 2.0 https://creativecommons.org/licenses/by/2.0/".

Download Full Size | PDF

2. Theory and method

In the following section, the theoretical framework for deriving the analytical solution is presented, which is based on the single scattering solution to the radiative transfer equation according to Chandrasekhar [32]. In detail, subsection 2.1 treats the isotropic point source with angular uniform radiation characteristic and subsection 2.2 the cone source with truncated angular reach.

2.1 Isotropic source

The single scattering radiance at position $\mathbf {r}$ in direction $\hat {\mathbf {s}}$ due to an isotropic point source located at the origin can be derived as [21]

$$L_1^{\textrm{iso}}(\mathbf{r}, \hat{\mathbf{s}}) = \dfrac{\mu_s}{4\pi} \int_0^\infty \dfrac{e^{-\mu_t(l+R)}}{R^2} f\left(\dfrac{r \mu -\ell}{R} \right) \,\mathrm{d}\ell,$$
where $f(\mu )$ is the phase function and $r$, $\mu$ and $R$ are defined as
$$r = \lvert \mathbf{r}\rvert,$$
$$\mu = \dfrac{\mathbf{r}\cdot \hat{\mathbf{s}}}{r} = \cos\alpha,$$
$$R = \lvert\mathbf{r} - \ell \hat{\mathbf{s}}\rvert = \sqrt{\ell^2+r^2-2\ell r\mu} \,.$$
Furthermore, we have $\mu _t=\mu _a+\mu _s$ as the total attenuation coefficient with $\mu _a$ and $\mu _s$ being respectively the absorption and scattering coefficient. Equation (1) is the above-mentioned air-light that we want to improve. This representation is very similar to the approach shown by Pegoraro and Parker [17] but different coordinates are applied. Due to the symmetry of the source, the angular dependency is only on the angle $\alpha$ between the position and the direction $\hat {\mathbf {s}}$, see a schematic view in Fig. 2. The integration path is depicted by the red line. This integral must be computed numerically. Unfortunately, for the most important region of $\mu$ near $1$, the integrand has a sharp peak at $l \approx \mu r$. This is shown in Fig. 3. Simple quadrature algorithms are problematic for such integrals and even more sophisticated ones may not be optimal, due to the sharp peak requiring thousands of evaluation points for small angles $\alpha$. Therefore, the integral is changed to an alternative form using the substitution
$$y = \arccos\left(\dfrac{r \cos\alpha -\ell}{\sqrt{\ell^2+r^2-2 \ell r\cos\alpha}}\right),$$
with
$$\dfrac{\partial y}{\partial \ell} = \dfrac{r\sin\alpha}{\ell^2+r^2-2 \ell r\cos\alpha} = \dfrac{r\sin\alpha}{R^2} \,.$$
The angle $y$ for this substitution is also depicted in Fig. 2. For the phase function argument, we then get
$$\dfrac{r\cos\alpha -\ell}{R} = \cos y \,.$$
Additionally, Mollweide’s formula [33] provides the relation
$$\dfrac{\ell+R}{r} = \dfrac{\cos\left(\frac{y-2\alpha}{2}\right)}{\sin\left(\frac{\pi-y}{2}\right)},$$
which leads to
$$\ell+R = r\cos\alpha + r\sin\alpha\tan\left(\frac{y}{2}\right) \,.$$
Plugging everything into the single scattering integral and adjusting the integral limits accordingly results in the representation
$$L_1^{\textrm{iso}}(r,\alpha) = \mu_s \dfrac{e^{-\mu_t r \cos\alpha}}{4\pi r\sin\alpha} \int_\alpha^\pi e^{-\mu_t r\sin\alpha \tan\left(\frac{y}{2}\right)} f\left(\cos y\right) \,\mathrm{d}y \,.$$
In contrast to the original integrand that exhibits a sharp peak for small angles $\alpha$, this integrand can be integrated even for very small angles. In addition, in the case of small angles, $L_1^{\textrm {iso}}$ can be approximated with
$$L_1^{\textrm{iso}}(r,\theta) \overset{\alpha \ll {1}^{{\circ}}}{\approx} \mu_s \dfrac{e^{-\mu_t r}}{4\pi r\sin\alpha} \int_0^\pi f\left(\cos y\right) \,\mathrm{d}y,$$
for faster computation. The integral $L_1^{\textrm {iso}}(r,\alpha )$ can be mapped to the standard interval $[-1,1]$ resulting in
$$L_1^{\textrm{iso}}(r,\alpha) = \mu_s \dfrac{e^{-\mu_t r \cos\alpha}}{4\pi r\sin\alpha} \frac{\pi-\alpha}{2} \int_{{-}1}^1 e^{-\mu_t r\sin\alpha \tan\left(\frac{\pi+\alpha}{4}+x\frac{\pi-\alpha}{4}\right)} f\left(\cos\left(\frac{\pi+\alpha}{2}+x\frac{\pi-\alpha}{2}\right)\right) \,\mathrm{d}x \,.$$
The integral in Eq. (12) can be solved by double exponential substitution methods [34,35]:
$$x = \tanh\left( 2\sinh(u)\right),$$
$$\mathrm{d}x = \dfrac{2\cosh(u)}{\cosh^2\left(2\sinh(u)\right)} \mathrm{d}u \,. $$
The effect of these substitutions is shown in Fig. 4. The integrals can then be computed by the simple trapezoidal rule with fast numerical convergence.

 figure: Fig. 2.

Fig. 2. Sketch illustrating the coordinates used in the solution for an isotropic point source located at the origin. $\alpha$ is the angle between the position $\mathbf {r}$ and the radiance direction $\hat {\mathbf {s}}$. The angle $y$ is used in the alternative representation.

Download Full Size | PDF

 figure: Fig. 3.

Fig. 3. Plot of the single scattering integrand for different angles using the Henyey-Greenstein phase function. The parameters are $\mu _{\mathrm {s}}={0.08}{\textrm {m}^{-1}}$, $\mu _{\mathrm {a}}={1{\times }10^{-5}}{\textrm {m}^{-1}}$, $g={0.8}$, $r={20}\textrm {m}$, $\mu _{\mathrm {s}} r={1.6}$.

Download Full Size | PDF

 figure: Fig. 4.

Fig. 4. Plot of the transformed single scattering integrand (Eq. (12)) for different angles using the Henyey-Greenstein phase function. The parameters are $\mu _{\mathrm {s}}={0.08}{\textrm {m}^{-1}}$, $\mu _{\mathrm {a}}=1{\times }10^{-5}{\textrm {m}^{-1}}$, $g={0.8}$, $r={20}\textrm {m}$.

Download Full Size | PDF

The derived representation for the isotropic source was verified using Monte Carlo simulations, see Fig. 5. Additionally, the Monte Carlo simulation considering all scattering events $L_{n\geq 1}^{\textrm {iso}}$ is plotted for comparison. A Henyey-Greenstein phase function [36] with an anisotropy factor of $g={0.8}$ was used. The optical properties are chosen to be typical for fog. The scattering coefficient $\mu _{\mathrm {s}}={0.08}{\textrm {m}^{-1}}$ corresponds to a visibility of approximately $\frac {4}{\mu _{\mathrm {s}}}={50}\textrm {m}$, the distance from the source to the detector is $r={20}\textrm {m}$ and, thus, the optical depth is $\tau =\mu _{\mathrm {s}} r={1.6}$. The absorption coefficient $\mu _{\mathrm {a}}={1{\times }10^{-5}}{\textrm {m}^{-1}}$ is very small, since it can be mostly neglected for fog in the visible wavelength range. These optical properties are chosen as an example and will be used further below. It should be noted that all the solutions were also tested with other values and no significant differences were found. For $\alpha \ \to 0$ the solution diverges to $\infty$, but it is integrable over a solid angle $\alpha \ \in \left [0,\,\alpha _c\right ]$, $\gamma \ \in \left [0,\,2 \pi \right ]$ respecting the Jacobian $\sin \alpha$, it does not depend on the azimuth $\gamma$. The inset shows the relative error between the derived solution and the Monte Carlo simulation of the single scattered radiance. In Fig. 6, the calculation was repeated for smaller angles. In the left image the distance is $r={20}\textrm {m}$ giving an optical depth of $\tau =\mu _{\mathrm {s}} r={1.6}$ and in the right image it is $r={1}\textrm {m}$ giving an optical depth of $\tau =\mu _{\mathrm {s}} r={0.08}$. It can be seen that for small angles and low optical depth the single scattering approximation is a very good approximation to the complete solutions to the RTE considering all scattering interactions. Looking at the relative errors, the statistical noise increases for very small angles due to the Monte Carlo simulation. An additional verification with other phase function can be found in the Appendix A.

 figure: Fig. 5.

Fig. 5. Single scattering comparison using the Henyey-Greenstein phase function with $\mu _{\mathrm {s}}={0.08}{\textrm {m}^{-1}}$, $\mu _{\mathrm {a}}={1{\times }10^{-5}}{\textrm {m}^{-1}}$, $g={0.8}$, $r={20}\textrm {m}$, for an isotropic point source at the origin. The inset shows the relative error between the derived solution and the Monte Carlo simulation of the single scattered radiance.

Download Full Size | PDF

 figure: Fig. 6.

Fig. 6. Single scattering comparison using the Henyey-Greenstein phase function with $\mu _{\mathrm {s}}={0.08}{\textrm {m}^{-1}}$, $\mu _{\mathrm {a}}={1{\times }10^{-5}}{\textrm {m}^{-1}}$, $g={0.8}$, for an isotropic point source at the origin for small angles. The insets show the relative error between the analytical solution and the Monte Carlo simulation of the single scattered radiance. Left: $r={20}\textrm {m}$, Right: $r={1}\textrm {m}$.

Download Full Size | PDF

2.2 Cone source

This section deals with the single-scattered radiance due to a cone source of the form

$$S({\mathbf{r}},{ \hat{\mathbf{s}}},\theta_0)=\frac{\delta({\mathbf{r}})}{2\pi(1-\mu_0)}\int_C\delta({ \hat{\mathbf{s}}}'-{ \hat{\mathbf{s}}})\,\mathrm{d}{ \hat{\mathbf{s}}}'^2=\frac{\delta({\mathbf{r}})}{2\pi(1-\mu_0)} \begin{cases} 1,\quad 0\leq \theta\leq\theta_0, \\ 0,\quad \theta_0<\theta\leq\pi, \end{cases}$$
where $C=\{{ \hat {\mathbf {s}}}'\in {\mathbb {S}}^2\,|\,0\leq \theta '\leq \theta _0,\,0\leq \phi '<2\pi \}$ with $\theta _0\in (0,\pi )$ being the opening angle, $\mu _0=\cos \theta _0$ and $\theta$ is the angle between the direction $\hat {\mathbf {s}}$ and the z-axis $\hat {\mathbf {z}}$: $\cos \theta = \hat {\mathbf {s}} \cdot \hat {\mathbf {z}}$. The source is normalized to unity. We additionally note on the limit
$$\mathop {\lim }_{\theta_0 \to 0} S({\mathbf{r}},{ \hat{\mathbf{s}}},\theta_0)= \mathop {\lim }_{\mu_0 \to 1} S({\mathbf{r}},{ \hat{\mathbf{s}}},\mu_0)= \delta({\mathbf{r}})\frac{\delta(1-\cos\theta)}{2\pi}=\delta({\mathbf{r}})\frac{\delta(1-{ \hat{\mathbf{s}}} \cdot { \hat{\mathbf{z}}})}{2\pi}=\delta({\mathbf{r}})\delta({ \hat{\mathbf{s}}}-{ \hat{\mathbf{z}}}),$$
corresponding with a unidirectional point source that radiates along the positive $z$-direction. The coordinates used within the derivations are illustrated in Fig. 7. The intended solution can be derived in different ways. In our recent publication [25], we have obtained the single-scattered radiance for the infinite and semi-infinite medium caused by the unidirectional beam $\delta ({\mathbf {r}}-{\mathbf {r}}_0)\delta ({ \hat {\mathbf {s}}}-{ \hat {\mathbf {s}}}')$ that radiates along an arbitrary direction ${ \hat {\mathbf {s}}}'$. Therefore, we can find the solution for the cone source via integration of Green’s function for the infinite space over $C$. The required representation has the form [25]
$$G_1({\mathbf{r}},{ \hat{\mathbf{s}}},{ \hat{\mathbf{s}}}')=\mu_s f({ \hat{\mathbf{s}}}\cdot { \hat{\mathbf{s}}}') \int_0^\infty \frac{e^{-\mu_t(\ell+R)}}{R^2} \delta\!\left(\frac{{\mathbf{ r}}-\ell { \hat{\mathbf{s}}}}{R}-{ \hat{\mathbf{s}}}'\right)\!\!\,\mathrm{d}\ell,$$
where $\delta (\cdot )$ denotes the Dirac delta function, $r$ and $R=R({\mathbf {r}},{ \hat {\mathbf {s}}},\ell )$ are defined above in Eqs. (2) and (4) and $\hat {\mathbf {r}}=\mathbf {r}/r$. Integration of (17) over $C$ under the use of
$$\int_C f({ \hat{\mathbf{s}}}\cdot { \hat{\mathbf{s}}}')\delta\!\left(\frac{{\mathbf{ r}}-\ell { \hat{\mathbf{s}}}}{R}-{ \hat{\mathbf{s}}}'\right)\!\mathrm{d}{ \hat{\mathbf{s}}}'^2=f\!\left(\frac{{\mathbf{ r}}\cdot{ \hat{\mathbf{s}}} -\ell }{R}\right)\!\Theta\!\left(\frac{z-\ell \cos\theta}{R}-\mu_0\right)$$
and normalization leads to the following single-scattered radiance
$$L_1^{\textrm{cone}}({\mathbf{r}},{ \hat{\mathbf{s}}})=\frac{1}{2\pi(1-\mu_0)}\int_C G_1({\mathbf{r}},{ \hat{\mathbf{s}}},{ \hat{\mathbf{s}}}')\,\mathrm{d}{ \hat{\mathbf{s}}}'^2$$
$$=\frac{\mu_s }{2\pi(1-\mu_0)} \int_0^\infty \frac{e^{-\mu_t(\ell+R)}}{R^2} f\!\left(\frac{{\mathbf{ r}}\cdot{ \hat{\mathbf{s}}} -\ell }{R}\right)\!\Theta\!\left(\frac{z-\ell \cos\theta}{R}-\mu_0\right)\!\!\,\mathrm{d}\ell,$$
with $\Theta (\cdot )$ being the Heaviside step function. Similar to the substitution (5), we perform a variable transformation $\ell :[0,\infty )\to [\vartheta _1,\pi )$ with the following monotonically increasing function
$$\ell \mapsto \vartheta(\ell) := {\rm{arccot}}\left(\frac{{\mathbf{r}}\cdot{ \hat{\mathbf{s}}}-\ell}{r\sqrt{1-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})^2}} \right)$$
and the corresponding inverse $\ell (\vartheta )={\mathbf {r}}\cdot { \hat {\mathbf {s}}}-r\sqrt {1-(\hat {\mathbf {r}}\cdot { \hat {\mathbf {s}}})^2}\cot \vartheta$. Furthermore, we have
$$\frac{\mathrm{d}\ell}{\mathrm{d}\vartheta}=\frac{r\sqrt{1-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})^2}}{\sin^2\vartheta}=\frac{R({\mathbf{r}},{ \hat{\mathbf{s}}},\ell)^2}{r\sqrt{1-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})^2}}$$
and $\cos \vartheta =({\mathbf { r}}\cdot { \hat {\mathbf {s}}} -\ell )/R$. The new lower and upper limit of integration are given by
$$\vartheta_1=\vartheta(0)={\rm{arccot}}\left(\frac{\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}}}{\sqrt{1-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})^2}} \right)=\arccos(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}}),$$
$$\vartheta_2=\vartheta(\infty)={\rm{arccot}}(-\infty)=\pi.$$
The resulting integral representation for the radiance in terms of the new variable is
$$L_1^{\textrm{cone}}({\mathbf{r}},{ \hat{\mathbf{s}}})=\frac{\mu_s\exp[-\mu_t ({\mathbf{r}}\cdot{ \hat{\mathbf{s}}})]}{2\pi (1-\mu_0)r\sqrt{1-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})^2}}\int_I f(\cos\vartheta)\exp\left(-\mu_t r \sqrt{1-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})^2}\tan\frac{\vartheta}{2}\right)$$
$$\times\Theta(\mu_c -\mu_0)\,\mathrm{d}\vartheta,$$
where $I=(\vartheta _1,\pi )$ and the function
$$\mu_c=\mu_c({\mathbf{r}},{ \hat{\mathbf{s}}},\vartheta)= \frac{(z/r-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})\cos\theta)\sin\vartheta+\cos\theta\sqrt{1-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})^2}\cos\vartheta}{\sqrt{1-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})^2}}.$$
In order to eliminate the step function in (25), we have to perform the integration over the restricted domain $D := \{\vartheta \in I\,|\,\mu _c({\mathbf {r}},{ \hat {\mathbf {s}}},\vartheta )-\mu _0 > 0\}$. In the possible case of $D=\emptyset$, we have $L_1^{\textrm {cone}}=0$. We start with the following conversion
$$\mu_c({\mathbf{r}},{ \hat{\mathbf{s}}},\vartheta)-\mu_0 > 0\iff $$
$$(z/r-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})\cos\theta)\sin\vartheta+\cos\theta\sqrt{1-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})^2}\cos\vartheta> \mu_0\sqrt{1-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})^2}.$$
To isolate the integration variable $\vartheta$, we convert the part left of the inequality sign into the form
$$(z/r-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})\cos\theta)\sin\vartheta+\cos\theta\sqrt{1-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})^2}\cos\vartheta=A\sin(\vartheta+\varphi),$$
where the amplitude $A$ and the phase shift $\varphi$ are given by
$$A=\sqrt{(z/r)^2+\cos^2\theta-2(z/r)(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})\cos\theta},$$
$$\varphi=\arg\{z/r-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})\cos\theta+i\cos\theta \sqrt{1-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})^2} \}\in(-\pi,\pi].$$
Thus, the inequality (28) can be written as
$$\sin(\vartheta+\varphi)> \mu_0\frac{\sqrt{1-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})^2}}{A}=: w.$$
Next, we perform a short case analysis. In this context, we make use of the property $(a,b)=\emptyset$ if $a\geq b$. Furthermore, it is sufficient to consider the inequality (33) for values $\vartheta +\varphi \in (-\pi,2\pi )$, because $\vartheta \in I$ and $\varphi \in (-\pi,\pi ]$.

 figure: Fig. 7.

Fig. 7. Sketch illustrating the coordinates used in the solution for the cone-shaped source. The source is located at the origin and pointing in $z$-direction $\hat {\mathbf {z}}$ in a cone-shaped manner with opening angle $\theta _0$. The detection position $\mathbf {r}$ is given in cylindrical coordinates $\rho$, $\phi _l$ and $z$. The angular variables of the radiance in direction $\hat {\mathbf {s}}$ are $\theta$ and $\phi$.

Download Full Size | PDF

Case 1: Opening angles with $\mu _0\geq 0$. In this situation, we have $w\geq 0$. For $w>1$, the inequality (33) cannot be satisfied so that $D=\emptyset$ and hence $L_1^{\textrm {cone}}=0$. If $0\leq w\leq 1$, we define $\beta := \arcsin w\in [0,\pi /2]$ and obtain $\sin (\vartheta +\varphi )>w$ for $\vartheta \in (\beta -\varphi,\pi -\beta -\varphi )$. The intersection with the integration interval results in

$$\vartheta\in(\beta-\varphi,\pi-\beta-\varphi)\cap (\vartheta_1,\pi)=D=(\max\{\beta-\varphi,\vartheta_1\},\min\{\pi-\beta-\varphi,\pi\}).$$

Case 2: Opening angles with $\mu _0< 0$. In this case, we have $w<0$. If $w<-1$, the inequality (33) is satisfied for all $\vartheta \in I$, yielding $D=(\vartheta _1,\pi )$. For $-1\leq w<0$, when $\beta \in [-\pi /2,0)$, a slightly longer analysis shows that

$$\vartheta\in((-\pi-\varphi,-\pi-\beta-\varphi)\cup (\beta-\varphi,\pi-\beta-\varphi) \cup(2\pi+\beta-\varphi,2\pi-\varphi))\cap(\vartheta_1,\pi).$$
The resulting integration domain becomes $D=D_1\cup D_2\cup D_3$, where $D_i$ $(i=1,2,3)$ are disjoint intervals defined by
$$D_1=(\max\{-\pi-\varphi,\vartheta_1\},\min\{-\pi-\beta-\varphi,\pi\}), $$
$$D_2 = (\max\{\beta-\varphi,\vartheta_1\},\min\{\pi-\beta-\varphi,\pi\}), $$
$$D_3 = (\max\{2\pi+\beta-\varphi,\vartheta_1\},\min\{2\pi-\varphi,\pi\}). $$
The integration over this union can therefore be carried out in the form $\int _D(\cdot )\,\mathrm {d}\vartheta =\sum \int _{D_i}(\cdot )\,\mathrm {d}\vartheta$. Both cases discussed above are illustrated in Fig. 8. The blue line segment displays exemplarily the solution of the inequality (33) for the first case, whereas the three red segments correspond with the three disjoint intervals which may appear in the second case. The final form for the single-scattered radiance due to a cone source without the use of the Heaviside step function becomes
$$L_1^{\textrm{cone}}({\mathbf{r}},{ \hat{\mathbf{s}}})=\frac{\mu_s\exp[-\mu_t ({\mathbf{r}}\cdot{ \hat{\mathbf{s}}})]}{2\pi (1-\mu_0)r\sqrt{1-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})^2}}\int_D f(\cos\vartheta)\exp\left(-\mu_t r \sqrt{1-(\hat {\mathbf{r}}\cdot{ \hat{\mathbf{s}}})^2}\tan\frac{\vartheta}{2}\right)\mathrm{d}\vartheta,$$
where the region of integration $D$ has been determined in the case analysis above. As shown above, the cone-shaped source includes the unidirectional point source in the limit $\theta _0\to {0}^{\circ }$. The resulting single-scattered radiance for this source type is given by [25]
$$L_1^{\textrm{uni}}(\mathbf{r}, \hat{\mathbf{s}})=\mathop {\lim }_{\theta_0 \to 0} L_1^{\textrm{cone}}({\mathbf{r}},{ \hat{\mathbf{s}}})=\mu_s f(\cos{\theta}) \frac{\exp(-\mu_t z -\mu_t \rho \tan(\theta/2))}{\rho \sin{\theta}} \delta(\phi_l-\phi) \Theta(z-\rho \cot{\theta}).$$
On the other hand, the corresponding solution for the isotropic source, when $\theta _0\to {180}^{\circ }$, can be found in (12). Again, we verified the derived single-scattered radiance via comparisons with the Monte Carlo method. The optical properties were the same as in the previous comparisons. Two detector positions $\mathbf {r}_1=\left (1,1,1\right )\textrm {m}$ and $\mathbf {r}_2=\left (1,1,-1\right )\textrm {m}$ and different opening angles between ${5}^{\circ }$ and ${180}^{\circ }$ were considered. In Fig. 9, it can be seen that the analytical results agree with the simulations. For smaller opening angles, the radiance becomes narrower and weaker. This is due to the fact that for the unidirectional source, there is no signal for $\phi \neq \phi _l$ due to the delta distribution in Eq. (40). For opening angles smaller than ${180}^{\circ }$, the single scattered radiance can go to zero, depending on the detector position.

 figure: Fig. 8.

Fig. 8. Illustration of possible solutions of the inequality (33). The blue segment corresponds with $\mu _0\geq 0$ and the red ones demonstrate a situation for $\mu _0<0$.

Download Full Size | PDF

 figure: Fig. 9.

Fig. 9. Single scattering comparison using a Henyey-Greenstein phase function with $\mu _{\mathrm {s}}={0.08}{\textrm {m}^{-1}}$, $\mu _{\mathrm {a}}={1{\times }10^{-5}}{\textrm {m}^{-1}}$, $g={0.8}$, $\chi =\lvert \phi _l-\phi \rvert ={30.5}^{\circ }$ for cone-shaped point sources at the origin pointing in $z$-direction. Left: $x={1}\textrm {m}$, $y={1}\textrm {m}$, $z={1}\textrm {m}$, Right: $x={1}\textrm {m}$, $y={1}\textrm {m}$, $z={-1}\textrm {m}$.

Download Full Size | PDF

3. Validity of the single scattering approximation

In this section we investigate the errors of the single scattering approximation compared to the full solution of the RTE considering all scattering interactions. The relative error between the single scattered radiance and the radiance considering all scattering events versus the optical depth $\tau =\mu _s r$ and the angle $\theta$ for $\mu _a={1{\times }10^{-5}}\;\textrm {m}^{-1}$ is shown in Fig. 10. Two types of phase functions were considered, a Henyey-Greenstein phase function with $g={0.8}$ (left) and a phase function for advection fog (right). It was calculated from a model particle size distribution applying Mie’s theory. Thus, it was assumed that the fog particles are spheres. The used phase function is shown in Fig. 11. We emphasize that this phase function is calculated from a particle size distribution that represents a fictitious fog. The particle size distribution is approximated using a modified gamma distribution like the simulation software MODTRAN, with the parameters for advection fog taken from Gebhart et al. [37]. It serves as an example for a general phase function generated from a particle size distribution with Mie’s theory.

 figure: Fig. 10.

Fig. 10. Relative error between the single scattering radiance and the radiance calculated using all scattering events for a Henyey-Greenstein phase function with $g=0.8$ (left) and advection fog [37] (right).

Download Full Size | PDF

 figure: Fig. 11.

Fig. 11. Phase function for advection fog at 550 nm calculated with Mie’s theory from the particle size distribution given by Gebhart et al. [37].

Download Full Size | PDF

The errors are simulated with a GPU-based Monte Carlo simulation using the lookup table based technique described by Naglič et al. [38]. As expected, for small angles and low scattering distances the errors decrease. Already for moderate optical depths of $\mu _{\mathrm {s}} r={2}$, the relative error is below ${10}\%$, even for the highly forward scattering advection fog.

4. Application of the solution: improved light sources in artificial adverse weather images

Here we show how to apply the solution derived above to improve the generation of artificial adverse weather images from real good weather images. The unscattered light of an isotropic light source in an infinitely extended medium is given by

$$I_0=P \frac{\exp(-\mu_t r)}{4 \pi r^2}\delta( \hat{\mathbf{s^\prime}}- \hat{\mathbf{r}}),$$
where $P$ is the source power. For a pixel with solid angle $\sigma$, the directly detected light is
$$Q = \int_\sigma I_0( \hat{\mathbf{s}^\prime}) \mathrm{d}\hat{\mathbf{s}^\prime}^2 = P \frac{\exp(-\mu_{\mathrm{t}} r)}{4 \pi r^2},$$
as long as $\hat {\mathbf {r}}$ is within the angular range of the pixel. From images without fog, $P$ can be estimated to be
$$P = 4 \pi r^2 Q(\mu_{\mathrm{t}}=0)\,.$$
With fog, the direct light can be calculated from the pixels without fog:
$$Q(\mu_{\mathrm{t}} \neq 0) = \exp(-\mu_t r) Q(\mu_{\mathrm{t}}=0) = t Q_{\textrm{clear}},$$
where $t=\exp (-\mu _t r)$, and $Q_{\textrm {clear}}=Q(\mu _{\mathrm {t}}=0)$ is the clear pixel value without fog. The scattered light is onto the pixel is
$$ V= \, P \int_{\sigma} R(\mu_s, \mu_a, r, g, \hat{\mathbf{r}}, \hat{\mathbf{s}^\prime}) \mathrm{d} \hat{\mathbf{s}^\prime}^2 $$
$$\approx \, P R(\mu_s, \mu_a, r, g, \hat{\mathbf{r}}, \hat{\mathbf{s}}) \sigma $$
$$ = \, Q_{\textrm{clear}} 4 \pi r^2 R(\mu_s, \mu_a, r, g, \hat{\mathbf{r}}, \hat{\mathbf{s}}) \sigma =: V_{\textrm{a}},$$
where $R(\mu _s, \mu _a, r, g, \hat {\mathbf {r}},\,\hat {\mathbf {s}})$ is the scattered radiance of a isotropic point source with a power of unity and $\hat {\mathbf {s}}$ is the average direction of the solid angle $\sigma$. Since for small distances or tenuous fog the scattered radiance can be approximated by the single scattered radiance, we set $R=L_1^{\textrm {iso}}$ from Eq. (12). For the images presented below, depth information for every pixel is available using a deep learning stereo depth retrieval algorithm [39]. So the distance from each light source to the camera is known. In previous works, the fog was added to the image by applying the following formula [4]:
$$A_{\textrm{foggy}} = t Q_{\textrm{clear}} + (1-t)M,$$
where $A_{\textrm {foggy}}$ is pixel value for the foggy image and $M$ is the global ambient component, which is retrieved applying the dark channel prior method [40] disallowing the simulation of active light sources in twilight conditions. Contrarily, our approach follows the procedure described by Shi et al. [7] in their Supplemental Material replacing the heuristic light source broadening with a physically accurate model. The process predicts a light source pixel with index $j$. Then each source was multiplied by the transmission coefficient $t$ to consider the attenuation of the unscattered component (Eq. (44)). Afterwards the single scattered light $V_{\textrm {neigh}}$ for each light source pixel onto the neighboring pixels was added. For each pixel $i$ of the image, the value was calculated using Eq. (47) by summing over all pixels $j$ identified as light sources
$$V_{\textrm{neigh},i}=\sum_{\substack{j \in \textrm{ light sources}\\ i\neq j }}V_{\textrm{a}}\left(\mathbf{r}_i, \hat{\mathbf{s}}_{j}\right)\,. $$
In the case $i=j$, the approximation of Eq. (47) cannot be used, since $L_1^{\textrm {iso}}\left ( \mathbf {r}_i, \hat {\mathbf {s}}= \hat {\mathbf {r}}_i\right )=\infty$, but the integral (45) over the solid angle $\sigma$ of a camera pixel can be computed
$$V_{\textrm{self}, i}=V\left( \mathbf{r}_i\right)\,.$$
Thus, the light on the camera pixels $i$ of the light sources in fog is
$$W_{\textrm{light},i} = \begin{cases} V_{\textrm{neigh},i} + (t Q_{\textrm{clear},i} + V_{\textrm{self}, i}) & \text{if}\; i\in \text{light sources}\\ V_{\textrm{neigh},i} & \text{else}\,. \end{cases} $$
The final image is calculated via
$$B_{\textrm{foggy}} = t Q_{\textrm{clear, no ls}} + (1-t)M + W_{\textrm{light}} $$
$$= t\left(Q_{\textrm{clear, no ls}} +Q_{\textrm{clear, ls}}\right) + (1-t)M + V_{\textrm{self, ls}} + V_{\textrm{neigh}} $$
$$ = t Q_{\textrm{clear}} + (1-t)M + V_{\textrm{self, ls}} + V_{\textrm{neigh}},$$
where $Y_{y\textrm {, ls}}$ is the subset of the image identified as light sources and $Y_{y\textrm {, no ls}}$ is the remaining image, $Y_y\in \{Q_{\textrm {clear}},\, V_{\textrm {self}} \}$. Comparing with Eq. (48), we find the same formula plus the single scattering terms ($P_{1\textrm {, ls}}$ and $P_{\textrm {neigh}}$). The solution is applied to an example image in Fig. 12 from the EU-project DENSE, for which depth information for the light sources exists. The applied advection fog is characterized by the same phase function already used above, see Fig. 11. The optical parameters are $\mu _{\mathrm {s}}={0.0287}{\textrm {m}^{-1}}$, $\mu _{\mathrm {a}}={3{\times }10^{-8}}{\textrm {m}^{-1}}$ which corresponds to a visibility of $\frac {4}{\mu _{\mathrm {s}}}={139}\textrm {m}$ and to an optical depth at the distance of the traffic lights ($\approx {35}\textrm {m}$) of $\tau \approx {1.0}$. They are directly derived from the particle distribution. One can clearly see the difference between the approaches. Previously, the effect of fog around the light sources was added by applying a Gaussian filter with a standard deviation of 12 pixels [7]. The distance to the light source, the fog density and the type of fog were not regarded. In contrast to the simple Gaussian filter, the single scattering solution takes the distance of the light source, the fog density and the fog type into account. As can be seen in the example, this leads to a very different appearance of the light source in a foggy environment. With Monte Carlo simulations we also calculated the appearance of the traffic lights considering all scattering events using $L_{n\geq 1}^{\textrm {iso}}$. This is shown in Fig. 12(g). Only the small section was calculating due to the high computation time of the Monte Carlo simulation. Comparing it to the single scattering solution (f), they show very similar results. This indicates that the single scattering solution in these conditions is already a very good approximation. An example with a different phase function is shown in the Appendix A.

 figure: Fig. 12.

Fig. 12. Fog is added to a clear image (a) using a Gaussian filter (c) and the single scattering solution for advection fog (e). Images (b) to (f) zoom in on the section with the traffic lights, marked with the red rectangles in the full images. In image (g) all scattering events are considered with Monte Carlo simulations.

Download Full Size | PDF

5. Summary and conclusion

In summary this work presents an analytical solution for the RTE in case of single scattering for an isotropic and a cone-shaped point source. Here, we assumed the location within a homogeneous infinitely extended scattering medium. To proof the findings, model predictions were validated using Monte Carlo simulations for small radiation angles and equivalent thin optical densities. Further, the solution of the isotropic point source was used to improve the appearance of simulated traffic light sources in dense foggy conditions. For cone-shaped light sources as street lanterns and known opening angles, the introduced solution derived in Eq. (54) could be used. Fog exists in different forms, depending on, e.g., the formation, location and environment, thus many densities and phase functions are possible. Depending on the size of the droplets and the applied wavelength the phase function can represent, both, strongly forward scattering as well as more uniformly scattering behavior. The droplet size distribution also changes during the fog life cycle [41]. We use the selected phase function as an example for a general phase function generated from a particle distribution and using Mie’s theory. The chosen particle distribution does not necessarily coincide with that of real fog, which can be much more complex. The wavelength is also only an example from the visible spectrum [42]. The presented method is not restricted to this particular fog type or wavelength, it can model any phase function. To generate robust datasets, we recommend simulating several fog types. We envision that these results can be used in future work to foster the development of deep learning methods robust to challenging illumination conditions as the one in Fig. 1, through an approach similar to [5,6]. In fact, the extension to challenging illumination conditions includes not only environmental conditions like advertising sign or street lighting, but more importantly also signal lights such as turn signals, traffic lights as shown in Fig. 12, stop lights, warning lights and blue lights of emergency vehicles, which are indispensable in complex urban traffic scenarios. All of those lightning conditions are currently modeled using time consuming Monte Carlo approaches prohibiting the widespread roll out of such of complex simulations due to time constraint issues. The here-proposed analytical solutions offer a reduction of computational time of some orders of magnitude, allowing to process millions of frames as required by state-of-the-art deep learning approaches [43,44]. We further indicate that the single scattering approximation limits the application of the method to certain circumstances. For very dense fog it may not be appropriate. In the future, the solution could be further expanded to include double scattering to reduce this limitation. Extending the solution to consider arbitrary angular radiation characteristics of the source will allow the modeling of vehicle headlamps with known cone illumination profile and differing characteristics for the US or European market.

Appendix

Since the phase function changes with the fog type and wavelength, we additionally show calculations for further phase functions, plotted in Fig. 13. They are again calculated using a modified gamma distribution to model the particle size distribution with parameters from [37]. As one can see, the advection fog is much more forward scattering than the convection fog. The single scattering radiance highly depends on the fog type as shown in Fig. 14. It is calculated using the analytical solution for all four phase functions from Fig. 13 for an isotropic point source at a distance of $r={30}\textrm {m}$. It agrees well to the Monte Carlo simulations that are shown for comparison. Figure 15 shows part of the images from Fig. 12 with an additional image (d), where the single scattering solution with the phase function from convection fog for a wavelength of 550 nm was used. A different fog type would also lead to a different visibility, but to compare only the influence of the phase function, we simulated the image with the same visibility of ${139}\textrm {m}$. Since the convection fog phase function represents a more isotropically scattering behavior, it has a very different effect on the light scattered from the traffic light compared to the advection fog phase function.

 figure: Fig. 13.

Fig. 13. Phase function for advection and convection fog at 550 nm and 1550 nm calculated with Mie’s theory from the particle size distribution given by Gebhart et al. [37].

Download Full Size | PDF

 figure: Fig. 14.

Fig. 14. Single scattering radiance using different fog types at a distance of $r={30}\textrm {m}$ and two different wavelengths, for an isotropic point source at the origin.

Download Full Size | PDF

 figure: Fig. 15.

Fig. 15. Comparison of the different fog types. Images (a) to (c) are already presented in Fig. 12. The additional image (d) show the same section, where the phase function of convection fog at 550 nm is used instead of the one from advection fog.

Download Full Size | PDF

Funding

Bundesministerium für Bildung und Forschung (16ME0344).

Acknowledgments

The research leading to these results is part of the AI-SEE project, which is a co-labelled PENTA and EURIPIDES 2 project endorsed by EUREKA. Co-funding is provided by the following national funding authorities: Austrian Research Promotion Agency (FFG), Business Finland, Federal Ministry of Education and Research (BMBF) and National Research Council of Canada Industrial Research Assistance Program (NRC-IRAP).

Disclosures

The authors declare that there are no conflicts of interest related to this article.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. A. Geiger, P. Lenz, and R. Urtasun, “Are we ready for autonomous driving? the kitti vision benchmark suite,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, (2012), pp. 3354–3361.

2. P. Sun, H. Kretzschmar, X. Dotiwalla, A. Chouard, V. Patnaik, P. Tsui, J. Guo, Y. Zhou, Y. Chai, B. Caine, V. Vasudevan, W. Han, J. Ngiam, H. Zhao, A. Timofeev, S. Ettinger, M. Krivokon, A. Gao, A. Joshi, Y. Zhang, J. Shlens, Z. Chen, and D. Anguelov, “Scalability in perception for autonomous driving: Waymo open dataset,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, (2020), pp. 2443–2451.

3. T.-Y. Lin, M. Maire, S. Belongie, J. Hays, P. Perona, D. Ramanan, P. Dollar, and C. L. Zitnick, “Coco: Commom objects in context,” in Proceedings of the European Conference on Computer Vision (ECCV), D. Fleet, T. Pajdla, B. Schiele, and T. Tuytelaars, eds. (Springer International PublishingCham, 2014), pp. 740–755.

4. M. Bijelic, F. Mannan, T. Gruber, W. Ritter, K. Dietmayer, and F. Heide, “Seeing Through Fog Without Seeing Fog: Deep Sensor Fusion in the Absence of Labeled Training Data,” arXiv:1902.08913 [cs] (2019).

5. S. S. Halder, J.-F. Lalonde, and R. d. Charette, “Physics-based rendering for improving robustness to rain,” in Proceedings of the International Conference on Computer Vision (ICCV), (2019).

6. C. Sakaridis, D. Dai, and L. Van Gool, “Semantic foggy scene understanding with synthetic data,” Int. J. Comput. Vis. 126(9), 973–992 (2018). [CrossRef]  

7. Z. Shi, E. Tseng, M. Bijelic, W. Ritter, and F. Heide, “ZeroScatter: Domain Transfer for Long Distance Imaging and Vision through Scattering Media,” in 2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), (2021), pp. 3475–3485.

8. J. C. Stewart, I. Kuščer, and N. J. McCormick, “Equivalence of special models in energy-dependent neutron transport and nongrey radiative transfer,” Ann. Phys. 40(2), 321–333 (1966). [CrossRef]  

9. E. Cerezo, F. Pérez, X. Pueyo, F. J. Seron, and F. X. Sillion, “A survey on participating media rendering techniques,” The Vis. Comput. 21(5), 303–328 (2005). [CrossRef]  

10. L. L. House and L. W. Avery, “The Monte Carlo technique applied to radiative transfer,” J. Quant. Spectrosc. Radiat. Transf. 9(12), 1579–1591 (1969). [CrossRef]  

11. S. A. Prahl, “A Monte Carlo model of light propagation in tissue,” (Berlin, Germany, 1989), p. 1030509.

12. S. Flock, M. Patterson, B. Wilson, and D. Wyman, “Monte Carlo modeling of light propagation in highly scattering tissues. I. Model predictions and comparison with diffusion theory,” IEEE Trans. Biomed. Eng. 36(12), 1162–1168 (1989). [CrossRef]  

13. C. Zoller and A. Kienle, “Fast and precise image generation of blood vessels embedded in skin,” J. Biomed. Opt. 24(01), 1 (2019). [CrossRef]  

14. U. M. Noebauer and S. A. Sim, “Monte Carlo radiative transfer,” Living Rev. Comput. Astrophys. 5(1), 1 (2019). [CrossRef]  

15. T. Ritschel, C. Dachsbacher, T. Grosch, and J. Kautz, “The State of the Art in Interactive Global Illumination,” Comput. Graph. Forum 31(1), 160–188 (2012). [CrossRef]  

16. T. Nishita, Y. Miyawaki, and E. Nakamae, “A shading model for atmospheric scattering considering luminous intensity distribution of light sources,” in Proceedings of the 14th annual conference on Computer graphics and interactive techniques, (Association for Computing MachineryNew York, NY, USA, 1987), SIGGRAPH ’87, pp. 303–310.

17. V. Pegoraro and S. G. Parker, “An Analytical Solution to Single Scattering in Homogeneous Participating Media,” Comput. Graph. Forum 28(2), 329–335 (2009). [CrossRef]  

18. V. Pegoraro, M. Schott, and S. G. Parker, “A Closed-Form Solution to Single Scattering for General Phase Functions and Light Distributions,” Comput. Graph. Forum 29(4), 1365–1374 (2010). [CrossRef]  

19. B. Z. Bentz, B. Z. Bentz, B. Z. Bentz, B. J. Redman, J. D. v. d. Laan, K. Westlake, A. Glen, A. L. Sanchez, and J. B. Wright, “Light transport with weak angular dependence in fog,” Opt. Express 29(9), 13231–13245 (2021). [CrossRef]  

20. M. L. Shendeleva, “Single-scattering solutions to radiative transfer in infinite turbid media,” J. Opt. Soc. Am. A 30(11), 2169–2174 (2013). [CrossRef]  

21. B. Sun, R. Ramamoorthi, S. G. Narasimhan, and S. K. Nayar, “A practical analytic single scattering model for real time rendering,” ACM Trans. Graph. 24(3), 1040–1049 (2005). [CrossRef]  

22. K. M. Case and P. F. Zweifel, Linear transport theory (Addison-Wesley, Reading, Mass., 1967).

23. J. J. Duderstadt and W. R. Martin, Transport theory (Wiley, 1979).

24. M. Machida, “How to Construct Three-Dimensional Transport Theory Using Rotated Reference Frames,” J. Comput. Theor. Transp. 45(7), 594–609 (2016). [CrossRef]  

25. A. Liemert, S. Geiger, and A. Kienle, “Solutions for single-scattered radiance in the semi-infinite medium based on radiative transport theory,” JOSA A 38(3), 405–411 (2021). [CrossRef]  

26. J. R. Frisvad, T. Hachisuka, and T. K. Kjeldsen, “Directional Dipole Model for Subsurface Scattering,” ACM Trans. Graph. 34(1), 1–12 (2014). [CrossRef]  

27. R. Habel, P. H. Christensen, and W. Jarosz, “Photon Beam Diffusion: A Hybrid Monte Carlo Method for Subsurface Scattering,” Comput. Graph. Forum 32(4), 27–37 (2013). [CrossRef]  

28. R. Ramamoorthi and P. Hanrahan, “A signal-processing framework for inverse rendering,” in Proceedings of the 28th annual conference on Computer graphics and interactive techniques, (Association for Computing Machinery, New York, NY, USA, 2001), SIGGRAPH ’01, pp. 117–128.

29. W. Jakob, E. d’Eon, O. Jakob, and S. Marschner, “A comprehensive framework for rendering layered materials,” ACM Trans. Graph. 33(4), 1–14 (2014). [CrossRef]  

30. H. W. Jensen and P. H. Christensen, “Efficient simulation of light transport in scenes with participating media using photon maps,” in Proceedings of the 25th annual conference on Computer graphics and interactive techniques, (Association for Computing Machinery, New York, NY, USA, 1998), SIGGRAPH ’98, pp. 311–320.

31. fusion-of-horizons, “cavernous,” https://www.flickr.com/photos/fusion_of_horizons/51839030680/ (2022).

32. S. Chandrasekhar, Radiative transfer (Dover Publications, 1960).

33. H. A. DeKleine, “Proof Without Words: Mollweide’s Equation,” Math. Mag. 61(5), 281 (1988). [CrossRef]  

34. H. Takahasi and M. Mori, “Double exponential formulas for numerical integration,” Publ. Res. Inst. for Math. Sci. 9(3), 721–741 (1973). [CrossRef]  

35. L. N. Trefethen and J. A. C. Weideman, “The Exponentially Convergent Trapezoidal Rule,” SIAM Rev. 56(3), 385–458 (2014). [CrossRef]  

36. L. G. Henyey and J. L. Greenstein, “Diffuse radiation in the Galaxy,” The Astrophys. J. 93, 70–83 (1941). [CrossRef]  

37. M. Gebhart, E. Leitgeb, S. Sheikh Muhammad, B. Flecker, C. Chlestil, M. Al Naboulsi, F. de Fornel, and H. Sizun, “Measurement of Light attenuation in dense fog conditions for FSO applications,” (San Diego, California, USA, 2005), p. 58910K.

38. P. Naglic, F. Pernuš, B. Likar, and M. Bürmen, “Lookup table-based sampling of the phase function for Monte Carlo simulations of light propagation in turbid media,” Biomed. Opt. Express 8(3), 1895–1910 (2017). [CrossRef]  

39. J.-R. Chang and Y.-S. Chen, “Pyramid stereo matching network,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, (2018), pp. 5410–5418.

40. K. He, J. Sun, and X. Tang, “Single image haze removal using dark channel prior,” in 2009 IEEE Conference on Computer Vision and Pattern Recognition, (2009), pp. 1956–1963.

41. M. Mazoyer, F. Burnet, and C. Denjean, “Experimental study on the evolution of droplet size distribution during the fog life cycle,” Atmos. Chem. Phys. 22(17), 11305–11321 (2022). [CrossRef]  

42. R. Nebuloni, “Empirical relationships between extinction coefficient and visibility in fog,” Appl. Opt. 44(18), 3795–3804 (2005). [CrossRef]  

43. T. Sun, M. Segu, J. Postels, Y. Wang, L. Van Gool, B. Schiele, F. Tombari, and F. Yu, “SHIFT: a synthetic driving dataset for continuous multi-task domain adaptation,” in Computer Vision and Pattern Recognition, (2022).

44. E. Wood, T. Baltrušaitis, C. Hewitt, S. Dziadzio, T. J. Cashman, and J. Shotton, “Fake it till you make it: face analysis in the wild using synthetic data alone,” in 2021 IEEE/CVF International Conference on Computer Vision (ICCV), (2021), pp. 3661–3671.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (15)

Fig. 1.
Fig. 1. Image of active light sources in foggy atmosphere, taken in Bucharest on January 9, 2022 by fusion-of-horizons [31]. The cone shape illumination pattern of street lanterns, signal and traffic light motivate the found approximation. License: "CC BY 2.0 https://creativecommons.org/licenses/by/2.0/".
Fig. 2.
Fig. 2. Sketch illustrating the coordinates used in the solution for an isotropic point source located at the origin. $\alpha$ is the angle between the position $\mathbf {r}$ and the radiance direction $\hat {\mathbf {s}}$. The angle $y$ is used in the alternative representation.
Fig. 3.
Fig. 3. Plot of the single scattering integrand for different angles using the Henyey-Greenstein phase function. The parameters are $\mu _{\mathrm {s}}={0.08}{\textrm {m}^{-1}}$, $\mu _{\mathrm {a}}={1{\times }10^{-5}}{\textrm {m}^{-1}}$, $g={0.8}$, $r={20}\textrm {m}$, $\mu _{\mathrm {s}} r={1.6}$.
Fig. 4.
Fig. 4. Plot of the transformed single scattering integrand (Eq. (12)) for different angles using the Henyey-Greenstein phase function. The parameters are $\mu _{\mathrm {s}}={0.08}{\textrm {m}^{-1}}$, $\mu _{\mathrm {a}}=1{\times }10^{-5}{\textrm {m}^{-1}}$, $g={0.8}$, $r={20}\textrm {m}$.
Fig. 5.
Fig. 5. Single scattering comparison using the Henyey-Greenstein phase function with $\mu _{\mathrm {s}}={0.08}{\textrm {m}^{-1}}$, $\mu _{\mathrm {a}}={1{\times }10^{-5}}{\textrm {m}^{-1}}$, $g={0.8}$, $r={20}\textrm {m}$, for an isotropic point source at the origin. The inset shows the relative error between the derived solution and the Monte Carlo simulation of the single scattered radiance.
Fig. 6.
Fig. 6. Single scattering comparison using the Henyey-Greenstein phase function with $\mu _{\mathrm {s}}={0.08}{\textrm {m}^{-1}}$, $\mu _{\mathrm {a}}={1{\times }10^{-5}}{\textrm {m}^{-1}}$, $g={0.8}$, for an isotropic point source at the origin for small angles. The insets show the relative error between the analytical solution and the Monte Carlo simulation of the single scattered radiance. Left: $r={20}\textrm {m}$, Right: $r={1}\textrm {m}$.
Fig. 7.
Fig. 7. Sketch illustrating the coordinates used in the solution for the cone-shaped source. The source is located at the origin and pointing in $z$-direction $\hat {\mathbf {z}}$ in a cone-shaped manner with opening angle $\theta _0$. The detection position $\mathbf {r}$ is given in cylindrical coordinates $\rho$, $\phi _l$ and $z$. The angular variables of the radiance in direction $\hat {\mathbf {s}}$ are $\theta$ and $\phi$.
Fig. 8.
Fig. 8. Illustration of possible solutions of the inequality (33). The blue segment corresponds with $\mu _0\geq 0$ and the red ones demonstrate a situation for $\mu _0<0$.
Fig. 9.
Fig. 9. Single scattering comparison using a Henyey-Greenstein phase function with $\mu _{\mathrm {s}}={0.08}{\textrm {m}^{-1}}$, $\mu _{\mathrm {a}}={1{\times }10^{-5}}{\textrm {m}^{-1}}$, $g={0.8}$, $\chi =\lvert \phi _l-\phi \rvert ={30.5}^{\circ }$ for cone-shaped point sources at the origin pointing in $z$-direction. Left: $x={1}\textrm {m}$, $y={1}\textrm {m}$, $z={1}\textrm {m}$, Right: $x={1}\textrm {m}$, $y={1}\textrm {m}$, $z={-1}\textrm {m}$.
Fig. 10.
Fig. 10. Relative error between the single scattering radiance and the radiance calculated using all scattering events for a Henyey-Greenstein phase function with $g=0.8$ (left) and advection fog [37] (right).
Fig. 11.
Fig. 11. Phase function for advection fog at 550 nm calculated with Mie’s theory from the particle size distribution given by Gebhart et al. [37].
Fig. 12.
Fig. 12. Fog is added to a clear image (a) using a Gaussian filter (c) and the single scattering solution for advection fog (e). Images (b) to (f) zoom in on the section with the traffic lights, marked with the red rectangles in the full images. In image (g) all scattering events are considered with Monte Carlo simulations.
Fig. 13.
Fig. 13. Phase function for advection and convection fog at 550 nm and 1550 nm calculated with Mie’s theory from the particle size distribution given by Gebhart et al. [37].
Fig. 14.
Fig. 14. Single scattering radiance using different fog types at a distance of $r={30}\textrm {m}$ and two different wavelengths, for an isotropic point source at the origin.
Fig. 15.
Fig. 15. Comparison of the different fog types. Images (a) to (c) are already presented in Fig. 12. The additional image (d) show the same section, where the phase function of convection fog at 550 nm is used instead of the one from advection fog.

Equations (54)

Equations on this page are rendered with MathJax. Learn more.

L 1 iso ( r , s ^ ) = μ s 4 π 0 e μ t ( l + R ) R 2 f ( r μ R ) d ,
r = | r | ,
μ = r s ^ r = cos α ,
R = | r s ^ | = 2 + r 2 2 r μ .
y = arccos ( r cos α 2 + r 2 2 r cos α ) ,
y = r sin α 2 + r 2 2 r cos α = r sin α R 2 .
r cos α R = cos y .
+ R r = cos ( y 2 α 2 ) sin ( π y 2 ) ,
+ R = r cos α + r sin α tan ( y 2 ) .
L 1 iso ( r , α ) = μ s e μ t r cos α 4 π r sin α α π e μ t r sin α tan ( y 2 ) f ( cos y ) d y .
L 1 iso ( r , θ ) α 1 μ s e μ t r 4 π r sin α 0 π f ( cos y ) d y ,
L 1 iso ( r , α ) = μ s e μ t r cos α 4 π r sin α π α 2 1 1 e μ t r sin α tan ( π + α 4 + x π α 4 ) f ( cos ( π + α 2 + x π α 2 ) ) d x .
x = tanh ( 2 sinh ( u ) ) ,
d x = 2 cosh ( u ) cosh 2 ( 2 sinh ( u ) ) d u .
S ( r , s ^ , θ 0 ) = δ ( r ) 2 π ( 1 μ 0 ) C δ ( s ^ s ^ ) d s ^ 2 = δ ( r ) 2 π ( 1 μ 0 ) { 1 , 0 θ θ 0 , 0 , θ 0 < θ π ,
lim θ 0 0 S ( r , s ^ , θ 0 ) = lim μ 0 1 S ( r , s ^ , μ 0 ) = δ ( r ) δ ( 1 cos θ ) 2 π = δ ( r ) δ ( 1 s ^ z ^ ) 2 π = δ ( r ) δ ( s ^ z ^ ) ,
G 1 ( r , s ^ , s ^ ) = μ s f ( s ^ s ^ ) 0 e μ t ( + R ) R 2 δ ( r s ^ R s ^ ) d ,
C f ( s ^ s ^ ) δ ( r s ^ R s ^ ) d s ^ 2 = f ( r s ^ R ) Θ ( z cos θ R μ 0 )
L 1 cone ( r , s ^ ) = 1 2 π ( 1 μ 0 ) C G 1 ( r , s ^ , s ^ ) d s ^ 2
= μ s 2 π ( 1 μ 0 ) 0 e μ t ( + R ) R 2 f ( r s ^ R ) Θ ( z cos θ R μ 0 ) d ,
ϑ ( ) := a r c c o t ( r s ^ r 1 ( r ^ s ^ ) 2 )
d d ϑ = r 1 ( r ^ s ^ ) 2 sin 2 ϑ = R ( r , s ^ , ) 2 r 1 ( r ^ s ^ ) 2
ϑ 1 = ϑ ( 0 ) = a r c c o t ( r ^ s ^ 1 ( r ^ s ^ ) 2 ) = arccos ( r ^ s ^ ) ,
ϑ 2 = ϑ ( ) = a r c c o t ( ) = π .
L 1 cone ( r , s ^ ) = μ s exp [ μ t ( r s ^ ) ] 2 π ( 1 μ 0 ) r 1 ( r ^ s ^ ) 2 I f ( cos ϑ ) exp ( μ t r 1 ( r ^ s ^ ) 2 tan ϑ 2 )
× Θ ( μ c μ 0 ) d ϑ ,
μ c = μ c ( r , s ^ , ϑ ) = ( z / r ( r ^ s ^ ) cos θ ) sin ϑ + cos θ 1 ( r ^ s ^ ) 2 cos ϑ 1 ( r ^ s ^ ) 2 .
μ c ( r , s ^ , ϑ ) μ 0 > 0
( z / r ( r ^ s ^ ) cos θ ) sin ϑ + cos θ 1 ( r ^ s ^ ) 2 cos ϑ > μ 0 1 ( r ^ s ^ ) 2 .
( z / r ( r ^ s ^ ) cos θ ) sin ϑ + cos θ 1 ( r ^ s ^ ) 2 cos ϑ = A sin ( ϑ + φ ) ,
A = ( z / r ) 2 + cos 2 θ 2 ( z / r ) ( r ^ s ^ ) cos θ ,
φ = arg { z / r ( r ^ s ^ ) cos θ + i cos θ 1 ( r ^ s ^ ) 2 } ( π , π ] .
sin ( ϑ + φ ) > μ 0 1 ( r ^ s ^ ) 2 A =: w .
ϑ ( β φ , π β φ ) ( ϑ 1 , π ) = D = ( max { β φ , ϑ 1 } , min { π β φ , π } ) .
ϑ ( ( π φ , π β φ ) ( β φ , π β φ ) ( 2 π + β φ , 2 π φ ) ) ( ϑ 1 , π ) .
D 1 = ( max { π φ , ϑ 1 } , min { π β φ , π } ) ,
D 2 = ( max { β φ , ϑ 1 } , min { π β φ , π } ) ,
D 3 = ( max { 2 π + β φ , ϑ 1 } , min { 2 π φ , π } ) .
L 1 cone ( r , s ^ ) = μ s exp [ μ t ( r s ^ ) ] 2 π ( 1 μ 0 ) r 1 ( r ^ s ^ ) 2 D f ( cos ϑ ) exp ( μ t r 1 ( r ^ s ^ ) 2 tan ϑ 2 ) d ϑ ,
L 1 uni ( r , s ^ ) = lim θ 0 0 L 1 cone ( r , s ^ ) = μ s f ( cos θ ) exp ( μ t z μ t ρ tan ( θ / 2 ) ) ρ sin θ δ ( ϕ l ϕ ) Θ ( z ρ cot θ ) .
I 0 = P exp ( μ t r ) 4 π r 2 δ ( s ^ r ^ ) ,
Q = σ I 0 ( s ^ ) d s ^ 2 = P exp ( μ t r ) 4 π r 2 ,
P = 4 π r 2 Q ( μ t = 0 ) .
Q ( μ t 0 ) = exp ( μ t r ) Q ( μ t = 0 ) = t Q clear ,
V = P σ R ( μ s , μ a , r , g , r ^ , s ^ ) d s ^ 2
P R ( μ s , μ a , r , g , r ^ , s ^ ) σ
= Q clear 4 π r 2 R ( μ s , μ a , r , g , r ^ , s ^ ) σ =: V a ,
A foggy = t Q clear + ( 1 t ) M ,
V neigh , i = j  light sources i j V a ( r i , s ^ j ) .
V self , i = V ( r i ) .
W light , i = { V neigh , i + ( t Q clear , i + V self , i ) if i light sources V neigh , i else .
B foggy = t Q clear, no ls + ( 1 t ) M + W light
= t ( Q clear, no ls + Q clear, ls ) + ( 1 t ) M + V self, ls + V neigh
= t Q clear + ( 1 t ) M + V self, ls + V neigh ,
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.