1 /*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 2012, 2011, 2010, 2009 University of Chicago
4 Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
5 Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
6
7 This library is free software; you can redistribute it and/or
8 modify it under the terms of the GNU Lesser General Public License
9 (LGPL) as published by the Free Software Foundation; either
10 version 2.1 of the License, or (at your option) any later version.
11 The terms of redistributing and/or modifying this software also
12 include exceptions to the LGPL that facilitate static linking.
13
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with this library; if not, write to Free Software Foundation, Inc.,
21 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22 */
23
24 #include "hoover.h"
25
26 /*
27 ** learned: if you're going to "simplify" code which computes some
28 ** floating point value within a loop using AFFINE() on the loop
29 ** control variable, by simply incrementing that value with the
30 ** correct amount iteration, BE SURE THAT THE INCREMENTING IS DONE in
31 ** every possible control path of the loop (wasn't incrementing ray
32 ** sample position if first sample wasn't inside the volume)
33 */
34
35 /*
36 ** _hooverLearnLengths()
37 **
38 ** This is where we enforce the constraint that the volume always fit
39 ** inside a cube with edge length 2, centered at the origin.
40 **
41 ** volHLen[i] is the HALF the length of the volume along axis i
42 **
43 ** NOTE: none of this comes into play if we have ctx->shape
44 */
45 void
_hooverLearnLengths(double volHLen[3],double voxLen[3],hooverContext * ctx)46 _hooverLearnLengths(double volHLen[3], double voxLen[3], hooverContext *ctx) {
47 double maxLen;
48 int numSamples[3], numElements[3];
49
50 ELL_3V_COPY(numSamples, ctx->volSize);
51 if (nrrdCenterNode == ctx->volCentering) {
52 numElements[0] = numSamples[0]-1;
53 numElements[1] = numSamples[1]-1;
54 numElements[2] = numSamples[2]-1;
55 } else {
56 numElements[0] = numSamples[0];
57 numElements[1] = numSamples[1];
58 numElements[2] = numSamples[2];
59 }
60 volHLen[0] = numElements[0]*ctx->volSpacing[0];
61 volHLen[1] = numElements[1]*ctx->volSpacing[1];
62 volHLen[2] = numElements[2]*ctx->volSpacing[2];
63 maxLen = AIR_MAX(volHLen[0], volHLen[1]);
64 maxLen = AIR_MAX(volHLen[2], maxLen);
65 volHLen[0] /= maxLen;
66 volHLen[1] /= maxLen;
67 volHLen[2] /= maxLen;
68 voxLen[0] = 2*volHLen[0]/numElements[0];
69 voxLen[1] = 2*volHLen[1]/numElements[1];
70 voxLen[2] = 2*volHLen[2]/numElements[2];
71 }
72
73 /*
74 ** _hooverExtraContext struct
75 **
76 ** Like hooverContext, this is READ-ONLY information which is not specific
77 ** to any thread.
78 ** Unlike hooverContext, it is solely for the benefit of the calculations
79 ** done in _hooverThreadBody.
80 **
81 ** No one outside hoover should need to know about this.
82 */
83 typedef struct {
84 double volHLen[3], /* length of x,y,z edges of volume bounding box */
85 voxLen[3], /* length of x,y,z edges of voxels */
86 uBase, uCap, /* uMin and uMax as seen on the near cutting plane */
87 vBase, vCap, /* analogous to uBase and uCap */
88 rayZero[3]; /* location of near plane, line of sight interxion */
89 } _hooverExtraContext;
90
91 _hooverExtraContext *
_hooverExtraContextNew(hooverContext * ctx)92 _hooverExtraContextNew(hooverContext *ctx) {
93 _hooverExtraContext *ec;
94
95 ec = (_hooverExtraContext *)calloc(1, sizeof(_hooverExtraContext));
96 if (ec) {
97 if (ctx->shape) {
98 ELL_3V_NAN_SET(ec->volHLen);
99 ELL_3V_NAN_SET(ec->voxLen);
100 } else {
101 _hooverLearnLengths(ec->volHLen, ec->voxLen, ctx);
102 }
103 ELL_3V_SCALE_ADD2(ec->rayZero,
104 1.0, ctx->cam->from,
105 ctx->cam->vspNeer, ctx->cam->N);
106 }
107 return ec;
108 }
109
110 _hooverExtraContext *
_hooverExtraContextNix(_hooverExtraContext * ec)111 _hooverExtraContextNix(_hooverExtraContext *ec) {
112
113 if (ec) {
114 free(ec);
115 }
116 return NULL;
117 }
118
119 /*
120 ** _hooverThreadArg struct
121 **
122 ** A pointer to this is passed to _hooverThreadBody. It contains all the
123 ** information which is not thread-specific, and all the thread-specific
124 ** information known at the level of hooverRender.
125 **
126 ** For simplicity sake, a pointer to a struct of this type is also
127 ** returned from _hooverThreadBody, so this is where we store an
128 ** error-signaling return value (errCode), and what function had
129 ** trouble (whichErr).
130 */
131 typedef struct {
132 /* ----------------------- input */
133 hooverContext *ctx;
134 _hooverExtraContext *ec;
135 void *render;
136 int whichThread;
137 /* ----------------------- output */
138 int whichErr;
139 int errCode;
140 } _hooverThreadArg;
141
142 void *
_hooverThreadBody(void * _arg)143 _hooverThreadBody(void *_arg) {
144 _hooverThreadArg *arg;
145 void *thread;
146 int ret, /* to catch return values from callbacks */
147 sampleI, /* which sample we're on */
148 inside, /* we're inside the volume */
149 vI, uI; /* integral coords in image */
150 double tmp,
151 mm, /* lowest position in index space, for all axes */
152 Mx, My, Mz, /* highest position in index space on each axis */
153 u, v, /* floating-point coords in image */
154 uvScale, /* how to scale (u,v) to go from image to
155 near plane, according to ortho or perspective */
156 lx, ly, lz, /* half edge-lengths of volume */
157 rayLen=0, /* length of segment formed by ray line intersecting
158 the near and far clipping planes */
159 rayT, /* current position along ray (world-space) */
160 rayDirW[3], /* unit-length ray direction (world-space) */
161 rayDirI[3], /* rayDirW transformed into index space;
162 not unit length, but a unit change in
163 world space along rayDirW translates to
164 this change in index space along rayDirI */
165 rayPosW[3], /* current ray location (world-space) */
166 rayPosI[3], /* current ray location (index-space) */
167 rayStartW[3], /* ray start on near plane (world-space) */
168 rayStartI[3], /* ray start on near plane (index-space) */
169 rayStep, /* distance between samples (world-space) */
170 vOff[3], uOff[3]; /* offsets in arg->ec->wU and arg->ec->wV
171 directions towards start of ray */
172
173 arg = (_hooverThreadArg *)_arg;
174 if ( (ret = (arg->ctx->threadBegin)(&thread,
175 arg->render,
176 arg->ctx->user,
177 arg->whichThread)) ) {
178 arg->errCode = ret;
179 arg->whichErr = hooverErrThreadBegin;
180 return arg;
181 }
182 if (arg->ctx->shape) {
183 lx = ly = lz = AIR_NAN;
184 if (nrrdCenterNode == arg->ctx->shape->center) {
185 mm = 0;
186 Mx = arg->ctx->shape->size[0]-1;
187 My = arg->ctx->shape->size[1]-1;
188 Mz = arg->ctx->shape->size[2]-1;
189 } else {
190 mm = -0.5;
191 Mx = arg->ctx->shape->size[0]-0.5;
192 My = arg->ctx->shape->size[1]-0.5;
193 Mz = arg->ctx->shape->size[2]-0.5;
194 }
195 } else {
196 lx = arg->ec->volHLen[0];
197 ly = arg->ec->volHLen[1];
198 lz = arg->ec->volHLen[2];
199 if (nrrdCenterNode == arg->ctx->volCentering) {
200 mm = 0;
201 Mx = arg->ctx->volSize[0]-1;
202 My = arg->ctx->volSize[1]-1;
203 Mz = arg->ctx->volSize[2]-1;
204 } else {
205 mm = -0.5;
206 Mx = arg->ctx->volSize[0]-0.5;
207 My = arg->ctx->volSize[1]-0.5;
208 Mz = arg->ctx->volSize[2]-0.5;
209 }
210 }
211
212 if (arg->ctx->cam->orthographic) {
213 ELL_3V_COPY(rayDirW, arg->ctx->cam->N);
214 if (arg->ctx->shape) {
215 double zeroW[3], zeroI[3];
216 ELL_3V_SET(zeroW, 0, 0, 0);
217 gageShapeWtoI(arg->ctx->shape, zeroI, zeroW);
218 gageShapeWtoI(arg->ctx->shape, rayDirI, rayDirW);
219 ELL_3V_SUB(rayDirI, rayDirI, zeroI);
220 } else {
221 rayDirI[0] = AIR_DELTA(-lx, rayDirW[0], lx, mm, Mx);
222 rayDirI[1] = AIR_DELTA(-ly, rayDirW[1], ly, mm, My);
223 rayDirI[2] = AIR_DELTA(-lz, rayDirW[2], lz, mm, Mz);
224 }
225 rayLen = arg->ctx->cam->vspFaar - arg->ctx->cam->vspNeer;
226 uvScale = 1.0;
227 } else {
228 uvScale = arg->ctx->cam->vspNeer/arg->ctx->cam->vspDist;
229 }
230
231 while (1) {
232 /* the work assignment is simply the next scanline to be rendered:
233 the result of all this is setting vI */
234 if (arg->ctx->workMutex) {
235 airThreadMutexLock(arg->ctx->workMutex);
236 }
237 vI = arg->ctx->workIdx;
238 if (arg->ctx->workIdx < arg->ctx->imgSize[1]) {
239 arg->ctx->workIdx += 1;
240 }
241 if (arg->ctx->workMutex) {
242 airThreadMutexUnlock(arg->ctx->workMutex);
243 }
244 if (vI == arg->ctx->imgSize[1]) {
245 /* we're done! */
246 break;
247 }
248
249 if (nrrdCenterCell == arg->ctx->imgCentering) {
250 v = uvScale*AIR_AFFINE(-0.5, vI, arg->ctx->imgSize[1]-0.5,
251 arg->ctx->cam->vRange[0],
252 arg->ctx->cam->vRange[1]);
253 } else {
254 v = uvScale*AIR_AFFINE(0.0, vI, arg->ctx->imgSize[1]-1.0,
255 arg->ctx->cam->vRange[0],
256 arg->ctx->cam->vRange[1]);
257 }
258 ELL_3V_SCALE(vOff, v, arg->ctx->cam->V);
259 for (uI=0; uI<arg->ctx->imgSize[0]; uI++) {
260 if (nrrdCenterCell == arg->ctx->imgCentering) {
261 u = uvScale*AIR_AFFINE(-0.5, uI, arg->ctx->imgSize[0]-0.5,
262 arg->ctx->cam->uRange[0],
263 arg->ctx->cam->uRange[1]);
264 } else {
265 u = uvScale*AIR_AFFINE(0.0, uI, arg->ctx->imgSize[0]-1.0,
266 arg->ctx->cam->uRange[0],
267 arg->ctx->cam->uRange[1]);
268 }
269 ELL_3V_SCALE(uOff, u, arg->ctx->cam->U);
270 ELL_3V_ADD3(rayStartW, uOff, vOff, arg->ec->rayZero);
271 if (arg->ctx->shape) {
272 gageShapeWtoI(arg->ctx->shape, rayStartI, rayStartW);
273 } else {
274 rayStartI[0] = AIR_AFFINE(-lx, rayStartW[0], lx, mm, Mx);
275 rayStartI[1] = AIR_AFFINE(-ly, rayStartW[1], ly, mm, My);
276 rayStartI[2] = AIR_AFFINE(-lz, rayStartW[2], lz, mm, Mz);
277 }
278 if (!arg->ctx->cam->orthographic) {
279 ELL_3V_SUB(rayDirW, rayStartW, arg->ctx->cam->from);
280 ELL_3V_NORM(rayDirW, rayDirW, tmp);
281 if (arg->ctx->shape) {
282 double zeroW[3], zeroI[3];
283 ELL_3V_SET(zeroW, 0, 0, 0);
284 gageShapeWtoI(arg->ctx->shape, zeroI, zeroW);
285 gageShapeWtoI(arg->ctx->shape, rayDirI, rayDirW);
286 ELL_3V_SUB(rayDirI, rayDirI, zeroI);
287 } else {
288 rayDirI[0] = AIR_DELTA(-lx, rayDirW[0], lx, mm, Mx);
289 rayDirI[1] = AIR_DELTA(-ly, rayDirW[1], ly, mm, My);
290 rayDirI[2] = AIR_DELTA(-lz, rayDirW[2], lz, mm, Mz);
291 }
292 rayLen = ((arg->ctx->cam->vspFaar - arg->ctx->cam->vspNeer)/
293 ELL_3V_DOT(rayDirW, arg->ctx->cam->N));
294 }
295 if ( (ret = (arg->ctx->rayBegin)(thread,
296 arg->render,
297 arg->ctx->user,
298 uI, vI, rayLen,
299 rayStartW, rayStartI,
300 rayDirW, rayDirI)) ) {
301 arg->errCode = ret;
302 arg->whichErr = hooverErrRayBegin;
303 return arg;
304 }
305
306 sampleI = 0;
307 rayT = 0;
308 while (1) {
309 ELL_3V_SCALE_ADD2(rayPosW, 1.0, rayStartW, rayT, rayDirW);
310 if (arg->ctx->shape) {
311 gageShapeWtoI(arg->ctx->shape, rayPosI, rayPosW);
312 } else {
313 ELL_3V_SCALE_ADD2(rayPosI, 1.0, rayStartI, rayT, rayDirI);
314 }
315 inside = (AIR_IN_CL(mm, rayPosI[0], Mx) &&
316 AIR_IN_CL(mm, rayPosI[1], My) &&
317 AIR_IN_CL(mm, rayPosI[2], Mz));
318 rayStep = (arg->ctx->sample)(thread,
319 arg->render,
320 arg->ctx->user,
321 sampleI, rayT,
322 inside,
323 rayPosW, rayPosI);
324 if (!AIR_EXISTS(rayStep)) {
325 /* sampling failed */
326 arg->errCode = 0;
327 arg->whichErr = hooverErrSample;
328 return arg;
329 }
330 if (!rayStep) {
331 /* ray decided to finish itself */
332 break;
333 }
334 /* else we moved to a new location along the ray */
335 rayT += rayStep;
336 if (!AIR_IN_CL(0, rayT, rayLen)) {
337 /* ray stepped outside near-far clipping region, its done. */
338 break;
339 }
340 sampleI++;
341 }
342
343 if ( (ret = (arg->ctx->rayEnd)(thread,
344 arg->render,
345 arg->ctx->user)) ) {
346 arg->errCode = ret;
347 arg->whichErr = hooverErrRayEnd;
348 return arg;
349 }
350 } /* end this scanline */
351 } /* end while(1) assignment of scanlines */
352
353 if ( (ret = (arg->ctx->threadEnd)(thread,
354 arg->render,
355 arg->ctx->user)) ) {
356 arg->errCode = ret;
357 arg->whichErr = hooverErrThreadEnd;
358 return arg;
359 }
360
361 /* returning NULL actually indicates that there was NOT an error */
362 return NULL;
363 }
364
365 typedef union {
366 _hooverThreadArg **h;
367 void **v;
368 } _htpu;
369
370 /*
371 ******** hooverRender()
372 **
373 ** because of the biff usage(), only one thread can call hooverRender(),
374 ** and no promises if the threads themselves call biff...
375 */
376 int
hooverRender(hooverContext * ctx,int * errCodeP,int * errThreadP)377 hooverRender(hooverContext *ctx, int *errCodeP, int *errThreadP) {
378 static const char me[]="hooverRender";
379 _hooverExtraContext *ec;
380 _hooverThreadArg args[HOOVER_THREAD_MAX];
381 _hooverThreadArg *errArg;
382 airThread *thread[HOOVER_THREAD_MAX];
383 _htpu u;
384
385 void *render;
386 int ret;
387 airArray *mop;
388 unsigned int threadIdx;
389
390 if (!( errCodeP && errThreadP )) {
391 biffAddf(HOOVER, "%s: got NULL int return pointer", me);
392 return hooverErrInit;
393 }
394
395 /* this calls limnCameraUpdate() */
396 if (hooverContextCheck(ctx)) {
397 biffAddf(HOOVER, "%s: problem detected in given context", me);
398 *errCodeP = 0;
399 *errThreadP = 0;
400 return hooverErrInit;
401 }
402
403 if (!(ec = _hooverExtraContextNew(ctx))) {
404 biffAddf(HOOVER, "%s: problem creating thread context", me);
405 *errCodeP = 0;
406 *errThreadP = 0;
407 return hooverErrInit;
408 }
409 mop = airMopNew();
410 airMopAdd(mop, ec, (airMopper)_hooverExtraContextNix, airMopAlways);
411 if ( (ret = (ctx->renderBegin)(&render, ctx->user)) ) {
412 *errCodeP = ret;
413 *errCodeP = 0;
414 *errThreadP = 0;
415 airMopError(mop);
416 return hooverErrRenderBegin;
417 }
418
419 for (threadIdx=0; threadIdx<ctx->numThreads; threadIdx++) {
420 args[threadIdx].ctx = ctx;
421 args[threadIdx].ec = ec;
422 args[threadIdx].render = render;
423 args[threadIdx].whichThread = threadIdx;
424 args[threadIdx].whichErr = hooverErrNone;
425 args[threadIdx].errCode = 0;
426 thread[threadIdx] = airThreadNew();
427 }
428 ctx->workIdx = 0;
429 if (1 < ctx->numThreads) {
430 ctx->workMutex = airThreadMutexNew();
431 } else {
432 ctx->workMutex = NULL;
433 }
434
435 /* (done): call airThreadStart() once per thread, passing the
436 address of a distinct (and appropriately intialized)
437 _hooverThreadArg to each. If return of airThreadStart() is
438 non-zero, put its return in *errCodeP, the number of the
439 problematic in *errThreadP, and return hooverErrThreadCreate.
440 Then call airThreadJoin() on all the threads, passing &errArg as
441 "retval". On non-zero return, set *errCodeP and *errThreadP,
442 and return hooverErrThreadJoin. If return of airThreadJoin() is
443 zero, but the errArg is non-NULL, then assume that this errArg
444 is actually just the passed _hooverThreadArg returned to us, and
445 from this copy errArg->errCode into *errCodeP, and return
446 errArg->whichErr */
447
448 if (1 < ctx->numThreads && !airThreadCapable) {
449 fprintf(stderr, "%s: WARNING: not multi-threaded; will do %d "
450 "\"threads\" serially !!!\n", me, ctx->numThreads);
451 }
452
453 for (threadIdx=0; threadIdx<ctx->numThreads; threadIdx++) {
454 if ((ret = airThreadStart(thread[threadIdx], _hooverThreadBody,
455 (void *) &args[threadIdx]))) {
456 *errCodeP = ret;
457 *errThreadP = threadIdx;
458 airMopError(mop);
459 return hooverErrThreadCreate;
460 }
461 }
462
463 for (threadIdx=0; threadIdx<ctx->numThreads; threadIdx++) {
464 u.h = &errArg;
465 if ((ret = airThreadJoin(thread[threadIdx], u.v))) {
466 *errCodeP = ret;
467 *errThreadP = threadIdx;
468 airMopError(mop);
469 return hooverErrThreadJoin;
470 }
471 if (errArg != NULL) {
472 *errCodeP = errArg->errCode;
473 *errThreadP = threadIdx;
474 return errArg->whichErr;
475 }
476 thread[threadIdx] = airThreadNix(thread[threadIdx]);
477 }
478
479 if (1 < ctx->numThreads) {
480 ctx->workMutex = airThreadMutexNix(ctx->workMutex);
481 }
482
483 if ( (ret = (ctx->renderEnd)(render, ctx->user)) ) {
484 *errCodeP = ret;
485 *errThreadP = -1;
486 return hooverErrRenderEnd;
487 }
488 render = NULL;
489 airMopOkay(mop);
490
491 *errCodeP = 0;
492 *errThreadP = 0;
493 return hooverErrNone;
494 }
495