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 1606 - (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 : glk 1605 // XXX - Hack specific to ../../data/vfrhand-nohip.nhdr"
69 :     // but for some reason still not matching what Diderot says
70 :     if(pos[0] < 2*0.9375 || pos[0] > (176-3)*0.9375) return false;
71 :     if(pos[1] < 2*0.9375 || pos[1] > (187-3)*0.9375) return false;
72 :     if(pos[2] < 2.0 || pos[2] > (190-3)) return false;
73 : jhr 1533 return true;
74 :     }
75 :    
76 :     STATIC_INLINE double lerp(double out_min, double out_max, double in_min, double in, double in_max)
77 :     {
78 :     return (out_min * (in_max - in) + out_max * (in - in_min)) / (in_max - in_min);
79 :     }
80 :    
81 :     STATIC_INLINE unsigned char ulerp(double in_min, double in, double in_max)
82 :     {
83 :     return (unsigned char)((0.0f * (in_max - in) + 255.0f * (in - in_min)) / (in_max - in_min));
84 :     }
85 :    
86 : glk 1600 int main (int argc, const char *argv[])
87 : jhr 1533 {
88 : glk 1600 const char *me = argv[0];
89 : jhr 1533 char *infile = "../../data/vfrhand-nohip.nhdr";
90 : glk 1605 char *outfile = "bmark-teem.nrrd";
91 : jhr 1533 double camEye[3] = {127.331, -1322.05, 272.53};
92 :     double camAt[3] = {63.0, 82.6536, 98.0};
93 :     double camUp[3] = {0.9987, 0.0459166, -0.0221267};
94 :     double camNear = -78.0;
95 :     double camFar = 78.0;
96 :     double camFOV = 5.0;
97 :     int imgResU = 480;
98 :     int imgResV = 345;
99 : glk 1605 double rayStep = 0.1;
100 : jhr 1533 double valOpacMin = 400.0; // 400.0 for skin, 1150.0 for bone
101 :     double valOpacMax = 700.0; // 700.0 for skin, 1450.0 for bone
102 : glk 1605
103 :    
104 : jhr 1533 double temp[3];
105 :     double temp2[3];
106 :     double temp3[3];
107 :     sub3(camAt, camEye, temp);
108 :     double camDist = mag3( temp );
109 :     double camVspNear = camDist + camNear;
110 :     double camVspFar = camDist + camFar;
111 :     double camN[3];
112 :     norm3(temp, camN);
113 :     double camU[3];
114 :     cross(camN, camUp, temp);
115 :     norm3(temp, camU);
116 :     double camV[3];
117 :     cross(camN, camU, camV);
118 : glk 1600 double camVmax = tan(camFOV*AIR_PI/360.0)*camDist;
119 : jhr 1533 double camUmax = camVmax*(double)(imgResU)/(double)(imgResV);
120 :     double lightVspDir[3] = {0.9, -1.0, -2.5};
121 :     scale_vec3(lightVspDir[0], camU, temp);
122 :     scale_vec3(lightVspDir[1], camV, temp2);
123 :     add3(temp, temp2, temp);
124 :     scale_vec3(lightVspDir[2], camN, temp2);
125 :     add3(temp, temp2, temp);
126 :     double lightDir[3];
127 :     norm3(temp, lightDir);
128 :     double phongKa = 0.05;
129 :     double phongKd = 0.80;
130 :     double phongKs = 0.20;
131 : glk 1600 double phongSp = 50.0;
132 : jhr 1533
133 :     // Read in the data
134 :     char* err;
135 :     Nrrd *nin;
136 :     int status;
137 : glk 1600 airArray *mop;
138 :     mop = airMopNew();
139 :    
140 : jhr 1533 nin = nrrdNew();
141 : glk 1600 airMopAdd(mop, nin, (airMopper)nrrdNuke, airMopAlways);
142 : jhr 1533
143 : glk 1605 // Get data for mip
144 : glk 1600 if (nrrdLoad(nin, infile, NULL)) {
145 :     airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
146 : glk 1605 fprintf(stderr, "%s: Trouble reading \"%s\":\n%s", me, infile, err);
147 : glk 1600 airMopError(mop);
148 : jhr 1533 return -1;
149 :     }
150 :    
151 :     gageContext *ctx;
152 :     gagePerVolume *pvl;
153 :     const double *val;
154 :     const double *grad;
155 : glk 1600 double kparm[NRRD_KERNEL_PARMS_NUM] = {0.0};
156 : jhr 1533
157 :     ctx = gageContextNew();
158 : glk 1600 airMopAdd(mop, ctx, (airMopper)gageContextNix, airMopAlways);
159 : glk 1605 if (!( pvl = gagePerVolumeNew(ctx, nin, gageKindScl) )) {
160 :     airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways);
161 :     fprintf(stderr, "%s: couldn't create gage context:\n%s", me, err);
162 : glk 1600 airMopError(mop);
163 : jhr 1533 return -1;
164 :     }
165 : glk 1600 if (gagePerVolumeAttach(ctx, pvl)
166 :     || gageKernelSet(ctx, gageKernel00, nrrdKernelBSpline3, kparm)
167 :     || gageKernelSet(ctx, gageKernel11, nrrdKernelBSpline3D, kparm)
168 :     || gageQueryItemOn(ctx, pvl, gageSclValue)
169 :     || gageQueryItemOn(ctx, pvl, gageSclGradVec)
170 :     || gageUpdate(ctx)) {
171 :     airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways);
172 : glk 1605 fprintf(stderr, "%s: Trouble setting up gage context:\n%s", me, err);
173 : glk 1600 airMopError(mop);
174 : jhr 1533 return -1;
175 :     }
176 :    
177 : glk 1600 if (!( (val = gageAnswerPointer(ctx, pvl, gageSclValue))
178 :     && (grad = gageAnswerPointer(ctx, pvl, gageSclGradVec)) )) {
179 : glk 1605 fprintf(stderr, "%s: Trouble getting answer pointer\n", me);
180 : glk 1600 airMopError(mop);
181 : jhr 1533 return -1;
182 :     }
183 :    
184 :     int ui, vi;
185 :     double rayU, rayV, rayN;
186 :     double rayVec[3];
187 : glk 1605 double toEye[3];
188 : jhr 1533 double rayPos[3];
189 :     double transp, gray;
190 :     double norm[3];
191 :     double alpha;
192 :     double ld;
193 :     double hd;
194 :     double mat;
195 : glk 1600
196 : glk 1605 // allocate output image
197 : glk 1600 Nrrd *nimg = nrrdNew();
198 :     airMopAdd(mop, nimg, (airMopper)nrrdNuke, airMopAlways);
199 :     if (nrrdAlloc_va(nimg, nrrdTypeDouble, 3,
200 :     AIR_CAST(size_t, 4),
201 :     AIR_CAST(size_t, imgResU),
202 : glk 1605 AIR_CAST(size_t, imgResV))) {
203 : glk 1600 airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
204 : glk 1605 fprintf(stderr, "%s: Trouble allocating image output:\n%s", me, err);
205 : glk 1600 airMopError(mop);
206 : jhr 1533 return -1;
207 :     }
208 : glk 1600 double *out_data = AIR_CAST(double *, nimg->data);
209 : glk 1606 /*
210 : glk 1605 fprintf(stderr, "hit return\n");
211 :     fgetc(stdin);
212 :     fprintf(stderr, "here we go\n");
213 : glk 1606 */
214 : glk 1600 double t0 = airTime(); // start timing
215 : jhr 1575 for(vi = 0; vi < imgResV; vi++) {
216 : jhr 1533 rayV = lerp(-camVmax, camVmax, -0.5, (double)(vi), (double)(imgResV) - 0.5);
217 : jhr 1575 for(ui = 0; ui < imgResU; ui++) {
218 : glk 1605 double len;
219 : jhr 1533 rayU = lerp(-camUmax, camUmax, -0.5, (double)(ui), (double)(imgResU) - 0.5);
220 : glk 1605 ELL_3V_SCALE_ADD3(rayVec, 1.0, camN, rayU/camDist, camU, rayV/camDist, camV);
221 :     ELL_3V_SCALE(toEye, -1, rayVec);
222 :     ELL_3V_NORM(toEye, toEye, len);
223 : jhr 1533 transp = 1;
224 :     gray = 0;
225 : glk 1605
226 : glk 1600 for (rayN = camVspNear; rayN <= camVspFar; rayN += rayStep) {
227 : glk 1605 ELL_3V_SCALE_ADD2(rayPos, 1.0, camEye, rayN, rayVec);
228 :    
229 : glk 1600 if (inside(rayPos)) {
230 : glk 1605 if (gageProbeSpace(ctx, rayPos[0], rayPos[1], rayPos[2],
231 :     AIR_FALSE /* index space */, AIR_TRUE /* clamp */)) {
232 :     fprintf(stderr, "%s: Trouble with gageProbe:\n%s", me, ctx->errStr);
233 : glk 1600 airMopError(mop); return -1;
234 : jhr 1533 }
235 : glk 1600 if (*val > valOpacMin) {
236 : glk 1605 ELL_3V_SCALE(temp2, -1, grad);
237 :     ELL_3V_NORM(norm, temp2, len);
238 : jhr 1533 alpha = lerp(0.0, 1.0, valOpacMin, *val, valOpacMax);
239 :     alpha = (alpha < 1.0) ? alpha : 1;
240 :     ld = dot3(norm, lightDir);
241 :     ld = (ld < 0) ? 0 : ld;
242 : glk 1605 ELL_3V_ADD2(temp2, toEye, lightDir);
243 :     ELL_3V_NORM(temp2, temp2, len);
244 : jhr 1533 hd = dot3(norm, temp2);
245 :     hd = (hd < 0) ? 0 : hd;
246 : glk 1605 mat = (phongKa +
247 :     phongKd * ld +
248 :     phongKs * pow(hd, phongSp));
249 : jhr 1533 gray = gray + transp * alpha * mat;
250 :     transp *= 1 - alpha;
251 :     }
252 :     }
253 : glk 1605 if (transp < 0.01) {
254 : jhr 1533 transp = 0;
255 :     break;
256 :     }
257 :     }
258 :    
259 : glk 1605 ELL_4V_SET(out_data + 4*(ui + imgResU*vi), gray, gray, gray, 1.0-transp);
260 :    
261 : jhr 1533 }
262 :     }
263 : jhr 1575
264 : glk 1600 double totalTime = airTime() - t0;
265 : jhr 1533 printf("usr=%f\n", totalTime);
266 :    
267 : glk 1605 if (nrrdSave(outfile, nimg, NULL)) {
268 : glk 1600 airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
269 : glk 1605 fprintf(stderr, "%s: Trouble writing output:\n%s", me, err);
270 : glk 1600 airMopError(mop);
271 : jhr 1533 return -1;
272 :     }
273 :    
274 : glk 1600 airMopOkay(mop);
275 : jhr 1533 return 0;
276 :    
277 :     }

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