/*

  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 "nrrd.h"

#include "privateNrrd.h"

/*

** learned: even when using doubles, because of limited floating point

** precision, you can get different results between quantizing

** unrescaled (value directly from nrrd, map domain set to nrrd range,

** as with early behavior of unu rmap) and rescaled (value from nrrd

** scaled to fit in existing map domain, as with unu imap -r) value,

** to the exact same index space.

*/

/*

** I won't try to support mapping individual values through a

** colormap, as with a function evaluation on a single passed value.

** That will be handled in an upcoming library...

*/

/*

** this identifies the different kinds of 1D maps, useful for the

** functions in this file only

*/

enum {

  kindLut=0,

  kindRmap=1,

  kindImap=2

};

double

_nrrdApplyDomainMin(const Nrrd *nmap, int ramps, int mapAxis) {

  double ret;

  AIR_UNUSED(ramps);

}

double

_nrrdApplyDomainMax(const Nrrd *nmap, int ramps, int mapAxis) {

  double ret;

  }

}

/*

** _nrrdApply1DSetUp()

**

** some error checking and initializing needed for 1D LUTS, regular,

** and irregular maps.  The intent is that if this succeeds, then

** there is no need for any further error checking.

**

** The only thing this function DOES is allocate the output nrrd, and

** set meta information.  The rest is just error checking.

**

** The given NrrdRange has to be fleshed out by the caller: it can't

** be NULL, and both range->min and range->max must exist.

*/

int

_nrrdApply1DSetUp(Nrrd *nout, const Nrrd *nin, const NrrdRange *range,

                  const Nrrd *nmap, int kind, int typeOut,

                  int rescale, int multi) {

  static const char me[]="_nrrdApply1DSetUp";

  char *mapcnt;

                                    "regular map",

                                    "irregular map"};

                                     "multi regular map",

                                     "multi irregular map"};

                                      /* wishful thinking */

                                    "rmap",

                                    "imap"};

                                     "mrmap",

                                     "mimap"}; /* wishful thinking */

  unsigned int ax, dim, entLen;

  double domMin, domMax;

  }

  }

  }

    }

    }

    /* HEY: consider adding an error check for range->min == range->max

       here; the code below now guards

       AIR_AFFINE(range->min, val, range->max, domMin, domMax)

       against producing NaNs, but maybe that's too clever */

  }

      }

      }

      /* need to make sure the relevant sizes match */

                   me, taxi, NRRD_DIM_MAX-1);

        }

                   "nin->axis[%d].size (%s) "

                   "!= nmap->axis[%d].size (%s): ", me,

                   mapAxis + 1 + ax,

          return 1;

        }

      }

    }

               me, mapAxis, domMin, domMax);

    }

    }

    }

    /* its an irregular map */

    }

    /* mapAxis has no meaning for irregular maps, but we'll pretend ... */

    copyMapAxis0 = AIR_TRUE;

  }

             "NRRD_DIM_MAX %d",

             me, nin->dim,

  }

  }

  }

  /*

  fprintf(stderr, "##%s: pre maybe alloc: nout->data = %p\n", me, nout->data);

  for (dim=0; dim<mapAxis + nin->dim; dim++) {

    fprintf(stderr, "    size[%d] = %d\n", d, (int)size[d]);

  }

  fprintf(stderr, "   nout->dim = %d; nout->type = %d = %s; sizes = %d,%d\n",

          nout->dim, nout->type,

          airEnumStr(nrrdType, nout->type));

  fprintf(stderr, "   typeOut = %d = %s\n", typeOut,

          airEnumStr(nrrdType, typeOut));

  */

  }

  /*

  fprintf(stderr, "   nout->dim = %d; nout->type = %d = %s\n",

          nout->dim, nout->type,

          airEnumStr(nrrdType, nout->type),

          nout->axis[0].size, nout->axis[1].size);

  for (d=0; d<nout->dim; d++) {

    fprintf(stderr, "    size[%d] = %d\n", d, (int)nout->axis[d].size);

  }

  fprintf(stderr, "##%s: post maybe alloc: nout->data = %p\n", me, nout->data);

  */

  }

                      NRRD_AXIS_INFO_SIZE_BIT);

  }

                        nin, "%s", mapcnt)) {

  }

                        NRRD_BASIC_INFO_DATA_BIT

                        | NRRD_BASIC_INFO_TYPE_BIT

                        | NRRD_BASIC_INFO_BLOCKSIZE_BIT

                        | NRRD_BASIC_INFO_DIMENSION_BIT

                        | NRRD_BASIC_INFO_CONTENT_BIT

                           ? 0

                           : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))) {

  }

}

/*

** _nrrdApply1DLutOrRegMap()

**

** the guts of nrrdApply1DLut and nrrdApply1DRegMap

**

** yikes, does NOT use biff, since we're only supposed to be called

** after copious error checking.

**

** FOR INSTANCE, this allows nout == nin, which could be a big

** problem if mapAxis == 1.

**

** we don't need a typeOut arg because nout has already been allocated

** as some specific type; we'll look at that.

**

** NOTE: non-existent values get passed through regular maps and luts

** "unchanged".  However, if the output type is integral, the results

** are probaby undefined.  HEY: there is currently no warning message

** or error handling based on nrrdStateDisallowIntegerNonExist, but

** there really should be.

*/

int

_nrrdApply1DLutOrRegMap(Nrrd *nout, const Nrrd *nin, const NrrdRange *range,

                        const Nrrd *nmap, int ramps, int rescale, int multi) {

  /* static const char me[]="_nrrdApply1DLutOrRegMap"; */

  char *inData, *outData, *mapData, *entData0, *entData1;

  size_t N, I;

  double (*inLoad)(const void *v), (*mapLup)(const void *v, size_t I),

    (*outInsert)(void *v, size_t I, double d),

    val, mapIdxFrac, domMin, domMax;

  unsigned int i, mapAxis, mapLen, mapIdx, entSize, entLen, inSize, outSize;

  }

                                       /* low end of map domain */

                                       /* high end of map domain */

            : 1);

    /* regular map */

      /* ...

      if (!(I % 100)) fprintf(stderr, "I = %d\n", (int)I);

      ... */

      /* ...

      fprintf(stderr, "##%s: val = \na% 31.15f --> ", me, val);

      ... */

               : domMin);

        /* ...

        fprintf(stderr, "\nb% 31.15f (min,max = %g,%g)--> ", val,

                range->min, range->max);

        ... */

      }

      /* ...

      fprintf(stderr, "\nc% 31.15f --> clamp(%g,%g), %d --> ",

              val, domMin, domMax, mapLen);

      ... */

        /* ...

        fprintf(stderr, "mapIdxFrac = \nd% 31.15f --> ",

                mapIdxFrac);

        ... */

        /* ...

        fprintf(stderr, "%s: (%d,\ne% 31.15f) --> \n",

                me, mapIdx, mapIdxFrac);

        ... */

          /* ...

          fprintf(stderr, "f% 31.15f\n", val);

          ... */

        }

      } else {

        /* copy non-existent values from input to output */

        }

      }

      }

    }

  } else {

    /* lookup table */

               : domMin);

      }

        }

      } else {

        /* copy non-existent values from input to output */

        }

      }

      }

    }

  }

}

/*

******** nrrdApply1DLut

**

** A "lut" is a simple lookup table: the data points are evenly spaced,

** with cell-centering assumed, and there is no interpolation except

** nearest neighbor.  The axis min and max are used to determine the

** range of values that can be mapped with the lut.

**

** Of the three kinds of 1-D maps, only luts can have output type block.

**

** If the lut nrrd is 1-D, then the output is a scalar nrrd with the

** same dimension as the input.  If the lut nrrd is 2-D, then each

** value in the input is mapped to a vector of values from the lut,

** which is always a scanline along axis 0.

**

** This allows lut length to be simply 1.

*/

int

nrrdApply1DLut(Nrrd *nout, const Nrrd *nin,

               const NrrdRange *_range, const Nrrd *nlut,

               int typeOut, int rescale) {

  static const char me[]="nrrdApply1DLut";

  NrrdRange *range;

  airArray *mop;

  }

  }

                        rescale, AIR_FALSE /* multi */)

                                 rescale, AIR_FALSE /* multi */)) {

  }

}

int

nrrdApplyMulti1DLut(Nrrd *nout, const Nrrd *nin,

                    const NrrdRange *_range, const Nrrd *nmlut,

                    int typeOut, int rescale) {

  static const char me[]="nrrdApplyMulti1DLut";

  NrrdRange *range;

  airArray *mop;

  }

  }

                        rescale, AIR_TRUE /* multi */)

                                 AIR_FALSE /* ramps */,

                                 rescale, AIR_TRUE /* multi */)) {

  }

}

/*

******** nrrdApply1DRegMap

**

** A regular map has data points evenly spaced, with node-centering and

** and linear interpolation assumed.  As with luts, the axis min and

** max determine the range of values that are mapped.  This function

** is used in nrrdHistoEq() and is the basis of the very popular "unu rmap".

**

** If the lut nrrd is 1-D, then the output is a scalar nrrd with the

** same dimension as the input.  If the lut nrrd is 2-D, then each

** value in the input is mapped to a linear weighting of vectors

** from the map; the vectors are the scanlines along axis 0.

**

** NB: this function makes NO provisions for non-existent input values.

** There won't be any memory errors, but the results are undefined.

**

** This allows map length to be simply 1.

*/

int

nrrdApply1DRegMap(Nrrd *nout, const Nrrd *nin,

                  const NrrdRange *_range, const Nrrd *nmap,

                  int typeOut, int rescale) {

  static const char me[]="nrrdApply1DRegMap";

  NrrdRange *range;

  airArray *mop;

  }

  }

                        rescale, AIR_FALSE /* multi */)

                                 rescale, AIR_FALSE /* multi */)) {

  }

}

int

nrrdApplyMulti1DRegMap(Nrrd *nout, const Nrrd *nin,

                       const NrrdRange *_range, const Nrrd *nmmap,

                       int typeOut, int rescale) {

  static const char me[]="nrrdApplyMulti1DRegMap";

  NrrdRange *range;

  airArray *mop;

  }

  }

                        rescale, AIR_TRUE /* multi */)

                                 rescale, AIR_TRUE /* multi */)) {

  }

}

/*

******** nrrd1DIrregMapCheck()

**

** return zero only for the valid forms of 1D irregular map.

** imap must be 2D, both sizes >= 2, non-block-type, no non-existent

** values in range.  If the first point's position is non-existent,

** than the first three points positions must be -inf, NaN, and +inf,

** and none of the other points locations can be non-existent, and

** they must increase monotonically.  There must be at least two

** points with existent positions.

*/

int

nrrd1DIrregMapCheck(const Nrrd *nmap) {

  static const char me[]="nrrd1DIrregMapCheck";

  double (*mapLup)(const void *v, size_t I);

  int i, entLen, mapLen, baseI;

  Nrrd *nrange;

  }

  }

  }

             me, nmap->dim);

  }

             me, entLen, mapLen);

  }

  }

  }

    baseI = 0;

    baseI = 3;

               "be >= 5 (not %d)", me, mapLen);

    }

               "of 1st three entries (%g:%d,%g:%d,%g:%d) not "

               "-inf, NaN, and +inf", me,

    }

  }

    }

  }

    }

  }

}

/*

******** nrrd1DIrregAclCheck()

**

** returns zero only on valid accelerators for 1D irregular mappings

*/

int

nrrd1DIrregAclCheck(const Nrrd *nacl) {

  static const char me[]="nrrd1DIrregAclCheck";

  }

  }

  }

  }

    return 1;

  }

}

/*

** _nrrd1DIrregMapDomain()

**

** ALLOCATES an array of doubles storing the existent control point

** locations, and sets its length in *poslenP.  If there are the three

** points with non-existent locations, these are ignored.

**

** Assumes that nrrd1DIrregMapCheck has been called on "nmap".

*/

double *

_nrrd1DIrregMapDomain(int *posLenP, int *baseIP, const Nrrd *nmap) {

  static const char me[]="_nrrd1DIrregMapDomain";

  int i, entLen, baseI, posLen;

  double *pos, (*mapLup)(const void *v, size_t I);

  }

  }

  }

  }

}

/*

** _nrrd1DIrregFindInterval()

**

** The hard part of doing 1D irregular mapping: given an array of

** control point locations, and a value, find which interval the value

** lies in.  The lowest and highest possible indices are given in

** "loI" and "hiI".  Results are undefined if these do not in fact

** bound the location of correct interval, or if loI > hiI, or if the

** query positon "p" is not in the domain vector "pos".  Intervals are

** identified by the integral index of the LOWER of the two control

** points spanning the interval.

**

** This imposes the same structure of half-open intervals that

** is done by airIndex().  That is, a value p is in interval i

** if pos[i] <= p < pos[i+1] for all but the last interval, and

** pos[i] <= p <= pos[i+1] for the last interval (in which case

** i == hiI)

*/

int

_nrrd1DIrregFindInterval(const double *pos, double p, int loI, int hiI) {

  int midI;

  /*

  fprintf(stderr, "##%s(%g): hi: %d/%g-%g | %d/%g-%g\n",

          "_nrrd1DIrregFindInterval", p,

          loI, pos[loI], pos[loI+1],

          hiI, pos[hiI], pos[hiI+1]);

  */

      /* p is between (pos[midI],pos[midI+1]): we're done */

      loI = hiI = midI;

      /* p is below interval midI: midI-1 is valid upper bound */

      /* p is above interval midI: midI+1 is valid lower bound */

    }

    /*

    fprintf(stderr, "##%s(%g): %d/%g-%g | %d/%g-%g | %d/%g-%g\n",

            "_nrrd1DIrregFindInterval", p,

            loI, pos[loI], pos[loI+1],

            midI, pos[midI], pos[midI+1],

            hiI, pos[hiI], pos[hiI+1]);

    */

  }

}

/*

******** nrrd1DIrregAclGenerate()

**

** Generates the "acl" that is used to speed up the action of

** nrrdApply1DIrregMap().  Basically, the domain of the map

** is quantized into "acllen" bins, and for each bin, the

** lowest and highest possible map interval is stored. This

** either obviates or speeds up the task of finding which

** interval contains a given value.

**

** Assumes that nrrd1DIrregMapCheck has been called on "nmap".

*/

int

nrrd1DIrregAclGenerate(Nrrd *nacl, const Nrrd *nmap, size_t aclLen) {

  static const char me[]="nrrd1DIrregAclGenerate";

  unsigned int ii;

  unsigned short *acl;

  double lo, hi, min, max, *pos;

  }

    return 1;

  }

                        AIR_CAST(size_t, 2), AIR_CAST(size_t, aclLen))) {

  }

  }

  }

}

/*

******** nrrdApply1DIrregMap()

**

** Linear interpolation between irregularly spaced control points.

** Obviously, the location of the control point has to be given

** explicitly.  The map nrrd must have dimension 2, and each

** control point is represented by a scanline along axis 0.  The

** first value is the position of the control point, and the remaining

** value(s) are linearly weighted according to the position of the

** input value among the control point locations.

**

** To allow "coloring" of non-existent values -inf, NaN, and +inf, if

** the very first value of the map (the location of the first control

** point) is non-existent, then the first three control point locations

** must be -inf, NaN, and +inf, in that order, and the information

** about these points will be used for corresponding input values.

** Doing this makes everything slower, however, because airFPClass_f()

** is called on every single value.

**

** This assumes that nrrd1DIrregMapCheck has been called on "nmap",

** and that nrrd1DIrregAclCheck has been called on "nacl" (if it is

** non-NULL).

*/

int

nrrdApply1DIrregMap(Nrrd *nout, const Nrrd *nin, const NrrdRange *_range,

                    const Nrrd *nmap, const Nrrd *nacl,

                    int typeOut, int rescale) {

  static const char me[]="nrrdApply1DIrregMap";

  size_t N, I;

    entSize, colSize, inSize, lo, hi, baseI;

  double val, *pos, domMin, domMax, mapIdxFrac,

    (*mapLup)(const void *v, size_t I),

    (*inLoad)(const void *v), (*outInsert)(void *v, size_t I, double d);

  char *inData, *outData, *entData0, *entData1;

  NrrdRange *range;

  airArray *mop;

  /*

  fprintf(stderr, "!%s: rescale = %d\n", me, rescale);

  */

  }

  }

                        kindImap, typeOut, rescale, AIR_FALSE /* multi */)) {

  }

  }

    acl = NULL;

    aclLen = 0;

  }

  }

  /*

  fprintf(stderr, "!%s: domMin, domMax = %g, %g\n", me, domMin, domMax);

  */

    /*

    fprintf(stderr, "!%s: (%d) val = % 31.15f\n", me, (int)I, val);

    */

      /* got a non-existent value */

        /* and we know how to deal with them */

        case airFP_NEG_INF:

          mapIdx = 0;

        case airFP_SNAN:

        case airFP_QNAN:

          mapIdx = 1;

        case airFP_POS_INF:

          mapIdx = 2;

        default:

          mapIdx = 0;

                  "not handled\n",

        }

        }

        continue;  /* we're done! (with this value) */

      } else {

        /* we don't know how to properly deal with this non-existent value:

           we use the first entry, and then fall through to code below */

        mapIdx = 0;

        mapIdxFrac = 0.0;

      }

      /* we have an existent value */