Annotation of /benchmarks/programs/vr-lite-cam/bmark-teem.c

Parent Directory | 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