/*! \file diderot.c
 *
 * \author John Reppy
 */

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

#include "Diderot/diderot.h"
#include <teem/nrrd.h>

/* Transformation matrix operations */

STATIC_INLINE Diderot_real_t DetM3x3 (
    Diderot_real_t a, Diderot_real_t b, Diderot_real_t c,
    Diderot_real_t d, Diderot_real_t e, Diderot_real_t f,
    Diderot_real_t g, Diderot_real_t h, Diderot_real_t i) 
{
   return ( (a)*(e)*(i) 
	  + (d)*(h)*(c) 
	  + (g)*(b)*(f) 
	  - (g)*(e)*(c) 
	  - (d)*(b)*(i) 
	  - (a)*(h)*(f)); 
}

static Diderot_real_t DetM4x4 (Diderot_Mat4x4_t m) 
{
#define M(i)	m[(i)>>2].r[(i)&3]
	return	(M( 0) * DetM3x3(M( 5), M( 6), M( 7), 
				 M( 9), M(10), M(11), 
				 M(13), M(14), M(15)) 						   
	       - M( 1) * DetM3x3(M( 4), M( 6), M( 7), 
				 M( 8), M(10), M(11), 
				 M(12), M(14), M(15)) 
	       + M( 2) * DetM3x3(M( 4), M( 5), M( 7), 
				 M( 8), M( 9), M(11), 
				 M(12), M(13), M(15)) 						   
	       - M( 3) * DetM3x3(M( 4), M( 5), M( 6), 
				 M( 8), M( 9), M(10), 
				 M(12), M(13), M(14))); 
#undef M

}

/*! \brief compute the inverse of \arg m, storing the result in \arg i.
 *  \param m the matrix to invert
 *  \param i the inverted matrix
 */
static void InvertM4x4 (Diderot_Mat4x4_t i, Diderot_Mat4x4_t m) 
{  
	Diderot_real_t scale = 1.0 / DetM4x4(m);

#define M(i)	m[(i)>>2].r[(i)&3]
#define I(j)	i[(j)>>2].r[(j)&3]
	I( 0) = DetM3x3(M( 5),M( 6),M( 7),
		       M( 9),M(10),M(11),
		       M(13),M(14),M(15)) * scale;
	
	I( 1) = -DetM3x3(M( 1),M( 2),M( 3),
			 M( 9),M(10),M(11),
			 M(13),M(14),M(15)) * scale;
					   
	I( 2) =  DetM3x3(M( 1),M( 2),M( 3),
			 M( 5),M( 6),M( 7),
			 M(13),M(14),M(15)) * scale;
					   
	I( 3) = -DetM3x3(M( 1),M( 2),M( 3),
			 M( 5),M( 6),M( 7),
			 M( 9),M(10),M(11)) * scale;
					   
	I( 4) = -DetM3x3(M( 4),M( 6),M( 7),
			 M( 8),M(10),M(11),
			 M(12),M(14),M(15)) * scale;
					   
	I( 5) =  DetM3x3(M( 0),M( 2),M( 3),
			 M( 8),M(10),M(11),
			 M(12),M(14),M(15)) * scale;
					   
	I( 6) = -DetM3x3(M( 0),M( 2),M( 3),
			 M( 4),M( 6),M( 7),
			 M(12),M(14),M(15)) * scale;
					   
	I( 7) =  DetM3x3(M( 0),M( 2),M( 3),
			 M( 4),M( 6),M( 7),
			 M( 8),M(10),M(11)) * scale;
					   
	I( 8) =  DetM3x3(M( 4),M( 5),M( 7),
			 M( 8),M( 9),M(11),
			 M(12),M(13),M(15)) * scale;
					   
	I( 9) = -DetM3x3(M( 0),M( 1),M( 3),
			 M( 8),M( 9),M(11),
			 M(12),M(13),M(15)) * scale;
					   
	I(10) =  DetM3x3(M( 0),M( 1),M( 3),
			 M( 4),M( 5),M( 7),
			 M(12),M(13),M(15)) * scale;
					   
	I(11) = -DetM3x3(M( 0),M( 1),M( 3),
			 M( 4),M( 5),M( 7),
			 M( 8),M( 9),M(11)) * scale;
					   
	I(12) = -DetM3x3(M( 4),M( 5),M( 6),
			 M( 8),M( 9),M(10),
			 M(12),M(13),M(14)) * scale;
					   
	I(13) =  DetM3x3(M( 0),M( 1),M( 2),
			 M( 8),M( 9),M(10),
			 M(12),M(13),M(14)) * scale;
					   
	I(14) = -DetM3x3(M( 0),M( 1),M( 2),
			 M( 4),M( 5),M( 6),
			 M(12),M(13),M(14)) * scale;
					   
	I(15) =  DetM3x3(M( 0),M( 1),M( 2),
			 M( 4),M( 5),M( 6),
			 M( 8),M( 9),M(10)) * scale;
#undef M
#undef I
} 

/*! \brief compute the inverse transpose of \arg m, storing the result in \arg i.
 *  \param m the matrix to invert
 *  \param i the inverted matrix
 */
static void InvertTransposeM4x4 (Diderot_Mat4x4_t i, Diderot_Mat4x4_t m) 
{  
	Diderot_real_t scale = 1.0 / DetM4x4(m);

#define M(i)	m[(i)>>2].r[(i)&3]
#define I(j)	i[(j)>>2].r[(j)&3]
	I( 0) = DetM3x3(M( 5),M( 6),M( 7),
		       M( 9),M(10),M(11),
		       M(13),M(14),M(15)) * scale;
	
	I( 4) = -DetM3x3(M( 1),M( 2),M( 3),
			 M( 9),M(10),M(11),
			 M(13),M(14),M(15)) * scale;
					   
	I( 8) =  DetM3x3(M( 1),M( 2),M( 3),
			 M( 5),M( 6),M( 7),
			 M(13),M(14),M(15)) * scale;
					   
	I(12) = -DetM3x3(M( 1),M( 2),M( 3),
			 M( 5),M( 6),M( 7),
			 M( 9),M(10),M(11)) * scale;
					   
	I( 1) = -DetM3x3(M( 4),M( 6),M( 7),
			 M( 8),M(10),M(11),
			 M(12),M(14),M(15)) * scale;
					   
	I( 5) =  DetM3x3(M( 0),M( 2),M( 3),
			 M( 8),M(10),M(11),
			 M(12),M(14),M(15)) * scale;
					   
	I( 9) = -DetM3x3(M( 0),M( 2),M( 3),
			 M( 4),M( 6),M( 7),
			 M(12),M(14),M(15)) * scale;
					   
	I( 7) =  DetM3x3(M( 0),M( 2),M( 3),
			 M( 4),M( 6),M( 7),
			 M( 8),M(10),M(11)) * scale;
					   
	I( 2) =  DetM3x3(M( 4),M( 5),M( 7),
			 M( 8),M( 9),M(11),
			 M(12),M(13),M(15)) * scale;
					   
	I( 6) = -DetM3x3(M( 0),M( 1),M( 3),
			 M( 8),M( 9),M(11),
			 M(12),M(13),M(15)) * scale;
					   
	I(10) =  DetM3x3(M( 0),M( 1),M( 3),
			 M( 4),M( 5),M( 7),
			 M(12),M(13),M(15)) * scale;
					   
	I(14) = -DetM3x3(M( 0),M( 1),M( 3),
			 M( 4),M( 5),M( 7),
			 M( 8),M( 9),M(11)) * scale;
					   
	I( 3) = -DetM3x3(M( 4),M( 5),M( 6),
			 M( 8),M( 9),M(10),
			 M(12),M(13),M(14)) * scale;
					   
	I( 7) =  DetM3x3(M( 0),M( 1),M( 2),
			 M( 8),M( 9),M(10),
			 M(12),M(13),M(14)) * scale;
					   
	I(11) = -DetM3x3(M( 0),M( 1),M( 2),
			 M( 4),M( 5),M( 6),
			 M(12),M(13),M(14)) * scale;
					   
	I(15) =  DetM3x3(M( 0),M( 1),M( 2),
			 M( 4),M( 5),M( 6),
			 M( 8),M( 9),M(10)) * scale;

#undef M
#undef I
} 

static void LoadTransformMatrix (Nrrd *nin, Diderot_Mat4x4_t m) 
{ 
#define M(i)	m[(i)>>2].r[(i)&3]
    int size =  nin->spaceDim; 

  // Image axis Scaling and Rotation 
    for (int i = 0; i < size; i++) { 
	NrrdAxisInfo axisInfo = nin->axis[i];
	for (int j = 0; j < size; j++) {
	    int k = (size + 1) * j + i;
	    M(k) = axisInfo.spaceDirection[j]; 
	}
      
      //Image Location 
	M((i * (size + 1)) + size) = nin->spaceOrigin[i];  
	
      //Bottom row of the Transform Matrix 
	M(((size + 1) * (size)) + i) = 0; 
    }
    M(((size + 1) * (size)) + size) = 1;

#undef M
}

static Nrrd *LoadNrrdFile (const char *filename) 
{
  /* create a nrrd; at this point this is just an empty container */
    Nrrd *nin = nrrdNew();

  /* read in the nrrd from file */
    if (nrrdLoad(nin, filename, NULL)) {
	char *err = biffGetDone(NRRD);
    	fprintf (stderr, "Diderot: trouble reading \"%s\":\n%s", filename, err);
    	free (err);
	return NULL;
    }
         
    return nin; 
  
}

/* load image data from Nrrd files */
//Status_t Diderot_LoadImage1D (Diderot_string_t name, Diderot_image1D_t *img);
//Status_t Diderot_LoadImage2D (Diderot_string_t name, Diderot_image2D_t *img);

Status_t Diderot_LoadImage3D (Diderot_string_t name, Diderot_image3D_t **imgOut)
{
    Nrrd *nin = LoadNrrdFile (name);
    Diderot_image3D_t *img = (Diderot_image3D_t *)malloc(sizeof(Diderot_image3D_t));

    img->dim = 3;
    img->size[0] = nin->axis[0].size;
    img->size[1] = nin->axis[1].size;
    img->size[2] = nin->axis[2].size;
    img->data = nin->data;

    LoadTransformMatrix (nin, img->m); 
    InvertM4x4 (img->m, img->mInv); 
    InvertTransposeM4x4 (img->m, img->mInvT); 

    *imgOut = img;
    return DIDEROT_OK;
}

/* functions to get input-parameter values */
Status_t Diderot_InputString (const char *name, const char **v, bool hasDflt)
{
    return DIDEROT_FAIL;
}

Status_t Diderot_Inputf (const char *name, float *v, bool hasDflt)
{
    float	f;

    while (true) {
	if (hasDflt)
	    printf("Enter value for %s (default %f): ", name, *v);
	else
	    printf("Enter value for %s: ", name);
	fflush (stdout);

	int n = scanf("%f\n", &f);
printf("n = %d\n", n);
	if (n == EOF)
	    return DIDEROT_FAIL;
	else if (n == 1) {
	    *v = f;
	    return DIDEROT_OK;;
	}
	else if (hasDflt)
	    return DIDEROT_OK;;
    }

}

Status_t Diderot_InputVec3f (const char *name, vec3f_t *v, bool hasDflt)
{
    return DIDEROT_FAIL;
}