Head

GCC Code Coverage Report

Directory:	./		Exec	Total	Coverage
File:	src/gage/deconvolve.c	Lines:	0	158	0.0 %
Date:	2017-05-26	Branches:	0	124	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"

int
gageDeconvolve(Nrrd *_nout, double *lastDiffP,
               const Nrrd *nin, const gageKind *kind,
               const NrrdKernelSpec *ksp, int typeOut,
               unsigned int maxIter, int saveAnyway,
               double step, double epsilon, int verbose) {
  static const char me[]="gageDeconvolve";
  gageContext *ctx[2];
  gagePerVolume *pvl[2];
  double *out[2], *val[2], alpha, (*lup)(const void *, size_t), meandiff=0;
  const double *ans[2];
  Nrrd *nout[2];
  airArray *mop;
  unsigned int sx, sy, sz, xi, yi, zi, anslen, thiz=0, last, inIdx, iter;
  int E, valItem;

  if (!(_nout && lastDiffP && nin && kind && ksp)) {
    biffAddf(GAGE, "%s: got NULL pointer", me);
    return 1;
  }
  if (!(nrrdTypeDefault == typeOut
        || !airEnumValCheck(nrrdType, typeOut))) {
    biffAddf(GAGE, "%s: typeOut %d not valid", me, typeOut);
    return 1;
  }
  if (!( maxIter >= 1 )) {
    biffAddf(GAGE, "%s: need maxIter >= 1 (not %u)", me, maxIter);
    return 1;
  }
  if (!( epsilon >= 0 )) {
    biffAddf(GAGE, "%s: need epsilon >= 0.0 (not %g)", me, epsilon);
    return 1;
  }

  /* this once changed from 0 to 1, but is unlikely to change again */
  valItem = 1;

  mop = airMopNew();
  for (iter=0; iter<2; iter++) {
    nout[iter] = nrrdNew();
    airMopAdd(mop, nout[iter], (airMopper)nrrdNuke, airMopAlways);
    if (nrrdConvert(nout[iter], nin, nrrdTypeDouble)) {
      biffMovef(GAGE, NRRD, "%s: couldn't allocate working buffer %u",
                me, iter);
      airMopError(mop); return 1;
    }
    ctx[iter] = gageContextNew();
    airMopAdd(mop, ctx[iter], (airMopper)gageContextNix, airMopAlways);
    E = 0;
    if (!E) E |= !(pvl[iter] = gagePerVolumeNew(ctx[iter], nout[iter], kind));
    if (!E) E |= gagePerVolumeAttach(ctx[iter], pvl[iter]);
    if (!E) E |= gageKernelSet(ctx[iter], gageKernel00,
                               ksp->kernel, ksp->parm);
    if (!E) E |= gageQueryItemOn(ctx[iter], pvl[iter], valItem);
    if (!E) E |= gageUpdate(ctx[iter]);
    if (E) {
      biffAddf(GAGE, "%s: trouble setting up context %u", me, iter);
      airMopError(mop); return 1;
    }
    out[iter] = AIR_CAST(double*, nout[iter]->data);
    ans[iter] = gageAnswerPointer(ctx[iter], pvl[iter], valItem);
  }

  anslen = kind->table[valItem].answerLength;
  lup = nrrdDLookup[nin->type];

  alpha = ksp->kernel->eval1_d(0.0, ksp->parm);
  sx = ctx[0]->shape->size[0];
  sy = ctx[0]->shape->size[1];
  sz = ctx[0]->shape->size[2];

  for (iter=0; iter<maxIter; iter++) {
    thiz = (iter+1) % 2;
    last = (iter+0) % 2;
    val[thiz] = out[thiz];
    val[last] = out[last];
    inIdx = 0;
    meandiff = 0;
    for (zi=0; zi<sz; zi++) {
      for (yi=0; yi<sy; yi++) {
        for (xi=0; xi<sx; xi++) {
          unsigned int ai;
          double in, aa;
          gageProbe(ctx[last], xi, yi, zi);
          for (ai=0; ai<anslen; ai++) {
            in = lup(nin->data, ai + anslen*inIdx);
            aa = ans[last][ai];
            val[thiz][ai] = val[last][ai] + step*(in - aa)/alpha;
            meandiff += 2*(in - aa)*(in - aa)/(DBL_EPSILON + in*in + aa*aa);
          }
          val[thiz] += anslen;
          val[last] += anslen;
          ++inIdx;
        }
      }
    }
    meandiff /= sx*sy*sz;
    if (verbose) {
      fprintf(stderr, "%s: iter %u meandiff = %g\n", me, iter, meandiff);
    }
    if (meandiff < epsilon) {
      /* we have indeed converged while iter < maxIter */
      break;
    }
  }
  if (iter == maxIter) {
    if (!saveAnyway) {
      biffAddf(GAGE, "%s: failed to converge in %u iterations, meandiff = %g",
               me, maxIter, meandiff);
      airMopError(mop); return 1;
    } else {
      if (verbose) {
        fprintf(stderr, "%s: at maxIter %u; err %g still > thresh %g\n", me,
                iter, meandiff, epsilon);
      }
    }
  }

  if (nrrdClampConvert(_nout, nout[thiz], (nrrdTypeDefault == typeOut
                                           ? nin->type
                                           : typeOut))) {
    biffAddf(GAGE, "%s: couldn't create output", me);
    airMopError(mop); return 1;
  }
  *lastDiffP = meandiff;

  airMopOkay(mop);
  return 0;
}

/*
*******************************
** all the following functionality should at some point be
** pushed down to nrrd . . .
*/

static void
deconvLine(double *line, size_t llen, const NrrdKernelSpec *ksp) {

  /* Add as many other parameters to this as you want,
     like number and location of poles, or whatever other
     buffers you think you need (they should be passed,
     not allocated and freed on a per-line basis) */

  /* comment these out when there is a real function body */
  AIR_UNUSED(line);
  AIR_UNUSED(llen);
  AIR_UNUSED(ksp);

  return;
}

static int
deconvTrivial(const NrrdKernelSpec *ksp) {
  int ret;

  /* HEY this will be much easier once kernels have a way of
     advertising whether or not they interpolate */
  if (1 == ksp->parm[0] &&
      (ksp->kernel == nrrdKernelHann ||
       ksp->kernel == nrrdKernelBlackman ||
       ksp->kernel == nrrdKernelBox ||
       ksp->kernel == nrrdKernelCheap ||
       ksp->kernel == nrrdKernelTent)) {
    ret = AIR_TRUE;
  } else {
    ret = AIR_FALSE;
  }
  return ret;
}

int
gageDeconvolveSeparableKnown(const NrrdKernelSpec *ksp) {
  int ret;

  if (!ksp) {
    ret = 0;
  } else if (deconvTrivial(ksp)
             || nrrdKernelBSpline3 == ksp->kernel
             || nrrdKernelBSpline5 == ksp->kernel) {
    ret = 1;
  } else {
    ret = 0;
  }
  return ret;
}

int
gageDeconvolveSeparable(Nrrd *nout, const Nrrd *nin,
                        const gageKind *kind,
                        const NrrdKernelSpec *ksp,
                        int typeOut) {
  static const char me[]="gageDeconvolveSeparable";
  double *line, (*lup)(const void *, size_t),
    (*ins)(void *, size_t, double);
  airArray *mop;
  size_t lineLen, sx, sy, sz, idx, ii, jj;
  unsigned int vi, valLen;

  if (!(nout && nin && kind && ksp)) {
    biffAddf(GAGE, "%s: got NULL pointer", me);
    return 1;
  }
  if (!(nrrdTypeDefault == typeOut
        || !airEnumValCheck(nrrdType, typeOut))) {
    biffAddf(GAGE, "%s: typeOut %d not valid", me, typeOut);
    return 1;
  }
  if (!gageDeconvolveSeparableKnown(ksp)) {
    biffAddf(GAGE, "%s: separable deconv not known for %s kernel",
             me, ksp->kernel->name);
    return 1;
  }
  if (gageKindVolumeCheck(kind, nin)) {
    biffAddf(GAGE, "%s: given volume doesn't fit %s kind",
             me, kind->name);
    return 1;
  }
  if (nrrdTypeDefault == typeOut
      ? nrrdCopy(nout, nin)
      : nrrdConvert(nout, nin, typeOut)) {
    biffMovef(GAGE, NRRD, "%s: problem allocating output", me);
    return 1;
  }
  if (deconvTrivial(ksp)) {
    /* if there's no real work for the deconvolution, then by
       copying the values we're already done; bye */
    return 0;
  }

  valLen = kind->valLen;
  sx = nin->axis[kind->baseDim + 0].size;
  sy = nin->axis[kind->baseDim + 1].size;
  sz = nin->axis[kind->baseDim + 2].size;
  lineLen = sx;
  lineLen = AIR_MAX(lineLen, sy);
  lineLen = AIR_MAX(lineLen, sz);
  lup = nrrdDLookup[nin->type];
  ins = nrrdDInsert[nout->type];

  mop = airMopNew();
  line = AIR_CALLOC(lineLen*valLen, double);
  airMopAdd(mop, line, airFree, airMopAlways);

  /* process along X scanlines */
  for (jj=0; jj<sy*sz; jj++) {
    /* xi = 0, yi = jj%sy, zi = jj/sy
       ==> xi + sx*(yi + sy*zi)
       == 0 + sx*(jj%sy + sy*(jj/sy)) == 0 + sx*jj */
    idx = 0 + valLen*(0 + sx*jj);
    for (ii=0; ii<sx; ii++) {
      for (vi=0; vi<valLen; vi++) {
        line[ii + sx*vi] = lup(nin->data, idx + vi + valLen*ii);
      }
    }
    for (vi=0; vi<valLen; vi++) {
      deconvLine(line + sx*vi, sx, ksp);
    }
    for (ii=0; ii<sx; ii++) {
      for (vi=0; vi<valLen; vi++) {
        ins(nout->data, idx + vi + valLen*ii, line[ii + sx*vi]);
      }
    }
  }

  /* process along Y scanlines */
  for (jj=0; jj<sx*sz; jj++) {
    /* xi = jj%sx, yi = 0, zi = jj/sx
       ==> xi + sx*(yi + sy*zi)
       == jj%sx + sx*(0 + sy*jj/sx) */
    idx = 0 + valLen*((jj%sx) + sx*(0 + sy*(jj/sx)));
    for (ii=0; ii<sy; ii++) {
      for (vi=0; vi<valLen; vi++) {
        line[ii + sy*vi] = lup(nin->data, idx + vi + valLen*sx*ii);
      }
    }
    for (vi=0; vi<valLen; vi++) {
      deconvLine(line + sy*vi, sy, ksp);
    }
    for (ii=0; ii<sx; ii++) {
      for (vi=0; vi<valLen; vi++) {
        ins(nout->data, idx + vi + valLen*sx*ii, line[ii + sy*vi]);
      }
    }
  }

  /* process along Z scanlines */
  for (jj=0; jj<sx*sy; jj++) {
    /* xi = jj%sx, yi = jj/sx, zi = 0
       ==> xi + sx*(yi + sy*zi)
       == jj%sx + sx*(jj/sx + sy*0)
       == jj%sx + sx*(jj/sx) == jj */
    idx = 0 + valLen*jj;
    for (ii=0; ii<sz; ii++) {
      for (vi=0; vi<valLen; vi++) {
        line[ii + sz*vi] = lup(nin->data, idx + vi + valLen*sx*sy*ii);
      }
    }
    for (vi=0; vi<valLen; vi++) {
      deconvLine(line + sz*vi, sz, ksp);
    }
    for (ii=0; ii<sx; ii++) {
      for (vi=0; vi<valLen; vi++) {
        ins(nout->data, idx + vi + valLen*sx*sy*ii, line[ii + sz*vi]);
      }
    }
  }

  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			int
28			gageDeconvolve(Nrrd _nout, double lastDiffP,
29			const Nrrd nin, const gageKind kind,
30			const NrrdKernelSpec *ksp, int typeOut,
31			unsigned int maxIter, int saveAnyway,
32			double step, double epsilon, int verbose) {
33			static const char me[]="gageDeconvolve";
34			gageContext *ctx[2];
35			gagePerVolume *pvl[2];
36			double out[2], val[2], alpha, (lup)(const void , size_t), meandiff=0;
37			const double *ans[2];
38			Nrrd *nout[2];
39			airArray *mop;
40			unsigned int sx, sy, sz, xi, yi, zi, anslen, thiz=0, last, inIdx, iter;
41			int E, valItem;
42
43			if (!(_nout && lastDiffP && nin && kind && ksp)) {
44			biffAddf(GAGE, "%s: got NULL pointer", me);
45			return 1;
46			}
47			if (!(nrrdTypeDefault == typeOut
48			\|\| !airEnumValCheck(nrrdType, typeOut))) {
49			biffAddf(GAGE, "%s: typeOut %d not valid", me, typeOut);
50			return 1;
51			}
52			if (!( maxIter >= 1 )) {
53			biffAddf(GAGE, "%s: need maxIter >= 1 (not %u)", me, maxIter);
54			return 1;
55			}
56			if (!( epsilon >= 0 )) {
57			biffAddf(GAGE, "%s: need epsilon >= 0.0 (not %g)", me, epsilon);
58			return 1;
59			}
60
61			/* this once changed from 0 to 1, but is unlikely to change again */
62			valItem = 1;
63
64			mop = airMopNew();
65			for (iter=0; iter<2; iter++) {
66			nout[iter] = nrrdNew();
67			airMopAdd(mop, nout[iter], (airMopper)nrrdNuke, airMopAlways);
68			if (nrrdConvert(nout[iter], nin, nrrdTypeDouble)) {
69			biffMovef(GAGE, NRRD, "%s: couldn't allocate working buffer %u",
70			me, iter);
71			airMopError(mop); return 1;
72			}
73			ctx[iter] = gageContextNew();
74			airMopAdd(mop, ctx[iter], (airMopper)gageContextNix, airMopAlways);
75			E = 0;
76			if (!E) E \|= !(pvl[iter] = gagePerVolumeNew(ctx[iter], nout[iter], kind));
77			if (!E) E \|= gagePerVolumeAttach(ctx[iter], pvl[iter]);
78			if (!E) E \|= gageKernelSet(ctx[iter], gageKernel00,
79			ksp->kernel, ksp->parm);
80			if (!E) E \|= gageQueryItemOn(ctx[iter], pvl[iter], valItem);
81			if (!E) E \|= gageUpdate(ctx[iter]);
82			if (E) {
83			biffAddf(GAGE, "%s: trouble setting up context %u", me, iter);
84			airMopError(mop); return 1;
85			}
86			out[iter] = AIR_CAST(double*, nout[iter]->data);
87			ans[iter] = gageAnswerPointer(ctx[iter], pvl[iter], valItem);
88			}
89
90			anslen = kind->table[valItem].answerLength;
91			lup = nrrdDLookup[nin->type];
92
93			alpha = ksp->kernel->eval1_d(0.0, ksp->parm);
94			sx = ctx[0]->shape->size[0];
95			sy = ctx[0]->shape->size[1];
96			sz = ctx[0]->shape->size[2];
97
98			for (iter=0; iter<maxIter; iter++) {
99			thiz = (iter+1) % 2;
100			last = (iter+0) % 2;
101			val[thiz] = out[thiz];
102			val[last] = out[last];
103			inIdx = 0;
104			meandiff = 0;
105			for (zi=0; zi<sz; zi++) {
106			for (yi=0; yi<sy; yi++) {
107			for (xi=0; xi<sx; xi++) {
108			unsigned int ai;
109			double in, aa;
110			gageProbe(ctx[last], xi, yi, zi);
111			for (ai=0; ai<anslen; ai++) {
112			in = lup(nin->data, ai + anslen*inIdx);
113			aa = ans[last][ai];
114			val[thiz][ai] = val[last][ai] + step*(in - aa)/alpha;
115			meandiff += 2(in - aa)(in - aa)/(DBL_EPSILON + inin + aaaa);
116			}
117			val[thiz] += anslen;
118			val[last] += anslen;
119			++inIdx;
120			}
121			}
122			}
123			meandiff /= sxsysz;
124			if (verbose) {
125			fprintf(stderr, "%s: iter %u meandiff = %g\n", me, iter, meandiff);
126			}
127			if (meandiff < epsilon) {
128			/* we have indeed converged while iter < maxIter */
129			break;
130			}
131			}
132			if (iter == maxIter) {
133			if (!saveAnyway) {
134			biffAddf(GAGE, "%s: failed to converge in %u iterations, meandiff = %g",
135			me, maxIter, meandiff);
136			airMopError(mop); return 1;
137			} else {
138			if (verbose) {
139			fprintf(stderr, "%s: at maxIter %u; err %g still > thresh %g\n", me,
140			iter, meandiff, epsilon);
141			}
142			}
143			}
144
145			if (nrrdClampConvert(_nout, nout[thiz], (nrrdTypeDefault == typeOut
146			? nin->type
147			: typeOut))) {
148			biffAddf(GAGE, "%s: couldn't create output", me);
149			airMopError(mop); return 1;
150			}
151			*lastDiffP = meandiff;
152
153			airMopOkay(mop);
154			return 0;
155			}
156
157			/*
158			*******************************
159			** all the following functionality should at some point be
160			** pushed down to nrrd . . .
161			*/
162
163			static void
164			deconvLine(double line, size_t llen, const NrrdKernelSpec ksp) {
165
166			/* Add as many other parameters to this as you want,
167			like number and location of poles, or whatever other
168			buffers you think you need (they should be passed,
169			not allocated and freed on a per-line basis) */
170
171			/* comment these out when there is a real function body */
172			AIR_UNUSED(line);
173			AIR_UNUSED(llen);
174			AIR_UNUSED(ksp);
175
176			return;
177			}
178
179			static int
180			deconvTrivial(const NrrdKernelSpec *ksp) {
181			int ret;
182
183			/* HEY this will be much easier once kernels have a way of
184			advertising whether or not they interpolate */
185			if (1 == ksp->parm[0] &&
186			(ksp->kernel == nrrdKernelHann \|\|
187			ksp->kernel == nrrdKernelBlackman \|\|
188			ksp->kernel == nrrdKernelBox \|\|
189			ksp->kernel == nrrdKernelCheap \|\|
190			ksp->kernel == nrrdKernelTent)) {
191			ret = AIR_TRUE;
192			} else {
193			ret = AIR_FALSE;
194			}
195			return ret;
196			}
197
198			int
199			gageDeconvolveSeparableKnown(const NrrdKernelSpec *ksp) {
200			int ret;
201
202			if (!ksp) {
203			ret = 0;
204			} else if (deconvTrivial(ksp)
205			\|\| nrrdKernelBSpline3 == ksp->kernel
206			\|\| nrrdKernelBSpline5 == ksp->kernel) {
207			ret = 1;
208			} else {
209			ret = 0;
210			}
211			return ret;
212			}
213
214			int
215			gageDeconvolveSeparable(Nrrd nout, const Nrrd nin,
216			const gageKind *kind,
217			const NrrdKernelSpec *ksp,
218			int typeOut) {
219			static const char me[]="gageDeconvolveSeparable";
220			double line, (lup)(const void *, size_t),
221			(ins)(void , size_t, double);
222			airArray *mop;
223			size_t lineLen, sx, sy, sz, idx, ii, jj;
224			unsigned int vi, valLen;
225
226			if (!(nout && nin && kind && ksp)) {
227			biffAddf(GAGE, "%s: got NULL pointer", me);
228			return 1;
229			}
230			if (!(nrrdTypeDefault == typeOut
231			\|\| !airEnumValCheck(nrrdType, typeOut))) {
232			biffAddf(GAGE, "%s: typeOut %d not valid", me, typeOut);
233			return 1;
234			}
235			if (!gageDeconvolveSeparableKnown(ksp)) {
236			biffAddf(GAGE, "%s: separable deconv not known for %s kernel",
237			me, ksp->kernel->name);
238			return 1;
239			}
240			if (gageKindVolumeCheck(kind, nin)) {
241			biffAddf(GAGE, "%s: given volume doesn't fit %s kind",
242			me, kind->name);
243			return 1;
244			}
245			if (nrrdTypeDefault == typeOut
246			? nrrdCopy(nout, nin)
247			: nrrdConvert(nout, nin, typeOut)) {
248			biffMovef(GAGE, NRRD, "%s: problem allocating output", me);
249			return 1;
250			}
251			if (deconvTrivial(ksp)) {
252			/* if there's no real work for the deconvolution, then by
253			copying the values we're already done; bye */
254			return 0;
255			}
256
257			valLen = kind->valLen;
258			sx = nin->axis[kind->baseDim + 0].size;
259			sy = nin->axis[kind->baseDim + 1].size;
260			sz = nin->axis[kind->baseDim + 2].size;
261			lineLen = sx;
262			lineLen = AIR_MAX(lineLen, sy);
263			lineLen = AIR_MAX(lineLen, sz);
264			lup = nrrdDLookup[nin->type];
265			ins = nrrdDInsert[nout->type];
266
267			mop = airMopNew();
268			line = AIR_CALLOC(lineLen*valLen, double);
269			airMopAdd(mop, line, airFree, airMopAlways);
270
271			/* process along X scanlines */
272			for (jj=0; jj<sy*sz; jj++) {
273			/* xi = 0, yi = jj%sy, zi = jj/sy
274			==> xi + sx(yi + syzi)
275			== 0 + sx(jj%sy + sy(jj/sy)) == 0 + sxjj /
276			idx = 0 + valLen(0 + sxjj);
277			for (ii=0; ii<sx; ii++) {
278			for (vi=0; vi<valLen; vi++) {
279			line[ii + sxvi] = lup(nin->data, idx + vi + valLenii);
280			}
281			}
282			for (vi=0; vi<valLen; vi++) {
283			deconvLine(line + sx*vi, sx, ksp);
284			}
285			for (ii=0; ii<sx; ii++) {
286			for (vi=0; vi<valLen; vi++) {
287			ins(nout->data, idx + vi + valLenii, line[ii + sxvi]);
288			}
289			}
290			}
291
292			/* process along Y scanlines */
293			for (jj=0; jj<sx*sz; jj++) {
294			/* xi = jj%sx, yi = 0, zi = jj/sx
295			==> xi + sx(yi + syzi)
296			== jj%sx + sx(0 + syjj/sx) */
297			idx = 0 + valLen((jj%sx) + sx(0 + sy*(jj/sx)));
298			for (ii=0; ii<sy; ii++) {
299			for (vi=0; vi<valLen; vi++) {
300			line[ii + syvi] = lup(nin->data, idx + vi + valLensx*ii);
301			}
302			}
303			for (vi=0; vi<valLen; vi++) {
304			deconvLine(line + sy*vi, sy, ksp);
305			}
306			for (ii=0; ii<sx; ii++) {
307			for (vi=0; vi<valLen; vi++) {
308			ins(nout->data, idx + vi + valLensxii, line[ii + sy*vi]);
309			}
310			}
311			}
312
313			/* process along Z scanlines */
314			for (jj=0; jj<sx*sy; jj++) {
315			/* xi = jj%sx, yi = jj/sx, zi = 0
316			==> xi + sx(yi + syzi)
317			== jj%sx + sx(jj/sx + sy0)
318			== jj%sx + sx(jj/sx) == jj /
319			idx = 0 + valLen*jj;
320			for (ii=0; ii<sz; ii++) {
321			for (vi=0; vi<valLen; vi++) {
322			line[ii + szvi] = lup(nin->data, idx + vi + valLensxsyii);
323			}
324			}
325			for (vi=0; vi<valLen; vi++) {
326			deconvLine(line + sz*vi, sz, ksp);
327			}
328			for (ii=0; ii<sx; ii++) {
329			for (vi=0; vi<valLen; vi++) {
330			ins(nout->data, idx + vi + valLensxsyii, line[ii + szvi]);
331			}
332			}
333			}
334
335			airMopOkay(mop);
336			return 0;
337			}