/*! \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
} 

/* 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_OK;
}

Status_t Diderot_Inputf (const char *name, float *v, bool hasDflt)
{
    return DIDEROT_OK;
}

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