Head

GCC Code Coverage Report

Directory:	./		Exec	Total	Coverage
File:	src/echo/renderEcho.c	Lines:	0	237	0.0 %
Date:	2017-05-26	Branches:	0	111	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 "echo.h"
#include "privateEcho.h"

int
echoThreadStateInit(int threadIdx, echoThreadState *tstate,
                    echoRTParm *parm, echoGlobalState *gstate) {
  static const char me[]="echoThreadStateInit";

  if (!(tstate && parm && gstate)) {
    biffAddf(ECHO, "%s: got NULL pointer", me);
    return 1;
  }
  /* tstate->thread set by echoThreadStateNew() */
  tstate->gstate = gstate;
  /* this will probably be over-written */
  tstate->verbose = gstate->verbose;
  tstate->threadIdx = threadIdx;
  if (nrrdMaybeAlloc_va(tstate->nperm, nrrdTypeInt, 2,
                        AIR_CAST(size_t, ECHO_JITTABLE_NUM),
                        AIR_CAST(size_t, parm->numSamples))) {
    biffMovef(ECHO, NRRD,
              "%s: couldn't allocate jitter permutation array", me);
    return 1;
  }
  nrrdAxisInfoSet_va(tstate->nperm, nrrdAxisInfoLabel,
                     "jittable", "sample");

  if (nrrdMaybeAlloc_va(tstate->njitt, echoPos_nt, 3,
                        AIR_CAST(size_t, 2),
                        AIR_CAST(size_t, ECHO_JITTABLE_NUM),
                        AIR_CAST(size_t, parm->numSamples))) {
    biffMovef(ECHO, NRRD, "%s: couldn't allocate jitter array", me);
    return 1;
  }
  nrrdAxisInfoSet_va(tstate->njitt, nrrdAxisInfoLabel,
                     "x,y", "jittable", "sample");

  tstate->permBuff = AIR_CAST(unsigned int *, airFree(tstate->permBuff));
  if (!(tstate->permBuff = AIR_CAST(unsigned int *,
                                    calloc(parm->numSamples, sizeof(int))))) {
    biffAddf(ECHO, "%s: couldn't allocate permutation buffer", me);
    return 1;
  }
  tstate->chanBuff = (echoCol_t *)airFree(tstate->chanBuff);
  if (!( tstate->chanBuff =
         (echoCol_t*)calloc(ECHO_IMG_CHANNELS * parm->numSamples,
                            sizeof(echoCol_t)) )) {
    biffAddf(ECHO, "%s: couldn't allocate img channel sample buffer", me);
    return 1;
  }

  airSrandMT_r(tstate->rst, AIR_CAST(unsigned int, (parm->seedRand
                                                    ? airTime()
                                                    : threadIdx)));
  tstate->returnPtr = NULL;

  return 0;
}

/*
******** echoJitterCompute()
**
**
*/
void
echoJitterCompute(echoRTParm *parm, echoThreadState *tstate) {
  echoPos_t *jitt, w;
  int s, i, j, xi, yi, n, N, *perm;

  N = parm->numSamples;
  n = (int)sqrt(N);
  w = 1.0/n;
  /* each row in perm[] is for one sample, for going through all jittables;
     each column is a different permutation of [0..parm->numSamples-1] */
  perm = (int *)tstate->nperm->data;
  for (j=0; j<ECHO_JITTABLE_NUM; j++) {
    airShuffle_r(tstate->rst, tstate->permBuff,
                 parm->numSamples, parm->permuteJitter);
    for (s=0; s<N; s++) {
      perm[j + ECHO_JITTABLE_NUM*s] = tstate->permBuff[s];
    }
  }
  jitt = (echoPos_t *)tstate->njitt->data;
  for (s=0; s<N; s++) {
    for (j=0; j<ECHO_JITTABLE_NUM; j++) {
      i = perm[j + ECHO_JITTABLE_NUM*s];
      xi = i % n;
      yi = i / n;
      switch(parm->jitterType) {
      case echoJitterNone:
        jitt[0 + 2*j] = 0.0;
        jitt[1 + 2*j] = 0.0;
        break;
      case echoJitterGrid:
        jitt[0 + 2*j] = NRRD_POS(nrrdCenterCell, -0.5, 0.5, n, xi);
        jitt[1 + 2*j] = NRRD_POS(nrrdCenterCell, -0.5, 0.5, n, yi);
        break;
      case echoJitterJitter:
        jitt[0 + 2*j] = (NRRD_POS(nrrdCenterCell, -0.5, 0.5, n, xi)
                         + w*(airDrandMT_r(tstate->rst) - 0.5));
        jitt[1 + 2*j] = (NRRD_POS(nrrdCenterCell, -0.5, 0.5, n, yi)
                         + w*(airDrandMT_r(tstate->rst) - 0.5));
        break;
      case echoJitterRandom:
        jitt[0 + 2*j] = airDrandMT_r(tstate->rst) - 0.5;
        jitt[1 + 2*j] = airDrandMT_r(tstate->rst) - 0.5;
        break;
      }
    }
    jitt += 2*ECHO_JITTABLE_NUM;
  }

  return;
}

/*
******** echoRTRenderCheck
**
** does all the error checking required of echoRTRender and
** everything that it calls
*/
int
echoRTRenderCheck(Nrrd *nraw, limnCamera *cam, echoScene *scene,
                  echoRTParm *parm, echoGlobalState *gstate) {
  static const char me[]="echoRTRenderCheck";
  int tmp;

  if (!(nraw && cam && scene && parm && gstate)) {
    biffAddf(ECHO, "%s: got NULL pointer", me);
    return 1;
  }
  if (limnCameraUpdate(cam)) {
    biffMovef(ECHO, LIMN, "%s: camera trouble", me);
    return 1;
  }
  if (scene->envmap) {
    if (limnEnvMapCheck(scene->envmap)) {
      biffMovef(ECHO, LIMN, "%s: environment map not valid", me);
      return 1;
    }
  }
  if (airEnumValCheck(echoJitter, parm->jitterType)) {
    biffAddf(ECHO, "%s: jitter method (%d) invalid", me, parm->jitterType);
    return 1;
  }
  if (!(parm->numSamples > 0)) {
    biffAddf(ECHO, "%s: # samples (%d) invalid", me, parm->numSamples);
    return 1;
  }
  if (!(parm->imgResU > 0 && parm->imgResV)) {
    biffAddf(ECHO, "%s: image dimensions (%dx%d) invalid", me,
             parm->imgResU, parm->imgResV);
    return 1;
  }
  if (!AIR_EXISTS(parm->aperture)) {
    biffAddf(ECHO, "%s: aperture doesn't exist", me);
    return 1;
  }

  switch (parm->jitterType) {
  case echoJitterNone:
  case echoJitterRandom:
    break;
  case echoJitterGrid:
  case echoJitterJitter:
    tmp = (int)sqrt(parm->numSamples);
    if (tmp*tmp != parm->numSamples) {
      biffAddf(ECHO,
               "%s: need a square # samples for %s jitter method (not %d)",
               me, airEnumStr(echoJitter, parm->jitterType), parm->numSamples);
      return 1;
    }
    break;
  }

  /* for the time being things are hard-coded to be r,g,b,a,time */
  if (ECHO_IMG_CHANNELS != 5) {
    biffAddf(ECHO, "%s: ECHO_IMG_CHANNELS != 5", me);
    return 1;
  }

  /* all is well */
  return 0;
}

void
echoChannelAverage(echoCol_t *img,
                   echoRTParm *parm, echoThreadState *tstate) {
  int s;
  echoCol_t R, G, B, A, T;

  R = G = B = A = T = 0;
  for (s=0; s<parm->numSamples; s++) {
    R += tstate->chanBuff[0 + ECHO_IMG_CHANNELS*s];
    G += tstate->chanBuff[1 + ECHO_IMG_CHANNELS*s];
    B += tstate->chanBuff[2 + ECHO_IMG_CHANNELS*s];
    A += tstate->chanBuff[3 + ECHO_IMG_CHANNELS*s];
    T += tstate->chanBuff[4 + ECHO_IMG_CHANNELS*s];
  }
  img[0] = R / parm->numSamples;
  img[1] = G / parm->numSamples;
  img[2] = B / parm->numSamples;
  img[3] = A / parm->numSamples;
  img[4] = T;

  return;
}

/*
******** echoRayColor
**
** This is called by echoRTRender and by the various color routines,
** following an intersection with non-phong non-light material (the
** things that require reflection or refraction rays).  As such, it is
** never called on shadow rays.
*/
void
echoRayColor(echoCol_t *chan, echoRay *ray,
             echoScene *scene, echoRTParm *parm, echoThreadState *tstate) {
  static const char me[]="echoRayColor";
  echoCol_t rgba[4];
  echoIntx intx;

  tstate->depth++;
  if (tstate->depth > parm->maxRecDepth) {
    /* we've exceeded the recursion depth, so no more rays for you */
    ELL_4V_SET(chan, parm->maxRecCol[0], parm->maxRecCol[1],
               parm->maxRecCol[2], 1.0);
    goto done;
  }

  intx.boxhits = 0;
  if (!echoRayIntx(&intx, ray, scene, parm, tstate)) {
    if (tstate->verbose) {
      fprintf(stderr, "%s%s: (nothing was hit)\n",_echoDot(tstate->depth), me);
    }
    /* ray hits nothing in scene */
    ELL_4V_SET_TT(chan, echoCol_t,
                  scene->bkgr[0], scene->bkgr[1], scene->bkgr[2],
                  (parm->renderBoxes
                   ? 1.0 - pow(1.0 - parm->boxOpac, intx.boxhits)
                   : 0.0));
    goto done;
  }

  if (tstate->verbose) {
    fprintf(stderr, "%s%s: hit a %d (%p) at (%g,%g,%g)\n"
            "%s    = %g along (%g,%g,%g)\n", _echoDot(tstate->depth), me,
            intx.obj->type, AIR_CAST(void*, intx.obj),
            intx.pos[0], intx.pos[1], intx.pos[2], _echoDot(tstate->depth),
            intx.t, ray->dir[0], ray->dir[1], ray->dir[2]);
  }
  echoIntxColor(rgba, &intx, scene, parm, tstate);
  ELL_4V_COPY(chan, rgba);
 done:
  tstate->depth--;
  return;
}

void *
_echoRTRenderThreadBody(void *_arg) {
  char done[20];
  int imgUi, imgVi,         /* integral pixel indices */
    samp;                   /* which sample are we doing */
  echoPos_t tmp0, tmp1,
    pixUsz, pixVsz,         /* U and V dimensions of a pixel */
    U[4], V[4], N[4],       /* view space basis (only first 3 elements used) */
    imgU, imgV,             /* floating point pixel center locations */
    eye[3],                 /* eye center before jittering */
    at[3],                  /* ray destination (pixel center post-jittering) */
    imgOrig[3];             /* image origin */
  double time0;
  echoRay ray;              /* (not a pointer) */
  echoThreadState *arg;
  echoCol_t *img, *chan;    /* current scanline of channel buffer array */
  Nrrd *nraw;               /* copies of arguments to echoRTRender . . . */
  limnCamera *cam;
  echoScene *scene;
  echoRTParm *parm;

  arg = (echoThreadState *)_arg;
  nraw = arg->gstate->nraw;
  cam = arg->gstate->cam;
  scene = arg->gstate->scene;
  parm = arg->gstate->parm;

  echoJitterCompute(arg->gstate->parm, arg);
  if (arg->gstate->verbose > 2) {
    nrrdSave("jitt.nrrd", arg->njitt, NULL);
  }

  /* set eye, U, V, N, imgOrig */
  ELL_3V_COPY(eye, arg->gstate->cam->from);
  ELL_4MV_ROW0_GET(U, cam->W2V);
  ELL_4MV_ROW1_GET(V, cam->W2V);
  ELL_4MV_ROW2_GET(N, cam->W2V);
  ELL_3V_SCALE_ADD2(imgOrig, 1.0, eye, cam->vspDist, N);

  /* determine size of a single pixel (based on cell-centering) */
  pixUsz = (cam->uRange[1] - cam->uRange[0])/(parm->imgResU);
  pixVsz = (cam->vRange[1] - cam->vRange[0])/(parm->imgResV);

  arg->depth = 0;
  ray.shadow = AIR_FALSE;
  arg->verbose = AIR_FALSE;

  while (1) {
    if (arg->gstate->workMutex) {
      airThreadMutexLock(arg->gstate->workMutex);
    }
    imgVi = arg->gstate->workIdx;
    if (arg->gstate->workIdx < parm->imgResV) {
      arg->gstate->workIdx += 1;
    }
    if (!(imgVi % 5)) {
      fprintf(stderr, "%s", airDoneStr(0, imgVi, parm->imgResV-1, done));
      fflush(stderr);
    }
    if (arg->gstate->workMutex) {
      airThreadMutexUnlock(arg->gstate->workMutex);
    }
    if (imgVi == parm->imgResV) {
      /* we're done! */
      break;
    }

    imgV = NRRD_POS(nrrdCenterCell, cam->vRange[0], cam->vRange[1],
                    parm->imgResV, imgVi);
    for (imgUi=0; imgUi<parm->imgResU; imgUi++) {
      imgU = NRRD_POS(nrrdCenterCell, cam->uRange[0], cam->uRange[1],
                      parm->imgResU, imgUi);
      img = ((echoCol_t *)nraw->data
             + ECHO_IMG_CHANNELS*(imgUi + parm->imgResU*imgVi));

      /* initialize things on first "scanline" */
      arg->jitt = (echoPos_t *)arg->njitt->data;
      chan = arg->chanBuff;

      /*
      arg->verbose = ( (48 == imgUi && 13 == imgVi)
                       || (49 == imgUi && 13 == imgVi) );
      */

      if (arg->verbose) {
        fprintf(stderr, "\n");
        fprintf(stderr, "-----------------------------------------------\n");
        fprintf(stderr, "----------------- (%3d, %3d) ------------------\n",
                imgUi, imgVi);
        fprintf(stderr, "-----------------------------------------------\n\n");
      }

      /* go through samples */
      for (samp=0; samp<parm->numSamples; samp++) {
        /* set ray.from[] */
        ELL_3V_COPY(ray.from, eye);
        if (parm->aperture) {
          tmp0 = parm->aperture*(arg->jitt[0 + 2*echoJittableLens]);
          tmp1 = parm->aperture*(arg->jitt[1 + 2*echoJittableLens]);
          ELL_3V_SCALE_ADD3(ray.from, 1, ray.from, tmp0, U, tmp1, V);
        }

        /* set at[] */
        tmp0 = imgU + pixUsz*(arg->jitt[0 + 2*echoJittablePixel]);
        tmp1 = imgV + pixVsz*(arg->jitt[1 + 2*echoJittablePixel]);
        ELL_3V_SCALE_ADD3(at, 1, imgOrig, tmp0, U, tmp1, V);

        /* do it! */
        ELL_3V_SUB(ray.dir, at, ray.from);
        ELL_3V_NORM(ray.dir, ray.dir, tmp0);
        ray.neer = 0.0;
        ray.faar = ECHO_POS_MAX;
        time0 = airTime();
        if (0) {
          memset(chan, 0, ECHO_IMG_CHANNELS*sizeof(echoCol_t));
        } else {
          echoRayColor(chan, &ray, scene, parm, arg);
        }
        chan[4] = AIR_CAST(echoCol_t, airTime() - time0);

        /* move to next "scanline" */
        arg->jitt += 2*ECHO_JITTABLE_NUM;
        chan += ECHO_IMG_CHANNELS;
      }
      echoChannelAverage(img, parm, arg);
      img += ECHO_IMG_CHANNELS;
      if (!parm->reuseJitter) {
        echoJitterCompute(parm, arg);
      }
    }
  }

  return _arg;
}


/*
******** echoRTRender
**
** top-level call to accomplish all (ray-tracing) rendering.  As much
** error checking as possible should be done here and not in the
** lower-level functions.
*/
int
echoRTRender(Nrrd *nraw, limnCamera *cam, echoScene *scene,
             echoRTParm *parm, echoGlobalState *gstate) {
  static const char me[]="echoRTRender";
  int tid, ret;
  airArray *mop;
  echoThreadState *tstate[ECHO_THREAD_MAX];

  if (echoRTRenderCheck(nraw, cam, scene, parm, gstate)) {
    biffAddf(ECHO, "%s: problem with input", me);
    return 1;
  }
  gstate->nraw = nraw;
  gstate->cam = cam;
  gstate->scene = scene;
  gstate->parm = parm;
  mop = airMopNew();
  if (nrrdMaybeAlloc_va(nraw, echoCol_nt, 3,
                        AIR_CAST(size_t, ECHO_IMG_CHANNELS),
                        AIR_CAST(size_t, parm->imgResU),
                        AIR_CAST(size_t, parm->imgResV))) {
    biffMovef(ECHO, NRRD, "%s: couldn't allocate output image", me);
    airMopError(mop); return 1;
  }
  airMopAdd(mop, nraw, (airMopper)nrrdNix, airMopOnError);
  nrrdAxisInfoSet_va(nraw, nrrdAxisInfoLabel,
                     "r,g,b,a,t", "x", "y");
  nrrdAxisInfoSet_va(nraw, nrrdAxisInfoMin,
                     AIR_NAN, cam->uRange[0], cam->vRange[0]);
  nrrdAxisInfoSet_va(nraw, nrrdAxisInfoMax,
                     AIR_NAN, cam->uRange[1], cam->vRange[1]);
  gstate->time = airTime();

  if (parm->numThreads > 1) {
    gstate->workMutex = airThreadMutexNew();
    airMopAdd(mop, gstate->workMutex,
              (airMopper)airThreadMutexNix, airMopAlways);
  } else {
    gstate->workMutex = NULL;
  }
  for (tid=0; tid<parm->numThreads; tid++) {
    if (!( tstate[tid] = echoThreadStateNew() )) {
      biffAddf(ECHO, "%s: failed to create thread state %d", me, tid);
      airMopError(mop); return 1;
    }
    if (echoThreadStateInit(tid, tstate[tid], parm, gstate)) {
      biffAddf(ECHO, "%s: failed to initialized thread state %d", me, tid);
      airMopError(mop); return 1;
    }
    airMopAdd(mop, tstate[tid], (airMopper)echoThreadStateNix, airMopAlways);
  }
  fprintf(stderr, "%s:       ", me);  /* prep for printing airDoneStr */
  gstate->workIdx = 0;
  for (tid=0; tid<parm->numThreads; tid++) {
    if (( ret = airThreadStart(tstate[tid]->thread, _echoRTRenderThreadBody,
                               (void *)(tstate[tid])) )) {
      biffAddf(ECHO, "%s: thread[%d] failed to start: %d", me, tid, ret);
      airMopError(mop); return 1;
    }
  }
  for (tid=0; tid<parm->numThreads; tid++) {
    if (( ret = airThreadJoin(tstate[tid]->thread,
                              (void **)(&(tstate[tid]->returnPtr))) )) {
      biffAddf(ECHO, "%s: thread[%d] failed to join: %d", me, tid, ret);
      airMopError(mop); return 1;
    }
  }

  gstate->time = airTime() - gstate->time;
  fprintf(stderr, "\n%s: time = %g\n", me, gstate->time);

  airMopOkay(mop);
  return 0;
}


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 "echo.h"
25			#include "privateEcho.h"
26
27			int
28			echoThreadStateInit(int threadIdx, echoThreadState *tstate,
29			echoRTParm parm, echoGlobalState gstate) {
30			static const char me[]="echoThreadStateInit";
31
32			if (!(tstate && parm && gstate)) {
33			biffAddf(ECHO, "%s: got NULL pointer", me);
34			return 1;
35			}
36			/* tstate->thread set by echoThreadStateNew() */
37			tstate->gstate = gstate;
38			/* this will probably be over-written */
39			tstate->verbose = gstate->verbose;
40			tstate->threadIdx = threadIdx;
41			if (nrrdMaybeAlloc_va(tstate->nperm, nrrdTypeInt, 2,
42			AIR_CAST(size_t, ECHO_JITTABLE_NUM),
43			AIR_CAST(size_t, parm->numSamples))) {
44			biffMovef(ECHO, NRRD,
45			"%s: couldn't allocate jitter permutation array", me);
46			return 1;
47			}
48			nrrdAxisInfoSet_va(tstate->nperm, nrrdAxisInfoLabel,
49			"jittable", "sample");
50
51			if (nrrdMaybeAlloc_va(tstate->njitt, echoPos_nt, 3,
52			AIR_CAST(size_t, 2),
53			AIR_CAST(size_t, ECHO_JITTABLE_NUM),
54			AIR_CAST(size_t, parm->numSamples))) {
55			biffMovef(ECHO, NRRD, "%s: couldn't allocate jitter array", me);
56			return 1;
57			}
58			nrrdAxisInfoSet_va(tstate->njitt, nrrdAxisInfoLabel,
59			"x,y", "jittable", "sample");
60
61			tstate->permBuff = AIR_CAST(unsigned int *, airFree(tstate->permBuff));
62			if (!(tstate->permBuff = AIR_CAST(unsigned int *,
63			calloc(parm->numSamples, sizeof(int))))) {
64			biffAddf(ECHO, "%s: couldn't allocate permutation buffer", me);
65			return 1;
66			}
67			tstate->chanBuff = (echoCol_t *)airFree(tstate->chanBuff);
68			if (!( tstate->chanBuff =
69			(echoCol_t)calloc(ECHO_IMG_CHANNELS parm->numSamples,
70			sizeof(echoCol_t)) )) {
71			biffAddf(ECHO, "%s: couldn't allocate img channel sample buffer", me);
72			return 1;
73			}
74
75			airSrandMT_r(tstate->rst, AIR_CAST(unsigned int, (parm->seedRand
76			? airTime()
77			: threadIdx)));
78			tstate->returnPtr = NULL;
79
80			return 0;
81			}
82
83			/*
84			******** echoJitterCompute()
85			**
86			**
87			*/
88			void
89			echoJitterCompute(echoRTParm parm, echoThreadState tstate) {
90			echoPos_t *jitt, w;
91			int s, i, j, xi, yi, n, N, *perm;
92
93			N = parm->numSamples;
94			n = (int)sqrt(N);
95			w = 1.0/n;
96			/* each row in perm[] is for one sample, for going through all jittables;
97			each column is a different permutation of [0..parm->numSamples-1] */
98			perm = (int *)tstate->nperm->data;
99			for (j=0; j<ECHO_JITTABLE_NUM; j++) {
100			airShuffle_r(tstate->rst, tstate->permBuff,
101			parm->numSamples, parm->permuteJitter);
102			for (s=0; s<N; s++) {
103			perm[j + ECHO_JITTABLE_NUM*s] = tstate->permBuff[s];
104			}
105			}
106			jitt = (echoPos_t *)tstate->njitt->data;
107			for (s=0; s<N; s++) {
108			for (j=0; j<ECHO_JITTABLE_NUM; j++) {
109			i = perm[j + ECHO_JITTABLE_NUM*s];
110			xi = i % n;
111			yi = i / n;
112			switch(parm->jitterType) {
113			case echoJitterNone:
114			jitt[0 + 2*j] = 0.0;
115			jitt[1 + 2*j] = 0.0;
116			break;
117			case echoJitterGrid:
118			jitt[0 + 2*j] = NRRD_POS(nrrdCenterCell, -0.5, 0.5, n, xi);
119			jitt[1 + 2*j] = NRRD_POS(nrrdCenterCell, -0.5, 0.5, n, yi);
120			break;
121			case echoJitterJitter:
122			jitt[0 + 2*j] = (NRRD_POS(nrrdCenterCell, -0.5, 0.5, n, xi)
123			+ w*(airDrandMT_r(tstate->rst) - 0.5));
124			jitt[1 + 2*j] = (NRRD_POS(nrrdCenterCell, -0.5, 0.5, n, yi)
125			+ w*(airDrandMT_r(tstate->rst) - 0.5));
126			break;
127			case echoJitterRandom:
128			jitt[0 + 2*j] = airDrandMT_r(tstate->rst) - 0.5;
129			jitt[1 + 2*j] = airDrandMT_r(tstate->rst) - 0.5;
130			break;
131			}
132			}
133			jitt += 2*ECHO_JITTABLE_NUM;
134			}
135
136			return;
137			}
138
139			/*
140			******** echoRTRenderCheck
141			**
142			** does all the error checking required of echoRTRender and
143			** everything that it calls
144			*/
145			int
146			echoRTRenderCheck(Nrrd nraw, limnCamera cam, echoScene *scene,
147			echoRTParm parm, echoGlobalState gstate) {
148			static const char me[]="echoRTRenderCheck";
149			int tmp;
150
151			if (!(nraw && cam && scene && parm && gstate)) {
152			biffAddf(ECHO, "%s: got NULL pointer", me);
153			return 1;
154			}
155			if (limnCameraUpdate(cam)) {
156			biffMovef(ECHO, LIMN, "%s: camera trouble", me);
157			return 1;
158			}
159			if (scene->envmap) {
160			if (limnEnvMapCheck(scene->envmap)) {
161			biffMovef(ECHO, LIMN, "%s: environment map not valid", me);
162			return 1;
163			}
164			}
165			if (airEnumValCheck(echoJitter, parm->jitterType)) {
166			biffAddf(ECHO, "%s: jitter method (%d) invalid", me, parm->jitterType);
167			return 1;
168			}
169			if (!(parm->numSamples > 0)) {
170			biffAddf(ECHO, "%s: # samples (%d) invalid", me, parm->numSamples);
171			return 1;
172			}
173			if (!(parm->imgResU > 0 && parm->imgResV)) {
174			biffAddf(ECHO, "%s: image dimensions (%dx%d) invalid", me,
175			parm->imgResU, parm->imgResV);
176			return 1;
177			}
178			if (!AIR_EXISTS(parm->aperture)) {
179			biffAddf(ECHO, "%s: aperture doesn't exist", me);
180			return 1;
181			}
182
183			switch (parm->jitterType) {
184			case echoJitterNone:
185			case echoJitterRandom:
186			break;
187			case echoJitterGrid:
188			case echoJitterJitter:
189			tmp = (int)sqrt(parm->numSamples);
190			if (tmp*tmp != parm->numSamples) {
191			biffAddf(ECHO,
192			"%s: need a square # samples for %s jitter method (not %d)",
193			me, airEnumStr(echoJitter, parm->jitterType), parm->numSamples);
194			return 1;
195			}
196			break;
197			}
198
199			/* for the time being things are hard-coded to be r,g,b,a,time */
200			if (ECHO_IMG_CHANNELS != 5) {
201			biffAddf(ECHO, "%s: ECHO_IMG_CHANNELS != 5", me);
202			return 1;
203			}
204
205			/* all is well */
206			return 0;
207			}
208
209			void
210			echoChannelAverage(echoCol_t *img,
211			echoRTParm parm, echoThreadState tstate) {
212			int s;
213			echoCol_t R, G, B, A, T;
214
215			R = G = B = A = T = 0;
216			for (s=0; s<parm->numSamples; s++) {
217			R += tstate->chanBuff[0 + ECHO_IMG_CHANNELS*s];
218			G += tstate->chanBuff[1 + ECHO_IMG_CHANNELS*s];
219			B += tstate->chanBuff[2 + ECHO_IMG_CHANNELS*s];
220			A += tstate->chanBuff[3 + ECHO_IMG_CHANNELS*s];
221			T += tstate->chanBuff[4 + ECHO_IMG_CHANNELS*s];
222			}
223			img[0] = R / parm->numSamples;
224			img[1] = G / parm->numSamples;
225			img[2] = B / parm->numSamples;
226			img[3] = A / parm->numSamples;
227			img[4] = T;
228
229			return;
230			}
231
232			/*
233			******** echoRayColor
234			**
235			** This is called by echoRTRender and by the various color routines,
236			** following an intersection with non-phong non-light material (the
237			** things that require reflection or refraction rays). As such, it is
238			** never called on shadow rays.
239			*/
240			void
241			echoRayColor(echoCol_t chan, echoRay ray,
242			echoScene scene, echoRTParm parm, echoThreadState *tstate) {
243			static const char me[]="echoRayColor";
244			echoCol_t rgba[4];
245			echoIntx intx;
246
247			tstate->depth++;
248			if (tstate->depth > parm->maxRecDepth) {
249			/* we've exceeded the recursion depth, so no more rays for you */
250			ELL_4V_SET(chan, parm->maxRecCol[0], parm->maxRecCol[1],
251			parm->maxRecCol[2], 1.0);
252			goto done;
253			}
254
255			intx.boxhits = 0;
256			if (!echoRayIntx(&intx, ray, scene, parm, tstate)) {
257			if (tstate->verbose) {
258			fprintf(stderr, "%s%s: (nothing was hit)\n",_echoDot(tstate->depth), me);
259			}
260			/* ray hits nothing in scene */
261			ELL_4V_SET_TT(chan, echoCol_t,
262			scene->bkgr[0], scene->bkgr[1], scene->bkgr[2],
263			(parm->renderBoxes
264			? 1.0 - pow(1.0 - parm->boxOpac, intx.boxhits)
265			: 0.0));
266			goto done;
267			}
268
269			if (tstate->verbose) {
270			fprintf(stderr, "%s%s: hit a %d (%p) at (%g,%g,%g)\n"
271			"%s = %g along (%g,%g,%g)\n", _echoDot(tstate->depth), me,
272			intx.obj->type, AIR_CAST(void*, intx.obj),
273			intx.pos[0], intx.pos[1], intx.pos[2], _echoDot(tstate->depth),
274			intx.t, ray->dir[0], ray->dir[1], ray->dir[2]);
275			}
276			echoIntxColor(rgba, &intx, scene, parm, tstate);
277			ELL_4V_COPY(chan, rgba);
278			done:
279			tstate->depth--;
280			return;
281			}
282
283			void *
284			_echoRTRenderThreadBody(void *_arg) {
285			char done[20];
286			int imgUi, imgVi, /* integral pixel indices */
287			samp; /* which sample are we doing */
288			echoPos_t tmp0, tmp1,
289			pixUsz, pixVsz, /* U and V dimensions of a pixel */
290			U[4], V[4], N[4], /* view space basis (only first 3 elements used) */
291			imgU, imgV, /* floating point pixel center locations */
292			eye[3], /* eye center before jittering */
293			at[3], /* ray destination (pixel center post-jittering) */
294			imgOrig[3]; /* image origin */
295			double time0;
296			echoRay ray; /* (not a pointer) */
297			echoThreadState *arg;
298			echoCol_t img, chan; /* current scanline of channel buffer array */
299			Nrrd nraw; / copies of arguments to echoRTRender . . . */
300			limnCamera *cam;
301			echoScene *scene;
302			echoRTParm *parm;
303
304			arg = (echoThreadState *)_arg;
305			nraw = arg->gstate->nraw;
306			cam = arg->gstate->cam;
307			scene = arg->gstate->scene;
308			parm = arg->gstate->parm;
309
310			echoJitterCompute(arg->gstate->parm, arg);
311			if (arg->gstate->verbose > 2) {
312			nrrdSave("jitt.nrrd", arg->njitt, NULL);
313			}
314
315			/* set eye, U, V, N, imgOrig */
316			ELL_3V_COPY(eye, arg->gstate->cam->from);
317			ELL_4MV_ROW0_GET(U, cam->W2V);
318			ELL_4MV_ROW1_GET(V, cam->W2V);
319			ELL_4MV_ROW2_GET(N, cam->W2V);
320			ELL_3V_SCALE_ADD2(imgOrig, 1.0, eye, cam->vspDist, N);
321
322			/* determine size of a single pixel (based on cell-centering) */
323			pixUsz = (cam->uRange[1] - cam->uRange[0])/(parm->imgResU);
324			pixVsz = (cam->vRange[1] - cam->vRange[0])/(parm->imgResV);
325
326			arg->depth = 0;
327			ray.shadow = AIR_FALSE;
328			arg->verbose = AIR_FALSE;
329
330			while (1) {
331			if (arg->gstate->workMutex) {
332			airThreadMutexLock(arg->gstate->workMutex);
333			}
334			imgVi = arg->gstate->workIdx;
335			if (arg->gstate->workIdx < parm->imgResV) {
336			arg->gstate->workIdx += 1;
337			}
338			if (!(imgVi % 5)) {
339			fprintf(stderr, "%s", airDoneStr(0, imgVi, parm->imgResV-1, done));
340			fflush(stderr);
341			}
342			if (arg->gstate->workMutex) {
343			airThreadMutexUnlock(arg->gstate->workMutex);
344			}
345			if (imgVi == parm->imgResV) {
346			/* we're done! */
347			break;
348			}
349
350			imgV = NRRD_POS(nrrdCenterCell, cam->vRange[0], cam->vRange[1],
351			parm->imgResV, imgVi);
352			for (imgUi=0; imgUi<parm->imgResU; imgUi++) {
353			imgU = NRRD_POS(nrrdCenterCell, cam->uRange[0], cam->uRange[1],
354			parm->imgResU, imgUi);
355			img = ((echoCol_t *)nraw->data
356			+ ECHO_IMG_CHANNELS(imgUi + parm->imgResUimgVi));
357
358			/* initialize things on first "scanline" */
359			arg->jitt = (echoPos_t *)arg->njitt->data;
360			chan = arg->chanBuff;
361
362			/*
363			arg->verbose = ( (48 == imgUi && 13 == imgVi)
364			\|\| (49 == imgUi && 13 == imgVi) );
365			*/
366
367			if (arg->verbose) {
368			fprintf(stderr, "\n");
369			fprintf(stderr, "-----------------------------------------------\n");
370			fprintf(stderr, "----------------- (%3d, %3d) ------------------\n",
371			imgUi, imgVi);
372			fprintf(stderr, "-----------------------------------------------\n\n");
373			}
374
375			/* go through samples */
376			for (samp=0; samp<parm->numSamples; samp++) {
377			/* set ray.from[] */
378			ELL_3V_COPY(ray.from, eye);
379			if (parm->aperture) {
380			tmp0 = parm->aperture(arg->jitt[0 + 2echoJittableLens]);
381			tmp1 = parm->aperture(arg->jitt[1 + 2echoJittableLens]);
382			ELL_3V_SCALE_ADD3(ray.from, 1, ray.from, tmp0, U, tmp1, V);
383			}
384
385			/* set at[] */
386			tmp0 = imgU + pixUsz(arg->jitt[0 + 2echoJittablePixel]);
387			tmp1 = imgV + pixVsz(arg->jitt[1 + 2echoJittablePixel]);
388			ELL_3V_SCALE_ADD3(at, 1, imgOrig, tmp0, U, tmp1, V);
389
390			/* do it! */
391			ELL_3V_SUB(ray.dir, at, ray.from);
392			ELL_3V_NORM(ray.dir, ray.dir, tmp0);
393			ray.neer = 0.0;
394			ray.faar = ECHO_POS_MAX;
395			time0 = airTime();
396			if (0) {
397			memset(chan, 0, ECHO_IMG_CHANNELS*sizeof(echoCol_t));
398			} else {
399			echoRayColor(chan, &ray, scene, parm, arg);
400			}
401			chan[4] = AIR_CAST(echoCol_t, airTime() - time0);
402
403			/* move to next "scanline" */
404			arg->jitt += 2*ECHO_JITTABLE_NUM;
405			chan += ECHO_IMG_CHANNELS;
406			}
407			echoChannelAverage(img, parm, arg);
408			img += ECHO_IMG_CHANNELS;
409			if (!parm->reuseJitter) {
410			echoJitterCompute(parm, arg);
411			}
412			}
413			}
414
415			return _arg;
416			}
417
418
419			/*
420			******** echoRTRender
421			**
422			** top-level call to accomplish all (ray-tracing) rendering. As much
423			** error checking as possible should be done here and not in the
424			** lower-level functions.
425			*/
426			int
427			echoRTRender(Nrrd nraw, limnCamera cam, echoScene *scene,
428			echoRTParm parm, echoGlobalState gstate) {
429			static const char me[]="echoRTRender";
430			int tid, ret;
431			airArray *mop;
432			echoThreadState *tstate[ECHO_THREAD_MAX];
433
434			if (echoRTRenderCheck(nraw, cam, scene, parm, gstate)) {
435			biffAddf(ECHO, "%s: problem with input", me);
436			return 1;
437			}
438			gstate->nraw = nraw;
439			gstate->cam = cam;
440			gstate->scene = scene;
441			gstate->parm = parm;
442			mop = airMopNew();
443			if (nrrdMaybeAlloc_va(nraw, echoCol_nt, 3,
444			AIR_CAST(size_t, ECHO_IMG_CHANNELS),
445			AIR_CAST(size_t, parm->imgResU),
446			AIR_CAST(size_t, parm->imgResV))) {
447			biffMovef(ECHO, NRRD, "%s: couldn't allocate output image", me);
448			airMopError(mop); return 1;
449			}
450			airMopAdd(mop, nraw, (airMopper)nrrdNix, airMopOnError);
451			nrrdAxisInfoSet_va(nraw, nrrdAxisInfoLabel,
452			"r,g,b,a,t", "x", "y");
453			nrrdAxisInfoSet_va(nraw, nrrdAxisInfoMin,
454			AIR_NAN, cam->uRange[0], cam->vRange[0]);
455			nrrdAxisInfoSet_va(nraw, nrrdAxisInfoMax,
456			AIR_NAN, cam->uRange[1], cam->vRange[1]);
457			gstate->time = airTime();
458
459			if (parm->numThreads > 1) {
460			gstate->workMutex = airThreadMutexNew();
461			airMopAdd(mop, gstate->workMutex,
462			(airMopper)airThreadMutexNix, airMopAlways);
463			} else {
464			gstate->workMutex = NULL;
465			}
466			for (tid=0; tid<parm->numThreads; tid++) {
467			if (!( tstate[tid] = echoThreadStateNew() )) {
468			biffAddf(ECHO, "%s: failed to create thread state %d", me, tid);
469			airMopError(mop); return 1;
470			}
471			if (echoThreadStateInit(tid, tstate[tid], parm, gstate)) {
472			biffAddf(ECHO, "%s: failed to initialized thread state %d", me, tid);
473			airMopError(mop); return 1;
474			}
475			airMopAdd(mop, tstate[tid], (airMopper)echoThreadStateNix, airMopAlways);
476			}
477			fprintf(stderr, "%s: ", me); /* prep for printing airDoneStr */
478			gstate->workIdx = 0;
479			for (tid=0; tid<parm->numThreads; tid++) {
480			if (( ret = airThreadStart(tstate[tid]->thread, _echoRTRenderThreadBody,
481			(void *)(tstate[tid])) )) {
482			biffAddf(ECHO, "%s: thread[%d] failed to start: %d", me, tid, ret);
483			airMopError(mop); return 1;
484			}
485			}
486			for (tid=0; tid<parm->numThreads; tid++) {
487			if (( ret = airThreadJoin(tstate[tid]->thread,
488			(void **)(&(tstate[tid]->returnPtr))) )) {
489			biffAddf(ECHO, "%s: thread[%d] failed to join: %d", me, tid, ret);
490			airMopError(mop); return 1;
491			}
492			}
493
494			gstate->time = airTime() - gstate->time;
495			fprintf(stderr, "\n%s: time = %g\n", me, gstate->time);
496
497			airMopOkay(mop);
498			return 0;
499			}