Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

[diderot] View of /branches/pure-cfg/src/lib/diderot.c
 [diderot] / branches / pure-cfg / src / lib / diderot.c

# View of /branches/pure-cfg/src/lib/diderot.c

Thu Mar 10 15:48:09 2011 UTC (8 years, 10 months ago) by jhr
File size: 8509 byte(s)
```  modified debug printing
```
```/*! \file diderot.c
*
* \author John Reppy
*/

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

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

static void PrintMat4x4 (const char *title, Diderot_Mat4x4_t m)
{
#define M(i)	m[(i)>>2].r[(i)&3]
printf ("%s @ %p:\n", title, m);
for (int i = 0;  i < 4;  i++) {
int j = 4*i;
printf ("  [ %6f  %6f  %6f  %6f ]\n", M(j+0), M(j+1), M(j+2), M(j+3));
}
#undef M

}

/* 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 m, Diderot_Mat4x4_t i)
{
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 m, Diderot_Mat4x4_t i)
{
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 */
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)
{
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;

printf("LoadImage \"%s\": dim = 3, axes = <%d, %d, %d>\n",
name, img->size[0], img->size[1], img->size[2]);
InvertM4x4 (img->m, img->mInv);
InvertTransposeM4x4 (img->m, img->mInvT);
PrintMat4x4("m", img->m);
PrintMat4x4("mInv", img->mInv);
PrintMat4x4("mInvT", 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)
{
char	buf[256];
float	f;

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

char *p = fgets(buf, sizeof(buf), stdin);
if (p == NULL)
return DIDEROT_FAIL;	// EOF
int n = sscanf(buf, "%f\n", &f);
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)
{
char	buf[256];
float	f1, f2, f3;

while (true) {
if (hasDflt)
printf("Enter value for %s (default %f %f %f): ",
name, ((union3f_t)*v).r[0], ((union3f_t)*v).r[1], ((union3f_t)*v).r[2]);
else
printf("Enter value for %s: ", name);
fflush (stdout);

char *p = fgets(buf, sizeof(buf), stdin);
if (p == NULL)
return DIDEROT_FAIL;	// EOF
int n = sscanf(buf, "%f %f %f\n", &f1, &f2, &f3);
if (n == 3) {
*v = vec3f(f1, f2, f3);
return DIDEROT_OK;;
}
else if (hasDflt)
return DIDEROT_OK;;
}

}
```