Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] View of /branches/lamont/test/MIP/mip_c.c
ViewVC logotype

View of /branches/lamont/test/MIP/mip_c.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2081 - (download) (as text) (annotate)
Mon Nov 5 23:26:06 2012 UTC (7 years ago) by lamonts
File size: 5501 byte(s)
Creating new developmented branch based on vis12
/*! \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_t	pos;
    Diderot_vec3_t	dir;
    Diderot_real_t	t;
    Diderot_real_t	maxval;
    Diderot_real_t	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; 
}

// thie C code for the actor's update method.  This is a general version that takes
// two self pointers; selfIn is the state at the beginning of the iteration and
// selfOut is the new state at the end of the step.  For this example, we could get
// away with just one self pointer, since there are no actor-actor interactions.
//
void RayCast_update (GlobalState_t *glob, RayCast_t *selfIn, RayCast_t *selfOut)
{
//    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_vec3_t dir = selfIn->dir;
    Diderot_vec3_t pos = selfIn->pos;
    pos = Diderot_AddV3(pos, Diderot_ScaleV3(glob->stepSz, dir.v));

    Diderot_Vec3_t imgPos = Diderot_TransformVec3 (pos, img->mInv);
    if (Diderot_Inside3D(imgPos, img)) {
	/* ?? probe ?? */
	Diderot_real_t maxval = selfIn->maxval;
	maxval = Diderot_Max(val, maxval);
	selfOut->maxval = maxval;
    }
    else {
	selfOut->maxval = selfIn->maxval;
    }

    Diderot_real_t t = selfIn->t;
    if (t > 20.0f) {
	selfOut->pos = pos;
	selfOut->t = t;
	RayCast_stabilize (glob, selfOut);
	return;
    }

    t = t + glob->stepSz;

    selfOut->pos = pos;
    selfOut->t = t;

}

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; 
} 

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0