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

SCM Repository

[diderot] Annotation of /trunk/test/teem/vr-lite-cam.c
ViewVC logotype

Annotation of /trunk/test/teem/vr-lite-cam.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : jhr 1299 /*! \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 : nseltzer 1294 #include <stdio.h>
12 :     #include <math.h>
13 :     #include <stdbool.h>
14 :     #include "teem/nrrd.h"
15 :     #include "teem/gage.h"
16 :    
17 : jhr 1299 #include <sys/time.h>
18 : nseltzer 1294
19 : jhr 1299 static double GetTime ()
20 : nseltzer 1294 {
21 : jhr 1299 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 : nseltzer 1294 return (vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] * vec2[2]);
33 :     }
34 :    
35 : jhr 1299 STATIC_INLINE double mag3(double vec[3])
36 : nseltzer 1294 {
37 :     return sqrt(dot3(vec, vec));
38 :     }
39 :    
40 : jhr 1299 STATIC_INLINE void norm3(double vec[3], double res[3])
41 : nseltzer 1294 {
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 : jhr 1299 STATIC_INLINE void add3(double vec1[3], double vec2[3], double res[3])
49 : nseltzer 1294 {
50 :     res[0] = vec1[0] + vec2[0];
51 :     res[1] = vec1[1] + vec2[1];
52 :     res[2] = vec1[2] + vec2[2];
53 :     }
54 :    
55 : jhr 1299 STATIC_INLINE void sub3(double vec1[3], double vec2[3], double res[3])
56 : nseltzer 1294 {
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 : jhr 1299 STATIC_INLINE void cross(double vec1[3], double vec2[3], double res[3])
64 : nseltzer 1294 {
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 : jhr 1299 STATIC_INLINE void scale_vec3(double scl, double vec[3], double res[3])
71 : nseltzer 1294 {
72 :     res[0] = scl * vec[0];
73 :     res[1] = scl * vec[1];
74 :     res[2] = scl * vec[2];
75 :     }
76 :    
77 : jhr 1299 STATIC_INLINE bool inside(double pos[3])
78 : nseltzer 1294 {
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 : jhr 1299 STATIC_INLINE double lerp(double out_min, double out_max, double in_min, double in, double in_max)
87 : nseltzer 1294 {
88 :     return (out_min * (in_max - in) + out_max * (in - in_min)) / (in_max - in_min);
89 :     }
90 :    
91 : jhr 1299 STATIC_INLINE unsigned char ulerp(double in_min, double in, double in_max)
92 : nseltzer 1294 {
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 : jhr 1299 char *outfile = "vr-lite-cam.pgm";
100 : nseltzer 1294 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 :     if(nin == NULL)
147 :     {
148 :     err = biffGetDone(NRRD);
149 :     fprintf(stderr, "Trouble allocating nrrd struct:\n%s", err);
150 :     free(err);
151 :     return -1;
152 :     }
153 :    
154 :    
155 :     status = nrrdLoad(nin, infile, NULL);
156 :     if (status)
157 :     {
158 :     err = biffGetDone(NRRD);
159 :     fprintf(stderr, "Trouble reading \"%s\":\n%s", infile, err);
160 :     free(err);
161 :     nrrdNix(nin);
162 :     return -1;
163 :     }
164 :    
165 :     // Get data for mip
166 :     gageContext *ctx;
167 :     gagePerVolume *pvl;
168 :     const double *val;
169 :     const double *grad;
170 :     double kparm[NRRD_KERNEL_PARMS_NUM];
171 :    
172 :     ctx = gageContextNew();
173 :     if(ctx == NULL)
174 :     {
175 :     err = biffGetDone(GAGE);
176 :     fprintf(stderr, "Trouble allocating gage context:\n%s", err);
177 :     free(err);
178 :     nrrdNuke(nin);
179 :     return -1;
180 :     }
181 :    
182 :     pvl = gagePerVolumeNew(ctx, nin, gageKindScl);
183 :     if(pvl == NULL)
184 :     {
185 :     err = biffGetDone(GAGE);
186 :     fprintf(stderr, "Trouble allocating gage PerVolume:\n%s", err);
187 :     free(err);
188 :     gageContextNix(ctx);
189 :     nrrdNuke(nin);
190 :     return -1;
191 :     }
192 :    
193 :     status = gagePerVolumeAttach(ctx, pvl);
194 :     if(status)
195 :     {
196 :     err = biffGetDone(GAGE);
197 :     fprintf(stderr, "Trouble attaching gage PerVolume:\n%s", err);
198 :     free(err);
199 :     gageContextNix(ctx);
200 :     nrrdNuke(nin);
201 :     return -1;
202 :     }
203 :    
204 :     // It's okay to use kparm unitialized here: gageKernelSet ignores it.
205 :     status = gageKernelSet(ctx, gageKernel00, nrrdKernelBSpline3, kparm);
206 :     if(status)
207 :     {
208 :     err = biffGetDone(GAGE);
209 :     fprintf(stderr, "Trouble setting kernel:\n%s", err);
210 :     free(err);
211 :     gageContextNix(ctx);
212 :     nrrdNuke(nin);
213 :     return -1;
214 :     }
215 :    
216 :     // It's okay to use kparm unitialized here: gageKernelSet ignores it.
217 :     status = gageKernelSet(ctx, gageKernel11, nrrdKernelBSpline3D, kparm);
218 :     if(status)
219 :     {
220 :     err = biffGetDone(GAGE);
221 :     fprintf(stderr, "Trouble setting kernel:\n%s", err);
222 :     free(err);
223 :     gageContextNix(ctx);
224 :     nrrdNuke(nin);
225 :     return -1;
226 :     }
227 :    
228 :     status = gageQueryItemOn(ctx, pvl, gageSclValue);
229 :     if(status)
230 :     {
231 :     err = biffGetDone(GAGE);
232 :     fprintf(stderr, "Trouble with gageQueryItemOn:\n%s", err);
233 :     free(err);
234 :     gageContextNix(ctx);
235 :     nrrdNuke(nin);
236 :     return -1;
237 :     }
238 :    
239 :     status = gageQueryItemOn(ctx, pvl, gageSclGradVec);
240 :     if(status)
241 :     {
242 :     err = biffGetDone(GAGE);
243 :     fprintf(stderr, "Trouble with gageQueryItemOn:\n%s", err);
244 :     free(err);
245 :     gageContextNix(ctx);
246 :     nrrdNuke(nin);
247 :     return -1;
248 :     }
249 :    
250 :     status = gageUpdate(ctx);
251 :     if(status)
252 :     {
253 :     err = biffGetDone(GAGE);
254 :     fprintf(stderr, "Trouble with gageUpdate:\n%s", err);
255 :     free(err);
256 :     gageContextNix(ctx);
257 :     nrrdNuke(nin);
258 :     return -1;
259 :     }
260 :    
261 :     val = gageAnswerPointer(ctx, pvl, gageSclValue);
262 :     if(val == NULL)
263 :     {
264 :     err = biffGetDone(GAGE);
265 :     fprintf(stderr, "Trouble getting answer pointer:\n%s", err);
266 :     free(err);
267 :     gageContextNix(ctx);
268 :     nrrdNuke(nin);
269 :     return -1;
270 :     }
271 :    
272 :     grad = gageAnswerPointer(ctx, pvl, gageSclGradVec);
273 :     if(grad == NULL)
274 :     {
275 :     err = biffGetDone(GAGE);
276 :     fprintf(stderr, "Trouble getting answer pointer:\n%s", err);
277 :     free(err);
278 :     gageContextNix(ctx);
279 :     nrrdNuke(nin);
280 :     return -1;
281 :     }
282 :    
283 :    
284 :     int ui, vi;
285 :     double rayU, rayV, rayN;
286 :     double rayVec[3];
287 :     double rayDir[3];
288 :     double rayPos[3];
289 :     double transp, gray;
290 :     double output[4];
291 :     double norm[3];
292 :     double alpha;
293 :     double ld;
294 :     double hd;
295 :     double mat;
296 :     bool set = false;
297 :     double max_gray = 0;
298 :     double min_gray = 0;
299 :     double *out_data = malloc(sizeof(double) * 4 * imgResU * imgResV);
300 :     if(out_data == NULL)
301 :     {
302 :     fprintf(stderr, "Trouble with malloc.\n");
303 :     gageContextNix(ctx);
304 :     nrrdNuke(nin);
305 :     return -1;
306 :     }
307 :    
308 : jhr 1299 double t0 = GetTime();
309 : nseltzer 1294 for(vi = 0; vi < imgResV; vi++)
310 :     {
311 :     rayV = lerp(-camVmax, camVmax, -0.5, (double)(vi), (double)(imgResV) - 0.5);
312 :     scale_vec3(rayV, camV, temp);
313 :     for(ui = 0; ui < imgResU; ui++)
314 :     {
315 :     rayU = lerp(-camUmax, camUmax, -0.5, (double)(ui), (double)(imgResU) - 0.5);
316 :     scale_vec3(rayU, camU, temp2);
317 :     add3(temp, temp2, temp2);
318 :     scale_vec3(camDist, camN, temp3);
319 :     add3(temp2, temp3, temp3);
320 :     scale_vec3(1 / camDist, temp3, rayVec);
321 :     norm3(rayVec, rayDir);
322 :     transp = 1;
323 :     gray = 0;
324 :     output[0] = 0;
325 :     output[1] = 0;
326 :     output[2] = 0;
327 :     output[3] = 0;
328 :     for(rayN = camVspNear; rayN <= camVspFar; rayN += rayStep)
329 :     {
330 :     scale_vec3(rayN, rayDir, temp2);
331 :     add3(temp2, camEye, rayPos);
332 :     if(inside(rayPos))
333 :     {
334 :     status = gageProbe(ctx, rayPos[0], rayPos[1], rayPos[2]);
335 :     if(status)
336 :     {
337 :     err = ctx->errStr;
338 :     fprintf(stderr, "Trouble with gageProbe:\n%s", err);
339 :     free(err);
340 :     free(out_data);
341 :     gageContextNix(ctx);
342 :     nrrdNuke(nin);
343 :     return -1;
344 :     }
345 :     if(*val > valOpacMin)
346 :     {
347 :     temp2[0] = *(grad + 0);
348 :     temp2[1] = *(grad + 1);
349 :     temp2[2] = *(grad + 2);
350 :     sub3(Zero, temp2, temp2);
351 :     norm3(temp2, norm);
352 :     alpha = lerp(0.0, 1.0, valOpacMin, *val, valOpacMax);
353 :     alpha = (alpha < 1.0) ? alpha : 1;
354 :     ld = dot3(norm, lightDir);
355 :     ld = (ld < 0) ? 0 : ld;
356 :     sub3(camEye, rayPos, temp2);
357 :     norm3(temp2, temp2);
358 :     add3(temp2, lightDir, temp2);
359 :     norm3(temp2, temp2);
360 :     hd = dot3(norm, temp2);
361 :     hd = (hd < 0) ? 0 : hd;
362 :     mat = phongKa +
363 :     phongKd * ((ld > 0) ? ld : 0) +
364 :     phongKs * ((hd > 0) ? pow(hd, phongSp) : 0);
365 :     gray = gray + transp * alpha * mat;
366 :     transp *= 1 - alpha;
367 :     }
368 :     }
369 :     if(transp < 0.01)
370 :     {
371 :     transp = 0;
372 :     break;
373 :     }
374 :     }
375 :     *(out_data + ui * imgResV * 4 + vi * 4 + 0) = gray;
376 :     *(out_data + ui * imgResV * 4 + vi * 4 + 1) = gray;
377 :     *(out_data + ui * imgResV * 4 + vi * 4 + 2) = gray;
378 :     *(out_data + ui * imgResV * 4 + vi * 4 + 3) = 1.0f - transp;
379 :     if(!set)
380 :     {
381 :     max_gray = gray;
382 :     min_gray = gray;
383 :     set = true;
384 :     }
385 :     max_gray = (gray > max_gray) ? gray : max_gray;
386 :     min_gray = (gray < min_gray) ? gray : min_gray;
387 :     }
388 :     }
389 :    
390 :     unsigned char *uc_out_data = malloc(sizeof(unsigned char) * imgResU * imgResV);
391 :     for(vi = 0; vi < imgResV; vi++)
392 :     {
393 :     for(ui = 0; ui < imgResU; ui++)
394 :     {
395 :     *(uc_out_data + vi * imgResU + ui) = ulerp(min_gray, *(out_data + ui * imgResV * 4 + vi * 4 + 0), max_gray);
396 :     }
397 :     }
398 : jhr 1299 double totalTime = GetTime() - t0;
399 : nseltzer 1294
400 : jhr 1299 printf("usr=%f\n", totalTime);
401 :    
402 : nseltzer 1294 // Write out data
403 :     Nrrd *nout;
404 :     nout = nrrdNew();
405 :     if(nout == NULL)
406 :     {
407 :     err = biffGetDone(NRRD);
408 :     fprintf(stderr, "Trouble allocating nrrd struct:\n%s", err);
409 :     free(err);
410 :     free(out_data);
411 :     free(uc_out_data);
412 :     gageContextNix(ctx);
413 :     nrrdNuke(nin);
414 :     return -1;
415 :     }
416 :    
417 :     status = nrrdWrap_va(nout, uc_out_data, nrrdTypeUChar, 2, imgResU, imgResV);
418 :     if(status)
419 :     {
420 :     err = biffGetDone(NRRD);
421 :     fprintf(stderr, "Trouble wrapping nrrd struct:\n%s", err);
422 :     free(err);
423 :     nrrdNix(nout);
424 :     free(out_data);
425 :     gageContextNix(ctx);
426 :     nrrdNuke(nin);
427 :     return -1;
428 :     }
429 :    
430 :     status = nrrdSave(outfile, nout, NULL);
431 :     if(status)
432 :     {
433 :     err = biffGetDone(NRRD);
434 :     fprintf(stderr, "Trouble saving nrrd struct:\n%s", err);
435 :     free(err);
436 :     nrrdNuke(nout);
437 :     free(out_data);
438 :     gageContextNix(ctx);
439 :     nrrdNuke(nin);
440 :     return -1;
441 :     }
442 :    
443 :    
444 :     free(out_data);
445 :     nrrdNuke(nout);
446 :     gageContextNix(ctx);
447 :     nrrdNuke(nin);
448 :     return 0;
449 :    
450 :     }

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