// Kevin Beason
// http://www.cs.fsu.edu/~beason/
// 4/16/02
// 
// This just reads in a dataset and shows an isosurface at a user-
// controlled isovalue.
// 
// The main innovation here is that now, instead of a function-defined
// data set, you can use one read in from a file.

#include <math.h>
#include <stdlib.h>
#include <iostream.h>
#include <fstream.h>
#include <string.h>
#include <Inventor/actions/SoWriteAction.h>
#include <Inventor/Xt/SoXt.h>
#include <Inventor/Xt/viewers/SoXtExaminerViewer.h>
#include <Inventor/nodes/SoCoordinate3.h>
#include <Inventor/nodes/SoCube.h>
#include <Inventor/nodes/SoCone.h>
#include <Inventor/nodes/SoMaterial.h>
#include <Inventor/nodes/SoSeparator.h>
#include <Inventor/nodes/SoSphere.h>
#include <Inventor/nodes/SoTransform.h>
#include <Inventor/actions/SoGLRenderAction.h>
#include <Inventor/nodes/SoPointLight.h>
#include <Inventor/sensors/SoTimerSensor.h>
#include <Inventor/nodes/SoComplexity.h>
#include <Inventor/nodes/SoCamera.h>
#include <Inventor/SoInput.h>
#include <Inventor/SoDB.h>
#include <Inventor/nodes/SoCylinder.h>
#include <Inventor/SbLinear.h>
#include <Inventor/nodes/SoShapeHints.h>
#include <Inventor/draggers/SoTranslate1Dragger.h>
#include <Inventor/engines/SoCalculator.h>
#include <Inventor/nodes/SoVertexProperty.h>
#include <Inventor/nodes/SoIndexedFaceSet.h>
#include <Inventor/fields/SoSubField.h>
#include <Inventor/SbLinear.h>

#include <stdio.h>
#include "MarchingCubes.cpp"
#include "intersect.cxx"
#include "scene.cxx"
#include "EmissiveSampleSet.cpp"

SoPointLight *light;
const int dim = 10000;
int numSamples = 0;
int rendSamplesMax = 100000;
int rendSamplesPerPass = 100;
int rendSamples = 0;
double myGamma = 0.15;
Distribution *sphereEmits[dim];
SoCamera *Camera;
double r, g, b;
Scene myScene;
EmissiveSampleSet *eset;
//EmissiveSampleSet *rend;
SoCylinder *myCyl;
SoMaterial *cylMat;
SoTransform *cylXform;
SoTransform *cylXform2;
SoIndexedFaceSet *faceSet;
SoVertexProperty *vprop;
SoTranslate1Dragger *dragger;

////////////////////////////////////////////////////////////////////////
// Begin code from Josh Grant
////////////////////////////////////////////////////////////////////////

void
setData (SFScalarField &sfield, int xDim, int yDim, int zDim)
{
  // set dimensions of mesh
  int index;
  int dims[3] = {xDim, yDim, zDim};
  dims[0] = 128;
  dims[1] = 128;
  dims[2] = 40;
  
  // allocate memory for data
  float *data = new float[dims[0]*dims[1]*dims[2]];

  // set the minimum and maximum bounds of the exponential function
  SbVec3f min(-1.0, -1.0, -1.0);
  SbVec3f max( 1.0,  1.0,  1.0);
  SbVec3f size = max - min;

  // compute the size of each voxel
  SbVec3f voxel, p;
  voxel[0] = (dims[0] > 1) ? size[0] / (float)(dims[0] - 1) : 0.0;
  voxel[1] = (dims[1] > 1) ? size[1] / (float)(dims[1] - 1) : 0.0;
  voxel[2] = (dims[2] > 1) ? size[2] / (float)(dims[2] - 1) : 0.0;

  // open file in binary
  char surfaceFile[] = "pumpkin4.data";
  FILE *dataFile = fopen(surfaceFile, "rb");
  if (!dataFile){
    cerr << "Could not open file: " << surfaceFile << endl;
    exit(1);
  }
  
  char a, b;
  int c, cmax=0, cmin=HUGE_VAL, numc=0;;
  double avgc=0.0;
  int valSize=1;
  unsigned char buffer[256*256*2];
  
  // iterate through each grid point
  index = 0;
  for (int k = 0; k < dims[2]; k++) {
    p[2] = min[2] + voxel[2]*k;

    // read in data
    cerr << "Reading in slice " << k << endl;
    fread(buffer, valSize, dims[0]*dims[1], dataFile);
    
    for (int j = 0; j < dims[1]; j++) {
      p[1] = min[1] + voxel[1]*j;
      for (int i = 0; i < dims[0]; i++) {
	p[0] = min[0] + voxel[0]*i;

	// compute data
	//c = buffer[j*dims[0]*valSize+i*valSize+0] << 8;
	//c = c | buffer[j*dims[0]*valSize+i*valSize+0];
	c = (unsigned)buffer[j*dims[0]*valSize+i*valSize+0];
	numc++;
	avgc+=(double)c;

	// diagnostics
	if      (c>cmax) cmax = c;
	else if (c<cmin) cmin = c;
	
	data[index] = (double)(c)/275.0;
	
	//data[index] = sin(p[0]*p[0]) + cos(p[1]*p[1]) + p[2]*p[2];
	
	index++;
      }
    }
  }

  avgc/=numc;
  cerr << "Max: " << cmax << "\nMin: " << cmin
       << "\nAvg:" << avgc << endl;
  fclose(dataFile);

  // set data, but don't copy the data since we already have a valid copy
  // allocated.  This is just an option provided by the SFScalarField class.
  // Setting it to true will have the same effect if just means the
  // SFScalarField class will allocate memory for the data and copy all the
  // values to its local space.
  sfield.setValue(dims, data, false);

}



////////////////////////////////////////////////////////////////////////
// End code from Josh Grant
////////////////////////////////////////////////////////////////////////

SoSeparator *
makeScene ()
  {

  SoSeparator *scene = new SoSeparator();
  int i, j, k;
  double x, y, z;
  //double l, w, h;
  double xMin=-1.0, xMax=1.0;
  double rad, angle;
  SbVec3f point, normal;
  myScene.readScene("scene.txt");
  
  // cylinder
  SoSeparator *cylSep = new SoSeparator();
  myCyl = new SoCylinder();
  cylMat = new SoMaterial();
  cylXform = new SoTransform();
  cylXform2 = new SoTransform();
  myCyl->radius.setValue(0.0);
  myCyl->height.setValue(0.0);
  cylMat->ambientColor.setValue(0.0, 0.0, 0.0);
  cylMat->diffuseColor.setValue(0.0, 0.0, 0.0);
  cylMat->specularColor.setValue(0.0, 0.0, 0.0);
  cylMat->emissiveColor.setValue(0.0, 0.0, 0.0);
  cylMat->shininess.setValue(0.0);
  cylMat->transparency.setValue(0.0);
  cylXform->translation.setValue(0.0, 0.0, 0.0);
  cylXform2->translation.setValue(0.0, 0.0, 0.0);
  cylSep->addChild(cylMat);
  cylSep->addChild(cylXform2);
  cylSep->addChild(cylXform);
  cylSep->addChild(myCyl);
  scene->addChild(cylSep);

  int sphereNum;
  eset = new EmissiveSampleSet();
  //rend = new EmissiveSampleSet();
  SoSeparator *esetSep = new SoSeparator();
  SoMaterial *esetMat = new SoMaterial();
  esetMat->ambientColor.setValue(0.0, 0.0, 0.0);
  esetMat->diffuseColor.setValue(0.0, 0.0, 0.0);
  esetMat->specularColor.setValue(0.0, 0.0, 0.0);
  esetMat->emissiveColor.setValue(0.0, 0.0, 0.0);
  esetMat->shininess.setValue(0.0);
  esetMat->transparency.setValue(0.0);
  esetSep->addChild(esetMat);
  esetSep->addChild(eset);
  eset->sampleRadius = 0.15;
  scene->addChild(esetSep);

  for (i=0; i<numSamples; i++){ // loop through angles
    
    //l = 100.0/numSamples; //drand48() *0.5;
    point = SbVec3f(drand48()*2.0-1.0,
		    drand48()*2.0-1.0,
		    drand48()*2.0-1.0);
    while (point.length()>1.0){
        point = SbVec3f(drand48()*2.0-1.0,
			drand48()*2.0-1.0,
			drand48()*2.0-1.0);
    }
    point.normalize();
    normal = point;
    sphereNum = (int)(drand48()*myScene.sphere_i);
    point*=myScene.spheres[sphereNum].radius;
    point+=myScene.spheres[sphereNum].center;
    
    // emissivesamplesset node
    eset->location.set1Value(i, point);
    eset->normal.set1Value(i, normal);
    eset->color.set1Value(i, SbColor(0.5, 0.5, 0.5));
    
    // sphere emittance
    sphereEmits[i] = myScene.spheres[sphereNum].emittance;

  }


////////////////////////////////////////////////////////////////////////  
// Begin Josh Grant's code
////////////////////////////////////////////////////////////////////////

  // create MarchingCubes engine and connect fields
  MarchingCubes *mcubes = new MarchingCubes();
  mcubes->ref();

  setData(mcubes->data, 15, 15, 15);

  // set isoValue, in this case this means find the surface of the sphere
  // where the radius is 1.0
  mcubes->isoValue = 1.0;

  // MarchingCubes gives the triangles with the vertices ordered clockwise
  SoShapeHints *hints = new SoShapeHints();
  hints->vertexOrdering.setValue(SoShapeHints::CLOCKWISE);
  scene->addChild(hints);

  // Create a VertexProperty because they are more efficient the
  // SoCoordinate3 nodes.
  vprop = new SoVertexProperty();
  vprop->ref();
  vprop->vertex.connectFrom(&mcubes->points);
  vprop->normal.connectFrom(&mcubes->normals);

  // connect the indexes from MarchingCubes
  SoSeparator *faceSep = new SoSeparator();
  SoMaterial *faceMat = new SoMaterial();
  faceMat->ambientColor.setValue(0.0, 0.0, 0.0);
  faceMat->diffuseColor.setValue(0.0, 1.0, 0.0);
  faceMat->specularColor.setValue(0.0, 0.0, 0.0);
  faceMat->emissiveColor.setValue(0.0, 0.3, 0.0);
  faceMat->shininess.setValue(0.0);
  faceMat->transparency.setValue(0.97);
  faceSet = new SoIndexedFaceSet();
  faceSet->vertexProperty = vprop;
  faceSet->coordIndex.connectFrom(&mcubes->indexes);
  faceSep->addChild(faceMat);
  faceSep->addChild(faceSet);
  scene->addChild(faceSep);

  // Translate the dragger by two in the y direction
  SoTransform *trans = new SoTransform();
  trans->translation.setValue(0, 2, 10);
  scene->addChild(trans);

  // Create the dragger and connect an SoCalculator to it.
  // Why?  This allows for a slick way to make the program interactive.
  // Whenever the dragger is moved the 'isoValue' field of the MarchingCubes
  // engine will be updated, as well as the isosurface it outputs.
  dragger = new SoTranslate1Dragger();
  SoCalculator *calc = new SoCalculator();
  // translate the dragger to <1, 0, 0>. 
  dragger->translation.setValue(1, 0, 0);
  // connect the dragger translation field to an input to the calculator
  calc->A.connectFrom(&dragger->translation);
  // Now set the expression for the calculator to evaluate.  This means
  // everytime the input vector field 'A' changes assign the x component of
  // the vector to the output variable 'a' of the calculator.  In terms of
  // this program, that means...everytime the dragger moves the current x
  // position of the dragger will be the new isoValue for the MarchingCubes
  // engine to compute.
  calc->expression.setValue("oa = A[0]");
  // connect the output of the engine to the isoValue field of MarchingCubes
  mcubes->isoValue.connectFrom(&calc->oa);

  scene->addChild(dragger);

////////////////////////////////////////////////////////////////////////  
// End Josh Grant's code
////////////////////////////////////////////////////////////////////////
  
  return scene;
  }

float PresentTime = 0.0;
float DTime = 0.015;
int moveIt=0;
double old_pos = 1.0;

void
animate(void *, SoSensor *)
{
  int i;
  double v;
  double r, g, b;
  SbVec3f x, y, point, normal;
  int sphereNum;

  return;
  
  // reset points if dragger moves
  if (old_pos != dragger->translation.getValue()[0]){
    rendSamples = 0;
    eset->location.setNum(numSamples);
  }
  old_pos = dragger->translation.getValue()[0];

  // move a sample
  moveIt++;
  if (moveIt>=numSamples)
    moveIt=0;
  point.setValue(1.0, 1.0, 1.0);
  while (point.length()>1.0){
    point = SbVec3f(drand48()*2.0-1.0,
		    drand48()*2.0-1.0,
		    drand48()*2.0-1.0);
  }
  point.normalize();
  normal = point;
  sphereNum = (int)(drand48()*myScene.sphere_i);
  point*=myScene.spheres[sphereNum].radius;
  point+=myScene.spheres[sphereNum].center;
  
  // emissivesamplesset node
  if (numSamples>0){
    eset->location.set1Value(moveIt, point);
    eset->normal.set1Value(moveIt, normal);
  }

  // sphere emittance
  sphereEmits[moveIt] = myScene.spheres[sphereNum].emittance;
  
  // set emittance
  for (i=0; i<numSamples; i++){
    
     x = (eset->normal.getValues(i))[0];
     y=Camera->position.getValue();
     y.normalize();
     v=x.dot(y);

     // calculate emission
     sphereEmits[i]->get_vals(M_PI/2.0-acos(v), r, g, b);
     eset->color.set1Value(i, SbColor(r, g, b));
     
  }

  // while rendMaxSamples not reached
  if (rendSamples<=rendSamplesMax){
  
    // shoot several rays
    for (i=0; i<rendSamplesPerPass; i++){
  
      // shoot a ray
      SbVec3f p0, dp, intersection, minIs;
      int minIsSn=0;
      double minDist = HUGE;
      dp.setValue(1.0, 1.0, 1.0);
      while (dp.length()>1.0){
	dp = SbVec3f(drand48()*2.0-1.0,
		     drand48()*2.0-1.0,
		     drand48()*2.0-1.0);
      }
      dp.normalize();
      p0.setValue(0.0, 3.0, 2.0);

      if (i==0){
	// draw cyclinder!!
	myCyl->height.setValue(100.0);
	myCyl->radius.setValue(0.1);
	cylMat->emissiveColor.setValue(0.5, 0.0, 0.5);
	cylMat->diffuseColor.setValue(0.5, 0.0, 0.5);
	cylMat->transparency.setValue(0.97);
	SbRotation sbr;
	sbr.setValue(SbVec3f(0.0, 1.0, 0.0), dp);
	cylXform->rotation = sbr;
	cylXform2->translation.setValue(p0+51.0*dp);
      }

      //////////////////////////////////////////////////////////////
      // Begin Triangle intersection test
      //////////////////////////////////////////////////////////////

      const int32_t *indices = faceSet->coordIndex.getValues(0);
      const SbVec3f *vertices = vprop->vertex.getValues(0);
      const SbVec3f *normals = vprop->normal.getValues(0);
      int numCoords = faceSet->coordIndex.getNum();
      int idx, idx0, idx1, idx2;
      double orig[3], dir[3], vert0[3], vert1[3], vert2[3], t, u, v, last_t;

      r = g = b = 0.0;
      last_t = HUGE;
      
      for (idx=0; idx<numCoords; idx+=4){

	orig[0] = p0[0];
	orig[1] = p0[1];
	orig[2] = p0[2];

	dir[0] = dp[0];
	dir[1] = dp[1];
	dir[2] = dp[2];

	idx0 = indices[idx+0];
	idx1 = indices[idx+1];
	idx2 = indices[idx+2];

	vert0[0] = vertices[idx0][0];
	vert0[1] = vertices[idx0][1];
	vert0[2] = vertices[idx0][2];

	vert1[0] = vertices[idx1][0];
	vert1[1] = vertices[idx1][1];
	vert1[2] = vertices[idx1][2];

	vert2[0] = vertices[idx2][0];
	vert2[1] = vertices[idx2][1];
	vert2[2] = vertices[idx2][2];


	if (intersect_triangle(orig, dir, vert0, vert1, vert2, &t, &u, &v)
	    && t>0.0){

	  if (t>last_t)
	    break;
	  last_t = t;
//int
//intersect_triangle(double orig[3], double dir[3],
//                   double vert0[3], double vert1[3], double vert2[3],
//                   double *t, double *u, double *v)

	  // draw emissive sample
	  // draw remote sample
	  eset->location.set1Value(numSamples+rendSamples, p0+dp*t);
	  eset->normal.set1Value(numSamples+rendSamples, normals[idx0]);

	  // calculate emission
	  x = normals[idx0];
	  x.normalize();
	  y=p0-(p0+dp*t);
	  y.normalize();
	  v=x.dot(y);

	  // set local sample color
	  //myScene.spheres[0].emittance->get_vals(M_PI/2.0-acos(v), r, g,b);
	  Scene::findDist("MY_EMITTANCE0").get_vals(M_PI/2.0-acos(v), r, g, b);
	  eset->color.set1Value(numSamples+rendSamples, SbColor(r, g, b));
	  rendSamples++;

	  // adjust cylinder
	  if (i==0){
	    // p0+50.0*dp
	    cylXform2->translation.setValue(( (t-1.0)/2.0+1.0 )*dp+p0);
	    myCyl->height.setValue(t-1.0);
	    cylMat->emissiveColor.setValue(0.0, 0.0, 0.0);
	    cylMat->diffuseColor.setValue(r, g, b);
	    cylMat->transparency.setValue(0.0);
	  }
	}

      }
      
      
      //////////////////////////////////////////////////////////////
      // End Triangle intersection test
      //////////////////////////////////////////////////////////////

      /*
      // find closest intersection
      int intFound = 0;
      for (sphereNum=0; sphereNum<myScene.sphere_i; sphereNum++){
	if (myScene.spheres[sphereNum].intersect(intersection, p0, dp)){
	  intFound = 1;
	  // if closest intersectoin
	  if ((intersection-p0).length()<minDist){
	    // save intersection info
	    minDist = (intersection-p0).length();
	    minIs = intersection;
	    minIsSn = sphereNum;
	  }
	}
      }

      if (intFound){
	// calculate emission
	x = minIs-myScene.spheres[minIsSn].center;
	x.normalize();
	y=p0-minIs;
	y.normalize();
	v=x.dot(y);

	// set local sample color
	myScene.spheres[minIsSn].emittance->get_vals(M_PI/2.0-acos(v), r, g, b);

	// draw remote sample
	eset->location.set1Value(numSamples+rendSamples, minIs);
	eset->normal.set1Value(numSamples+rendSamples, x);
	eset->color.set1Value(numSamples+rendSamples, SbColor(r, g, b));
	rendSamples++;

	// adjust cylinder
	if (i==0){
	  cylXform->translation.setValue((minIs.length()/2.0+1.0)*dp);
	  myCyl->height.setValue(minIs.length()-1.0);
	  cylMat->diffuseColor.setValue(r, g, b);
	}
      }
      else {
	// set local sample color
	r = g = b = 0.0;

	// set cylinder color
	if (i==0)
	  cylMat->diffuseColor.setValue(0.5, 0.5, 0.5);
      }
      */

      // draw local sample
      x = dp;
      x.normalize();
      eset->location.set1Value(numSamples+rendSamples, p0+x);
      eset->normal.set1Value(numSamples+rendSamples, x);
      eset->color.set1Value(numSamples+rendSamples, SbColor(r, g, b));
      rendSamples++;

    }

  }
  else
    cylMat->transparency.setValue(1.0);
  
  PresentTime += 3.14159/20.0; //DTime;
}

////////////////////////////////////////////////////////////////////////
//
// This code is from Dr. David C. Banks
//
////////////////////////////////////////////////////////////
//
// The main program gets a window from X, 
// puts an interactive 3D examinerViewer in it, 
// attaches the scene graph to the viewer,
// renders the scene graph,
// and hands control over to the Inventor main loop.

void
main(int argc, char* argv[])
{

  // read command line arguments
  int i, j;
  for (i=1; i<argc; i++){
    if (strcmp(argv[i], "-s") == 0){
      numSamples = atoi(argv[i+1]);
      if (numSamples>dim){
	numSamples = dim;
	cerr << "Too many samples. Using " << dim << ".\n";
      }
    }
    else if (strcmp(argv[i], "-g") == 0)
      myGamma = atof(argv[i+1]);
    else if (strcmp(argv[i], "-h") == 0){
      cerr << "\nUsage: " << argv[0]
	   << " [ -s samples ] [ -g gamma ]\n\n";
      exit(1);
    }
  }

    
  // Initialize Inventor and Xt
  Widget myWindow = SoXt::init(argv[0]);

  // Initialize new classes
  EmissiveSampleSet::initClass();
  SFScalarField::initClass();
  MarchingCubes::initClass();
  
  if (myWindow == NULL) exit(1);

  // SoGLRenderAction *rend = new SoGLRenderAction();
  //rend->setTransparencyType(SoGLRenderAction::SORTED_OBJECT_ADD);
  
  // Create a scene graph
  SoSeparator *root = makeScene();

  SoTimerSensor *sensor = new SoTimerSensor(animate, root);
  //  sensor->setInterval(SbTime(0.015)); // 70 frames per second
  sensor->setInterval(SbTime(0.015)); // 70 frames per second
  sensor->schedule();
  
  SoXtExaminerViewer *myViewer = new SoXtExaminerViewer(myWindow);
  myViewer->setSceneGraph(root);
  Camera = myViewer->getCamera();
  myViewer->setTitle("Isosurface Viewer");
  myViewer->setTransparencyType(SoGLRenderAction::SORTED_OBJECT_BLEND);
  myViewer->setBackgroundColor(SbColor(0.035, 0.035, 0.035));
  myViewer->show();

  SoWriteAction wa;  // Write out the scene graph just for demo purposes 
  //wa.apply(root);    // The scenegraph can be displayed using
		     // ivview, SceneViewer, or gview .
		     // It can also be displayed in a Web browser by
		     // changing its header from
		     //
		     // #Inventor V2.1 ascii
		     // 
		     // to
		     //
		     // #VRML V1.0 ascii
		     //
		     // to exploit the (deprecated) original version
		     // of the Virtual Reality Modelling Language
		     // file format.
		     // 
		     // Visit http://www.vrml.org/ (deprecated URL) or
		     // http://www.web3d.org/ for info on 
		     // VRML begat Web3D begat X3D.


  SoXt::show(myWindow);
  SoXt::mainLoop();
  }