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

SCM Repository

[diderot] Annotation of /benchmarks/programs/vr-lite-cam/bmark-teem.c
ViewVC logotype

Annotation of /benchmarks/programs/vr-lite-cam/bmark-teem.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : jhr 1533 /*! \file vr-lite-cam.c
2 :     *
3 :     * \author Nick Seltzer
4 :     */
5 :    
6 :     /*
7 :     * COPYRIGHT (c) 2011 The Diderot Project (http://diderot-language.cs.uchicago.edu)
8 :     * All rights reserved.
9 :     */
10 :    
11 :     #include <stdio.h>
12 :     #include <math.h>
13 :     #include <stdbool.h>
14 :     #include "teem/nrrd.h"
15 :     #include "teem/gage.h"
16 :    
17 :     #include <sys/time.h>
18 :    
19 :     static double GetTime ()
20 :     {
21 :     struct timeval tv;
22 :    
23 :     gettimeofday (&tv, 0);
24 :    
25 :     return (double)tv.tv_sec + 0.000001 * (double)tv.tv_usec;
26 :     }
27 :    
28 :     #define STATIC_INLINE static inline
29 :    
30 :     STATIC_INLINE double dot3(double vec1[3], double vec2[3])
31 :     {
32 :     return (vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] * vec2[2]);
33 :     }
34 :    
35 :     STATIC_INLINE double mag3(double vec[3])
36 :     {
37 :     return sqrt(dot3(vec, vec));
38 :     }
39 :    
40 :     STATIC_INLINE void norm3(double vec[3], double res[3])
41 :     {
42 :     double mag = mag3(vec);
43 :     res[0] = vec[0] / mag;
44 :     res[1] = vec[1] / mag;
45 :     res[2] = vec[2] / mag;
46 :     }
47 :    
48 :     STATIC_INLINE void add3(double vec1[3], double vec2[3], double res[3])
49 :     {
50 :     res[0] = vec1[0] + vec2[0];
51 :     res[1] = vec1[1] + vec2[1];
52 :     res[2] = vec1[2] + vec2[2];
53 :     }
54 :    
55 :     STATIC_INLINE void sub3(double vec1[3], double vec2[3], double res[3])
56 :     {
57 :     res[0] = vec1[0] - vec2[0];
58 :     res[1] = vec1[1] - vec2[1];
59 :     res[2] = vec1[2] - vec2[2];
60 :     }
61 :    
62 :     // Note res cannot be vec1 or vec2
63 :     STATIC_INLINE void cross(double vec1[3], double vec2[3], double res[3])
64 :     {
65 :     res[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
66 :     res[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
67 :     res[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
68 :     }
69 :    
70 :     STATIC_INLINE void scale_vec3(double scl, double vec[3], double res[3])
71 :     {
72 :     res[0] = scl * vec[0];
73 :     res[1] = scl * vec[1];
74 :     res[2] = scl * vec[2];
75 :     }
76 :    
77 :     STATIC_INLINE bool inside(double pos[3])
78 :     {
79 :     // XXX - Hack
80 :     if(pos[0] < -0.5 || pos[0] > 175.5) return false;
81 :     if(pos[1] < -0.5 || pos[1] > 186.5) return false;
82 :     if(pos[2] < -0.5 || pos[2] > 189.5) return false;
83 :     return true;
84 :     }
85 :    
86 :     STATIC_INLINE double lerp(double out_min, double out_max, double in_min, double in, double in_max)
87 :     {
88 :     return (out_min * (in_max - in) + out_max * (in - in_min)) / (in_max - in_min);
89 :     }
90 :    
91 :     STATIC_INLINE unsigned char ulerp(double in_min, double in, double in_max)
92 :     {
93 :     return (unsigned char)((0.0f * (in_max - in) + 255.0f * (in - in_min)) / (in_max - in_min));
94 :     }
95 :    
96 :     int main ()
97 :     {
98 :     char *infile = "../../data/vfrhand-nohip.nhdr";
99 :     char *outfile = "vr-lite-cam.pgm";
100 :     double Zero[3] = {0, 0, 0};
101 :     double camEye[3] = {127.331, -1322.05, 272.53};
102 :     double camAt[3] = {63.0, 82.6536, 98.0};
103 :     double camUp[3] = {0.9987, 0.0459166, -0.0221267};
104 :     double camNear = -78.0;
105 :     double camFar = 78.0;
106 :     double camFOV = 5.0;
107 :     int imgResU = 480;
108 :     int imgResV = 345;
109 :     double rayStep = 0.5;
110 :     double valOpacMin = 400.0; // 400.0 for skin, 1150.0 for bone
111 :     double valOpacMax = 700.0; // 700.0 for skin, 1450.0 for bone
112 :     double temp[3];
113 :     double temp2[3];
114 :     double temp3[3];
115 :     sub3(camAt, camEye, temp);
116 :     double camDist = mag3( temp );
117 :     double camVspNear = camDist + camNear;
118 :     double camVspFar = camDist + camFar;
119 :     double camN[3];
120 :     norm3(temp, camN);
121 :     double camU[3];
122 :     cross(camN, camUp, temp);
123 :     norm3(temp, camU);
124 :     double camV[3];
125 :     cross(camN, camU, camV);
126 :     double camVmax = tan(camFOV*3.1415926536/360.0)*camDist;
127 :     double camUmax = camVmax*(double)(imgResU)/(double)(imgResV);
128 :     double lightVspDir[3] = {0.9, -1.0, -2.5};
129 :     scale_vec3(lightVspDir[0], camU, temp);
130 :     scale_vec3(lightVspDir[1], camV, temp2);
131 :     add3(temp, temp2, temp);
132 :     scale_vec3(lightVspDir[2], camN, temp2);
133 :     add3(temp, temp2, temp);
134 :     double lightDir[3];
135 :     norm3(temp, lightDir);
136 :     double phongKa = 0.05;
137 :     double phongKd = 0.80;
138 :     double phongKs = 0.20;
139 :     double phongSp = 30.0;
140 :    
141 :     // Read in the data
142 :     char* err;
143 :     Nrrd *nin;
144 :     int status;
145 :     nin = nrrdNew();
146 : jhr 1575 if(nin == NULL) {
147 : jhr 1533 err = biffGetDone(NRRD);
148 :     fprintf(stderr, "Trouble allocating nrrd struct:\n%s", err);
149 :     free(err);
150 :     return -1;
151 :     }
152 :    
153 :    
154 :     status = nrrdLoad(nin, infile, NULL);
155 : jhr 1575 if (status) {
156 : jhr 1533 err = biffGetDone(NRRD);
157 :     fprintf(stderr, "Trouble reading \"%s\":\n%s", infile, err);
158 :     free(err);
159 :     nrrdNix(nin);
160 :     return -1;
161 :     }
162 :    
163 :     // Get data for mip
164 :     gageContext *ctx;
165 :     gagePerVolume *pvl;
166 :     const double *val;
167 :     const double *grad;
168 :     double kparm[NRRD_KERNEL_PARMS_NUM];
169 :    
170 :     ctx = gageContextNew();
171 : jhr 1575 if(ctx == NULL) {
172 : jhr 1533 err = biffGetDone(GAGE);
173 :     fprintf(stderr, "Trouble allocating gage context:\n%s", err);
174 :     free(err);
175 :     nrrdNuke(nin);
176 :     return -1;
177 :     }
178 :    
179 :     pvl = gagePerVolumeNew(ctx, nin, gageKindScl);
180 : jhr 1575 if(pvl == NULL) {
181 : jhr 1533 err = biffGetDone(GAGE);
182 :     fprintf(stderr, "Trouble allocating gage PerVolume:\n%s", err);
183 :     free(err);
184 :     gageContextNix(ctx);
185 :     nrrdNuke(nin);
186 :     return -1;
187 :     }
188 :    
189 :     status = gagePerVolumeAttach(ctx, pvl);
190 : jhr 1575 if(status) {
191 : jhr 1533 err = biffGetDone(GAGE);
192 :     fprintf(stderr, "Trouble attaching gage PerVolume:\n%s", err);
193 :     free(err);
194 :     gageContextNix(ctx);
195 :     nrrdNuke(nin);
196 :     return -1;
197 :     }
198 :    
199 :     // It's okay to use kparm unitialized here: gageKernelSet ignores it.
200 :     status = gageKernelSet(ctx, gageKernel00, nrrdKernelBSpline3, kparm);
201 : jhr 1575 if(status) {
202 : jhr 1533 err = biffGetDone(GAGE);
203 :     fprintf(stderr, "Trouble setting kernel:\n%s", err);
204 :     free(err);
205 :     gageContextNix(ctx);
206 :     nrrdNuke(nin);
207 :     return -1;
208 :     }
209 :    
210 :     // It's okay to use kparm unitialized here: gageKernelSet ignores it.
211 :     status = gageKernelSet(ctx, gageKernel11, nrrdKernelBSpline3D, kparm);
212 : jhr 1575 if(status) {
213 : jhr 1533 err = biffGetDone(GAGE);
214 :     fprintf(stderr, "Trouble setting kernel:\n%s", err);
215 :     free(err);
216 :     gageContextNix(ctx);
217 :     nrrdNuke(nin);
218 :     return -1;
219 :     }
220 :    
221 :     status = gageQueryItemOn(ctx, pvl, gageSclValue);
222 : jhr 1575 if(status) {
223 : jhr 1533 err = biffGetDone(GAGE);
224 :     fprintf(stderr, "Trouble with gageQueryItemOn:\n%s", err);
225 :     free(err);
226 :     gageContextNix(ctx);
227 :     nrrdNuke(nin);
228 :     return -1;
229 :     }
230 :    
231 :     status = gageQueryItemOn(ctx, pvl, gageSclGradVec);
232 : jhr 1575 if(status) {
233 : jhr 1533 err = biffGetDone(GAGE);
234 :     fprintf(stderr, "Trouble with gageQueryItemOn:\n%s", err);
235 :     free(err);
236 :     gageContextNix(ctx);
237 :     nrrdNuke(nin);
238 :     return -1;
239 :     }
240 :    
241 :     status = gageUpdate(ctx);
242 : jhr 1575 if(status) {
243 : jhr 1533 err = biffGetDone(GAGE);
244 :     fprintf(stderr, "Trouble with gageUpdate:\n%s", err);
245 :     free(err);
246 :     gageContextNix(ctx);
247 :     nrrdNuke(nin);
248 :     return -1;
249 :     }
250 :    
251 :     val = gageAnswerPointer(ctx, pvl, gageSclValue);
252 : jhr 1575 if(val == NULL) {
253 : jhr 1533 err = biffGetDone(GAGE);
254 :     fprintf(stderr, "Trouble getting answer pointer:\n%s", err);
255 :     free(err);
256 :     gageContextNix(ctx);
257 :     nrrdNuke(nin);
258 :     return -1;
259 :     }
260 :    
261 :     grad = gageAnswerPointer(ctx, pvl, gageSclGradVec);
262 : jhr 1575 if(grad == NULL) {
263 : jhr 1533 err = biffGetDone(GAGE);
264 :     fprintf(stderr, "Trouble getting answer pointer:\n%s", err);
265 :     free(err);
266 :     gageContextNix(ctx);
267 :     nrrdNuke(nin);
268 :     return -1;
269 :     }
270 :    
271 :     int ui, vi;
272 :     double rayU, rayV, rayN;
273 :     double rayVec[3];
274 :     double rayDir[3];
275 :     double rayPos[3];
276 :     double transp, gray;
277 :     double output[4];
278 :     double norm[3];
279 :     double alpha;
280 :     double ld;
281 :     double hd;
282 :     double mat;
283 :     bool set = false;
284 :     double max_gray = 0;
285 :     double min_gray = 0;
286 :     double *out_data = malloc(sizeof(double) * 4 * imgResU * imgResV);
287 : jhr 1575 if(out_data == NULL) {
288 : jhr 1533 fprintf(stderr, "Trouble with malloc.\n");
289 :     gageContextNix(ctx);
290 :     nrrdNuke(nin);
291 :     return -1;
292 :     }
293 :    
294 : jhr 1575 double t0 = GetTime(); // start timing
295 :     for(vi = 0; vi < imgResV; vi++) {
296 : jhr 1533 rayV = lerp(-camVmax, camVmax, -0.5, (double)(vi), (double)(imgResV) - 0.5);
297 :     scale_vec3(rayV, camV, temp);
298 : jhr 1575 for(ui = 0; ui < imgResU; ui++) {
299 : jhr 1533 rayU = lerp(-camUmax, camUmax, -0.5, (double)(ui), (double)(imgResU) - 0.5);
300 :     scale_vec3(rayU, camU, temp2);
301 :     add3(temp, temp2, temp2);
302 :     scale_vec3(camDist, camN, temp3);
303 :     add3(temp2, temp3, temp3);
304 :     scale_vec3(1 / camDist, temp3, rayVec);
305 :     norm3(rayVec, rayDir);
306 :     transp = 1;
307 :     gray = 0;
308 :     output[0] = 0;
309 :     output[1] = 0;
310 :     output[2] = 0;
311 :     output[3] = 0;
312 : jhr 1575 for(rayN = camVspNear; rayN <= camVspFar; rayN += rayStep) {
313 : jhr 1533 scale_vec3(rayN, rayDir, temp2);
314 :     add3(temp2, camEye, rayPos);
315 : jhr 1575 if(inside(rayPos)) {
316 : jhr 1533 status = gageProbe(ctx, rayPos[0], rayPos[1], rayPos[2]);
317 : jhr 1575 if(status) {
318 : jhr 1533 err = ctx->errStr;
319 :     fprintf(stderr, "Trouble with gageProbe:\n%s", err);
320 :     free(err);
321 :     free(out_data);
322 :     gageContextNix(ctx);
323 :     nrrdNuke(nin);
324 :     return -1;
325 :     }
326 : jhr 1575 if(*val > valOpacMin) {
327 : jhr 1533 temp2[0] = *(grad + 0);
328 :     temp2[1] = *(grad + 1);
329 :     temp2[2] = *(grad + 2);
330 :     sub3(Zero, temp2, temp2);
331 :     norm3(temp2, norm);
332 :     alpha = lerp(0.0, 1.0, valOpacMin, *val, valOpacMax);
333 :     alpha = (alpha < 1.0) ? alpha : 1;
334 :     ld = dot3(norm, lightDir);
335 :     ld = (ld < 0) ? 0 : ld;
336 :     sub3(camEye, rayPos, temp2);
337 :     norm3(temp2, temp2);
338 :     add3(temp2, lightDir, temp2);
339 :     norm3(temp2, temp2);
340 :     hd = dot3(norm, temp2);
341 :     hd = (hd < 0) ? 0 : hd;
342 :     mat = phongKa +
343 :     phongKd * ((ld > 0) ? ld : 0) +
344 :     phongKs * ((hd > 0) ? pow(hd, phongSp) : 0);
345 :     gray = gray + transp * alpha * mat;
346 :     transp *= 1 - alpha;
347 :     }
348 :     }
349 : jhr 1575 if(transp < 0.01) {
350 : jhr 1533 transp = 0;
351 :     break;
352 :     }
353 :     }
354 :     *(out_data + ui * imgResV * 4 + vi * 4 + 0) = gray;
355 :     *(out_data + ui * imgResV * 4 + vi * 4 + 1) = gray;
356 :     *(out_data + ui * imgResV * 4 + vi * 4 + 2) = gray;
357 :     *(out_data + ui * imgResV * 4 + vi * 4 + 3) = 1.0f - transp;
358 : jhr 1575 if(!set) {
359 : jhr 1533 max_gray = gray;
360 :     min_gray = gray;
361 :     set = true;
362 :     }
363 :     max_gray = (gray > max_gray) ? gray : max_gray;
364 :     min_gray = (gray < min_gray) ? gray : min_gray;
365 :     }
366 :     }
367 :    
368 :     unsigned char *uc_out_data = malloc(sizeof(unsigned char) * imgResU * imgResV);
369 : jhr 1575 for(vi = 0; vi < imgResV; vi++) {
370 :     for(ui = 0; ui < imgResU; ui++) {
371 : jhr 1533 *(uc_out_data + vi * imgResU + ui) = ulerp(min_gray, *(out_data + ui * imgResV * 4 + vi * 4 + 0), max_gray);
372 :     }
373 :     }
374 : jhr 1575
375 : jhr 1533 double totalTime = GetTime() - t0;
376 :     printf("usr=%f\n", totalTime);
377 :    
378 :     // Write out data
379 :     Nrrd *nout;
380 :     nout = nrrdNew();
381 : jhr 1575 if(nout == NULL) {
382 : jhr 1533 err = biffGetDone(NRRD);
383 :     fprintf(stderr, "Trouble allocating nrrd struct:\n%s", err);
384 :     free(err);
385 :     free(out_data);
386 :     free(uc_out_data);
387 :     gageContextNix(ctx);
388 :     nrrdNuke(nin);
389 :     return -1;
390 :     }
391 :    
392 :     status = nrrdWrap_va(nout, uc_out_data, nrrdTypeUChar, 2, imgResU, imgResV);
393 : jhr 1575 if(status) {
394 : jhr 1533 err = biffGetDone(NRRD);
395 :     fprintf(stderr, "Trouble wrapping nrrd struct:\n%s", err);
396 :     free(err);
397 :     nrrdNix(nout);
398 :     free(out_data);
399 :     gageContextNix(ctx);
400 :     nrrdNuke(nin);
401 :     return -1;
402 :     }
403 :    
404 :     status = nrrdSave(outfile, nout, NULL);
405 : jhr 1575 if(status) {
406 : jhr 1533 err = biffGetDone(NRRD);
407 :     fprintf(stderr, "Trouble saving nrrd struct:\n%s", err);
408 :     free(err);
409 :     nrrdNuke(nout);
410 :     free(out_data);
411 :     gageContextNix(ctx);
412 :     nrrdNuke(nin);
413 :     return -1;
414 :     }
415 :    
416 :    
417 :     free(out_data);
418 :     nrrdNuke(nout);
419 :     gageContextNix(ctx);
420 :     nrrdNuke(nin);
421 :     return 0;
422 :    
423 :     }

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