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

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