# Numerically integratable fog volume feathering along a ray

I haven't been able to come up with a good solution for this myself, so I thought I'd ask if anyone else has any good ideas on how to approach this problem. Here's more or less what I'm trying to do and what I've been able to research / think of so far:

I have implemented volumetric light rendering in my engine, and the way you feed atmospheric effects right now is either by a constant global fog (parametrized as optical length), a global height fog (pretty much the same as https://iquilezles.org/articles/fog/) or with a local volume that's parametrized using an AABB and density (which modules the optical length) and also the ability to do falloff in the same way as the height fog (so exponential, towards a user direction). One issue with the volume is that the cutoff at the edges of the volume is quite abrupt, so I'd like to add some feathering to it. An issue with this is that it needs to be integrateble along a ray segment in the same way as the height fog. It turns out that this problem is quite tricky.

Here's the code from https://iquilezles.org/articles/fog/ which computes the global height fog, this is pretty much the same as in my engine so let's use this instead of my re-written version of it. The fogAmount is essentially computing the optical length from rayOri towards rayDir for a known distance (e.g. surface hit). Note that all it needs are these paramters, and it can analytically figure out the optical length from just that (it looks at height difference change to be able to do it analytically).

vec3 applyFog( in vec3  rgb,      // original color of the pixel
in float distance, // camera to point distance
in vec3  rayOri,   // camera position
in vec3  rayDir )  // camera to point vector
{
float fogAmount = (a/b) * exp(-rayOri.y*b) * (1.0-exp( -distance*rayDir.y*b ))/rayDir.y;
vec3  fogColor  = vec3(0.5,0.6,0.7);
return mix( rgb, fogColor, fogAmount );
}


The problem I'm trying to solve is similar, but I want to be able to, analytically, figure out the optical density if you had a AABB with some amount of feathering at the edge. The function signature is the same, but instead of a and b as parameters (these define the height falloff) we'd just have a feather parameter that represents how far the feathering should be applied inside the volume. Here's a figure of how something like this looks like:

I've not been able to find a good representation to modulate the optical length to make it be numerically integrateable in the same way as Inigo Quilez's height fog. It's possible to find the distance from some point to the surface of an AABB like this:

float distance = length( max( max( aabb.min - point, 0.0f ), max( point - aabb.max, 0.0f ) );


And it would be possible to use that as input to e.g. smoothstep to modulate the optical length: smoothstep( 0.0f, feather, distance ) but there is no physical point to do this with, as this is a ray integral. So I guess my "real" question is: is this even possible to solve if this has the form of a ray integral? For height fog à la Inigo Quilez it's possible (it's mentioned in the article that it's possible to derive new such functions by solving the ray integral again with a different set of parameters) but is it also possible to do it in this case as well where we would want to have feathering on the edges of an AABB to change the optical length?

I realize that this is a bit of an odd question that's more exploratory than anything else... and hopefully I was able to explain what the problem was. Any ideas and help are much appreciated! I've run out of good ideas for doing it as a ray integral :-(

I have a gut feeling that this somehow involves computing rate of change based on the feathering and taking the incoming ray direction to modulate it (much like the height fog is doing)...

Thanks for the help!

The value of an integration is the area under the curve between the limits of integration. Since the density of the fog is constantly increasing starting at zero until it reaches the density of the constant fog then the value of the increasing fog is directly related to the value of the constant fog by 1/2 since the function we are integrating for both is the same (or directly related only changing by density) at any give (infinitesimally small) slice. Also angle doesn't matter since the rate of change of the fog is still constant only being changed by the angle.

But, this is only true where the view goes completely from zero to constant fog density. In regions where a fragment is partially inside the non-constant fog then the two are related by a ratio of the distance between the pixel and the edge of the constant region. For example if a pixel is half way through the non-constant fog then it's value is 1/2 of 1/2 or 1/4.

Unfortunately there is another case where the camera is outside the fog and there is a pixel whose view vector runs parallel to the constant fog and passes only through the non-constant fog. In this case the non-constant fog is essentially constant but at a lowered density. This region has a changing section that increases until the lowered density and then has constant density until it saturates or leaves the lowered density region through another changing density region. In this case the new constant density is related to the ratio of the distance from the constant density fog again. All other calculations remain the same.

Unfortunately there is another corner case where the view ray passes through the non-constant density region at a corner of the bounding region and the fragment is in a zero density region. In this case the fog calculation is equal to the same value as constant density region where the constant density is related to the distance of the constant bounding region from the view ray. (line to bounding box distance). It is equal to the constant calculation because it pass through an increasing section and a decreasing section of equal size. So it is 2 times 1/2 of constant.

There is yet another set of corner cases where the view starts inside the non-constant region and looks into either the corners, the constant region or runs parallel to the constant region.

Implementation boils down to doing line to box intersections where misses return the distance of the ray from the box and a time T along the ray where the ray is closest to the box or intersects with it. Then using the distance and T to classify the fragment into various corner cases. Then finally computing the fog as a sum of the various cases.

At the end of the day there are so many corner cases (no pun intended) that switching to a ray marched implementation becomes the more attractive solution. (at least using the method I described here)