# Nonlinear-density Halfspace Fog

by James W. Walker
25 May 2016

## Review

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.

## General Density

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)$$.

## Crossing the Plane

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.

## Unified Formula

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})) \;.$

## Nonlinear Density

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)\;.$

$\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)\;.$

## References

Lengyel, Eric, "Unified Distance Formulas for Halfspace Fog", Journal of Graphics Tools, vol. 12, number 2, pp 23-32, 2007.