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

SCM Repository

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

Annotation of /benchmarks/programs/illust-vr/bmark-teem.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : glk 1609 /*! \file bmark-teem.c
2 :     *
3 :     * Gordon Kindlmann
4 :     */
5 :    
6 :     #if 1
7 :     # define REAL float
8 :     # define ATAN2 atan2f
9 :     # define SQRT sqrtf
10 :     #else
11 :     # define REAL double
12 :     # define ATAN2 atan2
13 :     # define SQRT sqrt
14 :     #endif
15 :    
16 :     #define ELL_3M_FROBSQRD(m) \
17 :     (ELL_3V_DOT((m)+0, (m)+0) + \
18 :     ELL_3V_DOT((m)+3, (m)+3) + \
19 :     ELL_3V_DOT((m)+6, (m)+6))
20 :    
21 :    
22 :     #include <stdio.h>
23 :     #include <math.h>
24 :     #include <stdbool.h>
25 :     #include "teem/nrrd.h"
26 :     #include "teem/gage.h"
27 :    
28 :     #include "teem-defs.h"
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(REAL pos[3])
78 :     {
79 :     // XXX - Hack specific to ../../data/vfrhand-nohip.nhdr"
80 :     // but for some reason still not matching what Diderot says
81 :     if(pos[0] < 2*0.9375 || pos[0] > (176-3)*0.9375) return false;
82 :     if(pos[1] < 2*0.9375 || pos[1] > (187-3)*0.9375) return false;
83 :     if(pos[2] < 2.0 || pos[2] > (190-3)) return false;
84 :     return true;
85 :     }
86 :    
87 :     STATIC_INLINE double lerp(double out_min, double out_max, double in_min, double in, double in_max)
88 :     {
89 :     return (out_min * (in_max - in) + out_max * (in - in_min)) / (in_max - in_min);
90 :     }
91 :    
92 :     STATIC_INLINE unsigned char ulerp(double in_min, double in, double in_max)
93 :     {
94 :     return (unsigned char)((0.0f * (in_max - in) + 255.0f * (in - in_min)) / (in_max - in_min));
95 :     }
96 :    
97 :     STATIC_INLINE REAL
98 :     alpha(float *data, REAL val, REAL gmag) {
99 :     unsigned int vi, gi;
100 :     REAL vv, gg, vf, gf;
101 :    
102 :     /* XXX HEY HACK */
103 :     vv = val/7.8431372549019605 + 3;
104 :     gg = gmag/3.5294117647058822 + 3;
105 :     vi = vv; vf = vv - vi;
106 :     gi = gg; gf = gg - gi;
107 :     return (data[vi+0 + 262*(gi+0)]*(1-vf)*(1-gf) +
108 :     data[vi+1 + 262*(gi+0)]*vf*(1-gf) +
109 :     data[vi+0 + 262*(gi+1)]*(1-vf)*gf +
110 :     data[vi+1 + 262*(gi+1)]*vf*gf);
111 :     }
112 :    
113 :     STATIC_INLINE REAL
114 :     shade(float *data, REAL val) {
115 :     unsigned int vi;
116 :     REAL vv, vf;
117 :    
118 :     /* XXX HEY HACK */
119 :     vv = val/0.010050251256281407 + 103;
120 :     vi = vv; vf = vv - vi;
121 :     return (data[vi]*(1-vf) + data[vi+1]*vf);
122 :     }
123 :    
124 :     STATIC_INLINE REAL
125 :     silho(float *data, REAL val) {
126 :     unsigned int vi;
127 :     REAL vv, vf;
128 :    
129 :     /* XXX HEY HACK */
130 :     vv = val/0.010050251256281407 + 3;
131 :     vi = vv; vf = vv - vi;
132 :     return (data[vi]*(1-vf) + data[vi+1]*vf);
133 :     }
134 :    
135 :     STATIC_INLINE REAL
136 :     rival(float *data, REAL rad, REAL the) {
137 :     unsigned int ri, ti;
138 :     REAL rr, tt, rf, tf;
139 :    
140 :     /* XXX HEY HACK */
141 :     rr = rad/0.0050251256281407036 + 3;
142 :     tt = (the + 3.0*AIR_PI/4.0)/0.015786894472361809 + 3;
143 :     ri = rr; rf = rr - ri;
144 :     ti = tt; tf = tt - ti;
145 :     return (data[ri+0 + 206*(ti+0)]*(1-rf)*(1-tf) +
146 :     data[ri+1 + 206*(ti+0)]*rf*(1-tf) +
147 :     data[ri+0 + 206*(ti+1)]*(1-rf)*tf +
148 :     data[ri+1 + 206*(ti+1)]*rf*tf);
149 :     }
150 :    
151 :     STATIC_INLINE void
152 :     depth(float *data, REAL rgb[3], REAL dep) {
153 :     unsigned int di;
154 :     REAL dd, df;
155 :    
156 :     /* XXX HEY HACK */
157 :     dd = dep/0.0050251256281407036 + 3;
158 :     di = dd; df = dd - di;
159 :     ELL_3V_SCALE_ADD2(rgb,
160 :     1-df, data + 3*di,
161 :     df, data + 3*(di+1));
162 :     }
163 :    
164 :     int main (int argc, const char *argv[])
165 :     {
166 :     const char *me = argv[0];
167 :     char *inFile = "../../data/vfrhand-nohip-smooth.nrrd",
168 :     *alphaFile = "../../data/txf/alpha-bone.nrrd",
169 :     *shadeFile = "../../data/txf/shade.nrrd",
170 :     *silhoFile = "../../data/txf/silho.nrrd",
171 :     *rivalFile = "../../data/txf/ridgvall.nrrd",
172 :     *depthFile = "../../data/txf/depth.nrrd";
173 :     float *alphaData = NULL, *shadeData = NULL, *silhoData = NULL,
174 :     *rivalData = NULL, *depthData = NULL;
175 :     char *outfile = "bmark-teem.nrrd";
176 :     double camEye[3] = {127.331, -1322.05, 272.53};
177 :     double camAt[3] = {63.0, 82.6536, 98.0};
178 :     double camUp[3] = {0.9987, 0.0459166, -0.0221267};
179 :     double camNear = -78.0;
180 :     double camFar = 78.0;
181 :     double camFOV = 5.0;
182 :    
183 :     int imgResU = 640;
184 :     int imgResV = 480;
185 :     double rayStep = 0.15;
186 :     /*
187 :     int imgResU = 480;
188 :     int imgResV = 360;
189 :     double rayStep = 1;
190 :     */
191 :    
192 :     double temp[3];
193 :     double temp2[3];
194 :     double temp3[3];
195 :     sub3(camAt, camEye, temp);
196 :     double camDist = mag3( temp );
197 :     double camVspNear = camDist + camNear;
198 :     double camVspFar = camDist + camFar;
199 :     double camN[3];
200 :     norm3(temp, camN);
201 :     double camU[3];
202 :     cross(camN, camUp, temp);
203 :     norm3(temp, camU);
204 :     double camV[3];
205 :     cross(camN, camU, camV);
206 :     double camVmax = tan(camFOV*AIR_PI/360.0)*camDist;
207 :     double camUmax = camVmax*(double)(imgResU)/(double)(imgResV);
208 :     double lightVspDir[3] = {-2.0, -3.0, -2.0};
209 :     scale_vec3(lightVspDir[0], camU, temp);
210 :     scale_vec3(lightVspDir[1], camV, temp2);
211 :     add3(temp, temp2, temp);
212 :     scale_vec3(lightVspDir[2], camN, temp2);
213 :     add3(temp, temp2, temp);
214 :     double lightDir[3];
215 :     norm3(temp, lightDir);
216 :     double phongKa = 0.05;
217 :     double phongKd = 0.80;
218 :     double phongKs = 0.20;
219 :     double phongSp = 50.0;
220 :     #define sthick 0.4
221 :     double refStep = 1.5;
222 :    
223 :     // Read in the data
224 :     char* err;
225 :     Nrrd *nin, *nalpha, *nshade, *nsilho, *nrival, *ndepth;
226 :     int status;
227 :     airArray *mop;
228 :     mop = airMopNew();
229 :    
230 :     nin = nrrdNew();
231 :     airMopAdd(mop, nin, (airMopper)nrrdNuke, airMopAlways);
232 :     nalpha = nrrdNew();
233 :     airMopAdd(mop, nalpha, (airMopper)nrrdNuke, airMopAlways);
234 :     nshade = nrrdNew();
235 :     airMopAdd(mop, nshade, (airMopper)nrrdNuke, airMopAlways);
236 :     nsilho = nrrdNew();
237 :     airMopAdd(mop, nsilho, (airMopper)nrrdNuke, airMopAlways);
238 :     nrival = nrrdNew();
239 :     airMopAdd(mop, nrival, (airMopper)nrrdNuke, airMopAlways);
240 :     ndepth = nrrdNew();
241 :     airMopAdd(mop, ndepth, (airMopper)nrrdNuke, airMopAlways);
242 :    
243 :     if (nrrdLoad(nin, inFile, NULL)
244 :     || nrrdLoad(nalpha, alphaFile, NULL)
245 :     || nrrdLoad(nshade, shadeFile, NULL)
246 :     || nrrdLoad(nsilho, silhoFile, NULL)
247 :     || nrrdLoad(nrival, rivalFile, NULL)
248 :     || nrrdLoad(ndepth, depthFile, NULL)) {
249 :     airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
250 :     fprintf(stderr, "%s: Trouble reading:\n%s", me, err);
251 :     airMopError(mop);
252 :     return -1;
253 :     }
254 :     alphaData = AIR_CAST(float *, nalpha->data); // XXX HEY HACK
255 :     shadeData = AIR_CAST(float *, nshade->data); // XXX HEY HACK
256 :     silhoData = AIR_CAST(float *, nsilho->data); // XXX HEY HACK
257 :     rivalData = AIR_CAST(float *, nrival->data); // XXX HEY HACK
258 :     depthData = AIR_CAST(float *, ndepth->data); // XXX HEY HACK
259 :    
260 :     gageContext *ctx;
261 :     gagePerVolume *pvl;
262 :     const double *valP, *gradP, *hessP;
263 :     double kparm[NRRD_KERNEL_PARMS_NUM] = {0.0};
264 :    
265 :     ctx = gageContextNew();
266 :     airMopAdd(mop, ctx, (airMopper)gageContextNix, airMopAlways);
267 :     if (!( pvl = gagePerVolumeNew(ctx, nin, gageKindScl) )) {
268 :     airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways);
269 :     fprintf(stderr, "%s: couldn't create gage context:\n%s", me, err);
270 :     airMopError(mop);
271 :     return -1;
272 :     }
273 :     if (gagePerVolumeAttach(ctx, pvl)
274 :     || gageKernelSet(ctx, gageKernel00, nrrdKernelBSpline3, kparm)
275 :     || gageKernelSet(ctx, gageKernel11, nrrdKernelBSpline3D, kparm)
276 :     || gageKernelSet(ctx, gageKernel22, nrrdKernelBSpline3DD, kparm)
277 :     || gageQueryItemOn(ctx, pvl, gageSclValue)
278 :     || gageQueryItemOn(ctx, pvl, gageSclGradVec)
279 :     || gageQueryItemOn(ctx, pvl, gageSclHessian)
280 :     || gageUpdate(ctx)) {
281 :     airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways);
282 :     fprintf(stderr, "%s: Trouble setting up gage context:\n%s", me, err);
283 :     airMopError(mop);
284 :     return -1;
285 :     }
286 :    
287 :     if (!( (valP = gageAnswerPointer(ctx, pvl, gageSclValue))
288 :     && (gradP = gageAnswerPointer(ctx, pvl, gageSclGradVec))
289 :     && (hessP = gageAnswerPointer(ctx, pvl, gageSclHessian)) )) {
290 :     fprintf(stderr, "%s: Trouble getting answer pointer\n", me);
291 :     airMopError(mop);
292 :     return -1;
293 :     }
294 :    
295 :     int ui, vi;
296 :     REAL rayU, rayV, rayN;
297 :     REAL rayVec[3];
298 :     REAL vv[3];
299 :     REAL rayPos[3];
300 :     REAL rayTransp, rayRGB[3];
301 :     REAL norm[3];
302 :     REAL ld;
303 :     REAL hd;
304 :     REAL mat;
305 :    
306 :     // allocate output image
307 :     Nrrd *nimg = nrrdNew();
308 :     airMopAdd(mop, nimg, (airMopper)nrrdNuke, airMopAlways);
309 :     if (nrrdAlloc_va(nimg, nrrdTypeDouble, 3,
310 :     AIR_CAST(size_t, 4),
311 :     AIR_CAST(size_t, imgResU),
312 :     AIR_CAST(size_t, imgResV))) {
313 :     airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
314 :     fprintf(stderr, "%s: Trouble allocating image output:\n%s", me, err);
315 :     airMopError(mop);
316 :     return -1;
317 :     }
318 :     double *out_data = AIR_CAST(double *, nimg->data);
319 :     /*
320 :     fprintf(stderr, "%s: hit return\n", me);
321 :     fgetc(stdin);
322 :     fprintf(stderr, "%s: here we go\n", me);
323 :     */
324 :     double t0 = airTime(); // start timing
325 :     for (vi = 0; vi < imgResV; vi++) {
326 :     rayV = lerp(-camVmax, camVmax, -0.5, (double)(vi), (double)(imgResV) - 0.5);
327 :     for (ui = 0; ui < imgResU; ui++) {
328 :     REAL len;
329 :     rayU = lerp(-camUmax, camUmax, -0.5, (double)(ui), (double)(imgResU) - 0.5);
330 :     ELL_3V_SCALE_ADD3(rayVec, 1.0, camN, rayU/camDist, camU, rayV/camDist, camV);
331 :     ELL_3V_SCALE(vv, -1, rayVec);
332 :     ELL_3V_NORM(vv, vv, len);
333 :     rayTransp = 1;
334 :     ELL_3V_SET(rayRGB, 0, 0, 0);
335 :    
336 :     for (rayN = camVspNear; rayN <= camVspFar; rayN += rayStep) {
337 :     ELL_3V_SCALE_ADD2(rayPos, 1.0, camEye, rayN, rayVec);
338 :    
339 :     if (inside(rayPos)) {
340 :     if (gageProbeSpace(ctx, rayPos[0], rayPos[1], rayPos[2],
341 :     AIR_FALSE /* index space */, AIR_TRUE /* clamp */)) {
342 :     fprintf(stderr, "%s: Trouble with gageProbe:\n%s", me, ctx->errStr);
343 :     airMopError(mop); return -1;
344 :     }
345 :     REAL gmag = ELL_3V_LEN(gradP);
346 :     REAL aa = alpha(alphaData, AIR_CLAMP(0, *valP, 2000), AIR_CLAMP(0, gmag, 900));
347 :     REAL matRGB[3];
348 :     if (aa > 0.0) {
349 :     REAL norm[3], gray, gten[9], proj[9], matA[9], matB[9];
350 :     ELL_3V_SCALE(norm, -1.0/gmag, gradP);
351 :     aa = 1.0 - pow(1.0-aa, rayStep/refStep);
352 :     gray = shade(shadeData, ELL_3V_DOT(norm, lightDir));
353 :     ELL_3MV_OUTER(matA, norm, norm);
354 :     ELL_3M_SCALE(proj, -1, matA);
355 :     proj[0] += 1; proj[4] += 1; proj[8] += 1;
356 :    
357 :     ELL_3M_MUL(matB, hessP, proj);
358 :     ELL_3M_MUL(gten, proj, matB);
359 :     ELL_3M_SCALE(gten, -1.0/gmag, gten);
360 :     REAL kv = ELL_3MV_CONTR(gten, vv)/ELL_3MV_CONTR(proj, vv);
361 :     kv = AIR_CLAMP(0.0, kv, 1.0/sthick);
362 :     REAL vdn = ELL_3V_DOT(vv, norm);
363 :     vdn = vdn*(1.0 + 2.0*vdn*vdn);
364 :     REAL tmp = sthick*(kv - 1.0/sthick);
365 :     REAL sarg = vdn*vdn + tmp*tmp;
366 :     gray *= silho(silhoData, AIR_CLAMP(0, sarg, 2));
367 :    
368 :     REAL gtns = ELL_3M_FROBSQRD(gten);
369 :     REAL gttr = ELL_3M_TRACE(gten);
370 :     REAL disc = SQRT(2.0*gtns - gttr*gttr);
371 :     REAL k1 = (gttr + disc)/2.0;
372 :     REAL k2 = (gttr - disc)/2.0;
373 :     tmp = SQRT(gtns);
374 :     gray *= rival(rivalData, AIR_MIN(0.55, tmp)/0.55, ATAN2(k2,k1));
375 :    
376 :     depth(depthData, matRGB, lerp(0.0, 1.0, camVspNear, rayN, camVspFar));
377 :    
378 :     ELL_3V_SCALE_INCR(rayRGB, rayTransp*aa*gray, matRGB);
379 :     rayTransp *= 1 - aa;
380 :     }
381 :     }
382 :     if (rayTransp < 0.01) {
383 :     rayTransp = 0;
384 :     break;
385 :     }
386 :     }
387 :    
388 :     ELL_4V_SET(out_data + 4*(ui + imgResU*vi),
389 :     0.9*rayRGB[0], 0.9*rayRGB[1], 0.9*rayRGB[2], 1.0-rayTransp);
390 :    
391 :     }
392 :     }
393 :    
394 :     double totalTime = airTime() - t0;
395 :     printf("usr=%f\n", totalTime);
396 :    
397 :     if (nrrdSave(outfile, nimg, NULL)) {
398 :     airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
399 :     fprintf(stderr, "%s: Trouble writing output:\n%s", me, err);
400 :     airMopError(mop);
401 :     return -1;
402 :     }
403 :    
404 :     airMopOkay(mop);
405 :     return 0;
406 :    
407 :     }

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