/* Conversion of z-buffer values to "real" depth values
** by Wolfgang Stuerzlinger (stuerzlingerw AT-SIGN acm DOT org) 1997-
** Thanks to Gernot Schaufler for providing some essential hints/code segments
*/

#define SIZE 256                            // width & height of framebuffer

static float pcorrect[SIZE * SIZE];            // for perspective correction

// precomputation of perspective correction; has to be done only once
// because z-values are stored as computed in "orthogonal box-space"
// a plane normal to the viewer will have constant z-values for all points.
// for a perspective view this is not true.
// the computation below is for a fov of 90 degrees!

    int i,j;
    float step;
    Vector3f vec;

    step = 2. / SIZE;
    for (j = 0; j < SIZE;j++) {
        for (i = 0; i < SIZE; i++) {
            vec.Set(-1. + i * step + step * .5,            
                    -1. + j * step + step * .5, 1.);
            pcorrect[j * SIZE + i] = vec.Length();
            }
        }

// computation of depth values
// the explanation of the "magic" code: the code is the "inversion" of the
// depth.c source code provided by SGI in the example programs provided for
// the OpenGL SIGGRAPH 96 course.
// http://www.sgi.com/Technology/OpenGL/advanced/programs.html

    float cfar,cnear;                // input: near & far plane distance
    float fn,ztemp;
    float z[SIZE * SIZE];            // output: real depth values
    GL_FLOAT temp[SIZE * SIZE];

    // read the pixels from the frame buffer (SIZE * SIZE pixels)
    glReadPixels(rasterpos,0,SIZE,SIZE,GL_DEPTH_COMPONENT,GL_FLOAT,temp);

    // cnear & cfar are the near/far plane of the camera
    fn = cfar - cnear;

    // loop over all pixels
    for (i = 0; i < SIZE * SIZE; i++) {
        // and magic happens...
#ifdef ORIGINAL
        // direct "inversion" of depth.c
        ztemp = (2. * temp[i] - 1.) - (cfar + cnear) / fn;
        ztemp = -2. * cnear * cfar / ztemp / fn;
#else
        // algebraically simplified version
        ztemp = cfar * cnear / (temp[i] * (cnear - cfar) + cfar);
#endif

        // the z-value is now in "orthogonal space" -> convert to perspective
        z[i] = ztemp * pcorrect[i];
        }