/*! \file mip_c.c
 *
 * \author John Reppy
 *
 * This is a "by hand" translation of the vr-MIP.diderot example into parallel C code that
 * uses SSE instructions and pthreads.
 */

/*
 * COPYRIGHT (c) 2010 The Diderot Project (http://diderot-language.cs.uchicago.edu)
 * All rights reserved.
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <stdint.h>
#include <sys/sysctl.h>
#include <sys/stat.h>
#include <teem/nrrd.h>

/* SSE vector types */
typedef float float4 __attribute__ ((__vector_size__ (16)));
typedef double double2 __attribute__ ((__vector_size__ (16)));

/* union types for converting extracting vector data */
typedef union {
    float	f[4];
    float4	v;
} f4union;

typedef union {
    double	d[2];
    double2	v;
} d2union;

/* typedefs for Diderot types */
typedef int32_t Diderot_int;
typedef float Diderot_real;
typedef f4union Diderot_vec3;	// padded to fit in SSE register
typedef f4union Diderot_vec4;
typedef const char *Diderot_string;
typedef float 

/* global state representation */
typedef struct {
    Diderot_string	dataFile; // name of dataset
    Diderot_real	stepSz;	// size of steps
				// e.g. 0.1
    Diderot_vec3	eye;	// location of eye point
				// e.g. txs: (25,15,10)
                                //   vox1,x: (-8,2,2)
                                //   vox1,y: (2,-8,2)
                                //   vox1,z: (2,2,-8)
                                //  vfrhand: (127.331,-1322.05,272.53)
    Diderot_vec3	orig;	// location of pixel (0,0)
				// e.g. txs: (8.83877,2.5911,7.65275)  
				//   vox1,x: (0,3.4036,3.4036)
				//   vox1,y: (0.596402,0,3.4036)
				//   vox1,z: (3.4036,3.4036,0)
                                //  vfrhand: (122.835,17.7112,188.044)
    Diderot_vec3	cVec;	// vector between pixels horizontally
				// e.g. txs: (-0.0151831,0.0278357,0)
				//   vox1,x: (0,-0.014036,0)
				//   vox1,y: (0.014036,0,0)
				//   vox1,z: (-0.014036,0,0)
                                //  vfrhand: (-0.00403611,-0.029826,-0.244066)
    Diderot_vec3	rVec;	// vector between pixels vertically
				// e.g. txs: (0.0074887,0.00408474,-0.0305383)
				//   vox1,x: (0,0,-0.014036)
				//   vox1,y: (0,0,-0.014036)
				//   vox1,z: (0,-0.014036,0)
                                //  vfrhand: (-0.245595,-0.0112916,0.00544129)
    Diderot_image3Df	*img;
} GlobalState_t;

/* actor state representation */
typedef struct {
    Diderot_vec3	pos;
    Diderot_vec3	dir;
    Diderot_real	t;
    Diderot_real	maxval;
    Diderot_real	padding[2];
} RayCast_t;

void RayCast_init (GlobalState_t *glob, RayCast_t *self, int row, int col)
{
    self->pos = orig + real(row)*rVec + real(col)*cVec;
    self->dir = (pos - eye)/|pos - eye|;
    self->t = 0.0f;
    self->maxval = M_INF; 
}

void RayCast_update (GlobalState_t *glob, RayCast_t *self)
{
  // shadow state variables
    Diderot_vec3 pos, dir;
    Diderot_real t, maxval;

//    pos = pos + stepSz*dir;
//    if (inside (pos,F)) {
//      real val = F@pos;
//      maxval = max(val, maxval);
//    }
//    if (t > 20.0)
//	stabilize;
//    t = t + stepSz;

    Diderot_CopyV3 (&dir, &(glob->dir));
    Diderot_vec3 _t1;
    Diderot_ScaleV3 (&_t1, glob->stepSz,, &dir);
    Diderot_CopyV3 (&pos, &(glob->pos));
    Diderot_AddV3 (&pos, &pos, &_t);

/* ??? */
    if (/* ?? inside ?? */) {
	/* ?? probe ?? */
	maxval = self->maxval;
	maxval = Diderot_max(val, maxval);
	self->maxval = maxval;
    }

    t = self->t;
    if (t > 20.0f) {
	Diderot_CopyV3 (&(glob->pos), &pos);
	RayCast_stabilize (glob, self);
	return;
    }

    t = t + glob->stepSz;

    Diderot_CopyV3 (&(glob->pos), &pos);
    self->t = t;
    self->maxval = maxval;

}

void Global_init (GlobalState_t *glob)
{
    glob->stepSz = 0.1f;  // default value

    if ((Diderot_InputString("dataFile", &(glob->dataFile), false) != DIDEROT_OK)
    || (Diderot_InputReal("stepSz", &(glob->stepSz), true) != DIDEROT_OK)
    || (Diderot_InputVec3("eye", &(glob->eye), false) != DIDEROT_OK)
    || (Diderot_InputVec3("orig", &(glob->orig), false) != DIDEROT_OK)
    || (Diderot_InputVec3("cVec", &(glob->cVec), false) != DIDEROT_OK)
    || (Diderot_InputVec3("rVec", &(glob->rVec), false) != DIDEROT_OK))
      // error processing inputs
	exit (1);

  // image(3)[] img = load (dataFile);
    if (Diderot_LoadImage3D(glob->dataFile, &(glob->img)))
      // error loading image data
	exit (1);

}


int main (int argc, char ** argv) 
{
    GlobalState_t glob;

    glob.argc = argc;
    glob.argv = argv;
    Global_init (glob);

/* initial actors */

/* run simulation */

/* output */

    char *dataFile = "../../data/txs.nrrd"; 
    float transformMatrix[16]; 
    float inverseMatrix[16]; 

    if (argc == 2) {
	dataFile = argv[1];
    }

    float *out = (float *) malloc(sizeof(float) * (SIZE * SIZE)); 
    Nrrd * nin = loadNrrdFile (dataFile); 

    int sAxis[] = {nin->axis[0].size, nin->axis[1].size,nin->axis[2].size};  
    
    loadTransformMatrix(nin,transformMatrix); 
    invMatrix(transformMatrix,inverseMatrix); 
    

    exe_MIP_Kernel((float *)nin->data, (int)nrrdElementNumber(nin),inverseMatrix, sAxis, out);  
    
    return 0; 
}