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

SCM Repository

[diderot] Annotation of /trunk/test/MIP/mip_c.c
ViewVC logotype

Annotation of /trunk/test/MIP/mip_c.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3349 - (view) (download) (as text)

1 : jhr 434 /*! \file mip_c.c
2 :     *
3 :     * \author John Reppy
4 :     *
5 :     * This is a "by hand" translation of the vr-MIP.diderot example into parallel C code that
6 :     * uses SSE instructions and pthreads.
7 :     */
8 :    
9 :     /*
10 : jhr 3349 * This code is part of the Diderot Project (http://diderot-language.cs.uchicago.edu)
11 :     *
12 :     * COPYRIGHT (c) 2015 The University of Chicago
13 : jhr 434 * All rights reserved.
14 :     */
15 :    
16 :     #include <stdio.h>
17 :     #include <stdlib.h>
18 :     #include <stdbool.h>
19 :     #include <stdint.h>
20 :     #include <sys/sysctl.h>
21 :     #include <sys/stat.h>
22 :     #include <teem/nrrd.h>
23 :    
24 :     /* SSE vector types */
25 :     typedef float float4 __attribute__ ((__vector_size__ (16)));
26 :     typedef double double2 __attribute__ ((__vector_size__ (16)));
27 :    
28 :     /* union types for converting extracting vector data */
29 :     typedef union {
30 :     float f[4];
31 :     float4 v;
32 :     } f4union;
33 :    
34 :     typedef union {
35 :     double d[2];
36 :     double2 v;
37 :     } d2union;
38 :    
39 :     /* typedefs for Diderot types */
40 :     typedef int32_t Diderot_int;
41 :     typedef float Diderot_real;
42 :     typedef f4union Diderot_vec3; // padded to fit in SSE register
43 :     typedef f4union Diderot_vec4;
44 :     typedef const char *Diderot_string;
45 :     typedef float
46 :    
47 :     /* global state representation */
48 :     typedef struct {
49 :     Diderot_string dataFile; // name of dataset
50 :     Diderot_real stepSz; // size of steps
51 :     // e.g. 0.1
52 :     Diderot_vec3 eye; // location of eye point
53 :     // e.g. txs: (25,15,10)
54 :     // vox1,x: (-8,2,2)
55 :     // vox1,y: (2,-8,2)
56 :     // vox1,z: (2,2,-8)
57 :     // vfrhand: (127.331,-1322.05,272.53)
58 :     Diderot_vec3 orig; // location of pixel (0,0)
59 :     // e.g. txs: (8.83877,2.5911,7.65275)
60 :     // vox1,x: (0,3.4036,3.4036)
61 :     // vox1,y: (0.596402,0,3.4036)
62 :     // vox1,z: (3.4036,3.4036,0)
63 :     // vfrhand: (122.835,17.7112,188.044)
64 :     Diderot_vec3 cVec; // vector between pixels horizontally
65 :     // e.g. txs: (-0.0151831,0.0278357,0)
66 :     // vox1,x: (0,-0.014036,0)
67 :     // vox1,y: (0.014036,0,0)
68 :     // vox1,z: (-0.014036,0,0)
69 :     // vfrhand: (-0.00403611,-0.029826,-0.244066)
70 :     Diderot_vec3 rVec; // vector between pixels vertically
71 :     // e.g. txs: (0.0074887,0.00408474,-0.0305383)
72 :     // vox1,x: (0,0,-0.014036)
73 :     // vox1,y: (0,0,-0.014036)
74 :     // vox1,z: (0,-0.014036,0)
75 :     // vfrhand: (-0.245595,-0.0112916,0.00544129)
76 :     Diderot_image3Df *img;
77 :     } GlobalState_t;
78 :    
79 :     /* actor state representation */
80 :     typedef struct {
81 : jhr 446 Diderot_vec3_t pos;
82 :     Diderot_vec3_t dir;
83 :     Diderot_real_t t;
84 :     Diderot_real_t maxval;
85 :     Diderot_real_t padding[2];
86 : jhr 434 } RayCast_t;
87 :    
88 :     void RayCast_init (GlobalState_t *glob, RayCast_t *self, int row, int col)
89 :     {
90 :     self->pos = orig + real(row)*rVec + real(col)*cVec;
91 :     self->dir = (pos - eye)/|pos - eye|;
92 :     self->t = 0.0f;
93 :     self->maxval = M_INF;
94 :     }
95 :    
96 : jhr 447 // thie C code for the actor's update method. This is a general version that takes
97 :     // two self pointers; selfIn is the state at the beginning of the iteration and
98 :     // selfOut is the new state at the end of the step. For this example, we could get
99 :     // away with just one self pointer, since there are no actor-actor interactions.
100 :     //
101 :     void RayCast_update (GlobalState_t *glob, RayCast_t *selfIn, RayCast_t *selfOut)
102 : jhr 434 {
103 : jhr 438 // pos = pos + stepSz*dir;
104 :     // if (inside (pos,F)) {
105 :     // real val = F@pos;
106 :     // maxval = max(val, maxval);
107 :     // }
108 :     // if (t > 20.0)
109 :     // stabilize;
110 :     // t = t + stepSz;
111 : jhr 439
112 : jhr 447 Diderot_vec3_t dir = selfIn->dir;
113 :     Diderot_vec3_t pos = selfIn->pos;
114 : jhr 448 pos = Diderot_AddV3(pos, Diderot_ScaleV3(glob->stepSz, dir.v));
115 : jhr 439
116 : jhr 446 Diderot_Vec3_t imgPos = Diderot_TransformVec3 (pos, img->mInv);
117 :     if (Diderot_Inside3D(imgPos, img)) {
118 : jhr 439 /* ?? probe ?? */
119 : jhr 447 Diderot_real_t maxval = selfIn->maxval;
120 : jhr 446 maxval = Diderot_Max(val, maxval);
121 : jhr 447 selfOut->maxval = maxval;
122 : jhr 439 }
123 : jhr 447 else {
124 :     selfOut->maxval = selfIn->maxval;
125 :     }
126 : jhr 439
127 : jhr 447 Diderot_real_t t = selfIn->t;
128 : jhr 439 if (t > 20.0f) {
129 : jhr 447 selfOut->pos = pos;
130 : jhr 448 selfOut->t = t;
131 :     RayCast_stabilize (glob, selfOut);
132 : jhr 439 return;
133 :     }
134 :    
135 :     t = t + glob->stepSz;
136 :    
137 : jhr 447 selfOut->pos = pos;
138 :     selfOut->t = t;
139 : jhr 439
140 : jhr 434 }
141 :    
142 : jhr 437 void Global_init (GlobalState_t *glob)
143 :     {
144 : jhr 438 glob->stepSz = 0.1f; // default value
145 :    
146 :     if ((Diderot_InputString("dataFile", &(glob->dataFile), false) != DIDEROT_OK)
147 :     || (Diderot_InputReal("stepSz", &(glob->stepSz), true) != DIDEROT_OK)
148 :     || (Diderot_InputVec3("eye", &(glob->eye), false) != DIDEROT_OK)
149 :     || (Diderot_InputVec3("orig", &(glob->orig), false) != DIDEROT_OK)
150 :     || (Diderot_InputVec3("cVec", &(glob->cVec), false) != DIDEROT_OK)
151 :     || (Diderot_InputVec3("rVec", &(glob->rVec), false) != DIDEROT_OK))
152 : jhr 437 // error processing inputs
153 :     exit (1);
154 :    
155 :     // image(3)[] img = load (dataFile);
156 : jhr 438 if (Diderot_LoadImage3D(glob->dataFile, &(glob->img)))
157 :     // error loading image data
158 :     exit (1);
159 : jhr 437
160 :     }
161 :    
162 :    
163 : jhr 434 int main (int argc, char ** argv)
164 : jhr 437 {
165 :     GlobalState_t glob;
166 :    
167 :     glob.argc = argc;
168 :     glob.argv = argv;
169 :     Global_init (glob);
170 :    
171 :     /* initial actors */
172 :    
173 :     /* run simulation */
174 :    
175 :     /* output */
176 :    
177 : jhr 434 char *dataFile = "../../data/txs.nrrd";
178 :     float transformMatrix[16];
179 :     float inverseMatrix[16];
180 :    
181 :     if (argc == 2) {
182 :     dataFile = argv[1];
183 :     }
184 :    
185 :     float *out = (float *) malloc(sizeof(float) * (SIZE * SIZE));
186 : jhr 437 Nrrd * nin = loadNrrdFile (dataFile);
187 : jhr 434
188 :     int sAxis[] = {nin->axis[0].size, nin->axis[1].size,nin->axis[2].size};
189 :    
190 :     loadTransformMatrix(nin,transformMatrix);
191 :     invMatrix(transformMatrix,inverseMatrix);
192 :    
193 :    
194 :     exe_MIP_Kernel((float *)nin->data, (int)nrrdElementNumber(nin),inverseMatrix, sAxis, out);
195 :    
196 :     return 0;
197 :     }

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