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

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

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

/* 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 *img)
{
    Nrrd *nin = loadNrrdFile (name);

    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, transformMatrix); 
    InvertM (transformMatrix, inverseMatrix); 

}

/* functions to get input-parameter values */
extern Status_t Diderot_InputString (const char *, const char **, bool);
extern Status_t Diderot_InputReal (const char *, Diderot_real_t *, bool);
extern Status_t Diderot_InputVec3 (const char *, Diderot_vec3_t *v);

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

Diderot_real_t DetM4x4 (Matrix4x4_t m) 
{ 
	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])); 

} 
void InvertM4x4 (Matrix4x4_t i, Matrix4x4_t m) 
{  
	Diderot_real_t scale = 1.0 / DetM4x4(m);

	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;

}