by James W. Walker
25 May 2016
Traditional OpenGL fog assumes that the fog has a fixed density \(D\), and then the fractional transparency is given by a formula such as
\[e^{-D\; z}\]
where \(z\) is the distance from the camera to the fragment being rendered. Halfspace fog is fog that only exists on one side of a specified plane. Most often, the fog would be above or below a certain height.
The discussion of halfspace fog by Eric Lengyel (see references) uses a constant density or a linear density function
\[\delta(x) = a\; x\]
where \(x\) is the distance from the boundary plane and \(a\) is a constant.
Let's recall some of Lengyel's notation. Let \(\mathbf{F}\) be a 4-dimensional vector describing a halfspace, such that a homogeneous point \(\mathbf{X}\) is on the foggy side of the plane when \(\mathbf{F}•\mathbf{X} < 0\). We will further assume that the \(xyz\) part of \(\mathbf{F}\) is normalized, so that \(\mathbf{F}•\mathbf{X}\) is the signed distance between \(\mathbf{X}\) and the boundary plane. Let \(\mathbf{C}\) be the camera position, let \(\mathbf{P}\) be the location of a fragment being rendered, and let \(\mathbf{V} = \mathbf{C} - \mathbf{P}\). We can write a parametric equation of the line segment from \(\mathbf{P}\) to \(\mathbf{C}\) as \(\mathbf{Q}(t) = \mathbf{P} + t\,\mathbf{V}\) where \(t\) varies from 0 to 1. The differential distance can be expressed as \(ds = \left\|\mathbf{V}\right\| dt\). The function \(g(\mathbf{P})\) will represent the amount of fog encountered between \(\mathbf{P}\) and \(\mathbf{C}\), and can replace the product \(D\; z\) in the exponential fog formula.
Let's see how far we can get without getting specific about the density function \(\delta\). We will just assume that it has an antiderivative \(\Delta\).
First consider the case where \(\mathbf{P}\) and \(\mathbf{C}\) are both on the foggy side of the plane, i.e., \(\mathbf{F}•\mathbf{P} < 0\) and \(\mathbf{F}•\mathbf{C} < 0\). In that case,
\[g(\mathbf{P}) = \int_{P}^{C}dg = \int_{0}^{1}\delta(-\mathbf{F}•\mathbf{Q}(t))\, \left\|\mathbf{V}\right\|\, dt = \left\|\mathbf{V}\right\|\int_{0}^{1}\delta(-\mathbf{F}•\mathbf{Q}(t))\,dt\;.\]
Note that \(\mathbf{F}•\mathbf{Q}(t) = \mathbf{F}•\mathbf{P} + t\,\mathbf{F}•\mathbf{V}\). For convenience, let \(p = \mathbf{F}•\mathbf{P}\) and \(v = \mathbf{F}•\mathbf{V}\), so \(\mathbf{F}•\mathbf{Q}(t) = p + t\,v\). Thus, we need to evaluate
\[g(\mathbf{P}) = \left\|\mathbf{V}\right\|\int_{0}^{1}\delta(-(p + t\,v))\,dt\;.\]
Let's make the substitution \(u = p + t\,v\), so this becomes
\[g(\mathbf{P}) = \left\|\mathbf{V}\right\|\int_{p}^{p+v}\delta(-u)\,\frac{1}{v}\,du = \frac{\left\|\mathbf{V}\right\|}{v}\int_{p}^{p+v}\delta(-u)\,du\;.\]
Let \(c = p+v = \mathbf{F}•\mathbf{C}\), which produces
\[g(\mathbf{P}) = \left\|\mathbf{V}\right\|\int_{p}^{c}\delta(-u)\,\frac{1}{c-p}\,du = \frac{\left\|\mathbf{V}\right\|}{c-p}\int_{p}^{c}\delta(-u)\,du\;.\]
Substitute \(w\) for \(-u\), and we can rewrite this as
\[ g(\mathbf{P}) = \frac{-\left\|\mathbf{V}\right\|}{c-p}\int_{-p}^{-c}\delta(w)\,dw = \frac{-\left\|\mathbf{V}\right\|}{c-p}\,(\Delta(-c) - \Delta(-p)) \;. \]
It is possible that \(p=c\), in which case this equation involves a division by zero. However, that's actually not much of a problem, since it has a finite limit. In fact, if you define \(f(x) = \Delta(-x)\), we can say:
\[ \lim_{p\to c}\frac{\Delta(-c) - \Delta(-p)}{c-p} = \lim_{p\to c}\frac{f(c) - f(p)}{c-p} = f^\prime(c) = -\delta(-x) \]
We won't be far off if we just add a small epsilon to the denominator. In terms of GLSL built-in functions, we could replace \(c-p\) by \((\mathrm{abs}(c-p)+\epsilon)\,\mathrm{sign}(c-p)\).
If \(\mathbf{P}\) and \(\mathbf{C}\) are both on the non-foggy side of the plane, then of course the line segment encounters no fog and \(g(\mathbf{P}) = 0\). It remains to consider the cases where \(\mathbf{P}\) and \(\mathbf{C}\) are on different sides of the plane. In either case, the line segment encounters the plane where \(p + t\,v = 0\), so \(t = -p/v\). One can verify that \(-p/v\) is between 0 and 1. This parameter will take the place of one of the limits of integration.
If \(\mathbf{F}•\mathbf{P} < 0\) and \(\mathbf{F}•\mathbf{C} > 0\), the line starts in the fog, so \(t\) ranges from 0 to \(-p/v\). As a result, the upper limit on \(u\) and \(w\) becomes 0. The integral reduces to
\[ g(\mathbf{P}) = \frac{-\left\|\mathbf{V}\right\|}{c-p}\int_{-p}^{-c}\delta(w)\,dw = \frac{-\left\|\mathbf{V}\right\|}{c-p}\,(\Delta(0) - \Delta(-p)) \;. \]
Similarly, if \(\mathbf{F}•\mathbf{P} > 0\) and \(\mathbf{F}•\mathbf{C} < 0\), the lower limit on \(u\) and \(w\) becomes 0, and the integral becomes
\[ g(\mathbf{P}) = \frac{-\left\|\mathbf{V}\right\|}{c-p}\int_{-p}^{-c}\delta(w)\,dw = \frac{-\left\|\mathbf{V}\right\|}{c-p}\,(\Delta(-c) - \Delta(0)) \;. \]
In these cases, there can be no division by zero, since \(c\) and \(p\) are of opposite sign.
In the OpenGL shading language, we want to reduce branching, but there is a built-in minimum function. Let's define \(\bar{c} = \min(c, 0)\) and \(\bar{p} = \min(p, 0)\). Then we can handle all four cases with one formula:
\[ g(\mathbf{P}) = \frac{-\left\|\mathbf{V}\right\|}{c-p}\int_{-p}^{-c}\delta(w)\,dw = \frac{-\left\|\mathbf{V}\right\|}{c-p}\,(\Delta(-\bar{c}) - \Delta(-\bar{p})) \;. \]
To make the model more realistic than constant or linear halfspace fog, I want a density function that starts out with zero density at the boundary plane and smoothly approaches a fixed density \(D\) as you get farther from the boundary plane, with an adjustable derivative \(S\) (for boundary sharpness) at zero. I also want a function with an antiderivative that can be easily computed in a shader program.
One can design many such functions. For example, if you let
\[\delta(x) = D\left(1 - \left(\frac{S}{2\,D}\,x + 1\right)^{-2}\right)\]
which has the antiderivative
\[\Delta(x) = D\left(x + \frac{2\,D}{S}\left(\frac{S}{2\,D}\,x + 1\right)^{-1}\right)\]
and plug that into our unified formula and simplify, you obtain:
\[ g(\mathbf{P}) = \frac{\left\|\mathbf{V}\right\|\,D}{c-p} \left(\bar{c} - \bar{p} - \frac{2\,D}{S} \left(\left(1 - \frac{S}{2\,D}\,\bar{c}\right)^{-1} - \left(1 - \frac{S}{2\,D}\,\bar{p}\right)^{-1}\right)\right)\;. \]
Or if you like exponential functions, you can start with
\[\delta(x) = D\,(1-{\large e^{-\frac{S}{D}x}})\]
with the antiderivative
\[\Delta(x) = D\left(x + \frac{D}{S}\,{\large e^{-\frac{S}{D}\,x}}\right)\;.\]
Plugging this into our unified formula and simplifying a bit, we obtain
\[ \large g(\mathbf{P}) = \frac{\left\|\mathbf{V}\right\|\,D}{c-p} \left(\bar{c} - \bar{p} - \frac{D}{S}\left( e^{\frac{S}{D}\,\bar{c}} - e^{\frac{S}{D}\,\bar{p}}\right)\right)\;. \]
Lengyel, Eric, "Unified Distance Formulas for Halfspace Fog", Journal of Graphics Tools, vol. 12, number 2, pp 23-32, 2007.