/*

  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"

/*

** RETIRED: nrrdMinMaxSet()

**

** Sets nrrd->min and nrrd->max to the extremal (existent) values in

** the given nrrd, by calling the appropriate member of nrrdMinMaxExactFind[]

**

** calling this function will result in nrrd->hasNonExist being set

** (because of the nrrdFindMinMax[] functions)

**

** decided NOT to use biff, so that this is a more distinct alternative to

** nrrdMinMaxCleverSet().

void

nrrdMinMaxSet(Nrrd *nrrd) {

  NRRD_TYPE_BIGGEST _min, _max;

  if (nrrd) {

    if (!airEnumValCheck(nrrdType, nrrd->type)

        && nrrdTypeBlock != nrrd->type) {

      nrrdFindMinMax[nrrd->type](&_min, &_max, nrrd);

      nrrd->min = nrrdDLoad[nrrd->type](&_min);

      nrrd->max = nrrdDLoad[nrrd->type](&_max);

    } else {

      nrrd->min = nrrd->max = AIR_NAN;

    }

  }

  return;

}

*/

/*

** RETIRED: nrrdMinMaxCleverSet()

**

** basically a wrapper around nrrdMinMaxSet(), with bells + whistles:

** 1) will call nrrdMinMaxSet only when one of nrrd->min and nrrd->max

**    are non-existent, with the end result that only the non-existent

**    values are over-written

** 2) obeys the nrrdStateClever8BitMinMax global state to short-cut

**    finding min and max for 8-bit data.  Values for nrrd->min or

**    nrrd->max which were existent to start with are untouched.

** 3) reports error if there are no existent values in nrrd (AIR_EXISTS()

**    fails on every value)

**

** Like nrrdMinMaxSet(), this will always set nrrd->hasNonExist.

**

** Uses biff.

int

nrrdMinMaxCleverSet(Nrrd *nrrd) {

  static const char me[]="nrrdMinMaxCleverSet";

  double min, max;

  if (!nrrd) {

    biffAddf(NRRD, "%s: got NULL pointer", me);

    return 1;

  }

  if (airEnumValCheck(nrrdType, nrrd->type)) {

    biffAddf(NRRD, "%s: input nrrd has invalid type (%d)", me, nrrd->type);

    return 1;

  }

  if (nrrdTypeBlock == nrrd->type) {

    sprintf(err, "%s: can't find min/max of type %s", me,

            airEnumStr(nrrdType, nrrdTypeBlock));

  }

  if (AIR_EXISTS(nrrd->min) && AIR_EXISTS(nrrd->max)) {

    / * both of min and max already set, so we won't look for those, but

       we have to comply with stated behavior of always setting hasNonExist * /

    nrrdHasNonExistSet(nrrd);

    return 0;

  }

  if (nrrdStateClever8BitMinMax

      && (nrrdTypeChar == nrrd->type || nrrdTypeUChar == nrrd->type)) {

    if (nrrdTypeChar == nrrd->type) {

      if (!AIR_EXISTS(nrrd->min))

        nrrd->min = SCHAR_MIN;

      if (!AIR_EXISTS(nrrd->max))

        nrrd->max = SCHAR_MAX;

    } else {

      if (!AIR_EXISTS(nrrd->min))

        nrrd->min = 0;

      if (!AIR_EXISTS(nrrd->max))

        nrrd->max = UCHAR_MAX;

    }

    nrrdHasNonExistSet(nrrd);

    return 0;

  }

  / * at this point we need to find either min and/or max (at least

     one of them was non-existent on the way in) * /

  / * save incoming values in case they exist * /

  min = nrrd->min;

  max = nrrd->max;

  / * this will set nrrd->min, nrrd->max, and hasNonExist * /

  nrrdMinMaxSet(nrrd);

  if (!( AIR_EXISTS(nrrd->min) && AIR_EXISTS(nrrd->max) )) {

    biffAddf(NRRD, "%s: no existent values!", me);

    return 1;

  }

  / * re-enstate the existent incoming min and/or max values * /

  if (AIR_EXISTS(min))

    nrrd->min = min;

  if (AIR_EXISTS(max))

    nrrd->max = max;

  return 0;

}

*/

static int

clampRoundConvert(Nrrd *nout, const Nrrd *nin, int type,

                  int doClamp, int roundDir) {

  static const char me[]="clampRoundConvert";

  }

  }

  }

    /* there's a risk of non-existent values getting converted to

       non-sensical integral values */

               "non-existent values in input", me,

    }

  }

  /* if we're actually converting to the same type, just do a copy */

      /* nout == nin is allowed if the input and output type are

         of the same size, which will certainly be the case if the

         input and output types are identical, so there's actually

         no work to do */

    } else {

      }

    }

  } else {

    /* allocate space if necessary */

    /* MUST be nrrdMaybeAlloc_nva (not nrrd Alloc_nva) because we allow

       nout==nin if type sizes match */

    }

    /* call the appropriate converter */

                                                 doClamp, roundDir);

    }

    /* copy peripheral information */

    }

    /* the min and max have probably changed if there was a conversion

       to integral values, or to a lower precision representation */

                          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

                          | NRRD_BASIC_INFO_COMMENTS_BIT

                             ? 0

                             : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))) {

    }

  }

/*

******** nrrdConvert()

**

** copies values from one type of nrrd to another, without any

** transformation, except what you get with a cast.  The point is to

** make available on Nrrds the exact same behavior as you have in C

** with casts and assignments.

*/

int

nrrdConvert(Nrrd *nout, const Nrrd *nin, int type) {

  static const char me[]="nrrdConvert";

                        AIR_FALSE /* clamp */,

                        0 /* round */)) {

  }

/*

******** nrrdClampConvert()

**

** same as nrrdConvert, but with clamping to output value representation range

** so as to avoid wrap-around artifacts

**

** HEY: WARNING: may have loss of data when processing long long

** (either signed or unsigned)

*/

int

nrrdClampConvert(Nrrd *nout, const Nrrd *nin, int type) {

  static const char me[]="nrrdClampConvert";

                        AIR_TRUE  /* clamp */,

                        0 /* round */)) {

  }

}

/*

******** nrrdCastClampRound()

**

** For (maybe) doing rounding to integer, then (maybe) clamping

** and then casting.  The maybe-ness of rounding and clamping

** in this function is inconsistent with the certainty of clamping

** nrrdClampConvert, so this function name may be regretted.

**

** NOTE! Rounding is not performed when outType is for float

** or double!  That logic is implemented here.

**

** And warning same as above:

** HEY: WARNING: may have loss of data when processing long long

** (either signed or unsigned)

*/

int

nrrdCastClampRound(Nrrd *nout, const Nrrd *nin, int outType,

                   int doClamp, int roundDir) {

  static const char me[]="nrrdCastClampRound";

  }

}

/*

******** nrrdQuantize()

**

** convert values to 8, 16, or 32 bit unsigned quantities

** by mapping the value range delimited by the nrrd's min

** and max to the representable range

**

** NOTE: for the time being, this uses a "double" as the intermediate

** value holder, which may mean needless loss of precision

*/

int

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

             unsigned int bits) {

  static const char me[]="nrrdQuantize", func[]="quantize";

  double valIn, minIn, maxIn, eps;

  int type=nrrdTypeUnknown;

  unsigned char *outUC;

  unsigned short *outUS;

  unsigned int *outUI;

  airArray *mop;

  NrrdRange *range;

  }

  }

  /* determine nrrd type from number of bits */

  default:

  }

  }

  }

             "(NaN, +/-inf)", me);

  }

  /* allocate space if necessary */

  /* MUST be nrrdMaybeAlloc_nva (not nrrd Alloc_nva) because we allow

     nout==nin if type sizes match */

  }

  /* the skinny */

  case 8:

    }

    break;

  case 16:

    }

    break;

  case 32:

                          airIndexClampULL(minIn, valIn, maxIn+eps,

                                           AIR_ULLONG(1) << 32));

    }

    break;

  }

  /* set information in new volume */

  }

                        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

                        | NRRD_BASIC_INFO_OLDMIN_BIT

                        | NRRD_BASIC_INFO_OLDMAX_BIT

                        | NRRD_BASIC_INFO_COMMENTS_BIT

                           ? 0

                           : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))) {

  }

/*

** _nrrdTypeNumberOfValues[]

**

** only meaningful for integral values, and only correct for

** 32-bit values; tells the number of different integral values that

** can be represented by the type

**

** HEY: this used to be public, but that was stopped when it was clear

** that only the following function was using the array.  The reason

** for the array is questionable, and the implementation below

** should be re-evaluated.

*/

static const double

_nrrdTypeNumberOfValues[NRRD_TYPE_MAX+1] = {

  0,                         /* unknown */

  UCHAR_MAX+1,               /* char */

  UCHAR_MAX+1,               /* unsigned char */

  USHRT_MAX+1,               /* short */

  USHRT_MAX+1,               /* unsigned short */

  (double)UINT_MAX+1,        /* int */

  (double)UINT_MAX+1,        /* unsigned int */

  (double)NRRD_ULLONG_MAX+1, /* long long */

  (double)NRRD_ULLONG_MAX+1, /* unsigned long long */

  0,                         /* float */

  0,                         /* double */

  0                          /* punt */

};

/*

******** nrrdUnquantize()

**

** try to recover floating point values from a quantized nrrd,

** using the oldMin and oldMax values, if they exist.  If they

** don't exist, the output range will be 0.0 to 1.0.  However,

** because we're using NRRD_CELL_POS to recover values,

** the output values will never be exactly 0.0 to 1.0 (or oldMin

** to oldMax).  In unsigned char data, for instance, the value

** V will be mapped to:

** NRRD_CELL_POS(0.0, 1.0, 256, V) ==

** AIR_AFFINE(0, V + 0.5, 256, 0.0, 1.0) ==

** ((double)(1.0)-(0.0))*((double)(V+0.5)-(0))/((double)(256)-(0)) + (0.0) ==

** (1.0)*(V+0.5) / (256.0) + (0.0) ==

** (V+0.5)/256

** so a 0 will be mapped to 1/512 = 0.00195

*/

int

nrrdUnquantize(Nrrd *nout, const Nrrd *nin, int type) {

  static const char me[]="nrrdUnquantize", func[]="unquantize";

  float *outF;

  double *outD, minIn, numValIn, minOut, maxOut, valIn;

  }

  }

  }

  }

  }

  }

  }

    minOut = nin->oldMin;

    maxOut = nin->oldMax;

    minOut = 0.0;

    maxOut = 1.0;

  }

  case nrrdTypeFloat:

                          NRRD_CELL_POS(minOut, maxOut, numValIn, valIn));

    }

    break;

  case nrrdTypeDouble:

    }

    break;

  }

  /* set information in new volume */

  }

                        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

                        | NRRD_BASIC_INFO_OLDMIN_BIT

                        | NRRD_BASIC_INFO_OLDMAX_BIT

                        | NRRD_BASIC_INFO_COMMENTS_BIT

                           ? 0

                           : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))) {

  }

/*

** _nrrdHistoEqCompare()

**

** used by nrrdHistoEq in smart mode to sort the "steady" array

** in _descending_ order

*/

int

_nrrdHistoEqCompare(const void *a, const void *b) {

}

/*

******** nrrdHistoEq()

**

** performs histogram equalization on given nrrd, treating it as a

** big one-dimensional array.  The procedure is as follows:

** - create a histogram of nrrd (using "bins" bins)

** - integrate the histogram, and normalize and shift this so it is

**   a monotonically increasing function from min to max, where

**   (min,max) is the range of values in the nrrd

** - map the values in the nrrd through the adjusted histogram integral

**

** If the histogram of the given nrrd is already as flat as can be,

** the histogram integral will increase linearly, and the adjusted

** histogram integral should be close to the identity function, so

** the values shouldn't change much.

**

** If the nhistP arg is non-NULL, then it is set to point to

** the histogram that was used for calculation. Otherwise this

** histogram is deleted on return.

**

** This is all that is done normally, when "smart" is == 0.  In

** "smart" mode (activated by setting "smart" to something greater

** than 0), the histogram is analyzed during its creation to detect if

** there are a few bins which keep getting hit with the same value

** over and over.  It may be desirable to ignore these bins in the

** histogram integral because they may not contain any useful

** information, and so they should not effect how values are

** re-mapped.  The value of "smart" is the number of bins that will be

** ignored.  For instance, use the value 1 if the problem with naive

** histogram equalization is a large amount of background (which is

** exactly one fixed value).

*/

int

nrrdHistoEq(Nrrd *nout, const Nrrd *nin, Nrrd **nmapP,

            unsigned int bins, unsigned int smart, float amount) {

  static const char me[]="nrrdHistoEq", func[]="heq";

  Nrrd *nhist, *nmap;

  size_t I, num;

  airArray *mop;

  NrrdRange *range;

  unsigned bii;

  }

  }

  }

  /* we start by simply copying, because the only thing we're

     changing is the values themselves, and all peripheral

     information is unchanged by this value remapping */

    }

  }

  }

    }

    /* for "smart" mode, we have to some extra work while creating the

       histogram to look for bins incessantly hit with the exact same

       value */

    }

    /* allocate the respect, steady, and last arrays */

    }

    /* steady[0 + 2*bii] == how many times has bin bii seen the same value

       steady[1 + 2*bii] == bii (steady will be rearranged by qsort()) */

    }

    /* now create the histogram */

               "Min and max are equivalent (min,max = %g).", me, range->min);

    }

    min = range->min;

    max = range->max;

                               : 0);

        }

      }

    }

    /* now sort the steady array */

    /* we ignore some of the bins according to "smart" arg */

      /* printf("%s: disrespecting bin %d\n", me, steady[1+2*bii]); */

    }

  }

  }

  /* integrate the histogram then normalize it */

                                                 : 1);

    }

  }

  /* if we've done smart, the integral will have little flat spots

     where we ignored hits in the histogram.  That means the mapping

     will actually be lower at that value than it should be.  In

     truth, we should be using an irregular mapping for this, and the

     control points at the ignored bins should just be missing.  So we

     have to do this silliness to raise those control points in the

     regular map. */

    /* there are bins+1 control points, with indices 0 to bins.

       We'll fix control points 1 to bins-1.  ycoord[bii] is too low

       if hist[bii-1] was not respected (!respect[bii-1]) */

        /* lort and hirt will bracket the index of the bad control point

           with points corresponding either to respected bins or the

           endpoints of the histogram */

          ;

          ;

                                 ycoord[lort], ycoord[hirt]);

      }

    }

    /* the very last control point has to be handled differently */

    }

  }

  /* rescale the histogram integration to span the original

     value range, and affect the influence of "amount" */

                             AIR_AFFINE(0, bii, bins, min, max),

                             ycoord[bii]);

  }

  /* map the nrrd values through the normalized histogram integral */

  }

  /*

  for (I=0; I<num; I++) {

    val = nrrdDLookup[nin->type](nin->data, I);

    if (AIR_EXISTS(val)) {

      AIR_INDEX(min, val, max, bins, idx);

      val = AIR_AFFINE(xcoord[idx], val, xcoord[idx+1],

                       ycoord[idx], ycoord[idx+1]);

    }

    nrrdDInsert[nout->type](nout->data, I, val);

  }

  */

  /* if user is interested, set pointer to map nrrd,

     otherwise it will be nixed by airMop */

  }

  /* fiddling with content is the only thing we'll do */

  }

                        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

                        | NRRD_BASIC_INFO_COMMENTS_BIT

                           ? 0

                           : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))) {

  }

}