Head

GCC Code Coverage Report

Directory:	./		Exec	Total	Coverage
File:	src/gage/st.c	Lines:	0	153	0.0 %
Date:	2017-05-26	Branches:	0	94	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 "gage.h"
#include "privateGage.h"

#define GAGE_CACHE_LEN 1013

unsigned int
_gageHash(int x, int y, int z) {
  unsigned int h, g;
  unsigned char s[6];
  int i;

  s[0] = x & 0xff;
  s[1] = (x >> 8) & 0xff;
  s[2] = y & 0xff;
  s[3] = (y >> 8) & 0xff;
  s[4] = z & 0xff;
  s[5] = (z >> 8) & 0xff;

  h = 0;
  for (i=0; i<=5; i++) {
    h = (h << 4) + s[i];
    if ((g = h & 0xf0000000)) {
      h = h ^ (g >> 24);
      h = h ^ g;
    }
  }
  return h % GAGE_CACHE_LEN;
}

void
_gageCacheProbe(gageContext *ctx, double *grad,
                int *cc, double *gc,
                int x, int y, int z) {
  int hi;

  hi = _gageHash(x, y, z);
  if ( (cc[3*hi + 0] == x) &&
       (cc[3*hi + 1] == y) &&
       (cc[3*hi + 2] == z) ) {
    /* cache hit */
    ELL_3V_COPY(grad, gc + 3*hi);
  } else {
    /* cache miss */
    cc[3*hi + 0] = x;
    cc[3*hi + 1] = y;
    cc[3*hi + 2] = z;
    gageProbe(ctx, x, y, z);
    ELL_3V_COPY(gc + 3*hi, grad);
  }
  return ;
}

/*
******** gageStructureTensor()
**
** Computes volume of structure tensors.  Currently, only works on a scalar
** fields (for multi-scalar fields, just add structure tensors from each
** component scalar), and only uses the B-spline kernel for differentiation
** and derivative blurring.
**
** Note, if you want to use dsmp > 1, its your responsibility to give
** an appropriate iScale > 1, so that you don't undersample.
*/
int
gageStructureTensor(Nrrd *nout, const Nrrd *nin,
                    int dScale, int iScale, int dsmp) {
  static const char me[]="gageStructureTensor";
  NrrdKernelSpec *gk0, *gk1, *ik0;
  int E, rad, diam, osx, osy, osz, oxi, oyi, ozi,
    _ixi, _iyi, _izi, ixi, iyi, izi, wi, *coordCache;
  gageContext *ctx;
  gageQuery query;
  gagePerVolume *pvl;
  airArray *mop;
  double *grad, *ixw, *iyw, *izw, wght, sten[6], *gradCache, *out;
  double xs, ys, zs, ms;

  if (!(nout && nin)) {
    biffAddf(GAGE, "%s: got NULL pointer", me);
    return 1;
  }
  if (!( 3 == nin->dim && nrrdTypeBlock != nin->type )) {
    biffAddf(GAGE, "%s: nin isn't a 3D non-block type nrrd", me);
    return 1;
  }
  /*
  if (!( AIR_EXISTS(nin->axis[0].spacing)
         && AIR_EXISTS(nin->axis[1].spacing)
         && AIR_EXISTS(nin->axis[2].spacing) )) {
    biffAddf(GAGE, "%s: nin's axis 0,1,2 spacings don't all exist", me);
    return 1;
  }
  */
  if (!( dScale >= 1 && iScale >= 1 && dsmp >= 1 )) {
    biffAddf(GAGE, "%s: dScale (%d), iScale (%d), dsmp (%d) not all >= 1",
             me, dScale, iScale, dsmp);
    return 1;
  }

  mop = airMopNew();
  gk0 = nrrdKernelSpecNew();
  gk1 = nrrdKernelSpecNew();
  ik0 = nrrdKernelSpecNew();
  airMopAdd(mop, gk0, (airMopper)nrrdKernelSpecNix, airMopAlways);
  airMopAdd(mop, gk1, (airMopper)nrrdKernelSpecNix, airMopAlways);
  airMopAdd(mop, ik0, (airMopper)nrrdKernelSpecNix, airMopAlways);
  if ( nrrdKernelSpecParse(gk0, "cubic:1,0")
       || nrrdKernelSpecParse(gk1, "cubicd:1,0")
       || nrrdKernelSpecParse(ik0, "cubic:1,0")) {
    biffMovef(GAGE, NRRD, "%s: couldn't set up kernels", me);
    airMopError(mop); return 1;
  }
  /* manually set scale parameters */
  gk0->parm[0] = dScale;
  gk1->parm[0] = dScale;
  ik0->parm[0] = 1.0;       /* this is more complicated . . . */
  ctx = gageContextNew();
  airMopAdd(mop, ctx, (airMopper)gageContextNix, airMopAlways);
  gageParmSet(ctx, gageParmRenormalize, AIR_TRUE);
  E = 0;
  if (!E) E |= !(pvl = gagePerVolumeNew(ctx, nin, gageKindScl));
  if (!E) E |= gagePerVolumeAttach(ctx, pvl);
  if (!E) E |= gageKernelSet(ctx, gageKernel00, gk0->kernel, gk0->parm);
  if (!E) E |= gageKernelSet(ctx, gageKernel11, gk1->kernel, gk1->parm);
  if (!E) GAGE_QUERY_RESET(query);
  if (!E) GAGE_QUERY_ITEM_ON(query, gageSclGradVec);
  if (!E) E |= gageQuerySet(ctx, pvl, query);
  if (!E) E |= gageUpdate(ctx);
  if (E) {
    biffAddf(GAGE, "%s: ", me);
    airMopError(mop); return 1;
  }
  grad = _gageAnswerPointer(ctx, pvl, gageSclGradVec);

  xs = nin->axis[0].spacing;
  ys = nin->axis[1].spacing;
  zs = nin->axis[2].spacing;
  xs = AIR_EXISTS(xs) ? AIR_ABS(xs) : 1.0;
  ys = AIR_EXISTS(ys) ? AIR_ABS(ys) : 1.0;
  zs = AIR_EXISTS(zs) ? AIR_ABS(zs) : 1.0;
  ms = airCbrt(xs*ys*zs);
  /* ms is geometric mean of {xs,ys,zs} */
  xs /= ms;
  ys /= ms;
  zs /= ms;
  fprintf(stderr, "iScale = %d, xs, ys, zs = %g, %g, %g\n",
          iScale, xs, ys, xs);

  rad = 0;
  ik0->parm[0] = iScale/xs;
  rad = AIR_MAX(rad, AIR_ROUNDUP(ik0->kernel->support(ik0->parm)));
  ik0->parm[0] = iScale/ys;
  rad = AIR_MAX(rad, AIR_ROUNDUP(ik0->kernel->support(ik0->parm)));
  ik0->parm[0] = iScale/zs;
  rad = AIR_MAX(rad, AIR_ROUNDUP(ik0->kernel->support(ik0->parm)));
  diam = 2*rad + 1;
  ixw = AIR_CALLOC(diam, double);
  iyw = AIR_CALLOC(diam, double);
  izw = AIR_CALLOC(diam, double);
  airMopAdd(mop, ixw, airFree, airMopAlways);
  airMopAdd(mop, iyw, airFree, airMopAlways);
  airMopAdd(mop, izw, airFree, airMopAlways);
  if (!(ixw && iyw && izw)) {
    biffAddf(GAGE, "%s: couldn't allocate grad vector or weight buffers", me);
    airMopError(mop); return 1;
  }

  /* the only reason that it is thread-safe to cache gageProbe results,
     without having the cache hang directly off the gageContext, is that
     we're doing all the probing for one context in one shot- producing
     an entirely volume of structure tensors with one function call */
  gradCache = AIR_CALLOC(3*GAGE_CACHE_LEN, double);
  coordCache = AIR_CALLOC(3*GAGE_CACHE_LEN, int);
  airMopAdd(mop, gradCache, airFree, airMopAlways);
  airMopAdd(mop, coordCache, airFree, airMopAlways);
  if (!(gradCache && coordCache)) {
    biffAddf(GAGE, "%s: couldn't allocate caches", me);
    airMopError(mop); return 1;
  }
  for (ixi=0; ixi<GAGE_CACHE_LEN; ixi++) {
    coordCache[3*ixi + 0] = -1;
    coordCache[3*ixi + 1] = -1;
    coordCache[3*ixi + 2] = -1;
  }

  for (wi=-rad; wi<=rad; wi++) {
    ik0->parm[0] = iScale/xs;
    ixw[wi+rad] = ik0->kernel->eval1_d(wi, ik0->parm);
    ik0->parm[0] = iScale/ys;
    iyw[wi+rad] = ik0->kernel->eval1_d(wi, ik0->parm);
    ik0->parm[0] = iScale/zs;
    izw[wi+rad] = ik0->kernel->eval1_d(wi, ik0->parm);
    fprintf(stderr, "%d --> (%g,%g,%g) -> (%g,%g,%g)\n",
            wi, wi/xs, wi/ys, wi/zs, ixw[wi+rad], iyw[wi+rad], izw[wi+rad]);
  }

  osx = (nin->axis[0].size)/dsmp;
  osy = (nin->axis[1].size)/dsmp;
  osz = (nin->axis[2].size)/dsmp;
  if (nrrdMaybeAlloc_va(nout, nrrdTypeDouble, 4,
                        AIR_CAST(size_t, 7),
                        AIR_CAST(size_t, osx),
                        AIR_CAST(size_t, osy),
                        AIR_CAST(size_t, osz))) {
    biffAddf(GAGE, NRRD, "%s: couldn't allocate output", me);
    airMopError(mop); return 1;
  }
  airMopAdd(mop, nout, (airMopper)nrrdEmpty, airMopOnError);

  out = AIR_CAST(double *, nout->data);
  for (ozi=0; ozi<osz; ozi++) {
    fprintf(stderr, "%s: z = %d/%d\n", me, ozi+1, osz);
    for (oyi=0; oyi<osy; oyi++) {
      for (oxi=0; oxi<osx; oxi++) {

        sten[0] = sten[1] = sten[2] = sten[3] = sten[4] = sten[5] = 0;
        for (_izi=0; _izi<diam; _izi++) {
          izi = AIR_CLAMP(0, _izi - rad + ozi*dsmp,
                          (int)nin->axis[2].size-1);
          if (!izw[_izi]) continue;
          for (_iyi=0; _iyi<diam; _iyi++) {
            iyi = AIR_CLAMP(0, _iyi - rad + oyi*dsmp,
                            (int)nin->axis[1].size-1);
            if (!iyw[_iyi]) continue;
            for (_ixi=0; _ixi<diam; _ixi++) {
              ixi = AIR_CLAMP(0, _ixi - rad + oxi*dsmp,
                              (int)nin->axis[0].size-1);
              if (!ixw[_ixi]) continue;
              wght = ixw[_ixi]*iyw[_iyi]*izw[_izi];
              _gageCacheProbe(ctx, grad, coordCache, gradCache, ixi, iyi, izi);
              sten[0] += wght*grad[0]*grad[0];
              sten[1] += wght*grad[0]*grad[1];
              sten[2] += wght*grad[0]*grad[2];
              sten[3] += wght*grad[1]*grad[1];
              sten[4] += wght*grad[1]*grad[2];
              sten[5] += wght*grad[2]*grad[2];
            }
          }
        }
        out[0] = 1.0;
        out[1] = sten[0];
        out[2] = sten[1];
        out[3] = sten[2];
        out[4] = sten[3];
        out[5] = sten[4];
        out[6] = sten[5];
        out += 7;

      }
    }
  }

  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 "gage.h"
25			#include "privateGage.h"
26
27			#define GAGE_CACHE_LEN 1013
28
29			unsigned int
30			_gageHash(int x, int y, int z) {
31			unsigned int h, g;
32			unsigned char s[6];
33			int i;
34
35			s[0] = x & 0xff;
36			s[1] = (x >> 8) & 0xff;
37			s[2] = y & 0xff;
38			s[3] = (y >> 8) & 0xff;
39			s[4] = z & 0xff;
40			s[5] = (z >> 8) & 0xff;
41
42			h = 0;
43			for (i=0; i<=5; i++) {
44			h = (h << 4) + s[i];
45			if ((g = h & 0xf0000000)) {
46			h = h ^ (g >> 24);
47			h = h ^ g;
48			}
49			}
50			return h % GAGE_CACHE_LEN;
51			}
52
53			void
54			_gageCacheProbe(gageContext ctx, double grad,
55			int cc, double gc,
56			int x, int y, int z) {
57			int hi;
58
59			hi = _gageHash(x, y, z);
60			if ( (cc[3*hi + 0] == x) &&
61			(cc[3*hi + 1] == y) &&
62			(cc[3*hi + 2] == z) ) {
63			/* cache hit */
64			ELL_3V_COPY(grad, gc + 3*hi);
65			} else {
66			/* cache miss */
67			cc[3*hi + 0] = x;
68			cc[3*hi + 1] = y;
69			cc[3*hi + 2] = z;
70			gageProbe(ctx, x, y, z);
71			ELL_3V_COPY(gc + 3*hi, grad);
72			}
73			return ;
74			}
75
76			/*
77			******** gageStructureTensor()
78			**
79			** Computes volume of structure tensors. Currently, only works on a scalar
80			** fields (for multi-scalar fields, just add structure tensors from each
81			** component scalar), and only uses the B-spline kernel for differentiation
82			** and derivative blurring.
83			**
84			** Note, if you want to use dsmp > 1, its your responsibility to give
85			** an appropriate iScale > 1, so that you don't undersample.
86			*/
87			int
88			gageStructureTensor(Nrrd nout, const Nrrd nin,
89			int dScale, int iScale, int dsmp) {
90			static const char me[]="gageStructureTensor";
91			NrrdKernelSpec gk0, gk1, *ik0;
92			int E, rad, diam, osx, osy, osz, oxi, oyi, ozi,
93			_ixi, _iyi, _izi, ixi, iyi, izi, wi, *coordCache;
94			gageContext *ctx;
95			gageQuery query;
96			gagePerVolume *pvl;
97			airArray *mop;
98			double grad, ixw, iyw, izw, wght, sten[6], gradCache, out;
99			double xs, ys, zs, ms;
100
101			if (!(nout && nin)) {
102			biffAddf(GAGE, "%s: got NULL pointer", me);
103			return 1;
104			}
105			if (!( 3 == nin->dim && nrrdTypeBlock != nin->type )) {
106			biffAddf(GAGE, "%s: nin isn't a 3D non-block type nrrd", me);
107			return 1;
108			}
109			/*
110			if (!( AIR_EXISTS(nin->axis[0].spacing)
111			&& AIR_EXISTS(nin->axis[1].spacing)
112			&& AIR_EXISTS(nin->axis[2].spacing) )) {
113			biffAddf(GAGE, "%s: nin's axis 0,1,2 spacings don't all exist", me);
114			return 1;
115			}
116			*/
117			if (!( dScale >= 1 && iScale >= 1 && dsmp >= 1 )) {
118			biffAddf(GAGE, "%s: dScale (%d), iScale (%d), dsmp (%d) not all >= 1",
119			me, dScale, iScale, dsmp);
120			return 1;
121			}
122
123			mop = airMopNew();
124			gk0 = nrrdKernelSpecNew();
125			gk1 = nrrdKernelSpecNew();
126			ik0 = nrrdKernelSpecNew();
127			airMopAdd(mop, gk0, (airMopper)nrrdKernelSpecNix, airMopAlways);
128			airMopAdd(mop, gk1, (airMopper)nrrdKernelSpecNix, airMopAlways);
129			airMopAdd(mop, ik0, (airMopper)nrrdKernelSpecNix, airMopAlways);
130			if ( nrrdKernelSpecParse(gk0, "cubic:1,0")
131			\|\| nrrdKernelSpecParse(gk1, "cubicd:1,0")
132			\|\| nrrdKernelSpecParse(ik0, "cubic:1,0")) {
133			biffMovef(GAGE, NRRD, "%s: couldn't set up kernels", me);
134			airMopError(mop); return 1;
135			}
136			/* manually set scale parameters */
137			gk0->parm[0] = dScale;
138			gk1->parm[0] = dScale;
139			ik0->parm[0] = 1.0; /* this is more complicated . . . */
140			ctx = gageContextNew();
141			airMopAdd(mop, ctx, (airMopper)gageContextNix, airMopAlways);
142			gageParmSet(ctx, gageParmRenormalize, AIR_TRUE);
143			E = 0;
144			if (!E) E \|= !(pvl = gagePerVolumeNew(ctx, nin, gageKindScl));
145			if (!E) E \|= gagePerVolumeAttach(ctx, pvl);
146			if (!E) E \|= gageKernelSet(ctx, gageKernel00, gk0->kernel, gk0->parm);
147			if (!E) E \|= gageKernelSet(ctx, gageKernel11, gk1->kernel, gk1->parm);
148			if (!E) GAGE_QUERY_RESET(query);
149			if (!E) GAGE_QUERY_ITEM_ON(query, gageSclGradVec);
150			if (!E) E \|= gageQuerySet(ctx, pvl, query);
151			if (!E) E \|= gageUpdate(ctx);
152			if (E) {
153			biffAddf(GAGE, "%s: ", me);
154			airMopError(mop); return 1;
155			}
156			grad = _gageAnswerPointer(ctx, pvl, gageSclGradVec);
157
158			xs = nin->axis[0].spacing;
159			ys = nin->axis[1].spacing;
160			zs = nin->axis[2].spacing;
161			xs = AIR_EXISTS(xs) ? AIR_ABS(xs) : 1.0;
162			ys = AIR_EXISTS(ys) ? AIR_ABS(ys) : 1.0;
163			zs = AIR_EXISTS(zs) ? AIR_ABS(zs) : 1.0;
164			ms = airCbrt(xsyszs);
165			/* ms is geometric mean of {xs,ys,zs} */
166			xs /= ms;
167			ys /= ms;
168			zs /= ms;
169			fprintf(stderr, "iScale = %d, xs, ys, zs = %g, %g, %g\n",
170			iScale, xs, ys, xs);
171
172			rad = 0;
173			ik0->parm[0] = iScale/xs;
174			rad = AIR_MAX(rad, AIR_ROUNDUP(ik0->kernel->support(ik0->parm)));
175			ik0->parm[0] = iScale/ys;
176			rad = AIR_MAX(rad, AIR_ROUNDUP(ik0->kernel->support(ik0->parm)));
177			ik0->parm[0] = iScale/zs;
178			rad = AIR_MAX(rad, AIR_ROUNDUP(ik0->kernel->support(ik0->parm)));
179			diam = 2*rad + 1;
180			ixw = AIR_CALLOC(diam, double);
181			iyw = AIR_CALLOC(diam, double);
182			izw = AIR_CALLOC(diam, double);
183			airMopAdd(mop, ixw, airFree, airMopAlways);
184			airMopAdd(mop, iyw, airFree, airMopAlways);
185			airMopAdd(mop, izw, airFree, airMopAlways);
186			if (!(ixw && iyw && izw)) {
187			biffAddf(GAGE, "%s: couldn't allocate grad vector or weight buffers", me);
188			airMopError(mop); return 1;
189			}
190
191			/* the only reason that it is thread-safe to cache gageProbe results,
192			without having the cache hang directly off the gageContext, is that
193			we're doing all the probing for one context in one shot- producing
194			an entirely volume of structure tensors with one function call */
195			gradCache = AIR_CALLOC(3*GAGE_CACHE_LEN, double);
196			coordCache = AIR_CALLOC(3*GAGE_CACHE_LEN, int);
197			airMopAdd(mop, gradCache, airFree, airMopAlways);
198			airMopAdd(mop, coordCache, airFree, airMopAlways);
199			if (!(gradCache && coordCache)) {
200			biffAddf(GAGE, "%s: couldn't allocate caches", me);
201			airMopError(mop); return 1;
202			}
203			for (ixi=0; ixi<GAGE_CACHE_LEN; ixi++) {
204			coordCache[3*ixi + 0] = -1;
205			coordCache[3*ixi + 1] = -1;
206			coordCache[3*ixi + 2] = -1;
207			}
208
209			for (wi=-rad; wi<=rad; wi++) {
210			ik0->parm[0] = iScale/xs;
211			ixw[wi+rad] = ik0->kernel->eval1_d(wi, ik0->parm);
212			ik0->parm[0] = iScale/ys;
213			iyw[wi+rad] = ik0->kernel->eval1_d(wi, ik0->parm);
214			ik0->parm[0] = iScale/zs;
215			izw[wi+rad] = ik0->kernel->eval1_d(wi, ik0->parm);
216			fprintf(stderr, "%d --> (%g,%g,%g) -> (%g,%g,%g)\n",
217			wi, wi/xs, wi/ys, wi/zs, ixw[wi+rad], iyw[wi+rad], izw[wi+rad]);
218			}
219
220			osx = (nin->axis[0].size)/dsmp;
221			osy = (nin->axis[1].size)/dsmp;
222			osz = (nin->axis[2].size)/dsmp;
223			if (nrrdMaybeAlloc_va(nout, nrrdTypeDouble, 4,
224			AIR_CAST(size_t, 7),
225			AIR_CAST(size_t, osx),
226			AIR_CAST(size_t, osy),
227			AIR_CAST(size_t, osz))) {
228			biffAddf(GAGE, NRRD, "%s: couldn't allocate output", me);
229			airMopError(mop); return 1;
230			}
231			airMopAdd(mop, nout, (airMopper)nrrdEmpty, airMopOnError);
232
233			out = AIR_CAST(double *, nout->data);
234			for (ozi=0; ozi<osz; ozi++) {
235			fprintf(stderr, "%s: z = %d/%d\n", me, ozi+1, osz);
236			for (oyi=0; oyi<osy; oyi++) {
237			for (oxi=0; oxi<osx; oxi++) {
238
239			sten[0] = sten[1] = sten[2] = sten[3] = sten[4] = sten[5] = 0;
240			for (_izi=0; _izi<diam; _izi++) {
241			izi = AIR_CLAMP(0, _izi - rad + ozi*dsmp,
242			(int)nin->axis[2].size-1);
243			if (!izw[_izi]) continue;
244			for (_iyi=0; _iyi<diam; _iyi++) {
245			iyi = AIR_CLAMP(0, _iyi - rad + oyi*dsmp,
246			(int)nin->axis[1].size-1);
247			if (!iyw[_iyi]) continue;
248			for (_ixi=0; _ixi<diam; _ixi++) {
249			ixi = AIR_CLAMP(0, _ixi - rad + oxi*dsmp,
250			(int)nin->axis[0].size-1);
251			if (!ixw[_ixi]) continue;
252			wght = ixw[_ixi]iyw[_iyi]izw[_izi];
253			_gageCacheProbe(ctx, grad, coordCache, gradCache, ixi, iyi, izi);
254			sten[0] += wghtgrad[0]grad[0];
255			sten[1] += wghtgrad[0]grad[1];
256			sten[2] += wghtgrad[0]grad[2];
257			sten[3] += wghtgrad[1]grad[1];
258			sten[4] += wghtgrad[1]grad[2];
259			sten[5] += wghtgrad[2]grad[2];
260			}
261			}
262			}
263			out[0] = 1.0;
264			out[1] = sten[0];
265			out[2] = sten[1];
266			out[3] = sten[2];
267			out[4] = sten[3];
268			out[5] = sten[4];
269			out[6] = sten[5];
270			out += 7;
271
272			}
273			}
274			}
275
276			airMopOkay(mop);
277			return 0;
278			}