/***************************************************************************
 * vrend.cpp
 * Note: Images are oriented with pixel (0,0) in the top left corner.
 ***************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <iostream>
#include <GL/glut.h>
#include <limits.h>    /* This is where SHRT_MAX is defined */
#include <time.h>

#define TRUE  1
#define FALSE 0

#ifndef DBL_MAX
#define DBL_MAX         1.7976931348623157E+308
#endif
//#define PIE 3.14159265358979323846
#define PIE 3.141592653589
#define TWO_PIE (2*PIE)
//6.2831853072
#define DBL_TOL .0000001


#define FOV 1.047 //PI/3.0
#define XDIM 256
#define YDIM 256
#define BPP 3

#define NONE   0
#define ROTATE 1
#define ZOOM   2
#define CUT    3
#define TIME   4

#define FLAT    0
#define DEPTH   1
#define COLOR   2
#define DENSITY 3
#define NORMAL  4
#define NICE    5

#define ROOT2 1.41421356
#define ROOT3 1.73205080
#define INV_ROOT2 0.70710678
#define INV_ROOT3 0.57735027

class Volume;

//GLOBAL VARAIBLES
Volume *volume = NULL;
int g_delay;
int mode = NONE;
double p1x=0, p1y=0, p2x=0, p2y=0;
int g_slice_id = 0;
int g_x_prev = 0;
int g_y_prev = 0;
double g_scale = 1.0;
double g_frames_count = 0.0;
clock_t  g_clock_start = 0;


//Function prototypes
void display();
void reshape(int w, int h);

class Volume{
  //Description of volume
  //Faced normals
  //

  public:
    int needs_update;
    double eyex, eyey, eyez;
    double atx, aty, atz;
    double upx, upy, upz;
    double cx, cy, cz;
    double axis_x, axis_y, axis_z;
    double d_width, d_height, d_depth;
    int width, height, depth, im_width, im_height;
    double neg_dist;//Note this is a negative number.
    double depth_dist, depth_scale;
    double zmult;
    unsigned char *data;
    unsigned char *vr_data;
    unsigned char *im_data;
    int im_dim;
    unsigned char bgr, bgg, bgb;
    unsigned char render_type;
    unsigned char threshhold;
    unsigned char max_threshhold;
    int quality_rend;
    int show_skin;
    int curr_pic;

    //Cached values for the ray caster
    double *dirx;
    double *diry;
    double *dirz;
    Volume(int w, int h, int d, int xres, int yres, double z_multiplier){
      double mag;

      zmult = z_multiplier;
      im_width = xres;
      im_height = yres;
      im_dim = xres*yres;

      width = w;
      height = h;
      depth = d;

      d_width = (double)w;
      d_height = (double)h;
      d_depth = (double)d;
      d_depth /= z_multiplier;
      
      data = new unsigned char [width*height*depth];
      vr_data = new unsigned char [width*height*depth];

      im_data = new unsigned char[im_width*im_height*BPP];
      dirx = new double[im_width*im_height];
      diry = new double[im_width*im_height];
      dirz = new double[im_width*im_height];
      atx = cx = (double)width/2;
      aty = cy = (double)height/2;
      atz = cz = (double)depth/2;

      eyex = (double)width/2;
      eyey = (double)height/2;
      eyez = neg_dist = atz - (double)im_width/(2*tan(FOV*.5));

      upy = 1.0;
      upx = 0;
      upz = 0;

      //Start further away from the object
      atz -= width;
      eyez -= width;
      
      depth_scale =  ROOT3*(double)(width/3);
      depth_dist = fabs(eyez-cz) -  depth_scale;// ROOT3*(double)(width/2);
      if (depth_dist < 0.0) depth_dist = 0.0;
      
      //Set the bg color
      bgr = 0;
      bgg = 0;
      bgb = 0;
      
      axis_x = atx-eyex;
      axis_y = aty-eyey;
      axis_z = atz-eyez;
      mag = sqrt(axis_x*axis_x + axis_y*axis_y + axis_z*axis_z);
      if (mag < DBL_TOL) mag = 1.0;
      axis_x /= mag; axis_y /= mag; axis_z /= mag;

      curr_pic = 0;
      threshhold = 64;
      max_threshhold = 255;
      render_type = NICE;
      quality_rend = FALSE;
      needs_update = TRUE;
      show_skin = TRUE;
      calcDir();
      clearVrData(0);
    }
    ~Volume(){
      if (NULL != data){ delete [] data; }
      if (NULL != vr_data){ delete [] vr_data; }
      if (NULL != im_data){ delete [] im_data; }
      if (NULL != dirx){ delete [] dirx; }
      if (NULL != diry){ delete [] diry; }
      if (NULL != dirz){ delete [] dirz; }
    }
    void calcDir();
    void clearVrData(unsigned char val);
    void cut(double x1, double y1, double x2, double y2);
    void dilate();
    void erode();
    inline void forceUpdate(){
      needs_update = TRUE;
    }
    int read(char *filename);
    inline void rotate(double theta, double gamma);
    inline void setIntersectColor(int loc, double rayx, double rayy, double rayz);
    inline void setDepthColor(int x0, int y0, int z0, int x1, int y1, int z1, int loc, double idist);
    inline void setDataColor(int x0, int y0, int z0, int x1, int y1, int z1, int loc);
    inline void setDensityColor(int x0, int y0, int z0, int x1, int y1, int z1, int loc, double idist);
    inline void setNiceColor(int x0, int y0, int z0, int x1, int y1, int z1, int loc);
    inline void setNormalColor(int x0, int y0, int z0, int x1, int y1, int z1, int loc);
    inline void setRenderType(unsigned char rt){
      render_type = rt;
      needs_update = TRUE;
    }
    inline void setThresh(unsigned char t, unsigned char max_t){
      threshhold = t;
      max_threshhold = max_t;
      needs_update = TRUE;
    }
    void thresh(unsigned char t);
    inline void toggleQualityRend(){
      quality_rend = !quality_rend;
      needs_update = TRUE;
    }
    inline void toggleShowSkin(){
      show_skin = !show_skin;
      needs_update = TRUE;
    }
    void updateDisp();
    void writeImage();
    inline void zoom(double theta);
};

//VOLUME FUNCTIONS
//VOLUME FUNCTIONS
//VOLUME FUNCTIONS
//VOLUME FUNCTIONS

//Perform a 3D grayscale dilation
void Volume::dilate(){
  register int y_off = width;
  register int z_off = width*height;
  int loc, x_end, y_end, z_end;
  unsigned char max;
  unsigned char val;
  unsigned char *temp;

  z_end = height*width*depth - z_off;
  loc = z_off;//start point for loc
  while (loc<z_end){
    y_end = loc + z_off - y_off;
    loc += y_off;
    while (loc<y_end){
      x_end = loc + y_off - 1;
      loc ++;

      while (loc<x_end){
        //Temp for test.
        //data[loc] = 255;

          //Front (in x)
        max = data[loc + 1 + y_off - z_off];
        val = data[loc + 1 + y_off]; if (val > max) max = val;
        val = data[loc + 1 + y_off + z_off]; if (val > max) max = val;

        val = data[loc + 1 - z_off]; if (val > max) max = val;
        val = data[loc + 1 ]; if (val > max) max = val;
        val = data[loc + 1 + z_off]; if (val > max) max = val;

        val = data[loc + 1 - y_off - z_off]; if (val > max) max = val;
        val = data[loc + 1 - y_off]; if (val > max) max = val;
        val = data[loc + 1 - y_off + z_off]; if (val > max) max = val;


        //Middle
        val = data[loc + y_off - z_off]; if (val > max) max = val;
        val = data[loc + y_off]; if (val > max) max = val;
        val = data[loc + y_off + z_off]; if (val > max) max = val;

        val = data[loc - z_off]; if (val > max) max = val;
        val = data[loc ]; if (val > max) max = val; //original pixel location
        val = data[loc + z_off]; if (val > max) max = val;

        val = data[loc - y_off - z_off]; if (val > max) max = val;
        val = data[loc - y_off]; if (val > max) max = val;
        val = data[loc - y_off + z_off]; if (val > max) max = val;

        //Back
        val = data[loc - 1 + y_off - z_off]; if (val > max) max = val;
        val = data[loc - 1 + y_off]; if (val > max) max = val;
        val = data[loc - 1 + y_off + z_off]; if (val > max) max = val;

        val = data[loc - 1 - z_off]; if (val > max) max = val;
        val = data[loc - 1 ]; if (val > max) max = val;
        val = data[loc - 1 + z_off]; if (val > max) max = val;

        val = data[loc - 1 - y_off - z_off]; if (val > max) max = val;
        val = data[loc - 1 - y_off]; if (val > max) max = val;
        val = data[loc - 1 - y_off + z_off]; if (val > max) max = val;


        vr_data[loc] = max;

        loc++;
      }
      loc++;//put it at the end of the row
    }
    loc += y_off;//put it at the end of the page
  }
  //Swap the data pointers
  temp = data;
  data = vr_data;
  vr_data = temp;
  needs_update = TRUE;
}

void Volume::erode(){
  register int y_off = width;
  register int z_off = width*height;
  int loc, x_end, y_end, z_end;
  unsigned char min;
  unsigned char val;
  unsigned char *temp;

  z_end = height*width*depth - z_off;
  loc = z_off;//start point for loc
  while (loc<z_end){
    y_end = loc + z_off - y_off;
    loc += y_off;
    while (loc<y_end){
      x_end = loc + y_off - 1;
      loc ++;

      while (loc<x_end){
        //Temp for test.
        //data[loc] = 255;

          //Front (in x)
        min = data[loc + 1 + y_off - z_off];
        val = data[loc + 1 + y_off]; if (val < min) min = val;
        val = data[loc + 1 + y_off + z_off]; if (val < min) min = val;

        val = data[loc + 1 - z_off]; if (val < min) min = val;
        val = data[loc + 1 ]; if (val < min) min = val;
        val = data[loc + 1 + z_off]; if (val < min) min = val;

        val = data[loc + 1 - y_off - z_off]; if (val < min) min = val;
        val = data[loc + 1 - y_off]; if (val < min) min = val;
        val = data[loc + 1 - y_off + z_off]; if (val < min) min = val;


        //Middle
        val = data[loc + y_off - z_off]; if (val < min) min = val;
        val = data[loc + y_off]; if (val < min) min = val;
        val = data[loc + y_off + z_off]; if (val < min) min = val;

        val = data[loc - z_off]; if (val < min) min = val;
        val = data[loc ]; if (val < min) min = val; //original pixel location
        val = data[loc + z_off]; if (val < min) min = val;

        val = data[loc - y_off - z_off]; if (val < min) min = val;
        val = data[loc - y_off]; if (val < min) min = val;
        val = data[loc - y_off + z_off]; if (val < min) min = val;

        //Back
        val = data[loc - 1 + y_off - z_off]; if (val < min) min = val;
        val = data[loc - 1 + y_off]; if (val < min) min = val;
        val = data[loc - 1 + y_off + z_off]; if (val < min) min = val;

        val = data[loc - 1 - z_off]; if (val < min) min = val;
        val = data[loc - 1 ]; if (val < min) min = val;
        val = data[loc - 1 + z_off]; if (val < min) min = val;

        val = data[loc - 1 - y_off - z_off]; if (val < min) min = val;
        val = data[loc - 1 - y_off]; if (val < min) min = val;
        val = data[loc - 1 - y_off + z_off]; if (val < min) min = val;


        vr_data[loc] = min;

        loc++;
      }
      loc++;//put it at the end of the row
    }
    loc += y_off;//put it at the end of the page
  }
  //Swap the data pointers
  temp = data;
  data = vr_data;
  vr_data = temp;
  needs_update = TRUE;
}

void Volume::thresh(unsigned char t){
  int i, dim;

  dim = height*width*depth;
  for (i=0; i<dim; i++){
    if (data[i] < t)
      data[i] = 0;
  }
  needs_update = TRUE;
}

/*
 * Calculates all the directions for each ray
 */
void Volume::calcDir(){
  int i, j, iw, jh, im_loc, loc;
  double nx, ny, nz, ux, uy, uz, vx, vy, vz, mag;
  double dx, dy, dz;
  int halfwidth, halfheight;


  //Define new coordinate system based on eye and at points.
  nx = atx - eyex;
  ny = aty - eyey;
  nz = atz - eyez;
  mag = sqrt(nx*nx + ny*ny + nz*nz);
  nx/=mag;
  ny/=mag;
  nz/=mag;

  //cross normal with up to get u (horizontal)
  //Up = 0, 1, 0
  ux =   ny*upz - upy*nz;
  uy = - nx*upz + upx*nz;
  uz =   nx*upy - upx*ny;
  mag = sqrt(ux*ux + uy*uy + uz*uz);
  ux/=mag;
  uy/=mag;
  uz/=mag;

  //cross u with normal to get v (vertical)
  vx =   uy*nz - ny*uz;
  vy = - ux*nz + nx*uz;
  vz =   ux*ny - nx*uy;

  mag = sqrt(vx*vx + vy*vy + vz*vz);
  vx/=mag;
  vy/=mag;
  vz/=mag;


  halfwidth = im_width/2;
  halfheight = im_height/2;



  //Do ray casting for the volume data

  for (j=0; j<im_height; j++){
    jh = halfheight - j;
    loc = j*im_width;
    for (i=0; i<im_width; i++){
      iw = i - halfwidth;
      im_loc = i +loc;
      //Calculate the ray.
       
      dx = atx + iw*ux - jh*vx - eyex;
      dy = aty + iw*uy - jh*vy - eyey;
      dz = atz + iw*uz - jh*vz - eyez;
      mag = sqrt(dx*dx + dy*dy + dz*dz);
      dirx[im_loc] = dx / mag;
      diry[im_loc] = dy / mag;
      dirz[im_loc] = dz / mag;
    }
  }
  needs_update = TRUE;
}

void Volume::clearVrData(unsigned char val){
  int i, dim;

  dim = height*width*depth;
  for (i=0; i<dim; i++){
    vr_data[i] = val;
  }
  needs_update = TRUE;
}

//Lobs off the top half of the volume based on the plane created with this line.
void Volume::cut(double x1, double y1, double x2, double y2){
  int halfwidth, halfheight;
  double new_x1, new_y1, new_z1, new_x2, new_y2, new_z2;
  double ax, ay, az, bx, by, bz, cx, cy, cz;
  double nx, ny, nz, mag;
  double ux, uy, uz, vx, vy, vz;
  double d;
  int i, j, k, wh;
  

  //Take the cross product between these two points and the eye point.
  //Must transform the points into world space.
  halfwidth = im_width/2;
  halfheight = im_height/2;

  bx = atx - eyex;
  by = aty - eyey;
  bz = atz - eyez;
  mag = sqrt(bx*bx + by*by + bz*bz);
  if (mag < DBL_TOL) mag = 1.0;
  bx /= mag;
  by /= mag;
  bz /= mag;

  ux =   by*upz - upy*bz;
  uy = - bx*upz + upx*bz;
  uz =   bx*upy - upx*by;
  mag = sqrt(ux*ux + uy*uy + uz*uz);
  if (mag < DBL_TOL) mag = 1.0;
  ux/=mag;
  uy/=mag;
  uz/=mag;

  //cross u with normal to get v (vertical)
  vx =   uy*bz - by*uz;
  vy = - ux*bz + bx*uz;
  vz =   ux*by - bx*uy;
  mag = sqrt(vx*vx + vy*vy + vz*vz);
  if (mag < DBL_TOL) mag = 1.0;
  vx/=mag;
  vy/=mag;
  vz/=mag;

  x1 = x1 - halfwidth;
  y1 = y1 - halfheight;
  x2 = x2 - halfwidth;
  y2 = y2 - halfheight;

  new_x1 = atx + x1*ux - y1*vx;
  new_y1 = aty + x1*uy - y1*vy;
  new_z1 = atz + x1*uz - y1*vz;

  new_x2 = atx + x2*ux - y2*vx;
  new_y2 = aty + x2*uy - y2*vy;
  new_z2 = atz + x2*uz - y2*vz;

  ax = new_x2 - new_x1;
  ay = new_y2 - new_y1;
  az = new_z2 - new_z1;

  cx = new_x2 - eyex;
  cy = new_y2 - eyey;
  cz = new_z2 - eyez;

  //Plane normal
  nx =  cy*az - ay*cz;
  ny = -cx*az + ax*cz;
  nz =  cx*ay - ax*cy;

  mag = sqrt(nx*nx + ny*ny + nz*nz);
  if (mag < DBL_TOL) mag = 1.0;
  nx/=mag;
  ny/=mag;
  nz/=mag;
  //Offset to plane
  d = nx*new_x1 + ny*new_y1 + nz*new_z1;

  wh = width*height;

//  printf("nx: %f ny: %f nz: %f\n", nx, ny, nz);

  printf("Cutting data...  ");
  for (i=0; i<width; i++){
    for (j=0; j<height; j++){
      for (k=0; k<depth; k++){
        if (nx*(double)i + ny*(double)j + nz*(double)k - d > 0)
          data[i + j*width + k*wh] = 0;
      }
    }
  }
  printf("Finished.\n");

  needs_update = TRUE;
}

int Volume::read(char *filename){
  FILE *fp;
  size_t num_read;

  fp = fopen(filename, "rb");
  num_read = fread(data, 1, height * width *depth, fp );
  if (num_read != (size_t) (width*height*depth)){
    printf("Error reading file %s. num_read = %d\n", filename, num_read);
    return FALSE;
  }
  needs_update = TRUE;
  return TRUE;
}

void Volume::setIntersectColor(int im_loc,
                               double dx, double dy, double dz){ 
  double t1;
  double enter_t;
  int enter_px=0, enter_py=0, enter_pz=0;
  double px, py, pz;
  int exit_px=0, exit_py=0, exit_pz=0;
  int enter_face;

  /*
double b0f0nx =  0.0, b0f0ny =  0.0, b0f0nz = -1.0;
double b0f1nx =  0.0, b0f1ny =  1.0, b0f1nz =  0.0;
double b0f2nx =  1.0, b0f2ny =  0.0, b0f2nz =  0.0;
double b0f3nx =  0.0, b0f3ny = -1.0, b0f3nz =  0.0;
double b0f4nx = -1.0, b0f4ny =  0.0, b0f4nz =  0.0;
double b0f5nx =  0.0, b0f5ny =  0.0, b0f5nz =  1.0;

t = (D - dot(normal, eye) / dot(normal, direction);
D = distance from plane to origin.

NOTE:  Box sits with bottom left corner at origin in space
*/

  enter_face = 6;
  enter_t = 1000000;
    //Face 0 front
  t1 = eyez/-dz;
  if (t1 > 0.0){
    px = eyex + dx*t1; py = eyey + dy*t1; pz = 0.0;
    if (px < d_width && py < d_height && px >= 0.0 && py >= 0.0){
      enter_t = t1;
      enter_face = 0;
      enter_px = (int)px;
      enter_py = (int)py;
      enter_pz = (int)(pz*zmult);
    }
  }
    //Face 1 top
  t1 = (d_height - eyey)/dy;
  if (t1 > 0.0){// if (t1 < enter_t){
    px = eyex + dx*t1; py = d_height; pz = eyez + dz*t1;
    if (px < d_width && pz < d_depth && px >= 0.0 && pz >= 0.0){
      //We hit a side, not determine enter and exit points.
      if (t1 < enter_t){
        //This is the entering face
        exit_px = enter_px;
        exit_py = enter_py;
        exit_pz = enter_pz;

        enter_t = t1;
        enter_face = 1;
        enter_px = (int)px;
        enter_py = (int)py;
        enter_pz = (int)(pz*zmult);
      }
      else{
        //This is the exiting face.
        exit_px = (int)px;
        exit_py = (int)py;
        exit_pz = (int)(pz*zmult);
      }
    }
  }
    //Face 2 right side
  t1 = (d_width - eyex )/dx;
  if (t1 > 0.0){
    px = d_width; py = eyey + dy*t1; pz = eyez + dz*t1;
    if (pz < d_depth && py < d_height && pz >= 0.0 && py >= 0.0){
      if (t1 < enter_t){
        exit_px = enter_px;
        exit_py = enter_py;
        exit_pz = enter_pz;

        enter_t = t1;
        enter_face = 2;
        enter_px = (int)px;
        enter_py = (int)py;
        enter_pz = (int)(pz*zmult);
      }
      else{
        exit_px = (int)px;
        exit_py = (int)py;
        exit_pz = (int)(pz*zmult);
      }
    }
  }
    //Face 3 bottom
  t1 = eyey/-dy;
  if (t1 > 0.0){
    px = eyex + dx*t1; py = 0.0; pz = eyez + dz*t1;
    if (px < d_width && pz < d_depth && px >= 0.0 && pz >= 0.0){
      if (t1 < enter_t){
        exit_px = enter_px;
        exit_py = enter_py;
        exit_pz = enter_pz;

        enter_t = t1;
        enter_face = 3;
        enter_px = (int)px;
        enter_py = (int)py;
        enter_pz = (int)(pz*zmult);
      }
      else{
        exit_px = (int)px;
        exit_py = (int)py;
        exit_pz = (int)(pz*zmult);
      }
    }
  }
    //Face 4 left side
  t1 = eyex/-dx;
  if (t1 > 0.0){
    px = 0.0; py = eyey + dy*t1; pz = eyez + dz*t1;
    if (pz < d_depth && py < d_height && pz >= 0.0 && py >= 0.0){
      if (t1 < enter_t){
        exit_px = enter_px;
        exit_py = enter_py;
        exit_pz = enter_pz;

        enter_t = t1;
        enter_face = 4;
        enter_px = (int)px;
        enter_py = (int)py;
        enter_pz = (int)(pz*zmult);
      }
      else{
        exit_px = (int)px;
        exit_py = (int)py;
        exit_pz = (int)(pz*zmult);
      }
    }
  }
    //Face 5 back
  t1 = (d_depth - eyez)/dz;
  if (t1 > 0.0){
    px = eyex + dx*t1; py = eyey + dy*t1; pz = d_depth;
    
    if (px < d_width && py < d_height && px >= 0.0 && py >= 0.0){
      if (t1 < enter_t){
        exit_px = enter_px;
        exit_py = enter_py;
        exit_pz = enter_pz;

        enter_t = t1;
        enter_face = 5;
        enter_px = (int)px;
        enter_py = (int)py;
        enter_pz = (int)(pz*zmult);
      }
      else{
        exit_px = (int)px;
        exit_py = (int)py;
        exit_pz = (int)(pz*zmult);
      }
    }
  }
  

  switch (render_type){
    case FLAT:
      switch (enter_face){
        case 0:
          im_data[im_loc] = 255;
          im_data[im_loc+1] =
                 data[enter_px + enter_py*width + g_slice_id*width*height];
//                 data[exit_px + exit_py*width + g_slice_id*width*height];
          im_data[im_loc+2] =
                 data[enter_px + enter_py*width + g_slice_id*width*height];
//                 data[exit_px + exit_py*width + g_slice_id*width*height];
          break;
        case 1:
          im_data[im_loc] = 0;
          im_data[im_loc+1] = 255;
          im_data[im_loc+2] = 0; break; case 2:
          im_data[im_loc] = 0;
          im_data[im_loc+1] = 0;
          im_data[im_loc+2] = 255;
          break;
        case 3:
          im_data[im_loc] = 255;
          im_data[im_loc+1] = 255;
          im_data[im_loc+2] = 0;
          break;
        case 4:
          im_data[im_loc] = 0;
          im_data[im_loc+1] = 255;
          im_data[im_loc+2] = 255;
          break;
        case 5:
          im_data[im_loc] = 255;
          im_data[im_loc+1] = 0;
          im_data[im_loc+2] = 255;
          break;
        case 6:
          im_data[im_loc] = bgr;
          im_data[im_loc+1] = bgg;
          im_data[im_loc+2] = bgb;
          break;
      }
      break;

    case DEPTH:
      if (enter_face > 5){
        im_data[im_loc] = bgr;
        im_data[im_loc+1] = bgg;
        im_data[im_loc+2] = bgb;
        return;
      }
      setDepthColor(enter_px, enter_py, enter_pz,
                     exit_px, exit_py, exit_pz, im_loc, enter_t);
      break;
      
    case COLOR:
      if (enter_face > 5){
        im_data[im_loc] = bgr;
        im_data[im_loc+1] = bgg;
        im_data[im_loc+2] = bgb;
        return;
      }
      setDataColor(enter_px, enter_py, enter_pz,
               exit_px, exit_py, exit_pz, im_loc);
      break;

    case DENSITY:
      if (enter_face > 5){
        im_data[im_loc] = bgr;
        im_data[im_loc+1] = bgg;
        im_data[im_loc+2] = bgb;
        return;
      }
      setDensityColor(enter_px, enter_py, enter_pz,
               exit_px, exit_py, exit_pz, im_loc, enter_t);
      break;

    case NORMAL:
      if (enter_face > 5){
        im_data[im_loc] = bgr;
        im_data[im_loc+1] = bgg;
        im_data[im_loc+2] = bgb;
        return;
      }
      setNormalColor(enter_px, enter_py, enter_pz,
               exit_px, exit_py, exit_pz, im_loc);
      break;

    case NICE:
      if (enter_face > 5){
        im_data[im_loc] = bgr;
        im_data[im_loc+1] = bgg;
        im_data[im_loc+2] = bgb;
        return;
      }
      setNiceColor(enter_px, enter_py, enter_pz,
               exit_px, exit_py, exit_pz, im_loc);
      break;

    default:
      break;
  }

}

void Volume::setDataColor(int x0, int y0, int z0, int x1, int y1, int z1, int loc){
  int dx, dy, dz;
  int magx, magy, magz;
  int flag;
  unsigned char ddist;
  unsigned char val = 0;

  dx = x1 - x0;
  dy = y1 - y0;
  dz = z1 - z0;
  magx = abs(dx);
  magy = abs(dy);
  magz = abs(dz);

  if (magx > magy){
    if (magx > magz){
      //Mag x is the biggest
      if (magx == 0 ) flag = 3;
      else flag = 0;
    }
    else{
      //Mag z is the biggest
      if (magz == 0 ) flag = 3;
      else flag = 2;
    }
  }
  else if (magy > magz){
    //Mag y is the biggest
    if (magy == 0 ) flag = 3;
    else flag = 1;
  }
  else{
    //Mag z is the biggest
    if (magz == 0 ) flag = 3;
    else flag = 2;
  }


  ddist = 255;


  switch (flag){
    case 0://X biggest
      dy = (dy<<16)/magx;
      dz = (dz<<16)/magx;
      dx = (dx<<16)/magx;
      x0 <<= 16;
      y0 <<= 16;
      z0 <<= 16;
      x1 <<= 16;
      while (x0!=x1){
        //plot/check point
        val = data[(x0>>16) + (y0>>16)*width + (z0>>16)*width*height];
        if (val >= threshhold && val <= max_threshhold){
          im_data[loc] = val;
          im_data[loc+1] = val;
          im_data[loc+2] = val;
          return;
        }
        x0+=dx;
        y0+=dy;
        z0+=dz;
        if (ddist > 0) ddist --;
      }
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      im_data[loc+2] = bgb;
      break;
    case 1://Y biggest
      dy = (dy<<16)/magy;
      dz = (dz<<16)/magy;
      dx = (dx<<16)/magy;
      x0 <<= 16;
      y0 <<= 16;
      z0 <<= 16;
      y1 <<= 16;
      while (y0!=y1){
        //plot/check point
        val = data[(x0>>16) + (y0>>16)*width + (z0>>16)*width*height];
        if (val >= threshhold && val <= max_threshhold){
          im_data[loc] = val;
          im_data[loc+1] = val;
          im_data[loc+2] = val;
          return;
        }
        x0+=dx;
        y0+=dy;
        z0+=dz;
        if (ddist > 0) ddist --;
      }
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      im_data[loc+2] = bgb;
      break;
    case 2://Z biggest
      dy = (dy<<16)/magz;
      dz = (dz<<16)/magz;
      dx = (dx<<16)/magz;
      x0 <<= 16;
      y0 <<= 16;
      z0 <<= 16;
      z1 <<= 16;
      while (z0!=z1){
        //plot/check point
        val = data[(x0>>16) + (y0>>16)*width + (z0>>16)*width*height];
        if (val >= threshhold && val <= max_threshhold){
          im_data[loc] = val;
          im_data[loc+1] = val;
          im_data[loc+2] = val;
          return;
        }
        x0+=dx;
        y0+=dy;
        z0+=dz;
        if (ddist > 0) ddist --;
      }
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      im_data[loc+2] = bgb;
      break;
    case 3:
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      im_data[loc+2] = bgb;
      //Do nothing
      break;
  }

}

void Volume::setDensityColor(int x0, int y0, int z0, int x1, int y1, int z1, int loc, double idist){
  int dx, dy, dz;
  int magx, magy, magz;
  int flag;
  unsigned char ddist;
  unsigned int val = 0;
  unsigned int den_t;
  double temp;
  unsigned char tval = 0;
  unsigned char addin = 0;
  int check_tval = show_skin;

  den_t = ((unsigned int) threshhold*width/2);
 // Image img(256, 256, 6);

  dx = x1 - x0;
  dy = y1 - y0;
  dz = z1 - z0;
  magx = abs(dx);
  magy = abs(dy);
  magz = abs(dz);
//  printf("magx: %d magy: %d magz: %d\n", magx, magy, magz);

//  img.clear(128);
  if (magx > magy){
    if (magx > magz){
      //Mag x is the biggest
      if (magx == 0 ) flag = 3;
      else flag = 0;
    }
    else{
      //Mag z is the biggest
      if (magz == 0 ) flag = 3;
      else flag = 2;
    }
  }
  else if (magy > magz){
    //Mag y is the biggest
    if (magy == 0 ) flag = 3;
    else flag = 1;
  }
  else{
    //Mag z is the biggest
    if (magz == 0 ) flag = 3;
    else flag = 2;
  }


//  ddist = 255;


  switch (flag){
    case 0://X biggest
//      printf("magx: %d magy: %d magz: %d\n", magx, magy, magz);
      dy = (dy<<16)/magx;
      dz = (dz<<16)/magx;
      dx = (dx<<16)/magx;
      x0 <<= 16;
      y0 <<= 16;
      z0 <<= 16;
      x1 <<= 16;
      im_data[loc+2] = bgb;
      while (x0!=x1){
        //plot/check point
        tval = data[(x0>>16) + (y0>>16)*width + (z0>>16)*width*height];
        val += (unsigned int) tval;
        if (check_tval && tval > threshhold){
          temp = (idist-depth_dist)/depth_scale;
          if (temp > 1.0) addin = 0;
          else if (temp < 0.0) addin = 255;
          else{
            addin = (unsigned char)(255.0 - 255.0*temp);
          }
          im_data[loc+2] = addin;
          check_tval = FALSE;
        }

        if (val >= den_t){
          temp = (idist-depth_dist)/depth_scale;
          if (temp > 1.0) ddist = 0;
          else if (temp < 0.0) ddist = 255;
          else{
            ddist = (unsigned char)(255.0 - 255.0*temp);
          }

          im_data[loc] = ddist;
          im_data[loc+1] = ddist;
          return;
        }
        x0+=dx;
        y0+=dy;
        z0+=dz;
        idist++;
      }
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      break;


    case 1://Y biggest
      dy = (dy<<16)/magy;
      dz = (dz<<16)/magy;
      dx = (dx<<16)/magy;
      x0 <<= 16;
      y0 <<= 16;
      z0 <<= 16;
      y1 <<= 16;
      im_data[loc+2] = bgb;
      while (y0!=y1){
        //plot/check point
        tval = data[(x0>>16) + (y0>>16)*width + (z0>>16)*width*height];
        val += (unsigned int) tval;
        if (check_tval && tval > threshhold){
          temp = (idist-depth_dist)/depth_scale;
          if (temp > 1.0) addin = 0;
          else if (temp < 0.0) addin = 255;
          else{
            addin = (unsigned char)(255.0 - 255.0*temp);
          }
          im_data[loc+2] = addin;
          check_tval = FALSE;
        }

        if (val >= den_t){
          temp = (idist-depth_dist)/depth_scale;
          if (temp > 1.0) ddist = 0;
          else if (temp < 0.0) ddist = 255;
          else{
            ddist = (unsigned char)(255.0 - 255.0*temp);
          }

          im_data[loc] = ddist;
          im_data[loc+1] = ddist;
          return;
        }
        x0+=dx;
        y0+=dy;
        z0+=dz;
        idist++;
      }
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      break;


    case 2://Z biggest
      dy = (dy<<16)/magz;
      dz = (dz<<16)/magz;
      dx = (dx<<16)/magz;
      x0 <<= 16;
      y0 <<= 16;
      z0 <<= 16;
      z1 <<= 16;
      im_data[loc+2] = bgb;
      while (z0!=z1){
        //plot/check point
        tval = data[(x0>>16) + (y0>>16)*width + (z0>>16)*width*height];
        val += (unsigned int) tval;
        if (check_tval && tval > threshhold){
          temp = (idist-depth_dist)/depth_scale;
          if (temp > 1.0) addin = 0;
          else if (temp < 0.0) addin = 255;
          else{
            addin = (unsigned char)(255.0 - 255.0*temp);
          }
          im_data[loc+2] = addin;
          check_tval = FALSE;
        }

        if (val >= den_t){
          temp = (idist-depth_dist)/depth_scale;
          if (temp > 1.0) ddist = 0;
          else if (temp < 0.0) ddist = 255;
          else{
            ddist = (unsigned char)(255.0 - 255.0*temp);
          }

          im_data[loc] = ddist;
          im_data[loc+1] = ddist;
          return;
        }
        x0+=dx;
        y0+=dy;
        z0+=dz;
        idist++;
      }
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      break;
      

    case 3:
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      im_data[loc+2] = bgb;
      //Do nothing
      break;
  }

}


void Volume::setDepthColor(int x0, int y0, int z0, int x1, int y1, int z1, int loc, double idist){
  int dx, dy, dz;
  int magx, magy, magz;
  int flag;
  unsigned char ddist;
  unsigned char val = 0;
  double temp;

  dx = x1 - x0;
  dy = y1 - y0;
  dz = z1 - z0;
  magx = abs(dx);
  magy = abs(dy);
  magz = abs(dz);

//  img.clear(128);
  if (magx > magy){
    if (magx > magz){
      //Mag x is the biggest
      if (magx == 0 ) flag = 3;
      else flag = 0;
    }
    else{
      //Mag z is the biggest
      if (magz == 0 ) flag = 3;
      else flag = 2;
    }
  }
  else if (magy > magz){
    //Mag y is the biggest
    if (magy == 0 ) flag = 3;
    else flag = 1;
  }
  else{
    //Mag z is the biggest
    if (magz == 0 ) flag = 3;
    else flag = 2;
  }


  ddist = 255;


  switch (flag){
    case 0://X biggest
      dy = (dy<<16)/magx;
      dz = (dz<<16)/magx;
      dx = (dx<<16)/magx;
      x0 <<= 16;
      y0 <<= 16;
      z0 <<= 16;
      x1 <<= 16;
      while (x0!=x1){
        //plot/check point
        val = data[(x0>>16) + (y0>>16)*width + (z0>>16)*width*height];
        if (val >= threshhold && val <= max_threshhold){

          temp = (idist-depth_dist)/depth_scale;
          if (temp > 1.0) ddist = 0;
          else if (temp < 0.0) ddist = 255;
          else{
            ddist = (unsigned char)(255.0 - 255.0*temp);
          }

          im_data[loc] =   ddist;
          im_data[loc+1] = ddist;
          im_data[loc+2] = ddist;
          return;
        }
        x0+=dx;
        y0+=dy;
        z0+=dz;

        idist ++;
      }
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      im_data[loc+2] = bgb;
      break;
    case 1://Y biggest
      dy = (dy<<16)/magy;
      dz = (dz<<16)/magy;
      dx = (dx<<16)/magy;
      x0 <<= 16;
      y0 <<= 16;
      z0 <<= 16;
      y1 <<= 16;
      while (y0!=y1){
        //plot/check point
        val = data[(x0>>16) + (y0>>16)*width + (z0>>16)*width*height];
        if (val >= threshhold && val <= max_threshhold){

          temp = (idist-depth_dist)/depth_scale;
          if (temp > 1.0) ddist = 0;
          else if (temp < 0.0) ddist = 255;
          else{
            ddist = (unsigned char)(255.0 - 255.0*temp);
          }

          im_data[loc] = ddist;
          im_data[loc+1] = ddist;
          im_data[loc+2] = ddist;
          return;
        }
        x0+=dx;
        y0+=dy;
        z0+=dz;
        //if (ddist > 0) ddist --;
        idist ++;
      }
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      im_data[loc+2] = bgb;
      break;
    case 2://Z biggest
      dy = (dy<<16)/magz;
      dz = (dz<<16)/magz;
      dx = (dx<<16)/magz;
      x0 <<= 16;
      y0 <<= 16;
      z0 <<= 16;
      z1 <<= 16;
      while (z0!=z1){
        //plot/check point
        val = data[(x0>>16) + (y0>>16)*width + (z0>>16)*width*height];

        if (val >= threshhold && val <= max_threshhold){

          temp = (idist-depth_dist)/depth_scale;
          if (temp > 1.0) ddist = 0;
          else if (temp < 0.0) ddist = 255;
          else{
            ddist = (unsigned char)(255.0 - 255.0*temp);
          }

          im_data[loc] = ddist;
          im_data[loc+1] = ddist;
          im_data[loc+2] = ddist;
          return;
        }
        x0+=dx;
        y0+=dy;
        z0+=dz;
        //if (ddist > 0) ddist --;
        idist ++;
      }
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      im_data[loc+2] = bgb;
      break;
    case 3:
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      im_data[loc+2] = bgb;
      //Do nothing
      break;
  }

}

void Volume::setNiceColor(int x0, int y0, int z0, int x1, int y1, int z1, int loc){
  int dx, dy, dz;
  int magx, magy, magz;
  int flag;
  double nx, ny, nz, mag;
  unsigned char ddist;
  unsigned char val = 0;
  unsigned char lval = 0;
  register int base_loc;
  register int y_off = width;
  register int z_off = width*height;

  //The base pixel is v14
  //Front face (posative x) = v1 -> v9
  //Back face (negative x)  = v19 -> v27
  double v1, v2, v3, v4, v5, v6, v7, v8, v9, v10, v11, v12, v13, v15, v16,
         v17, v18, v19, v20, v21, v22, v23, v24, v25, v26, v27;

  dx = x1 - x0;
  dy = y1 - y0;
  dz = z1 - z0;
  magx = abs(dx);
  magy = abs(dy);
  magz = abs(dz);

  if (magx > magy){
    if (magx > magz){
      //Mag x is the biggest
      if (magx == 0 ) flag = 3;
      else flag = 0;
    }
    else{
      //Mag z is the biggest
      if (magz == 0 ) flag = 3;
      else flag = 2;
    }
  }
  else if (magy > magz){
    //Mag y is the biggest
    if (magy == 0 ) flag = 3;
    else flag = 1;
  }
  else{
    //Mag z is the biggest
    if (magz == 0 ) flag = 3;
    else flag = 2;
  }


  ddist = 255;

  switch (flag){
    case 0://X biggest
//      printf("magx: %d magy: %d magz: %d\n", magx, magy, magz);
      dy = (dy<<16)/magx;
      dz = (dz<<16)/magx;
      dx = (dx<<16)/magx;
      x0 <<= 16;
      y0 <<= 16;
      z0 <<= 16;
      x1 <<= 16;
      while (x0!=x1){
        //plot/check point
        base_loc = (x0>>16) + (y0>>16)*width + (z0>>16)*width*height;
        val = data[base_loc];
        if (val >= threshhold && val <= max_threshhold){

          //Front (in x)
          v1 = INV_ROOT3*(double)data[base_loc + 1 + y_off - z_off];
          v2 = INV_ROOT2*(double)data[base_loc + 1 + y_off];
          v3 = INV_ROOT3*(double)data[base_loc + 1 + y_off + z_off];

          v4 = INV_ROOT2*(double)data[base_loc + 1 - z_off];
          v5 =       (double)data[base_loc + 1 ];
          v6 = INV_ROOT2*(double)data[base_loc + 1 + z_off];

          v7 = INV_ROOT3*(double)data[base_loc + 1 - y_off - z_off];
          v8 = INV_ROOT2*(double)data[base_loc + 1 - y_off];
          v9 = INV_ROOT3*(double)data[base_loc + 1 - y_off + z_off];


          //Middle
          v10 = INV_ROOT2*(double)data[base_loc + y_off - z_off];
          v11 =       (double)data[base_loc + y_off];
          v12 = INV_ROOT2*(double)data[base_loc + y_off + z_off];

          v13 =       (double)data[base_loc - z_off];
//          v14 =       (double)data[base_loc ]; //original pixel location
          v15 =       (double)data[base_loc + z_off];

          v16 = INV_ROOT2*(double)data[base_loc - y_off - z_off];
          v17 =       (double)data[base_loc - y_off];
          v18 = INV_ROOT2*(double)data[base_loc - y_off + z_off];

          //Back
          v19 = INV_ROOT3*(double)data[base_loc - 1 + y_off - z_off];
          v20 = INV_ROOT2*(double)data[base_loc - 1 + y_off];
          v21 = INV_ROOT3*(double)data[base_loc - 1 + y_off + z_off];

          v22 = INV_ROOT2*(double)data[base_loc - 1 - z_off];
          v23 =       (double)data[base_loc - 1 ];
          v24 = INV_ROOT2*(double)data[base_loc - 1 + z_off];

          v25 = INV_ROOT3*(double)data[base_loc - 1 - y_off - z_off];
          v26 = INV_ROOT2*(double)data[base_loc - 1 - y_off];
          v27 = INV_ROOT3*(double)data[base_loc - 1 - y_off + z_off];

          nx = v1 + v2 + v3 + v4 + v5 + v6 + v7 + v8 + v9
              - v19 - v20 - v21 - v22 - v23 - v24 - v25 - v26 - v27;

          ny = v1 + v2 + v3 + v10 + v11 + v12 + v19 + v20 + v21
              - v7 - v8 - v9 - v16 - v17 - v18 - v25 - v26 - v27;

          nz = v21 + v12 + v3 + v24 + v15 + v6 + v27 + v18 + v9
              - v19 - v10 - v1 - v22 - v13 - v4 - v25 - v16 - v7;

          mag = sqrt(nx*nx + ny*ny + nz*nz);
          if (mag<DBL_TOL) mag = 1.0;
          nx/=mag;
          ny/=mag;
          nz/=mag;
          //lval =  (unsigned char)(fabs(nz)*128 + fabs(nx*axis_x + ny*axis_y + nz*axis_z)*128);//(unsigned char)(fabs(nz)*255);
          lval =  (unsigned char)(fabs(nx*axis_x + ny*axis_y + nz*axis_z)*255);
          im_data[loc] = lval;
          im_data[loc+1] =  lval;
          im_data[loc+2] =  lval;

          return;
        }
        x0+=dx;
        y0+=dy;
        z0+=dz;
        if (ddist > 0) ddist --;
      }
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      im_data[loc+2] = bgb;
      break;



    case 1://Y biggest
      dy = (dy<<16)/magy;
      dz = (dz<<16)/magy;
      dx = (dx<<16)/magy;
      x0 <<= 16;
      y0 <<= 16;
      z0 <<= 16;
      y1 <<= 16;
      while (y0!=y1){
        //plot/check point
        base_loc = (x0>>16) + (y0>>16)*width + (z0>>16)*width*height;
        val = data[base_loc];
 //       data[x0>>16 + (y0>>16)*width + (z0>>16)*width*height] = 255;
        if (val >= threshhold && val <= max_threshhold){

          //Front (in x)
          v1 = INV_ROOT3*(double)data[base_loc + 1 + y_off - z_off];
          v2 = INV_ROOT2*(double)data[base_loc + 1 + y_off];
          v3 = INV_ROOT3*(double)data[base_loc + 1 + y_off + z_off];

          v4 = INV_ROOT2*(double)data[base_loc + 1 - z_off];
          v5 =       (double)data[base_loc + 1 ];
          v6 = INV_ROOT2*(double)data[base_loc + 1 + z_off];

          v7 = INV_ROOT3*(double)data[base_loc + 1 - y_off - z_off];
          v8 = INV_ROOT2*(double)data[base_loc + 1 - y_off];
          v9 = INV_ROOT3*(double)data[base_loc + 1 - y_off + z_off];


          //Middle
          v10 = INV_ROOT2*(double)data[base_loc + y_off - z_off];
          v11 =       (double)data[base_loc + y_off];
          v12 = INV_ROOT2*(double)data[base_loc + y_off + z_off];

          v13 =       (double)data[base_loc - z_off];
//          v14 =       (double)data[base_loc ]; //original pixel location
          v15 =       (double)data[base_loc + z_off];

          v16 = INV_ROOT2*(double)data[base_loc - y_off - z_off];
          v17 =       (double)data[base_loc - y_off];
          v18 = INV_ROOT2*(double)data[base_loc - y_off + z_off];

          //Back
          v19 = INV_ROOT3*(double)data[base_loc - 1 + y_off - z_off];
          v20 = INV_ROOT2*(double)data[base_loc - 1 + y_off];
          v21 = INV_ROOT3*(double)data[base_loc - 1 + y_off + z_off];

          v22 = INV_ROOT2*(double)data[base_loc - 1 - z_off];
          v23 =       (double)data[base_loc - 1 ];
          v24 = INV_ROOT2*(double)data[base_loc - 1 + z_off];

          v25 = INV_ROOT3*(double)data[base_loc - 1 - y_off - z_off];
          v26 = INV_ROOT2*(double)data[base_loc - 1 - y_off];
          v27 = INV_ROOT3*(double)data[base_loc - 1 - y_off + z_off];

          nx = v1 + v2 + v3 + v4 + v5 + v6 + v7 + v8 + v9
              - v19 - v20 - v21 - v22 - v23 - v24 - v25 - v26 - v27;

          ny = v1 + v2 + v3 + v10 + v11 + v12 + v19 + v20 + v21
              - v7 - v8 - v9 - v16 - v17 - v18 - v25 - v26 - v27;

          nz = v21 + v12 + v3 + v24 + v15 + v6 + v27 + v18 + v9
              - v19 - v10 - v1 - v22 - v13 - v4 - v25 - v16 - v7;

          mag = sqrt(nx*nx + ny*ny + nz*nz);
          if (mag<DBL_TOL) mag = 1.0;
          nx/=mag;
          ny/=mag;
          nz/=mag;
          lval =  (unsigned char)(fabs(nx*axis_x + ny*axis_y + nz*axis_z)*255);
          im_data[loc] = lval;
          im_data[loc+1] = lval;
          im_data[loc+2] = lval;

          return;
        }
        x0+=dx;
        y0+=dy;
        z0+=dz;
        if (ddist > 0) ddist --;
      }
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      im_data[loc+2] = bgb;
      break;
    case 2://Z biggest
      dy = (dy<<16)/magz;
      dz = (dz<<16)/magz;
      dx = (dx<<16)/magz;
      x0 <<= 16;
      y0 <<= 16;
      z0 <<= 16;
      z1 <<= 16;
      while (z0!=z1){
        //plot/check point
        base_loc = (x0>>16) + (y0>>16)*width + (z0>>16)*width*height;
        val = data[base_loc];

        if (val >= threshhold && val <= max_threshhold){
          
          //Front (in x)
          v1 = INV_ROOT3*(double)data[base_loc + 1 + y_off - z_off];
          v2 = INV_ROOT2*(double)data[base_loc + 1 + y_off];
          v3 = INV_ROOT3*(double)data[base_loc + 1 + y_off + z_off];

          v4 = INV_ROOT2*(double)data[base_loc + 1 - z_off];
          v5 =       (double)data[base_loc + 1 ];
          v6 = INV_ROOT2*(double)data[base_loc + 1 + z_off];

          v7 = INV_ROOT3*(double)data[base_loc + 1 - y_off - z_off];
          v8 = INV_ROOT2*(double)data[base_loc + 1 - y_off];
          v9 = INV_ROOT3*(double)data[base_loc + 1 - y_off + z_off];


          //Middle
          v10 = INV_ROOT2*(double)data[base_loc + y_off - z_off];
          v11 =       (double)data[base_loc + y_off];
          v12 = INV_ROOT2*(double)data[base_loc + y_off + z_off];

          v13 =       (double)data[base_loc - z_off];
//          v14 =       (double)data[base_loc ]; //original pixel location
          v15 =       (double)data[base_loc + z_off];

          v16 = INV_ROOT2*(double)data[base_loc - y_off - z_off];
          v17 =       (double)data[base_loc - y_off];
          v18 = INV_ROOT2*(double)data[base_loc - y_off + z_off];

          //Back
          v19 = INV_ROOT3*(double)data[base_loc - 1 + y_off - z_off];
          v20 = INV_ROOT2*(double)data[base_loc - 1 + y_off];
          v21 = INV_ROOT3*(double)data[base_loc - 1 + y_off + z_off];

          v22 = INV_ROOT2*(double)data[base_loc - 1 - z_off];
          v23 =       (double)data[base_loc - 1 ];
          v24 = INV_ROOT2*(double)data[base_loc - 1 + z_off];

          v25 = INV_ROOT3*(double)data[base_loc - 1 - y_off - z_off];
          v26 = INV_ROOT2*(double)data[base_loc - 1 - y_off];
          v27 = INV_ROOT3*(double)data[base_loc - 1 - y_off + z_off];

          nx = v1 + v2 + v3 + v4 + v5 + v6 + v7 + v8 + v9
              - v19 - v20 - v21 - v22 - v23 - v24 - v25 - v26 - v27;

          ny = v1 + v2 + v3 + v10 + v11 + v12 + v19 + v20 + v21
              - v7 - v8 - v9 - v16 - v17 - v18 - v25 - v26 - v27;

          nz = v21 + v12 + v3 + v24 + v15 + v6 + v27 + v18 + v9
              - v19 - v10 - v1 - v22 - v13 - v4 - v25 - v16 - v7;

          mag = sqrt(nx*nx + ny*ny + nz*nz);
          if (mag<DBL_TOL) mag = 1.0;
          nx/=mag;
          ny/=mag;
          nz/=mag;
          lval =  (unsigned char)(fabs(nx*axis_x + ny*axis_y + nz*axis_z)*255);
          im_data[loc] = lval;
          im_data[loc+1] = lval;
          im_data[loc+2] = lval;

          return;
        }
        x0+=dx;
        y0+=dy;
        z0+=dz;
        if (ddist > 0) ddist --;
      }
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      im_data[loc+2] = bgb;
      break;
    case 3:
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      im_data[loc+2] = bgb;
      //Do nothing
      break;
  }

}

void Volume::setNormalColor(int x0, int y0, int z0, int x1, int y1, int z1, int loc){
  int dx, dy, dz;
  int magx, magy, magz;
  int flag;
  double nx, ny, nz, mag;
  unsigned char ddist;
  unsigned char val = 0;
  register int base_loc;
  register int y_off = width;
  register int z_off = width*height;

  dx = x1 - x0;
  dy = y1 - y0;
  dz = z1 - z0;
  magx = abs(dx);
  magy = abs(dy);
  magz = abs(dz);

  if (magx > magy){
    if (magx > magz){
      //Mag x is the biggest
      if (magx == 0 ) flag = 3;
      else flag = 0;
    }
    else{
      //Mag z is the biggest
      if (magz == 0 ) flag = 3;
      else flag = 2;
    }
  }
  else if (magy > magz){
    //Mag y is the biggest
    if (magy == 0 ) flag = 3;
    else flag = 1;
  }
  else{
    //Mag z is the biggest
    if (magz == 0 ) flag = 3;
    else flag = 2;
  }


  ddist = 255;


  switch (flag){
    case 0://X biggest
//      printf("magx: %d magy: %d magz: %d\n", magx, magy, magz);
      dy = (dy<<16)/magx;
      dz = (dz<<16)/magx;
      dx = (dx<<16)/magx;
      x0 <<= 16;
      y0 <<= 16;
      z0 <<= 16;
      x1 <<= 16;
      while (x0!=x1){
        //plot/check point
        base_loc = (x0>>16) + (y0>>16)*width + (z0>>16)*width*height;
        val = data[base_loc];
        if (val >= threshhold && val <= max_threshhold){


          //Calculate the normal to the surface here.
          nx = (double)data[base_loc + 1]-
               (double)data[base_loc - 1];
          ny = (double)data[base_loc + y_off]-
               (double)data[base_loc - y_off];
          nz = (double)data[base_loc + z_off]-
               (double)data[base_loc - z_off];

          mag = sqrt(nx*nx + ny*ny + nz*nz);
          nx/=mag;
          ny/=mag;
          nz/=mag;
          im_data[loc] = (unsigned char)(fabs(nz)*255);//ddist;
          im_data[loc+1] = (unsigned char)(fabs(nz)*255);//ddist;
          im_data[loc+2] = (unsigned char)(fabs(nz)*255);//ddist;

          return;
        }
        x0+=dx;
        y0+=dy;
        z0+=dz;
        if (ddist > 0) ddist --;
      }
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      im_data[loc+2] = bgb;
      break;
    case 1://Y biggest
      dy = (dy<<16)/magy;
      dz = (dz<<16)/magy;
      dx = (dx<<16)/magy;
      x0 <<= 16;
      y0 <<= 16;
      z0 <<= 16;
      y1 <<= 16;
      while (y0!=y1){
        //plot/check point
        base_loc = (x0>>16) + (y0>>16)*width + (z0>>16)*width*height;
        val = data[base_loc];
        if (val >= threshhold && val <= max_threshhold){


          //Calculate the normal to the surface here.
          nx = (double)data[base_loc + 1]-
               (double)data[base_loc - 1];
          ny = (double)data[base_loc + y_off]-
               (double)data[base_loc - y_off];
          nz = (double)data[base_loc + z_off]-
               (double)data[base_loc - z_off];

          mag = sqrt(nx*nx + ny*ny + nz*nz);
          nx/=mag;
          ny/=mag;
          nz/=mag;
          im_data[loc] = (unsigned char)(fabs(nz)*255);//ddist;
          im_data[loc+1] = (unsigned char)(fabs(nz)*255);//ddist;
          im_data[loc+2] = (unsigned char)(fabs(nz)*255);//ddist;

          return;
        }
        x0+=dx;
        y0+=dy;
        z0+=dz;
        if (ddist > 0) ddist --;
      }
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      im_data[loc+2] = bgb;
      break;
    case 2://Z biggest
      dy = (dy<<16)/magz;
      dz = (dz<<16)/magz;
      dx = (dx<<16)/magz;
      x0 <<= 16;
      y0 <<= 16;
      z0 <<= 16;
      z1 <<= 16;
      while (z0!=z1){
        //plot/check point
        base_loc = (x0>>16) + (y0>>16)*width + (z0>>16)*width*height;
        val = data[base_loc];
        if (val >= threshhold && val <= max_threshhold){


          //Calculate the normal to the surface here.
          nx = (double)data[base_loc + 1]-
               (double)data[base_loc - 1];
          ny = (double)data[base_loc + y_off]-
               (double)data[base_loc - y_off];
          nz = (double)data[base_loc + z_off]-
               (double)data[base_loc - z_off];

          mag = sqrt(nx*nx + ny*ny + nz*nz);
          nx/=mag;
          ny/=mag;
          nz/=mag;
          im_data[loc] = (unsigned char)(fabs(nz)*255);//ddist;
          im_data[loc+1] = (unsigned char)(fabs(nz)*255);//ddist;
          im_data[loc+2] = (unsigned char)(fabs(nz)*255);//ddist;

          return;
        }
        x0+=dx;
        y0+=dy;
        z0+=dz;
        if (ddist > 0) ddist --;
      }
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      im_data[loc+2] = bgb;
      break;
    case 3:
      im_data[loc] = bgr;
      im_data[loc+1] = bgg;
      im_data[loc+2] = bgb;
      //Do nothing
      break;
  }

}

void Volume::rotate(double theta, double gamma){
  int i;
  int dim;
  double ct, st, cg, sg;
  double rx, ry, rz;

  ct = cos(theta);
  st = sin(theta);
  cg = cos(gamma);
  sg = sin(gamma);
  
  //Rotate the from (or eye) point
  rx = eyex - cx;
  ry = eyey - cy;
  rz = eyez - cz;
  eyex = ct*rx + st*rz + cx;
  eyey = st*sg*rx + cg*ry - ct*sg*rz + cy;
  eyez = -cg*st*rx + sg*ry + cg*ct*rz + cz;

  //Rotate the at point
  rx = atx - cx;
  ry = aty - cy;
  rz = atz - cz;
  atx = ct*rx + st*rz + cx;
  aty = st*sg*rx + cg*ry - ct*sg*rz + cy;
  atz = -cg*st*rx + sg*ry + cg*ct*rz + cz;

  //Rotate the up vector
  rx = upx;
  ry = upy;
  rz = upz;
  upx = ct*rx + st*rz;
  upy = st*sg*rx + cg*ry - ct*sg*rz;
  upz = -cg*st*rx + sg*ry + cg*ct*rz;

  dim = im_width*im_height;

  for (i=0; i<dim; i++){
    rx = dirx[i];
    ry = diry[i];
    rz = dirz[i];

    /*
    dirx[i] = ct*cg*rx - sg*ry + st*cg*rz;
    diry[i] = ct*sg*rx + cg*ry + st*sg*rz;
    dirz[i] = -st*rx           + ct*rz;
    */

    dirx[i] = ct*rx + st*rz;
    diry[i] = st*sg*rx + cg*ry - ct*sg*rz;
    dirz[i] = -cg*st*rx + sg*ry + cg*ct*rz;
    /*
    dirx[i] = ct*cg*rx - sg*ry + st*cg*rz;
    diry[i] = ct*sg*rx + cg*ry + st*sg*rz;
    dirz[i] = -st*rx           + ct*rz;
    */
  }
  needs_update = TRUE;

}

/*
 * Updates the disp image based on the volume data and the current viewing
 * options.
 */
void Volume::updateDisp(){

  int i, j, loc;
  int loc0, loc1, loc2, loc3, loc4, loc5, loc6, loc7;
  int xoff, yoff, xstep;
  unsigned char r, r_25, r_50;
  unsigned char g, g_25, g_50;
  unsigned char b, b_25, b_50;
  double mag;

  if (!needs_update) return;

  axis_x = atx-eyex;
  axis_y = aty-eyey;
  axis_z = atz-eyez;
  mag = sqrt(axis_x*axis_x + axis_y*axis_y + axis_z*axis_z);
  if (mag < DBL_TOL) mag = 1.0;
  axis_x /= mag; axis_y /= mag; axis_z /= mag;


  //Do ray casting for the volume data
  if (quality_rend){
    for (i=0, j=0; i<im_dim; i++, j+=BPP){
      setIntersectColor(j, dirx[i], diry[i], dirz[i]);
    }

  }
  else{
    for (i=1, j=3; i<im_dim; i+=2, j+=(BPP<<1)){
      setIntersectColor(j, dirx[i], diry[i], dirz[i]);
    }

    yoff = im_width*BPP;
    xoff = BPP;
    xstep = BPP<<1;

    //Now fill in the gaps
    loc =  xoff + yoff;//the cursor position (1,1)
    loc0 = xstep + yoff;//right 2, 1
    loc1 = xstep;//top right 2, 0 
    loc2 = xoff;//top 1, 0 
    loc3 = 0;//top left 0, 0 
    loc4 = yoff;//left 0, 1 
    loc5 = yoff+yoff;//bot left 0, 2 
    loc6 = xoff + yoff+yoff;//bot 1, 2 
    loc7 = xstep + yoff+yoff;//bot right 2, 2 

    while (loc<im_width*im_height*BPP-xoff){
      for (i=0; i<im_width-1; i+=2){
        r = im_data[loc++];
        g = im_data[loc++];
        b = im_data[loc++];
        r_25 = r>>2;
        r_50 = r>>1;
        g_25 = g>>2;
        g_50 = g>>1;
        b_25 = b>>2;
        b_50 = b>>1;

        //Add parts of this to neighboring pixels
        im_data[loc0++] = r_50;
        im_data[loc0++] = g_50;
        im_data[loc0++] = b_50;
        im_data[loc1++] += r_25;
        im_data[loc1++] += g_25;
        im_data[loc1++] += b_25;
        im_data[loc2++] += r_50;
        im_data[loc2++] += g_50;
        im_data[loc2++] += b_50;
        im_data[loc3++] += r_25;
        im_data[loc3++] += g_25;
        im_data[loc3++] += b_25;
        im_data[loc4++] += r_50;
        im_data[loc4++] += g_50;
        im_data[loc4++] += b_50;
        im_data[loc5++] += r_25;
        im_data[loc5++] += g_25;
        im_data[loc5++] += b_25;
        im_data[loc6++] = r_50;
        im_data[loc6++] = g_50;
        im_data[loc6++] = b_50;
        im_data[loc7++] = r_25;
        im_data[loc7++] = g_25;
        im_data[loc7++] = b_25;

        //skip the next pixel
        loc  += xoff;
        loc0 += xoff;
        loc1 += xoff;
        loc2 += xoff;
        loc3 += xoff;
        loc4 += xoff;
        loc5 += xoff;
        loc6 += xoff;
        loc7 += xoff;
      }
      //skip a row
      loc  += yoff;
      loc0 += yoff;
      loc1 += yoff;
      loc2 += yoff;
      loc3 += yoff;
      loc4 += yoff;
      loc5 += yoff;
      loc6 += yoff;
      loc7 += yoff;
    }
  }
  needs_update = FALSE;

}

//Saves out an image of what is on the screen.
void Volume::writeImage(){
  int i;
  FILE   *fpt;
  unsigned char *save_line;
  char filename[255];

  sprintf(filename, "pic%03d.ppm", curr_pic);
  printf("saving %s... ", filename);

  fpt = fopen( filename, "wb" );
  if ( !fpt ){
    printf("Error: Could not create file: %s\n", filename);
    return;
  }

  fprintf( fpt, "P6\n" );
  //if ( HeaderInfo!=NULL ) fprintf( F, "# %s\n", HeaderInfo );
  //else fprintf( F, "# Grayscale Image\n" );
  fprintf( fpt, "%d %d\n255\n", im_width, im_height );
  for ( i=0; i<im_height; i++ ) {
    save_line = &im_data[i*im_width*BPP];
    if ( fwrite( (void*)save_line, 1, (size_t)(im_width*BPP), fpt ) !=
         (unsigned int)(im_width*BPP)  ){
      puts("fail to write");
      return;
    }
  }
  printf("done\n");
  fclose( fpt );


  curr_pic++;
}

void Volume::zoom(double theta){
  double rx, ry, rz, mag;

  //First the eye point
  rx = eyex - cx;
  ry = eyey - cy;
  rz = eyez - cz;

  mag = sqrt(rx*rx + ry*ry + rz*rz);
  if (mag + theta < DBL_TOL) return;
  rx /= mag; ry /= mag; rz /= mag;

  mag += theta;

  depth_dist = mag - depth_scale;
  if (depth_dist < 0.0) depth_dist = 0.0;
  eyex = rx*mag + cx;
  eyey = ry*mag + cy;
  eyez = rz*mag + cz;

  //Now the at point. Note: neg_dist is negative.
  atx = rx*neg_dist + eyex;
  aty = ry*neg_dist + eyey;
  atz = rz*neg_dist + eyez;

  needs_update = TRUE;
}





//OPENGL FUNCTIONS
//OPENGL FUNCTIONS
//OPENGL FUNCTIONS
//OPENGL FUNCTIONS
//OPENGL FUNCTIONS

/******************************************************************************
 * Purpose:  used by OpenGL to display the image.
 *****************************************************************************/
void display (void) {
  double g_time_count;

  volume->updateDisp();
  glRasterPos2i( 0, 0 );
  glDrawPixels(volume->im_width, volume->im_height, GL_RGB, GL_UNSIGNED_BYTE, volume->im_data);
  if (CUT == mode){
    //Draw a line
    glColor3f(0.0, 1.0, 0.0);
    glBegin(GL_LINES);
      glVertex2f(p1x, volume->im_height - p1y);
      glVertex2f(p2x, volume->im_height - p2y);
    glEnd();
    glFlush();
    glutSwapBuffers();
  }
  else if (TIME == mode){
    glFlush();//finish rendering it all and dump it out
    glutSwapBuffers();//since we're using the GLUT_DOUBLE
    if (g_frames_count < 10.0){
      g_frames_count ++;
//      printf("g_frames_count: %f\n", g_frames_count);
      volume->forceUpdate();
      display();
    }
    else{
//      printf("start: %d, clock(): %d diff: %d\n", (int)g_clock_start, (int)clock(), (int)(clock()-g_clock_start));
      g_time_count = (double)(clock() - g_clock_start); 
      printf("fps: %f\n", CLOCKS_PER_SEC*g_frames_count/g_time_count);
      g_frames_count = 0.0;
      mode = NONE;
    }
  }
  else{
    glFlush();//finish rendering it all and dump it out
    glutSwapBuffers();//since we're using the GLUT_DOUBLE
  }

  /*
  wait_time = clock() + g_delay;
  while(wait_time > clock())
    ;
  glutPostRedisplay();
  */
}

void reshape(int w, int h){
  int zoomx, zoomy;

  if (w < h){
    g_scale = (double)w/(double)volume->im_width;
  }
  else{
    g_scale = (double)h/(double)volume->im_height;
  }

  zoomx = (int)(g_scale*(double)volume->im_width);
  zoomy = (int)(g_scale*(double)volume->im_height);
  glViewport(0, 0, (GLsizei) zoomx,  (GLsizei) zoomy);
  /*
  glMatrixMode ( GL_PROJECTION );
  glLoadIdentity();
  gluOrtho2D( 0.0, xres, 0.0, yres);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
  */

  glPixelZoom(g_scale, g_scale);

}

/******************************************************************************
 * Purpose:  initialize the opengl stuff
 *****************************************************************************/
void init () {
  glClearColor (0.0, 0.0, 0.0, 0.0);
  glMatrixMode (GL_PROJECTION);
  glLoadIdentity();
  gluOrtho2D(0.0, volume->im_width, 0.0, volume->im_height);
  glColor3f ( 1.0, 1.0, 1.0);
  glClear (GL_COLOR_BUFFER_BIT );//need to clear the buffer to do the z-buffer
  glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
//  glPixelZoom((float)XDIM/(float)IM_XDIM, (float)YDIM/(float)IM_YDIM);
  g_delay = (int)( (double)CLOCKS_PER_SEC/30.0);
}

/******************************************************************************
 * Purpose:  called by opengl to handle mouse events
 *****************************************************************************/
void mouseButton(int button, int state, int x, int y) {
  /*
  if (x<0) x = 0; if (x>volume->im_width-1) x = volume->im_width-1;
  if (y<0) y=0; if (y>volume->im_height-1) y = volume->im_height-1;
  */
//  y = volume->height - y;//switch the coordiates
  //if (CUT == mode) return;

  switch (button) {
    case GLUT_LEFT_BUTTON:
      if (state == GLUT_DOWN) {
        if (CUT == mode){
          p1x = (double)x/g_scale;
          p1y = (double)y/g_scale;
        }
        else{
          mode = ROTATE;
          g_x_prev = x;
          g_y_prev = y;
        }
      }
      if (state == GLUT_UP) {
        if (CUT == mode){
          p2x = (double)x/g_scale;
          p2y = (double)y/g_scale;
        }
        else{
          mode = NONE;
        }
      }
      break;

    case GLUT_RIGHT_BUTTON:
      if (state == GLUT_DOWN) {
        mode = ZOOM;
        g_x_prev = x;
        g_y_prev = y;
      }
      if (state == GLUT_UP) {
        mode = NONE;
      }
      break;
  }
}

/******************************************************************************
 * Purpose:  called by opengl to handle mouse events
 *****************************************************************************/
void mouseMotion(int x, int y) {
//  y = volume->height - y;//switch the coordiates
//  if (x<0) x = 0; if (x>volume->im_width-1) x = volume->im_width-1;
//  if (y<0) y=0; if (y>volume->im_height-1) y = volume->im_height-1;
  double theta, gamma;

  switch (mode){
    case NONE:
      break;
    case ROTATE:
      theta = 0.01*(double)(x - g_x_prev);
      gamma = 0.01*(double)(y - g_y_prev);
      volume->rotate(theta, gamma);
      g_x_prev = x;
      g_y_prev = y;
      glutPostRedisplay();
      break;

    case ZOOM:
      theta = (double)(y - g_y_prev);
      volume->zoom(theta);
      g_x_prev = x;
      g_y_prev = y;
      glutPostRedisplay();
      break;

    case CUT:
      p2x = (double)x/g_scale;
      p2y = (double)y/g_scale;
      glutPostRedisplay();
      break;

    default:
      break;
  }

}

/******************************************************************************
 * Purpose:  called by opengl to handle keyboard events
 *****************************************************************************/
void keyboard(unsigned char key, int x, int y) {
  if (x<0) x = 0; if (x>volume->im_width-1) x = volume->im_width-1;
  if (y<0) y=0; if (y>volume->im_height-1) y = volume->im_height-1;

  switch ( key ) {

    case 'a':
      g_slice_id++;
      if (g_slice_id >= volume->depth) g_slice_id = 0;
      glutPostRedisplay();
      break;

    case 'b':
      volume->toggleQualityRend();
      glutPostRedisplay();
      break;

    case 'c':
      if (CUT == mode){
        volume->cut(p1x, p1y, p2x, p2y);
      }
      else{
        volume->setRenderType(COLOR);
      }
      glutPostRedisplay();
      break;

    case 'd':
      volume->setRenderType(DEPTH);
      glutPostRedisplay();
      break;

    case 'e':
      volume->setRenderType(DENSITY);
      glutPostRedisplay();
      break;

    case 'f':
      volume->setRenderType(FLAT);
      glutPostRedisplay();
      break;

    case 'g':
      volume->setRenderType(NICE);
      glutPostRedisplay();
      break;

    case 'h':
      printf(" KEYS:\n");
      printf(" 'a' - ADVANCE the slice in Flat shade.\n");
      printf(" 'b' - Toggle BETTER rendering quality.\n");
      printf(" 'c' - Use the COLOR from the data set image.\n");
      printf(" 'd' - DEPTH shade.\n");
      printf(" 'e' - DENSITY shade.\n");
      printf(" 'f' - FLAT shade.\n");
      printf(" 'g' - GOOD nice normal shade.\n");
      printf(" 'h' - HELP menu.\n");
      printf(" 'i' - DILATE the data.\n");
      printf(" 'j' - THRESHOLD the data.\n");
      printf(" 'k' - Calculate FRAMES per second.\n");
      printf(" 's' - Show SKIN in the density shading.\n");
      printf(" 't' - Increment the THRESHHOLD.\n");
      printf(" 'T' - Decrement the THRESHHOLD.\n");
      printf(" 'n' - NORMAL shade\n");
      printf(" 'p' - Save out a PICTURE of the view\n");
      printf(" 'q' - QUIT\n");
      printf(" 'u' - ERODE the data\n");
      printf(" 'w' - Add WINDOWING of the threshhold by 16.\n");
      printf(" 'x' - Make a CROSS SECTION CUT in the model.\n");
      break;

    case 'i':
      volume->dilate();
      glutPostRedisplay();
      break;

    case 'j':
      volume->thresh(volume->threshhold);
      glutPostRedisplay();
      break;

    case 'k':
      mode = TIME;
      g_frames_count = 0.0;
      g_clock_start = clock();
      glutPostRedisplay();
      break;

    case 'n':
      volume->setRenderType(NORMAL);
      glutPostRedisplay();
      break;

    case 'p':
      volume->writeImage();
      break;

    case 'q':
      exit (0);
      break;
      
    case 'u':
      volume->erode();
      glutPostRedisplay();
      break;

    case 's':
      volume->toggleShowSkin();
      glutPostRedisplay();
      break;

    case 't':
      volume->setThresh( volume->threshhold+1, 255);
      glutPostRedisplay();
      break;

    case 'T':
      volume->setThresh( volume->threshhold-1, 255);
      glutPostRedisplay();
      break;

    case 'w':
      if (volume->threshhold < 240){
        volume->setThresh( volume->threshhold, volume->threshhold+16);
      }
      else{
        volume->setThresh( volume->threshhold, 255);
      }
      glutPostRedisplay();
      break;

    case 'x':
      if (CUT != mode)
        mode = CUT;
      else
        mode = NONE;
      glutPostRedisplay();
      break;
  }
}


//MAIN
//MAIN
//MAIN
//MAIN
//MAIN
/******************************************************************************
 * Purpose:  main... The head hancho function.
 *****************************************************************************/
int main(int argc, char** argv) {
  char name[255];
  /*
  int dummy;
  dummy = -8;
  printf("dummy: %d\n", dummy); 
  printf("dummy<<3: %d\n", dummy<<3); 
  printf("dummy>>1: %d\n", dummy>>1); 
  */
  if (argc < 2 || argc >3) {
    if (argc == 1){
      //-head_mr
      volume = new Volume(256, 256, 256, XDIM, YDIM, .02/.015);//.02, .02, .015
      volume->read("data/Head_256x256x256.raw");
    }
    else{
      printf("  Error:\nIncorrect command line arguments.  The format is:\n");
      printf(">vrend [data_file]\n");
      exit (1);
    }
  }
  else if (strcmp("-head_sm", argv[1]) == 0){
    volume = new Volume(256, 256, 106, XDIM, YDIM, .435547);//.435547, .435547, 1
    volume->read("data/CT_Head_small.raw");
  }
  else if (strcmp("-engine", argv[1]) == 0){
    volume = new Volume(256, 256, 110, XDIM, YDIM, 1.0);//1, 1, 1
    volume->read("data/Engine.raw");
  }
  else if (strcmp("-teddy", argv[1]) == 0){
    volume = new Volume(128, 128, 62, XDIM, YDIM, 2.8/5);//2.8, 2.8, 5
    volume->read("data/Teddybear.raw");
  }
  else if (strcmp("-tbone", argv[1]) == 0){
    volume = new Volume(512, 512, 72, XDIM, YDIM, .11718/.5);//.11718, .11718, .5
    volume->read("data/tbone.raw");
  }
  else if (strcmp("-bonsai", argv[1]) == 0){
    volume = new Volume(256, 256, 128, XDIM, YDIM, .585938/.5);//.585938, .585938, .5
    volume->read("data/Bonsai.raw");
  }
  else if (strcmp("-head_mr", argv[1]) == 0){
    volume = new Volume(256, 256, 256, XDIM, YDIM, .02/.015);//.02, .02, .015
    volume->read("data/Head_256x256x256.raw");
  }
  else{
    printf("Flag '%s' not supported\n", argv[1]);
    printf("Valid flags are:\n");
    printf("-head_sm\n");
    printf("-engine\n");
    printf("-teddy\n");
    printf("-tbone\n");
    printf("-bonsai\n");
    printf("-head_mr\n");
    exit(0);
  }

  glutInit(&argc, argv);
    //GLUT_DOUBLE used for double buffering
  glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE );
  glutInitWindowSize (volume->im_width, volume->im_height);
  glutInitWindowPosition (0, 0);
  strcpy(name, argv[0]);//name of the program
  strcat(name, " Volume Renderer");
  glutCreateWindow (name);
  init();
      // The names of the functions that you use for displaying and events go
      // as the parameters of this functions below
  glutDisplayFunc(display);
  glutReshapeFunc(reshape);
  glutMouseFunc(mouseButton);
  glutMotionFunc(mouseMotion);
  glutKeyboardFunc(keyboard);
  glutMainLoop();

  return 0;
}
//EOF