Head

GCC Code Coverage Report

Directory:	./		Exec	Total	Coverage
File:	src/gage/update.c	Lines:	162	234	69.2 %
Date:	2017-05-26	Branches:	103	200	51.5 %


/*
  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 "gage.h"
#include "privateGage.h"

int
_gagePvlFlagCheck(gageContext *ctx, int pvlFlag) {
  int ret;
  unsigned int pvlIdx;

  ret = AIR_FALSE;
  for (pvlIdx=0; pvlIdx<ctx->pvlNum; pvlIdx++) {
    ret |= ctx->pvl[pvlIdx]->flag[pvlFlag];
  }
  return ret;
}

void
_gagePvlFlagDown(gageContext *ctx, int pvlFlag) {
  unsigned int pvlIdx;

  for (pvlIdx=0; pvlIdx<ctx->pvlNum; pvlIdx++) {
    ctx->pvl[pvlIdx]->flag[pvlFlag] = AIR_FALSE;
  }
}

/*
** One could go from all the pvls' queries to the context's needD in
** one shot, but doing it in two steps (as below) seems a little clearer,
** and it means that pvl->needD isn't needlessly re-computed for
** pvl's whose query hasn't changed.
*/

/*
** for each pvl: pvl's query --> pvl's needD
*/
void
_gagePvlNeedDUpdate(gageContext *ctx) {
  static const char me[]="_gagePvlNeedDUpdate";
  gagePerVolume *pvl;
  int que, needD[GAGE_DERIV_MAX+1];
  unsigned int pvlIdx, di;

  if (ctx->verbose) fprintf(stderr, "%s: hello\n", me);
  for (pvlIdx=0; pvlIdx<ctx->pvlNum; pvlIdx++) {
    pvl = ctx->pvl[pvlIdx];
    if (pvl->flag[gagePvlFlagQuery]) {
      GAGE_DV_SET(needD, 0, 0, 0);
      que = pvl->kind->itemMax+1;
      do {
        que--;
        if (GAGE_QUERY_ITEM_TEST(pvl->query, que)) {
          needD[pvl->kind->table[que].needDeriv] = 1;
        }
      } while (que);
      if (!GAGE_DV_EQUAL(needD, pvl->needD)) {
        if (ctx->verbose) {
          fprintf(stderr, "%s: updating pvl[%d]'s needD to (", me, pvlIdx);
          for (di=0; di<=GAGE_DERIV_MAX; di++) {
            fprintf(stderr, "%s%d", di ? "," : "", needD[di]);
          }
          fprintf(stderr, "\n");
        }
        GAGE_DV_COPY(pvl->needD, needD);
        pvl->flag[gagePvlFlagNeedD] = AIR_TRUE;
      }
    }
  }
  if (ctx->verbose) fprintf(stderr, "%s: bye\n", me);
  return;
}

/*
** all pvls' needD --> ctx's needD
*/
void
_gageNeedDUpdate(gageContext *ctx) {
  static const char me[]="_gageNeedDUpdate";
  gagePerVolume *pvl;
  int needD[GAGE_DERIV_MAX+1];
  unsigned int pvlIdx, di;

  if (ctx->verbose) fprintf(stderr, "%s: hello\n", me);
  GAGE_DV_SET(needD, 0, 0, 0);
  for (pvlIdx=0; pvlIdx<ctx->pvlNum; pvlIdx++) {
    pvl = ctx->pvl[pvlIdx];
    for (di=0; di<=GAGE_DERIV_MAX; di++) {
      needD[di] |= pvl->needD[di];
    }
  }
  if (!GAGE_DV_EQUAL(needD, ctx->needD)) {
    if (ctx->verbose) {
      fprintf(stderr, "%s: updating ctx's needD to (", me);
      for (di=0; di<=GAGE_DERIV_MAX; di++) {
        fprintf(stderr, "%s%d", di ? "," : "", needD[di]);
      }
      fprintf(stderr, "\n");
    }
    GAGE_DV_COPY(ctx->needD, needD);
    ctx->flag[gageCtxFlagNeedD] = AIR_TRUE;
  }
  if (ctx->verbose) fprintf(stderr, "%s: bye\n", me);

  return;
}

/*
** ctx's needD & k3pack --> needK
*/
static int
_gageNeedKUpdate(gageContext *ctx) {
  static const char me[]="_gageNeedKUpdate";
  int kernIdx, needK[GAGE_KERNEL_MAX+1], change;
  unsigned int di;

  if (ctx->verbose) fprintf(stderr, "%s: hello\n", me);
  for (kernIdx=gageKernelUnknown+1; kernIdx<gageKernelLast; kernIdx++) {
    needK[kernIdx] = AIR_FALSE;
  }
  if (!ctx->parm.k3pack) {
    biffAddf(GAGE, "%s: sorry, only 3-pack filtering implemented now", me);
    return 1;
  }
  di = 0;
  if (ctx->needD[di]) {
    needK[gageKernel00] = AIR_TRUE;
  }
  di += 1;
  if (ctx->needD[di]) {
    needK[gageKernel00] = needK[gageKernel11] = AIR_TRUE;
  }
  di += 1;
  if (ctx->needD[di]) {
    needK[gageKernel00] = needK[gageKernel11] =
      needK[gageKernel22] = AIR_TRUE;
  }
  if (GAGE_DERIV_MAX != di) {
    biffAddf(GAGE, "%s: sorry, code needs updating for GAGE_DERIV_MAX %u",
             me, GAGE_DERIV_MAX);
    return 1;
  }
  change = AIR_FALSE;
  for (kernIdx=gageKernelUnknown+1; kernIdx<gageKernelLast; kernIdx++) {
    change |= (needK[kernIdx] != ctx->needK[kernIdx]);
  }
  if (change) {
    if (ctx->verbose) {
      fprintf(stderr, "%s: changing needK to (", me);
      for (kernIdx=gageKernelUnknown+1; kernIdx<gageKernelLast; kernIdx++) {
        fprintf(stderr, "%s%d", kernIdx > gageKernelUnknown+1 ? "," : "",
                needK[kernIdx]);
      }
      fprintf(stderr, ")\n");
    }
    for (kernIdx=gageKernelUnknown+1; kernIdx<gageKernelLast; kernIdx++) {
      ctx->needK[kernIdx] = needK[kernIdx];
    }
    ctx->flag[gageCtxFlagNeedK] = AIR_TRUE;
  }
  if (ctx->verbose) fprintf(stderr, "%s: bye\n", me);

  return 0;
}

/*
** ctx's ksp[] & needK --> radius
**
*/
int
_gageRadiusUpdate(gageContext *ctx) {
  static const char me[]="_gageRadiusUpdate";
  unsigned int kernIdx, radius;
  double maxRad, rad;
  NrrdKernelSpec *ksp;

  if (ctx->verbose) fprintf(stderr, "%s: hello\n", me);
  maxRad = 0;
  for (kernIdx=gageKernelUnknown+1; kernIdx<gageKernelLast; kernIdx++) {
    if (ctx->needK[kernIdx]) {
      ksp = ctx->ksp[kernIdx];
      if (!ksp) {
        biffAddf(GAGE, "%s: need kernel %s but it hasn't been set",
                 me, airEnumStr(gageKernel, kernIdx));
        return 1;
      }
      rad = ksp->kernel->support(ksp->parm);
      maxRad = AIR_MAX(maxRad, rad);
      if (ctx->verbose) {
        fprintf(stderr, "%s: k[%s]=%s -> rad = %g -> maxRad = %g\n", me,
                airEnumStr(gageKernel, kernIdx), ksp->kernel->name,
                rad, maxRad);
      }
    }
  }
  radius = AIR_ROUNDUP(maxRad);
  /* In case either kernels have tiny supports (less than 0.5), or if
     we in fact don't need any kernels, then we need to do this to
     ensure that we generate a valid radius */
  radius = AIR_MAX(radius, 1);
  if (ctx->parm.stackUse && (nrrdKernelHermiteScaleSpaceFlag
                             == ctx->ksp[gageKernelStack]->kernel)) {
    if (ctx->verbose) {
      fprintf(stderr, "%s: hermite on stack: bumping radius %d --> %d\n",
              me, radius, radius+1);
    }
    radius += 1;
  }
  if (radius != ctx->radius) {
    if (ctx->verbose) {
      fprintf(stderr, "%s: changing radius from %d to %d\n",
              me, ctx->radius, radius);
    }
    ctx->radius = radius;
    ctx->flag[gageCtxFlagRadius] = AIR_TRUE;
  }
  if (ctx->verbose) fprintf(stderr, "%s: bye\n", me);

  return 0;
}

int
_gageCacheSizeUpdate(gageContext *ctx) {
  static const char me[]="_gageCacheSizeUpdate";
  int fd;
  gagePerVolume *pvl;
  unsigned int pvlIdx;

  if (ctx->verbose) fprintf(stderr, "%s: hello (radius = %d)\n",
                            me, ctx->radius);
  if (!( ctx->radius > 0 )) {
    biffAddf(GAGE, "%s: have bad radius %d", me, ctx->radius);
    return 1;
  }
  fd = 2*ctx->radius;
  ctx->fsl = (double *)airFree(ctx->fsl);
  ctx->fw = (double *)airFree(ctx->fw);
  ctx->off = (unsigned int *)airFree(ctx->off);
  ctx->fsl = (double *)calloc(fd*3, sizeof(double));
  ctx->fw = (double *)calloc(fd*3*(GAGE_KERNEL_MAX+1), sizeof(double));
  ctx->off = (unsigned int *)calloc(fd*fd*fd, sizeof(unsigned int));
  if (!(ctx->fsl && ctx->fw && ctx->off)) {
    biffAddf(GAGE, "%s: couldn't allocate filter caches for fd=%d", me, fd);
    return 1;
  }
  /* note that all pvls get caches allocated for the same size, even if
     their queries don't involve the largest-size kernels. This is actually
     the feature that enabled the stack functionality to be easily added. */
  for (pvlIdx=0; pvlIdx<ctx->pvlNum; pvlIdx++) {
    pvl = ctx->pvl[pvlIdx];
    pvl->iv3 = (double *)airFree(pvl->iv3);
    pvl->iv2 = (double *)airFree(pvl->iv2);
    pvl->iv1 = (double *)airFree(pvl->iv1);
    pvl->iv3 = (double *)calloc(fd*fd*fd*pvl->kind->valLen, sizeof(double));
    pvl->iv2 = (double *)calloc(fd*fd*pvl->kind->valLen, sizeof(double));
    pvl->iv1 = (double *)calloc(fd*pvl->kind->valLen, sizeof(double));
    if (!(pvl->iv3 && pvl->iv2 && pvl->iv1)) {
      biffAddf(GAGE, "%s: couldn't allocate pvl[%d]'s value caches for fd=%d",
               me, pvlIdx, fd);
      return 1;
    }
  }
  if (ctx->verbose) fprintf(stderr, "%s: bye\n", me);

  return 0;
}

void
_gageOffValueUpdate(gageContext *ctx) {
  static const char me[]="_gageOffValueUpdate";
  int fd, i, j, k;
  unsigned int sx, sy;

  if (ctx->verbose) fprintf(stderr, "%s: hello\n", me);

  sx = ctx->shape->size[0];
  sy = ctx->shape->size[1];
  fd = 2*ctx->radius;
  /* HEY: look into special casing this for small fd */
  for (k=0; k<fd; k++) {
    for (j=0; j<fd; j++) {
      for (i=0; i<fd; i++) {
        ctx->off[i + fd*(j + fd*k)] = i + sx*(j + sy*k);
      }
    }
  }
  /* no flags to set for further action */
  if (ctx->verbose) fprintf(stderr, "%s: bye\n", me);

  return;
}

/*
******** gageUpdate()
**
** call just before probing begins.
*/
int
gageUpdate(gageContext *ctx) {
  static const char me[]="gageUpdate";
  unsigned int pi;
  int i, haveQuery;

  if (!( ctx )) {
    biffAddf(GAGE, "%s: got NULL pointer", me);
    return 1;
  }
  if (0 == ctx->pvlNum) {
    biffAddf(GAGE, "%s: context has no attached pervolumes", me);
    return 1;
  }
  haveQuery = AIR_FALSE;
  for (pi=0; pi<ctx->pvlNum; pi++) {
    haveQuery |= GAGE_QUERY_NONZERO(ctx->pvl[pi]->query);
  }
  if (!haveQuery) {
    biffAddf(GAGE, "%s: no query item set in %s", me,
             (ctx->pvlNum == 1
              ? "the pervolume"
              : "any of the pervolumes"));
    return 1;
  }

  /* HEY: shouldn't there be some more logic/state for this? */
  if (ctx->parm.stackUse) {
    if (!ctx->ksp[gageKernelStack]) {
      biffAddf(GAGE, "%s: can't do stack without ksp[%s]", me,
               airEnumStr(gageKernel, gageKernelStack));
      return 1;
    }
    if (!( 2 <= ctx->pvlNum )) {
      biffAddf(GAGE, "%s: need at least 2 pervolumes for stack", me);
      return 1;
    }
    for (pi=1; pi<ctx->pvlNum; pi++) {
      if (!( ctx->pvl[0]->kind == ctx->pvl[pi]->kind )) {
        biffAddf(GAGE, "%s: pvl[%u] kind (%s) != pvl[0] kind (%s)", me,
                 pi, ctx->pvl[pi]->kind->name, ctx->pvl[0]->kind->name);
        return 1;
      }
    }
  }

  /* start traversing the whole update graph . . . */
  if (ctx->verbose) {
    fprintf(stderr, "%s: hello ____________________ \n", me);
    fprintf(stderr, "    context flags:");
    for (i=gageCtxFlagUnknown+1; i<gageCtxFlagLast; i++) {
      fprintf(stderr, " %d=%d", i, ctx->flag[i]);
    }
    fprintf(stderr, "\n");
    fprintf(stderr, "    pvl flags:");
    for (i=gagePvlFlagUnknown+1; i<gagePvlFlagLast; i++) {
      fprintf(stderr, " %d=%d", i, _gagePvlFlagCheck(ctx, i));
    }
    fprintf(stderr, "\n");
  }
  if (_gagePvlFlagCheck(ctx, gagePvlFlagQuery)) {
    _gagePvlNeedDUpdate(ctx);
    _gagePvlFlagDown(ctx, gagePvlFlagQuery);
  }
  if (_gagePvlFlagCheck(ctx, gagePvlFlagNeedD)) {
    _gageNeedDUpdate(ctx);
    _gagePvlFlagDown(ctx, gagePvlFlagNeedD);
  }
  if (ctx->flag[gageCtxFlagNeedD] || ctx->flag[gageCtxFlagK3Pack]) {
    if (_gageNeedKUpdate(ctx)) {
      biffAddf(GAGE, "%s: trouble", me); return 1;
    }
    ctx->flag[gageCtxFlagNeedD] = AIR_FALSE;
    ctx->flag[gageCtxFlagK3Pack] = AIR_FALSE;
  }
  if (ctx->flag[gageCtxFlagKernel] || ctx->flag[gageCtxFlagNeedK]) {
    if (_gageRadiusUpdate(ctx)) {
      biffAddf(GAGE, "%s: trouble", me); return 1;
    }
    ctx->flag[gageCtxFlagKernel] = AIR_FALSE;
    ctx->flag[gageCtxFlagNeedK] = AIR_FALSE;
  }
  if (ctx->flag[gageCtxFlagRadius]
      /* HEY HEY HEY: this is a total hack: right now its possible for a
         new pvl to have unallocated iv3,iv2,iv1, if it was attached to a
         context which had already been probing, as was the case with
         _tenRegisterDoit.  So, with this hack we reallocate ALL caches
         just because a new pervolume was attached . . . */
      || _gagePvlFlagCheck(ctx, gagePvlFlagVolume)) {
    if (_gageCacheSizeUpdate(ctx)) {
      biffAddf(GAGE, "%s: trouble", me); return 1;
    }
  }
  if (ctx->flag[gageCtxFlagRadius]
      || ctx->flag[gageCtxFlagShape]
      /* see above; following flags that triggered _gageCacheSizeUpdate(ctx) */
      || _gagePvlFlagCheck(ctx, gagePvlFlagVolume)) {
    _gageOffValueUpdate(ctx);
    ctx->flag[gageCtxFlagShape] = AIR_FALSE;
  }
  ctx->flag[gageCtxFlagRadius] = AIR_FALSE;

  /* chances are, something above has invalidated the state maintained
     during successive calls to gageProbe() */
  gagePointReset(&ctx->point);

  for (pi=0; pi<ctx->pvlNum; pi++) {
    if (ctx->pvl[pi]->kind->pvlDataUpdate) {
      if (ctx->pvl[pi]->kind->pvlDataUpdate(ctx->pvl[pi]->kind,
                                            ctx,
                                            ctx->pvl[pi],
                                            ctx->pvl[pi]->data)) {
        biffAddf(GAGE, "%s: pvlDataUpdate(pvl[%u]) failed", me, pi);
        return 1;
      }
    }
  }

  if (ctx->verbose > 3 && ctx->stackPos) {
    fprintf(stderr, "%s: pvlNum = %u -> stack of %u [0,%u]\n", me,
            ctx->pvlNum, ctx->pvlNum-1, ctx->pvlNum-2);
    for (pi=0; pi<ctx->pvlNum-1; pi++) {
      fprintf(stderr, "%s: stackPos[%u] = %g\n", me, pi, ctx->stackPos[pi]);
    }
  }
  if (ctx->verbose) fprintf(stderr, "%s: bye ^^^^^^^^^^^^^^^^^^^ \n", me);

  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 "gage.h"
25			#include "privateGage.h"
26
27			int
28			_gagePvlFlagCheck(gageContext *ctx, int pvlFlag) {
29			int ret;
30			unsigned int pvlIdx;
31
32			ret = AIR_FALSE;
33	✓✓	3906	for (pvlIdx=0; pvlIdx<ctx->pvlNum; pvlIdx++) {
34		1494	ret \|= ctx->pvl[pvlIdx]->flag[pvlFlag];
35			}
36		306	return ret;
37			}
38
39			void
40			_gagePvlFlagDown(gageContext *ctx, int pvlFlag) {
41			unsigned int pvlIdx;
42
43	✓✓	2986	for (pvlIdx=0; pvlIdx<ctx->pvlNum; pvlIdx++) {
44		1094	ctx->pvl[pvlIdx]->flag[pvlFlag] = AIR_FALSE;
45			}
46		266	}
47
48			/*
49			** One could go from all the pvls' queries to the context's needD in
50			** one shot, but doing it in two steps (as below) seems a little clearer,
51			** and it means that pvl->needD isn't needlessly re-computed for
52			** pvl's whose query hasn't changed.
53			*/
54
55			/*
56			** for each pvl: pvl's query --> pvl's needD
57			*/
58			void
59			_gagePvlNeedDUpdate(gageContext *ctx) {
60			static const char me[]="_gagePvlNeedDUpdate";
61			gagePerVolume *pvl;
62		266	int que, needD[GAGE_DERIV_MAX+1];
63			unsigned int pvlIdx, di;
64
65	✗✓	133	if (ctx->verbose) fprintf(stderr, "%s: hello\n", me);
66	✓✓	1360	for (pvlIdx=0; pvlIdx<ctx->pvlNum; pvlIdx++) {
67		547	pvl = ctx->pvl[pvlIdx];
68	✓✗	547	if (pvl->flag[gagePvlFlagQuery]) {
69		547	GAGE_DV_SET(needD, 0, 0, 0);
70		547	que = pvl->kind->itemMax+1;
71		547	do {
72		21559	que--;
73	✓✓	21559	if (GAGE_QUERY_ITEM_TEST(pvl->query, que)) {
74		1457	needD[pvl->kind->table[que].needDeriv] = 1;
75		1457	}
76	✓✓	21559	} while (que);
77	✓✓✗✓ ✗✗	747	if (!GAGE_DV_EQUAL(needD, pvl->needD)) {
78	✗✓	547	if (ctx->verbose) {
79			fprintf(stderr, "%s: updating pvl[%d]'s needD to (", me, pvlIdx);
80			for (di=0; di<=GAGE_DERIV_MAX; di++) {
81			fprintf(stderr, "%s%d", di ? "," : "", needD[di]);
82			}
83			fprintf(stderr, "\n");
84			}
85		547	GAGE_DV_COPY(pvl->needD, needD);
86		547	pvl->flag[gagePvlFlagNeedD] = AIR_TRUE;
87		547	}
88			}
89			}
90	✗✓	133	if (ctx->verbose) fprintf(stderr, "%s: bye\n", me);
91			return;
92		133	}
93
94			/*
95			** all pvls' needD --> ctx's needD
96			*/
97			void
98			_gageNeedDUpdate(gageContext *ctx) {
99			static const char me[]="_gageNeedDUpdate";
100			gagePerVolume *pvl;
101		266	int needD[GAGE_DERIV_MAX+1];
102			unsigned int pvlIdx, di;
103
104	✗✓	133	if (ctx->verbose) fprintf(stderr, "%s: hello\n", me);
105		133	GAGE_DV_SET(needD, 0, 0, 0);
106	✓✓	1360	for (pvlIdx=0; pvlIdx<ctx->pvlNum; pvlIdx++) {
107		547	pvl = ctx->pvl[pvlIdx];
108	✓✓	4376	for (di=0; di<=GAGE_DERIV_MAX; di++) {
109		1641	needD[di] \|= pvl->needD[di];
110			}
111			}
112	✓✓✗✓ ✗✗	153	if (!GAGE_DV_EQUAL(needD, ctx->needD)) {
113	✗✓	133	if (ctx->verbose) {
114			fprintf(stderr, "%s: updating ctx's needD to (", me);
115			for (di=0; di<=GAGE_DERIV_MAX; di++) {
116			fprintf(stderr, "%s%d", di ? "," : "", needD[di]);
117			}
118			fprintf(stderr, "\n");
119			}
120		133	GAGE_DV_COPY(ctx->needD, needD);
121		133	ctx->flag[gageCtxFlagNeedD] = AIR_TRUE;
122		133	}
123	✗✓	133	if (ctx->verbose) fprintf(stderr, "%s: bye\n", me);
124
125			return;
126		133	}
127
128			/*
129			** ctx's needD & k3pack --> needK
130			*/
131			static int
132			_gageNeedKUpdate(gageContext *ctx) {
133			static const char me[]="_gageNeedKUpdate";
134		266	int kernIdx, needK[GAGE_KERNEL_MAX+1], change;
135			unsigned int di;
136
137	✗✓	133	if (ctx->verbose) fprintf(stderr, "%s: hello\n", me);
138	✓✓	2128	for (kernIdx=gageKernelUnknown+1; kernIdx<gageKernelLast; kernIdx++) {
139		931	needK[kernIdx] = AIR_FALSE;
140			}
141	✗✓	133	if (!ctx->parm.k3pack) {
142			biffAddf(GAGE, "%s: sorry, only 3-pack filtering implemented now", me);
143			return 1;
144			}
145			di = 0;
146	✓✗	133	if (ctx->needD[di]) {
147		133	needK[gageKernel00] = AIR_TRUE;
148		133	}
149			di += 1;
150	✓✓	133	if (ctx->needD[di]) {
151		111	needK[gageKernel00] = needK[gageKernel11] = AIR_TRUE;
152		111	}
153			di += 1;
154	✓✓	133	if (ctx->needD[di]) {
155		111	needK[gageKernel00] = needK[gageKernel11] =
156		111	needK[gageKernel22] = AIR_TRUE;
157		111	}
158	✗✓	133	if (GAGE_DERIV_MAX != di) {
159			biffAddf(GAGE, "%s: sorry, code needs updating for GAGE_DERIV_MAX %u",
160			me, GAGE_DERIV_MAX);
161			return 1;
162			}
163			change = AIR_FALSE;
164	✓✓	2128	for (kernIdx=gageKernelUnknown+1; kernIdx<gageKernelLast; kernIdx++) {
165		931	change \|= (needK[kernIdx] != ctx->needK[kernIdx]);
166			}
167	✓✗	133	if (change) {
168	✗✓	133	if (ctx->verbose) {
169			fprintf(stderr, "%s: changing needK to (", me);
170			for (kernIdx=gageKernelUnknown+1; kernIdx<gageKernelLast; kernIdx++) {
171			fprintf(stderr, "%s%d", kernIdx > gageKernelUnknown+1 ? "," : "",
172			needK[kernIdx]);
173			}
174			fprintf(stderr, ")\n");
175			}
176	✓✓	2128	for (kernIdx=gageKernelUnknown+1; kernIdx<gageKernelLast; kernIdx++) {
177		931	ctx->needK[kernIdx] = needK[kernIdx];
178			}
179		133	ctx->flag[gageCtxFlagNeedK] = AIR_TRUE;
180		133	}
181	✗✓	133	if (ctx->verbose) fprintf(stderr, "%s: bye\n", me);
182
183		133	return 0;
184		133	}
185
186			/*
187			** ctx's ksp[] & needK --> radius
188			**
189			*/
190			int
191			_gageRadiusUpdate(gageContext *ctx) {
192			static const char me[]="_gageRadiusUpdate";
193			unsigned int kernIdx, radius;
194			double maxRad, rad;
195			NrrdKernelSpec *ksp;
196
197	✗✓	266	if (ctx->verbose) fprintf(stderr, "%s: hello\n", me);
198			maxRad = 0;
199	✓✓	2128	for (kernIdx=gageKernelUnknown+1; kernIdx<gageKernelLast; kernIdx++) {
200	✓✓	931	if (ctx->needK[kernIdx]) {
201		355	ksp = ctx->ksp[kernIdx];
202	✗✓	355	if (!ksp) {
203			biffAddf(GAGE, "%s: need kernel %s but it hasn't been set",
204			me, airEnumStr(gageKernel, kernIdx));
205			return 1;
206			}
207		355	rad = ksp->kernel->support(ksp->parm);
208		355	maxRad = AIR_MAX(maxRad, rad);
209	✗✓	355	if (ctx->verbose) {
210			fprintf(stderr, "%s: k[%s]=%s -> rad = %g -> maxRad = %g\n", me,
211			airEnumStr(gageKernel, kernIdx), ksp->kernel->name,
212			rad, maxRad);
213			}
214			}
215			}
216		133	radius = AIR_ROUNDUP(maxRad);
217			/* In case either kernels have tiny supports (less than 0.5), or if
218			we in fact don't need any kernels, then we need to do this to
219			ensure that we generate a valid radius */
220		133	radius = AIR_MAX(radius, 1);
221	✗✓✗✗ ✗✗	133	if (ctx->parm.stackUse && (nrrdKernelHermiteScaleSpaceFlag
222			== ctx->ksp[gageKernelStack]->kernel)) {
223			if (ctx->verbose) {
224			fprintf(stderr, "%s: hermite on stack: bumping radius %d --> %d\n",
225			me, radius, radius+1);
226			}
227			radius += 1;
228			}
229	✓✓	133	if (radius != ctx->radius) {
230	✗✓	113	if (ctx->verbose) {
231			fprintf(stderr, "%s: changing radius from %d to %d\n",
232			me, ctx->radius, radius);
233			}
234		113	ctx->radius = radius;
235		113	ctx->flag[gageCtxFlagRadius] = AIR_TRUE;
236		113	}
237	✗✓	133	if (ctx->verbose) fprintf(stderr, "%s: bye\n", me);
238
239		133	return 0;
240		133	}
241
242			int
243			_gageCacheSizeUpdate(gageContext *ctx) {
244			static const char me[]="_gageCacheSizeUpdate";
245			int fd;
246			gagePerVolume *pvl;
247			unsigned int pvlIdx;
248
249	✗✓	266	if (ctx->verbose) fprintf(stderr, "%s: hello (radius = %d)\n",
250			me, ctx->radius);
251	✗✓	133	if (!( ctx->radius > 0 )) {
252			biffAddf(GAGE, "%s: have bad radius %d", me, ctx->radius);
253			return 1;
254			}
255		133	fd = 2*ctx->radius;
256		133	ctx->fsl = (double *)airFree(ctx->fsl);
257		133	ctx->fw = (double *)airFree(ctx->fw);
258		133	ctx->off = (unsigned int *)airFree(ctx->off);
259		133	ctx->fsl = (double )calloc(fd3, sizeof(double));
260		133	ctx->fw = (double )calloc(fd3*(GAGE_KERNEL_MAX+1), sizeof(double));
261		133	ctx->off = (unsigned int )calloc(fdfd*fd, sizeof(unsigned int));
262	✓✗✓✗ ✗✓	399	if (!(ctx->fsl && ctx->fw && ctx->off)) {
263			biffAddf(GAGE, "%s: couldn't allocate filter caches for fd=%d", me, fd);
264			return 1;
265			}
266			/* note that all pvls get caches allocated for the same size, even if
267			their queries don't involve the largest-size kernels. This is actually
268			the feature that enabled the stack functionality to be easily added. */
269	✓✓	1360	for (pvlIdx=0; pvlIdx<ctx->pvlNum; pvlIdx++) {
270		547	pvl = ctx->pvl[pvlIdx];
271		547	pvl->iv3 = (double *)airFree(pvl->iv3);
272		547	pvl->iv2 = (double *)airFree(pvl->iv2);
273		547	pvl->iv1 = (double *)airFree(pvl->iv1);
274		547	pvl->iv3 = (double )calloc(fdfdfdpvl->kind->valLen, sizeof(double));
275		547	pvl->iv2 = (double )calloc(fdfd*pvl->kind->valLen, sizeof(double));
276		547	pvl->iv1 = (double )calloc(fdpvl->kind->valLen, sizeof(double));
277	✓✗✓✗ ✗✓	1641	if (!(pvl->iv3 && pvl->iv2 && pvl->iv1)) {
278			biffAddf(GAGE, "%s: couldn't allocate pvl[%d]'s value caches for fd=%d",
279			me, pvlIdx, fd);
280			return 1;
281			}
282			}
283	✗✓	133	if (ctx->verbose) fprintf(stderr, "%s: bye\n", me);
284
285		133	return 0;
286		133	}
287
288			void
289			_gageOffValueUpdate(gageContext *ctx) {
290			static const char me[]="_gageOffValueUpdate";
291			int fd, i, j, k;
292			unsigned int sx, sy;
293
294	✗✓	266	if (ctx->verbose) fprintf(stderr, "%s: hello\n", me);
295
296		133	sx = ctx->shape->size[0];
297		133	sy = ctx->shape->size[1];
298		133	fd = 2*ctx->radius;
299			/* HEY: look into special casing this for small fd */
300	✓✓	2062	for (k=0; k<fd; k++) {
301	✓✓	20876	for (j=0; j<fd; j++) {
302	✓✓	287096	for (i=0; i<fd; i++) {
303		134008	ctx->off[i + fd(j + fdk)] = i + sx(j + syk);
304			}
305			}
306			}
307			/* no flags to set for further action */
308	✗✓	133	if (ctx->verbose) fprintf(stderr, "%s: bye\n", me);
309
310			return;
311		133	}
312
313			/*
314			******** gageUpdate()
315			**
316			** call just before probing begins.
317			*/
318			int
319			gageUpdate(gageContext *ctx) {
320			static const char me[]="gageUpdate";
321			unsigned int pi;
322			int i, haveQuery;
323
324	✗✓	266	if (!( ctx )) {
325			biffAddf(GAGE, "%s: got NULL pointer", me);
326			return 1;
327			}
328	✗✓	133	if (0 == ctx->pvlNum) {
329			biffAddf(GAGE, "%s: context has no attached pervolumes", me);
330			return 1;
331			}
332			haveQuery = AIR_FALSE;
333	✓✓	1360	for (pi=0; pi<ctx->pvlNum; pi++) {
334		547	haveQuery \|= GAGE_QUERY_NONZERO(ctx->pvl[pi]->query);
335			}
336	✗✓	133	if (!haveQuery) {
337			biffAddf(GAGE, "%s: no query item set in %s", me,
338			(ctx->pvlNum == 1
339			? "the pervolume"
340			: "any of the pervolumes"));
341			return 1;
342			}
343
344			/* HEY: shouldn't there be some more logic/state for this? */
345	✗✓	133	if (ctx->parm.stackUse) {
346			if (!ctx->ksp[gageKernelStack]) {
347			biffAddf(GAGE, "%s: can't do stack without ksp[%s]", me,
348			airEnumStr(gageKernel, gageKernelStack));
349			return 1;
350			}
351			if (!( 2 <= ctx->pvlNum )) {
352			biffAddf(GAGE, "%s: need at least 2 pervolumes for stack", me);
353			return 1;
354			}
355			for (pi=1; pi<ctx->pvlNum; pi++) {
356			if (!( ctx->pvl[0]->kind == ctx->pvl[pi]->kind )) {
357			biffAddf(GAGE, "%s: pvl[%u] kind (%s) != pvl[0] kind (%s)", me,
358			pi, ctx->pvl[pi]->kind->name, ctx->pvl[0]->kind->name);
359			return 1;
360			}
361			}
362			}
363
364			/* start traversing the whole update graph . . . */
365	✗✓	133	if (ctx->verbose) {
366			fprintf(stderr, "%s: hello ____________________ \n", me);
367			fprintf(stderr, " context flags:");
368			for (i=gageCtxFlagUnknown+1; i<gageCtxFlagLast; i++) {
369			fprintf(stderr, " %d=%d", i, ctx->flag[i]);
370			}
371			fprintf(stderr, "\n");
372			fprintf(stderr, " pvl flags:");
373			for (i=gagePvlFlagUnknown+1; i<gagePvlFlagLast; i++) {
374			fprintf(stderr, " %d=%d", i, _gagePvlFlagCheck(ctx, i));
375			}
376			fprintf(stderr, "\n");
377			}
378	✓✗	133	if (_gagePvlFlagCheck(ctx, gagePvlFlagQuery)) {
379		133	_gagePvlNeedDUpdate(ctx);
380		133	_gagePvlFlagDown(ctx, gagePvlFlagQuery);
381		133	}
382	✓✗	133	if (_gagePvlFlagCheck(ctx, gagePvlFlagNeedD)) {
383		133	_gageNeedDUpdate(ctx);
384		133	_gagePvlFlagDown(ctx, gagePvlFlagNeedD);
385		133	}
386	✗✓✗✗	133	if (ctx->flag[gageCtxFlagNeedD] \|\| ctx->flag[gageCtxFlagK3Pack]) {
387	✗✓	133	if (_gageNeedKUpdate(ctx)) {
388			biffAddf(GAGE, "%s: trouble", me); return 1;
389			}
390		133	ctx->flag[gageCtxFlagNeedD] = AIR_FALSE;
391		133	ctx->flag[gageCtxFlagK3Pack] = AIR_FALSE;
392		133	}
393	✗✓✗✗	133	if (ctx->flag[gageCtxFlagKernel] \|\| ctx->flag[gageCtxFlagNeedK]) {
394	✗✓	133	if (_gageRadiusUpdate(ctx)) {
395			biffAddf(GAGE, "%s: trouble", me); return 1;
396			}
397		133	ctx->flag[gageCtxFlagKernel] = AIR_FALSE;
398		133	ctx->flag[gageCtxFlagNeedK] = AIR_FALSE;
399		133	}
400	✓✗	153	if (ctx->flag[gageCtxFlagRadius]
401			/* HEY HEY HEY: this is a total hack: right now its possible for a
402			new pvl to have unallocated iv3,iv2,iv1, if it was attached to a
403			context which had already been probing, as was the case with
404			_tenRegisterDoit. So, with this hack we reallocate ALL caches
405			just because a new pervolume was attached . . . */
406	✓✓	153	\|\| _gagePvlFlagCheck(ctx, gagePvlFlagVolume)) {
407	✗✓	133	if (_gageCacheSizeUpdate(ctx)) {
408			biffAddf(GAGE, "%s: trouble", me); return 1;
409			}
410			}
411	✓✗	153	if (ctx->flag[gageCtxFlagRadius]
412	✓✓	153	\|\| ctx->flag[gageCtxFlagShape]
413			/* see above; following flags that triggered _gageCacheSizeUpdate(ctx) */
414	✓✗	40	\|\| _gagePvlFlagCheck(ctx, gagePvlFlagVolume)) {
415		133	_gageOffValueUpdate(ctx);
416		133	ctx->flag[gageCtxFlagShape] = AIR_FALSE;
417		133	}
418		133	ctx->flag[gageCtxFlagRadius] = AIR_FALSE;
419
420			/* chances are, something above has invalidated the state maintained
421			during successive calls to gageProbe() */
422		133	gagePointReset(&ctx->point);
423
424	✓✓	1360	for (pi=0; pi<ctx->pvlNum; pi++) {
425	✓✓	547	if (ctx->pvl[pi]->kind->pvlDataUpdate) {
426	✗✓	16	if (ctx->pvl[pi]->kind->pvlDataUpdate(ctx->pvl[pi]->kind,
427			ctx,
428			ctx->pvl[pi],
429		16	ctx->pvl[pi]->data)) {
430			biffAddf(GAGE, "%s: pvlDataUpdate(pvl[%u]) failed", me, pi);
431			return 1;
432			}
433			}
434			}
435
436	✗✓✗✗	133	if (ctx->verbose > 3 && ctx->stackPos) {
437			fprintf(stderr, "%s: pvlNum = %u -> stack of %u [0,%u]\n", me,
438			ctx->pvlNum, ctx->pvlNum-1, ctx->pvlNum-2);
439			for (pi=0; pi<ctx->pvlNum-1; pi++) {
440			fprintf(stderr, "%s: stackPos[%u] = %g\n", me, pi, ctx->stackPos[pi]);
441			}
442			}
443	✗✓	133	if (ctx->verbose) fprintf(stderr, "%s: bye ^^^^^^^^^^^^^^^^^^^ \n", me);
444
445		133	return 0;
446		133	}