/*======================================================================
 * CGLS: OpenGL Shadows: TEST PROGRAM
 *----------------------------------------------------------------------
 * This program is copyright Michael McCool, 1998.   However, this 
 * program, in whole or part, may be used for any purpose, commercial or 
 * educational, as long as this copyright notice is retained and a
 * note describing the usage is sent to the following:
 *
 * Contact: Michael McCool
 *          mmccool@cgl.uwaterloo.ca
 *          http://www.cgl.uwaterloo.ca/~mmccool/
 *
 * Usage information is printed when the program is run.
 *
 * Note: The program may crash if a light view window is opened and
 *       resized and then deleted.   So if you open one, leave it alone.
 *       Some silly memory-reallocation thing I haven't had time to track 
 *       down...
 *
 * Computer Graphics Lab
 * Department of Computer Science
 * University of Waterloo
 * Waterloo, Ontario
 * N2L 3G1 Canada
 *======================================================================*/
#include <math.h>
#include <stdio.h>
#include <GL/glut.h>

/* use RLE encoding of edge maps.   The default; this define is only
 * here so I can compile it without and time the difference
 */
#define RLE 1

/* Forward references */
void update_menus();

void
print_sep (
    FILE* fp,
    int c
) {
    int i;
    for (i=0; i<79; i++) {
        fputc(c,fp);
    }
    fputc('\n',fp);
}

/************************************************************************
 * Error Checking
 *
 * Check if an OpenGL error has occured, and if so, print it in English
 * and die.
 ************************************************************************/

void
check_error(char *w) {
    int ecode, err;

    err = 0;
    while ((ecode = glGetError()) != GL_NO_ERROR) {
        fprintf(stderr,"%s: OpenGL ERROR [%d] %s\n",
            w,ecode,gluErrorString(ecode));
        err = 1;
    }
    if (err) exit(1);
}


/************************************************************************
 * State
 *
 * Here are most of the state variables for the interface and the 
 * values they can take.
 ************************************************************************/

/* Window identifiers
 */
int main_window;
int light_window;

/* Light window status
 */
int light_window_active;

/* Resolution of main window
 */
int xr = 400;
int yr = 400;

/* Resolution of light window
 */
int lxr = 400;
int lyr = 400;

/* Allocated size of shadow map
 */
int map_xr = -1;
int map_yr = -1;

/* Depth of shadow map 
 */
typedef GLint map_z_type;
#define MAP_Z_TYPE GL_INT
#define MAP_Z_FORMAT "%d"
#define MAX_MAP_Z 0x7fffffffL
#define MIN_MAP_Z 0x80000000L

map_z_type** map_z = NULL;

/* Edge detection thresholds
 * (set automatically after detecting depth precision).
 */
map_z_type th_H;
map_z_type th_L;
map_z_type th_s;

/* Macro to convert a depth map sample into an OpenGL vertex.
 *
 * We assume maxxed out z values are the background at infinity.
 * HOWEVER, if real points at infinity are used, when they are
 * mapped through the projective transformation to the light
 * source, they can wrap around behind the light source.   Really
 * what we should do is analyse the inverse projective transformation
 * and pick "safe" large z values, but in the meantime we just
 * use fudge factors.
 */
#define WZ 1e-2
#define ZINF 1e8

#define VERTEX(i,j)                                                 \
    if (map_z[j][i] == MAX_MAP_Z) {                                 \
        glVertex4d(WZ*(i), WZ*(j), ZINF, WZ);                       \
    } else {                                                        \
        glVertex3i((GLint)i,                                        \
                   (GLint)j,                                        \
                   (GLint)map_z[j][i]);                             \
    }

/* Should optimize this a bit more someday...
 */
#define TRI(i0,j0, i1,j1, i2,j2)                                    \
    VERTEX(i0,j0)                                                   \
    VERTEX(i1,j1)                                                   \
    VERTEX(i2,j2)                                                   

#define glTRI(i0,j0, i1,j1, i2,j2)                                  \
    glBegin(GL_TRIANGLES);                                          \
        TRI(i0,j0, i1,j1, i2,j2)                                    \
    glEnd();

/* Macro to evaluate a sample against a plane equation.
 *
 * We have to take the hack above into account...
 * p is the plane equation, the result is returned in q.
 */
#define PLANE(p,i,j,q)                                              \
    {                                                               \
        GLdouble x, y, z, w;                                        \
        if (map_z[j][i] == MAX_MAP_Z) {                             \
            x = WZ*(i); y = WZ*(j); z = ZINF; w = WZ;               \
        } else {                                                    \
            x = i; y = j; z = map_z[j][i]; w = 1;                   \
        }                                                           \
        q = p[0]*x + p[1]*y + p[2]*z + p[3]*w;                      \
    }

/* Macro to generate a cap vertex.
 *
 * We intersect a ray passing through the sample with the cap plane.
 * Note that we DON'T do a depth test; we assume it's been done
 * already.
 */
#define CAP_VERTEX(p,i,j)                                                   \
    if (p[2] > 0) {                                                         \
        glVertex4d(p[2]*(i),p[2]*(j),-(p[0]*(i) + p[1]*(j) + p[3]),p[2]);   \
    } else {                                                                \
        glVertex4d(-p[2]*(i),-p[2]*(j),(p[0]*(i) + p[1]*(j) + p[3]),-p[2]); \
    }

/* Macro to generate a clip vertex.
 *
 * These are vertices generated on the cap between samples, due to
 * clipping with the cap plane.   The macro takes the indices of the 
 * two samples and the plane equation generates the cap vertex 
 * under OpenGL.
 *
 * The two sample points should be on opposite sides of the clipping
 * plane, or strange things are likely to happen.
 */
#define CLIP_VERTEX(p,i0,j0,i1,j1)                                      \
    {                                                                   \
        GLdouble x0, x1, y0, y1, z0, z1, w0, w1, p0, p1, t0, t1;        \
        if (map_z[j0][i0] == MAX_MAP_Z) {                               \
            x0 = WZ*(i0); y0 = WZ*(j0); z0 = ZINF; w0 = WZ;             \
        } else {                                                        \
            x0 = i0; y0 = j0; z0 = map_z[j0][i0]; w0 = 1.0;             \
        }                                                               \
        if (map_z[j1][i1] == MAX_MAP_Z) {                               \
            x1 = WZ*(i1); y1 = WZ*(j1); z1 = ZINF; w1 = WZ;             \
        } else {                                                        \
            x1 = i1; y1 = j1; z1 = map_z[j1][i1]; w1 = 1.0;             \
        }                                                               \
        p0 = p[0]*x0 + p[1]*y0 + p[2]*z0 + p[3]*w0;                     \
        p1 = p[0]*x1 + p[1]*y1 + p[2]*z1 + p[3]*w1;                     \
        t0 = p1 / (p1 - p0); t1 = 1.0 - t0;                             \
        glVertex4d(t0*x0+t1*x1,t0*y0+t1*y1,t0*z0+t1*z1,t0*w0+t1*w1);    \
    }

/* Edge and cap codes
 * 
 * EDGE CODES:
 *
 * One of 16 codes defining various configurations of 
 * edges, organized as a combination of 4 bits.  
 *
 * CAP CODES:
 *
 * One of 16 codes defining various configurations of 
 * samples above/below the camera near clipping plane,
 * organized as a combination of 4 bits.
 * 
 * INDEX_CODES:
 *
 * For run-length encoding of edge maps (this increases
 * reconstruction speed by a significant factor, i.e. up to
 * 10x)
 *
 * Here's a map (the corners are depth map samples):
 *
 *    01 H1 11
 *         
 *    V0    V1     
 *         
 *    00 H0 10
 */
#ifdef RLE

typedef GLuint map_code_type;

#define MAP_CODE_TYPE GL_UNSIGNED_INT
#define MAP_CODE_FORMAT "0x%x"

#else

typedef GLubyte map_code_type;

#define MAP_CODE_TYPE GL_UNSIGNED_BYTE
#define MAP_CODE_FORMAT "0x%x"

#endif

#define EDGE_H0 (1L << 2)
#define EDGE_H1 (1L << 3)
#define EDGE_V0 (1L << 0)
#define EDGE_V1 (1L << 1)

#define EDGE_MASK 0x0000000fL

#define EDGE(v0,v1,h0,h1)                       \
    (                                           \
        ((v0) ? EDGE_V0 : 0L) |                 \
        ((v1) ? EDGE_V1 : 0L) |                 \
        ((h0) ? EDGE_H0 : 0L) |                 \
        ((h1) ? EDGE_H1 : 0L)                   \
    )

#define PEDGE(fp,ec)                            \
    fprintf(fp,"/%c%c%c%c/",                    \
        (((ec) & EDGE_V0) ? '<' : '.'),         \
        (((ec) & EDGE_V1) ? '>' : '.'),         \
        (((ec) & EDGE_H0) ? 'V' : '.'),         \
        (((ec) & EDGE_H1) ? 'A' : '.')          \
    )

/* Capping codes packed into same storage as edge codes
 */
#define CAP_00  (1L << 4)
#define CAP_01  (1L << 5)
#define CAP_10  (1L << 6)
#define CAP_11  (1L << 7)

#define CAP_MASK 0x000000f0L

#define CAP(s00,s01,s10,s11)               \
    (                                           \
        ((s00) ? CAP_00 : 0) |             \
        ((s01) ? CAP_01 : 0) |             \
        ((s10) ? CAP_10 : 0) |             \
        ((s11) ? CAP_11 : 0)               \
    )

/* Reconstruction category codes
 *
 * 3-bit codes describing the triangles that were generated in
 * the reconstruction phase.
 *
 * U, D, and V are mutually exclusive.  They should be read as
 * "up", "down", and "vector".
 *
 * CAT_U_01:
 *
 *    *   * 
 *       /| 
 *      / | 
 *     /  |
 *    *---* 
 *
 * CAT_U_10:
 *
 *    *---* 
 *    |  /  
 *    | /   
 *    |/   
 *    *   * 
 *
 * CAT_U_10:
 *
 *    *---* 
 *    |  /| 
 *    | / | 
 *    |/  |
 *    *---* 
 *
 * CAT_D_11:
 *
 *    *   * 
 *    |\      
 *    | \   
 *    |  \ 
 *    *---* 
 *
 * CAT_D_10:
 *
 *    *---* 
 *     \  | 
 *      \ |  
 *       \|
 *    *   * 
 *
 * CAT_D_11:
 *
 *    *---* 
 *    |\  | 
 *    | \ |  
 *    |  \|
 *    *---* 
 */
#define CAT_NONE ((0L << 10) | (0L << 9) | (0L << 8))
#define CAT_U_01 ((0L << 10) | (1L << 9) | (0L << 8))
#define CAT_U_10 ((1L << 10) | (0L << 9) | (0L << 8))
#define CAT_U_11 ((1L << 10) | (1L << 9) | (0L << 8))
#define CAT_D_01 ((0L << 10) | (1L << 9) | (1L << 8))
#define CAT_D_10 ((1L << 10) | (0L << 9) | (1L << 8))
#define CAT_D_11 ((1L << 10) | (1L << 9) | (1L << 8))
#define CAT_VEC  ((0L << 10) | (0L << 9) | (1L << 8))

#define CAT_MASK 0x00000700

#define INDEX_MASK 0xfffff800
#define INDEX_SHIFT 11

map_code_type** map_code = NULL;

/* Flag to indicate that edges and shadow volume need to be
 * reacquired
 */
int edge_update = 1;

/* Epsilons to hide shadow map grunginess.
 */
GLdouble light_ze;
GLdouble eye_ze;
GLdouble line_ze;

/* Number of subdivisions in models
 */
#define SUBD_LOW     10
#define SUBD_MEDIUM  30
#define SUBD_HIGH    60
int subd = SUBD_MEDIUM;

/* View control mode (main window)
 */
#define VIEW_RESET_SCENE            0
#define VIEW_TUMBLE_SCENE           1
#define VIEW_TUMBLE_LIGHT           2
#define VIEW_TUMBLE_SCENE_AND_LIGHT (VIEW_TUMBLE_SCENE|VIEW_TUMBLE_LIGHT)
int view_mode = VIEW_TUMBLE_SCENE;

/* Scene selection
 */
#define NSCENES 11

#define SCENE_START_DEMO        0
#define SCENE_END_DEMO          5

#define SCENE_N                 8

#define SCENE_SPHERE_AND_CUBE_1 0
#define SCENE_SPHERE_AND_CUBE_2 1
#define SCENE_CUBES             2
#define SCENE_GRID              3
#define SCENE_SPHERES           4
#define SCENE_CONE_AND_TORUS    5
#define SCENE_TEAPOT            6
#define SCENE_TEAPOTS           7
char* scene_name[] = {
    "Sphere and Cube 1",
    "Sphere and Cube 2",
    "Cubes",
    "Grid",
    "Spheres",
    "Cone and Torus",
    "Teapot",
    "Teapots"
};

#define SCENE_SHADOW_VOLUME     10

static int main_scene_base = -1;
static int light_scene_base = -1;
static int scene = SCENE_SPHERE_AND_CUBE_1;

/* Following matrix code shamelessly stolen from Mesa and modified
 * for our purposes
 */

/* 4x4 Matrix Type
 */
typedef GLdouble Mat4x4[16];

/* 2x2 Matrix Type
 */
typedef GLdouble Mat2x2[4];

/* Macro to access elements of 4x4 matrix
 * Note that OpenGL matrices are column major, not row major as
 * usual in C.
 */
#define M4(i,j) ((i) + 4*(j))

/* Macro to access elements of 2x2 matrix
 */
#define M2(i,j) ((i) + 2*(j))

/* Identity matrix
 */
Mat4x4 Identity = {
   1.0, 0.0, 0.0, 0.0,
   0.0, 1.0, 0.0, 0.0,
   0.0, 0.0, 1.0, 0.0,
   0.0, 0.0, 0.0, 1.0
};

Mat4x4 eye_xform = {
    1.0, 0.0, 0.0, 0.0,
    0.0, 1.0, 0.0, 0.0,
    0.0, 0.0, 1.0, 0.0,
    0.0, 0.0, 0.0, 1.0
};

Mat4x4 light_xform = {
    1.0, 0.0, 0.0, 0.0,
    0.0, 1.0, 0.0, 0.0,
    0.0, 0.0, 1.0, 0.0,
    0.0, 0.0, 0.0, 1.0
};

/* Demo mode control
 */
int demo_mode = 0;
int old_demo_mode = 0;
int demo_delay = 1;

/* Light state
 */
GLfloat light_position[] = { 0.0, 0.0, 30.0, 1.0 };
GLfloat light_power[] = { 1.0, 1.0, 1.0, 1.0 };
GLfloat light_fov = 30.0;
GLfloat light_near = 5.0;
GLfloat light_far = 45.0;

/* Light window mode
 */
#define LW_SHADED         0
#define LW_DEPTH          1
#define LW_EDGE           2
#define LW_HW_EDGE        3
int light_window_mode = LW_SHADED;

/* Eye state
 */
GLfloat eye_position[] = { 0.0, 0.0, 15.0, 1.0 };
GLfloat eye_lookat[] = { 0.0, 0.0, 0.0, 1.0 };
GLfloat eye_fov = 45.0;
GLfloat eye_near = 2.0;
GLfloat eye_far = 30.0;

/* Ambient lighting state
 */
int ambient_mode = 0;
GLfloat ambient[] = { 0.0, 0.0, 0.0, 1.0 };

/* Display mode (main window)
 */
#define DISPLAY_N                           12

#define DISPLAY_SCENE                       0
#define DISPLAY_DARK_SHADOWS                1
#define DISPLAY_AMBIENT_SHADOWS             2
#define DISPLAY_COMPOSITED_SHADOWS          3
#define DISPLAY_SHADOW_VOLUME               4
#define DISPLAY_OUTLINED_SHADOW_VOLUME      5
#define DISPLAY_SHADOW_VOLUME_POINTS        6
#define DISPLAY_SHADOW_VOLUME_MESH          7
#define DISPLAY_SHADOW_VOLUME_NS            8
#define DISPLAY_OUTLINED_SHADOW_VOLUME_NS   9
#define DISPLAY_SHADOW_VOLUME_POINTS_NS     10
#define DISPLAY_SHADOW_VOLUME_MESH_NS       11
char* display_mode_name[] = {
    "No Shadowing",
    "Dark Shadows",
    "Ambient Shadows",
    "Composited Shadows",
    "Shadow Volume",
    "Outlined Shadow Volume",
    "Shadow Volume Points",
    "Shadow Volume Mesh",
    "Shadow Volume, No Scene",
    "Outlined Shadow Volume, No Scene",
    "Shadow Volume Points, No Scene",
    "Shadow Volume Mesh, No Scene"
};

int display_mode = DISPLAY_SCENE;

/* Edge detector mode
 */
#define EDGE_N             2

#define EDGE_SIMPLE        0
#define EDGE_CANNY         1
char* edge_mode_name[] = {
    "Gradient Threshold",
    "Canny"
};

int edge_mode = EDGE_CANNY;

/* Shadow volume reconstruction mode
 */
#define RECON_N           6

#define RECON_ALL         0
#define RECON_QUADS       1
#define RECON_QUADS_FLIP  2
#define RECON_TRIS        3
#define RECON_TRIS_FLIP   4
#define RECON_VECTORIZE   5
char* recon_mode_name[] = {
    "All (edge detection disabled)",
    "Quads",
    "Flipped Quads",
    "Triangles",
    "Flipped Triangles",
    "Vectorize",
};

int recon_mode = RECON_TRIS_FLIP;

/* Shadow volume capping mode
 */
#define CAP_N             2

#define CAP_OFF           0
#define CAP_VIEW_NEAR     1
char* cap_mode_name[] = {
    "Off",
    "View Near"
};

int cap_mode = CAP_VIEW_NEAR;

/* Lighting mode
 */
#define LIGHT_POINT        0
#define LIGHT_DIRECTION    1
#define LIGHT_SPOT         2
int light_mode = LIGHT_POINT;

/* Main menu indices
 */ 
#define M_SCENE          1
#define M_DISPLAY        2
#define M_EDGE           3
#define M_RECON          4
#define M_CAP            5
#define M_AMBIENT        6
#define M_SUBD           7
#define M_VIEW           8
#define M_LIGHT          9
#define M_DEMO           10
#define M_LIGHT_WINDOW   11
#define M_QUIT           12

#define LM_CLOSE         1
#define LM_QUIT          2

/* Menu identifiers
 */
int main_menu_id;

int scene_menu_id;
int display_menu_id;
int edge_menu_id;
int recon_menu_id;
int cap_menu_id;
int ambient_menu_id;
int subd_menu_id;
int view_menu_id;
int light_menu_id;
int demo_menu_id;

int light_window_menu_id;
int light_window_mode_menu_id;

/* Flags to manage menu updates
 */
int menu_state = GLUT_MENU_NOT_IN_USE;
int menu_change = 0;

/************************************************************************
 * Timing
 *
 * Some stuff to time durations of various things, to obtain performance
 * numbers.   This stuff is only guaranteed to work on SGIs, hence the
 * ifdef.    Real time AND cpu time is evaluated, but to maximize real
 * time performance the system should not be otherwise loaded.  
 *
 * These functions evaluate both the variance of the measurements taken
 * and the variance that could have been due to quantization.   Only when
 * the number of trials is such that the variance of the result is 
 * on the same order as the quantization variance can we trust the
 * result.
 ************************************************************************/

int timer_report = 0;
int timer_frames = 0;

#ifdef TIMERS
#include <time.h>

typedef struct {
    double  s;
    double  e;
    double  sq;

    int     n;
} Timer;

/* Initialize a timer
 */
void
timer_init (
    Timer* t
) {
    t->s = 0.0;
    t->e = 0.0;
    t->sq = 0.0;
    t->n = 0;
}

/* Start a timer
 */
void
timer_start(
    Timer* t
) {
    struct timespec p;

    clock_gettime(CLOCK_REALTIME,&p);
    t->s = p.tv_sec + 1e-9*p.tv_nsec;
}

/* Stop a timer (and increment number of measurements)
 */
void
timer_stop (
    Timer* t
) {
    struct timespec p;
    double s, d;

    /* So timings for graphics renderings will be stage-accurate
     */
    glFinish();

    clock_gettime(CLOCK_REALTIME,&p);
    s = p.tv_sec + 1e-9*p.tv_nsec;
    d = s - t->s;
    t->e += d;
    t->sq += d*d;
    t->n += 1;
}

/* Return the timer resolution
 */
double
timer_res () {
    struct timespec p;
    clock_getres(CLOCK_REALTIME,&p);
    return p.tv_sec + 1e-9*p.tv_nsec;
}

void
timer_title (
    FILE* fp
) {
    fprintf(fp,"%20s ","PHASE");
    fprintf(fp,"%8s ","N");
    fprintf(fp,"%10s ","Total (s)");
    fprintf(fp,"%10s ","Mean (s)");
    fprintf(fp,"%10s ","Std Err (s)");
    fprintf(fp,"%10s ","Mean Rate (Hz)");
    fprintf(fp,"\n");
}

/* Print stats for a particular timer (stop timer first)
 */
void
timer_print (
    FILE* fp,
    char* msg,
    Timer* t
) {
    double m, v;
    struct timespec p;
    clock_getres(CLOCK_REALTIME,&p);

    fprintf(fp,"%20s ",msg);
    if (t->n > 0) {
        fprintf(fp,"%8d ",t->n);
        fprintf(fp,"%10.5f ",t->e);
        m = t->e/t->n;
        fprintf(fp,"%10.5f ",m);
        if (t->n > 1) {
            v = (t->sq - (t->n)*m*m)/(t->n-1);
            fprintf(fp,"%10.5f ",sqrt(v));
        } else {
            fprintf(fp,"%10s ","UNKNOWN");
        }
        fprintf(fp,"%10.5f ",1.0/m);
    } else {
        fprintf(fp,"NO INFORMATION");
    }
    fprintf(fp,"\n");
}

#else

typedef int Timer;

void 
timer_init (
    Timer* t
) {}

void
timer_start (
    Timer* t
) {}

void
timer_stop (
    Timer* t
) {}

double
timer_res () { }

void
timer_title (
    FILE* fp
) { }

void
timer_print (
    FILE* fp,
    char* msg,
    Timer* t
) {}

#endif

/* Timers for various stages of algorithm
 */
Timer timer_render_scene;
Timer timer_render_depth;
Timer timer_read_depth;
Timer timer_edge_detect;
Timer timer_recon_volume;
Timer timer_render_volume;
Timer timer_cap;
Timer timer_render_shadow;
Timer timer_frame;

void
init_timers () {
    timer_init(&timer_render_scene);
    timer_init(&timer_render_depth);
    timer_init(&timer_read_depth);
    timer_init(&timer_edge_detect);
    timer_init(&timer_recon_volume);
    timer_init(&timer_render_volume);
    timer_init(&timer_cap);
    timer_init(&timer_render_shadow);
    timer_init(&timer_frame);
}

void
print_timers (  
    FILE* fp
) {
    print_sep(fp,'-');
    fprintf(fp,"Timer resolution: %fs\n",timer_res());
    fprintf(fp,"Main window resolution: %d x %d\n",xr,yr);
    if (light_window_active) {
        fprintf(fp,"Light window resolution: %d x %d\n",lxr,lyr);
    }
    fprintf(fp,"Shadow map resolution: %d x %d\n",map_xr,map_yr);
    fprintf(fp,"Display mode: %s\n",display_mode_name[display_mode]);
    fprintf(fp,"Edge mode: %s\n",edge_mode_name[edge_mode]);
    fprintf(fp,"Recon mode: %s\n",recon_mode_name[recon_mode]);
    fprintf(fp,"Cap mode: %s\n",cap_mode_name[cap_mode]);
    fprintf(fp,"Scene: %s\n",scene_name[scene]);
    timer_title(fp);
    timer_print(fp,"Render Scene",&timer_render_scene);
    timer_print(fp,"Render Depth",&timer_render_depth);
    timer_print(fp,"Read Depth",&timer_read_depth);
    timer_print(fp,"Edge Detection",&timer_edge_detect);
    timer_print(fp,"Reconstruction",&timer_recon_volume);
    timer_print(fp,"Stencil Volume",&timer_render_volume);
    timer_print(fp,"Cap",&timer_cap);
    timer_print(fp,"Render Shadows",&timer_render_shadow);
    timer_print(fp,"Total Frame",&timer_frame);
    fflush(fp);
}

/************************************************************************
 * RGBA Colours
 ************************************************************************/

GLfloat black[] =  { 0.0, 0.0, 0.0, 1.0 };
GLfloat grey[] =   { 0.5, 0.5, 0.5, 1.0 };
GLfloat white[] =  { 1.0, 1.0, 1.0, 1.0 };
GLfloat red[] =    { 1.0, 0.0, 0.0, 1.0 };
GLfloat green[] =  { 0.0, 1.0, 0.0, 1.0 };
GLfloat blue[] =   { 0.0, 0.0, 1.0, 1.0 };
GLfloat yellow[] = { 1.0, 1.0, 0.0, 1.0 };

GLfloat transparent_black[] =  { 0.0, 0.0, 0.0, 0.5 };

/************************************************************************
 * Matrix Manipulation (in OpenGL column-major format)
 ************************************************************************/

/* Transform a vector.
 */
void
xform (
    GLdouble y[4],
    Mat4x4 A,
    GLdouble x[4]
) {
    int i, j;
    for (i=0; i<4; i++) {
        y[i] = A[M4(i,0)] * x[0];
        for (j=1; j<4; j++) {
            y[i] += A[M4(i,j)] * x[j];
        }
    }
}

/* Multiply two 4x4 matrices.
 *
 * C = A * B.
 * N.B. The output matrix must be different from the input.
 */
void
mult (
    Mat4x4 C,
    Mat4x4 A,
    Mat4x4 B
) {
    int i, j, k;
    for (i=0; i<4; i++) {
        for (j=0; j<4; j++) {
            C[M4(i,j)] = A[M4(i,0)] * B[M4(0,j)];
            for (k=1; k<4; k++) {
                C[M4(i,j)] += A[M4(i,k)] * B[M4(k,j)];
            }
        }
    }
}

/* Transpose a matrix.
 */
void
transpose (
    Mat4x4 At,
    Mat4x4 A
) {
    int i, j;
    for (i=0; i<4; i++) {
        for (j=0; j<4; j++) {
            At[M4(i,j)] = A[M4(j,i)];
        }
    }
}

/* Compute the inverse of a 4x4 matrix.
 *
 * From an algorithm by V. Strassen, 1969, _Numerishe Mathematik_, vol. 13,
 * pp. 354-356.  60 multiplies, 24 additions, 10 subtractions, 8 negations, 
 * 2 divisions, 48 assignments, _0_ branches
 *
 * This implementation by Scott McCaskill.
 */
void 
invert ( 
    Mat4x4 C,   /* C is the 4x4 destination matrix */
    Mat4x4 A    /* A is the 4x4 source matrix (to be inverted) */
) {
   Mat2x2 r1, r2, r3, r4, r5, r6, r7;
   GLdouble one_over_det;

   /*
    * a11 is the 2x2 matrix in the upper left quadrant of A
    * a12 is the 2x2 matrix in the upper right quadrant of A
    * a21 is the 2x2 matrix in the lower left quadrant of A
    * a22 is the 2x2 matrix in the lower right quadrant of A
    * similarly, cXX are the 2x2 quadrants of the destination matrix
    */

   /* R1 = inverse( a11 ) */
   one_over_det = 1.0 / 
        ( ( A[M4(0,0)] * A[M4(1,1)] ) - ( A[M4(1,0)] * A[M4(0,1)] ) );
   r1[M2(0,0)] = one_over_det *  A[M4(1,1)];
   r1[M2(0,1)] = one_over_det * -A[M4(0,1)];
   r1[M2(1,0)] = one_over_det * -A[M4(1,0)];
   r1[M2(1,1)] = one_over_det *  A[M4(0,0)];

   /* R2 = a21 x R1 */
   r2[M2(0,0)] = A[M4(2,0)] * r1[M2(0,0)] + A[M4(2,1)] * r1[M2(1,0)];
   r2[M2(0,1)] = A[M4(2,0)] * r1[M2(0,1)] + A[M4(2,1)] * r1[M2(1,1)];
   r2[M2(1,0)] = A[M4(3,0)] * r1[M2(0,0)] + A[M4(3,1)] * r1[M2(1,0)];
   r2[M2(1,1)] = A[M4(3,0)] * r1[M2(0,1)] + A[M4(3,1)] * r1[M2(1,1)];

   /* R3 = R1 x a12 */
   r3[M2(0,0)] = r1[M2(0,0)] * A[M4(0,2)] + r1[M2(0,1)] * A[M4(1,2)];
   r3[M2(0,1)] = r1[M2(0,0)] * A[M4(0,3)] + r1[M2(0,1)] * A[M4(1,3)];
   r3[M2(1,0)] = r1[M2(1,0)] * A[M4(0,2)] + r1[M2(1,1)] * A[M4(1,2)];
   r3[M2(1,1)] = r1[M2(1,0)] * A[M4(0,3)] + r1[M2(1,1)] * A[M4(1,3)];

   /* R4 = a21 x R3 */
   r4[M2(0,0)] = A[M4(2,0)] * r3[M2(0,0)] + A[M4(2,1)] * r3[M2(1,0)];
   r4[M2(0,1)] = A[M4(2,0)] * r3[M2(0,1)] + A[M4(2,1)] * r3[M2(1,1)];
   r4[M2(1,0)] = A[M4(3,0)] * r3[M2(0,0)] + A[M4(3,1)] * r3[M2(1,0)];
   r4[M2(1,1)] = A[M4(3,0)] * r3[M2(0,1)] + A[M4(3,1)] * r3[M2(1,1)];

   /* R5 = R4 - a22 */
   r5[M2(0,0)] = r4[M2(0,0)] - A[M4(2,2)];
   r5[M2(0,1)] = r4[M2(0,1)] - A[M4(2,3)];
   r5[M2(1,0)] = r4[M2(1,0)] - A[M4(3,2)];
   r5[M2(1,1)] = r4[M2(1,1)] - A[M4(3,3)];

   /* R6 = inverse( R5 ) */
   one_over_det = 1.0 / 
        ( ( r5[M2(0,0)] * r5[M2(1,1)] ) - ( r5[M2(1,0)] * r5[M2(0,1)] ) );
   r6[M2(0,0)] = one_over_det *  r5[M2(1,1)];
   r6[M2(0,1)] = one_over_det * -r5[M2(0,1)];
   r6[M2(1,0)] = one_over_det * -r5[M2(1,0)];
   r6[M2(1,1)] = one_over_det *  r5[M2(0,0)];

   /* c12 = R3 x R6 */
   C[M4(0,2)] = r3[M2(0,0)] * r6[M2(0,0)] + r3[M2(0,1)] * r6[M2(1,0)];
   C[M4(0,3)] = r3[M2(0,0)] * r6[M2(0,1)] + r3[M2(0,1)] * r6[M2(1,1)];
   C[M4(1,2)] = r3[M2(1,0)] * r6[M2(0,0)] + r3[M2(1,1)] * r6[M2(1,0)];
   C[M4(1,3)] = r3[M2(1,0)] * r6[M2(0,1)] + r3[M2(1,1)] * r6[M2(1,1)];

   /* c21 = R6 x R2 */
   C[M4(2,0)] = r6[M2(0,0)] * r2[M2(0,0)] + r6[M2(0,1)] * r2[M2(1,0)];
   C[M4(2,1)] = r6[M2(0,0)] * r2[M2(0,1)] + r6[M2(0,1)] * r2[M2(1,1)];
   C[M4(3,0)] = r6[M2(1,0)] * r2[M2(0,0)] + r6[M2(1,1)] * r2[M2(1,0)];
   C[M4(3,1)] = r6[M2(1,0)] * r2[M2(0,1)] + r6[M2(1,1)] * r2[M2(1,1)];

   /* R7 = R3 x c21 */
   r7[M2(0,0)] = r3[M2(0,0)] * C[M4(2,0)] + r3[M2(0,1)] * C[M4(3,0)];
   r7[M2(0,1)] = r3[M2(0,0)] * C[M4(2,1)] + r3[M2(0,1)] * C[M4(3,1)];
   r7[M2(1,0)] = r3[M2(1,0)] * C[M4(2,0)] + r3[M2(1,1)] * C[M4(3,0)];
   r7[M2(1,1)] = r3[M2(1,0)] * C[M4(2,1)] + r3[M2(1,1)] * C[M4(3,1)];

   /* c11 = R1 - R7 */
   C[M4(0,0)] = r1[M2(0,0)] - r7[M2(0,0)];
   C[M4(0,1)] = r1[M2(0,1)] - r7[M2(0,1)];
   C[M4(1,0)] = r1[M2(1,0)] - r7[M2(1,0)];
   C[M4(1,1)] = r1[M2(1,1)] - r7[M2(1,1)];

   /* c22 = -R6 */
   C[M4(2,2)] = -r6[M2(0,0)];
   C[M4(2,3)] = -r6[M2(0,1)];
   C[M4(3,2)] = -r6[M2(1,0)];
   C[M4(3,3)] = -r6[M2(1,1)];
}

void
print_mat (
    FILE* fp,
    Mat4x4 M
) {
    int i, j;
    for (i=0; i<4; i++) {
        fprintf(fp,"[ ");
        for (j=0; j<4; j++) {
            fprintf(fp,"%10.10f ",M[M4(i,j)]);
        }
        fprintf(fp,"]\n");
    }
    fprintf(fp,"{ ");
    for (i=0; i<16; i++) {
        fprintf(fp,"%10.10f ",M[i]);
    }
    fprintf(fp,"}\n");
}

void
print_vec (
    FILE* fp,
    GLdouble v[4]
) {
    int i;
    fprintf(fp,"[ ");
    for (i=0; i<4; i++) {
        fprintf(fp,"%10.10f ",v[i]);
    }
    fprintf(fp,"] => {");
    for (i=0; i<3; i++) {
        fprintf(fp,"%10.10f ",v[i]/v[3]);
    }
    fprintf(fp,"}\n");
}

/************************************************************************
 * Scene Definitions
 ************************************************************************/

/* Features common to all scenes.
 */
void
scene_common() {

    glMatrixMode(GL_MODELVIEW);
    glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,black);

    /* Ground Plane (drawn as squashed cube, actually)
     */
    glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,white);
    glMaterialfv(GL_FRONT_AND_BACK,GL_SPECULAR,black);
    glPushMatrix();
        glScaled(10.0, 0.04, 10.0);
        glutSolidCube(1.0);
    glPopMatrix();

    /* Occluders (outside camera radius!)
     */
    glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,blue);
    glMaterialfv(GL_FRONT_AND_BACK,GL_SPECULAR,black);

    glPushMatrix();
        glTranslated(0.0, 20.0, 0.0);
        glutSolidCube(0.2);
    glPopMatrix();
    glPushMatrix();
        glTranslated(0.22, 20.0, 0.22);
        glutSolidCube(0.2);
    glPopMatrix();
    glPushMatrix();
        glTranslated(-0.22, 20.0, 0.22);
        glutSolidCube(0.2);
    glPopMatrix();
    glPushMatrix();
        glTranslated(-0.22, 20.0, -0.22);
        glutSolidCube(0.2);
    glPopMatrix();
    glPushMatrix();
        glTranslated(0.22, 20.0, -0.22);
        glutSolidCube(0.2);
    glPopMatrix();

    glPushMatrix();
	  glRotated(180.0,1.0,0.0,0.0);
	  glPushMatrix();
	      glTranslated(0.0, 20.0, 0.0);
	      glutSolidCube(0.2);
	  glPopMatrix();
	  glPushMatrix();
	      glTranslated(0.22, 20.0, 0.22);
	      glutSolidCube(0.2);
	  glPopMatrix();
	  glPushMatrix();
	      glTranslated(-0.22, 20.0, 0.22);
	      glutSolidCube(0.2);
	  glPopMatrix();
	  glPushMatrix();
	      glTranslated(-0.22, 20.0, -0.22);
	      glutSolidCube(0.2);
	  glPopMatrix();
	  glPushMatrix();
	      glTranslated(0.22, 20.0, -0.22);
	      glutSolidCube(0.2);
	  glPopMatrix();
    glPopMatrix();
}

/* Set up scene display lists.
 * Needs to be called at beginning of program and if the subdivision
 * level changes.
 */
void
init_scenes(
    int *scene_base_p
) {
    int i, j, k;

    if (*scene_base_p < 0) {
        *scene_base_p = glGenLists(NSCENES);
    }

    glEnable(GL_NORMALIZE);

    /* Sphere and cube 1.
     */
    glNewList(*scene_base_p + SCENE_SPHERE_AND_CUBE_1, GL_COMPILE);

        scene_common();

        /* Sphere 
         */
        glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,red);
        glMateriali(GL_FRONT_AND_BACK,GL_SHININESS,100);
        glMaterialfv(GL_FRONT_AND_BACK,GL_SPECULAR,white);
        glPushMatrix();
            glTranslated(1.5, 3, 1.5);
            glutSolidSphere(1.0, subd, subd);
        glPopMatrix();
        glPushMatrix();
	    glRotated(180.0,1.0,0.0,0.0);
            glTranslated(1.5, 3, 1.5);
            glutSolidSphere(1.0, subd, subd);
        glPopMatrix();

        /* Cube 
         */
        glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,green);
        glMateriali(GL_FRONT_AND_BACK,GL_SHININESS,0);
        glMaterialfv(GL_FRONT_AND_BACK,GL_SPECULAR,black);
        glPushMatrix();
            glTranslated(-1.0, 1.0, -1.0);
            glutSolidCube(2.0);
        glPopMatrix();
        glPushMatrix();
	    glRotated(180.0,1.0,0.0,0.0);
            glTranslated(-1.0, 1.0, -1.0);
            glutSolidCube(2.0);
        glPopMatrix();
    glEndList();

    /* Sphere and cube 2.
     */
    glNewList(*scene_base_p + SCENE_SPHERE_AND_CUBE_2, GL_COMPILE);

        scene_common();

        /* Sphere 
         */
        glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,red);
        glMateriali(GL_FRONT_AND_BACK,GL_SHININESS,100);
        glMaterialfv(GL_FRONT_AND_BACK,GL_SPECULAR,white);
        glPushMatrix();
            glTranslated(1.9, 1.2, -0.3);
            glutSolidSphere(1.2, subd, subd);
        glPopMatrix();
        glPushMatrix();
	    glRotated(180.0,1.0,0.0,0.0);
            glTranslated(1.9, 1.2, -0.3);
            glutSolidSphere(1.2, subd, subd);
        glPopMatrix();

        /* Cube 
         */
        glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,green);
        glMateriali(GL_FRONT_AND_BACK,GL_SHININESS,0);
        glMaterialfv(GL_FRONT_AND_BACK,GL_SPECULAR,black);
        glPushMatrix();
            glTranslated(-1.0, 1.0, -1.0);
            glutSolidCube(2.0);
        glPopMatrix();
        glPushMatrix();
	    glRotated(180.0,1.0,0.0,0.0);
            glTranslated(-1.0, 1.0, -1.0);
            glutSolidCube(2.0);
        glPopMatrix();
    glEndList();

    /* Cubes.
     */
    glNewList(*scene_base_p + SCENE_CUBES, GL_COMPILE);

        scene_common();

        /* Cube Grid
         */
        glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,green);
        glMateriali(GL_FRONT_AND_BACK,GL_SHININESS,10);
        glMaterialfv(GL_FRONT_AND_BACK,GL_SPECULAR,white);
        for (i=0; i<5; i++) {
            for (j=0; j<5; j++) {
                for (k=0; k<5; k++) {
                    glPushMatrix();
                        glTranslated(i-2, k+1, j-2);
                        glutSolidCube(0.5);
                    glPopMatrix();
                    glPushMatrix();
	    	        glRotated(180.0,1.0,0.0,0.0);
                        glTranslated(i-2, k+1, j-2);
                        glutSolidCube(0.5);
                    glPopMatrix();
                }
            }
        }
    glEndList();

    /* Grid.
     */
    glNewList(*scene_base_p + SCENE_GRID, GL_COMPILE);

        scene_common();

        /* Cube Grid
         */
        glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,green);
        glMateriali(GL_FRONT_AND_BACK,GL_SHININESS,10);
        glMaterialfv(GL_FRONT_AND_BACK,GL_SPECULAR,white);
        for (i=0; i<30; i++) {
            for (k=0; k<30; k++) {
                glPushMatrix();
                    glTranslated(3.0*(i-15)/15.0, 3.0*(k+1)/15.0, 0);
                    glutSolidCube(0.1);
                glPopMatrix();
                glPushMatrix();
	    	    glRotated(180.0,1.0,0.0,0.0);
                    glTranslated(3.0*(i-15)/15.0, 3.0*(k+1)/15.0, 0);
                    glutSolidCube(0.1);
                glPopMatrix();
            }
        }
    glEndList();

    /* Spheres.
     */
    glNewList(*scene_base_p + SCENE_SPHERES, GL_COMPILE);

        scene_common();

        /* Sphere Grid
         */
        glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,red);
        glMateriali(GL_FRONT_AND_BACK,GL_SHININESS,100);
        glMaterialfv(GL_FRONT_AND_BACK,GL_SPECULAR,white);
        for (i=0; i<5; i++) {
            for (j=0; j<5; j++) {
                for (k=0; k<5; k++) {
                    glPushMatrix();
                        glTranslated(i-2, k+1, j-2);
                        glutSolidSphere(0.4,subd/2,subd/2);
                    glPopMatrix();
                    glPushMatrix();
	    	        glRotated(180.0,1.0,0.0,0.0);
                        glTranslated(i-2, k+1, j-2);
                        glutSolidSphere(0.4,subd/2,subd/2);
                    glPopMatrix();
                }
            }
        }
    glEndList();

    /* Cone and Torus
     */
    glNewList(*scene_base_p + SCENE_CONE_AND_TORUS, GL_COMPILE);

        scene_common();

        /* Cone   
         */
        glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,red);
        glMateriali(GL_FRONT_AND_BACK,GL_SHININESS,100);
        glMaterialfv(GL_FRONT_AND_BACK,GL_SPECULAR,white);
        glPushMatrix();
            glRotated(-90.0, 1.0, 0.0, 0.0);
            glutSolidCone(2.0,4.0,subd,subd);
        glPopMatrix();
        glPushMatrix();
            glRotated(90.0, 1.0, 0.0, 0.0);
            glutSolidCone(2.0,4.0,subd,subd);
        glPopMatrix();

        /* Torus 
         */
        glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,green);
        glMateriali(GL_FRONT_AND_BACK,GL_SHININESS,64);
        glMaterialfv(GL_FRONT_AND_BACK,GL_SPECULAR,white);
        glPushMatrix();
            glTranslated(0.0, 1.5, 0.0);
            glRotated(75.0, 1.0, 0.0, 0.0);
            glutSolidTorus(0.5, 2.0, subd, subd);
        glPopMatrix();
        glPushMatrix();
            glRotated(180.0, 1.0, 0.0, 0.0);
            glTranslated(0.0, 1.5, 0.0);
            glRotated(75.0, 1.0, 0.0, 0.0);
            glutSolidTorus(0.5, 2.0, subd, subd);
        glPopMatrix();
    glEndList();

    /* Teapot
     */
    glNewList(*scene_base_p + SCENE_TEAPOT, GL_COMPILE);

        scene_common();

        glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,yellow);
        glMateriali(GL_FRONT_AND_BACK,GL_SHININESS,70);
        glMaterialfv(GL_FRONT_AND_BACK,GL_SPECULAR,white);
	glFrontFace(GL_CW);
	glDisable(GL_CULL_FACE);
        glPushMatrix();
            glTranslated(0.0, 1.5, 0.0);
	    glutSolidTeapot(1.5);
        glPopMatrix();
        glPushMatrix();
            glRotated(180.0, 1.0, 0.0, 0.0);
            glTranslated(0.0, 1.5, 0.0);
	    glutSolidTeapot(1.5);
        glPopMatrix();
	glFrontFace(GL_CCW);
	glEnable(GL_CULL_FACE);
    glEndList();

    /* Teapots
     */
    glNewList(*scene_base_p + SCENE_TEAPOTS, GL_COMPILE);

        scene_common();

        glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,yellow);
        glMateriali(GL_FRONT_AND_BACK,GL_SHININESS,70);
        glMaterialfv(GL_FRONT_AND_BACK,GL_SPECULAR,white);
	glFrontFace(GL_CW);
	glDisable(GL_CULL_FACE);
        for (i=0; i<5; i+=2) {
            for (j=0; j<5; j+=2) {
	      glPushMatrix();
                  glTranslated(i-2, 1.5, j-2);
		  glRotatef(45, 0, 1, 0);
		  glutSolidTeapot(0.8);
	      glPopMatrix();
	      glPushMatrix();
		  glRotated(180.0, 1.0, 0.0, 0.0);
                  glTranslated(i-2, 1.5, j-2);
		  glRotatef(45, 0, 1, 0);
		  glutSolidTeapot(0.8);
	      glPopMatrix();
	   }
	}
	glFrontFace(GL_CCW);
	glEnable(GL_CULL_FACE);
    glEndList();

    check_error("init_scenes");
}

/* Render scene.
 * Draws chosen scene with current rotation.
 */
void
render_scene(
    int scene_base
) {
    glMatrixMode(GL_MODELVIEW);
    glPushMatrix();
        glCallList(scene+scene_base);
    glPopMatrix();

    check_error("render_scene");
}

/* Set up initial light and eye/scene transformations
 */
void
init_xforms(
) {
    glMatrixMode(GL_MODELVIEW);
    glPushMatrix();
        glLoadIdentity();
        glRotated(50,1.0,0.0,0.0);
        glGetDoublev(GL_MODELVIEW_MATRIX,eye_xform);

        glLoadIdentity();
        glRotated(50,0.0,1.0,0.0);
        glRotated(-60,1.0,0.0,0.0);
        glGetDoublev(GL_MODELVIEW_MATRIX,light_xform);
    glPopMatrix();

    check_error("init_xforms");
}

/************************************************************************
 * Shadow Mesh Geometry
 ************************************************************************/

void
get_light_proj (
    Mat4x4 C,
    int xr,
    int yr,
    GLdouble ze
) {
    Mat4x4 T, M, iM, P, iP, iD, A, B;

    /* Set up translation matrix to move volume slightly away
     * from light source, and over half a pixel
     */
    T[0] = 1.0;
    T[1] = 0.0;
    T[2] = 0.0;
    T[3] = 0.0;

    T[4] = 0.0;
    T[5] = 1.0;
    T[6] = 0.0;
    T[7] = 0.0;

    T[8] = 0.0;
    T[9] = 0.0;
    T[10] = 1.0;
    T[11] = 0.0;

    T[12] = 0.5;
    T[13] = 0.5;
    T[14] = MAX_MAP_Z * ze; 
    T[15] = 1.0;

    /* Compute inverse of NDC -> DC transformation
     */
    iD[0] = 2.0 / (GLdouble)xr;
    iD[1] = 0.0;
    iD[2] = 0.0;
    iD[3] = 0.0;

    iD[4] = 0.0;
    iD[5] = 2.0 / (GLdouble)yr;
    iD[6] = 0.0;
    iD[7] = 0.0;

    iD[8] = 0.0;
    iD[9] = 0.0;
    iD[10] = 2.0 / (GLdouble)MAX_MAP_Z;
    iD[11] = 0.0;

    iD[12] = -1.0;
    iD[13] = -1.0;
    iD[14] = -1.0;
    iD[15] = 1.0;

    /* Get projection matrix and compute its inverse.
     */
    glGetDoublev(GL_PROJECTION_MATRIX,P);
    invert(iP,P);

    /* Get modelview matrix and compute its inverse.
     */
    glGetDoublev(GL_MODELVIEW_MATRIX,M);
    invert(iM,M);

    /* Combine inverse matrices to get matrix to transform
     * back into world space
     */
    mult(A,iD,T);
    mult(B,iP,A);
    mult(C,iM,B);

    check_error("get_light_proj");
}

    
/* Reconstruct everything into a shadow map, using
 * triangle strips for maximum performance.
 */
void
recon_all () {
    int ix, iy;

    /* Reconstruct everything into a shadow map, using
     * triangle strips for maximum performance.
     */
    for (iy=0; iy < map_yr-1; iy++) {
        glBegin(GL_TRIANGLE_STRIP);
            VERTEX(0,iy+1);
            VERTEX(0,iy);
            map_code[iy][0] = CAT_U_11;
            for (ix=1; ix < map_xr; ix++) {
                map_code[iy][ix] = CAT_U_11;
                VERTEX(ix,iy+1);
                VERTEX(ix,iy);
            }
        glEnd();
    }
}

                /* "Quads" are actually reconstructed as triangles to have 
                 * some control over how the quads are split, which is 
                 * important to know when we generate the cap.
                 *
                 * Here we just split them one way; use a "flipped"
                 * mode to test for the best way.
                 */
void
recon_quads () {
    int ix, iy;
    map_code_type c;
#ifdef RLE
    map_code_type next_ix;
#endif

    glBegin(GL_TRIANGLES);
        for (iy=1; iy<map_yr-2; iy++) {
#ifdef RLE
            c = map_code[iy][0];
            next_ix = (c >> INDEX_SHIFT);
            ix = 0;
            while (next_ix) {
                ix += next_ix;
                c = map_code[iy][ix];
                next_ix = (c >> INDEX_SHIFT);
#else
            for (ix=1; ix<map_xr-2; ix++) {
                c = map_code[iy][ix];
#endif
        
                /* If there is an edge mark on any boundary of
                 * a cell, generate a quad covering that cell
                 */
                if (EDGE_MASK & c) {
                    map_code[iy][ix] = c | CAT_U_11;

                    VERTEX(ix,iy);
                    VERTEX(ix+1,iy);
                    VERTEX(ix+1,iy+1);

                    VERTEX(ix,iy);
                    VERTEX(ix+1,iy+1);
                    VERTEX(ix,iy+1);
                }
            }
        }
    glEnd();
}

/* We actually reconstruct quads as triangles to have some
 * control over how the quads are split, which is 
 * important to know when we generate the cap.
 *
 * Here we split quads in the way that maximizes the
 * difference in the orientation of the two triangles.
 */
void
recon_quads_flip () {
    int ix, iy;
    map_code_type c;
#ifdef RLE
    map_code_type next_ix;
#endif
    map_z_type dz0, dz1;

    glBegin(GL_TRIANGLES);

        for (iy=1; iy<map_yr-2; iy++) {
#ifdef RLE
            c = map_code[iy][0];
            next_ix = (c >> INDEX_SHIFT);
            ix = 0;
            while (next_ix) {
                ix += next_ix;
                c = map_code[iy][ix];
                next_ix = (c >> INDEX_SHIFT);
#else
            for (ix=1; ix<map_xr-2; ix++) {
                c = map_code[iy][ix];
#endif
                /* If there is an edge mark on any boundary of
                 * a cell, generate a quad covering that cell
                 */
                if (EDGE_MASK & c) {

                    /* Compute magnitudes of diagonal deltas
                     */
                    dz0 = map_z[iy+1][ix+1] - map_z[iy][ix];
                    if (dz0 < 0) dz0 = -dz0;
                    dz1 = map_z[iy+1][ix] - map_z[iy][ix+1];
                    if (dz1 < 0) dz1 = -dz1;

                    /* Pick decomposition that cuts across 
                     * maximum delta z
                     */
                    if (dz0 < dz1) {
                        /* Use upwards decomposition
                         */
                        map_code[iy][ix] = c | CAT_U_11;

                        VERTEX(ix,iy);
                        VERTEX(ix+1,iy);
                        VERTEX(ix+1,iy+1);

                        VERTEX(ix,iy);
                        VERTEX(ix+1,iy+1);
                        VERTEX(ix,iy+1);
                    } else {
                        /* Use downwards decomposition
                         */
                        map_code[iy][ix] = c | CAT_D_11;

                        VERTEX(ix,iy);
                        VERTEX(ix+1,iy);
                        VERTEX(ix,iy+1);

                        VERTEX(ix,iy+1);
                        VERTEX(ix+1,iy);
                        VERTEX(ix+1,iy+1);
                    }
                }
            }
        }
    glEnd();
} 

/* Generate independent triangles using table lookup.
 * 
 * We do not test for flip direction.
 */
void
recon_tris () {
    int ix, iy;
    map_code_type c;
#ifdef RLE
    map_code_type next_ix;
#endif

    glBegin(GL_TRIANGLES);
        for (iy=1; iy<map_yr-2; iy++) {
#ifdef RLE
            c = map_code[iy][0];
            next_ix = (c >> INDEX_SHIFT);
            ix = 0;
            while (next_ix) {
                ix += next_ix;
                c = map_code[iy][ix];
                next_ix = (c >> INDEX_SHIFT);
#else
            for (ix=1; ix<map_xr-2; ix++) {
                c = map_code[iy][ix];
#endif

                /* For each possible edge configuration, 
                 * generate zero, one, or two triangles.
                 * Also overwrite edge map with reconstruction
                 * category (just edge code for now...)
                 */
                switch (EDGE_MASK & c) {
                    case EDGE(0,0,0,0): 
                        break;

                    case EDGE(1,0,1,0): 
                        map_code[iy][ix] = c | CAT_D_01;

                        VERTEX(ix,iy);
                        VERTEX(ix+1,iy);
                        VERTEX(ix,iy+1);
                        break;

                    case EDGE(0,1,1,0): 
                        map_code[iy][ix] = c | CAT_U_01;

                        VERTEX(ix,iy);
                        VERTEX(ix+1,iy);
                        VERTEX(ix+1,iy+1);
                        break;

                    case EDGE(1,0,0,1): 
                        map_code[iy][ix] = c | CAT_U_10;

                        VERTEX(ix,iy);
                        VERTEX(ix+1,iy+1);
                        VERTEX(ix,iy+1);
                        break;

                    case EDGE(0,1,0,1): 
                        map_code[iy][ix] = c | CAT_D_10;

                        VERTEX(ix,iy+1);
                        VERTEX(ix+1,iy);
                        VERTEX(ix+1,iy+1);
                        break;

                    case EDGE(1,0,0,0): 
                    case EDGE(0,1,0,0): 
                    case EDGE(0,0,1,0): 
                    case EDGE(0,0,0,1): 

                    case EDGE(1,1,0,0): 
                    case EDGE(1,1,1,0): 
                    case EDGE(1,1,0,1): 
                    case EDGE(0,0,1,1): 
                    case EDGE(1,0,1,1): 
                    case EDGE(0,1,1,1): 
                    case EDGE(1,1,1,1): 
                        map_code[iy][ix] = c | CAT_U_11;

                        VERTEX(ix,iy);
                        VERTEX(ix+1,iy+1);
                        VERTEX(ix,iy+1);

                        VERTEX(ix,iy);
                        VERTEX(ix+1,iy);
                        VERTEX(ix+1,iy+1);
                        break;
                }
            }
        }
    glEnd();
} 

/* Generate independent triangles using table lookup.
 * 
 * We test for flip direction to improve the quality of the
 * mesh when quads are rendered.   
 */
void
recon_tris_flip () {
    int ix, iy;
    map_code_type c;
#ifdef RLE
    map_code_type next_ix;
#endif
    map_z_type dz0, dz1;

    glBegin(GL_TRIANGLES);
        for (iy=1; iy<map_yr-2; iy++) {
#ifdef RLE
            c = map_code[iy][0];
            next_ix = (c >> INDEX_SHIFT);
            ix = 0;
            while (next_ix) {
                ix += next_ix;
                c = map_code[iy][ix];
                next_ix = (c >> INDEX_SHIFT);
#else
            for (ix=1; ix<map_xr-2; ix++) {
                c = map_code[iy][ix];
#endif

                /* For each possible edge configuration, 
                 * generate zero, one, or two triangles.
                 * Also overwrite edge map with reconstruction
                 * category (just edge code for now...)
                 */
                switch (EDGE_MASK & c) {
                    case EDGE(0,0,0,0): 
                        break;

                    case EDGE(1,0,1,0): 
                        map_code[iy][ix] = c | CAT_D_01;

                        VERTEX(ix,iy);
                        VERTEX(ix+1,iy);
                        VERTEX(ix,iy+1);
                        break;

                    case EDGE(0,1,1,0): 
                        map_code[iy][ix] = c | CAT_U_01;

                        VERTEX(ix,iy);
                        VERTEX(ix+1,iy);
                        VERTEX(ix+1,iy+1);
                        break;

                    case EDGE(1,0,0,1): 
                        map_code[iy][ix] = c | CAT_U_10;

                        VERTEX(ix,iy);
                        VERTEX(ix+1,iy+1);
                        VERTEX(ix,iy+1);
                        break;

                    case EDGE(0,1,0,1): 
                        map_code[iy][ix] = c | CAT_D_10;

                        VERTEX(ix,iy+1);
                        VERTEX(ix+1,iy);
                        VERTEX(ix+1,iy+1);
                        break;

                    case EDGE(1,0,0,0): 
                    case EDGE(0,1,0,0): 
                    case EDGE(0,0,1,0): 
                    case EDGE(0,0,0,1): 

                    case EDGE(1,1,0,0): 
                    case EDGE(1,1,1,0): 
                    case EDGE(1,1,0,1): 
                    case EDGE(0,0,1,1): 
                    case EDGE(1,0,1,1): 
                    case EDGE(0,1,1,1): 
                    case EDGE(1,1,1,1): 

                        /* Compute magnitudes of diagonal deltas
                         */
                        dz0 = map_z[iy+1][ix+1] - map_z[iy][ix];
                        if (dz0 < 0) dz0 = -dz0;
                        dz1 = map_z[iy+1][ix] - map_z[iy][ix+1];
                        if (dz1 < 0) dz1 = -dz1;

                        /* Pick decomposition that cuts across 
                         * maximum delta z
                         */
                        if (dz0 < dz1) {
                            /* Use upwards decomposition
                             */
                            map_code[iy][ix] = c | CAT_U_11;

                            VERTEX(ix,iy);
                            VERTEX(ix+1,iy);
                            VERTEX(ix+1,iy+1);

                            VERTEX(ix,iy);
                            VERTEX(ix+1,iy+1);
                            VERTEX(ix,iy+1);
                        } else {
                            /* Use downwards decomposition
                             */
                            map_code[iy][ix] = c | CAT_D_11;

                            VERTEX(ix,iy);
                            VERTEX(ix+1,iy);
                            VERTEX(ix,iy+1);

                            VERTEX(ix,iy+1);
                            VERTEX(ix+1,iy);
                            VERTEX(ix+1,iy+1);
                        }
                        break;
                }
            }
        }
    glEnd();
} 

/* WARNING!   The following vectorization code is incomplete and buggy.
 */

/* Manage storage for convex hulls used during vectorization.
 */
typedef struct sHullPoint {
    int x;
    int y;
} HullPoint;

typedef struct sHull {
    HullPoint* p;
    int i;
    int m;
} Hull;

#define HULL_BLOCK 32

/* Extend a convex hull with a new point.   
 */
static void
extend_hull (
    Hull* A,
    Hull* B,
    int x2,
    int y2, 
    int s
) {
    int i, x0, y0, x1, y1;

    /* Allocate more memory for hull point storage if needed
     */
    if (A->i == A->m) {
        if (A->p) {
            A->m *= 2;
            A->p = (HullPoint *)realloc(A->p,(A->m)*sizeof(HullPoint));
        } else {
            A->m = HULL_BLOCK;
            A->p = (HullPoint *)malloc(HULL_BLOCK*sizeof(HullPoint));
        }
    }

    /* Add new point to hull
     */
    switch (A->i) {
        case 0:
        case 1:
            A->p[A->i].x = x2;
            A->p[A->i].y = y2;
            A->i += 1;
            break;
        default: 
            /* Scan back through old hull points, detecting and deleting
             * points that are inside the new hull
             */
            i = A->i;
            if (s > 0) {
                do {
                    x0 = A->p[i-1].x;
                    y0 = A->p[i-1].y;

                    x1 = A->p[i-2].x;
                    y1 = A->p[i-2].y;

                    i--;
                } while (0 < (x1-x0)*(y2-y1) - (y1-y0)*(x2-x1) && i >= 2);
            } else {
                do {
                    x0 = A->p[i-1].x;
                    y0 = A->p[i-1].y;

                    x1 = A->p[i-2].x;
                    y1 = A->p[i-2].y;

                    i--;
                } while (0 > (x1-x0)*(y2-y1) - (y1-y0)*(x2-x1) && i >= 2);
            }

            /* Store the new point, implicitly deleting
             * any points scanned over.
             */
            A->p[i].x = x2;
            A->p[i].y = y2;

            A->i = i+1;
    }

    /* Check if the last segment added interferes with the last
     * segment of a sister hull.
     */
}


typedef struct sFacet {
    map_code_type code;   /* allowable steps */
    map_code_type prim;   /* primary step direction step */
    map_code_type last;   /* last step direction */
    int error_area;       /* area of error */
    int len;
    int ex0;
    int ey0;
    int ex1;
    int ey1;
    int x0;
    int y0;
    int x1;
    int y1;
    int x2;
    int y2;
    int x3;
    int y3;
} Facet;

/* This is an implied parameter of the functions below; only one facet is
 * generated at a time.
 */
static Facet facet;

/* Print facet for debugging purposes
 */
static void
print_facet () {
    fprintf(stderr,"[");
    PEDGE(stderr,facet.code);
    fprintf(stderr," (%d,%d),(%d,%d),(%d,%d),(%d,%d)]",
        facet.x0, facet.y0, 
        facet.x1, facet.y1,
        facet.x2, facet.y2,
        facet.x3, facet.y3
    );
}

/* Make facet inactive.
 */
static void
clear_facet () {
    facet.code = 0;
}

/* Start vectorization facet on entry to new cell.
 * The vector thus started lies along the entry edge and is
 * designed to extend into the current cell.
 */
static void
start_facet (
    map_code_type entry,
    int x, int y
) {
    /* Compute starting positions for each possible entry and constraints
     * on future directions of edge steps
     */
    facet.code = EDGE_MASK & ~entry;
    switch (entry) {
        case EDGE(0,0,0,1):
            facet.x0 = x+1;
            facet.y0 = y+1;
            facet.ex0 = 2*x + 1;
            facet.ey0 = 2*y + 2;
            facet.x1 = x;
            facet.y1 = y+1;
            break;
        case EDGE(0,0,1,0):
            facet.x0 = x;
            facet.y0 = y;
            facet.ex0 = 2*x + 1;
            facet.ey0 = 2*y;
            facet.x1 = x+1;
            facet.y1 = y;
            break;
        case EDGE(0,1,0,0):
            facet.x0 = x+1;
            facet.y0 = y;
            facet.ex0 = 2*x + 2;
            facet.ey0 = 2*y + 1;
            facet.x1 = x+1;
            facet.y1 = y+1;
            break;
        case EDGE(1,0,0,0):
            facet.x0 = x;
            facet.y0 = y+1;
            facet.ex0 = 2*x;
            facet.ey0 = 2*y + 1;
            facet.x1 = x;
            facet.y1 = y;
            break;
        default:
            fprintf(stderr,
                "Bad entry in start_facet: 0x%x\n",
                entry);
            exit(1);
    }
    facet.len = 0;
    facet.error_area = 0;
}

/* End facet at entry to current cell.     
 * The entry code and coordinates given are both those
 * of the CURRENT cell.    The facet thereby terminated does not
 * enter the cell, but lies along the entry edge.
 */
static void
end_facet (
    map_code_type entry,
    int x, int y
) {
    /* Can't terminate facets that aren't open 
     */
    if (!facet.code) return;

    /* Figure out endpoint info
     */
    switch (entry) {
        case EDGE(0,0,0,1):
            facet.x2 = x;
            facet.y2 = y+1;
            facet.x3 = x+1;
            facet.y3 = y+1;
            break;
        case EDGE(0,0,1,0):
            facet.x2 = x+1;
            facet.y2 = y;
            facet.x3 = x;
            facet.y3 = y;
            break;
        case EDGE(0,1,0,0):
            facet.x2 = x+1;
            facet.y2 = y+1;
            facet.x3 = x+1;
            facet.y3 = y;
            break;
        case EDGE(1,0,0,0):
            facet.x2 = x;
            facet.y2 = y;
            facet.x3 = x;
            facet.y3 = y+1;
            break;
        default:
            fprintf(stderr,
                "Bad move in end_facet: 0x%x\n",
                entry);
            exit(1);
    }

    /* DEBUG
    print_facet();
     */

    /* Render facet.
     */
    glBegin(GL_TRIANGLE_STRIP);
        VERTEX(facet.x0, facet.y0);
        VERTEX(facet.x1, facet.y1);
        VERTEX(facet.x3, facet.y3);
        VERTEX(facet.x2, facet.y2);
    glEnd();
}

#define LINE_TOL 10

/* Attempt to extend facet.   The coordinates and edge code given
 * are those of the CURRENT cell.    The entry code given is the
 * entry code to this cell.
 */
static int
extend_facet (
    map_code_type c,            /* edge code of current cell */
    map_code_type entry,        /* entry code into current cell */
    int x, int y
) {
    int extended;
    
    /* Always mark cell so not visited again
     */
    map_code[y][x] = (c | CAT_VEC);

    /* Only extend facet if active, otherwise start a new one
     * upon this entry.
     */
    if (!facet.code) {
        start_facet(entry,x,y);
        return 0;
    }

    /* Test if current vector should be extended or not
     */
    facet.code &= ~entry;
    extended = (facet.code) ? 1 : 0;

    if (extended && entry == facet.last) {
        /* Repeated entry detected.   This is only permitted for the 
         * ``primary'' step direction.
         */ 
        if (facet.prim) {
            if (facet.prim != entry) { 
                extended = 0;
            }
        } else {
            facet.prim = entry;
        }
    }
    facet.last = entry;

    extended = 1;

    if (extended) {
        int ex, ey;

        /* Update error area term
         */
        switch (entry) {
            case EDGE(0,0,0,1):
                ex = 2*x + 1;
                ey = 2*y + 2;
                break;
            case EDGE(0,0,1,0):
                ex = 2*x + 1;
                ey = 2*y;
                break;
            case EDGE(0,1,0,0):
                ex = 2*x + 2;
                ey = 2*y + 1;
                break;
            case EDGE(1,0,0,0):
                ex = 2*x;
                ey = 2*y + 1;
                break;
            default:
                fprintf(stderr,
                    "Bad entry in extend_facet: 0x%x\n",
                    entry);
                exit(1);
        }
        if (facet.len == 0) {
            facet.ex1 = ex;
            facet.ey1 = ey;
        } else {
            facet.error_area += (facet.ex1 - facet.ex0)*(ey - facet.ey1)
                              - (facet.ey1 - facet.ey0)*(ex - facet.ex1);

            if (facet.error_area > LINE_TOL || facet.error_area < -LINE_TOL) {
                extended = 0;
            } else {
                facet.ex1 = ex;
                facet.ey1 = ey;
            }
        }
    }

    /* Count length of facet
     */
    facet.len++;

    /* Either extend or terminate facet and start a new one
     */
    if (!extended) {
        end_facet(entry,x,y);
        start_facet(entry,x,y);
    }

    return extended;
}

static void
quad_flip (
    map_code_type c,
    int x, int y
) {
    map_z_type dz0, dz1;

    dz0 = map_z[y+1][x+1] - map_z[y][x];
    if (dz0 < 0) dz0 = -dz0;
    dz1 = map_z[y+1][x] - map_z[y][x+1];
    if (dz1 < 0) dz1 = -dz1;

    if (dz0 < dz1) {
        glBegin(GL_TRIANGLES);
            TRI(x,y, x+1,y, x+1,y+1);
            TRI(x,y, x+1,y+1, x,y+1);
        glEnd();
        map_code[y][x] = (c | CAT_U_11);
    } else {
        glBegin(GL_TRIANGLES);
            TRI(x,y, x+1,y, x,y+1);
            TRI(x,y+1, x+1,y, x+1,y+1);
        glEnd();
        map_code[y][x] = (c | CAT_D_11);
    }
}

/* Walk around contours in the mesh and vectorize the
 * edges, producing a cleaner shadow volume with fewer polygons.
 */
void
recon_vectorize () {
    int ix, iy;
#ifdef RLE
    map_code_type next_ix;
#endif
    int x, y;
    int done, cont;
    map_code_type c, cn, entry, egress;

    for (iy=1; iy<map_yr-2; iy++) {
#ifdef RLE
        c = map_code[iy][0];
        next_ix = (c >> INDEX_SHIFT);
        ix = 0;
        while (next_ix) {
            ix += next_ix;
            c = map_code[iy][ix];
            next_ix = (c >> INDEX_SHIFT);
#else
        for (ix=1; ix<map_xr-2; ix++) {
            c = map_code[iy][ix];
#endif
            /* Found edge mark; as long as is not 
             * ``used up'', walk around contour, segmenting
             * it into vectors to be rendered.
             */
            if ((CAT_MASK & c) == CAT_NONE) {
                x = ix;
                y = iy;
                entry = 0;
                done = 0;
                clear_facet();

                do {
                    /* Find all edge codes which are not the entry edge
                     */
                    egress = c & ~entry;

                    switch (EDGE_MASK & c) {
                        case EDGE(1,0,1,0): 

                            if (!entry) {
                                entry = EDGE(0,0,1,0);
                                egress = EDGE(1,0,0,0);
                            }
                            extend_facet(c,entry,x,y);

                            if ((EDGE(1,0,0,0) & egress) && x > 2) {
                                cn = map_code[y][x-1];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(0,1,0,0);
                                    x--;
                                    break;
                                }
                            }  

                            if ((EDGE(0,0,1,0) & egress) && (y > 2)) {
                                cn = map_code[y-1][x];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(0,0,0,1);
                                    y--;
                                    break;
                                }
                            }  

                            if (EDGE(0,0,1,0) == entry) {
                                end_facet(EDGE(0,1,0,0),x-1,y);
                            } else {
                                end_facet(EDGE(0,0,0,1),x,y-1);
                            }

                            done = 1;
                            break;

                        case EDGE(0,1,1,0): 

                            if (!entry) {
                                entry = EDGE(0,0,1,0);
                                egress = EDGE(0,1,0,0);
                            }
                            extend_facet(c,entry,x,y);

                            if ((EDGE(0,1,0,0) & egress) && (x < map_xr-2)) {
                                cn = map_code[y][x+1];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(1,0,0,0);
                                    x++;
                                    break;
                                }
                            }  

                            if ((EDGE(0,0,1,0) & egress) && (y > 2)) {
                                cn = map_code[y-1][x];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(0,0,0,1);
                                    y--;
                                    break;
                                }
                            }  

                            if (EDGE(0,0,1,0) == entry) {
                                end_facet(EDGE(1,0,0,0),x+1,y);
                            } else {
                                end_facet(EDGE(0,0,0,1),x,y-1);
                            }

                            done = 1;
                            break;

                        case EDGE(1,0,0,1): 

                            if (!entry) {
                                entry = EDGE(0,0,0,1);
                                egress = EDGE(1,0,0,0);
                            }
                            extend_facet(c,entry,x,y);

                            if ((EDGE(1,0,0,0) & egress) && x > 2) {
                                cn = map_code[y][x-1];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(0,1,0,0);
                                    x--;
                                    break;
                                }
                            } 

                            if ((EDGE(0,0,0,1) & egress) && (y < map_yr-2)) {
                                cn = map_code[y+1][x];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(0,0,1,0);
                                    y++;
                                    break;
                                }
                            } 

                            if (EDGE(0,0,0,1) == entry) {
                                end_facet(EDGE(0,1,0,0),x-1,y);
                            } else {
                                end_facet(EDGE(0,0,1,0),x,y+1);
                            }

                            done = 1;
                            break;

                        case EDGE(0,1,0,1): 

                            if (!entry) {
                                entry = EDGE(0,0,0,1);
                                egress = EDGE(0,1,0,0);
                            }
                            extend_facet(c,entry,x,y);

                            if ((EDGE(0,1,0,0) & egress) && (x < map_xr-2)) {
                                cn = map_code[y][x+1];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(1,0,0,0);
                                    x++;
                                    break;
                                }
                            }  

                            if ((EDGE(0,0,0,1) & egress) && (y < map_yr-2)) {
                                cn = map_code[y+1][x];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(0,0,1,0);
                                    y++;
                                    break;
                                }
                            } 

                            if (EDGE(0,0,0,1) == entry) {
                                end_facet(EDGE(1,0,0,0),x+1,y);
                            } else {
                                end_facet(EDGE(0,0,1,0),x,y+1);
                            }

                            done = 1;
                            break;

                        case EDGE(1,0,0,0): 

                            if (!entry) {
                                entry = EDGE(0,1,0,0);
                                egress = EDGE(1,0,0,0);
                            }
                            extend_facet(c,entry,x,y);

                            if ((EDGE(1,0,0,0) & egress) && x > 2) {
                                cn = map_code[y][x-1];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(0,1,0,0);
                                    x--;
                                    break;
                                }
                            } 

                            if (EDGE(0,1,0,0) == entry) {
                                end_facet(EDGE(0,1,0,0),x+1,y);
                            } else {
                                end_facet(EDGE(1,0,0,0),x-1,y);
                            }

                            done = 1;
                            break;

                        case EDGE(0,1,0,0): 

                            if (!entry) {
                                entry = EDGE(1,0,0,0);
                                egress = EDGE(0,1,0,0);
                            }
                            extend_facet(c,entry,x,y);

                            if ((EDGE(0,1,0,0) & egress) && (x < map_xr-2)) {
                                cn = map_code[y][x+1];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(1,0,0,0);
                                    x++;
                                    break;
                                }
                            }  

                            if (EDGE(1,0,0,0) == entry) {
                                end_facet(EDGE(1,0,0,0),x+1,y);
                            } else {
                                end_facet(EDGE(0,1,0,0),x-1,y);
                            }

                            done = 1;
                            break;

                        case EDGE(0,0,1,0): 

                            if (!entry) {
                                entry = EDGE(0,0,0,1);
                                egress = EDGE(0,0,1,0);
                            }
                            extend_facet(c,entry,x,y);

                            if ((EDGE(0,0,1,0) & egress) && (y > 2)) {
                                cn = map_code[y-1][x];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(0,0,0,1);
                                    y--;
                                    break;
                                }
                            }  

                            if (EDGE(0,0,0,1) == entry) {
                                end_facet(EDGE(0,0,0,1),x,y-1);
                            } else {
                                end_facet(EDGE(0,0,1,0),x,y+1);
                            }

                            done = 1;
                            break;

                        case EDGE(0,0,0,1): 

                            if (!entry) {
                                entry = EDGE(0,0,1,0);
                                egress = EDGE(0,0,0,1);
                            }
                            extend_facet(c,entry,x,y);

                            if ((EDGE(0,0,0,1) & egress) && (y < map_yr-2)) {
                                cn = map_code[y+1][x];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(0,0,1,0);
                                    y++;
                                    break;
                                }
                            } 

                            if (EDGE(0,0,1,0) == entry) {
                                end_facet(EDGE(0,0,1,0),x,y+1);
                            } else {
                                end_facet(EDGE(0,0,0,1),x,y-1);
                            }

                            done = 1;
                            break;

                        case EDGE(1,1,0,0): 

                            if (!entry) {
                                entry = EDGE(0,1,0,0);
                                egress = EDGE(1,0,0,0);
                            }
                            extend_facet(c,entry,x,y);

                            if ((EDGE(1,0,0,0) & egress) && x > 2) {
                                cn = map_code[y][x-1];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(0,1,0,0);
                                    x--;
                                    break;
                                }
                            } 

                            if ((EDGE(0,1,0,0) & egress) && (x < map_xr-2)) {
                                cn = map_code[y][x+1];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(1,0,0,0);
                                    x++;
                                    break;
                                }
                            }  

                            if (EDGE(0,1,0,0) == entry) {
                                end_facet(EDGE(0,1,0,0),x-1,y);
                            } else {
                                end_facet(EDGE(1,0,0,0),x+1,y);
                            }

                            done = 1;
                            break;

                        case EDGE(0,0,1,1): 

                            if (!entry) {
                                entry = EDGE(0,0,0,1);
                                egress = EDGE(0,0,1,0);
                            }
                            extend_facet(c,entry,x,y);

                            if ((EDGE(0,0,1,0) & egress) && (y > 2)) {
                                cn = map_code[y-1][x];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(0,0,0,1);
                                    y--;
                                    break;
                                }
                            }  

                            if ((EDGE(0,0,0,1) & egress) && (y < map_yr-2)) {
                                cn = map_code[y+1][x];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(0,0,1,0);
                                    y++;
                                    break;
                                }
                            } 

                            if (EDGE(0,0,0,1) == entry) {
                                end_facet(EDGE(0,0,0,1),x,y-1);
                            } else {
                                end_facet(EDGE(0,0,1,0),x,y+1);
                            }

                            done = 1;
                            break;

                        case EDGE(1,1,1,0): 
                        case EDGE(1,1,0,1): 
                        case EDGE(1,0,1,1): 
                        case EDGE(0,1,1,1): 
                        case EDGE(1,1,1,1): 

                            /* Terminate vectorization at 3 or 4-way
                             * branches and reconstruct cell with a quad
                             */
                            end_facet(entry,x,y);
                            quad_flip(c,x,y);

                            /* Pick exit and try to restart vectorization
                             */
                            if ((EDGE(1,0,0,0) & egress) && x > 2) {
                                cn = map_code[y][x-1];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(0,1,0,0);
                                    x--;
                                    break;
                                }
                            } 
                            if ((EDGE(0,1,0,0) & egress) && (x < map_xr-2)) {
                                cn = map_code[y][x+1];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(1,0,0,0);
                                    x++;
                                    break;
                                }
                            }  
                            if ((EDGE(0,0,1,0) & egress) && (y > 2)) {
                                cn = map_code[y-1][x];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(0,0,0,1);
                                    y--;
                                    break;
                                }
                            }  
                            if ((EDGE(0,0,0,1) & egress) && (y < map_yr-2)) {
                                cn = map_code[y+1][x];
                                if ((CAT_MASK & cn) == CAT_NONE) {
                                    c = cn;
                                    entry = EDGE(0,0,1,0);
                                    y++;
                                    break;
                                }
                            } 

                            done = 1;
                            break;
                        default:
                            fprintf(stderr,
                                "Bad edge code: 0x%x\n",EDGE_MASK & c);
                            exit(1);
                    }
                } while (!done);
            }
        }
    }
} 

/* Reconstruct the shadow geometry in world space.
 * The matrix L should be a transformation from light device space
 * into world space, and sm should be a display list id to bind the
 * shadow volume to.
 */ 
void
recon_volume (
    Mat4x4 L,
    int sm
) {
    timer_start(&timer_recon_volume);
    glNewList(sm, GL_COMPILE);
        glMatrixMode(GL_MODELVIEW);
        glPushMatrix();
            glMultMatrixd(L);

            switch (recon_mode) {
                case RECON_ALL: 
                    recon_all();
                    break;
                case RECON_QUADS: 
                    recon_quads();
                    break;
                case RECON_QUADS_FLIP:  
                    recon_quads_flip();
                    break;
                case RECON_TRIS: 
                    recon_tris();
                    break;
                case RECON_TRIS_FLIP: 
                    recon_tris_flip();
                    break;
                case RECON_VECTORIZE: 
                    recon_vectorize();
                    break;
            } 
        glPopMatrix();
    glEndList();
    check_error("recon_volume");
    timer_stop(&timer_recon_volume);
}

void
get_view_proj (
    Mat4x4 P,
    Mat4x4 iP,
    Mat4x4 M,
    Mat4x4 iM,
    Mat4x4 V,
    Mat4x4 iV,
    int xr,
    int yr
) {
    /* Get projection matrix and compute its inverse.
     */
    glGetDoublev(GL_PROJECTION_MATRIX,P);
    invert(iP,P);

    /* Get modelview matrix and compute its inverse.
     */
    glGetDoublev(GL_MODELVIEW_MATRIX,M);
    invert(iM,M);

    /* Combine forward matrices to get matrix to transform
     * from world space into the view NDC space
     */
    mult(V,P,M);

    /* Combine inverse matrices to get matrix to transform
     * from view DCS back into world space
     */
    mult(iV,iM,iP);

    check_error("get_view_proj");
}

/* Generate a cap.
 *
 * The matrix L is the light device to world coordinate system
 * transformation.   This function should be called when the viewing
 * transformation is available in the pipeline.    
 */
void
cap (
    Mat4x4 L,
    int xr, 
    int yr
) {
    Mat4x4 oldP;
    Mat4x4 V, iV, P, iP, M, iM, T, Tt, iT;
    GLdouble v[4];
    GLdouble q[4], p[4];
    GLdouble p_ll[4], p_lr[4], p_ul[4], p_ur[4];
    GLdouble x, y, t;
    GLint xi, yi, lo_x, lo_y, hi_x, hi_y;
    map_code_type c;
    int full;

    if (cap_mode == CAP_OFF) return;

    timer_start(&timer_cap);

    /* Get projection matrix, do a pretranslation in z to avoid
     * cap clipping, and multiply the P matrix back in, then
     * reget it and compute its inverse. 
     */
    glGetDoublev(GL_PROJECTION_MATRIX,oldP);
    /* avoid pushing stack, usually not very deep; just restore at end */

    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    glTranslated(0.0, 0.0, 1.0);   /* put near plane at z=0 */
    glMultMatrixd(oldP);
    glGetDoublev(GL_PROJECTION_MATRIX,P);

    invert(iP,P);

    /* Get modelview matrix and compute its inverse.
     */
    glGetDoublev(GL_MODELVIEW_MATRIX,M);
    invert(iM,M);

    /* Combine forward matrices to get matrix to transform
     * from world space into the view NDC space
     */
    mult(V,P,M);

    /* Combine inverse matrices to get matrix to transform
     * from view DCS back into world space
     */
    mult(iV,iM,iP);

    /* Get the viewing transformation and its inverse;
     * compute the overall transformation from the light DCS
     * to the view NDCS.   The V matrix is the product of 
     * P (the projection matrix) and M (the modelview matrix)
     */
    mult(T,V,L);

    /* Transform near plane equation.   Note that normal is set up
     * to point INTO the scene, AWAY the eye, so values with a 
     * negative distance from the plane equation are CLIPPED.
     * Note: transform a plane by the inverse transpose,
     * an inverse of the inverse is just the matrix...
     */
    transpose(Tt,T);
    q[0] = 0;
    q[1] = 0;
    q[2] = 1;
    q[3] = 0;
    xform(p,Tt,q);

    /* Transform the corners of the near plane into the light DCS
     */
    invert(iT,T);
    v[0] = -1;
    v[1] = -1;
    v[2] =  0;
    v[3] =  1;
    xform(p_ll,iT,v);

    v[0] =  1;
    v[1] = -1;
    v[2] =  0;
    v[3] =  1;
    xform(p_lr,iT,v);

    v[0] = -1;
    v[1] =  1;
    v[2] =  0;
    v[3] =  1;
    xform(p_ul,iT,v);

    v[0] =  1;
    v[1] =  1;
    v[2] =  0;
    v[3] =  1;
    xform(p_ur,iT,v);

    /* Fit an integer-valued bounding box around the intersections
     * of rays passing through the eye near plane with the 
     * light near plane.
     */
    x = p_ll[0]/p_ll[3];
    xi = (GLint)(x - 0.5);
    lo_x = xi;
    xi = (GLint)(x + 1.5);
    hi_x = xi;

    x = p_lr[0]/p_lr[3];
    xi = (GLint)(x - 0.5);
    if (xi < lo_x) lo_x = xi;
    xi = (GLint)(x + 1.6);
    if (xi > hi_x) hi_x = xi;

    x = p_ul[0]/p_ul[3];
    xi = (GLint)(x - 0.5);
    if (xi < lo_x) lo_x = xi;
    xi = (GLint)(x + 1.6);
    if (xi > hi_x) hi_x = xi;

    x = p_ur[0]/p_ur[3];
    xi = (GLint)(x - 0.5);
    if (xi < lo_x) lo_x = xi;
    xi = (GLint)(x + 1.6);
    if (xi > hi_x) hi_x = xi;

    y = p_ll[1]/p_ll[3];
    yi = (GLint)(y - 0.5);
    lo_y = yi;
    yi = (GLint)(y + 1.6);
    hi_y = yi;

    y = p_lr[1]/p_lr[3];
    yi = (GLint)(y - 0.5);
    if (yi < lo_y) lo_y = yi;
    yi = (GLint)(y + 1.6);
    if (yi > hi_y) hi_y = yi;

    y = p_ul[1]/p_ul[3];
    yi = (GLint)(y - 0.5);
    if (yi < lo_y) lo_y = yi;
    yi = (GLint)(y + 1.6);
    if (yi > hi_y) hi_y = yi;

    y = p_ur[1]/p_ur[3];
    yi = (GLint)(y - 0.5);
    if (yi < lo_y) lo_y = yi;
    yi = (GLint)(y + 1.6);
    if (yi > hi_y) hi_y = yi;

    /* Clip minmax coords to active area of light map
     */
    if (lo_x < 1) lo_x = 1;
    if (hi_x > xr - 1) hi_x = xr - 1;
    if (lo_y < 1) lo_y = 1;
    if (hi_y > yr - 1) hi_y = yr - 1;

    /* Don't do cap if can't see near plane in light view
     */
    if (lo_x < hi_x && lo_y < hi_y) {

        if (p[2] < 0.0) {
            p[0] = -p[0];
            p[1] = -p[1];
            p[2] = -p[2];
            p[3] = -p[3];
        }

        /* Apply the plane equation to all samples in the region
         * and update the cap codes (not changing the cat codes!)
         */
        full = 1;
        for (yi = lo_y; yi <= hi_y; yi++) {
            for (xi = lo_x; xi <= hi_x; xi++) {
                /* Clear any existing cap codes in region
                 */
                map_code[yi][xi] &= ~CAP_MASK;

                /* Do plane test and set appropriate cap codes
                 */
                PLANE(p,xi,yi,t);
                if (t < 0) {
                    map_code[yi][xi]     |= CAP_00;
                    map_code[yi-1][xi]   |= CAP_01;
                    map_code[yi][xi-1]   |= CAP_10;
                    map_code[yi-1][xi-1] |= CAP_11;
                } else {
                    full = 0;
                }
            }
        }

        /* If ALL the samples we test are behind the near plane,
         * the near plane is fully in shadow and we can just render 
         * a single quad to invert the entire stencil buffer.
         */
        if (full) {

            glMatrixMode(GL_MODELVIEW);
            glPushMatrix();
            glLoadIdentity();

            glMatrixMode(GL_PROJECTION);
            glLoadIdentity();

            /* Note: let clipper reduce to cover all edges for sure...
             */
            glBegin(GL_QUADS);
                glVertex3i(-2,-2,0);
                glVertex3i( 2,-2,0);
                glVertex3i( 2, 2,0);
                glVertex3i(-2, 2,0);
            glEnd();

            glMatrixMode(GL_MODELVIEW);
            glPopMatrix();

        } else {

            /* Reconstruction of the cap is done in the light
             * device coordinate system, so we have to transform
             * back out...   but V already in effect...
             */
            glMatrixMode(GL_MODELVIEW);
            glPushMatrix();
            glMultMatrixd(L);

                /* Scan all cap codes and generate cap polygons.
                 */
                for (yi = lo_y; yi < hi_y; yi++) {
                    for (xi = lo_x; xi < hi_x; xi++) {

                        c = map_code[yi][xi];

                        switch (CAP_MASK & c) {
                            case CAP(0,0,0,0):
                                break;
                            case CAP(1,1,1,1):

                                glBegin(GL_QUADS);

                                CAP_VERTEX(p,xi,yi+1);
                                CAP_VERTEX(p,xi,yi);

                                /* Try to grab an entire run at once, and
                                 * render as a single quad.
                                 */
                                xi++;
                                c = map_code[yi][xi];
                                while (
                                    (xi < hi_x) &&
                                    CAP(1,1,1,1) == (CAP_MASK & c)) {

                                    xi++;
                                    c = map_code[yi][xi];
                                }

                                CAP_VERTEX(p,xi,yi);
                                CAP_VERTEX(p,xi,yi+1);

                                glEnd();

                                /* Need to back up one...
                                 */
                                xi--;
                                break;
                            case CAP(0,0,0,1):
                                switch (CAT_MASK & c) {
                                    case CAT_NONE:
                                    case CAT_D_01:
                                    case CAT_D_10:
                                    case CAT_D_11:
                                        glBegin(GL_TRIANGLES);
                                        CAP_VERTEX(p,xi+1,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi+1,yi);
                                        glEnd();
                                        break;
                                    case CAT_U_01:
                                    case CAT_U_10:
                                    case CAT_U_11:
                                        glBegin(GL_TRIANGLES);
                                        CAP_VERTEX(p,xi+1,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi,yi);

                                        CAP_VERTEX(p,xi+1,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi,yi);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi+1,yi);
                                        glEnd();
                                        break;
                                    default:
                                        fprintf(stderr,
                                            "Unexpected cat code 0x%x\n", c);
                                        exit(1);
                                }
                                break;
                            case CAP(0,0,1,0):
                                switch (CAT_MASK & c) {
                                    case CAT_D_01:
                                    case CAT_D_10:
                                    case CAT_D_11:
                                        glBegin(GL_TRIANGLES);
                                        CAP_VERTEX(p,xi+1,yi);
                                        CLIP_VERTEX(p,xi+1,yi,xi,yi);
                                        CLIP_VERTEX(p,xi+1,yi,xi,yi+1);

                                        CAP_VERTEX(p,xi+1,yi);
                                        CLIP_VERTEX(p,xi+1,yi,xi,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi,xi+1,yi+1);
                                        glEnd();
                                        break;
                                    case CAT_NONE:
                                    case CAT_U_01:
                                    case CAT_U_10:
                                    case CAT_U_11:
                                        glBegin(GL_TRIANGLES);
                                        CAP_VERTEX(p,xi+1,yi);
                                        CLIP_VERTEX(p,xi+1,yi,xi,yi);
                                        CLIP_VERTEX(p,xi+1,yi,xi+1,yi+1);
                                        glEnd();
                                        break;
                                    default:
                                        fprintf(stderr,
                                            "Unexpected cat code: 0x%x\n", c);
                                        exit(1);
                                }
                                break;
                            case CAP(0,1,0,0):
                                switch (CAT_MASK & c) {
                                    case CAT_D_01:
                                    case CAT_D_10:
                                    case CAT_D_11:
                                        glBegin(GL_TRIANGLES);
                                        CAP_VERTEX(p,xi,yi+1);
                                        CLIP_VERTEX(p,xi,yi+1,xi+1,yi);
                                        CLIP_VERTEX(p,xi,yi+1,xi+1,yi+1);

                                        CAP_VERTEX(p,xi,yi+1);
                                        CLIP_VERTEX(p,xi,yi+1,xi,yi);
                                        CLIP_VERTEX(p,xi,yi+1,xi+1,yi);
                                        glEnd();
                                        break;
                                    case CAT_NONE:
                                    case CAT_U_01:
                                    case CAT_U_10:
                                    case CAT_U_11:
                                        glBegin(GL_TRIANGLES);
                                        CAP_VERTEX(p,xi,yi+1);
                                        CLIP_VERTEX(p,xi,yi+1,xi,yi);
                                        CLIP_VERTEX(p,xi,yi+1,xi+1,yi+1);
                                        glEnd();
                                        break;
                                    default:
                                        fprintf(stderr,
                                            "Unexpected cat code: 0x%x\n", c);
                                        exit(1);
                                }
                                break;
                            case CAP(1,0,0,0):
                                switch (CAT_MASK & c) {
                                    case CAT_NONE:
                                    case CAT_D_01:
                                    case CAT_D_10:
                                    case CAT_D_11:
                                        glBegin(GL_TRIANGLES);
                                        CAP_VERTEX(p,xi,yi);
                                        CLIP_VERTEX(p,xi,yi,xi+1,yi);
                                        CLIP_VERTEX(p,xi,yi,xi,yi+1);
                                        glEnd();
                                        break;
                                    case CAT_U_01:
                                    case CAT_U_10:
                                    case CAT_U_11:
                                        glBegin(GL_TRIANGLES);
                                        CAP_VERTEX(p,xi,yi);
                                        CLIP_VERTEX(p,xi,yi,xi+1,yi);
                                        CLIP_VERTEX(p,xi,yi,xi+1,yi+1);

                                        CAP_VERTEX(p,xi,yi);
                                        CLIP_VERTEX(p,xi,yi,xi+1,yi+1);
                                        CLIP_VERTEX(p,xi,yi,xi,yi+1);
                                        glEnd();
                                        break;
                                    default:
                                        fprintf(stderr,
                                            "Unexpected cat code: 0x%x\n", c);
                                        exit(1);
                                }
                                break;
                            case CAP(1,1,1,0):
                                switch (CAT_MASK & c) {
                                    case CAT_NONE:
                                    case CAT_D_01:
                                    case CAT_D_10:
                                    case CAT_D_11:
                                        glBegin(GL_POLYGON);
                                        CAP_VERTEX(p,xi,yi);
                                        CAP_VERTEX(p,xi,yi+1);
                                        CLIP_VERTEX(p,xi,yi+1,xi+1,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi,xi+1,yi+1);
                                        CAP_VERTEX(p,xi+1,yi);
                                        glEnd();
                                        break;
                                    case CAT_U_01:
                                    case CAT_U_10:
                                    case CAT_U_11:
                                        glBegin(GL_QUADS);
                                        CAP_VERTEX(p,xi,yi);
                                        CAP_VERTEX(p,xi,yi+1);
                                        CLIP_VERTEX(p,xi,yi+1,xi+1,yi+1);
                                        CLIP_VERTEX(p,xi,yi,xi+1,yi+1);

                                        CAP_VERTEX(p,xi,yi);
                                        CLIP_VERTEX(p,xi,yi,xi+1,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi,xi+1,yi+1);
                                        CAP_VERTEX(p,xi+1,yi);
                                        glEnd();
                                        break;
                                    default:
                                        fprintf(stderr,
                                            "Unexpected cat code: 0x%x\n", c);
                                        exit(1);
                                }
                                break;
                            case CAP(1,1,0,1):
                                switch (CAT_MASK & c) {
                                    case CAT_D_01:
                                    case CAT_D_10:
                                    case CAT_D_11:
                                        glBegin(GL_QUADS);
                                        CAP_VERTEX(p,xi,yi+1);
                                        CAP_VERTEX(p,xi,yi);
                                        CLIP_VERTEX(p,xi,yi,xi+1,yi);
                                        CLIP_VERTEX(p,xi,yi+1,xi+1,yi);

                                        CAP_VERTEX(p,xi,yi+1);
                                        CLIP_VERTEX(p,xi,yi+1,xi+1,yi);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi+1,yi);
                                        CAP_VERTEX(p,xi+1,yi+1);
                                        glEnd();
                                        break;
                                    case CAT_NONE:
                                    case CAT_U_01:
                                    case CAT_U_10:
                                    case CAT_U_11:
                                        glBegin(GL_POLYGON);
                                        CAP_VERTEX(p,xi,yi+1);
                                        CAP_VERTEX(p,xi,yi);
                                        CLIP_VERTEX(p,xi,yi,xi+1,yi);
                                        CLIP_VERTEX(p,xi+1,yi,xi+1,yi+1);
                                        CAP_VERTEX(p,xi+1,yi+1);
                                        glEnd();
                                        break;
                                    default:
                                        fprintf(stderr,
                                            "Unexpected cat code: 0x%x\n", c);
                                        exit(1);
                                }
                                break;
                            case CAP(1,0,1,1):
                                switch (CAT_MASK & c) {
                                    case CAT_D_01:
                                    case CAT_D_10:
                                    case CAT_D_11:
                                        glBegin(GL_QUADS);
                                        CAP_VERTEX(p,xi,yi);
                                        CLIP_VERTEX(p,xi,yi,xi,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi,xi,yi+1);
                                        CAP_VERTEX(p,xi+1,yi);

                                        CAP_VERTEX(p,xi+1,yi);
                                        CLIP_VERTEX(p,xi+1,yi,xi,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi,yi+1);
                                        CAP_VERTEX(p,xi+1,yi+1);
                                        glEnd();
                                        break;
                                    case CAT_NONE:
                                    case CAT_U_01:
                                    case CAT_U_10:
                                    case CAT_U_11:
                                        glBegin(GL_POLYGON);
                                        CAP_VERTEX(p,xi,yi);
                                        CLIP_VERTEX(p,xi,yi,xi,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi,yi+1);
                                        CAP_VERTEX(p,xi+1,yi+1);
                                        CAP_VERTEX(p,xi+1,yi);
                                        glEnd();
                                        break;
                                    default:
                                        fprintf(stderr,
                                            "Unexpected cat code: 0x%x\n",c);
                                        exit(1);
                                }
                                break;
                            case CAP(0,1,1,1):
                                switch (CAT_MASK & c) {
                                    case CAT_NONE:
                                    case CAT_D_01:
                                    case CAT_D_10:
                                    case CAT_D_11:
                                        glBegin(GL_POLYGON);
                                        CLIP_VERTEX(p,xi,yi,xi,yi+1);
                                        CAP_VERTEX(p,xi,yi+1);
                                        CAP_VERTEX(p,xi+1,yi+1);
                                        CAP_VERTEX(p,xi+1,yi);
                                        CLIP_VERTEX(p,xi+1,yi,xi,yi);
                                        glEnd();
                                        break;
                                    case CAT_U_01:
                                    case CAT_U_10:
                                    case CAT_U_11:
                                        glBegin(GL_QUADS);
                                        CAP_VERTEX(p,xi,yi+1);
                                        CAP_VERTEX(p,xi+1,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi,yi);
                                        CLIP_VERTEX(p,xi,yi+1,xi,yi);

                                        CAP_VERTEX(p,xi+1,yi+1);
                                        CAP_VERTEX(p,xi+1,yi);
                                        CLIP_VERTEX(p,xi+1,yi,xi,yi);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi,yi);
                                        glEnd();
                                        break;
                                    default:
                                        fprintf(stderr,
                                            "Unexpected cat code: 0x%x\n", c);
                                        exit(1);
                                }
                                break;
                            case CAP(0,0,1,1):
                                switch (CAT_MASK & c) {
                                    case CAT_NONE:
                                    case CAT_D_01:
                                    case CAT_D_10:
                                    case CAT_D_11:
                                        glBegin(GL_TRIANGLES);
                                        CAP_VERTEX(p,xi+1,yi);
                                        CLIP_VERTEX(p,xi+1,yi,xi,yi);
                                        CLIP_VERTEX(p,xi+1,yi,xi,yi+1);
                                        glEnd();

                                        glBegin(GL_QUADS);
                                        CAP_VERTEX(p,xi+1,yi);
                                        CLIP_VERTEX(p,xi+1,yi,xi,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi,yi+1);
                                        CAP_VERTEX(p,xi+1,yi+1);
                                        glEnd();
                                        break;
                                    case CAT_U_01:
                                    case CAT_U_10:
                                    case CAT_U_11:
                                        glBegin(GL_TRIANGLES);
                                        CAP_VERTEX(p,xi+1,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi,yi);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi,yi+1);
                                        glEnd();

                                        glBegin(GL_QUADS);
                                        CAP_VERTEX(p,xi+1,yi);
                                        CLIP_VERTEX(p,xi+1,yi,xi,yi);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi,yi);
                                        CAP_VERTEX(p,xi+1,yi+1);
                                        glEnd();
                                        break;
                                    default:
                                        fprintf(stderr,
                                            "Unexpected cat code: 0x%x\n", c);
                                        exit(1);
                                }
                                break;
                            case CAP(0,1,0,1):
                                switch (CAT_MASK & c) {
                                    case CAT_NONE:
                                    case CAT_D_01:
                                    case CAT_D_10:
                                    case CAT_D_11:
                                        glBegin(GL_TRIANGLES);
                                        CAP_VERTEX(p,xi,yi+1);
                                        CLIP_VERTEX(p,xi,yi+1,xi+1,yi);
                                        CLIP_VERTEX(p,xi,yi+1,xi,yi);
                                        glEnd();

                                        glBegin(GL_QUADS);
                                        CAP_VERTEX(p,xi,yi+1);
                                        CLIP_VERTEX(p,xi,yi+1,xi+1,yi);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi+1,yi);
                                        CAP_VERTEX(p,xi+1,yi+1);
                                        glEnd();
                                        break;
                                    case CAT_U_01:
                                    case CAT_U_10:
                                    case CAT_U_11:
                                        glBegin(GL_TRIANGLES);
                                        CAP_VERTEX(p,xi+1,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi,yi);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi+1,yi);
                                        glEnd();

                                        glBegin(GL_QUADS);
                                        CAP_VERTEX(p,xi,yi+1);
                                        CLIP_VERTEX(p,xi,yi+1,xi,yi);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi,yi);
                                        CAP_VERTEX(p,xi+1,yi+1);
                                        glEnd();
                                        break;
                                    default:
                                        fprintf(stderr,
                                            "Unexpected cat code: 0x%x\n", c);
                                        exit(1);
                                }
                                break;
                            case CAP(1,0,1,0):
                                switch (CAT_MASK & c) {
                                    case CAT_NONE:
                                    case CAT_D_01:
                                    case CAT_D_10:
                                    case CAT_D_11:
                                        glBegin(GL_TRIANGLES);
                                        CAP_VERTEX(p,xi+1,yi);
                                        CLIP_VERTEX(p,xi+1,yi,xi,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi,xi+1,yi+1);
                                        glEnd();

                                        glBegin(GL_QUADS);
                                        CAP_VERTEX(p,xi,yi);
                                        CLIP_VERTEX(p,xi,yi,xi,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi,xi,yi+1);
                                        CAP_VERTEX(p,xi+1,yi);
                                        glEnd();
                                        break;
                                    case CAT_U_01:
                                    case CAT_U_10:
                                    case CAT_U_11:
                                        glBegin(GL_TRIANGLES);
                                        CAP_VERTEX(p,xi,yi);
                                        CLIP_VERTEX(p,xi,yi,xi,yi+1);
                                        CLIP_VERTEX(p,xi,yi,xi+1,yi+1);
                                        glEnd();

                                        glBegin(GL_QUADS);
                                        CAP_VERTEX(p,xi,yi);
                                        CLIP_VERTEX(p,xi,yi,xi+1,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi,xi+1,yi+1);
                                        CAP_VERTEX(p,xi+1,yi);
                                        glEnd();
                                        break;
                                    default:
                                        fprintf(stderr,
                                            "Unexpected cat code: 0x%x\n", c);
                                        exit(1);
                                }
                                break;
                            case CAP(1,1,0,0):
                                switch (CAT_MASK & c) {
                                    case CAT_NONE:
                                    case CAT_D_01:
                                    case CAT_D_10:
                                    case CAT_D_11:
                                        glBegin(GL_TRIANGLES);
                                        CAP_VERTEX(p,xi,yi+1);
                                        CLIP_VERTEX(p,xi,yi+1,xi+1,yi+1);
                                        CLIP_VERTEX(p,xi,yi+1,xi+1,yi);
                                        glEnd();

                                        glBegin(GL_QUADS);
                                        CAP_VERTEX(p,xi,yi+1);
                                        CLIP_VERTEX(p,xi,yi+1,xi+1,yi);
                                        CLIP_VERTEX(p,xi,yi,xi+1,yi);
                                        CAP_VERTEX(p,xi,yi);
                                        glEnd();
                                        break;
                                    case CAT_U_01:
                                    case CAT_U_10:
                                    case CAT_U_11:
                                        glBegin(GL_TRIANGLES);
                                        CAP_VERTEX(p,xi,yi);
                                        CLIP_VERTEX(p,xi,yi,xi+1,yi);
                                        CLIP_VERTEX(p,xi,yi,xi+1,yi+1);
                                        glEnd();

                                        glBegin(GL_QUADS);
                                        CAP_VERTEX(p,xi,yi);
                                        CLIP_VERTEX(p,xi,yi,xi+1,yi+1);
                                        CLIP_VERTEX(p,xi,yi+1,xi+1,yi+1);
                                        CAP_VERTEX(p,xi,yi+1);
                                        glEnd();
                                        break;
                                    default:
                                        fprintf(stderr,
                                            "Unexpected cat code: 0x%x\n", c);
                                        exit(1);
                                }
                                break;
                            case CAP(0,1,1,0):
                                switch (CAT_MASK & c) {
                                    case CAT_NONE:
                                    case CAT_D_01:
                                    case CAT_D_10:
                                    case CAT_D_11:
                                        glBegin(GL_QUADS);
                                        CAP_VERTEX(p,xi,yi+1);
                                        CAP_VERTEX(p,xi+1,yi);
                                        CLIP_VERTEX(p,xi+1,yi,xi+1,yi+1);
                                        CLIP_VERTEX(p,xi,yi+1,xi+1,yi+1);

                                        CAP_VERTEX(p,xi,yi+1);
                                        CAP_VERTEX(p,xi+1,yi);
                                        CLIP_VERTEX(p,xi+1,yi,xi,yi);
                                        CLIP_VERTEX(p,xi,yi+1,xi,yi);
                                        glEnd();
                                        break;
                                    case CAT_U_01:
                                    case CAT_U_10:
                                    case CAT_U_11:
                                        glBegin(GL_TRIANGLES);
                                        CAP_VERTEX(p,xi,yi+1);
                                        CLIP_VERTEX(p,xi,yi+1,xi+1,yi+1);
                                        CLIP_VERTEX(p,xi,yi+1,xi,yi);

                                        CAP_VERTEX(p,xi+1,yi);
                                        CLIP_VERTEX(p,xi+1,yi,xi,yi);
                                        CLIP_VERTEX(p,xi+1,yi,xi+1,yi+1);
                                        glEnd();
                                        break;
                                    default:
                                        fprintf(stderr,
                                            "Unexpected cat code: 0x%x\n", c);
                                        exit(1);
                                }
                                break;
                            case CAP(1,0,0,1):
                                switch (CAT_MASK & c) {
                                    case CAT_NONE:
                                    case CAT_D_01:
                                    case CAT_D_10:
                                    case CAT_D_11:
                                        glBegin(GL_TRIANGLES);
                                        CAP_VERTEX(p,xi+1,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi+1,yi);

                                        CAP_VERTEX(p,xi,yi);
                                        CLIP_VERTEX(p,xi,yi,xi+1,yi);
                                        CLIP_VERTEX(p,xi,yi,xi,yi+1);
                                        glEnd();
                                        break;
                                    case CAT_U_01:
                                    case CAT_U_10:
                                    case CAT_U_11:
                                        glBegin(GL_QUADS);
                                        CAP_VERTEX(p,xi,yi);
                                        CAP_VERTEX(p,xi+1,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi,yi+1);
                                        CLIP_VERTEX(p,xi,yi,xi,yi+1);

                                        CAP_VERTEX(p,xi,yi);
                                        CAP_VERTEX(p,xi+1,yi+1);
                                        CLIP_VERTEX(p,xi+1,yi+1,xi+1,yi);
                                        CLIP_VERTEX(p,xi,yi,xi+1,yi);
                                        glEnd();
                                        break;
                                    default:
                                        fprintf(stderr,
                                            "Unexpected cat code: 0x%x\n", c);
                                        exit(1);
                                }
                                break;
                            default:
                                fprintf(stderr, 
                                    "Unexpected cap code: 0x%x\n", c);
                                exit(1);
                        }
                    }
                }
            glMatrixMode(GL_MODELVIEW);
            glPopMatrix();
        }
    }

    /* `pop' back to old projection w/o using (too-shallow) stack */
    glMatrixMode(GL_PROJECTION);
    glLoadMatrixd(oldP);

    check_error("cap");
    timer_stop(&timer_cap);
}

/************************************************************************
 * Utility Functions
 ************************************************************************/

/* The following thresholds are used by the edge detection process.
 * We set the thresholds as fractions of the depth resolution to be
 * somewhat independent of the hardware's capabilities, i.e. the
 * thresholds are geometric distances.   However, if your scenes
 * are more complex than the ones used here, these thresholds may 
 * not be correct.
 */
void
set_thresholds() {
    GLint i;
    long m;

    glGetIntegerv(GL_DEPTH_BITS,&i);
    m = ~((~0uL) << (i-1));
    th_H = (m >> 10) << (33 - i);
    th_L = (m >> 18) << (33 - i);
    th_s = (m >> 15) << (33 - i);

    check_error("set_thresholds");
}

/* The following espilons are used to hide shadow map grunginess.
 * They are functions of one of the thresholds and the depth map
 * precision.
 */
void
set_epsilons () {
    GLint i;
    long m;

    glGetIntegerv(GL_DEPTH_BITS,&i);
    m = ~((~0uL) << (i-1));

    light_ze = 1e-5;
    eye_ze = 5e-4;
    line_ze = 8e-4;

    check_error("set_epsilons");
}

void 
update_window_title()
{
    char title[256];

    /* Reflect size of shadow map in title 
     */
    if (light_window_active) {
        sprintf(title,"shadow [%dx%d] <- [%dx%d]",xr,yr,lxr,lyr);
    } else {
        sprintf(title,"shadow [%dx%d]",xr,yr);
    }
    glutSetWindowTitle(title); 

    check_error("update_window_title");
}

void 
post_edges() {
    edge_update = 1;
}

void 
post_redisplay()
{
    glutSetWindow(light_window); 
    glutPostRedisplay(); 
    glutSetWindow(main_window);
    glutPostRedisplay();
}

void 
map_realloc (
    int xxr,
    int yyr
) {
    int i;
    map_z_type **old_z;
    map_z_type *old_z_body;
    map_code_type **old_code;
    map_code_type *old_code_body;

    /* Avoid doing unnecessary work
     */
    if (map_xr == xxr && map_yr == yyr) return;

    /* Make SURE recompute edges if resize happened for any reason
     */
    edge_update = 1;

    /* Reallocate map storage: z
     */
    old_z = map_z;
    if (old_z) {
        old_z_body = old_z[0];
    } else {
        old_z_body = NULL;
    }
    map_z =
        (map_z_type **)realloc(old_z,yyr*sizeof(map_z_type *));
    map_z[0] = 
        (map_z_type *)realloc(old_z_body,xxr*yyr*sizeof(map_z_type));
    for (i=1; i<yyr; i++) {
        map_z[i] = map_z[0] + i*xxr;
    }

    /* Reallocate edge code map
     */
    old_code = map_code;
    if (old_code) {
        old_code_body = old_code[0];
    } else {
        old_code_body = NULL;
    }
    map_code = 
        (map_code_type **)realloc(old_code,yyr*sizeof(map_code_type *));
    map_code[0] = 
        (map_code_type *)realloc(old_code_body,xxr*yyr*sizeof(map_code_type));
    for (i=1; i<yyr; i++) {
        map_code[i] = map_code[0] + i*xxr;
    }

    /* Record new size of map
     */
    map_xr = xxr;
    map_yr = yyr;
}

/* A simple threshold-of-the-gradient operator, using 
 * horizontal and vertical first differences.
 */
void 
edge_simple () {
    int xi, yi;
    map_z_type dv, dh;
    map_z_type *z00_p;
    map_z_type *z10_p;
    map_z_type *z01_p;
    map_code_type c;
    map_code_type *c00_p;
    map_code_type *c10_p;
    map_code_type *c01_p;
#ifdef RLE
    int last_xi;
#endif

    for (xi=0; xi<map_xr; xi++) {
        map_code[0][xi] = 0;
    }

    for (yi=1; yi<map_yr-1; yi++) {

        map_code[yi][0] = 0;

        for (xi=1; xi<map_xr-1; xi++) {

            /* Compute magnitude of vertical first difference
             */
            dv = map_z[yi+1][xi] - map_z[yi][xi];

            /* Compute magnitude of horizontal first difference
             */
            dh = map_z[yi][xi+1] - map_z[yi][xi];

            /* Compute edge codes 
             */
            c = 0;
            if (dh > th_H || dh < -th_H) {
                map_code[yi-1][xi] |= EDGE_H1;
                c |= EDGE_H0;
            } 
            if (dv > th_H || dv < -th_H) {
                map_code[yi][xi-1] |= EDGE_V1;
                c |= EDGE_V0;
            } 

            /* Write edge code.   Note that this also clears upper
             * bits, even if there is no edge.
             */
            map_code[yi][xi] = c;
        }

#ifdef RLE
        /* Rescan finished scanlines to run-length encode edge marks
         */
        last_xi = 0;
        for (xi=1; xi<map_xr-1; xi++) {
            if (EDGE_MASK & map_code[yi-1][xi]) {
                map_code[yi-1][last_xi] |= 
                    ((xi - last_xi) << INDEX_SHIFT);
                last_xi = xi;
            }
        }
#endif
    }
}

/* Candidate edges that pass a thresholded-magnitude gradient test
 * are subjected to a zero-crossing of the second derivative test
 * (this is not really a full Canny edge detector)
 */
void 
edge_canny () {
    int xi, yi;
    map_z_type dv, dh;
    long s0, s1;
    map_code_type c;
#ifdef RLE
    int last_xi;
#endif

    for (xi=0; xi<map_xr; xi++) {
        map_code[0][xi] = 0;
    }

    for (yi=1; yi<map_yr-2; yi++) {

        map_code[yi][0] = 0;

        for (xi=1; xi<map_xr-2; xi++) {

            /* Check for vertical edge
             */
            if (map_z[yi+1][xi] > map_z[yi][xi]) {

                /* Increasing first vertical difference
                 */
                dv = map_z[yi+1][xi] - map_z[yi][xi];

                if (dv > th_H) {
                    c = EDGE_V0;
                    map_code[yi][xi-1] |= EDGE_V1;
                } else if (dv > th_L) {

                    s0 = (long)map_z[yi+1][xi]
                       - ((long)map_z[yi][xi] << 1)
                       + (long)map_z[yi-1][xi];

                    s1 = (long)map_z[yi+2][xi]  
                       - ((long)map_z[yi+1][xi] << 1)
                       + (long)map_z[yi][xi];

                    if (s1 < -th_s && s0 > th_s) {
                        c = EDGE_V0;
                        map_code[yi][xi-1] |= EDGE_V1;
                    } else {
                        c = 0;
                    } 
                } else {
                    c = 0;
                }

            } else {

                /* Decreasing first vertical difference
                 */
                dv = map_z[yi][xi] - map_z[yi+1][xi];

                if (dv > th_H) {
                    c = EDGE_V0;
                    map_code[yi][xi-1] |= EDGE_V1;
                } else if (dv > th_L) {

                    s0 = (long)map_z[yi+1][xi] 
                       - ((long)map_z[yi][xi] << 1)
                       + (long)map_z[yi-1][xi];

                    s1 = (long)map_z[yi+2][xi]
                       - ((long)map_z[yi+1][xi] << 1)
                       + (long)map_z[yi][xi];

                    if (s1 > th_s && s0 < -th_s) {
                        c = EDGE_V0;
                        map_code[yi][xi-1] |= EDGE_V1;
                    } else {
                        c = 0;
                    }
                } else {
                    c = 0;
                }

            }

            /* Check for horizontal edge
             */
            if (map_z[yi][xi+1] > map_z[yi][xi]) {

                /* Increasing first horizontal difference
                 */
                dh = map_z[yi][xi+1] - map_z[yi][xi];

                if (dh > th_H) {
                    c |= EDGE_H0;
                    map_code[yi-1][xi] |= EDGE_H1;
                } else if (dh > th_L) {

                    s0 = (long)map_z[yi][xi+1] 
                       - ((long)map_z[yi][xi] << 1)
                       + (long)map_z[yi][xi-1];

                    s1 = (long)map_z[yi][xi+2] 
                       - ((long)map_z[yi][xi+1] << 1)
                       + (long)map_z[yi][xi];

                    if (s1 < -th_s && s0 > th_s) {
                        c |= EDGE_H0;
                        map_code[yi-1][xi] |= EDGE_H1;
                    } 
                }
            } else {

                /* Increasing first horizontal difference
                 */
                dh = map_z[yi][xi] - map_z[yi][xi+1];

                if (dh > th_H) {
                    c |= EDGE_H0;
                    map_code[yi-1][xi] |= EDGE_H1;
                } else if (dh > th_L) {

                    s0 = (long)map_z[yi][xi+1] 
                       - ((long)map_z[yi][xi] << 1)
                       + (long)map_z[yi][xi-1];

                    s1 = (long)map_z[yi][xi+2] 
                       - ((long)map_z[yi][xi+1] << 1)
                       + (long)map_z[yi][xi];

                    if (s1 > th_s && s0 < -th_s) {
                        c |= EDGE_H0;
                        map_code[yi-1][xi] |= EDGE_H1;
                    } 
                }
            }

            /* Record edge code
             */
            map_code[yi][xi] = c;
        }

#ifdef RLE
        /* Rescan finished scanlines to run-length encode edge marks
         */
        last_xi = 0;
        for (xi=1; xi<map_xr-2; xi++) {
            if (EDGE_MASK & map_code[yi-1][xi]) {
                map_code[yi-1][last_xi] |= 
                    ((xi - last_xi) << INDEX_SHIFT);
                last_xi = xi;
            }
        }
#endif
    }
}

/* The null edge detector; just marks everything as an edge, so a 
 * full mesh shadow volume will be generated.   For visual and 
 * performance comparison.
 */
void 
edge_all () {
    int xi, yi;
    map_z_type d;

    for (yi=0; yi<map_yr; yi++) {
        for (xi=0; xi<map_xr-1; xi++) {
            map_code[yi][xi] = (1 << INDEX_SHIFT) | EDGE(1,1,1,1);
        }
        map_code[yi][xi] = EDGE(1,1,1,1);
    }
}

void
configure_lightview(
    int xr, 
    int yr
) {
    Mat4x4 inv_light_xform;

    /* Set up perspective projection for light
     */
    glViewport(0, 0, xr, yr);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluPerspective(light_fov, 1.0, light_near, light_far);

    /* Load light position and orientation.
     */
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    gluLookAt(
        light_position[0], light_position[1], light_position[2],  
        0.0, 0.0, 0.0,       
        0, 1, 0    /* up */
    );
    invert(inv_light_xform,light_xform);
    glMultMatrixd(inv_light_xform);

    check_error("configure_lightview");
}

void
compute_edges () {

    timer_start(&timer_render_depth);

    /* Set correct depth/stencil modes
     */
    glEnable(GL_DEPTH_TEST);
    glDepthMask(GL_TRUE);
    glStencilFunc(GL_ALWAYS,0,0);

    /* Clear depth buffer and render scene
     */
    glClear(GL_DEPTH_BUFFER_BIT);
    glMatrixMode(GL_MODELVIEW);
    if (light_window_active) {
        render_scene(light_scene_base);
    } else {
        render_scene(main_scene_base);
    }

    /* Flush OpenGL pipeline (is this necessary?)
     */
    glFinish();
    timer_stop(&timer_render_depth);

    /* Don't seem to be able to read back depth buffer without
     * disabling depth test... 
     */
    timer_start(&timer_read_depth);
    glDisable(GL_DEPTH_TEST);
    glDepthMask(GL_FALSE);

    /* Get depth values
     */
    glPixelTransferi(GL_DEPTH_SCALE,1.0);
    glPixelTransferi(GL_DEPTH_BIAS,0);
    glReadPixels(0,0,map_xr,map_yr,GL_DEPTH_COMPONENT,MAP_Z_TYPE,map_z[0]);
    timer_stop(&timer_read_depth);

    /* Do edge detection
     */
    if (recon_mode != RECON_ALL) {
        timer_start(&timer_edge_detect);
        switch (edge_mode) {
            case EDGE_SIMPLE:
                edge_simple();
                break;
            case EDGE_CANNY:
                edge_canny();
                break;
        }
        timer_stop(&timer_edge_detect);
    }

    check_error("compute_edges");
}

Mat4x4 L;

int
acquire_map (
    int map_win,
    int mxr, 
    int myr,
    int recon_win,
    int recon_id
) {
    map_realloc(mxr,myr);
    if (!edge_update) {
        return 0;
    }

    glutSetWindow(map_win);
    configure_lightview(mxr,myr);
    compute_edges();
    get_light_proj(L,mxr,myr,light_ze);

    glutSetWindow(recon_win);
    recon_volume(L, recon_id);

    edge_update = 0;

    check_error("acquire_map");
    return 1;
}

void 
tb_rotate(
    GLdouble rx, GLdouble ry, GLdouble rz, 
    GLdouble scale 
) {
    int i;
    GLdouble r, ir, angle;

    /* Find the length of the vector.
     */
    r = sqrt(rx*rx + ry*ry + rz*rz);

    /* If the vector has zero length do nothing.
     */
    if (r > -1e-6 && r < 1e-6) {
        return;
    }

    /* Normalize the rotation axis vector.
     */
    ir = 1.0 / r;
    rx *= ir;
    ry *= ir;
    rz *= ir;

    /* Compute the rotation angle
     */
    angle = r * scale * 180.0;  

    /* Apply OpenGL rotation
     */
    glRotated(angle, rx, ry, rz);

    check_error("tb_rotate");
}

/* Computes a vector for a virtual trackball.
 */
void 
tb_gen_vec(
    int xr, int yr,
    int x, int y, 
    GLdouble *vx_p,
    GLdouble *vy_p, 
    GLdouble *vz_p
) {
    int xo, yo;
    GLdouble vx, vy, vz, vz_sq;
    GLdouble r, f;

    /* Compute center of window.
     */
    xo = xr / 2;
    yo = yr / 2;

    /* Compute the radius of the trackball sphere
     */
    r = (xr > yr) ? yr : xr;
    r /= 2.0/0.9;  

    /* Calculate components of rotation axis
     */
    vx = (x - xo) / r; 
    vy = (y - yo) / r; 
    vz_sq = (1.0 - vx*vx - vy*vy);

    /* Adjust if moving around outside virtual trackball.
     */
    if (vz_sq < 0.0) {
       f = sqrt(1.0 - vz_sq);
       vx /= f;
       vy /= f;
       vz = 0.0;
    } else {
       vz = sqrt(vz_sq);
    }

    /* Store result
     */
    *vx_p = vx;
    *vy_p = vy;
    *vz_p = vz;
}

/* Generate an incremental rotation matrix using the cross product,
 * and update the old vector.
 */
void 
tb_inc_vec(
    int xr, int yr, 
    int x, int y, 
    GLdouble *vx_p,
    GLdouble *vy_p, 
    GLdouble *vz_p,
    GLdouble *rx_p,
    GLdouble *ry_p, 
    GLdouble *rz_p
) {
    GLdouble vx, vy, vz;
    GLdouble nvx, nvy, nvz;

    /* Compute current vector.
     */
    tb_gen_vec(xr,yr,x,y,&nvx,&nvy,&nvz);

    /* Get old vector.
     */
    vx = *vx_p;
    vy = *vy_p;
    vz = *vz_p;

    /* Generate incremental rotation vector using cross product.
     */
    *rx_p = vy*nvz - nvy*vz;
    *ry_p = vz*nvx - nvz*vx; 
    *rz_p = vx*nvy - nvx*vy; 

    /* Update `last vector' for next time.
     */
    *vx_p = nvx;
    *vy_p = nvy;
    *vz_p = nvz; 
}

/************************************************************************
 * Demo
 *
 * This function is activated by a timer, and puts the system through
 * its paces automatically, so the program can operate in ``attract''
 * mode.   The timer must be reset in the display callback.
 ************************************************************************/

int demo_ticks = 0;

GLdouble demo_rx_scene, demo_ry_scene, demo_rz_scene;
GLdouble demo_angle_scene;

GLdouble demo_rx_light, demo_ry_light, demo_rz_light;
GLdouble demo_angle_light;

void
demo (
    int state
) {
    GLdouble r;

    if (demo_mode) {
        if (999 == demo_ticks % 1000) {
            print_timers(stdout);
            init_timers();
            scene++;
            if (scene > SCENE_END_DEMO) {
                scene = SCENE_START_DEMO;
            }
        }

        if (0 == demo_ticks % 76) {
            demo_rx_scene = drand48();
            demo_ry_scene = drand48();
            demo_rz_scene = drand48();
            r = sqrt(demo_rx_scene*demo_rx_scene +
                     demo_ry_scene*demo_ry_scene +
                     demo_rz_scene*demo_rz_scene);
            demo_rx_scene /= r;
            demo_ry_scene /= r;
            demo_rz_scene /= r;
            demo_angle_scene = 1.0 + 3.0*drand48();
        }

        if (0 == demo_ticks % 89) {

            demo_rx_light = drand48();
            demo_ry_light = drand48();
            demo_rz_light = drand48();
            r = sqrt(demo_rx_light*demo_rx_light +
                     demo_ry_light*demo_ry_light +
                     demo_rz_light*demo_rz_light);
            demo_rx_light /= r;
            demo_ry_light /= r;
            demo_rz_light /= r;
            demo_angle_light = 1.0 + 3.0*drand48();
        }

        glMatrixMode(GL_MODELVIEW);
        glPushMatrix();
            glLoadMatrixd(eye_xform);
            glRotated(demo_angle_scene,
                demo_rx_scene,demo_ry_scene,demo_rz_scene);
            glGetDoublev(GL_MODELVIEW_MATRIX,eye_xform);

            glLoadMatrixd(light_xform);
            glRotated(demo_angle_light,
                demo_rx_light,demo_ry_light,demo_rz_light);
            glGetDoublev(GL_MODELVIEW_MATRIX,light_xform);
        glPopMatrix();

        demo_ticks++;

        post_redisplay();
        post_edges();

        check_error("demo");
    }
}

/************************************************************************
 * Event Callbacks
 ************************************************************************/

/* Resize event for main window
 */
void 
reshape(
    int width, 
    int height
) {
    xr = width;
    yr = height;
    if (!light_window_active) post_edges();

    /* Algorithms very size-dependent, so clear timers if change size
     * of window
     */
    init_timers();
}

/* Resize event for light window
 */
void 
reshape_light(
    int width, 
    int height
) {
    lxr = width;
    lyr = height;
    post_edges();
    post_redisplay();

    /* Algorithms very size-dependent, so clear timers if change size
     * of window
     */
    init_timers();
}

/* Print out list of keyboard commands
 */
void
keyboard_help () {
    fprintf(stdout,"Keypress Commands:\n");
    fprintf(stdout,"\tH, h, ?, /: list keypress commands\n");
    fprintf(stdout,"\tESC, Q, q: quit\n");
    fprintf(stdout,"\tC, c: close light view\n");
    fprintf(stdout,"\tA, a: demo mode on (animate)\n");
    fprintf(stdout,"\tZ, z: demo off\n");
    fprintf(stdout,"\tE, e: edge map\n");
    fprintf(stdout,"\tD, d: depth map\n");
    fprintf(stdout,"\tS, s: shaded light view\n");
    fprintf(stdout,"\t1, 2, 3, etc: change scene\n");
    fprintf(stdout,"\to, O: no shadows\n");
    fprintf(stdout,"\tp, P: shadow volume\n");
    fprintf(stdout,"\t[: dark shadows (auto ambient off)\n");
    fprintf(stdout,"\t]: ambient shadows (auto ambient on)\n");
    fprintf(stdout,"\t\\: composited shadows\n");
    fprintf(stdout,"\tt: print timing report\n");
    fprintf(stdout,"\tT: timing reports on/off\n");
    fprintf(stdout,"\tr, R: clear timing info\n");
}

/* Keypress events
 */
void
keyboard(
    unsigned char key, 
    int x, 
    int y 
) {
    switch (key) {
        case '?': 
        case '/': 
        case 'h':
        case 'H':
            keyboard_help();
            break;
        case '\x1B': 
        case 'q':
        case 'Q':
            exit(0);
            break;
        case 'c':
        case 'C':
            light_window_active = 0;
            glutSetWindow(light_window);
            glutHideWindow();
            post_edges();
            post_redisplay();
            break;
        case 'a':
        case 'A':
            menu_change = 1;
            demo_mode = 1;
            glutTimerFunc(demo_delay,demo,0);
            break;
        case 'z':
        case 'Z':
            menu_change = 1;
            demo_mode = 0;
            break;
        case 'e':
        case 'E':
            menu_change = 1;
            light_window_active = 1;
            glutSetWindow(light_window);
            glutShowWindow();
            light_window_mode = LW_EDGE;
            post_redisplay();
            post_edges();
            update_menus();
            break;
        case 'd':
        case 'D':
            menu_change = 1;
            light_window_active = 1;
            glutSetWindow(light_window);
            glutShowWindow();
            light_window_mode = LW_DEPTH;
            post_redisplay();
            post_edges();
            update_menus();
            break;
        case 's':
        case 'S':
            menu_change = 1;
            light_window_active = 1;
            glutSetWindow(light_window);
            glutShowWindow();
            light_window_mode = LW_SHADED;
            post_redisplay();
            post_edges();
            update_menus();
            break;
        case 't':
            print_timers(stderr);
            break;
        case 'T':
            timer_report = !timer_report;
            timer_frames = 0;
            break;
        case 'r':
        case 'R':
            init_timers();
            break;
        case '1':
            scene = SCENE_SPHERE_AND_CUBE_1;
            init_timers();
            post_redisplay();
            post_edges();
            update_menus();
            break;
        case '2':
            scene = SCENE_SPHERE_AND_CUBE_2;
            init_timers();
            post_redisplay();
            post_edges();
            update_menus();
            break;
        case '3':
            scene = SCENE_CUBES;
            init_timers();
            post_redisplay();
            post_edges();
            update_menus();
            break;
        case '4':
            scene = SCENE_GRID;
            init_timers();
            post_redisplay();
            post_edges();
            update_menus();
            break;
        case '5':
            scene = SCENE_SPHERES;
            init_timers();
            post_redisplay();
            post_edges();
            update_menus();
            break;
        case '6':
            scene = SCENE_CONE_AND_TORUS;
            init_timers();
            post_redisplay();
            post_edges();
            update_menus();
            break;
        case '7':
            scene = SCENE_TEAPOT;
            init_timers();
            post_redisplay();
            post_edges();
            update_menus();
            break;
        case '8':
            scene = SCENE_TEAPOTS;
            init_timers();
            post_redisplay();
            post_edges();
            update_menus();
            break;
        case 'o':
        case 'O':
            display_mode = DISPLAY_SCENE;
            init_timers();
            post_redisplay();
            post_edges();
            update_menus();
            break;
        case 'p':
        case 'P':
            display_mode = DISPLAY_OUTLINED_SHADOW_VOLUME;
            init_timers();
            post_redisplay();
            post_edges();
            update_menus();
            break;
        case '[':
            display_mode = DISPLAY_DARK_SHADOWS;
            init_timers();
            ambient[0] = 0.0;
            ambient[1] = 0.0;
            ambient[2] = 0.0;
            ambient_mode = 0;
            post_redisplay();
            post_edges();
            update_menus();
            break;
        case ']':
            display_mode = DISPLAY_AMBIENT_SHADOWS;
            init_timers();
            ambient[0] = 0.2;
            ambient[1] = 0.2;
            ambient[2] = 0.2;
            ambient_mode = 1;
            post_redisplay();
            post_edges();
            update_menus();
            break;
        case '\\':
            display_mode = DISPLAY_COMPOSITED_SHADOWS;
            init_timers();
            post_redisplay();
            post_edges();
            update_menus();
            break;
        default:
            break;
    }
    update_window_title();
}

static GLdouble vx, vy, vz;
int motion_active = 0;

/* Mouse button events
 */
void
mouse(
    int button, 
    int state, 
    int x, 
    int y
) {
    /* Flip coordinate systems
     */
    y = yr - y;

    /* Change motion mode
     */
    if (GLUT_LEFT_BUTTON == button) {
        if (state == GLUT_DOWN) {
            motion_active = 1;
            tb_gen_vec(xr,yr,x,y,&vx,&vy,&vz);
            old_demo_mode = demo_mode;
            demo_mode = 0;
        } else {
            motion_active = 0;
            demo_mode = old_demo_mode;
            demo(0);
        }
    }
}

/* Mouse motion events
 */
void
motion(
    int x, 
    int y
) {
    GLdouble rx, ry, rz;
    Mat4x4 iL, iS;

    if (motion_active) {

        /* Flip coordinate systems
         */
        y = yr - y;

        /* Do appropriate state manipulation 
         */
        tb_inc_vec(xr,yr,x,y,&vx,&vy,&vz,&rx,&ry,&rz);
        glMatrixMode(GL_MODELVIEW);
        switch (view_mode) {
            case VIEW_TUMBLE_SCENE:
                glPushMatrix();
                    glLoadIdentity();
                    tb_rotate(rx,ry,rz,0.1);
                    glMultMatrixd(eye_xform);
                    glGetDoublev(GL_MODELVIEW_MATRIX,eye_xform);
                glPopMatrix();
                glPushMatrix();
                    glLoadIdentity();
                    glMultMatrixd(light_xform);
                    tb_rotate(-rx,-ry,-rz,0.1);
                    glGetDoublev(GL_MODELVIEW_MATRIX,light_xform);
                glPopMatrix();
                post_edges();
                break;
            case VIEW_TUMBLE_LIGHT:
                glPushMatrix();
                    glLoadIdentity();
                    tb_rotate(rx,ry,rz,0.1);
                    glMultMatrixd(eye_xform);
                    glMultMatrixd(light_xform);
                    {   Mat4x4 A, B;
                        glGetDoublev(GL_MODELVIEW_MATRIX,B);
                        invert(A,eye_xform);
                        mult(light_xform,A,B);
                    }
                glPopMatrix();
                post_edges();
                break;
            case VIEW_TUMBLE_SCENE_AND_LIGHT:
                glPushMatrix();
                    glLoadIdentity();
                    tb_rotate(rx,ry,rz,0.1);
                    glMultMatrixd(eye_xform);
                    glGetDoublev(GL_MODELVIEW_MATRIX,eye_xform);
                glPopMatrix();
                break;
            default:
                break;
        }

        /* Update rendering in both windows
         */
        update_window_title();
        post_redisplay();
    }
    check_error("motion");
}

static GLdouble lvx, lvy, lvz;
int light_motion_active = 0;

/* Mouse button events in light view
 */
void
mouse_light(
    int button, 
    int state, 
    int x, 
    int y
) {
    /* Flip coordinate systems
     */
    y = lyr - y;

    /* Change motion mode
     */
    if (GLUT_LEFT_BUTTON == button) {
        if (state == GLUT_DOWN) {
            light_motion_active = 1;
            tb_gen_vec(lxr,lyr,x,y,&lvx,&lvy,&lvz);
            old_demo_mode = demo_mode;
            demo_mode = 0;
        } else {
            light_motion_active = 0;
            demo_mode = old_demo_mode;
            demo(0);
        }
    }
}

/* Mouse motion events
 */
void
motion_light(
    int x, 
    int y
) {
    GLdouble rx, ry, rz;

    if (light_motion_active) {

        /* Flip coordinate systems
         */
        y = lyr - y;

        /* Do appropriate state manipulation 
         */
        tb_inc_vec(lxr,lyr,x,y,&lvx,&lvy,&lvz,&rx,&ry,&rz);
        glMatrixMode(GL_MODELVIEW);
        glPushMatrix();
            glLoadMatrixd(light_xform);
            tb_rotate(-rx,-ry,-rz,0.1);
            glGetDoublev(GL_MODELVIEW_MATRIX,light_xform);
        glPopMatrix();

        /* Update rendering in both windows
         */
        update_window_title();
        post_redisplay();
        post_edges();

        check_error("motion_light");
    }
}

/************************************************************************
 * Display Callback
 * 
 * This is really the most important part of this program.   It provides
 * an example of how the CGLS library can be used.
 ************************************************************************/

/* Redisplay event handler for main window.
 */
void
display_scene() {
    GLdouble ze = 0.0001;
    Mat4x4 M, Pt, P;
    GLdouble p[4];
    int xxr, yyr;

    timer_start(&timer_frame);

    /* Get edge map and generate shadow mesh if necessary
     */
    switch (display_mode) {
        case DISPLAY_SCENE:
            break;
        case DISPLAY_SHADOW_VOLUME:
        case DISPLAY_OUTLINED_SHADOW_VOLUME:
        case DISPLAY_SHADOW_VOLUME_POINTS: 
        case DISPLAY_SHADOW_VOLUME_MESH:
        case DISPLAY_SHADOW_VOLUME_NS:
        case DISPLAY_OUTLINED_SHADOW_VOLUME_NS:
        case DISPLAY_SHADOW_VOLUME_POINTS_NS: 
        case DISPLAY_SHADOW_VOLUME_MESH_NS:
        case DISPLAY_DARK_SHADOWS:
        case DISPLAY_AMBIENT_SHADOWS:
        case DISPLAY_COMPOSITED_SHADOWS:
            if (light_window_active) {
                xxr = lxr;
                yyr = lyr;
                acquire_map(light_window,lxr,lyr,main_window,
                            main_scene_base + SCENE_SHADOW_VOLUME);
            } else {
                xxr = xr;
                yyr = yr;
                acquire_map(main_window,xr,yr,main_window,
                            main_scene_base + SCENE_SHADOW_VOLUME);
            }
            break;
    }
    glutSetWindow(main_window);


    /* Set up perspective projection for camera
     */
    glViewport(0, 0, xr, yr);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluPerspective(eye_fov, (GLfloat)xr / (GLfloat)yr, eye_near, eye_far);

    /* Load camera position and orientation.
     */
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    gluLookAt(
        eye_position[0], eye_position[1], eye_position[2],  
        eye_lookat[0], eye_lookat[1], eye_lookat[2],       
        0, 1, 0    /* up */
    );

    /* Set up scene transformation
     */
    glMultMatrixd(eye_xform);

    /* Set up lighting
     */
    glEnable(GL_LIGHTING);
    glPushMatrix();
        glMultMatrixd(light_xform);

	{
	  GLfloat light_direction[] = { 0.0, 0.0, 0.0, 1.0 };
	  light_direction[0] = -light_position[0];
	  light_direction[1] = -light_position[1];
	  light_direction[2] = -light_position[2];

	  switch (light_mode) {
	  case LIGHT_DIRECTION:
	    light_position[3] = 0.0;
	    glLightf(GL_LIGHT0, GL_SPOT_EXPONENT, 0);
	    glLightf(GL_LIGHT0, GL_SPOT_CUTOFF, 180);
	    break;
	  case LIGHT_SPOT:
	    light_position[3] = 1.0;
	    glLightf(GL_LIGHT0, GL_SPOT_EXPONENT, 60);
	    glLightf(GL_LIGHT0, GL_SPOT_CUTOFF, 12);
	    break;
	  case LIGHT_POINT:
	    light_position[3] = 1.0;
	    glLightf(GL_LIGHT0, GL_SPOT_EXPONENT, 0);
	    glLightf(GL_LIGHT0, GL_SPOT_CUTOFF, 180);
	    break;
	  }

	  glLightfv(GL_LIGHT0, GL_SPOT_DIRECTION, light_direction);
	  glLightfv(GL_LIGHT0, GL_POSITION, light_position);
	  glLightfv(GL_LIGHT0, GL_DIFFUSE, light_power);
	  glLightfv(GL_LIGHT0, GL_SPECULAR, light_power);
	  glLightfv(GL_LIGHT0, GL_AMBIENT, black);
	  glEnable(GL_LIGHT0);
	}

        glLightModelfv(GL_LIGHT_MODEL_AMBIENT, ambient);
        glLightModeli(GL_LIGHT_MODEL_LOCAL_VIEWER, 1);
        glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, 1);
    glPopMatrix();

    /* Set correct depth/stencil modes
     */
    glEnable(GL_DEPTH_TEST);
    glDepthMask(GL_TRUE);
    glDepthFunc(GL_LEQUAL);
    glStencilFunc(GL_ALWAYS,0,0);

    /* Clear buffers and render scene
     */
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT | GL_STENCIL_BUFFER_BIT);

    glCullFace(GL_BACK);
    glEnable(GL_CULL_FACE);
    switch (display_mode) {
        case DISPLAY_SHADOW_VOLUME_NS:
        case DISPLAY_OUTLINED_SHADOW_VOLUME_NS:
        case DISPLAY_SHADOW_VOLUME_POINTS_NS: 
        case DISPLAY_SHADOW_VOLUME_MESH_NS:
            break;
        default:
            timer_start(&timer_render_scene);
            render_scene(main_scene_base);
            timer_stop(&timer_render_scene);
            break;
    }

    /* Render stencil if called for
     */
    switch (display_mode) {
        case DISPLAY_DARK_SHADOWS:
        case DISPLAY_AMBIENT_SHADOWS:
        case DISPLAY_COMPOSITED_SHADOWS:
            timer_start(&timer_render_volume);

            /* Set up stencil buffer to invert 
             */
            glDepthMask(GL_FALSE);
            glEnable(GL_DEPTH_TEST);

            glDisable(GL_LIGHTING);
            glColorMask(GL_FALSE, GL_FALSE, GL_FALSE, GL_FALSE);

            glEnable(GL_STENCIL_TEST);
            glStencilFunc(GL_ALWAYS, 1, 1);
            glStencilOp(GL_KEEP, GL_KEEP, GL_INVERT);

            /* Render shadow volume into stencil buffer
             */
            glMatrixMode(GL_PROJECTION);
            glPushMatrix();

                /* Shift polygons away from eye
                 */
                glGetDoublev(GL_PROJECTION_MATRIX,M);
                glLoadIdentity();
                glTranslated(0,0,eye_ze);
                glMultMatrixd(M);

                /* Render both front and back faces
                 */
                glDisable(GL_CULL_FACE);
                glCallList(main_scene_base + SCENE_SHADOW_VOLUME);

                timer_stop(&timer_render_volume);

                /* Render caps to correct parity
                 */
                glDisable(GL_DEPTH_TEST);
                cap(L,xxr,yyr);

            glMatrixMode(GL_PROJECTION);
            glPopMatrix();

            break;
        default:
            break;
    }

    /* Render shadows or shadow volume
     */
    switch (display_mode) {
        case DISPLAY_SCENE:
            break;
        case DISPLAY_DARK_SHADOWS:
            timer_start(&timer_render_shadow);

            /* Render black polygon over scene
             */
            glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE);
            glStencilFunc(GL_EQUAL, 1, 1);
            glStencilOp(GL_KEEP, GL_KEEP, GL_KEEP);
            glDisable(GL_LIGHTING);
            glColor3fv(black);

            glMatrixMode(GL_MODELVIEW);
            glPushMatrix();
                glLoadIdentity();
                glMatrixMode(GL_PROJECTION);
                glPushMatrix();
                    glLoadIdentity();
                    glOrtho(0,1,0,1,-1.0,1.0);
                    glBegin(GL_QUADS);
                        glVertex2i(0,0);
                        glVertex2i(1,0);
                        glVertex2i(1,1);
                        glVertex2i(0,1);
                    glEnd();
                glMatrixMode(GL_PROJECTION);
                glPopMatrix();
            glMatrixMode(GL_MODELVIEW);
            glPopMatrix();

            timer_stop(&timer_render_shadow);
            break;
        case DISPLAY_COMPOSITED_SHADOWS:
            timer_start(&timer_render_shadow);

            /* Composite black polygon into scene
             */
            glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE);
            glStencilFunc(GL_EQUAL, 1, 1);
            glStencilOp(GL_KEEP, GL_KEEP, GL_KEEP);
            glDisable(GL_LIGHTING);
            glColor4fv(transparent_black);

            glMatrixMode(GL_MODELVIEW);
            glPushMatrix();
                glLoadIdentity();
                glMatrixMode(GL_PROJECTION);
                glPushMatrix();
                    glLoadIdentity();
                    glOrtho(0,1,0,1,-1.0,1.0);
                    glEnable(GL_BLEND);
                    glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);
                    glBegin(GL_QUADS);
                        glVertex2i(0,0);
                        glVertex2i(1,0);
                        glVertex2i(1,1);
                        glVertex2i(0,1);
                    glEnd();
                    glDisable(GL_BLEND);
                glMatrixMode(GL_PROJECTION);
                glPopMatrix();
            glMatrixMode(GL_MODELVIEW);
            glPopMatrix();

            timer_stop(&timer_render_shadow);
            break;
        case DISPLAY_AMBIENT_SHADOWS:
            timer_start(&timer_render_shadow);

            /* Rerender scene with stencil and only ambient illumination
             */
            glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE);
            glDepthMask(GL_TRUE);
            glDepthFunc(GL_LEQUAL);
            glEnable(GL_DEPTH_TEST);
            glEnable(GL_LIGHTING);
            glDepthFunc(GL_LEQUAL);
            glStencilFunc(GL_EQUAL, 1, 1);
            glStencilOp(GL_KEEP, GL_KEEP, GL_KEEP);

            /* Change to ambient-only illumination
             */
            glLightfv(GL_LIGHT0, GL_DIFFUSE, black);
            glLightfv(GL_LIGHT0, GL_SPECULAR, black);

            /* Render scene
             */
            glMatrixMode(GL_MODELVIEW);
            glCullFace(GL_BACK);
            glEnable(GL_CULL_FACE);
            render_scene(main_scene_base);

            timer_stop(&timer_render_shadow);
            break;
        case DISPLAY_SHADOW_VOLUME:
        case DISPLAY_SHADOW_VOLUME_NS:
            timer_start(&timer_render_shadow);

            glDisable(GL_LIGHTING);
            glDepthMask(GL_TRUE);
            glEnable(GL_DEPTH_TEST);
            glDepthFunc(GL_LESS);
            glEnable(GL_CULL_FACE);

            glMatrixMode(GL_PROJECTION);
            glPushMatrix();

                /* Shift polygons away from eye
                 */
                glGetDoublev(GL_PROJECTION_MATRIX,M);
                glLoadIdentity();
                glTranslated(0,0,eye_ze);
                glMultMatrixd(M);

                /* Front faces in yellow
                 */
                glColor3d(1.0,1.0,0.0);
                glCullFace(GL_BACK);
                glCallList(main_scene_base + SCENE_SHADOW_VOLUME);

                /* Back faces in blue
                 */
                glColor3d(0.0,0.0,1.0);
                glCullFace(GL_FRONT);
                glCallList(main_scene_base + SCENE_SHADOW_VOLUME);

                /* Caps in semitransparent red
                 */
                glEnable(GL_BLEND);
                glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);
                glCullFace(GL_BACK);
                glDisable(GL_CULL_FACE);
                glDisable(GL_DEPTH_TEST);
                glDepthMask(GL_FALSE);
                glColor4d(1.0,0.0,0.0,0.5);
                cap(L,xxr,yyr);
                glDisable(GL_BLEND);

            glMatrixMode(GL_PROJECTION);
            glPopMatrix();

            timer_stop(&timer_render_shadow);
            break;
        case DISPLAY_OUTLINED_SHADOW_VOLUME:
        case DISPLAY_OUTLINED_SHADOW_VOLUME_NS:
            timer_start(&timer_render_shadow);

            glDisable(GL_LIGHTING);
            glDepthMask(GL_TRUE);
            glEnable(GL_DEPTH_TEST);
            glDepthFunc(GL_LESS);
            glEnable(GL_CULL_FACE);

            glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);

            glMatrixMode(GL_PROJECTION);
            glPushMatrix();

                /* Shift polygons away from eye
                 */
                glGetDoublev(GL_PROJECTION_MATRIX,M);
                glLoadIdentity();
                glTranslated(0,0,line_ze);
                glMultMatrixd(M);

                /* Front faces in yellow
                 */
                glColor3d(1.0,1.0,0.0);
                glCullFace(GL_BACK);
                glCallList(main_scene_base + SCENE_SHADOW_VOLUME);

                /* Back faces in blue
                 */
                glColor3d(0.0,0.0,1.0);
                glCullFace(GL_FRONT);
                glCallList(main_scene_base + SCENE_SHADOW_VOLUME);

                /* Caps in semitransparent red
                 */
                glEnable(GL_BLEND);
                glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);
                glDisable(GL_CULL_FACE);
                glDisable(GL_DEPTH_TEST);
                glDepthMask(GL_FALSE);
                glColor4d(1.0,0.0,0.0,0.5);
                cap(L,xxr,yyr);
                glDisable(GL_BLEND);

            glMatrixMode(GL_PROJECTION);
            glPopMatrix();

            glEnable(GL_DEPTH_TEST);
            glEnable(GL_CULL_FACE);
            glDepthMask(GL_TRUE);

            glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);

            glMatrixMode(GL_PROJECTION);
            glPushMatrix();

                /* Shift polygons away from eye
                 */
                glGetDoublev(GL_PROJECTION_MATRIX,M);
                glLoadIdentity();
                glTranslated(0,0,eye_ze*0.9);
                glMultMatrixd(M);

                /* Front faces with red edging
                 */
                glColor3d(1.0,0.0,0.0);
                glCullFace(GL_BACK);
                glCallList(main_scene_base + SCENE_SHADOW_VOLUME);

                /* Back faces with green edging
                 */
                glColor3d(0.0,1.0,0.0);
                glCullFace(GL_FRONT);
                glCallList(main_scene_base + SCENE_SHADOW_VOLUME);

                /* Caps with opaque blue edging
                 */
                glCullFace(GL_BACK);
                glDisable(GL_CULL_FACE);
                glDisable(GL_DEPTH_TEST);
                glColor3d(0.0,0.0,1.0);
                cap(L,xxr,yyr);

            glMatrixMode(GL_PROJECTION);
            glPopMatrix();

            glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);

            timer_stop(&timer_render_shadow);
            break;
        case DISPLAY_SHADOW_VOLUME_POINTS:
        case DISPLAY_SHADOW_VOLUME_POINTS_NS:
            timer_start(&timer_render_shadow);

            glDisable(GL_LIGHTING);
            glDepthMask(GL_TRUE);
            glEnable(GL_DEPTH_TEST);
            glDepthFunc(GL_LEQUAL);
            glEnable(GL_CULL_FACE);

            glPolygonMode(GL_FRONT_AND_BACK,GL_POINT);

            glMatrixMode(GL_PROJECTION);
            glPushMatrix();

                /* Back faces with green edging
                 */
                glColor3d(0.0,1.0,0.0);
                glCullFace(GL_FRONT);
                glCallList(main_scene_base + SCENE_SHADOW_VOLUME);

                /* Front faces with red edging
                 */
                glColor3d(1.0,0.0,0.0);
                glCullFace(GL_BACK);
                glCallList(main_scene_base + SCENE_SHADOW_VOLUME);

                /* Caps with opaque blue edging
                 */
                glCullFace(GL_BACK);
                glDisable(GL_CULL_FACE);
                glDisable(GL_DEPTH_TEST);
                glColor3d(0.0,0.0,1.0);
                cap(L,xxr,yyr);

            glMatrixMode(GL_PROJECTION);
            glPopMatrix();

            glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);

            timer_stop(&timer_render_shadow);
            break;
        case DISPLAY_SHADOW_VOLUME_MESH:
        case DISPLAY_SHADOW_VOLUME_MESH_NS:
            timer_start(&timer_render_shadow);

            glDisable(GL_LIGHTING);
            glDepthMask(GL_TRUE);
            glEnable(GL_DEPTH_TEST);
            glDepthFunc(GL_LEQUAL);
            glEnable(GL_CULL_FACE);

            glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);

            glMatrixMode(GL_PROJECTION);
            glPushMatrix();

                /* Back faces with green edging
                 */
                glColor3d(0.0,1.0,0.0);
                glCullFace(GL_FRONT);
                glCallList(main_scene_base + SCENE_SHADOW_VOLUME);

                /* Front faces with red edging
                 */
                glColor3d(1.0,0.0,0.0);
                glCullFace(GL_BACK);
                glCallList(main_scene_base + SCENE_SHADOW_VOLUME);

                /* Caps with opaque blue edging
                 */
                glCullFace(GL_BACK);
                glDisable(GL_CULL_FACE);
                glDisable(GL_DEPTH_TEST);
                glColor3d(0.0,0.0,1.0);
                cap(L,xxr,yyr);

            glMatrixMode(GL_PROJECTION);
            glPopMatrix();

            glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);

            timer_stop(&timer_render_shadow);
            break;
    }

    /* Show rendering.
     */
    glutSwapBuffers();

    timer_stop(&timer_frame);

    /* Trigger next frame of demo here, to avoid overrunning
     * rendering latency
     */
    if (demo_mode) {
        glutTimerFunc(demo_delay,demo,0);
    }

    check_error("display_scene");
    if (timer_report && (timer_frames % 10 == 0)) {
        print_timers(stderr);
    }
    timer_frames++;
}

/* Redisplay event handler for light view window.
 * Note that we don't update the edge detection stuff here,
 * which may result in some extra work.   However, if the light
 * view is active, it will be used to render the shadow map,
 * so the resolution can be controlled.    
 */
void
display_lightview() {
    int edges_updated;

    configure_lightview(lxr,lyr);

    /* Set up rendering mode
     */
    switch (light_window_mode) {
        case LW_SHADED:
            glEnable(GL_LIGHTING);
            glDisable(GL_FOG);

            glMatrixMode(GL_MODELVIEW);
            glPushMatrix();

                glLoadIdentity();
                glLightfv(GL_LIGHT0, GL_POSITION, light_position);
                glLightfv(GL_LIGHT0, GL_DIFFUSE, light_power);
                glLightfv(GL_LIGHT0, GL_SPECULAR, light_power);
                glLightfv(GL_LIGHT0, GL_AMBIENT, black);
                glEnable(GL_LIGHT0);

                glLightModelfv(GL_LIGHT_MODEL_AMBIENT, ambient);
                glLightModeli(GL_LIGHT_MODEL_LOCAL_VIEWER, 1);
                glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, 1);

            glPopMatrix();

            glEnable(GL_DEPTH_TEST);
            glDepthMask(GL_TRUE);
            glStencilFunc(GL_ALWAYS,0,0);

            /* Clear buffers and render scene
             */
            glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
            glMatrixMode(GL_MODELVIEW);
            glCullFace(GL_BACK);
            glEnable(GL_CULL_FACE);

            render_scene(light_scene_base);

            break;

        case LW_DEPTH:
            glDisable(GL_LIGHTING);
            glColor3d(1.0,1.0,1.0);
            glEnable(GL_FOG);
            glFogf(GL_FOG_MODE,GL_LINEAR);
            glFogf(GL_FOG_START,light_near);
            glFogf(GL_FOG_END,light_far);
            glFogfv(GL_FOG_COLOR,black);

            glEnable(GL_DEPTH_TEST);
            glDepthMask(GL_TRUE);
            glStencilFunc(GL_ALWAYS,0,0);

            /* Clear buffers and render scene
             */
            glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
            glMatrixMode(GL_MODELVIEW);
            glCullFace(GL_BACK);
            glEnable(GL_CULL_FACE);

            render_scene(light_scene_base);

            break;

        case LW_EDGE:

            glDisable(GL_LIGHTING);
            glDisable(GL_FOG);

            /* Update edge map (also builds shadow volume
             * in main window, which is a waste of time if not
             * rendering shadows, but most of the time we will be,
             * except in stupid demos like this)
             */
            acquire_map(light_window,lxr,lyr,main_window,
                        main_scene_base + SCENE_SHADOW_VOLUME);
            glutSetWindow(light_window);

            /* Draw edge map on image
             */
            glMatrixMode(GL_MODELVIEW);
            glPushMatrix();

                glLoadIdentity();

                glMatrixMode(GL_PROJECTION);
                glPushMatrix();

                    glLoadIdentity();
                    glOrtho(0,lxr,0,lyr,-1.0,1.0);
                    glMatrixMode(GL_MODELVIEW);

                    glColorMask(GL_TRUE,GL_TRUE,GL_TRUE,GL_TRUE);
                    glClear(GL_COLOR_BUFFER_BIT);

                    {
                        int i;
                        GLfloat m[256];

                        for (i=0; i<256; i++) {
                            if (EDGE_MASK & i) {
                                    m[i] = 1.0;
                            } else {
                                    m[i] = 0.0;
                            }
                        }

                        glPixelMapfv(GL_PIXEL_MAP_I_TO_R,256,m);
                        glPixelMapfv(GL_PIXEL_MAP_I_TO_B,256,m);

                        for (i=16; i<256; i++) {
                            m[i] = 1.0;
                        }
                        glPixelMapfv(GL_PIXEL_MAP_I_TO_G,256,m);

                        for (i=0; i<16; i++) {
                            m[i] = 1.0;
                        }
                        glPixelMapfv(GL_PIXEL_MAP_I_TO_A,256,m);

                    }
                    glPixelStorei(GL_UNPACK_ALIGNMENT,sizeof(map_code_type));
                    glRasterPos2i(0,0);
                    glDrawPixels(map_xr,map_yr,
                                 GL_COLOR_INDEX,MAP_CODE_TYPE,map_code[0]);

                glMatrixMode(GL_PROJECTION);
                glPopMatrix();

            glMatrixMode(GL_MODELVIEW);
            glPopMatrix();

            break;

        case LW_HW_EDGE:

            glDepthMask(GL_TRUE);
            glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE);
            glStencilMask(~0);
            glClearStencil(0);

            /* Clear buffers
             */
            glClear(GL_COLOR_BUFFER_BIT | 
                    GL_DEPTH_BUFFER_BIT |
                    GL_STENCIL_BUFFER_BIT 
            );

            /* Render depth (also show depth-cued image in background)
             */
            glDisable(GL_LIGHTING);

            glColor3d(1.0, 0.0, 0.0);
            glEnable(GL_FOG);
            glFogf(GL_FOG_MODE,GL_LINEAR);
            glFogf(GL_FOG_START,light_near);
            glFogf(GL_FOG_END,light_far);
            glFogfv(GL_FOG_COLOR,black);
            glCullFace(GL_BACK);
            glEnable(GL_CULL_FACE);
            glDepthFunc(GL_LESS);
            glEnable(GL_DEPTH_TEST);
            glDisable(GL_STENCIL_TEST);
            glStencilMask(0);
            glStencilOp(GL_ZERO, GL_ZERO, GL_ZERO);

            glMatrixMode(GL_MODELVIEW);
            render_scene(light_scene_base);

            glDisable(GL_FOG);
            glDisable(GL_CULL_FACE);
            glColorMask(GL_FALSE, GL_FALSE, GL_FALSE, GL_FALSE);

            /* Copy shifted and biased versions of the depth
             * buffer to do edge detection (results in colour and
             * stencil buffers).
             */
            glMatrixMode(GL_MODELVIEW);
            glPushMatrix();
                glLoadIdentity();
                glMatrixMode(GL_PROJECTION);
                    glPushMatrix();
                    glLoadIdentity();
                    glOrtho(0,lxr, 0,lyr, -1,1);
                    glMatrixMode(GL_MODELVIEW);

                    glDepthMask(GL_FALSE);
                    glEnable(GL_DEPTH_TEST);

                    glEnable(GL_STENCIL_TEST);
                    glStencilFunc(GL_ALWAYS,3,3);
                    glStencilOp(GL_KEEP, GL_KEEP, GL_REPLACE);

                /*
                    glPixelTransferf(GL_DEPTH_SCALE, 1-1e-2);
                    glPixelTransferf(GL_DEPTH_BIAS,  1e-2);
                */
                    glPixelTransferf(GL_DEPTH_BIAS,  
                        (GLfloat)th_L/(GLfloat)MAX_MAP_Z);
                    glDepthFunc(GL_LESS);

                    glStencilMask(1);
                    glRasterPos2i(1,0);
                    glCopyPixels(0,0,lxr-2,lyr-1,GL_DEPTH);

                    glStencilMask(2);
                    glRasterPos2i(0,1);
                    glCopyPixels(0,0,lxr-1,lyr-2,GL_DEPTH);

                /*
                    glPixelTransferf(GL_DEPTH_SCALE, 1+1e-2);
                    glPixelTransferf(GL_DEPTH_BIAS,  -1e-2);
                */
                    glPixelTransferf(GL_DEPTH_BIAS,  
                        -(GLfloat)th_L/(GLfloat)MAX_MAP_Z);
                    glDepthFunc(GL_GREATER);

                    glStencilMask(1);
                    glRasterPos2i(1,0);
                    glCopyPixels(0,0,lxr-2,lyr-1,GL_DEPTH);

                    glStencilMask(2);
                    glRasterPos2i(0,1);
                    glCopyPixels(0,0,lxr-1,lyr-2,GL_DEPTH);

                    /* Render polygon over scene where stencil set
                     */
                    glStencilFunc(GL_NOTEQUAL,0,3);
                    glStencilOp(GL_KEEP,GL_KEEP,GL_KEEP);

                    glDisable(GL_DEPTH_TEST);

                    glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE);
                    glColor3d(0.0, 1.0, 0.0);

                    glBegin(GL_QUADS);
                        glVertex2i(0,0);
                        glVertex2i(lxr,0);
                        glVertex2i(lxr,lyr);
                        glVertex2i(0,lyr);
                    glEnd();

                glMatrixMode(GL_PROJECTION);
                glPopMatrix();
            glMatrixMode(GL_MODELVIEW);
            glPopMatrix();

            glPixelTransferi(GL_DEPTH_SCALE,1);
            glPixelTransferi(GL_DEPTH_BIAS,0);
            glDisable(GL_STENCIL_TEST);
            glDepthFunc(GL_LESS);
            glColorMask(GL_TRUE,GL_TRUE,GL_TRUE,GL_TRUE);
            break;
    }

    glutSwapBuffers();
    check_error("display_lightview");
}

/************************************************************************
 * Menu Update
 *
 * Changes the menu entries to reflect the current status of various
 * parameters.
 ************************************************************************/

void
update_menus() {
    char buf[1000];
    if (menu_state == GLUT_MENU_IN_USE) return;
    if (!menu_change) return;

    glutSetWindow(main_window);
    glutSetMenu(main_menu_id);

    sprintf(buf,"SCENE: %s", scene_name[scene]);
    glutChangeToSubMenu(M_SCENE, buf, scene_menu_id);

    sprintf(buf,"DISPLAY: %s", display_mode_name[display_mode]);
    glutChangeToSubMenu(M_DISPLAY, buf, display_menu_id);

    sprintf(buf,"EDGE: %s", edge_mode_name[edge_mode]);
    glutChangeToSubMenu(M_EDGE, buf, edge_menu_id);

    sprintf(buf,"RECON: %s", recon_mode_name[recon_mode]);
    glutChangeToSubMenu(M_RECON, buf, recon_menu_id);

    sprintf(buf,"CAP: %s", cap_mode_name[cap_mode]);
    glutChangeToSubMenu(M_CAP, buf, cap_menu_id);

    switch (ambient_mode) {
        case 0:
            glutChangeToSubMenu(M_AMBIENT, "AMBIENT: Off", 
                                ambient_menu_id);
            break;
        case 1:
            glutChangeToSubMenu(M_AMBIENT, "AMBIENT: On", 
                                ambient_menu_id);
            break;
        case 2:
            glutChangeToSubMenu(M_AMBIENT, "AMBIENT: Bright", 
                                ambient_menu_id);
            break;
    }
    switch (subd) {
        case SUBD_LOW:
            glutChangeToSubMenu(M_SUBD, "SUBD: Low", 
                                subd_menu_id);
            break;
        case SUBD_MEDIUM:
            glutChangeToSubMenu(M_SUBD, "SUBD: Medium", 
                                subd_menu_id);
            break;
        case SUBD_HIGH:
            glutChangeToSubMenu(M_SUBD, "SUBD: High", 
                                subd_menu_id);
            break;
    }
    switch (view_mode) {
        case VIEW_TUMBLE_SCENE:
            glutChangeToSubMenu(M_VIEW, "VIEW: Tumble Scene", 
                                view_menu_id);
            break;
        case VIEW_TUMBLE_LIGHT:
            glutChangeToSubMenu(M_VIEW, "VIEW: Tumble Light", 
                                view_menu_id);
            break;
        case VIEW_TUMBLE_SCENE_AND_LIGHT:
            glutChangeToSubMenu(M_VIEW, "VIEW: Tumble Scene and Light", 
                                view_menu_id);
            break;
    }
    switch (light_mode) {
        case LIGHT_POINT:
            glutChangeToSubMenu(M_LIGHT, "LIGHT: Point Source", 
                                light_menu_id);
            break;
        case LIGHT_DIRECTION:
            glutChangeToSubMenu(M_LIGHT, "LIGHT: Directional Source", 
                                light_menu_id);
            break;
        case LIGHT_SPOT:
            glutChangeToSubMenu(M_LIGHT, "LIGHT: Spotlight Source", 
                                light_menu_id);
            break;
    }
    switch (demo_mode) {
        case 0:
            glutChangeToSubMenu(M_DEMO, "DEMO: Off", 
                                demo_menu_id);
            break;
        case 1:
            glutChangeToSubMenu(M_DEMO, "DEMO: On", 
                                demo_menu_id);
            break;
    }
    menu_change = 0;
}

void
menu_state_func(
    int status
) {
    menu_state = status;
}

/************************************************************************
 * Menu Callbacks
 ************************************************************************/

/* Scene selection sub-menu.
 */
void
scene_menu(
    int selection
) {
    menu_change = 1;
    scene = selection;
    post_redisplay();
    post_edges();
    update_menus();
    init_timers();
}

/* Display mode sub-menu.
 */
void
display_menu(
    int selection
) {
    menu_change = 1;
    display_mode = selection;
    post_redisplay();
    update_menus();
    init_timers();
}

/* Edge detection operator sub-menu.
 */
void
edge_menu(
    int selection
) {
    menu_change = 1;
    edge_mode = selection;
    post_redisplay();
    post_edges();
    update_menus();
    init_timers();
}

/* Reconstruction technique sub-menu.
 */
void
recon_menu(
    int selection
) {
    menu_change = 1;
    recon_mode = selection;
    post_redisplay();
    post_edges();
    update_menus();
    init_timers();
}

/* Capping mode sub-menu.
 */
void
cap_menu(
    int selection
) {
    menu_change = 1;
    cap_mode = selection;
    post_redisplay();
    post_edges();
    update_menus();
    init_timers();
}

/* Ambient on/off sub-menu.
 */
void
ambient_menu(
    int selection
) {
    switch (selection) {
        case 0:
            menu_change = 1;
            ambient[0] = 0.0;
            ambient[1] = 0.0;
            ambient[2] = 0.0;
            ambient_mode = 0;
            glutPostRedisplay();
            break;
        case 1:
            menu_change = 1;
            ambient[0] = 0.2;
            ambient[1] = 0.2;
            ambient[2] = 0.2;
            ambient_mode = 1;
            glutPostRedisplay();
            break;
        case 2:
            menu_change = 1;
            ambient[0] = 0.5;
            ambient[1] = 0.5;
            ambient[2] = 0.5;
            ambient_mode = 2;
            glutPostRedisplay();
            break;
        default:
            break;
    }
    update_menus();
}

/* Subdivision sub-menu.
 */
void
subd_menu(
    int selection
) {
    menu_change = 1;
    subd = selection;
    glutSetWindow(light_window);
    init_scenes(&light_scene_base);
    glutSetWindow(main_window);
    init_scenes(&main_scene_base);
    post_redisplay();
    post_edges();
    update_menus();
}

/* Tumble mode sub-menu.
 */
void
view_menu(
    int selection
) {
    menu_change = 1;
    switch (selection) {
        case VIEW_RESET_SCENE:
            init_xforms();
            break;
        default:
            view_mode = selection;
            update_menus();
            break;
    }
    post_redisplay();
}

/* Lighting model sub-menu.
 */
void
light_menu(
    int selection
) {
    menu_change = 1;
    light_mode = selection;
    post_redisplay();
    update_menus();
}

/* Demo on/off sub-menu.
 */
void
demo_menu(
    int selection
) {
    switch (selection) {
        case 0:
            menu_change = 1;
            demo_mode = 0;
            break;
        case 1:
            menu_change = 1;
            demo_mode = 1;
            glutTimerFunc(demo_delay,demo,0);
            break;
        default:
            break;
    }
}

/* Main menu (and status display).
 */
void
main_menu(
    int selection
) {
    switch (selection) {
        case M_LIGHT_WINDOW:
            light_window_active = 1;
            glutSetWindow(light_window);
            glutShowWindow();
            post_redisplay();
            post_edges();
            break;
        case M_QUIT:
            exit(0);
            break;
        default:
            break;
    }
}

/* Light window menu
 */
void
light_window_menu(
    int selection
) {
    switch (selection) {
        case LM_CLOSE:
            light_window_active = 0;
            glutHideWindow();
            post_edges();
            post_redisplay();
            break;
        case LM_QUIT:
            exit(0);
            break;
        default:
            break;
    }
}

/* Light window mode selector
 */
void
light_window_mode_menu(
    int selection
) {
    light_window_mode = selection;
    post_redisplay();
}

/* Initialise menus
 */
void
init_menus() {
    int i;

    glutSetWindow(main_window);

    scene_menu_id = glutCreateMenu(scene_menu);
    for (i=0; i < SCENE_N; i++) glutAddMenuEntry(scene_name[i],i); 

    display_menu_id = glutCreateMenu(display_menu);
    for (i=0; i < DISPLAY_N; i++) glutAddMenuEntry(display_mode_name[i],i);

    edge_menu_id = glutCreateMenu(edge_menu);
    for (i=0; i < EDGE_N; i++) glutAddMenuEntry(edge_mode_name[i],i);

    recon_menu_id = glutCreateMenu(recon_menu);
    for (i=0; i < RECON_N; i++) glutAddMenuEntry(recon_mode_name[i],i);

    cap_menu_id = glutCreateMenu(cap_menu);
    for (i=0; i < CAP_N; i++) glutAddMenuEntry(cap_mode_name[i],i);

    ambient_menu_id = glutCreateMenu(ambient_menu);
    glutAddMenuEntry("Off", 0);
    glutAddMenuEntry("On", 1);
    glutAddMenuEntry("Bright", 2);

    subd_menu_id = glutCreateMenu(subd_menu);
    glutAddMenuEntry("Low",    
        SUBD_LOW);
    glutAddMenuEntry("Medium", 
        SUBD_MEDIUM);
    glutAddMenuEntry("High",   
        SUBD_HIGH);

    view_menu_id = glutCreateMenu(view_menu);
    glutAddMenuEntry("RESET SCENE",    
        VIEW_RESET_SCENE);
    glutAddMenuEntry("Tumble Scene",    
        VIEW_TUMBLE_SCENE);
    glutAddMenuEntry("Tumble Light",             
        VIEW_TUMBLE_LIGHT);
    glutAddMenuEntry("Tumble Scene and Light",   
        VIEW_TUMBLE_SCENE_AND_LIGHT);

    light_menu_id = glutCreateMenu(light_menu);
    glutAddMenuEntry("Point Source",    
        LIGHT_POINT);
    glutAddMenuEntry("Directional Source",             
        LIGHT_DIRECTION);
    glutAddMenuEntry("Spotlight Source",   
        LIGHT_SPOT);

    demo_menu_id = glutCreateMenu(demo_menu);
    glutAddMenuEntry("Off", 0);
    glutAddMenuEntry("On", 1);

    main_menu_id = glutCreateMenu(main_menu);
    glutAddSubMenu("SCENE", scene_menu_id);
    glutAddSubMenu("DISPLAY", display_menu_id);
    glutAddSubMenu("EDGE", edge_menu_id);
    glutAddSubMenu("RECON", recon_menu_id);
    glutAddSubMenu("CAP", cap_menu_id);
    glutAddSubMenu("AMBIENT", ambient_menu_id);
    glutAddSubMenu("SUBD", subd_menu_id);
    glutAddSubMenu("VIEW", view_menu_id);
    glutAddSubMenu("LIGHT", light_menu_id);
    glutAddSubMenu("DEMO", demo_menu_id);
    glutAddMenuEntry("CREATE LIGHT VIEW", M_LIGHT_WINDOW);
    glutAddMenuEntry("QUIT", M_QUIT);

    glutMenuStateFunc(menu_state_func);
    menu_change = 1;
    update_menus();

    glutAttachMenu(GLUT_RIGHT_BUTTON);

    glutSetWindow(light_window);

    light_window_mode_menu_id = glutCreateMenu(light_window_mode_menu);
    glutAddMenuEntry("Shaded", LW_SHADED);
    glutAddMenuEntry("Depth", LW_DEPTH);
    glutAddMenuEntry("Edges", LW_EDGE);
    glutAddMenuEntry("Hardware Edges", LW_HW_EDGE);

    light_window_menu_id = glutCreateMenu(light_window_menu);
    glutAddSubMenu("MODE:", light_window_mode_menu_id);
    glutAddMenuEntry("CLOSE", LM_CLOSE);
    glutAddMenuEntry("QUIT", LM_QUIT);

    glutAttachMenu(GLUT_RIGHT_BUTTON);
}

/************************************************************************
 * MAIN
 ************************************************************************/

int 
main(int argc, char** argv)
{
    /* Identification
     */
    print_sep(stdout,'=');
    fprintf(stdout,"Shadow Volume Reconstruction\n");
    fprintf(stdout,"TEST/DEMO PROGRAM\n");
    print_sep(stdout,'-');
    fprintf(stdout,"Contact:\n");
    fprintf(stdout,"\tMichael McCool\n");
    fprintf(stdout,"\tComputer Graphics Lab\n");
    fprintf(stdout,"\tUniversity of Waterloo\n");
    fprintf(stdout,"\tWaterloo, Ontario N2L 3G1\n");
    fprintf(stdout,"\tCanada\n");
    fprintf(stdout,"\tmmccool@cgl.uwaterloo.ca\n");
    fprintf(stdout,"\thttp://www.cgl.uwaterloo.ca/~mmccool/\n");

    /* Initialize GLUT
     */
    glutInitDisplayMode(GLUT_RGB | GLUT_DEPTH | GLUT_DOUBLE | GLUT_STENCIL);
    glutInitWindowSize(xr,yr);
    glutInit(&argc, argv);

    /* Create main window
     */
    main_window = glutCreateWindow("shadow");

    /* Register callbacks with main window
     */
    glutReshapeFunc(reshape);
    glutDisplayFunc(display_scene);
    glutKeyboardFunc(keyboard);
    glutMotionFunc(motion);
    glutMouseFunc(mouse);

    /* Build scene display lists
     */
    init_scenes(&main_scene_base);

    /* Create light view window
     */
    light_window = glutCreateWindow("shadow: light view");

    /* Register callbacks with light view window
     */
    glutReshapeFunc(reshape_light);
    glutKeyboardFunc(keyboard);
    glutDisplayFunc(display_lightview);
    glutMotionFunc(motion_light);
    glutMouseFunc(mouse_light);

    /* Build scene display lists (in light view window context)
     */
    init_scenes(&light_scene_base);

    /* Hide light view window until needed.
     */
    glutHideWindow();
    light_window_active = 0;

    /* Construct menus
     */
    init_menus();

    /* Set up initial transformations
     */
    init_xforms();

    /* Print out hardware resolutions
     */
    print_sep(stdout,'-');
    fprintf(stdout,"Hardware Configuration:\n");
    {
        GLint i;

        glGetIntegerv(GL_RED_BITS,&i);
        fprintf(stdout,"\tGL_RED_BITS = %d\n",i);
        glGetIntegerv(GL_GREEN_BITS,&i);
        fprintf(stdout,"\tGL_GREEN_BITS = %d\n",i);
        glGetIntegerv(GL_BLUE_BITS,&i);
        fprintf(stdout,"\tGL_BLUE_BITS = %d\n",i);
        glGetIntegerv(GL_ALPHA_BITS,&i);
        fprintf(stdout,"\tGL_ALPHA_BITS = %d\n",i);
        glGetIntegerv(GL_DEPTH_BITS,&i);
        fprintf(stdout,"\tGL_DEPTH_BITS = %d\n",i);
        glGetIntegerv(GL_STENCIL_BITS,&i);
        fprintf(stdout,"\tGL_STENCIL_BITS = %d\n",i);
        glGetIntegerv(GL_MAX_LIGHTS,&i);
        fprintf(stdout,"\tGL_MAX_LIGHTS = %d\n",i);
    }

    /* Determine depth thresholds based on resolution of depth buffer
     */
    set_thresholds();
    print_sep(stdout,'-');
    fprintf(stdout,"Edge Detector Thresholds:\n");
    {
        GLint i;

        fprintf(stdout,"\tth_H = %d\n",th_H);
        fprintf(stdout,"\tth_L = %d\n",th_L);
        fprintf(stdout,"\tth_s = %d\n",th_s);
    }

    /* Determine epsilons based on resolution of depth buffer
     */
    set_epsilons();
    print_sep(stdout,'-');
    fprintf(stdout,"Shadow Map Bias (Relative):\n");
    {
        GLint i;

        fprintf(stdout,"\tlight_ze = %f\n",light_ze);
        fprintf(stdout,"\teye_ze = %f\n",eye_ze);
    }

    print_sep(stdout,'-');
    keyboard_help();
    fflush(stdout);

    /* Initialize timers
     */
    init_timers();

    /* Enter GLUT event processing loop.
     */
    glutSetWindow(main_window);
    update_window_title();
    glutMainLoop();
    return 0;
}
