Why does this gl_FragDepth calculation work?

Question

I am raytracing in GLSL. My gl_FragDepth calculation wasn't right, so I did some hunting and found this solution (P is the world-space coordinate of the pixel and pv is proj*view):

    vec4 Pclip = pv * vec4 (P, 1);

    float ndc_depth = Pclip.z / Pclip.w;

    gl_FragDepth = (
        (gl_DepthRange.diff * ndc_depth) +
        gl_DepthRange.near + gl_DepthRange.far) / 2.0;

It works. But why?

Okay, so I need to divide by w and scale it to the expected z-range. That makes sense qualitatively, but the calculation surprised me.

If you told me to "divide by w and scale to gl_DepthRange" I would have written

float ndc_depth = Pclip.z / Pclip.w;

gl_FragDepth = gl_DepthRange.diff * ndc_depth + glDepthRange.near;

At this point I'd like to point out that this depth calculation looks correct with the raytraced geometry which is intersecting with my normal triangle models, so I guess OpenGL is performing the same depth-scaling implicitly for those triangles as well.

Three questions:

1. Why add gl_DepthRange.far to the calculation?
2. Why divide by 2?
3. I set the pv matrix in C++ with glm::perspective(...)*glm::lookAt(...) so how the hell does gl_DepthRange know what the near and far planes are?

Answer 1

What you are missing is, that in OpenGL's NDC space (i.e. clip space after division by w) all 3 coordinates are in the range $[-1,1]$. So ndc_depth is in the range $[-1,1]$, while your computation assumes a range of $[0,1]$.

So let's take your computation but add an initial step to map $[-1,1]$ to $[0,1]$:

float ndc_depth = Pclip.z / Pclip.w;
float win_depth = (ndc_depth + 1.0) / 2.0;
gl_FragDepth = gl_DepthRange.diff * win_depth + gl_DepthRange.near;

If we do some transformations (and realizing that gl_DepthRange.diff is just convenience for gl_DepthRange.far - gl_DepthRange.near), we get

gl_FragDepth = gl_DepthRange.diff * (ndc_depth + 1.0) / 2.0 + gl_DepthRange.near;

gl_FragDepth = (gl_DepthRange.diff * ndc_depth + gl_DepthRange.diff) / 2.0 
               + gl_DepthRange.near;

gl_FragDepth = (gl_DepthRange.diff * ndc_depth + gl_DepthRange.diff + 
                2.0 * gl_DepthRange.near) / 2.0;

gl_FragDepth = (gl_DepthRange.diff * ndc_depth + 
                gl_DepthRange.far - gl_DepthRange.near + 2.0 * gl_DepthRange.near) / 2.0;

gl_FragDepth = (gl_DepthRange.diff * ndc_depth + 
                gl_DepthRange.far + gl_DepthRange.near) / 2.0;

So, to answer your specfic questions:

1. That's just a result from simplifying the calculation based on the fact that gl_DepthRange.diff is nothing else than gl_DepthRange.far - gl_DepthRange.near.

2. Because the depth range in NDC space is $[-1,1]$, which has a width of $2$.

3. It doesn't. These values have nothing to do with the near and far values you used to construct your projection matrix, in fact you could use any arbitrary projection matrix (or none at all). They are the depth range values set with glDepthRange, which have a totally different meaning (okay, let's say slightly different but related).

   Basically, and special situations notwithstanding, the projection matrix defines where you near and far clipping planes are in eye-space, i.e. which eye-space range is mapped to the $[-1,1]$ NDC depth range. Then the depth range defines what these NDC near/far values of -1 and 1 are mapped to in the actual depth buffer. So the names "near"/"far" do relate to the near and far clipping planes. But both transformations are basically independent of each other, using the general NDC (i.e. the normalized device space) as a mediator between both stages. This way glDepthRange more or less just extends what glViewport does onto the z-dimension.