/*************************************************************************
 * main.c -- the sample program that takes an environment file in NETCDF,
 *           multiple time-step volume inputfiles either in NetCDF format 
 *           or on raw format, renders and generates output images in the
 *           format specified in the environment file.
 * 
 *	(C) Copyright 1998, California Institute of Technology
 *	ALL RIGHTS RESERVED.
 *      U.S. Government Sponsorship Acknowledged under NAS7-1260
 *
 *
 *  Author:  P. Li
 *  Creation Date: 1/7/98
 *  Revision History: NONE
 *
 *************************************************************************/
static char SccsId[] = "@(#)main.c	1.3 3/16/98";

#include <unistd.h>
#include <stdio.h>
#include <parvox.h>               /* parvox header file */

/* Constants */
#define MAX_FILENAME 80
extern long _MPP_MY_PE;
			      
/* Global variable for input arguments */
struct THE {
    int  index;                   /* variable index, if the input file contains
				     more than one variable */
    char input_file[MAX_FILENAME+1]; /* input file name */
    int t_start, t_stop, t_inc;	  /* start, stop, and increment timestep */
    int filetype;                 /* input file type, VL_NETCDF or VL_RAW */
    char environfile[MAX_FILENAME+1]; /* environment file name */
    char render_type[12];         /* render type, "vol", "iso", or "slice" */
} This, *The;

main(int argc, char** argv)
{
    char   c;
    extern int optind;
    extern char* optarg;
    int i, j, k, n, cnt;
    vlvolume 	*vol;               /* defined in vlvoxel.h */
    vlvariable 	*variable;          /* defined in vlvoxel.h */
    int offset, nid, time_step;
    vlmeta_voxel4d	voxel4d;    /* defined in vlmeta.h */
    vlmeta_environ      *environ;   /* defined in vlmeta.h */
    int index;
    int render_type;
   
    The = &This;

    /* Default values for the input argument structure */
    The->t_start = 0;
    The->t_stop = 0;
    The->t_inc = 1;
    The->index = 0;
    The->filetype = VL_NETCDF;
    strcpy(The->render_type,"vol");

    /******************************************************************* 
       USAGE: main  [-s begin,end,increment_timestep]
		    [-v index ]      - variable index to be rendered
		    [-t "vol"|"slice"|"iso"] - render type
		    [-r] 
		    environfile
		    inputfile		    
     *******************************************************************/  
			 
    while ((c = getopt(argc, argv, "v:s:t:r")) != '\377')
    {
	switch (c) 
	{	
	case 's':
	  sscanf(optarg, "%d,%d,%d", &The->t_start,&The->t_stop,&The->t_inc);
	  break;
	case 'v':
	  The->index = atoi(optarg);
	  break;
	case 't':
	  strcpy(The->render_type,optarg);
	  break;
	case 'r':
	  The->filetype = VL_RAW;
	  break;
	default:
	  break;
	}
    }
    if (argc > (optind+1)) {
        strncpy(The->environfile, argv[optind++], MAX_FILENAME);
	strncpy(The->input_file, argv[optind++], MAX_FILENAME);
    } else
    {
	vlerror("Insuffient Arguments");
    }

    /* initiate the render context */
    vlinit();

    /* read in the environment file and set up the rendering parameters,
       vlread_environ() calls vlinit_fb() implicitly */
    environ = vlread_environ(The->environfile);

    if (!environ) {
      vlerror("Error reading environment file %s\n", The->environfile);
    }

    /* read and distribute the multiple input data files based on the block 
       dimension defined in the environment file.  */
    /* dimensions here (or blk_per_pe) are specified by user */

    vol = vlread_block(The->input_file, The->t_start, 
		       The->t_stop, The->t_inc,
			environ->block_per_pe, environ->block_idim,
			environ->block_jdim, environ->block_kdim, 
		       (void *)&voxel4d, The->filetype);

    if (vol == NULL)
      vlerror("Error reading input file %s.\n", The->input_file);

    /* right now, only one variable per volume is implemented */
    index = The->index; 
    variable = &vol->vars[index];

    /* the following input-data related information are stored in environ 
       data structure, and has to be set explicitly after the input file is
       read in */

    /* 1. set the material property for all the variables */
    for (i = 0; i < environ->num_vars; i++) {
      vlmaterial(&vol->vars[i], VL_AMBIENT, environ->material[i].ka,0,0);
      vlmaterial(&vol->vars[i], VL_DIFFUSE, environ->material[i].kd,0,0);
      vlmaterial(&vol->vars[i], VL_SPECULAR, environ->material[i].ks,0,0);
      vlmaterial(&vol->vars[i], VL_SHINENESS, environ->material[i].shine,0,0);
    }
    
    /* 2. set nodata color, one for all the variables in the volume */
    vlset_nodata_color(vol, vol->nodata, 
		       environ->nodata_color.r, 
		       environ->nodata_color.g,
		       environ->nodata_color.b);

    /* check if proper definition is set for certain render type */
    /* i.e., at least one colormap and one opacitymap is set for volume,
             at least one slice defined for slice,
	     at least one surface defined for isosurface.
	     */
    if (!strncmp(The->render_type, "vol", 3)) {
      if (environ->num_colormap == 0 || environ->num_opacitymap == 0)
	vlerror("No colormap or opacity map is defined in the environ file");
      else 
	render_type = VL_VOLUME;
    } else if (!strncmp(The->render_type, "slice", 5)) {
      if (environ->num_slice[index] == 0)
	vlerror("No slices are defined in the environment file");
      else
	render_type = VL_SLICE;
    } else if (!strncmp(The->render_type, "iso", 3)) {
      if (environ->num_iso[index] == 0)
	vlerror("No isosurfaces are defined in the environment file");
      else
	render_type = VL_ISO;
    }
    
    /* 3. set render type for each time step of the specific variable to 
       be rendered */
    for (j = The->t_start; j <= The->t_stop; j+=The->t_inc)
      vlset_render_type(vol, variable, render_type, j-The->t_start);
      
    /* 4. Based on the render type, we need to classify the data sets 
       accordingly */

    switch (render_type) {
    case VL_VOLUME:
    case VL_SLICE:
      /* Volume and slice rendering requires colormap and opacity map been set
	 for each variable. */
      for (i = 0; i < environ->num_vars; i++) {
	vol->vars[i].colormapptr = environ->colormap[environ->map[i][0]];
	vol->vars[i].omapptr = environ->opacity[environ->map[i][1]];
      }
      /* Input data has to be normalized */
      for (j = The->t_start; j <= The->t_stop; j+=The->t_inc)
	vlnormalize_data(vol, variable, j-The->t_start);

      /* Pre-multiple the opacity to the color map for volume rendering */
      /* This is a necessary call to set up some ParVox internal stuff, we
	 need to make it transparent to user in our next release.  For now,
	 be sure to call either vlset_opmap() for volume or vlcopy_opmap()
	 for slice rendering after the colormap and opacity map pointers are
	 set */
      if (render_type == VL_VOLUME)
	vlset_opmap(variable, vol);
      else 
	vlcopy_opmap(variable,vol);

      /* Classification call for slice rendering */
      if (render_type == VL_SLICE) {
	for (j = The->t_start; j <= The->t_stop; j+=The->t_inc) 
	  vlset_slices(vol, variable, j-The->t_start, environ->slices, 
		       environ->num_slice[index]);
      }
      break;
    case VL_ISO:
      /* Enable surface first */
      for (i = 0; i < environ->num_iso[index]; i++) {
	vlsurface *surf = &(environ->surface[index][i]);

	/* if isosurface color is not set, pick it from the color table */
	/* the GUI doesn't have the information for the surface color,
	   therefore, the GUI didnot save the color */
	if (surf->col.r < 0) {
	  if (environ->colormap)
	    variable->colormapptr = environ->colormap[environ->map[index][0]];
	  vlpick_iso_color(variable, surf->threshold, &(surf->col));
	}
	vlenable_surface(variable, i);
      }
      /* Now, shade the input data with the isosurface information.  */
      for (i = The->t_start; i <= The->t_stop; i+=The->t_inc) {
	vliso_shade(vol, variable, i-The->t_start, environ->surface[index], 
		      environ->num_iso[index] );
      }
      break;
    }

    /* set the volume and the variable to be rendered. 
     */
    vlset_new_volume(vol);
    vlenable_variable(vol, index);
     
    /* Do the rendering,  you have to do matrix stuff for each time step, then 
       call vlflush().  Note that the matrix transformation information is 
       stored in the environment file in three forms, the transformation 
       matrix, the scaling factors, and the x, y tranlation for the convenience 
       of restoring the GUI setup. Therefore, we have to build the modelview
       matrix explicitly in the render loop using these three sets of data.
       Also note that the projection matrix, determined by vlortho(), is bounded 
       by the input data dimension and is not stored in the environment file.
    */

    cnt = 0;
    /* Render each time step */
    for (time_step = The->t_start; time_step <= The->t_stop;
    				   time_step+=The->t_inc, cnt++)
    {
        float max, depth, ratio;
	long t1, t2;

	/* aspect ratio */
	ratio = (float)environ->width/(float)environ->height;

#ifdef TIMING           /* Cray code for Timing */
	barrier();	/* sync up all procs to start at once */
	t1 = clock();
#endif

	/* Set up the modelview matrix */
	vlmatrix_mode(VL_MODELVIEW);
	vlload_identity();
	/* 1. Traslate to the center of the volume, taking into account of
	   the spacing and the scaling factors, spacing is defined in the
	   volume itself, and the scaling is defined in the environment
	   file */
	vltranslate(- vol->idim * .5 * vol->spacing.x*environ->scale[0],
		    - vol->jdim * .5 * vol->spacing.y*environ->scale[1],
		    - vol->kdim * .5 * vol->spacing.z*environ->scale[2]);

	/* 2, now, multiple the modelview matrix by the transformation matrix
	   defined in the environment file */
	vlmult_matrix(environ->trans_matrix);
	
	/* 3, last, traslate the matrix by the traslation vector defined in
	   the environment file */
	vltranslate(environ->translate[0], environ->translate[1], 0.0);
	/* 4. finally, select a look-at point at (0,0,5), gaze (0,0,0) and
	   up vection (0, -1,0).  Note that this is hard coded in netparvox,
	   it will results in the same viewing position as you see 
	   interactively */
	vllook_at(0, 0, 5, 0, 0, 0, 0, -1, 0);

	/* Set up the projection matrix */
	vlmatrix_mode(VL_PROJECTION);
	vlload_identity();
	max = (vol->idim > vol->jdim)? vol->idim: vol->jdim;
	depth = (max > vol->kdim)? max: vol->kdim;

	/* Define ortho, again, this is hard coded in netparvox.  The following
	   line will give you the same viewing position as you see 
	   interactively. */
	vlortho(-(max * .75) * ratio, max * .75 * ratio,
		    -max * .75, max * .75, .1, depth * 10);

	/* After everything is set, do it */
	vlflush(vol, time_step - The->t_start);

#ifdef TIMING	 /* Cray specific code for timing */
	barrier();
	t2 = rtclock();
	vl0printf("Done rendering in %f seconds.\n", (float)(t2-t1)/150000000.0);
#endif
    }
}