Head

GCC Code Coverage Report

Directory:	./		Exec	Total	Coverage
File:	src/hoover/rays.c	Lines:	0	226	0.0 %
Date:	2017-05-26	Branches:	0	112	0.0 %


/*
  Teem: Tools to process and visualize scientific data and images             .
  Copyright (C) 2013, 2012, 2011, 2010, 2009  University of Chicago
  Copyright (C) 2008, 2007, 2006, 2005  Gordon Kindlmann
  Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998  University of Utah

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public License
  (LGPL) as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  The terms of redistributing and/or modifying this software also
  include exceptions to the LGPL that facilitate static linking.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public License
  along with this library; if not, write to Free Software Foundation, Inc.,
  51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
*/

#include "hoover.h"

/*
** learned: if you're going to "simplify" code which computes some
** floating point value within a loop using AFFINE() on the loop
** control variable, by simply incrementing that value with the
** correct amount iteration, BE SURE THAT THE INCREMENTING IS DONE in
** every possible control path of the loop (wasn't incrementing ray
** sample position if first sample wasn't inside the volume)
*/

/*
** _hooverLearnLengths()
**
** This is where we enforce the constraint that the volume always fit
** inside a cube with edge length 2, centered at the origin.
**
** volHLen[i] is the HALF the length of the volume along axis i
**
** NOTE: none of this comes into play if we have ctx->shape
*/
void
_hooverLearnLengths(double volHLen[3], double voxLen[3], hooverContext *ctx) {
  double maxLen;
  int numSamples[3], numElements[3];

  ELL_3V_COPY(numSamples, ctx->volSize);
  if (nrrdCenterNode == ctx->volCentering) {
    numElements[0] = numSamples[0]-1;
    numElements[1] = numSamples[1]-1;
    numElements[2] = numSamples[2]-1;
  } else {
    numElements[0] = numSamples[0];
    numElements[1] = numSamples[1];
    numElements[2] = numSamples[2];
  }
  volHLen[0] = numElements[0]*ctx->volSpacing[0];
  volHLen[1] = numElements[1]*ctx->volSpacing[1];
  volHLen[2] = numElements[2]*ctx->volSpacing[2];
  maxLen = AIR_MAX(volHLen[0], volHLen[1]);
  maxLen = AIR_MAX(volHLen[2], maxLen);
  volHLen[0] /= maxLen;
  volHLen[1] /= maxLen;
  volHLen[2] /= maxLen;
  voxLen[0] = 2*volHLen[0]/numElements[0];
  voxLen[1] = 2*volHLen[1]/numElements[1];
  voxLen[2] = 2*volHLen[2]/numElements[2];
}

/*
** _hooverExtraContext struct
**
** Like hooverContext, this is READ-ONLY information which is not specific
** to any thread.
** Unlike hooverContext, it is solely for the benefit of the calculations
** done in _hooverThreadBody.
**
** No one outside hoover should need to know about this.
*/
typedef struct {
  double volHLen[3],     /* length of x,y,z edges of volume bounding box */
    voxLen[3],           /* length of x,y,z edges of voxels */
    uBase, uCap,         /* uMin and uMax as seen on the near cutting plane */
    vBase, vCap,         /* analogous to uBase and uCap */
    rayZero[3];          /* location of near plane, line of sight interxion */
} _hooverExtraContext;

_hooverExtraContext *
_hooverExtraContextNew(hooverContext *ctx) {
  _hooverExtraContext *ec;

  ec = (_hooverExtraContext *)calloc(1, sizeof(_hooverExtraContext));
  if (ec) {
    if (ctx->shape) {
      ELL_3V_NAN_SET(ec->volHLen);
      ELL_3V_NAN_SET(ec->voxLen);
    } else {
      _hooverLearnLengths(ec->volHLen, ec->voxLen, ctx);
    }
    ELL_3V_SCALE_ADD2(ec->rayZero,
                      1.0, ctx->cam->from,
                      ctx->cam->vspNeer, ctx->cam->N);
  }
  return ec;
}

_hooverExtraContext *
_hooverExtraContextNix(_hooverExtraContext *ec) {

  if (ec) {
    free(ec);
  }
  return NULL;
}

/*
** _hooverThreadArg struct
**
** A pointer to this is passed to _hooverThreadBody.  It contains all the
** information which is not thread-specific, and all the thread-specific
** information known at the level of hooverRender.
**
** For simplicity sake, a pointer to a struct of this type is also
** returned from _hooverThreadBody, so this is where we store an
** error-signaling return value (errCode), and what function had
** trouble (whichErr).
*/
typedef struct {
  /* ----------------------- input */
  hooverContext *ctx;
  _hooverExtraContext *ec;
  void *render;
  int whichThread;
  /* ----------------------- output */
  int whichErr;
  int errCode;
} _hooverThreadArg;

void *
_hooverThreadBody(void *_arg) {
  _hooverThreadArg *arg;
  void *thread;
  int ret,               /* to catch return values from callbacks */
    sampleI,             /* which sample we're on */
    inside,              /* we're inside the volume */
    vI, uI;              /* integral coords in image */
  double tmp,
    mm,                  /* lowest position in index space, for all axes */
    Mx, My, Mz,          /* highest position in index space on each axis */
    u, v,                /* floating-point coords in image */
    uvScale,             /* how to scale (u,v) to go from image to
                            near plane, according to ortho or perspective */
    lx, ly, lz,          /* half edge-lengths of volume */
    rayLen=0,            /* length of segment formed by ray line intersecting
                            the near and far clipping planes */
    rayT,                /* current position along ray (world-space) */
    rayDirW[3],          /* unit-length ray direction (world-space) */
    rayDirI[3],          /* rayDirW transformed into index space;
                            not unit length, but a unit change in
                            world space along rayDirW translates to
                            this change in index space along rayDirI */
    rayPosW[3],          /* current ray location (world-space) */
    rayPosI[3],          /* current ray location (index-space) */
    rayStartW[3],        /* ray start on near plane (world-space) */
    rayStartI[3],        /* ray start on near plane (index-space) */
    rayStep,             /* distance between samples (world-space) */
    vOff[3], uOff[3];    /* offsets in arg->ec->wU and arg->ec->wV
                            directions towards start of ray */

  arg = (_hooverThreadArg *)_arg;
  if ( (ret = (arg->ctx->threadBegin)(&thread,
                                      arg->render,
                                      arg->ctx->user,
                                      arg->whichThread)) ) {
    arg->errCode = ret;
    arg->whichErr = hooverErrThreadBegin;
    return arg;
  }
  if (arg->ctx->shape) {
    lx = ly = lz = AIR_NAN;
    if (nrrdCenterNode == arg->ctx->shape->center) {
      mm = 0;
      Mx = arg->ctx->shape->size[0]-1;
      My = arg->ctx->shape->size[1]-1;
      Mz = arg->ctx->shape->size[2]-1;
    } else {
      mm = -0.5;
      Mx = arg->ctx->shape->size[0]-0.5;
      My = arg->ctx->shape->size[1]-0.5;
      Mz = arg->ctx->shape->size[2]-0.5;
    }
  } else {
    lx = arg->ec->volHLen[0];
    ly = arg->ec->volHLen[1];
    lz = arg->ec->volHLen[2];
    if (nrrdCenterNode == arg->ctx->volCentering) {
      mm = 0;
      Mx = arg->ctx->volSize[0]-1;
      My = arg->ctx->volSize[1]-1;
      Mz = arg->ctx->volSize[2]-1;
    } else {
      mm = -0.5;
      Mx = arg->ctx->volSize[0]-0.5;
      My = arg->ctx->volSize[1]-0.5;
      Mz = arg->ctx->volSize[2]-0.5;
    }
  }

  if (arg->ctx->cam->orthographic) {
    ELL_3V_COPY(rayDirW, arg->ctx->cam->N);
    if (arg->ctx->shape) {
      double zeroW[3], zeroI[3];
      ELL_3V_SET(zeroW, 0, 0, 0);
      gageShapeWtoI(arg->ctx->shape, zeroI, zeroW);
      gageShapeWtoI(arg->ctx->shape, rayDirI, rayDirW);
      ELL_3V_SUB(rayDirI, rayDirI, zeroI);
    } else {
      rayDirI[0] = AIR_DELTA(-lx, rayDirW[0], lx, mm, Mx);
      rayDirI[1] = AIR_DELTA(-ly, rayDirW[1], ly, mm, My);
      rayDirI[2] = AIR_DELTA(-lz, rayDirW[2], lz, mm, Mz);
    }
    rayLen = arg->ctx->cam->vspFaar - arg->ctx->cam->vspNeer;
    uvScale = 1.0;
  } else {
    uvScale = arg->ctx->cam->vspNeer/arg->ctx->cam->vspDist;
  }

  while (1) {
    /* the work assignment is simply the next scanline to be rendered:
       the result of all this is setting vI */
    if (arg->ctx->workMutex) {
      airThreadMutexLock(arg->ctx->workMutex);
    }
    vI = arg->ctx->workIdx;
    if (arg->ctx->workIdx < arg->ctx->imgSize[1]) {
      arg->ctx->workIdx += 1;
    }
    if (arg->ctx->workMutex) {
      airThreadMutexUnlock(arg->ctx->workMutex);
    }
    if (vI == arg->ctx->imgSize[1]) {
      /* we're done! */
      break;
    }

    if (nrrdCenterCell == arg->ctx->imgCentering) {
      v = uvScale*AIR_AFFINE(-0.5, vI, arg->ctx->imgSize[1]-0.5,
                             arg->ctx->cam->vRange[0],
                             arg->ctx->cam->vRange[1]);
    } else {
      v = uvScale*AIR_AFFINE(0.0, vI, arg->ctx->imgSize[1]-1.0,
                             arg->ctx->cam->vRange[0],
                             arg->ctx->cam->vRange[1]);
    }
    ELL_3V_SCALE(vOff, v, arg->ctx->cam->V);
    for (uI=0; uI<arg->ctx->imgSize[0]; uI++) {
      if (nrrdCenterCell == arg->ctx->imgCentering) {
        u = uvScale*AIR_AFFINE(-0.5, uI, arg->ctx->imgSize[0]-0.5,
                               arg->ctx->cam->uRange[0],
                               arg->ctx->cam->uRange[1]);
      } else {
        u = uvScale*AIR_AFFINE(0.0, uI, arg->ctx->imgSize[0]-1.0,
                               arg->ctx->cam->uRange[0],
                               arg->ctx->cam->uRange[1]);
      }
      ELL_3V_SCALE(uOff, u, arg->ctx->cam->U);
      ELL_3V_ADD3(rayStartW, uOff, vOff, arg->ec->rayZero);
      if (arg->ctx->shape) {
        gageShapeWtoI(arg->ctx->shape, rayStartI, rayStartW);
      } else {
        rayStartI[0] = AIR_AFFINE(-lx, rayStartW[0], lx, mm, Mx);
        rayStartI[1] = AIR_AFFINE(-ly, rayStartW[1], ly, mm, My);
        rayStartI[2] = AIR_AFFINE(-lz, rayStartW[2], lz, mm, Mz);
      }
      if (!arg->ctx->cam->orthographic) {
        ELL_3V_SUB(rayDirW, rayStartW, arg->ctx->cam->from);
        ELL_3V_NORM(rayDirW, rayDirW, tmp);
        if (arg->ctx->shape) {
          double zeroW[3], zeroI[3];
          ELL_3V_SET(zeroW, 0, 0, 0);
          gageShapeWtoI(arg->ctx->shape, zeroI, zeroW);
          gageShapeWtoI(arg->ctx->shape, rayDirI, rayDirW);
          ELL_3V_SUB(rayDirI, rayDirI, zeroI);
        } else {
          rayDirI[0] = AIR_DELTA(-lx, rayDirW[0], lx, mm, Mx);
          rayDirI[1] = AIR_DELTA(-ly, rayDirW[1], ly, mm, My);
          rayDirI[2] = AIR_DELTA(-lz, rayDirW[2], lz, mm, Mz);
        }
        rayLen = ((arg->ctx->cam->vspFaar - arg->ctx->cam->vspNeer)/
                  ELL_3V_DOT(rayDirW, arg->ctx->cam->N));
      }
      if ( (ret = (arg->ctx->rayBegin)(thread,
                                       arg->render,
                                       arg->ctx->user,
                                       uI, vI, rayLen,
                                       rayStartW, rayStartI,
                                       rayDirW, rayDirI)) ) {
        arg->errCode = ret;
        arg->whichErr = hooverErrRayBegin;
        return arg;
      }

      sampleI = 0;
      rayT = 0;
      while (1) {
        ELL_3V_SCALE_ADD2(rayPosW, 1.0, rayStartW, rayT, rayDirW);
        if (arg->ctx->shape) {
          gageShapeWtoI(arg->ctx->shape, rayPosI, rayPosW);
        } else {
          ELL_3V_SCALE_ADD2(rayPosI, 1.0, rayStartI, rayT, rayDirI);
        }
        inside = (AIR_IN_CL(mm, rayPosI[0], Mx) &&
                  AIR_IN_CL(mm, rayPosI[1], My) &&
                  AIR_IN_CL(mm, rayPosI[2], Mz));
        rayStep = (arg->ctx->sample)(thread,
                                     arg->render,
                                     arg->ctx->user,
                                     sampleI, rayT,
                                     inside,
                                     rayPosW, rayPosI);
        if (!AIR_EXISTS(rayStep)) {
          /* sampling failed */
          arg->errCode = 0;
          arg->whichErr = hooverErrSample;
          return arg;
        }
        if (!rayStep) {
          /* ray decided to finish itself */
          break;
        }
        /* else we moved to a new location along the ray */
        rayT += rayStep;
        if (!AIR_IN_CL(0, rayT, rayLen)) {
          /* ray stepped outside near-far clipping region, its done. */
          break;
        }
        sampleI++;
      }

      if ( (ret = (arg->ctx->rayEnd)(thread,
                                     arg->render,
                                     arg->ctx->user)) ) {
        arg->errCode = ret;
        arg->whichErr = hooverErrRayEnd;
        return arg;
      }
    }  /* end this scanline */
  } /* end while(1) assignment of scanlines */

  if ( (ret = (arg->ctx->threadEnd)(thread,
                                    arg->render,
                                    arg->ctx->user)) ) {
    arg->errCode = ret;
    arg->whichErr = hooverErrThreadEnd;
    return arg;
  }

  /* returning NULL actually indicates that there was NOT an error */
  return NULL;
}

typedef union {
  _hooverThreadArg **h;
  void **v;
} _htpu;

/*
******** hooverRender()
**
** because of the biff usage(), only one thread can call hooverRender(),
** and no promises if the threads themselves call biff...
*/
int
hooverRender(hooverContext *ctx, int *errCodeP, int *errThreadP) {
  static const char me[]="hooverRender";
  _hooverExtraContext *ec;
  _hooverThreadArg args[HOOVER_THREAD_MAX];
  _hooverThreadArg *errArg;
  airThread *thread[HOOVER_THREAD_MAX];
  _htpu u;

  void *render;
  int ret;
  airArray *mop;
  unsigned int threadIdx;

  if (!( errCodeP && errThreadP )) {
    biffAddf(HOOVER, "%s: got NULL int return pointer", me);
    return hooverErrInit;
  }

  /* this calls limnCameraUpdate() */
  if (hooverContextCheck(ctx)) {
    biffAddf(HOOVER, "%s: problem detected in given context", me);
    *errCodeP = 0;
    *errThreadP = 0;
    return hooverErrInit;
  }

  if (!(ec = _hooverExtraContextNew(ctx))) {
    biffAddf(HOOVER, "%s: problem creating thread context", me);
    *errCodeP = 0;
    *errThreadP = 0;
    return hooverErrInit;
  }
  mop = airMopNew();
  airMopAdd(mop, ec, (airMopper)_hooverExtraContextNix, airMopAlways);
  if ( (ret = (ctx->renderBegin)(&render, ctx->user)) ) {
    *errCodeP = ret;
    *errCodeP = 0;
    *errThreadP = 0;
    airMopError(mop);
    return hooverErrRenderBegin;
  }

  for (threadIdx=0; threadIdx<ctx->numThreads; threadIdx++) {
    args[threadIdx].ctx = ctx;
    args[threadIdx].ec = ec;
    args[threadIdx].render = render;
    args[threadIdx].whichThread = threadIdx;
    args[threadIdx].whichErr = hooverErrNone;
    args[threadIdx].errCode = 0;
    thread[threadIdx] = airThreadNew();
  }
  ctx->workIdx = 0;
  if (1 < ctx->numThreads) {
    ctx->workMutex = airThreadMutexNew();
  } else {
    ctx->workMutex = NULL;
  }

  /* (done): call airThreadStart() once per thread, passing the
     address of a distinct (and appropriately intialized)
     _hooverThreadArg to each.  If return of airThreadStart() is
     non-zero, put its return in *errCodeP, the number of the
     problematic in *errThreadP, and return hooverErrThreadCreate.
     Then call airThreadJoin() on all the threads, passing &errArg as
     "retval". On non-zero return, set *errCodeP and *errThreadP,
     and return hooverErrThreadJoin. If return of airThreadJoin() is
     zero, but the errArg is non-NULL, then assume that this errArg
     is actually just the passed _hooverThreadArg returned to us, and
     from this copy errArg->errCode into *errCodeP, and return
     errArg->whichErr */

  if (1 < ctx->numThreads && !airThreadCapable) {
    fprintf(stderr, "%s: WARNING: not multi-threaded; will do %d "
            "\"threads\" serially !!!\n", me, ctx->numThreads);
  }

  for (threadIdx=0; threadIdx<ctx->numThreads; threadIdx++) {
    if ((ret = airThreadStart(thread[threadIdx], _hooverThreadBody,
                              (void *) &args[threadIdx]))) {
      *errCodeP = ret;
      *errThreadP = threadIdx;
      airMopError(mop);
      return hooverErrThreadCreate;
    }
  }

  for (threadIdx=0; threadIdx<ctx->numThreads; threadIdx++) {
    u.h = &errArg;
    if ((ret = airThreadJoin(thread[threadIdx], u.v))) {
      *errCodeP = ret;
      *errThreadP = threadIdx;
      airMopError(mop);
      return hooverErrThreadJoin;
    }
    if (errArg != NULL) {
      *errCodeP = errArg->errCode;
      *errThreadP = threadIdx;
      return errArg->whichErr;
    }
    thread[threadIdx] = airThreadNix(thread[threadIdx]);
  }

  if (1 < ctx->numThreads) {
    ctx->workMutex = airThreadMutexNix(ctx->workMutex);
  }

  if ( (ret = (ctx->renderEnd)(render, ctx->user)) ) {
    *errCodeP = ret;
    *errThreadP = -1;
    return hooverErrRenderEnd;
  }
  render = NULL;
  airMopOkay(mop);

  *errCodeP = 0;
  *errThreadP = 0;
  return hooverErrNone;
}


Generated by: GCOVR (Version 3.3)

Line	Branch	Exec	Source
1			/*
2			Teem: Tools to process and visualize scientific data and images .
3			Copyright (C) 2013, 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
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 *
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 *
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 *
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
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			}