/*

  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 "teem/meet.h"

/*

** Tests:

** ... lots of gage stuff ...

**

** The main point of all this is to make sure of two (or maybe 2.5) separate

** things about gage, the testing of which requires so much in common that one

** program might as well do them all. The tasks are:

**

** 1) ensure that values and their derivatives (where the gageKind supports

** it) are correctly handled in the multi-value gageKinds (gageKindVec,

** tenGageKind, tenDwiGageKind), relative to the gageKindScl ground-truth: the

** per-component values and derivatives have to match those reconstructed from

** sets of scalar volumes of the components.

** ==> Also, that gageContextCopy works even on dynamic kinds (like DWIs)

**

** 1.5) that there is expected consistency between the scalar, vector, tensor,

** and DWI properties of a related set of volumes.  So the test starts with a

** diffusion tensor field and generates DWIs, scalars (norm squared of the

** tensor), and vectors (gradient of norm squared) from this.

**

** 2) that scale-space reconstruction works: that sets of pre-blurred volumes

** can be generated and saved via the utilities in meet, that the smart

** hermite-spline based scale-space interpolation is working (for all kinds),

** and that gageContextCopy works on scale-space contexts

*/

#define BKEY "probeSS"

#define KIND_NUM 4

#define KI_SCL 0

#define KI_VEC 1

#define KI_TEN 2

#define KI_DWI 3

/*

static const char *kindStr[KIND_NUM] = {

  "scalar",

  "vector",

  "tensor",

  "  DWI "

};

*/

#define HALTON_BASE 100

static unsigned int volSize[3] = {45, 46, 47};

#define NRRD_NEW(nn, mm)                                        \

  (nn) = nrrdNew();                                             \

       airMopAdd((mm), (nn), (airMopper)nrrdNuke, airMopAlways)

/* the weird function names of the various local functions here (created in a

   rather adhoc and organic way) should be read in usual Teem order: from

   right to left */

static int

engageGenTensor(gageContext *gctx, Nrrd *ncten,

                /* out/in */

                double noiseStdv, unsigned int seed,

                unsigned int sx, unsigned int sy, unsigned int sz) {

  static const char me[]="engageGenTensor";

  hestParm *hparm;

  airArray *smop;

  Nrrd *nclean;

  NrrdIter *narg0, *narg1;

  /*   0     1     2    3     4     5     6     7     8     9 */

    {"-s", NULL, NULL, NULL, "-v", "0", "-r", "40", "-o", NULL,

     "-ev", "0.00086", "0.00043", "0.00021", "-bg", "0.003", "-b", "3",

     NULL};

  int helixArgc;

  gagePerVolume *pvl;

  /* NOTE: this is currently the only place where a unrrduCmd

     is called from within C code; it was educational to get working.

     Learned:

     * hest does NOT tolerate having empty or NULL elements of

     its argv[]!  More error checking for this in hest is needed.

     * the "const char **argv" type is not very convenient to

     set up in a dynamic way; the per-element setting done below

     is certainly awkward

  */

  helixArgc = AIR_CAST(int, sizeof(helixArgv)/sizeof(char *)) - 1;

    /* error already went to stderr, not to any biff container */

  }

  }

  /* add some noise to tensor value; no, this isn't really physical;

     since we're adding noise to the tensor and then simulating DWIs,

     rather than adding noise to DWIs and then estimating tensor,

     but for the purposes of gage testing its fine */

                            narg0, narg1)) {

  }

  /* wrap in gage context */

  }

/*

** takes in ncten, measures "S" to nscl, copies that to nsclCopy,

** wraps those in gctx and gctxComp

*/

static int

engageGenScalar(gageContext *gctx, Nrrd *nscl,

                gageContext *gctxComp, Nrrd *nsclCopy,

                /* out/in */

                const Nrrd *ncten) {

  static const char me[]="engageGenScalar";

  gagePerVolume *pvl;

  }

  }

  /* wrap both in gage context */

  }

/*

** Makes a vector volume by measuring the gradient

** Being the gradient of the given scalar volume is just to make

** something vaguely interesting, and to test consistency with

** the the same gradient as measured in the tensor field.

*/

static int

engageGenVector(gageContext *gctx, Nrrd *nvec,

                /* out/in */

                const Nrrd *nscl) {

  static const char me[]="engageGenVector";

  Nrrd *ntmp;

  airArray *smop;

  float *vec, *scl;

  size_t sx, sy, sz, xi, yi, zi, Px, Mx, Py, My, Pz, Mz;

  double dnmX, dnmY, dnmZ;

  gagePerVolume *pvl;

  }

             AIR_CAST(ptrdiff_t, sx-1),

             AIR_CAST(ptrdiff_t, sy-1),

             AIR_CAST(ptrdiff_t, sz-1));

  /* we do axinsert and pad in order to keep all the per-axis info */

                     nrrdBoundaryPad, 0.0)) {

  }

#define INDEX(xj, yj, zj) (xj + sx*(yj + sy*zj))

      }

    }

  }

#undef INDEX

  /* wrap in gage context */

  }

/*

** make a DWI volume by simulating DWIs from given tensor

** this includes generating a new gradient set

*/

static int

genDwi(Nrrd *ndwi, Nrrd *ngrad,

       /* out/in */

       unsigned int gradNum, double bval, const Nrrd *ncten) {

  static const char me[]="genDwi";

  tenGradientParm *gparm;

  tenExperSpec *espec;

  Nrrd *nten, *nb0;

  NrrdIter *narg0, *narg1;

  airArray *smop;

                                       bval,

                                                ngrad->axis[1].size))) {

  }

             ncten->axis[0].size-1,

             ncten->axis[1].size-1,

             ncten->axis[2].size-1,

             ncten->axis[3].size-1);

  }

                            narg0, narg1)) {

  }

                       nb0, nten, AIR_TRUE /* keyValueSet */)) {

  }

int

engageMopDiceVector(gageContext *gctx, Nrrd *nvecComp[3],

                    /* out/in */

                    airArray *mop, const Nrrd* nvec) {

  static const char me[]="engageMopDiceVector";

  gagePerVolume *pvl;

  unsigned int ci;

  }

    }

    }

  }

int

engageMopDiceTensor(gageContext *gctx, Nrrd *nctenComp[7],

                    /* out/in */

                    airArray *mop, const Nrrd* ncten) {

  static const char me[]="engageMopDiceTensor";

  gagePerVolume *pvl;

  unsigned int ci;

                     AIR_TRUE /* useBiff */)) {

  }

    }

    }

  }

int

engageMopDiceDwi(gageContext *gctx, Nrrd ***ndwiCompP,

                 /* out/in */

                 airArray *mop, const Nrrd* ndwi) {

  static const char me[]="mopDiceDwi";

  Nrrd **ndwiComp;

  size_t dwiNum;

  gagePerVolume *pvl;

  unsigned int ci;

  }

  }

  }

    }

    }

  }

typedef struct {

  char name[AIR_STRLEN_SMALL];

  const double **aptr; /* array of pointers to (const) answers */

  gageItemSpec *ispec; /* array of gageItemSpecs (not pointers to them) */

  unsigned int *alen;  /* array of answer lengths */

  unsigned int *sidx;  /* array of index (within san) of copied answer */

  unsigned int anum;   /* lenth of aptr, ispec, alen, sidx arrays */

  double *san;         /* single buffer for copy + concating answers */

  unsigned int slen;   /* length of single buffer (sum of all alen[i]) */

} multiAnswer;

static void

multiAnswerInit(multiAnswer *man) {

  /* HEY sloppy: not resetting name */

static multiAnswer*

multiAnswerNew(char *name) {

  multiAnswer *man;

}

void

multiAnswerLenSet(multiAnswer *man, unsigned int anum) {

  /*

  static const char me[]="multiAnswerLenSet";

  fprintf(stderr, "!%s: %s hello -> answer number = %u\n", me,

          man->name, anum);

  */

  /* don't know answer lengths yet */

}

static multiAnswer*

multiAnswerNix(multiAnswer *man) {

}

static void

multiAnswerAdd(multiAnswer *man, unsigned int ansIdx,

               const gageContext *gctx, const gagePerVolume *pvl,

               unsigned int item) {

  /*

  static const char me[]="multiAnswerAdd";

  fprintf(stderr, "!%s: %s hello (aidx = %u)\n", me, man->name, ansIdx);

  */

  /* HEY hack: doing this tally and allocation only on the last time

     we're called, presuming calls with increasing ansIdx */

    unsigned int ai;

    /* fprintf(stderr, "!%s: %s hello aidx reached anum-1 = %u\n", me, man->name, man->anum-1); */

    }

    /*

    fprintf(stderr, "!%s: (%s) slen = %u\n", "multiAnswerAdd",

            pvl->kind->name, man->slen);

    */

}

static void

multiAnswerCollect(multiAnswer *man) {

  unsigned int ai;

  /*

  static const char me[]="multiAnswerCollect";

  fprintf(stderr, "%s: %s hi\n", me, man->name);

  */

    /*

    fprintf(stderr, "!%s: (%s/%s) ai=%u/%u  to %p+%u=%p for %u doubles\n", "multiAnswerCollect",

            man->ispec[ai].kind->name,

            airEnumStr(man->ispec[ai].kind->enm, man->ispec[ai].item),

            ai, man->anum, man->san, man->sidx[ai], man->san + man->sidx[ai],

            man->alen[ai]);

    */

           man->aptr[ai],

           man->alen[ai]*sizeof(double));

  }

  return;

static int

multiAnswerCompare(multiAnswer *man1, multiAnswer *man2) {

  static const char me[]="multiAnswerCompare";

  unsigned int si, slen;

#if 1

             man1->slen, man2->slen);

  }

  slen = man1->slen;

#else

  slen = AIR_MIN(man1->slen, man2->slen);

#endif

      /* HEY should track down which part of which answer,

         in man1 and man2, is different, which was the

         purpose of recording ispec and sidx */

               man1->san[si], man2->san[si]);

    }

  }

/*

** setting up gageContexts for the first of the two tasks listed above: making

** sure per-component information is handled correctly.  NOTE: The combination

** of function calls made here is very atypical for a Teem-using program

*/

static int

updateQueryKernelSetTask1(gageContext *gctxComp[KIND_NUM],

                          gageContext *gctx[KIND_NUM], int gSetup,

                          multiAnswer *manComp[KIND_NUM],

                          multiAnswer *man[KIND_NUM],

                          NrrdKernel *kpack[3],

                          double support) {

  static const char me[]="updateQueryKernelSetTask1";

  unsigned int kindIdx;

  if (4 != KIND_NUM) {

    biffAddf(BKEY, "%s: sorry, confused: KIND_NUM %u != 4",

             me, KIND_NUM);

    return 1;

  }

                  me, kindIdx);

      }

    }

    /* Note that the contexts for the diced-up volumes of components are always

       of kind gageKindScl, and all the items are from the gageScl* enum */

      unsigned int pvi;

      /*

        fprintf(stderr, "!%s: gctx[%u] has %u pvls\n", me, kindIdx, gctx[kindIdx]->pvlNum);

        fprintf(stderr, "!%s: gctxComp[%u] has %u pvls\n", me, kindIdx, gctxComp[kindIdx]->pvlNum);

      */

                    me, kindIdx, pvi);

        }

      }

                  me, kindIdx);

      }

    /* For the original contexts, we have to use the kind-specific items that

       correspond to the values and derivatives */

    }

                  me, kindIdx);

      }

    }

  } /* if (gSetup) */

  }

    unsigned int pvi;

                     gageSclValue);

    }

                       gageSclGradVec);

      }

                       gageSclHessian);

      }

    }

    case KI_SCL:

    case KI_VEC:

    case KI_TEN:

    case KI_DWI:

    }

  }

static int

probeTask1(gageContext *gctxComp[KIND_NUM],

           gageContext *gctx[KIND_NUM],

           multiAnswer *manComp[KIND_NUM],

           multiAnswer *man[KIND_NUM],

           unsigned int probeNum, NrrdKernel *kpack[3],

           double ksupport) {

  static const char me[]="probeTask1";

  unsigned int kindIdx, probeIdx, errNumMax, tenErrNum, sclErrNum;

    tenDiff[7], tenAvg[7], tenErr, tenErrMax,

    vecDiff[3], vecAvg[3], vecErr, vecErrMax, vecErrNum,

    sclDiff, sclAvg, sclErr, sclErrMax, errNumFrac;

  const double *dwiTenEstP,

    *tenTenP, *tenTenNormP, *tenTenNormGradP,

    *sclSclP, *sclGradP, *vecVecP;

  ELL_3V_SET(minp, ksupport, ksupport, ksupport);

  /* this is all for task 1.5 */

  tenErrNum = sclErrNum = 0;

  vecErrNum = 0.0;

  errNumFrac = 0.02;

    /* not actually here for any derivatives, mainly to check on

       tensor esimation (in addition to usual check of multivariate as

       set of scalars) */

    tenErrMax = 0.0002;

    sclErrMax = AIR_NAN;

    vecErrMax = AIR_NAN;

    sclErrMax = AIR_NAN;

    tenErrMax = 0.14;  /* honestly, not sure how meaningful this test

                          is, given how significant we're allowing the

                          error to be. might be more meaningful if

                          this was with comparison to a pre-computed

                          non-linear least-squares fit, instead of the

                          (log-based) linear least-squares.  */

    vecErrMax = 0.0016;

    sclErrMax = 0.0008;

  } else {

  }

#define PROBE(gctx, what, pos)                                          \

      if (gageProbeSpace((gctx), (pos)[0], (pos)[1], (pos)[2],          \

                         AIR_TRUE /* indexSpace */,                     \

                         AIR_FALSE /* clamp */)) {                      \

        biffAddf(BKEY, "%s: probe (%s) error(%d): %s", me,              \

                 what, (gctx)->errNum, (gctx)->errStr);                 \

        return 1;                                                       \

      }

#undef PROBE

                 probeIdx, kindIdx);

      }

    }

    /* see if probed tensors mostly agree with

       tensors estimated from probed DWIs */

                   me, probeIdx, tenErr, tenErrMax, tenErrNum, errNumMax);

        }

      }

    }

    /* see if invariant gradient learned from tensor volume agrees with

       volume of pre-computed invariant gradients, and with

       gradient of pre-computed invariants */

                   me, probeIdx, vecErr, vecErrMax, vecErrNum, errNumMax);

        }

      }

                   me, probeIdx, vecErr, vecErrMax, vecErrNum, errNumMax);

        }

      }

    }

    /* see if invariant learned from tensor volume agrees with

       volume of precomputed invariants */

                   me, probeIdx, sclErr, sclErrMax, sclErrNum, errNumMax);

        }

      }

    }

  }

static const char *testInfo = "for testing gage";

int

main(int argc, const char **argv) {

  const char *me;

  char *err = NULL;

  hestParm *hparm;

  airArray *mop;

  gageKind *dwikind = NULL;

    *gctx[KIND_NUM], *gctxCopy[KIND_NUM];

    /* these are the volumes that are used in gctxComp[] */

    *nsclCopy, *nvecComp[3], *nctenComp[7], **ndwiComp,

    *ngrad;  /* need access to list of gradients used to make DWIs;

                (this is not the gradient of a scalar field) */

    gradNum = 10; /* small number so that missing one will produce

                     a big reconstruction error */

  /* learn things from hest */

             "kernel support");

             "# of probes");

             "kernel to use; can be: box, cos, or ctmr");

                 AIR_TRUE, AIR_TRUE, AIR_TRUE);

               nrrdKernelBoxSupportDebug,

               nrrdKernelZero,

               nrrdKernelZero);

               nrrdKernelCos4SupportDebug,

               nrrdKernelCos4SupportDebugD,

               nrrdKernelCos4SupportDebugDD);

               nrrdKernelCatmullRomSupportDebug,

               nrrdKernelCatmullRomSupportDebugD,

               nrrdKernelCatmullRomSupportDebugDD);

  } else {

  }

  /* This was a tricky bug: Adding gageContextNix(gctx) to the mop

     as soon as a gctx is created makes sense.  But, in the first

     version of this code, gageContextNix was added to the mop

     BEFORE tenDwiGageKindNix was added, which meant that gageContextNix

     was being called AFTER tenDwiGageKindNix.  However, that meant

     that existing pvls had their pvl->kind being freed from under them,

     but gagePerVolumeNix certainly needs to look at and possibly call

     pvl->kind->pvlDataNix in order clean up pvl->data.  *Therefore*,

     we have to add tenDwiGageKindNix to the mop prior to gageContextNix,