/*

  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"

/*

** sets the filter sample location (fsl) array based on fractional

** probe location ctx->point->frac

**

** One possible rare surpise: if a filter is not continuous with 0

** at the end of its support, and if the sample location is at the

** highest possible point (xi == N-2, xf = 1.0), then the filter

** weights may not be the desired ones.  Forward differencing (via

** nrrdKernelForwDiff) is a good example of this.

*/

void

_gageFslSet(gageContext *ctx) {

  int fr, i;

  double *fslx, *fsly, *fslz;

  double xf, yf, zf;

  case 1:

  case 2:

  default:

    /* filter radius bigger than 2 */

    }

    break;

  }

  return;

/*

** renormalize weights of a reconstruction kernel with

** constraint: the sum of the weights must equal the continuous

** integral of the kernel

*/

void

_gageFwValueRenormalize(gageContext *ctx, int wch) {

  double integral, sumX, sumY, sumZ, *fwX, *fwY, *fwZ;

  int i, fd;

  sumX = sumY = sumZ = 0;

  }

  }

  return;

/*

** renormalize weights of a derivative kernel with

** constraint: the sum of the weights must be zero, but

** sign of individual weights must be preserved

*/

void

_gageFwDerivRenormalize(gageContext *ctx, int wch) {

  double negX, negY, negZ, posX, posY, posZ, fixX, fixY, fixZ,

    *fwX, *fwY, *fwZ;

  int i, fd;

  negX = negY = negZ = 0;

  posX = posY = posZ = 0;

  }

  /* fix is the sqrt() of factor by which the positive values

     are too big.  negative values are scaled up by fix;

     positive values are scaled down by fix */

  }

  }

  return;

void

_gageFwSet(gageContext *ctx, unsigned int sidx, double sfrac) {

  int kidx;

  unsigned int fd;

      continue;

    }

    /* we evaluate weights for all three axes with one call */

  }

        continue;

      }

      case gageKernel00:

      case gageKernel10:

      case gageKernel20:

      default:

      }

    }

    }

  }

    unsigned int jj;

    double scl, norm, *fwX, *fwY, *fwZ;

                     ctx->stackPos[sidx],

                     ctx->stackPos[sidx+1]);

#if 0

    double (*dgeval)(double x, const double *parm),

      dgparm[2] = {0, 3};

    dgeval = nrrdKernelDiscreteGaussian->eval1_d;

    dgparm[0] = scl;

    /* from Eq. (120) in T. Lindeberg. "Feature Detection with Automatic

       Scale Selection." International Journal of Computer Vision,

       1998, 30, 77-116 */

    /* 0.7978845608 ~= sqrt(2)/sqrt(pi) */

    norm = 0.7978845608/(dgeval(0.0, dgparm) + dgeval(1.0, dgparm));

#endif

    /* really simple; no lindeberg normalization, possible bias */

    kidx = gageKernel11;

    }

    kidx = gageKernel22;

    }

  }

  return;

/*

** _gageLocationSet

**

** updates probe location in general context, and things which

** depend on it:

** fsl, fw

**

** (_xi,_yi,_zi) is *index* space position in the volume

** _si is the index-space position in the stack, the value is ignored

** if there is no stack behavior

**

** does NOT use biff, but returns 1 on error and 0 if all okay

** Currently only error is probing outside volume, which sets

** ctx->errNum and sprints message into ctx->errStr.

*/

int

_gageLocationSet(gageContext *ctx,

                 double xif, double yif, double zif, double sif) {

  unsigned int top[3],  /* "top" x, y, z: highest valid index in volume */

    idx[4];

  int sdiff;      /* computed integral positions in volume */

  double frac[4], min, max[3];

  /* **** bounds checking **** */

    min = 0;

    min = -0.5;

  }

              "bounds [%g,%g]x[%g,%g]x[%g,%g]",

              me, xif, yif, zif,

              airEnumStr(nrrdCenter, ctx->shape->center),

              min, max[0], min, max[1], min, max[2]);

    }

  }

                "bounds [0,%u]", me, sif,

                airEnumStr(nrrdCenter, nrrdCenterNode), ctx->pvlNum-2);

      }

    }

  }

  /* **** computing integral and fractional sample locations **** */

  /* Thu Jan 14 19:46:53 CST 2010: detected that along the low edge

     (next to sample 0) in cell centered, the rounding behavior of

     AIR_CAST(unsigned int, xif), namely [-0.5,0] --> 0, meant that

     the low edge was not treated symmetrically with the high edge.

     This motivated the change from using idx to store the lower

     corner of the containing voxel, to the upper corner.  So, the new

     "idx" is always +1 of what the old idx was.  Code here and in

     ctx.c (since idx is saved into ctx->point.idx) has been changed

     accordingly */

             AIR_CAST(unsigned int, xif+1), /* +1: see above */

             AIR_CAST(unsigned int, yif+1),

             AIR_CAST(unsigned int, zif+1));

            "        --> idx %u %u %u\n",

            min, max[0], max[1], max[2], idx[0], idx[1], idx[2]);

  }

  /* these can only can kick in for node-centered, because that's when

     max[] has an integral value */

  }

             xif - (AIR_CAST(float, idx[0])-1),

             yif - (AIR_CAST(float, idx[1])-1),

             zif - (AIR_CAST(float, idx[2])-1));

    idx[3] = 0;

    frac[3] = 0;

    sdiff = AIR_FALSE;

  }

            "        pos (% 15.7f,% 15.7f,% 15.7f,% 15.7f) \n"

            "        -> i(%5d,%5d,%5d,%5d) \n"

            "         + f(% 15.7f,% 15.7f,% 15.7f,% 15.7f) \n",

            frac[0], frac[1], frac[2], frac[3]);

  }

  /* **** compute *spatial* fsl and fw ****

     these have to be reconsidered if anything changes about the

     fractional spatial position, or (if no fractional spatial change),

     movement along scale AND using normalization based on scale */

    /* We don't yet record the scale position in ctx->point because

       that's done below while setting stackFsl and stackFw. So, have

       to pass stack pos info to _gageFwSet() */

    /* these may take some time (especially if using renormalization),

       hence the conditional above */

  /* **** compute *stack* fsl and fw ****  */

  }

    double sum;

    unsigned int ii, nnz;

    NrrdKernelSpec *sksp;

    /* node-centered sampling of stack indices from 0 to ctx->pvlNum-2 */

    /* HEY: we really are quite far from implementing arbitrary

       nrrdBoundary behaviors here, and centering is stuck on node! */

    /* HEY: honestly, the whole idea that it still makes sense to do

       low-level operations in index space, when the world-space locations

       of the samples can be non-uniform, is a little suspect.  This is

       all legit for nrrdKernelTent and nrrdKernelHermiteScaleSpaceFlag,

       but is pretty fishy otherwise */

      }

    }

      }

    }

    /* compute stackFwNonZeroNum whether or not parm.stackNormalizeRecon! */

    nnz = 0;

      sum = 0;

      }

                  "can't do stack reconstruction", me);

        }

      }

      }

        }

      }

    } else {

      }

                  "can't do stack reconstruction", me);

        }

      }

    }

  }