/home/carolina/workspace/spvolren_cpu/main.c

Vá para a documentação deste arquivo.

/*
 * Copyright (c) 2005  Institute for Visualization and Interactive
 * Systems, University of Stuttgart, Germany
 *
 * This source code is distributed as part of the single-pass volume
 * rendering project. Details about this project can be found on the
 * project web page at http://www.vis.uni-stuttgart.de/eng/research/
 * fields/current/spvolren. This file may be distributed, modified,
 * and used free of charge as long as this copyright notice is
 * included in its original form. Commercial use is strictly
 * prohibited.
 *
 * Filename: main.c
 * 
 */
#include "carol_glsl.h"

#define DEBUG_OPENGL                            1

#define VOLUME_SHADERS_DIR                      "./glsl_volumeShaders"
#define SHADERS_DIR                             "./shaders"
#define BACKGROUND_IMAGE                        "backgrounds/chess.ppm"
#define SPHEREMAP                               "spheremaps/grace_probe.ppm"

#define FORCE_POWER_OF_TWO_TEXTURE      0
#define NORMALIZE_VOLUME_DATA           1
#define WRITE_FLOAT_GRADIENTS           0

#define SOBEL                           1
#define GRAD_FILTER_SIZE    5
#define SIGMA2                          5.0

#define MOUSE_SCALE         10.0

#define FRAG_PROG_EXT           ".fp"
#define AA_FRAG_PROG_EXT        ".fpa"
#define VOL_FILE_EXT            ".dat"
#define SETTINGS_EXT            ".stg"
#define GRADIENTS_EXT           ".grd"

#define MAX_FRAMES              50

#define WINDOW_WIDTH            512
#define WINDOW_HEIGHT           512

#define FOVY                            60.0
#define NEAR_CLIP                       0.1
#define FAR_CLIP                        3.0

#define SLICE_STEP                      .1
#define MAX_SLICE_THICKNESS 10.0
#define MIN_SLICE_THICKNESS SLICE_STEP

#define STEPSIZE_STEP           .001
#define MAX_STEPSIZE            1.0
#define MIN_STEPSIZE            0.0

#define GRAD_SCALE_STEP         .01
#define GRAD_OFFSET_STEP        .1

#define TEX_COORD_SCALE_STEP    .05

#define ISO_VALUE_STEP          .005

#define SCATTERING_STEP         .5

#define LIGHT_POS_STEP          1.0

#define AA_THRESHOLD_STEP       0.005

#define TEXTURE_TRANSFERFUNCTION        "TRANSFERFUNCTION"
#define TEXTURE_OPTICALDENSITY          "OPTICALDENSITY"
#define TEXTURE_BACKGROUND              "BACKGROUND"
#define TEXTURE_CLIPVOLUME              "CLIPVOLUME"
#define TEXTURE_SCATTERING              "SCATTERING"
#define TEXTURE_SPHEREMAP               "SPHEREMAP"
#define TEXTURE_VOLUME                  "VOLUME"

#ifndef SQR
#define SQR(a) ((a) * (a))
#endif
#define MAX(a,b) ((a) > (b) ? (a) : (b))
#define MIN(a,b) ((a) < (b) ? (a) : (b))

/* Determine the current framerate continously */
#define CONTINOUS_FPS           0

#if DEBUG_OPENGL == 1
#define CHECK_FOR_OGL_ERROR()                                           \
do {                                                                    \
    GLenum err;                                                         \
        err = glGetError();                                             \
        if (err != GL_NO_ERROR)                                         \
        {                                                               \
            fprintf(stderr, "%s(%d) glError: %s\n",                     \
                    __FILE__, __LINE__, gluErrorString(err));           \
        }                                                               \
} while(0)
#else 
#define CHECK_FOR_OGL_ERROR() 
#endif

int getNextPowerOfTwo(int n)
{
        int i;

        i = 1;
        while (i < n) {
                i *= 2;
        }

        return i;
}

unsigned char getVoxel8(int x, int y, int z)
{
        return ((unsigned char*)g.volData)[z * g.numSlices[0] * g.numSlices[1] + 
                                                                           y * g.numSlices[0] + 
                                                                           x];
}

unsigned short getVoxel16(int x, int y, int z)
{
        return ((unsigned short*)g.volData)[z * g.numSlices[0] * g.numSlices[1] + 
                                                                                y * g.numSlices[0] + 
                                                                                x];
}

float getVoxel(int x, int y, int z, DataType dataType)
{
        switch (dataType) {
                case DATRAW_UCHAR:
                        return (float)getVoxel8(x, y, z);
                        break;
                case DATRAW_USHORT:
                        return (float)getVoxel16(x, y, z);
                        break;
                default:
                        fprintf(stderr, "Unsupported data type\n");
                        exit(1);
                        break;
        }
        return 0.0;
}

char *getGradientsFilename(void)
{
    char *filename;

    if (! (filename = (char *)malloc(strlen(g.basename) +
                                     strlen(GRADIENTS_EXT) + 1))) {
        fprintf(stderr, "not enough memory for filename\n");
        exit(1);
    }

    strcpy(filename, g.basename);
    strcat(filename, GRADIENTS_EXT);

    return filename;
}

int loadGradients(void *gradients, int *sizes, DataType dataType)
{
    char *filename;
        int size;
    FILE *fp;

    filename = getGradientsFilename();
    if (! (fp = fopen(filename, "rb"))) {
        fprintf(stderr, "no pre-computed gradients found\n");
        return 0;
    }

        size = 3 * sizes[0] * sizes[1] * sizes[2]
                   * getDataTypeSize(dataType);
    if (fread(gradients, size, 1, fp) != 1) {
        fprintf(stderr, "reading gradients failed\n");
        exit(1);
    }

    fclose(fp);

    return 1;
}

void saveGradients(void *gradients, int *sizes, DataType dataType)
{
    char *filename;
        int size;
    FILE *fp;

    filename = getGradientsFilename();
    if (! (fp = fopen(filename, "wb"))) {
        perror("cannot open gradients file for writing");
        exit(1);
    }

        size = 3 * sizes[0] * sizes[1] * sizes[2]
                   * getDataTypeSize(dataType);
        
    if (fwrite(gradients, size, 1, fp) != 1) {
        fprintf(stderr, "writing gradients failed\n");
        exit(1);
    }

    fclose(fp);
}

#if WRITE_FLOAT_GRADIENTS == 1
char *getFloatGradientsFilename(void)
{
        char floatExt[] = "_float";
    char *filename;

    if (! (filename = (char *)malloc(strlen(g.basename) +
                                                                         strlen(floatExt) +
                                     strlen(GRADIENTS_EXT) + 1))) {
        fprintf(stderr, "not enough memory for filename\n");
        exit(1);
    }

    strcpy(filename, g.basename);
        strcat(filename, floatExt);
    strcat(filename, GRADIENTS_EXT);

    return filename;
}

void saveFloatGradients(float *gradients, int *sizes)
{
    char *filename;
    FILE *fp;

    filename = getFloatGradientsFilename();
    if (! (fp = fopen(filename, "wb"))) {
        perror("cannot open gradients file for writing");
        exit(1);
    }

    if (fwrite(gradients, 3 * sizes[0] * sizes[1] * sizes[2] * sizeof(float), 
                           1, fp) != 1) {
        fprintf(stderr, "writing float gradients failed\n");
        exit(1);
    }

    fclose(fp);
}
#endif

void computeGradients(float *gradients, int *sizes, DataType dataType)
{
    int i, j, k, dir, di, vdi, idz, idy, idx;
    float *gp;

    static int weights[][3][3][3] = {
        {{{-1, -3, -1},
          {-3, -6, -3},
          {-1, -3, -1}},
         {{ 0,  0,  0},
          { 0,  0,  0},
          { 0,  0,  0}},
         {{ 1,  3,  1},
          { 3,  6,  3},
          { 1,  3,  1}}},
        {{{-1, -3, -1},
          { 0,  0,  0},
          { 1,  3,  1}},
         {{-3, -6, -3},
          { 0,  0,  0},
          { 3,  6,  3}},
         {{-1, -3, -1},
          { 0,  0,  0},
          { 1,  3,  1}}},
        {{{-1,  0,  1},
          {-3,  0,  3},
          {-1,  0,  1}},
         {{-3,  0,  3},
          {-6,  0,  6},
          {-3,  0,  3}},
         {{-1,  0,  1},
          {-3,  0,  3},
          {-1,  0,  1}}}
    };

        fprintf(stderr, "computing gradients ... may take a while\n");

    di = 0;
    vdi = 0;
        gp = gradients;
    for (idz = 0; idz < sizes[2]; idz++) {
        for (idy = 0; idy < sizes[1]; idy++) {
            for (idx = 0; idx < sizes[0]; idx++) {
#if SOBEL == 1
                if (idx > 0 && idx < sizes[0] - 1 &&
                    idy > 0 && idy < sizes[1] - 1 &&
                    idz > 0 && idz < sizes[2] - 1) {

                    for (dir = 0; dir < 3; dir++) {
                        gp[dir] = 0.0;
                        for (i = -1; i < 2; i++) {
                            for (j = -1; j < 2; j++) {
                                for (k = -1; k < 2; k++) {
                                    gp[dir] += weights[dir][i + 1]
                                                           [j + 1]
                                                           [k + 1] *
                                               getVoxel(idx + i,
                                                      idy + j,
                                                      idz + k,
                                                                                                          dataType);
                                }
                            }
                        }

                                                gp[dir] /= 2.0 * g.sliceDists[dir];
                    }
                } else {
                    /* X-direction */
                    if (idx < 1) {
                        gp[0] = (getVoxel(idx + 1, idy, idz, dataType) -
                                 getVoxel(idx, idy, idz, dataType))/
                                (g.sliceDists[0]);
                    } else {
                        gp[0] = (getVoxel(idx, idy, idz, dataType) -
                                 getVoxel(idx - 1, idy, idz, dataType))/
                                (g.sliceDists[0]);
                    }

                    /* Y-direction */
                    if (idy < 1) {
                        gp[1] = (getVoxel(idx, idy + 1, idz, dataType) -
                                   getVoxel(idx, idy, idz, dataType))/
                                   (g.sliceDists[1]);
                    } else {
                        gp[1] = (getVoxel(idx, idy, idz, dataType) -
                                   getVoxel(idx, idy - 1, idz, dataType))/
                                   (g.sliceDists[1]);
                    }

                    /* Z-direction */
                    if (idz < 1) {
                        gp[2] = (getVoxel(idx, idy, idz + 1, dataType) -
                                 getVoxel(idx, idy, idz, dataType))/
                                (g.sliceDists[2]);
                    } else {
                        gp[2] = (getVoxel(idx, idy, idz, dataType) -
                                 getVoxel(idx, idy, idz - 1, dataType))/
                                (g.sliceDists[2]);
                    }
                }
#else
                /* X-direction */
                if (idx < 1) {
                    gp[0] = (getVoxel(idx + 1, idy, idz, dataType) -
                             getVoxel(idx, idy, idz, dataType))/
                            (g.sliceDists[0]);
                } else if (idx > g.numSlices[0] - 1) {
                    gp[0] = (getVoxel(idx, idy, idz, dataType) -
                             getVoxel(idx - 1, idy, idz, dataType))/
                            (g.sliceDists[0]);
                } else {
                    gp[0] = (getVoxel(idx + 1, idy, idz, dataType) -
                             getVoxel(idx - 1, idy, idz, dataType))/
                            (2.0 * g.sliceDists[0]);
                }

                /* Y-direction */
                if (idy < 1) {
                    gp[1] = (getVoxel(idx, idy + 1, idz, dataType) -
                             getVoxel(idx, idy, idz, dataType))/
                            (g.sliceDists[1]);
                } else if (idy > g.numSlices[1] - 1) {
                    gp[1] = (getVoxel(idx, idy, idz, dataType) -
                             getVoxel(idx, idy - 1, idz, dataType))/
                            (g.sliceDists[1]);
                } else {
                    gp[1] = (getVoxel(idx, idy + 1, idz, dataType) -
                             getVoxel(idx, idy - 1, idz, dataType))/
                            (2.0 * g.sliceDists[1]);
                }

                /* Z-direction */
                if (idz < 1) {
                    gp[2] = (getVoxel(idx, idy, idz + 1, dataType) -
                             getVoxel(idx, idy, idz, dataType))/
                            (g.sliceDists[2]);
                } else if (idz > g.numSlices[2] - 1) {
                    gp[2] = (getVoxel(idx, idy, idz, dataType) -
                             getVoxel(idx, idy, idz - 1, dataType))/
                            (g.sliceDists[2]);
                } else {
                    gp[2] = (getVoxel(idx, idy, idz + 1, dataType) -
                             getVoxel(idx, idy, idz - 1, dataType))/
                            (2.0 * g.sliceDists[2]);
                }
#endif
                gp += 3;
            }
        }
    }

}

void filterGradients(float *gradients, int *sizes)
{
    int i, j, k, idz, idy, idx, gi, ogi, filterWidth, n, borderDist[3];
    float sum, *filteredGradients, ****filter;

        fprintf(stderr, "filtering gradients ... may also take a while\n");

    if (! (filteredGradients = (float *)malloc(sizes[0] * sizes[1] * sizes[2]
                                                                                           * 3 * sizeof(float)))) {
        fprintf(stderr, "not enough memory for filtered gradients\n");
        exit(1);
    }

        /* Allocate storage for filter kernels */
        if (! (filter = (float ****)malloc((GRAD_FILTER_SIZE/2 + 1) * 
                                                                           sizeof(float ***)))) {
                fprintf(stderr, "failed to allocate gradient filter\n");
                exit(1);
        }
        
        for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
                if (! (filter[i] = (float ***)malloc((GRAD_FILTER_SIZE) * 
                                                                                         sizeof(float **)))) {
                        fprintf(stderr, "failed to allocate gradient filter\n");
                        exit(1);
                }
        }
        for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
                for (j = 0; j < GRAD_FILTER_SIZE; j++) {
                        if (! (filter[i][j] = (float **)malloc((GRAD_FILTER_SIZE) * 
                                                                                                   sizeof(float *)))) {
                                fprintf(stderr, "failed to allocate gradient filter\n");
                                exit(1);
                        }
                }
        }
        for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
                for (j = 0; j < GRAD_FILTER_SIZE; j++) {
                        for (k = 0; k < GRAD_FILTER_SIZE; k++) {
                                if (! (filter[i][j][k] = (float *)malloc((GRAD_FILTER_SIZE) *
                                                                                                                 sizeof(float)))) {
                                        fprintf(stderr, "failed to allocate gradient filter\n");
                                        exit(1);
                                }
                        }
                }
        }

    filterWidth = GRAD_FILTER_SIZE/2;
        
        /* Compute the filter kernels */
    for (n = 0; n <= filterWidth; n++) {
        sum = 0.0f;
        for (k = -filterWidth; k <= filterWidth; k++) {
            for (j = -filterWidth; j <= filterWidth; j++) {
                for (i = -filterWidth; i <= filterWidth; i++) {
                    sum += (filter[n][filterWidth + k]
                                                                         [filterWidth + j]
                                                                         [filterWidth + i] =
                                                        exp(-(SQR(i) + SQR(j) + SQR(k)) / SIGMA2));
                }
            }
        }
        for (k = -filterWidth; k <= filterWidth; k++) {
            for (j = -filterWidth; j <= filterWidth; j++) {
                for (i = -filterWidth; i <= filterWidth; i++) {
                    filter[n][filterWidth + k]
                                                         [filterWidth + j]
                                                         [filterWidth + i] /= sum;
                }
            }
        }
    }

    gi = 0;
    /* Filter the gradients */
    for (idz = 0; idz < sizes[2]; idz++) {
        for (idy = 0; idy < sizes[1]; idy++) {
            for (idx = 0; idx < sizes[0]; idx++) {
                borderDist[0] = MIN(idx, sizes[0] - idx - 1);
                borderDist[1] = MIN(idy, sizes[1] - idy - 1);
                borderDist[2] = MIN(idz, sizes[2] - idz - 1);

                filterWidth = MIN(GRAD_FILTER_SIZE/2,
                                  MIN(MIN(borderDist[0], 
                                                                                  borderDist[1]),
                                                                                  borderDist[2]));

                for (n = 0; n < 3; n++) {
                    filteredGradients[gi] = 0.0;
                    for (k = -filterWidth; k <= filterWidth; k++) {
                        for (j = -filterWidth; j <= filterWidth; j++) {
                            for (i = -filterWidth; i <= filterWidth; i++) {
                                ogi = (((idz + k) * sizes[1]  + (idy + j)) *
                                                                          sizes[0] + (idx + i)) * 3 + 
                                                                          n;
                                filteredGradients[gi] += filter[filterWidth]
                                                                                                           [filterWidth + k]
                                                                                                           [filterWidth + j]
                                                                                                           [filterWidth + i] * 
                                                                                                           gradients[ogi];
                            }
                        }
                    }
                    gi++;
                }
            }
        }
    }

        /* Replace the original gradients by the filtered gradients */
        memcpy(gradients, filteredGradients, 
                   sizes[0] * sizes[1] * sizes[2] * 3 * sizeof(float));

    free(filteredGradients);

        /* free storage of filter kernel(s) */
        for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
                for (j = 0; j < GRAD_FILTER_SIZE; j++) {
                        for (k = 0; k < GRAD_FILTER_SIZE; k++) {
                                free(filter[i][j][k]);
                        }
                }
        }
        for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
                for (j = 0; j < GRAD_FILTER_SIZE; j++) {
                        free(filter[i][j]);
                }
        }
        for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
                free(filter[i]);
        }
        free(filter);
}

void quantize8(float *grad, unsigned char *data)
{
        float len;
        int i;

        len = sqrt(SQR(grad[0]) + SQR(grad[1]) + SQR(grad[2]));

        if (len < EPS) {
                grad[0] = grad[1] = grad[2] = 0.0;
        } else {
                grad[0] /= len;
                grad[1] /= len;
                grad[2] /= len;
        }

        for (i = 0; i < 3; i++) {
                data[i] = (unsigned char)((grad[i] + 1.0)/2.0 * UCHAR_MAX);
        }
}

void quantize16(float *grad, unsigned short *data)
{
        float len;
        int i;

        len = sqrt(SQR(grad[0]) + SQR(grad[1]) + SQR(grad[2]));

        if (len < EPS) {
                grad[0] = grad[1] = grad[2] = 0.0;
        } else {
                grad[0] /= len;
                grad[1] /= len;
                grad[2] /= len;
        }

        for (i = 0; i < 3; i++) {
                data[i] = (unsigned short)((grad[i] + 1.0)/2.0 * USHRT_MAX);
        }
}

void quantizeGradients(float *gradientsIn, void *gradientsOut,
                                           int *sizes, DataType dataType)
{
        int idx, idy, idz, di;
        
    di = 0;
    for (idz = 0; idz < sizes[2]; idz++) {
        for (idy = 0; idy < sizes[1]; idy++) {
            for (idx = 0; idx < sizes[0]; idx++) {
                                switch (dataType) {
                                case DATRAW_UCHAR:
                                        quantize8(&gradientsIn[di], 
                                                          &((unsigned char*)gradientsOut)[di]);
                                        break;
                                case DATRAW_USHORT:
                                        quantize16(&gradientsIn[di], 
                                                           &((unsigned short*)gradientsOut)[di]);
                                        break;
                                default:
                                        fprintf(stderr, "unsupported data type\n");
                                        break;
                                }
                di += 3;
            }
        }
    }
}

void addTexture(GLuint texId, GLenum texTarget, const char *texName)
{
        Texture *texture;

        g.numTextures++;
        if (! (g.textures = (Texture *)realloc(g.textures, g.numTextures *
                                                                                   sizeof(Texture)))) {
                fprintf(stderr, "not enough memory for textures\n");
                exit(1);
        }
        texture = &g.textures[g.numTextures - 1];

        strcpy(texture->name, texName);
        texture->target = texTarget;
        texture->id = texId;
}

Texture *getTexture(const char *texName)
{
        Texture *texture;
        int i;
        
        for (i = 0; i < g.numTextures; i++) {
                texture = g.textures + i;
                if (!strcmp(texture->name, texName)) {
                        return texture;
                }
        }
        return NULL;
}

void loadVolumeTexture8()
{
        unsigned char *voldata = (unsigned char*)g.volData;
        unsigned char *data, *gradients, *gp, min, max;
        int di, vdi, idz, idy, idx,  haveGradients;
        GLuint texId;
        
    if (! (data = (unsigned char *)calloc(g.volTexSizes[0] *
                                          g.volTexSizes[1] *
                                          g.volTexSizes[2] * 4,
                                          sizeof(unsigned char)))) {
        fprintf(stderr, "not enough memory for volume texture\n");
        exit(1);
    }

    if (! (gradients = (unsigned char *)malloc(g.numSlices[0] *
                                               g.numSlices[1] *
                                               g.numSlices[2] * 
                                                                                           3 * sizeof(unsigned char)))) {
        fprintf(stderr, "not enough memory for gradients\n");
        exit(1);
    }

    if (! (haveGradients = loadGradients(gradients, g.numSlices, DATRAW_UCHAR))) {
                float *tempGradients;

                if (! (tempGradients = (float *)malloc(g.numSlices[0] *
                                                                                           g.numSlices[1] * 
                                                                                           g.numSlices[2] * 
                                                                                           3 * sizeof(float)))) {
                        fprintf(stderr, "not enough memory for temporary gradients\n");
                        exit(1);
                }
        computeGradients(tempGradients, g.numSlices, DATRAW_UCHAR);
        filterGradients(tempGradients, g.numSlices);
#if WRITE_FLOAT_GRADIENTS == 1
        saveFloatGradients(tempGradients, g.numSlices);
#endif
                quantizeGradients(tempGradients, gradients, g.numSlices, DATRAW_UCHAR);
        saveGradients(gradients, g.numSlices, DATRAW_UCHAR);

                free(tempGradients);
    }

        /* find min and max of scalar values */
        di = 0;
        min = UCHAR_MAX;
        max = 0;
        for (idz = 0; idz < g.numSlices[2]; idz++) {
        for (idy = 0; idy < g.numSlices[1]; idy++) {
            for (idx = 0; idx < g.numSlices[0]; idx++) {
                                if (voldata[di] > max) {
                                        max = voldata[di];
                                }
                                if (voldata[di] < min) {
                                        min = voldata[di];
                                }
                                di++;
                        }
                }
        }
        printf("data range: %d - %d\n", (int)min, (int)max);

    di = 0;
    vdi = 0;
    gp = gradients;
        // Pack the gradients and scalar values into a single texture
    for (idz = 0; idz < g.volTexSizes[2]; idz++) {
        for (idy = 0; idy < g.volTexSizes[1]; idy++) {
            for (idx = 0; idx < g.volTexSizes[0]; idx++) {
                if (idx < g.numSlices[0] &&
                    idy < g.numSlices[1] &&
                    idz < g.numSlices[2]) {
                        memcpy(&data[di], gp, 3);
                        gp += 3;
#if NORMALIZE_VOLUME_DATA == 1
                        data[di + 3] = (unsigned char) 
                                                        ((voldata[vdi++] - min)/(float)(max - min) * UCHAR_MAX);
#else
                                                data[di + 3] = voldata[vdi++];
#endif
                }
                                di += 4;
            }
        }
    }

    free(gradients);

    glGenTextures(1, &texId);
    glBindTexture(GL_TEXTURE_3D, texId);
    glTexImage3DEXT(GL_TEXTURE_3D, 0, GL_RGBA, g.volTexSizes[0],
                    g.volTexSizes[1], g.volTexSizes[2], 0, GL_RGBA,
                    GL_UNSIGNED_BYTE, data);
        addTexture(texId, GL_TEXTURE_3D, TEXTURE_VOLUME);

    glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
        free(data);
}

void loadVolumeTexture16()
{
        unsigned short *voldata = (unsigned short*)g.volData;
        unsigned short *data, *gradients, *gp, min, max;
        int di, vdi, idz, idy, idx,  haveGradients;
        GLuint texId;

    if (! (data = (unsigned short *)calloc(g.volTexSizes[0] *
                                           g.volTexSizes[1] *
                                           g.volTexSizes[2] * 4,
                                           sizeof(unsigned short)))) {
        fprintf(stderr, "not enough memory for volume texture\n");
        exit(1);
    }

    if (! (gradients = (unsigned short *)malloc(g.numSlices[0] *
                                                g.numSlices[1] *
                                                g.numSlices[2] * 3
                                                                                           * sizeof(unsigned short)))) {
        fprintf(stderr, "not enough memory for gradients\n");
        exit(1);
    }

    if (! (haveGradients = loadGradients(gradients, g.numSlices, DATRAW_USHORT))) {
                float *tempGradients;

                if (! (tempGradients = (float *)malloc(g.numSlices[0] *
                                                                                           g.numSlices[1] * 
                                                                                           g.numSlices[2] * 
                                                                                           3 * sizeof(float)))) {
                        fprintf(stderr, "not enough memory for temporary gradients\n");
                        exit(1);
                }
        computeGradients(tempGradients, g.numSlices, DATRAW_USHORT);
        filterGradients(tempGradients, g.numSlices);
#if WRITE_FLOAT_GRADIENTS == 1
        saveFloatGradients(tempGradients, g.numSlices);
#endif
                quantizeGradients(tempGradients, gradients, g.numSlices, DATRAW_USHORT);
        saveGradients(gradients, g.numSlices, DATRAW_USHORT);

                free(tempGradients);
    }

        /* find min and max of scalar values */
        di = 0;
        min = USHRT_MAX;
        max = 0;
        for (idz = 0; idz < g.numSlices[2]; idz++) {
        for (idy = 0; idy < g.numSlices[1]; idy++) {
            for (idx = 0; idx < g.numSlices[0]; idx++) {
                                if (voldata[di] > max) {
                                        max = voldata[di];
                                } 
                                if (voldata[di] < min) {
                                        min = voldata[di];
                                }
                                di++;
                        }
                }
        }
        printf("data range: %d - %d\n", (int)min, (int)max);

        /* Pack the gradients and scalar values into a single texture */
    di = 0;
    vdi = 0;
    gp = gradients;
    for (idz = 0; idz < g.volTexSizes[2]; idz++) {
        for (idy = 0; idy < g.volTexSizes[1]; idy++) {
            for (idx = 0; idx < g.volTexSizes[0]; idx++) {
                if (idx < g.numSlices[0] &&
                    idy < g.numSlices[1] &&
                    idz < g.numSlices[2]) {
                        memcpy(&data[di], gp, 6);
                        gp += 3;
#if NORMALIZE_VOLUME_DATA == 1
                        data[di + 3] = (unsigned short) 
                                                        ((voldata[vdi++] - min)/(float)(max - min) * USHRT_MAX);
#else
                                                data[di + 3] = voldata[vdi++];
#endif
                }
                                di += 4;
            }
        }
    }

    free(gradients);
    
    glGenTextures(1, &texId);
    glBindTexture(GL_TEXTURE_3D, texId);
        glTexImage3DEXT(GL_TEXTURE_3D, 0, GL_RGBA16, g.volTexSizes[0],
                    g.volTexSizes[1], g.volTexSizes[2], 0, GL_RGBA,
                    GL_UNSIGNED_SHORT, data);

        addTexture(texId, GL_TEXTURE_3D, TEXTURE_VOLUME);

    glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);

        free(data);
}

void loadVolumeTexture(char *filename)
{
    int i, numComponents, maxVolTexSize;
        float volSizes[3], maxVolSize;
    DataType dataType;

        readData(filename, g.numSlices, g.sliceDists, (void **)&g.volData,
                         &dataType, &numComponents);

    for (i = 0; i < 3; i++) {
#if FORCE_POWER_OF_TWO_TEXTURE == 1
        g.volTexSizes[i] = getNextPowerOfTwo(g.numSlices[i]);
#else
                g.volTexSizes[i] = g.numSlices[i];
#endif
    }

        fprintf(stderr, "volume: %i %i %i\n",
            g.volTexSizes[0],
            g.volTexSizes[1],
            g.volTexSizes[2]);

        if (numComponents != 1) {
        fprintf(stderr, "invalid volume data set\n");
        exit(1);
    }

        switch(dataType) {
        case DATRAW_UCHAR:
                loadVolumeTexture8();
                break;
        case DATRAW_USHORT:
                loadVolumeTexture16();
                break;
        default:
                fprintf(stderr, "Only 8bit and 16bit data format supported\n");
                exit(1);
        }

    for (i = 0; i < 3; i++) {
        volSizes[i] = g.numSlices[i] * g.sliceDists[i];
    }

    maxVolSize = MAX(MAX(volSizes[0], volSizes[1]), volSizes[2]);
    maxVolTexSize = MAX(MAX(g.volTexSizes[0], g.volTexSizes[1]), 
                                                g.volTexSizes[2]);

    for (i = 0; i < 3; i++) {
        g.scaleFactors[i] = maxVolSize/(g.volTexSizes[i] * g.sliceDists[i]);
        g.extents[i] = (g.numSlices[i] * g.sliceDists[i])/maxVolSize;
                g.scaleFactorsInv[i] = 1.0/g.scaleFactors[i];
        g.center[i] = g.extents[i]/2.0;
    }
    g.scaleFactorsInv[3] = 1.0;
    g.scaleFactors[3] = 0.0;
    g.center[3] = 0.0;
}

void updateTransferFunction(void)
{
        int i, j, index;

        index = 0;
        for (j = 0; j < NUM_TE_ENTRIES; j++) {
                /* First four arrays are RGBA */
                for (i = 0; i < 4; i++) {
                        g.tfFunc[index++] = g.dataTE[i][j];
                }
        }
}

void loadTETexture(void)
{
        static int initialized = 0;
        static GLuint texId;

        if (! initialized) {
                glGenTextures(1, &texId);
                addTexture(texId, GL_TEXTURE_2D, TEXTURE_TRANSFERFUNCTION);
                initialized = 1;
        }

    glBindTexture(GL_TEXTURE_2D, texId);
        glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, NUM_TE_ENTRIES, NUM_TE_ENTRIES, 
                                 0, GL_RGBA, GL_UNSIGNED_BYTE, g.preIntTable);

    glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
}

void loadClipTexture(char* filename)
{
        int numComponents, i, sizes[3], volTexSizes[3];
        unsigned char *temp, *data, noClipData = 255;
        DataType dataType;
        float dists[3];
        GLuint texId;
#if FORCE_POWER_OF_TWO_TEXTURE == 1             
        int idx, idy, idz;
#endif

        if (! filename) {
                data = &noClipData;
                for (i = 0; i < 3; i++) {
                        volTexSizes[i] = 1;
                }
        } else {
                readData(filename, sizes, dists, (void**)&temp, &dataType, &numComponents);

                if (dataType != DATRAW_UCHAR || numComponents != 1) {
                        fprintf(stderr, "invalid volume data set\n");
                        exit(1);
                }

                for (i = 0; i < 3; i++) {
#if FORCE_POWER_OF_TWO_TEXTURE == 1             
                        volTexSizes[i] = getNextPowerOfTwo(sizes[i]);
#else
                        volTexSizes[i] = sizes[i];
#endif
                }

                if (! (data = (unsigned char *)malloc(volTexSizes[0] *
                                                                                          volTexSizes[1] *
                                                                                          volTexSizes[2]))) {
                        fprintf(stderr, "not enough memory for clip volume data\n");
                        exit(1);
                }

#if FORCE_POWER_OF_TWO_TEXTURE == 1             
                i = 0;
                for (idz = 0; idz < volTexSizes[2]; idz++) {
                        for (idy = 0; idy < volTexSizes[1]; idy++) {
                                for (idx = 0; idx < volTexSizes[0]; idx++) {
                                        if (idx < sizes[0] && idy < sizes[1] && idz < sizes[2]) {
                                                data[i] = temp[(idz * sizes[1] + idy) *
                                                                           sizes[0] + idx];
                                        } else {
                                                data[i] = 0;
                                        }
                                        i++;
                                }
                        }
                }
#else
                memcpy(data, temp, sizes[0] * sizes[1] * sizes[2]);
#endif
        }

        glGenTextures(1, &texId);
    glBindTexture(GL_TEXTURE_3D, texId);
    glTexImage3DEXT(GL_TEXTURE_3D, 0, GL_LUMINANCE, volTexSizes[0],
                    volTexSizes[1], volTexSizes[2], 0, GL_LUMINANCE,
                    GL_UNSIGNED_BYTE, data);

        addTexture(texId, GL_TEXTURE_3D, TEXTURE_CLIPVOLUME);

    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);

        if (filename) {
        free(data);
                free(temp);
        }
}

void loadBackgroundTexture(void)
{
        static int initialized = 0;
        static PPMFile image;
        static GLuint texId;

        if (! initialized) {
                glGenTextures(1, &texId);
                addTexture(texId, GL_TEXTURE_2D, TEXTURE_BACKGROUND);
                ppmRead(BACKGROUND_IMAGE, &image);
                initialized = 1;
        }

    glBindTexture(GL_TEXTURE_2D, texId);

        if (g.u.backgroundImage) {
        fprintf(stderr, "%i\n", g.u.backgroundImage);
                glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, image.width, image.height, 
                                         0, GL_RGB, GL_UNSIGNED_BYTE, image.data);
        } else {
                unsigned char rgb[3];
                memset(rgb, g.u.backgroundGrayVal, 3);
                glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, 1, 1, 
                                         0, GL_RGB, GL_UNSIGNED_BYTE, rgb);
        }

    glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
}

void loadSphereMapTexture(void)
{
        PPMFile image;
        GLuint texId;

        glGenTextures(1, &texId);
        ppmRead(SPHEREMAP, &image);

    glBindTexture(GL_TEXTURE_2D, texId);
        glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, image.width, image.height, 
                                 0, GL_RGB, GL_UNSIGNED_BYTE, image.data);

        addTexture(texId, GL_TEXTURE_2D, TEXTURE_SPHEREMAP);

    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
}

void vertex(float x, float y, float z)
{
        glMultiTexCoord3fARB(GL_TEXTURE0_ARB, x, y, z);
        glVertex3f(x, y, z);
}

void vertexv(double *coords)
{
        glMultiTexCoord3dvARB(GL_TEXTURE0_ARB, coords);
        glVertex3dv(coords);
}

void drawQuads(void)
{
        glBegin(GL_QUADS);
                /* Back side */
                glNormal3f(0.0, 0.0, -1.0);
                glMultiTexCoord4fARB(GL_TEXTURE5_ARB, 0.0, 0.0, -1.0, 0.0);
                vertex(0.0, 0.0, 0.0);
                vertex(0.0, g.extents[1], 0.0);
                vertex(g.extents[0], g.extents[1], 0.0);
                vertex(g.extents[0], 0.0, 0.0);

                /* Front side */
                glNormal3f(0.0, 0.0, 1.0);
                glMultiTexCoord4fARB(GL_TEXTURE5_ARB, 0.0, 0.0, 1.0, 0.0);
                vertex(0.0, 0.0, g.extents[2]);
                vertex(g.extents[0], 0.0, g.extents[2]);
                vertex(g.extents[0], g.extents[1], g.extents[2]);
                vertex(0.0, g.extents[1], g.extents[2]);

                /* Top side */
                glNormal3f(0.0, 1.0, 0.0);
                glMultiTexCoord4fARB(GL_TEXTURE5_ARB, 0.0, 1.0, 0.0, 0.0);
                vertex(0.0, g.extents[1], 0.0);
                vertex(0.0, g.extents[1], g.extents[2]);
                vertex(g.extents[0], g.extents[1], g.extents[2]);
                vertex(g.extents[0], g.extents[1], 0.0);

                /* Bottom side */
                glNormal3f(0.0, -1.0, 0.0);
                glMultiTexCoord4fARB(GL_TEXTURE5_ARB, 0.0, -1.0, 0.0, 0.0);
                vertex(0.0, 0.0, 0.0);
                vertex(g.extents[0], 0.0, 0.0);
                vertex(g.extents[0], 0.0, g.extents[2]);
                vertex(0.0, 0.0, g.extents[2]);

                /* Left side */
                glNormal3f(-1.0, 0.0, 0.0);
                glMultiTexCoord4fARB(GL_TEXTURE5_ARB, -1.0, 0.0, 0.0, 0.0);
                vertex(0.0, 0.0, 0.0);
                vertex(0.0, 0.0, g.extents[2]);
                vertex(0.0, g.extents[1], g.extents[2]);
                vertex(0.0, g.extents[1], 0.0);

                /* Right side */
                glNormal3f(1.0, 0.0, 0.0);
                glMultiTexCoord4fARB(GL_TEXTURE5_ARB, 1.0, 0.0, 0.0, 0.0);
                vertex(g.extents[0], 0.0, 0.0);
                vertex(g.extents[0], g.extents[1], 0.0);
                vertex(g.extents[0], g.extents[1], g.extents[2]);
                vertex(g.extents[0], 0.0, g.extents[2]);
        glEnd();
}

void drawLight(void)
{
        glDisable(GL_LIGHTING);
        
        glColor4f(1.0, 1.0, 0.0, 1.0);
        
        glDisable(GL_CULL_FACE);

        glPushMatrix();
        glLoadIdentity();
        glTranslatef(0.0, 0.0, -g.u.camZ);
        glBegin(GL_LINES);
        glVertex3f(0.0, 0.0, 0.0);
        glVertex3f(g.u.lightPos[0], g.u.lightPos[1], g.u.lightPos[2]);
        glEnd();
        glTranslatef(g.u.lightPos[0], g.u.lightPos[1], g.u.lightPos[2]);
        glutSolidSphere(.1, 8, 8);
        glPopMatrix();
}

void drawWireframe(void)
{
        glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);

        glDisable(GL_CULL_FACE);

        glColor3f(2.0, 3.0, 5.0);
        drawQuads();

        glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
}

void setLight(void)
{
        GLfloat transLight[4];
        GLfloat invMVM[16];
        Vector3 axis;
        float angle;

        Quaternion_getAngleAxis(g.u.camRot, &angle, &axis);

        glPushMatrix();
        glLoadIdentity();
        glTranslatef(g.center[0], g.center[1], g.center[2]);
        glRotatef(-angle * 180.0 / M_PI, axis.x, axis.y, axis.z);
        glRotatef(-g.animatedAngle, 0.0, 1.0, 0.0);
        glTranslatef(0.0, 0.0, g.u.camZ);
        glTranslatef(-g.u.translate[0], -g.u.translate[1], -g.u.translate[2]);

        glGetFloatv(GL_MODELVIEW_MATRIX, invMVM);


        glPopMatrix();

        
        transLight[0] = invMVM[0] * g.u.lightPos[0] + 
                                        invMVM[4] * g.u.lightPos[1] + 
                                        invMVM[8] * g.u.lightPos[2] + 
                                        invMVM[12];
        transLight[1] = invMVM[1] * g.u.lightPos[0] + 
                                        invMVM[5] * g.u.lightPos[1] + 
                                        invMVM[9] * g.u.lightPos[2] + 
                                        invMVM[13];
        transLight[2] = invMVM[2] * g.u.lightPos[0] + 
                                        invMVM[6] * g.u.lightPos[1] + 
                                        invMVM[10] * g.u.lightPos[2] + 
                                        invMVM[14];
        transLight[3] = 1.0f;

        glPushMatrix();
        glLoadIdentity();
        glLightfv(GL_LIGHT0, GL_POSITION, transLight);
        glPopMatrix();
}

void raycastBackground(GLenum destinationBuffer)
{
        GLdouble modelview[16], projection[16], vertizes[4][3];
        GLint viewport[4];
        GLdouble x, y, z;
        int i;

        float texelHeight = 1.0/(float)g.windowHeight;
        float texelWidth = 1.0/(float)g.windowWidth;

        activateARBProg(g.currentFragProg);

        glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
        glEnable(GL_CULL_FACE);

        glProgramLocalParameter4fARB(GL_FRAGMENT_PROGRAM_ARB, 1, 
                                                                 g.u.texCoordScale, 0.0, 0.0, 0.0);

        glProgramLocalParameter4fvARB(GL_FRAGMENT_PROGRAM_ARB, 2, 
                                                                  g.center);

        glProgramLocalParameter4fARB(GL_FRAGMENT_PROGRAM_ARB, 3,
                                                                 g.extents[0] * g.scaleFactors[0], 
                                                                 g.extents[1] * g.scaleFactors[1], 
                                                                 g.extents[2] * g.scaleFactors[2], 
                                                                 0.0);

        glProgramLocalParameter4fvARB(GL_FRAGMENT_PROGRAM_ARB, 5, 
                                                                  g.scaleFactors);
        
        glProgramLocalParameter4fvARB(GL_FRAGMENT_PROGRAM_ARB, 6,
                                                                  g.scaleFactorsInv);

        glProgramLocalParameter4fARB(GL_FRAGMENT_PROGRAM_ARB, 7, 
                                                                 (1.0 - texelWidth)/(float)g.windowWidth, 
                                                                 (1.0 - texelHeight)/(float)g.windowHeight, 
                                                                 0.5 * texelWidth, 0.5 * texelHeight);

        glGetIntegerv(GL_VIEWPORT, viewport);
        glGetDoublev(GL_MODELVIEW_MATRIX, modelview);
        glGetDoublev(GL_PROJECTION_MATRIX, projection);

        for (i = 0; i < 4; i++) {
                x = ((i&1)>>0) * viewport[2];
                y = ((i&2)>>1) * viewport[3];
                z = .5;

                gluUnProject(x, y, z, modelview, projection, viewport,
                                         &vertizes[i][0], &vertizes[i][1], &vertizes[i][2]);
        }

        glBegin(GL_QUADS);
                vertexv(vertizes[0]);
                vertexv(vertizes[1]);
                vertexv(vertizes[3]);
                vertexv(vertizes[2]);
        glEnd();

        deactivateARBProg(g.currentFragProg);

        glDrawBuffer(GL_BACK);
}

void renderVolume(void)
{
        float texelHeight = 1.0/(float)g.windowHeight;
        float texelWidth = 1.0/(float)g.windowWidth;

        glProgramLocalParameter4fARB(GL_FRAGMENT_PROGRAM_ARB, 1,
                                                                 g.u.texCoordScale, g.u.numIterations, 
                                                                 g.u.isoValue, 0);

        glProgramLocalParameter4fvARB(GL_FRAGMENT_PROGRAM_ARB, 2, 
                                                                  g.center);

        glProgramLocalParameter4fARB(GL_FRAGMENT_PROGRAM_ARB, 3,
                                                                 g.extents[0] * g.scaleFactors[0], 
                                                                 g.extents[1] * g.scaleFactors[1], 
                                                                 g.extents[2] * g.scaleFactors[2], 
                                                                 0.0);

        glProgramLocalParameter4fARB(GL_FRAGMENT_PROGRAM_ARB, 4,
                                                                 g.u.scatteringScale, g.u.clipIsoValue, 
                                                                 0.0, 0.0);

        glProgramLocalParameter4fvARB(GL_FRAGMENT_PROGRAM_ARB, 5, 
                                                                  g.scaleFactors);

        glProgramLocalParameter4fvARB(GL_FRAGMENT_PROGRAM_ARB, 6,
                                                                  g.scaleFactorsInv);

        glProgramLocalParameter4fARB(GL_FRAGMENT_PROGRAM_ARB, 0, g.u.stepSize,
                                                                 g.u.gradScale, g.u.gradOffset, -5.0);

        glProgramLocalParameter4fARB(GL_FRAGMENT_PROGRAM_ARB, 7, 
                                                                 (1.0 - texelWidth)/(float)g.windowWidth, 
                                                                 (1.0 - texelHeight)/(float)g.windowHeight, 
                                                                 0.5 * texelWidth, 0.5 * texelHeight);
        drawQuads();
}

void raycastVolume()
{
        activateARBProg(g.currentFragProg);

        glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
    glDisable(GL_DEPTH_TEST);
        
        glEnable(GL_CULL_FACE);
        
        renderVolume();

        deactivateARBProg(g.currentFragProg);
}

void display(void)
{
        static int frameCount = 0;
        static double t0, t1;
        Vector3 axis;
        float angle;

        if (!g.initialized) {
                glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
                glViewport(0, 0, g.windowWidth, g.windowHeight);


                glMatrixMode(GL_PROJECTION);
                glLoadIdentity();
#if 1
                gluPerspective(FOVY, (float)g.windowWidth/(float)g.windowHeight, 
                                           NEAR_CLIP, FAR_CLIP);
#else
                gluPerspective(FOVY, (float)g.windowWidth/(float)g.windowHeight, 
                                                g.u.camZ - 1.0, g.u.camZ + 1.0);
#endif
                resizeTE();
                CHECK_FOR_OGL_ERROR();
                g.initialized = 1;
        }

        CHECK_FOR_OGL_ERROR();

        if (frameCount == 0) {
                t0 = timer();
        }

        glMatrixMode(GL_PROJECTION);
        glLoadIdentity();
#if 1
        gluPerspective(FOVY, (float)g.windowWidth/(float)g.windowHeight, 
                                   NEAR_CLIP, FAR_CLIP);
#else
        gluPerspective(FOVY, (float)g.windowWidth/(float)g.windowHeight, 
                                        g.u.camZ - 1.0, g.u.camZ + 1.0);
#endif
    glMatrixMode(GL_MODELVIEW);
        glLoadIdentity();

        glTranslatef(g.u.translate[0], g.u.translate[1], g.u.translate[2]);
        glTranslatef(0.0, 0.0, -g.u.camZ);
        glRotatef(g.animatedAngle, 0.0, 1.0, 0.0);
        Quaternion_getAngleAxis(g.u.camRot, &angle, &axis);
        glRotatef(angle * 180.0 / M_PI, axis.x, axis.y, axis.z);

        glTranslatef(-g.center[0], -g.center[1], -g.center[2]);

        setLight();

        CHECK_FOR_OGL_ERROR();

    glDisable(GL_DEPTH_TEST);
    glClear(GL_DEPTH_BUFFER_BIT);

    //@carol add
    raycastBackground(GL_BACK);
    //raycastVolume();
    GLSL_raycastVolume();

    if (g.u.wireframe) {
                drawWireframe();
        }

        if (g.u.drawLight) {
                drawLight();
        }

        /* Draw transfer function editor */
        drawTE();

        glutSwapBuffers();

        frameCount++;

        if (frameCount == MAX_FRAMES) {
                float fps;

                t1 = timer();
                fps = (float)MAX_FRAMES/(t1 - t0) * 1000.0;

                frameCount = 0;
                
#if CONTINOUS_FPS == 0
                if (fps < g.minFps) {
                        g.minFps = fps;
                }

                if (fps > g.maxFps) {
                        g.maxFps = fps;
                }

                g.fpsSum += fps;
                g.animatedFrames++;
#else
                fprintf(stdout, "%f fps\n", fps);
#endif
        }
}

void idle(void)
{

        if (g.animated) {
                g.animatedAngle += 1.0;
                if (g.animatedAngle > 359.5) {
#if CONTINOUS_FPS == 0          
                        fprintf(stdout, "minFps: %f\n", g.minFps);
                        fprintf(stdout, "avgFps: %f\n", g.fpsSum/(float)g.animatedFrames);
                        fprintf(stdout, "maxFps: %f\n\n", g.maxFps);
                        g.animatedFrames = 0;
                        g.minFps = FLT_MAX;
                        g.fpsSum = 0.0;
                        g.maxFps = 0.0;
#endif                  
                        g.animatedAngle = 0.0;
                }
        }
        glutPostRedisplay();
}

void updatePreIntTable(void)
{
        updateTransferFunction();

        calcPreIntTable(NUM_TE_ENTRIES, g.u.sliceThickness, g.tfFunc, 
                                        g.preIntTable);

        loadTETexture();
}

void updateOptDensityTexture(void)
{
        static int initialized = 0;
        static GLuint texId;

        if (! initialized) {
                glGenTextures(1, &texId);
                addTexture(texId, GL_TEXTURE_1D, TEXTURE_OPTICALDENSITY);
                initialized = 1;
        }

        memcpy(g.optDensity, g.dataTE[4], NUM_TE_ENTRIES);

    glBindTexture(GL_TEXTURE_1D, texId);
        glTexImage1D(GL_TEXTURE_1D, 0, GL_LUMINANCE, NUM_TE_ENTRIES, 
                                 0, GL_LUMINANCE, GL_UNSIGNED_BYTE, g.optDensity);

    glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
    glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

    glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
}

void initScatteringTexture(void)
{
        unsigned char line[100];
        unsigned char *data;
        int numEntries;
        GLuint texId;
        FILE *fp;

        if (! (fp = fopen("colortables/colortable_wax.tab", "rb"))) {
                fprintf(stderr, "opening scattering colortable failed\n");
                exit(1);
        }

        if (! fgets((char *)line, sizeof(line), fp)) {
                fprintf(stderr, "invalid colortable\n");
                exit(1);
        }

        if (sscanf((const char *)line, "%i\n", &numEntries) != 1) {
                fprintf(stderr, "invalid colortable\n");
                exit(1);
        }

        if (! (data = (unsigned char *)malloc(numEntries * 3))) {
                fprintf(stderr, "not enough memory for scattering colortable\n");
                exit(1);
        }

        if (fread(data, numEntries * 3, 1, fp) != 1) {
                fprintf(stderr, "reading scattering colortable failed\n");
                exit(1);
        }

        fclose(fp);

        glGenTextures(1, &texId);
    glBindTexture(GL_TEXTURE_1D, texId);
        glTexImage1D(GL_TEXTURE_1D, 0, GL_RGB, numEntries, 0, GL_RGB,
                                 GL_UNSIGNED_BYTE, data);

        addTexture(texId, GL_TEXTURE_1D, TEXTURE_SCATTERING);
        
    glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

    glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
}

void printARBProgs(int numProgs, ARBProgram *progs)
{
        ARBProgram *prog;
        int i;

        for (i = 0; i < numProgs && i < 12; i++) {
                prog = &progs[i];
                fprintf(stderr, "\tF%i: %s (%i bytes)\n", 
                                i + 1, prog->filename, prog->size);

#if 0
                for (j = 0; j < prog->numTextureObjects; j++) {
                        fprintf(stderr, "\t\ttexture unit %i: %s\n", 
                                        prog->textureObjects[j].texUnit,
                                        prog->textureObjects[j].texture.name);
                }
#endif
        }

        fprintf(stderr, "\n");
}

void loadVolumeShaders()
{
        loadARBProgs(VOLUME_SHADERS_DIR, FRAG_PROG_EXT,
                                 g.numTextures, g.textures,
                                 &g.numFragProgs, &g.fragProgs);
        if (g.numFragProgs <= 0) {
                fprintf(stderr, "No shaders found\n");
                exit(1);
        }
        printARBProgs(g.numFragProgs, g.fragProgs);
        g.currentFragProg = &g.fragProgs[0];
}

void freeVolumeShaders(void)
{
        
        ARBProgram *prog;
        int i;
        
        for (i = 0; i < g.numFragProgs; i++) {
                prog = &g.fragProgs[i];
                //@carol add
                glDeleteProgram(prog->id);
                //glDeleteProgramsARB(1, &prog->id);
                free(prog->filename);
                free(prog->prog);
        }
        free(g.fragProgs);
        
}

char *getSettingsFilename(void)
{
        char *filename;

        if (! (filename = (char *)malloc(strlen(g.basename) + 1 +
                                strlen(g.currentFragProg->filename) - strlen(FRAG_PROG_EXT) +
                                + strlen(SETTINGS_EXT) + 1))) {
                fprintf(stderr, "not enough memory for filename\n");
                exit(1);
        }

        strcpy(filename, g.basename);
        strcat(filename, "_");
        strcat(filename, g.currentFragProg->filename);
        *strrchr(filename, '.') = 0;
        strcat(filename, SETTINGS_EXT);

        return filename;
}

void saveSettings(void)
{
        char *filename;
        FILE *fp;
        int i;

        filename = getSettingsFilename();

        fprintf(stderr, "saving settings to file %s\n", filename);

        if (! (fp = fopen(filename, "wb"))) {
                perror("cannot open settings file for writing");
                return;
        }

        if (fwrite(&g.u, sizeof(UserSettings), 1, fp) != 1) {
                fprintf(stderr, "saving user settings failed\n");
                exit(1);
        }

        for (i = 0; i < NUM_TE; i++) {
                if (fwrite(g.dataTE[i], NUM_TE_ENTRIES * sizeof(DataTypeTE), 
                                   1, fp) != 1) {
                        fprintf(stderr, "saving transfer functions failed\n");
                        exit(1);
                }
        }

        fclose(fp);

        free(filename);
}

void loadSettings(void)
{
        char *filename;
        FILE *fp;
        int i;

        filename = getSettingsFilename();

        fprintf(stderr, "loading settings from file %s\n", filename);

        if (! (fp = fopen(filename, "rb"))) {
                perror("cannot open settings file for reading");
                /* Do not exit if there are any problems */
        } else {
                if (fread(&g.u, sizeof(UserSettings), 1, fp) != 1) {
                        fprintf(stderr, "loading user settings failed\n");
                        exit(1);
                }

                for (i = 0; i < NUM_TE; i++) {
                        if (fread(g.dataTE[i], NUM_TE_ENTRIES * sizeof(DataTypeTE), 
                                          1, fp) != 1) {
                                fprintf(stderr, "loading transfer functions failed\n");
                                exit(1);
                        }
                }
                updateOptDensityTexture();
                loadBackgroundTexture();
                updatePreIntTable();

                fclose(fp);
        }

        free(filename);

        glutPostRedisplay();
}

void printSettings(void)
{
        Vector3 axis;
        float angle;

        printf("step size ........................: %g\n", g.u.stepSize);
        printf("slice thickness ..................: %g\n", g.u.sliceThickness);
        printf("gradient scale ...................: %g\n", g.u.gradScale);
        printf("gradient offset ..................: %g\n", g.u.gradOffset);
        printf("texture coordinate scale .........: %g\n", g.u.texCoordScale);
        printf("isovalue of volume ...............: %g\n", g.u.isoValue);
        printf("isovalue of clipping volume ......: %g\n", g.u.clipIsoValue);
        printf("isovalue step ....................: %g\n", g.isoValueStep);
        printf("number of iterations per loop ....: %i\n", g.u.numIterations);
        printf("wireframe mode ...................: %s\n", 
                   g.u.wireframe ? "yes" : "no");
        printf("light show mode ..................: %s\n", 
                   g.u.drawLight ? "yes" : "no");
        printf("scattering scale .................: %g\n", g.u.scatteringScale);
        printf("light source position ............: %.2g %.2g %.2g\n", 
                   g.u.lightPos[0], g.u.lightPos[1], g.u.lightPos[2]);
        printf("object translation ...............: %.2g %.2g %.2g\n", 
                   g.u.translate[0], g.u.translate[1], g.u.translate[2]);
        printf("camera z-value ...................: %g\n", g.u.camZ);
        Quaternion_getAngleAxis(g.u.camRot, &angle, &axis);
        printf("camera rotation ..................: (%.2g, %.2g, %.2g) "
                   "at %g degrees\n", axis.x, axis.y, axis.z, angle);
        printf("background image mode ............: %s\n", 
                   g.u.backgroundImage ? "yes" : "no");
        printf("background gray value ............: %i\n", g.u.backgroundGrayVal);
        printf("\n");
}

void key(unsigned char k, int x, int y)
{
    int already_handled;
    already_handled = keyTE(k, x, y);

    if (already_handled) {
                updateOptDensityTexture();
        } else {
        switch (k) {
            case 27:
            case 'q':
                exit(0);
                break;
                        case '!':
                                g.initialized = 0;
                                break;
            case 'w':
                g.u.wireframe = ! g.u.wireframe;
                break;
            case 't':
                toggleHidden();
                break;
                        case 'p':
                                updatePreIntTable();
                                break;
                        case '+':
                                g.u.sliceThickness += SLICE_STEP;
                                if (g.u.sliceThickness > MAX_SLICE_THICKNESS) {
                                        g.u.sliceThickness = MAX_SLICE_THICKNESS;
                                }
                                break;
                        case '-':
                                g.u.sliceThickness -= SLICE_STEP;
                                if (g.u.sliceThickness < MIN_SLICE_THICKNESS) {
                                        g.u.sliceThickness = MIN_SLICE_THICKNESS;
                                }
                                break;
                        case '>':
                                g.u.stepSize += STEPSIZE_STEP;
                                if (g.u.stepSize > MAX_STEPSIZE) {
                                        g.u.stepSize = MAX_STEPSIZE;
                                }
                                break;
                        case '<':
                                g.u.stepSize -= STEPSIZE_STEP;
                                if (g.u.stepSize < MIN_STEPSIZE) {
                                        g.u.stepSize = MIN_STEPSIZE;
                                }
                                break;
                        case ']':
                                g.u.stepSize *= 2.0;
                                if (g.u.stepSize > MAX_STEPSIZE) {
                                        g.u.stepSize = MAX_STEPSIZE;
                                }
                                break;
                        case '[':
                                g.u.stepSize /= 2.0;
                                if (g.u.stepSize < MIN_STEPSIZE) {
                                        g.u.stepSize = MIN_STEPSIZE;
                                }
                                break;
                        case 's':
                                g.u.gradScale += GRAD_SCALE_STEP;
                                fprintf(stderr, "gradient: scale %f, offset %f\n", 
                                                g.u.gradScale, g.u.gradOffset);
                                break;
                        case 'x':
                                g.u.gradScale -= GRAD_SCALE_STEP;
                                break;
                        case 'd':
                                g.u.gradOffset += GRAD_OFFSET_STEP;
                                break;
                        case 'c':
                                g.u.gradOffset -= GRAD_OFFSET_STEP;
                                break;
                        case 'f':
                                g.u.texCoordScale += TEX_COORD_SCALE_STEP;
                                break;
                        case 'v':
                                g.u.texCoordScale -= TEX_COORD_SCALE_STEP;
                                break;
                        case 'l':
                                g.u.drawLight = ! g.u.drawLight;
                                break;
                        case 'j':
                                if (glutGetModifiers() & GLUT_ACTIVE_ALT) {
                                        g.isoValueStep *= 2.0;
                                        if (g.isoValueStep > 1.0) {
                                                g.isoValueStep = 1.0;
                                        }
                                } else {
                                        g.u.isoValue += g.isoValueStep;
                                        if (g.u.isoValue > 1.0) {
                                                g.u.isoValue = 1.0;
                                        }
                                }                                       
                                break;
                        case 'm':
                                if (glutGetModifiers() & GLUT_ACTIVE_ALT) {
                                        g.isoValueStep /= 2.0;                                  
                                } else {
                                        g.u.isoValue -= g.isoValueStep;
                                        if (g.u.isoValue < 0.0) {
                                                g.u.isoValue = 0.0;
                                        } 
                                } 
                                break;
                        case 'J':
                                g.u.clipIsoValue += g.isoValueStep;
                                if (g.u.clipIsoValue > 1.0) {
                                        g.u.clipIsoValue = 1.0;
                                }
                                break;
                        case 'M':
                                g.u.clipIsoValue -= g.isoValueStep;
                                if (g.u.clipIsoValue < 0.0) {
                                        g.u.clipIsoValue = 0.0;
                                }
                                break;
                        case 'r':
                                freeVolumeShaders();
                                loadVolumeShaders();
                                break;
                        case '.':
                                g.u.scatteringScale += SCATTERING_STEP;
                                break;
                        case ',':
                                g.u.scatteringScale -= SCATTERING_STEP;
                                if (g.u.scatteringScale < 0.0) {
                                        g.u.scatteringScale = 0.0;
                                }
                                break;
                        case '1':
                        case '2':
                        case '3':
                        case '4':
                        case '5':
                        case '6': {
                                int coord = (k - '1')/2;
                                if ((k - '1')%2) {
                                        g.u.lightPos[coord] += LIGHT_POS_STEP;
                                } else {
                                        g.u.lightPos[coord] -= LIGHT_POS_STEP;
                                }
                                break;
                        }
                        case 'S':
                                saveSettings();
                                break;
                        case 'L':
                                loadSettings();
                                break;
                        case 'k':
                                g.u.backgroundImage = ! g.u.backgroundImage;
                                loadBackgroundTexture();
                                break;
                        case '}':
                                if (g.u.backgroundGrayVal < 255) {
                                        g.u.backgroundGrayVal++;
                                        loadBackgroundTexture();
                                }
                                break;
                        case '{':
                                if (g.u.backgroundGrayVal > 0) {
                                        g.u.backgroundGrayVal--;
                                        loadBackgroundTexture();
                                }
                                break;
                        case 'i':
                                printSettings();
                                break;
                        case ' ':
                                g.animated = ! g.animated;
                                g.animatedFrames = 0;
                                g.minFps = FLT_MAX;
                                g.fpsSum = 0.0;
                                g.maxFps = 0.0;
                                g.animatedAngle = 0.0;
                                break;
        }
    }

    glutPostRedisplay();
}

void special(int key, int x, int y)
{
        int num;

        if (key >= GLUT_KEY_F1 && key <= GLUT_KEY_F12) {
                num = key - GLUT_KEY_F1;

                if (num < g.numFragProgs) {
                        g.currentFragProg = &g.fragProgs[num];
                        fprintf(stderr, "activated fragment program '%s'\n",
                                        g.currentFragProg->filename);

#if 0
                        loadSettings();
#endif
                        glutPostRedisplay();
                }
        }
}

void mouseMouseInteract(int button, int state, int x, int y)
{
        int modifiers = glutGetModifiers();

        g.mousePosLast[0] = x;
        g.mousePosLast[1] = y;

        
        switch (button) {
                case GLUT_LEFT_BUTTON:
                        if (modifiers & GLUT_ACTIVE_CTRL) {
                                g.mouseMode = MOUSE_MOVE_LIGHT_XY;
                        }
                        else {
                                g.mouseMode = MOUSE_ROTATE;
                        }
                        break;
                case GLUT_MIDDLE_BUTTON:
                        if (modifiers & GLUT_ACTIVE_CTRL) {
                        }
                        else {
                                g.mouseMode = MOUSE_TRANSLATE;
                        }
                        break;  
                case GLUT_RIGHT_BUTTON:
                        if (modifiers & GLUT_ACTIVE_CTRL) {
                                g.mouseMode = MOUSE_MOVE_LIGHT_Z;
                        }
            else {
                                g.mouseMode = MOUSE_DOLLY;
                        }
                        break;
                default:
                        break;
        }
}

void motionMouseInteract(int x, int y)
{
        Vector3 mouseVector, view, rotAxis;
        Quaternion rot;

        float dx = MOUSE_SCALE * (x - g.mousePosLast[0]) / (float)g.windowWidth;
        float dy = MOUSE_SCALE * (g.mousePosLast[1] - y) / (float)g.windowHeight;

        mouseVector.x = dx;
        mouseVector.y = dy;
        mouseVector.z = 0.0;

        view.x = 0.0;
        view.y = 0.0;
        view.z = -1.0;

        rotAxis = Vector3_cross(mouseVector, view);
        Vector3_normalize(&rotAxis);

        switch (g.mouseMode) {
                case MOUSE_DOLLY:
                        g.u.camZ += dy;
                        break;
                case MOUSE_ROTATE:
                        rot = Quaternion_fromAngleAxis(SQR(dx) + SQR(dy), rotAxis);
                        g.u.camRot = Quaternion_mult(rot, g.u.camRot);
                        Quaternion_normalize(&g.u.camRot);
                        break;
                case MOUSE_TRANSLATE:
                        g.u.translate[0] += mouseVector.x;              
                        g.u.translate[1] += mouseVector.y;              
                        g.u.translate[2] += mouseVector.z;              
                        break;
                case MOUSE_MOVE_LIGHT_XY:
                        g.u.lightPos[0] += dx;
                        g.u.lightPos[1] += dy;
                        fprintf(stderr, "lightPos= %f, %f, %f\n", g.u.lightPos[0], g.u.lightPos[1], g.u.lightPos[2]);
                        break;
                case MOUSE_MOVE_LIGHT_Z:
                        g.u.lightPos[2] -= dy;
                        fprintf(stderr, "lightPos= %f, %f, %f\n", g.u.lightPos[0], g.u.lightPos[1], g.u.lightPos[2]);
                        break;
                default:
                        break;
        }

        g.mousePosLast[0] = x;
        g.mousePosLast[1] = y;
}

void mouse(int button, int state, int x, int y)
{
    int already_handled;

    already_handled = mouseTE(x, y);

    if (! already_handled)  {
                mouseMouseInteract(button, state, x, y);
        }

    glutPostRedisplay();
}

void motion(int x, int y)
{
    int already_handled;

    already_handled = motionTE(x, y);

    if (already_handled) {
                updateOptDensityTexture();
        } else {
                motionMouseInteract(x, y);
        }

    glutPostRedisplay();
}

void init(char *volfilename, char *clipfilename)
{
        g.initialized = 0;
        g.numTextures = 0;
        g.textures = NULL;
        g.mouseMode = MOUSE_ROTATE;    
    g.u.camZ = 2.0;
    g.u.camRot.x = g.u.camRot.y = g.u.camRot.z = 0.0; 
        g.u.camRot.w = 1.0;
        g.windowWidth = WINDOW_WIDTH; 
        g.windowHeight = WINDOW_HEIGHT;
        g.u.translate[0] = 0.0;
        g.u.translate[1] = 0.0;
        g.u.translate[2] = 0.0;
        g.u.wireframe = 0;
        g.u.drawLight = 0;
        g.u.stepSize = 1.0/128.0;/*.02;*/ /*1.0/512.0;*/
        g.u.sliceThickness = 1.0;
        g.u.gradOffset = .9;
        g.u.gradScale = .2;
        g.u.texCoordScale = .45;
        /* max. # of iterations supported for fragment programs on NV40 */
        g.u.numIterations = 255; 
        g.u.isoValue = 0.5;
        g.isoValueStep = ISO_VALUE_STEP;
        g.u.clipIsoValue = 0.5;
        g.u.scatteringScale = 6.0;
        g.u.lightPos[0] = 2.0;
        g.u.lightPos[1] = 2.0;
        g.u.lightPos[2] = 2.0;
        g.u.lightPos[3] = 1.0;
        g.u.backgroundImage = 0;
        g.u.backgroundGrayVal = 127;
        g.animatedAngle = 0.0;
        g.animated = 0;
        g.minFps = FLT_MAX;
        g.fpsSum = 0.0;
        g.maxFps = 0.0;

        if (! (g.basename = strdup(volfilename))) {
                fprintf(stderr, "cannot duplicate filename;\n");
        }
        if (endsWith(g.basename, VOL_FILE_EXT)) {
                *strrchr(g.basename, '.') = 0;
        } 
        fprintf(stderr, "%s\n", g.basename);

        if (! (g.preIntTable = (unsigned char *)malloc(NUM_TE_ENTRIES *
                                                                                                   NUM_TE_ENTRIES * 4))) {
                fprintf(stderr, "not enough memory for preintegrated table\n");
                exit(1);
        }

        if (! (g.tfFunc = (unsigned char *)malloc(4 * NUM_TE_ENTRIES))) {
                fprintf(stderr, "not enough memory for transfer function data\n");
                exit(1);
        }

        if (! (g.optDensity = (unsigned char *)malloc(NUM_TE_ENTRIES))) {
                fprintf(stderr, "not enough memory for optical density data\n");
                exit(1);
        }
}

void initOpenGL(char *volfilename, char *clipfilename)
{
        GLfloat lightSpecular[] = {1.0, 1.0, 1.0, 1.0};
        GLfloat lightAmbient[] = {0.3, 0.3, 0.3, 1.0};
        GLfloat lightDiffuse[] = {0.7, 0.7, 0.7, 1.0};

#if 0
        fprintf(stdout, "%s\n", glGetString(GL_EXTENSIONS));
#endif
    //@carol delete
    /*if (!glh_init_extensions(
          "GL_EXT_texture3D "
          "GL_ARB_fragment_program "
          "GL_ARB_multitexture "
                  "GL_ARB_vertex_program "
          "GL_NV_occlusion_query "
                  "GL_NV_fragment_program2 "
                  "GL_NV_vertex_program3 "
        ))
    {
        printf("Unable to load the following extension(s): %s\n\nExiting...\n", 
                glh_get_unsupported_extensions());
        exit(-1);
    }*/

    initTE();
        g.dataTE = getDataTE();
    loadVolumeTexture(volfilename);
        loadClipTexture(clipfilename);
        loadBackgroundTexture();
        loadSphereMapTexture();
        updateOptDensityTexture();
        initScatteringTexture();
        updatePreIntTable();

  
        loadVolumeShaders();
        GLSL_init();

    glClearColor(1.0f, 0.0f, 0.0f, 0.0f);
    glDisable(GL_DEPTH_TEST);

        glLightfv(GL_LIGHT0, GL_AMBIENT, lightAmbient);
        glLightfv(GL_LIGHT0, GL_DIFFUSE, lightDiffuse);
        glLightfv(GL_LIGHT0, GL_SPECULAR, lightSpecular);
        
        CHECK_FOR_OGL_ERROR();
}

void resize(int w, int h)
{
    if (h == 0) {
        return;
        }

        g.windowWidth = w;
        g.windowHeight = h;
        
        g.initialized = 0;

    glutPostRedisplay();
}

int main(int argc, char *argv[])
{
        if ( argc != 2 && argc != 3) {
                fprintf(stderr, "usage: spvolren <volfilename.dat>\
                                                 [<clipfilename.dat>]\n");
                exit(1);
        }

        if (argc == 2) {
                init(argv[1], NULL);
        } else {
                init(argv[1], argv[2]);
        }

    glutInit(&argc, argv);
        glutInitWindowPosition(
                (720 - WINDOW_WIDTH)/2, 
                (576 - WINDOW_HEIGHT)/2);
    glutInitWindowSize(WINDOW_WIDTH, WINDOW_HEIGHT);
    glutInitDisplayMode(GLUT_DOUBLE | GLUT_DEPTH | GLUT_RGBA | GLUT_ALPHA);
    glutCreateWindow("Single Pass Volume Raycasting Demo");

        if (argc == 2) {
        initOpenGL(argv[1], NULL);
        } else { 
                initOpenGL(argv[1], argv[2]);
        }

        glutDisplayFunc(display);
    glutMotionFunc(motion);
    glutMouseFunc(mouse);
    glutIdleFunc(idle);
    glutKeyboardFunc(key);
        glutSpecialFunc(special);
    glutReshapeFunc(resize);

        glutMainLoop();

    return 0;
}

---