Head

GCC Code Coverage Report

Directory:	./		Exec	Total	Coverage
File:	src/gage/sclfilter.c	Lines:	225	229	98.3 %
Date:	2017-05-26	Branches:	216	252	85.7 %


/*
  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"

/*
** The functions gageScl3PFilter2, gageScl3PFilter4, and
** gageScl3PFilterN weren't really intended to be public, but they
** have become so because they implement the brains of the lowest
** level of dot-products needed to do the convolution-based
** reconstruction, and, to do the the transforms that affect gradient
** and hessian values.  These functions are very useful for the
** implementation of other gageKinds, in which non-scalar values are
** filtered per-component.
*/

/*
** Teem trivia: "3P" stands for "3-pack", the kind of kernel set used (k00,
** k11, and k22).  6-pack filtering is still unimplemented. The name
** "3-pack" persists even with the possibility of higher derivatives, and
** what would be called 10-pack filtering (ability to specify and use all
** 10 kinds of spatial kernels) is also still unimplemented.
*/

#define X 0
#define Y 1
#define Z 2

void
gageScl3PFilter2(gageShape *shape,
                 double *ivX, double *ivY, double *ivZ,
                 double *fw0, double *fw1, double *fw2,
                 double *val, double *gvec, double *hess,
                 const int *needD) {
  int doV, doD1, doD2;
  doV = needD[0];
  doD1 = needD[1];
  doD2 = needD[2];

  /* fw? + 2*?
       |     |
       |     +- along which axis (0:x, 1:y, 2:z)
       |
       + what information (0:value, 1:1st deriv, 2:2nd deriv)

     ivX: 3D cube cache of original volume values
          (its scanlines are along the X axis)
     ivY: 2D square cache of intermediate filter results
          (its scanlines are along the Y axis)
     ivZ: 1D linear cache of intermediate filter results
          (it is a scanline along the Z axis)
  */

#define DOT_2(a, b) ((a)[0]*(b)[0] + (a)[1]*(b)[1])
#define VL_2(i, axis) DOT_2(fw0 + (axis)*2, iv##axis + i*2)
#define D1_2(i, axis) DOT_2(fw1 + (axis)*2, iv##axis + i*2)
#define D2_2(i, axis) DOT_2(fw2 + (axis)*2, iv##axis + i*2)

  /* x0 */
  ivY[0] = VL_2(0,X); /* interpolate values of 0th scanline along X axis */
  ivY[1] = VL_2(1,X);
  ivY[2] = VL_2(2,X);
  ivY[3] = VL_2(3,X);
  /* x0y0 */
  ivZ[0] = VL_2(0,Y);
  ivZ[1] = VL_2(1,Y);
  /* x0y0z0 */
  if (doV) {
    *val = VL_2(0,Z);                          /* f */
  }

  if (!( doD1 || doD2 ))
    return;

  /* x0y0z1 */
  if (doD1) {
    gvec[2] = D1_2(0,Z);                       /* g_z */
  }
  if (doD2) {
    /* actually, there is no possible way in which it makes sense to
       try to measure a second derivative with only two samples, so
       all this "if (doD2)" code is basically bogus, but we'll keep it
       around for generality . . . */
    /* x0y0z2 */
    hess[8] = D2_2(0,Z);                       /* h_zz */
  }
  /* x0y1 */
  ivZ[0] = D1_2(0,Y);
  ivZ[1] = D1_2(1,Y);
  /* x0y1z0 */
  if (doD1) {
    gvec[1] = VL_2(0,Z);                       /* g_y */
  }
  if (doD2) {
    /* x0y1z1 */
    hess[5] = hess[7] = D1_2(0,Z);             /* h_yz */
    /* x0y2 */
    ivZ[0] = D2_2(0,Y);
    ivZ[1] = D2_2(1,Y);
    /* x0y2z0 */
    hess[4] = VL_2(0,Z);                       /* h_yy */
  }
  /* x1 */
  ivY[0] = D1_2(0,X);
  ivY[1] = D1_2(1,X);
  ivY[2] = D1_2(2,X);
  ivY[3] = D1_2(3,X);
  /* x1y0 */
  ivZ[0] = VL_2(0,Y);
  ivZ[1] = VL_2(1,Y);
  /* x1y0z0 */
  if (doD1) {
    gvec[0] = VL_2(0,Z);                       /* g_x */
  }

  ell_3mv_mul_d(gvec, shape->ItoWSubInvTransp, gvec);

  if (!doD2)
    return;

  /* x1y0z1 */
  hess[2] = hess[6] = D1_2(0,Z);               /* h_xz */
  /* x1y1 */
  ivZ[0] = D1_2(0,Y);
  ivZ[1] = D1_2(1,Y);
  /* x1y1z0 */
  hess[1] = hess[3] = VL_2(0,Z);               /* h_xy */
  /* x2 */
  ivY[0] = D2_2(0,X);
  ivY[1] = D2_2(1,X);
  ivY[2] = D2_2(2,X);
  ivY[3] = D2_2(3,X);
  /* x2y0 */
  ivZ[0] = VL_2(0,Y);
  ivZ[1] = VL_2(1,Y);
  /* x2y0z0 */
  hess[0] = VL_2(0,Z);                         /* h_xx */

  if (1) {
    double matA[9];
    ELL_3M_MUL(matA, shape->ItoWSubInvTransp, hess);
    ELL_3M_MUL(hess, matA, shape->ItoWSubInv);
  }

  return;
}

void
gageScl3PFilter4(gageShape *shape,
                 double *ivX, double *ivY, double *ivZ,
                 double *fw0, double *fw1, double *fw2,
                 double *val, double *gvec, double *hess,
                 const int *needD) {
  int doV, doD1, doD2;
  doV = needD[0];
  doD1 = needD[1];
  doD2 = needD[2];

  /* fw? + 4*?
       |     |
       |     +- along which axis (0:x, 1:y, 2:z)
       |
       + what information (0:value, 1:1st deriv, 2:2nd deriv)

     ivX: 3D cube cache of original volume values
          (its scanlines are along the X axis)
     ivY: 2D square cache of intermediate filter results
          (its scanlines are along the Y axis)
     ivZ: 1D linear cache of intermediate filter results
          (it is a scanline along the Z axis)
  */

#define DOT_4(a,b) ((a)[0]*(b)[0]+(a)[1]*(b)[1]+(a)[2]*(b)[2]+(a)[3]*(b)[3])
#define VL_4(i, axis) DOT_4(fw0 + (axis)*4, iv##axis + i*4)
#define D1_4(i, axis) DOT_4(fw1 + (axis)*4, iv##axis + i*4)
#define D2_4(i, axis) DOT_4(fw2 + (axis)*4, iv##axis + i*4)

  /* x0 */
  ivY[ 0] = VL_4( 0,X);
  ivY[ 1] = VL_4( 1,X);
  ivY[ 2] = VL_4( 2,X);
  ivY[ 3] = VL_4( 3,X);
  ivY[ 4] = VL_4( 4,X);
  ivY[ 5] = VL_4( 5,X);
  ivY[ 6] = VL_4( 6,X);
  ivY[ 7] = VL_4( 7,X);
  ivY[ 8] = VL_4( 8,X);
  ivY[ 9] = VL_4( 9,X);
  ivY[10] = VL_4(10,X);
  ivY[11] = VL_4(11,X);
  ivY[12] = VL_4(12,X);
  ivY[13] = VL_4(13,X);
  ivY[14] = VL_4(14,X);
  ivY[15] = VL_4(15,X);
  /*
  */
  /* x0y0 */
  ivZ[ 0] = VL_4( 0,Y);
  ivZ[ 1] = VL_4( 1,Y);
  ivZ[ 2] = VL_4( 2,Y);
  ivZ[ 3] = VL_4( 3,Y);
  /* x0y0z0 */
  if (doV) {
    *val = VL_4( 0,Z);                          /* f */
  }

  if (!( doD1 || doD2 ))
    return;

  /* x0y0z1 */
  if (doD1) {
    gvec[2] = D1_4( 0,Z);                       /* g_z */
  }
  if (doD2) {
    /* x0y0z2 */
    hess[8] = D2_4( 0,Z);                       /* h_zz */
  }
  /* x0y1 */
  ivZ[ 0] = D1_4( 0,Y);
  ivZ[ 1] = D1_4( 1,Y);
  ivZ[ 2] = D1_4( 2,Y);
  ivZ[ 3] = D1_4( 3,Y);
  /* x0y1z0 */
  if (doD1) {
    gvec[1] = VL_4( 0,Z);                       /* g_y */
  }
  if (doD2) {
    /* x0y1z1 */
    hess[5] = hess[7] = D1_4( 0,Z);             /* h_yz */
    /* x0y2 */
    ivZ[ 0] = D2_4( 0,Y);
    ivZ[ 1] = D2_4( 1,Y);
    ivZ[ 2] = D2_4( 2,Y);
    ivZ[ 3] = D2_4( 3,Y);
    /* x0y2z0 */
    hess[4] = VL_4( 0,Z);                       /* h_yy */
  }
  /* x1 */
  ivY[ 0] = D1_4( 0,X);
  ivY[ 1] = D1_4( 1,X);
  ivY[ 2] = D1_4( 2,X);
  ivY[ 3] = D1_4( 3,X);
  ivY[ 4] = D1_4( 4,X);
  ivY[ 5] = D1_4( 5,X);
  ivY[ 6] = D1_4( 6,X);
  ivY[ 7] = D1_4( 7,X);
  ivY[ 8] = D1_4( 8,X);
  ivY[ 9] = D1_4( 9,X);
  ivY[10] = D1_4(10,X);
  ivY[11] = D1_4(11,X);
  ivY[12] = D1_4(12,X);
  ivY[13] = D1_4(13,X);
  ivY[14] = D1_4(14,X);
  ivY[15] = D1_4(15,X);
  /* x1y0 */
  ivZ[ 0] = VL_4( 0,Y);
  ivZ[ 1] = VL_4( 1,Y);
  ivZ[ 2] = VL_4( 2,Y);
  ivZ[ 3] = VL_4( 3,Y);
  /* x1y0z0 */
  if (doD1) {
    gvec[0] = VL_4( 0,Z);                       /* g_x */
  }

  ell_3mv_mul_d(gvec, shape->ItoWSubInvTransp, gvec);

  if (!doD2)
    return;

  /* x1y0z1 */
  hess[2] = hess[6] = D1_4( 0,Z);               /* h_xz */
  /* x1y1 */
  ivZ[ 0] = D1_4( 0,Y);
  ivZ[ 1] = D1_4( 1,Y);
  ivZ[ 2] = D1_4( 2,Y);
  ivZ[ 3] = D1_4( 3,Y);
  /* x1y1z0 */
  hess[1] = hess[3] = VL_4( 0,Z);               /* h_xy */
  /* x2 */
  ivY[ 0] = D2_4( 0,X);
  ivY[ 1] = D2_4( 1,X);
  ivY[ 2] = D2_4( 2,X);
  ivY[ 3] = D2_4( 3,X);
  ivY[ 4] = D2_4( 4,X);
  ivY[ 5] = D2_4( 5,X);
  ivY[ 6] = D2_4( 6,X);
  ivY[ 7] = D2_4( 7,X);
  ivY[ 8] = D2_4( 8,X);
  ivY[ 9] = D2_4( 9,X);
  ivY[10] = D2_4(10,X);
  ivY[11] = D2_4(11,X);
  ivY[12] = D2_4(12,X);
  ivY[13] = D2_4(13,X);
  ivY[14] = D2_4(14,X);
  ivY[15] = D2_4(15,X);
  /* x2y0 */
  ivZ[ 0] = VL_4( 0,Y);
  ivZ[ 1] = VL_4( 1,Y);
  ivZ[ 2] = VL_4( 2,Y);
  ivZ[ 3] = VL_4( 3,Y);
  /* x2y0z0 */
  hess[0] = VL_4( 0,Z);                         /* h_xx */

  if (1) {
    double matA[9];
    ELL_3M_MUL(matA, shape->ItoWSubInvTransp, hess);
    ELL_3M_MUL(hess, matA, shape->ItoWSubInv);
  }

  return;
}

void
gageScl3PFilter6(gageShape *shape,
                 double *ivX, double *ivY, double *ivZ,
                 double *fw0, double *fw1, double *fw2,
                 double *val, double *gvec, double *hess,
                 const int *needD) {
  int i, j;
  double T;
  int doV, doD1, doD2;
  doV = needD[0];
  doD1 = needD[1];
  doD2 = needD[2];

#define fd 6
#include "scl3pfilterbody.c"
#undef fd

  return;
}

void
gageScl3PFilter8(gageShape *shape,
                 double *ivX, double *ivY, double *ivZ,
                 double *fw0, double *fw1, double *fw2,
                 double *val, double *gvec, double *hess,
                 const int *needD) {
  int i, j;
  double T;
  int doV, doD1, doD2;
  doV = needD[0];
  doD1 = needD[1];
  doD2 = needD[2];

#define fd 8
#include "scl3pfilterbody.c"
#undef fd

  return;
}

void
gageScl3PFilterN(gageShape *shape, int fd,
                 double *ivX, double *ivY, double *ivZ,
                 double *fw0, double *fw1, double *fw2,
                 double *val, double *gvec, double *hess,
                 const int *needD) {
  int i, j;
  double T;
  int doV, doD1, doD2;
  doV = needD[0];
  doD1 = needD[1];
  doD2 = needD[2];

#include "scl3pfilterbody.c"

  return;
}

void
_gageSclFilter(gageContext *ctx, gagePerVolume *pvl) {
  char me[]="_gageSclFilter";
  int fd;
  double *fw00, *fw11, *fw22;
  gageScl3PFilter_t *filter[5] = {NULL, gageScl3PFilter2, gageScl3PFilter4,
                                  gageScl3PFilter6, gageScl3PFilter8};

  fd = 2*ctx->radius;
  if (!ctx->parm.k3pack) {
    fprintf(stderr, "!%s: sorry, 6-pack filtering not implemented\n", me);
    return;
  }
  fw00 = ctx->fw + fd*3*gageKernel00;
  fw11 = ctx->fw + fd*3*gageKernel11;
  fw22 = ctx->fw + fd*3*gageKernel22;
  /* perform the filtering */
  if (fd <= 8) {
    filter[ctx->radius](ctx->shape, pvl->iv3, pvl->iv2, pvl->iv1,
                        fw00, fw11, fw22,
                        pvl->directAnswer[gageSclValue],
                        pvl->directAnswer[gageSclGradVec],
                        pvl->directAnswer[gageSclHessian],
                        pvl->needD);
  } else {
    gageScl3PFilterN(ctx->shape, fd,
                     pvl->iv3, pvl->iv2, pvl->iv1,
                     fw00, fw11, fw22,
                     pvl->directAnswer[gageSclValue],
                     pvl->directAnswer[gageSclGradVec],
                     pvl->directAnswer[gageSclHessian],
                     pvl->needD);
  }

  return;
}

#undef X
#undef Y
#undef Z


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			/*
28			** The functions gageScl3PFilter2, gageScl3PFilter4, and
29			** gageScl3PFilterN weren't really intended to be public, but they
30			** have become so because they implement the brains of the lowest
31			** level of dot-products needed to do the convolution-based
32			** reconstruction, and, to do the the transforms that affect gradient
33			** and hessian values. These functions are very useful for the
34			** implementation of other gageKinds, in which non-scalar values are
35			** filtered per-component.
36			*/
37
38			/*
39			** Teem trivia: "3P" stands for "3-pack", the kind of kernel set used (k00,
40			** k11, and k22). 6-pack filtering is still unimplemented. The name
41			** "3-pack" persists even with the possibility of higher derivatives, and
42			** what would be called 10-pack filtering (ability to specify and use all
43			** 10 kinds of spatial kernels) is also still unimplemented.
44			*/
45
46			#define X 0
47			#define Y 1
48			#define Z 2
49
50			void
51	✓✓✓✓ ✓✓✓✓ ✓✓✓✓	2654223002	gageScl3PFilter2(gageShape *shape,
52			double ivX, double ivY, double *ivZ,
53	✓✓✓✓ ✓✓✓✓ ✓✓✓✓	167472314	double fw0, double fw1, double *fw2,
54			double val, double gvec, double *hess,
55	✓✗✓✗ ✓✗	392749	const int *needD) {
56	✓✓✓✓ ✓✓	11196483	int doV, doD1, doD2;
57		1690125	doV = needD[0];
58		648688	doD1 = needD[1];
59	✓✓✓✓ ✓✓	1041437	doD2 = needD[2];
60		193200
61			/* fw? + 2*?
62			\| \|
63	✓✗✓✗ ✓✗	199549	\| +- along which axis (0:x, 1:y, 2:z)
64	✓✓✓✓ ✓✓	5659683	\|
65		199549	+ what information (0:value, 1:1st deriv, 2:2nd deriv)
66	✓✗✓✗ ✓✗	199549
67			ivX: 3D cube cache of original volume values
68	✓✓✓✓ ✓✓	5659683	(its scanlines are along the X axis)
69		199549	ivY: 2D square cache of intermediate filter results
70			(its scanlines are along the Y axis)
71	✓✓✓✓ ✓✓✓✓ ✓✓✓✓	84309914	ivZ: 1D linear cache of intermediate filter results
72			(it is a scanline along the Z axis)
73	✓✗✓✗ ✓✗	199549	*/
74	✓✓✓✓ ✓✓	5659683
75		199549	#define DOT_2(a, b) ((a)[0](b)[0] + (a)[1](b)[1])
76	✓✗✓✗ ✓✗	199549	#define VL_2(i, axis) DOT_2(fw0 + (axis)2, iv##axis + i2)
77			#define D1_2(i, axis) DOT_2(fw1 + (axis)2, iv##axis + i2)
78	✓✓✓✓ ✓✓	5659683	#define D2_2(i, axis) DOT_2(fw2 + (axis)2, iv##axis + i2)
79		199549
80			/* x0 */
81	✓✓✓✓ ✓✓✓✓ ✓✓✓✓	82428084	ivY[0] = VL_2(0,X); /* interpolate values of 0th scanline along X axis */
82		648688	ivY[1] = VL_2(1,X);
83	✓✓✓✓ ✓✓	6108822	ivY[2] = VL_2(2,X);
84		848237	ivY[3] = VL_2(3,X);
85			/* x0y0 */
86	✓✓✓✓ ✓✓✓✓ ✓✓✓✓	1333189290	ivZ[0] = VL_2(0,Y);
87		648688	ivZ[1] = VL_2(1,Y);
88	✓✓✓✓ ✓✓✓✓ ✓✓✓✓	84309914	/* x0y0z0 */
89	✓✗	648688	if (doV) {
90	✓✗✓✗ ✓✗	848237	val = VL_2(0,Z); / f */
91	✓✓✓✓ ✓✓	6308371	}
92		199549
93	✓✓	648688	if (!( doD1 \|\| doD2 ))
94		345509	return;
95
96	✗✓✗✓ ✗✓	199549	/* x0y0z1 */
97	✓✗	502728	if (doD1) {
98		502728	gvec[2] = D1_2(0,Z); /* g_z */
99		502728	}
100	✓✗✓✓ ✓✓✓✓	6162411	if (doD2) {
101		199549	/* actually, there is no possible way in which it makes sense to
102			try to measure a second derivative with only two samples, so
103	✓✓✓✓ ✓✓✓✓ ✓✓✓✓	81779396	all this "if (doD2)" code is basically bogus, but we'll keep it
104			around for generality . . . */
105	✓✓✓✓ ✓✓	5460134	/* x0y0z2 */
106		702277	hess[8] = D2_2(0,Z); /* h_zz */
107		502728	}
108	✓✓✓✓ ✓✓✓✓ ✓✓✓✓	1332540602	/* x0y1 */
109		502728	ivZ[0] = D1_2(0,Y);
110	✓✓✓✓ ✓✓✓✓ ✓✓✓✓	82282124	ivZ[1] = D1_2(1,Y);
111			/* x0y1z0 */
112	✓✗✓✓ ✓✓✓✓	5962862	if (doD1) {
113		502728	gvec[1] = VL_2(0,Z); /* g_y */
114		502728	}
115	✓✗	502728	if (doD2) {
116		199549	/* x0y1z1 */
117		702277	hess[5] = hess[7] = D1_2(0,Z); /* h_yz */
118			/* x0y2 */
119		502728	ivZ[0] = D2_2(0,Y);
120		502728	ivZ[1] = D2_2(1,Y);
121			/* x0y2z0 */
122		502728	hess[4] = VL_2(0,Z); /* h_yy */
123		502728	}
124			/* x1 */
125		502728	ivY[0] = D1_2(0,X);
126		502728	ivY[1] = D1_2(1,X);
127		502728	ivY[2] = D1_2(2,X);
128		502728	ivY[3] = D1_2(3,X);
129			/* x1y0 */
130		502728	ivZ[0] = VL_2(0,Y);
131		502728	ivZ[1] = VL_2(1,Y);
132			/* x1y0z0 */
133	✓✗	502728	if (doD1) {
134		502728	gvec[0] = VL_2(0,Z); /* g_x */
135		502728	}
136
137		502728	ell_3mv_mul_d(gvec, shape->ItoWSubInvTransp, gvec);
138
139	✗✓	502728	if (!doD2)
140			return;
141
142			/* x1y0z1 */
143		502728	hess[2] = hess[6] = D1_2(0,Z); /* h_xz */
144			/* x1y1 */
145		502728	ivZ[0] = D1_2(0,Y);
146		502728	ivZ[1] = D1_2(1,Y);
147			/* x1y1z0 */
148		502728	hess[1] = hess[3] = VL_2(0,Z); /* h_xy */
149			/* x2 */
150		502728	ivY[0] = D2_2(0,X);
151		502728	ivY[1] = D2_2(1,X);
152		502728	ivY[2] = D2_2(2,X);
153		502728	ivY[3] = D2_2(3,X);
154			/* x2y0 */
155		502728	ivZ[0] = VL_2(0,Y);
156		502728	ivZ[1] = VL_2(1,Y);
157			/* x2y0z0 */
158		502728	hess[0] = VL_2(0,Z); /* h_xx */
159
160			if (1) {
161			double matA[9];
162		502728	ELL_3M_MUL(matA, shape->ItoWSubInvTransp, hess);
163		502728	ELL_3M_MUL(hess, matA, shape->ItoWSubInv);
164			}
165
166		502728	return;
167		648688	}
168
169			void
170			gageScl3PFilter4(gageShape *shape,
171			double ivX, double ivY, double *ivZ,
172			double fw0, double fw1, double *fw2,
173			double val, double gvec, double *hess,
174			const int *needD) {
175			int doV, doD1, doD2;
176		954176	doV = needD[0];
177		477088	doD1 = needD[1];
178		477088	doD2 = needD[2];
179
180			/* fw? + 4*?
181			\| \|
182			\| +- along which axis (0:x, 1:y, 2:z)
183			\|
184			+ what information (0:value, 1:1st deriv, 2:2nd deriv)
185
186			ivX: 3D cube cache of original volume values
187			(its scanlines are along the X axis)
188			ivY: 2D square cache of intermediate filter results
189			(its scanlines are along the Y axis)
190			ivZ: 1D linear cache of intermediate filter results
191			(it is a scanline along the Z axis)
192			*/
193
194			#define DOT_4(a,b) ((a)[0](b)[0]+(a)[1](b)[1]+(a)[2](b)[2]+(a)[3](b)[3])
195			#define VL_4(i, axis) DOT_4(fw0 + (axis)4, iv##axis + i4)
196			#define D1_4(i, axis) DOT_4(fw1 + (axis)4, iv##axis + i4)
197			#define D2_4(i, axis) DOT_4(fw2 + (axis)4, iv##axis + i4)
198
199			/* x0 */
200		477088	ivY[ 0] = VL_4( 0,X);
201		477088	ivY[ 1] = VL_4( 1,X);
202		477088	ivY[ 2] = VL_4( 2,X);
203		477088	ivY[ 3] = VL_4( 3,X);
204		477088	ivY[ 4] = VL_4( 4,X);
205		477088	ivY[ 5] = VL_4( 5,X);
206		477088	ivY[ 6] = VL_4( 6,X);
207		477088	ivY[ 7] = VL_4( 7,X);
208		477088	ivY[ 8] = VL_4( 8,X);
209		477088	ivY[ 9] = VL_4( 9,X);
210		477088	ivY[10] = VL_4(10,X);
211		477088	ivY[11] = VL_4(11,X);
212		477088	ivY[12] = VL_4(12,X);
213		477088	ivY[13] = VL_4(13,X);
214		477088	ivY[14] = VL_4(14,X);
215		477088	ivY[15] = VL_4(15,X);
216			/*
217			*/
218			/* x0y0 */
219		477088	ivZ[ 0] = VL_4( 0,Y);
220		477088	ivZ[ 1] = VL_4( 1,Y);
221		477088	ivZ[ 2] = VL_4( 2,Y);
222		477088	ivZ[ 3] = VL_4( 3,Y);
223			/* x0y0z0 */
224	✓✗	477088	if (doV) {
225		477088	val = VL_4( 0,Z); / f */
226		477088	}
227
228	✓✓	477088	if (!( doD1 \|\| doD2 ))
229		129460	return;
230
231			/* x0y0z1 */
232	✓✗	347628	if (doD1) {
233		347628	gvec[2] = D1_4( 0,Z); /* g_z */
234		347628	}
235	✓✗	347628	if (doD2) {
236			/* x0y0z2 */
237		347628	hess[8] = D2_4( 0,Z); /* h_zz */
238		347628	}
239			/* x0y1 */
240		347628	ivZ[ 0] = D1_4( 0,Y);
241		347628	ivZ[ 1] = D1_4( 1,Y);
242		347628	ivZ[ 2] = D1_4( 2,Y);
243		347628	ivZ[ 3] = D1_4( 3,Y);
244			/* x0y1z0 */
245	✓✗	347628	if (doD1) {
246		347628	gvec[1] = VL_4( 0,Z); /* g_y */
247		347628	}
248	✓✗	347628	if (doD2) {
249			/* x0y1z1 */
250		347628	hess[5] = hess[7] = D1_4( 0,Z); /* h_yz */
251			/* x0y2 */
252		347628	ivZ[ 0] = D2_4( 0,Y);
253		347628	ivZ[ 1] = D2_4( 1,Y);
254		347628	ivZ[ 2] = D2_4( 2,Y);
255		347628	ivZ[ 3] = D2_4( 3,Y);
256			/* x0y2z0 */
257		347628	hess[4] = VL_4( 0,Z); /* h_yy */
258		347628	}
259			/* x1 */
260		347628	ivY[ 0] = D1_4( 0,X);
261		347628	ivY[ 1] = D1_4( 1,X);
262		347628	ivY[ 2] = D1_4( 2,X);
263		347628	ivY[ 3] = D1_4( 3,X);
264		347628	ivY[ 4] = D1_4( 4,X);
265		347628	ivY[ 5] = D1_4( 5,X);
266		347628	ivY[ 6] = D1_4( 6,X);
267		347628	ivY[ 7] = D1_4( 7,X);
268		347628	ivY[ 8] = D1_4( 8,X);
269		347628	ivY[ 9] = D1_4( 9,X);
270		347628	ivY[10] = D1_4(10,X);
271		347628	ivY[11] = D1_4(11,X);
272		347628	ivY[12] = D1_4(12,X);
273		347628	ivY[13] = D1_4(13,X);
274		347628	ivY[14] = D1_4(14,X);
275		347628	ivY[15] = D1_4(15,X);
276			/* x1y0 */
277		347628	ivZ[ 0] = VL_4( 0,Y);
278		347628	ivZ[ 1] = VL_4( 1,Y);
279		347628	ivZ[ 2] = VL_4( 2,Y);
280		347628	ivZ[ 3] = VL_4( 3,Y);
281			/* x1y0z0 */
282	✓✗	347628	if (doD1) {
283		347628	gvec[0] = VL_4( 0,Z); /* g_x */
284		347628	}
285
286		347628	ell_3mv_mul_d(gvec, shape->ItoWSubInvTransp, gvec);
287
288	✗✓	347628	if (!doD2)
289			return;
290
291			/* x1y0z1 */
292		347628	hess[2] = hess[6] = D1_4( 0,Z); /* h_xz */
293			/* x1y1 */
294		347628	ivZ[ 0] = D1_4( 0,Y);
295		347628	ivZ[ 1] = D1_4( 1,Y);
296		347628	ivZ[ 2] = D1_4( 2,Y);
297		347628	ivZ[ 3] = D1_4( 3,Y);
298			/* x1y1z0 */
299		347628	hess[1] = hess[3] = VL_4( 0,Z); /* h_xy */
300			/* x2 */
301		347628	ivY[ 0] = D2_4( 0,X);
302		347628	ivY[ 1] = D2_4( 1,X);
303		347628	ivY[ 2] = D2_4( 2,X);
304		347628	ivY[ 3] = D2_4( 3,X);
305		347628	ivY[ 4] = D2_4( 4,X);
306		347628	ivY[ 5] = D2_4( 5,X);
307		347628	ivY[ 6] = D2_4( 6,X);
308		347628	ivY[ 7] = D2_4( 7,X);
309		347628	ivY[ 8] = D2_4( 8,X);
310		347628	ivY[ 9] = D2_4( 9,X);
311		347628	ivY[10] = D2_4(10,X);
312		347628	ivY[11] = D2_4(11,X);
313		347628	ivY[12] = D2_4(12,X);
314		347628	ivY[13] = D2_4(13,X);
315		347628	ivY[14] = D2_4(14,X);
316		347628	ivY[15] = D2_4(15,X);
317			/* x2y0 */
318		347628	ivZ[ 0] = VL_4( 0,Y);
319		347628	ivZ[ 1] = VL_4( 1,Y);
320		347628	ivZ[ 2] = VL_4( 2,Y);
321		347628	ivZ[ 3] = VL_4( 3,Y);
322			/* x2y0z0 */
323		347628	hess[0] = VL_4( 0,Z); /* h_xx */
324
325			if (1) {
326			double matA[9];
327		347628	ELL_3M_MUL(matA, shape->ItoWSubInvTransp, hess);
328		347628	ELL_3M_MUL(hess, matA, shape->ItoWSubInv);
329			}
330
331		347628	return;
332		477088	}
333
334			void
335			gageScl3PFilter6(gageShape *shape,
336			double ivX, double ivY, double *ivZ,
337			double fw0, double fw1, double *fw2,
338			double val, double gvec, double *hess,
339			const int *needD) {
340			int i, j;
341			double T;
342			int doV, doD1, doD2;
343		29786	doV = needD[0];
344		14893	doD1 = needD[1];
345		14893	doD2 = needD[2];
346
347			#define fd 6
348			#include "scl3pfilterbody.c"
349			#undef fd
350
351		8893	return;
352		14893	}
353
354			void
355			gageScl3PFilter8(gageShape *shape,
356			double ivX, double ivY, double *ivZ,
357			double fw0, double fw1, double *fw2,
358			double val, double gvec, double *hess,
359			const int *needD) {
360			int i, j;
361			double T;
362			int doV, doD1, doD2;
363		331056	doV = needD[0];
364		165528	doD1 = needD[1];
365		165528	doD2 = needD[2];
366
367			#define fd 8
368			#include "scl3pfilterbody.c"
369			#undef fd
370
371		83628	return;
372		165528	}
373
374			void
375			gageScl3PFilterN(gageShape *shape, int fd,
376			double ivX, double ivY, double *ivZ,
377			double fw0, double fw1, double *fw2,
378			double val, double gvec, double *hess,
379			const int *needD) {
380			int i, j;
381			double T;
382			int doV, doD1, doD2;
383		424656	doV = needD[0];
384		212328	doD1 = needD[1];
385		212328	doD2 = needD[2];
386
387			#include "scl3pfilterbody.c"
388
389		107028	return;
390		212328	}
391
392			void
393			_gageSclFilter(gageContext ctx, gagePerVolume pvl) {
394		1789650	char me[]="_gageSclFilter";
395			int fd;
396			double fw00, fw11, *fw22;
397		894825	gageScl3PFilter_t *filter[5] = {NULL, gageScl3PFilter2, gageScl3PFilter4,
398			gageScl3PFilter6, gageScl3PFilter8};
399
400		894825	fd = 2*ctx->radius;
401	✗✓	894825	if (!ctx->parm.k3pack) {
402			fprintf(stderr, "!%s: sorry, 6-pack filtering not implemented\n", me);
403			return;
404			}
405		894825	fw00 = ctx->fw + fd3gageKernel00;
406		894825	fw11 = ctx->fw + fd3gageKernel11;
407		894825	fw22 = ctx->fw + fd3gageKernel22;
408			/* perform the filtering */
409	✓✓	894825	if (fd <= 8) {
410		1629594	filter[ctx->radius](ctx->shape, pvl->iv3, pvl->iv2, pvl->iv1,
411			fw00, fw11, fw22,
412		814797	pvl->directAnswer[gageSclValue],
413		814797	pvl->directAnswer[gageSclGradVec],
414		814797	pvl->directAnswer[gageSclHessian],
415		814797	pvl->needD);
416		814797	} else {
417		160056	gageScl3PFilterN(ctx->shape, fd,
418		80028	pvl->iv3, pvl->iv2, pvl->iv1,
419			fw00, fw11, fw22,
420		80028	pvl->directAnswer[gageSclValue],
421		80028	pvl->directAnswer[gageSclGradVec],
422		80028	pvl->directAnswer[gageSclHessian],
423		80028	pvl->needD);
424			}
425
426		894825	return;
427		894825	}
428
429			#undef X
430			#undef Y
431			#undef Z