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 1602 - (view) (download) (as text)

1 : jhr 1602 /*! \file bmark-teem.c
2 : jhr 1533 *
3 : glk 1600 * \author Nick Seltzer & Gordon Kindlmann
4 : jhr 1533 */
5 :    
6 : jhr 1602 /*
7 :     * COPYRIGHT (c) 2011 The Diderot Project (http://diderot-language.cs.uchicago.edu)
8 :     * All rights reserved.
9 :     */
10 :    
11 : jhr 1533 #include <stdio.h>
12 :     #include <math.h>
13 :     #include <stdbool.h>
14 :     #include "teem/nrrd.h"
15 :     #include "teem/gage.h"
16 :    
17 : jhr 1586 #include "teem-defs.h"
18 : jhr 1533
19 :     STATIC_INLINE double dot3(double vec1[3], double vec2[3])
20 :     {
21 :     return (vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] * vec2[2]);
22 :     }
23 :    
24 :     STATIC_INLINE double mag3(double vec[3])
25 :     {
26 :     return sqrt(dot3(vec, vec));
27 :     }
28 :    
29 :     STATIC_INLINE void norm3(double vec[3], double res[3])
30 :     {
31 :     double mag = mag3(vec);
32 :     res[0] = vec[0] / mag;
33 :     res[1] = vec[1] / mag;
34 :     res[2] = vec[2] / mag;
35 :     }
36 :    
37 :     STATIC_INLINE void add3(double vec1[3], double vec2[3], double res[3])
38 :     {
39 :     res[0] = vec1[0] + vec2[0];
40 :     res[1] = vec1[1] + vec2[1];
41 :     res[2] = vec1[2] + vec2[2];
42 :     }
43 :    
44 :     STATIC_INLINE void sub3(double vec1[3], double vec2[3], double res[3])
45 :     {
46 :     res[0] = vec1[0] - vec2[0];
47 :     res[1] = vec1[1] - vec2[1];
48 :     res[2] = vec1[2] - vec2[2];
49 :     }
50 :    
51 :     // Note res cannot be vec1 or vec2
52 :     STATIC_INLINE void cross(double vec1[3], double vec2[3], double res[3])
53 :     {
54 :     res[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
55 :     res[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
56 :     res[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
57 :     }
58 :    
59 :     STATIC_INLINE void scale_vec3(double scl, double vec[3], double res[3])
60 :     {
61 :     res[0] = scl * vec[0];
62 :     res[1] = scl * vec[1];
63 :     res[2] = scl * vec[2];
64 :     }
65 :    
66 :     STATIC_INLINE bool inside(double pos[3])
67 :     {
68 :     // XXX - Hack
69 :     if(pos[0] < -0.5 || pos[0] > 175.5) return false;
70 :     if(pos[1] < -0.5 || pos[1] > 186.5) return false;
71 :     if(pos[2] < -0.5 || pos[2] > 189.5) return false;
72 :     return true;
73 :     }
74 :    
75 :     STATIC_INLINE double lerp(double out_min, double out_max, double in_min, double in, double in_max)
76 :     {
77 :     return (out_min * (in_max - in) + out_max * (in - in_min)) / (in_max - in_min);
78 :     }
79 :    
80 :     STATIC_INLINE unsigned char ulerp(double in_min, double in, double in_max)
81 :     {
82 :     return (unsigned char)((0.0f * (in_max - in) + 255.0f * (in - in_min)) / (in_max - in_min));
83 :     }
84 :    
85 : glk 1600 int main (int argc, const char *argv[])
86 : jhr 1533 {
87 : glk 1600 const char *me = argv[0];
88 : jhr 1533 char *infile = "../../data/vfrhand-nohip.nhdr";
89 :     char *outfile = "vr-lite-cam.pgm";
90 :     double camEye[3] = {127.331, -1322.05, 272.53};
91 :     double camAt[3] = {63.0, 82.6536, 98.0};
92 :     double camUp[3] = {0.9987, 0.0459166, -0.0221267};
93 :     double camNear = -78.0;
94 :     double camFar = 78.0;
95 :     double camFOV = 5.0;
96 :     int imgResU = 480;
97 :     int imgResV = 345;
98 :     double rayStep = 0.5;
99 :     double valOpacMin = 400.0; // 400.0 for skin, 1150.0 for bone
100 :     double valOpacMax = 700.0; // 700.0 for skin, 1450.0 for bone
101 :     double temp[3];
102 :     double temp2[3];
103 :     double temp3[3];
104 :     sub3(camAt, camEye, temp);
105 :     double camDist = mag3( temp );
106 :     double camVspNear = camDist + camNear;
107 :     double camVspFar = camDist + camFar;
108 :     double camN[3];
109 :     norm3(temp, camN);
110 :     double camU[3];
111 :     cross(camN, camUp, temp);
112 :     norm3(temp, camU);
113 :     double camV[3];
114 :     cross(camN, camU, camV);
115 : glk 1600 double camVmax = tan(camFOV*AIR_PI/360.0)*camDist;
116 : jhr 1533 double camUmax = camVmax*(double)(imgResU)/(double)(imgResV);
117 :     double lightVspDir[3] = {0.9, -1.0, -2.5};
118 :     scale_vec3(lightVspDir[0], camU, temp);
119 :     scale_vec3(lightVspDir[1], camV, temp2);
120 :     add3(temp, temp2, temp);
121 :     scale_vec3(lightVspDir[2], camN, temp2);
122 :     add3(temp, temp2, temp);
123 :     double lightDir[3];
124 :     norm3(temp, lightDir);
125 :     double phongKa = 0.05;
126 :     double phongKd = 0.80;
127 :     double phongKs = 0.20;
128 : glk 1600 double phongSp = 50.0;
129 : jhr 1533
130 :     // Read in the data
131 :     char* err;
132 :     Nrrd *nin;
133 :     int status;
134 : glk 1600 airArray *mop;
135 :     mop = airMopNew();
136 :    
137 : jhr 1533 nin = nrrdNew();
138 : glk 1600 airMopAdd(mop, nin, (airMopper)nrrdNuke, airMopAlways);
139 : jhr 1533
140 : glk 1600 if (nrrdLoad(nin, infile, NULL)) {
141 :     airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
142 : jhr 1533 fprintf(stderr, "Trouble reading \"%s\":\n%s", infile, err);
143 : glk 1600 airMopError(mop);
144 : jhr 1533 return -1;
145 :     }
146 :    
147 :     // Get data for mip
148 :     gageContext *ctx;
149 :     gagePerVolume *pvl;
150 :     const double *val;
151 :     const double *grad;
152 : glk 1600 double kparm[NRRD_KERNEL_PARMS_NUM] = {0.0};
153 :     int E;
154 : jhr 1533
155 :     ctx = gageContextNew();
156 : glk 1600 airMopAdd(mop, ctx, (airMopper)gageContextNix, airMopAlways);
157 : jhr 1533 pvl = gagePerVolumeNew(ctx, nin, gageKindScl);
158 : glk 1600 if (!( ctx && pvl )) {
159 :     fprintf(stderr, "%s: couldn't allocate gage ctx or pvl", me);
160 :     airMopError(mop);
161 : jhr 1533 return -1;
162 :     }
163 : glk 1600 if (gagePerVolumeAttach(ctx, pvl)
164 :     || gageKernelSet(ctx, gageKernel00, nrrdKernelBSpline3, kparm)
165 :     || gageKernelSet(ctx, gageKernel11, nrrdKernelBSpline3D, kparm)
166 :     || gageQueryItemOn(ctx, pvl, gageSclValue)
167 :     || gageQueryItemOn(ctx, pvl, gageSclGradVec)
168 :     || gageUpdate(ctx)) {
169 :     airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways);
170 :     fprintf(stderr, "Trouble setting up gage context:\n%s", err);
171 :     airMopError(mop);
172 : jhr 1533 return -1;
173 :     }
174 :    
175 : glk 1600 if (!( (val = gageAnswerPointer(ctx, pvl, gageSclValue))
176 :     && (grad = gageAnswerPointer(ctx, pvl, gageSclGradVec)) )) {
177 :     fprintf(stderr, "Trouble getting answer pointer\n");
178 :     airMopError(mop);
179 : jhr 1533 return -1;
180 :     }
181 :    
182 :     int ui, vi;
183 :     double rayU, rayV, rayN;
184 :     double rayVec[3];
185 :     double rayDir[3];
186 :     double rayPos[3];
187 :     double transp, gray;
188 :     double output[4];
189 :     double norm[3];
190 :     double alpha;
191 :     double ld;
192 :     double hd;
193 :     double mat;
194 :     bool set = false;
195 :     double max_gray = 0;
196 :     double min_gray = 0;
197 : glk 1600
198 :     // allocate image and output
199 :     Nrrd *nimg = nrrdNew();
200 :     airMopAdd(mop, nimg, (airMopper)nrrdNuke, airMopAlways);
201 :     Nrrd *nout = nrrdNew();
202 :     airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);
203 :     if (nrrdAlloc_va(nimg, nrrdTypeDouble, 3,
204 :     AIR_CAST(size_t, 4),
205 :     AIR_CAST(size_t, imgResU),
206 :     AIR_CAST(size_t, imgResV))
207 :     || nrrdAlloc_va(nout, nrrdTypeUChar, 2,
208 :     AIR_CAST(size_t, imgResU),
209 :     AIR_CAST(size_t, imgResV))) {
210 :     airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
211 :     fprintf(stderr, "Trouble allocating image and output:\n%s", err);
212 :     airMopError(mop);
213 : jhr 1533 return -1;
214 :     }
215 : glk 1600 double *out_data = AIR_CAST(double *, nimg->data);
216 :     unsigned char *uc_out_data = AIR_CAST(unsigned char *, nout->data);
217 : jhr 1533
218 : glk 1600 double t0 = airTime(); // start timing
219 : jhr 1575 for(vi = 0; vi < imgResV; vi++) {
220 : jhr 1533 rayV = lerp(-camVmax, camVmax, -0.5, (double)(vi), (double)(imgResV) - 0.5);
221 :     scale_vec3(rayV, camV, temp);
222 : jhr 1575 for(ui = 0; ui < imgResU; ui++) {
223 : jhr 1533 rayU = lerp(-camUmax, camUmax, -0.5, (double)(ui), (double)(imgResU) - 0.5);
224 :     scale_vec3(rayU, camU, temp2);
225 :     add3(temp, temp2, temp2);
226 :     scale_vec3(camDist, camN, temp3);
227 :     add3(temp2, temp3, temp3);
228 :     scale_vec3(1 / camDist, temp3, rayVec);
229 :     norm3(rayVec, rayDir);
230 :     transp = 1;
231 :     gray = 0;
232 :     output[0] = 0;
233 :     output[1] = 0;
234 :     output[2] = 0;
235 :     output[3] = 0;
236 : glk 1600 for (rayN = camVspNear; rayN <= camVspFar; rayN += rayStep) {
237 : jhr 1533 scale_vec3(rayN, rayDir, temp2);
238 :     add3(temp2, camEye, rayPos);
239 : glk 1600 if (inside(rayPos)) {
240 :     if (gageProbe(ctx, rayPos[0], rayPos[1], rayPos[2])) {
241 : glk 1587 fprintf(stderr, "Trouble with gageProbe:\n%s", ctx->errStr);
242 : glk 1600 airMopError(mop); return -1;
243 : jhr 1533 }
244 : glk 1600 if (*val > valOpacMin) {
245 :     temp2[0] = -*(grad + 0);
246 :     temp2[1] = -*(grad + 1);
247 :     temp2[2] = -*(grad + 2);
248 : jhr 1533 norm3(temp2, norm);
249 :     alpha = lerp(0.0, 1.0, valOpacMin, *val, valOpacMax);
250 :     alpha = (alpha < 1.0) ? alpha : 1;
251 :     ld = dot3(norm, lightDir);
252 :     ld = (ld < 0) ? 0 : ld;
253 :     sub3(camEye, rayPos, temp2);
254 :     norm3(temp2, temp2);
255 :     add3(temp2, lightDir, temp2);
256 :     norm3(temp2, temp2);
257 :     hd = dot3(norm, temp2);
258 :     hd = (hd < 0) ? 0 : hd;
259 :     mat = phongKa +
260 :     phongKd * ((ld > 0) ? ld : 0) +
261 :     phongKs * ((hd > 0) ? pow(hd, phongSp) : 0);
262 :     gray = gray + transp * alpha * mat;
263 :     transp *= 1 - alpha;
264 :     }
265 :     }
266 : jhr 1575 if(transp < 0.01) {
267 : jhr 1533 transp = 0;
268 :     break;
269 :     }
270 :     }
271 :     *(out_data + ui * imgResV * 4 + vi * 4 + 0) = gray;
272 :     *(out_data + ui * imgResV * 4 + vi * 4 + 1) = gray;
273 :     *(out_data + ui * imgResV * 4 + vi * 4 + 2) = gray;
274 :     *(out_data + ui * imgResV * 4 + vi * 4 + 3) = 1.0f - transp;
275 : jhr 1575 if(!set) {
276 : jhr 1533 max_gray = gray;
277 :     min_gray = gray;
278 :     set = true;
279 :     }
280 :     max_gray = (gray > max_gray) ? gray : max_gray;
281 :     min_gray = (gray < min_gray) ? gray : min_gray;
282 :     }
283 :     }
284 :    
285 : jhr 1575 for(vi = 0; vi < imgResV; vi++) {
286 :     for(ui = 0; ui < imgResU; ui++) {
287 : jhr 1533 *(uc_out_data + vi * imgResU + ui) = ulerp(min_gray, *(out_data + ui * imgResV * 4 + vi * 4 + 0), max_gray);
288 :     }
289 :     }
290 : jhr 1575
291 : glk 1600 double totalTime = airTime() - t0;
292 : jhr 1533 printf("usr=%f\n", totalTime);
293 :    
294 : glk 1600 if (nrrdSave(outfile, nout, NULL)) {
295 :     airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
296 :     fprintf(stderr, "Trouble writing output:\n%s", err);
297 :     airMopError(mop);
298 : jhr 1533 return -1;
299 :     }
300 :    
301 : glk 1600 airMopOkay(mop);
302 : jhr 1533 return 0;
303 :    
304 :     }

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