/*

  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"

/*

******** nrrdInvertPerm()

**

** given an array (p) which represents a permutation of n elements,

** compute the inverse permutation ip.  The value of this function

** is not its core functionality, but all the error checking it

** provides.

*/

int

nrrdInvertPerm(unsigned int *invp, const unsigned int *pp, unsigned int nn) {

  static const char me[]="nrrdInvertPerm";

  int problem;

  unsigned int ii;

  }

  /* use the given array "invp" as a temp buffer for validity checking */

               "%s: permutation element #%d == %d out of bounds [0,%d]",

               me, ii, pp[ii], nn-1);

    }

  }

  /* for some reason when this code was written (revision 2700 Sun Jul

     3 04:18:33 2005 UTC) it was decided that all problems with the

     permutation would be reported with a pile of error messages in

     biff; rather than bailing at the first problem.  Not clear if

     this is a good idea. */

  problem = AIR_FALSE;

               me, ii, invp[ii]);

      problem = AIR_TRUE;

    }

  }

  }

  /* the skinny */

  }

}

/*

******** nrrdAxesInsert

**

** like reshape, but preserves axis information on old axes, and

** this is only for adding a "stub" axis with length 1.  All other

** axis attributes are initialized as usual.

*/

int

nrrdAxesInsert(Nrrd *nout, const Nrrd *nin, unsigned int axis) {

  static const char me[]="nrrdAxesInsert", func[]="axinsert";

  unsigned int ai;

  }

             me, axis, nin->dim);

  }

             me, NRRD_DIM_MAX);

  }

                                 ? 0

                                 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT)))) {

    }

  }

                      NRRD_AXIS_INFO_NONE);

  }

  /* the ONLY thing we can say about the new axis is its size */

    /* except maybe the kind */

  }

  /* all basic info has already been copied by nrrdCopy() above */

/*

******** nrrdAxesPermute

**

** changes the scanline ordering of the data in a nrrd

**

** The basic means by which data is moved around is with memcpy().

** The goal is to call memcpy() as few times as possible, on memory

** segments as large as possible.  Currently, this is done by

** detecting how many of the low-index axes are left untouched by

** the permutation- this constitutes a "scanline" which can be

** copied around as a unit.  For permuting the y and z axes of a

** matrix-x-y-z order matrix volume, this optimization produced a

** factor of 5 speed up (exhaustive multi-platform tests, of course).

**

** The axes[] array determines the permutation of the axes.

** axis[i] = j means: axis i in the output will be the input's axis j

** (axis[i] answers: "what do I put here", from the standpoint of the output,

** not "where do I put this", from the standpoint of the input)

*/

int

nrrdAxesPermute(Nrrd *nout, const Nrrd *nin, const unsigned int *axes) {

  static const char me[]="nrrdAxesPermute", func[]="permute";

    lineSize,                /* size of block of memory which can be

                                moved contiguously from input to output,

                                thought of as a "scanline" */

    numLines,                /* how many "scanlines" there are to permute */

    szIn[NRRD_DIM_MAX], *lszIn,

    szOut[NRRD_DIM_MAX], *lszOut,

    cIn[NRRD_DIM_MAX],

    cOut[NRRD_DIM_MAX];

  char *dataIn, *dataOut;

    ai,                      /* running index along dimensions */

    lowPax,                  /* lowest axis which is "p"ermutated */

    ldim,                    /* nin->dim - lowPax */

    ip[NRRD_DIM_MAX+1],      /* inverse of permutation in "axes" */

    laxes[NRRD_DIM_MAX+1];   /* copy of axes[], but shifted down by lowPax

                                elements, to remove i such that i == axes[i] */

  airArray *mop;

  }

  /* we don't actually need ip[], computing it is for error checking */

  }

  /* this shouldn't actually be necessary .. */

  }

    ;

  lowPax = ai;

  /* allocate output by initial copy */

    }

    }

  }

    /* if lowPax == nin->dim, then we were given the identity permutation, so

       there's nothing to do other than the copy already done.  Otherwise,

       here we are (actually, lowPax < nin->dim-1) */

    }

    }

    /* the skinny */

    lineSize = 1;

    }

    }

      /* in our representation of the coordinates of the start of the

         scanlines that we're copying, we are not even storing all the

         zeros in the coordinates prior to lowPax, and when we go to

         a linear index for the memcpy(), we multiply by lineSize */

      }

    }

    /* set content */

    }

    }

                            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))) {

      }

    }

  }

}

/*

******** nrrdShuffle

**

** rearranges hyperslices of a nrrd along a given axis according to

** given permutation.  This could be used to on a 4D array,

** representing a 3D volume of vectors, to re-order the vector

** components.

**

** the given permutation array must allocated for at least as long as

** the input nrrd along the chosen axis.  perm[j] = i means that the

** value at position j in the _new_ array should come from position i

** in the _old_array.  The standpoint is from the new, looking at

** where to find the values amid the old array (perm answers "what do

** I put here", not "where do I put this").  This allows multiple

** positions in the new array to copy from the same old position, and

** insures that there is an source for all positions along the new

** array.

*/

int

nrrdShuffle(Nrrd *nout, const Nrrd *nin, unsigned int axis,

            const size_t *perm) {

  static const char me[]="nrrdShuffle", func[]="shuffle";

  /* Sun Feb 8 13:13:58 CST 2009: There was a memory bug here caused

     by using the same buff1[NRRD_DIM_MAX*30] declaration that had

     worked fine for nrrdAxesPermute and nrrdReshape, but does NOT

     work here because now samples along an axes are re-ordered, not

     axes, so its often not allocated for long enough to hold the

     string that's printed to it. Ideally there'd be another argument

     that says whether to document the shuffle in the content string,

     which would mean an API change.  Or, we can use a secret

     heuristic (or maybe later a nrrdState variable) for determining

     when an axis is short enough to make documenting the shuffle

     interesting.  This is useful since functions like nrrdFlip()

     probably do *not* need the shuffle (the sample reversal) to be

     documented for long axes */

#define LONGEST_INTERESTING_AXIS 42

  unsigned int ai, ldim, len;

    cIn[NRRD_DIM_MAX+1], cOut[NRRD_DIM_MAX+1];

  char *dataIn, *dataOut;

  }

  }

  }

      return 1;

    }

  }

  /* this shouldn't actually be necessary .. */

  }

  /* set information in new volume */

  }

  }

  /* the min and max along the shuffled axis are now meaningless */

  /* do the safe thing first */

  /* try cleverness */

      /* these kinds have no intrinsic ordering */

    }

  }

  /* the skinny */

  lineSize = 1;

  }

  }

  /* Set content. The LONGEST_INTERESTING_AXIS hack avoids the

     previous array out-of-bounds bug */

    }

    }

  } else {

    }

  }

                        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))) {

  }

#undef LONGEST_INTERESTING_AXIS

}

/* ---- BEGIN non-NrrdIO */

/*

******** nrrdAxesSwap()

**

** for when you just want to switch the order of two axes, without

** going through the trouble of creating the permutation array

** needed to call nrrdAxesPermute()

*/

int

nrrdAxesSwap(Nrrd *nout, const Nrrd *nin, unsigned int ax1, unsigned int ax2) {

  static const char me[]="nrrdAxesSwap", func[]="swap";

  }

  }

  }

  }

  /* basic info already copied by nrrdAxesPermute */

}

/*

******** nrrdFlip()

**

** reverse the order of slices along the given axis.

** Actually, just a wrapper around nrrdShuffle() (with some

** extra setting of axis info)

*/

int

nrrdFlip(Nrrd *nout, const Nrrd *nin, unsigned int axis) {

  static const char me[]="nrrdFlip", func[]="flip";

  size_t *perm, si;

  airArray *mop;

  unsigned int axisIdx;

  }

  }

  }

  }

  /* nrrdBasicInfoCopy called by nrrdShuffle() */

  }

                    NRRD_AXIS_INFO_SIZE_BIT

                    | NRRD_AXIS_INFO_KIND_BIT);

  /* HEY: (Tue Jan 18 00:28:26 EST 2005) there's a basic question to

     be answered here: do we want to keep the "location" of the

     samples fixed, while changing their ordering, or do want to flip

     the location of the samples?  In the former, the position

     information has to be flipped to cancel the flipping of the the

     sample order, so that samples maintain location.  In the latter,

     the position information is copied verbatim from the original.  */

  /* (Tue Sep 13 09:59:12 EDT 2005) answer: we keep the "location" of

     the samples fixed, while changing their ordering.  This is the

     low-level thing to do, so for a nrrd function, its the right thing

     to do.  You don't need a nrrd function to simply manipulate

     per-axis meta-information */

  /* HEY: Fri Jan 14 02:53:30 EST 2005: isn't spacing supposed to be

     the step from one sample to the next?  So its a signed quantity.

     If min and max can be flipped (so min > max), then spacing can

     be negative, right?  */

  /* HEY: Fri Jan 14 02:53:30 EST 2005: but not thickness */

  /* need to set general orientation info too */

  }

  /* modify origin only if we flipped a spatial axis */

                          1.0,

  }

}

/*

**

** NOTE: this seems to destroy all space/orientation info.  What

** should be done?

*/

int

nrrdJoin(Nrrd *nout, const Nrrd *const *nin, unsigned int ninNum,

         unsigned int axis, int incrDim) {

  static const char me[]="nrrdJoin";

    permute[NRRD_DIM_MAX], ipermute[NRRD_DIM_MAX];

  char *dataPerm;

  Nrrd *ntmpperm,    /* axis-permuted version of output */

    **ninperm;

  airArray *mop;

  /* error checking */

  }

  }

    }

    }

  }

  }

  }

             me, diffdim);

  }

             me, axis, maxdim);

  }

  /* figure out dimension of output (outdim) */

    /* case A: (example) 2D slices and 3D slabs are being joined

       together to make a bigger 3D volume */

    outdim = maxdim;

    /* diffdim == 0 */

      /* case B: this is like the old "stitch": a bunch of equal-sized

         slices of dimension N are being stacked together to make an

         N+1 dimensional volume, which is essentially just the result of

         concatenating the memory of individual inputs */

      /* case C: axis < maxdim; maxdim == mindim */

      /* case C1 (!incrDim): a bunch of N-D slabs are being joined

         together to make a bigger N-D volume.  The axis along which

         they are being joined could be any of existing axes (from 0

         to maxdim-1) */

      /* case C2 (incrDim): this is also a "stitch", but the new axis

         created by the stitching is inserted into the existing

         axes. (ex: stitch 3 PGMs (R, G, B) together into a PPM (with

         color on axis zero) */

    }

  }

             me, outdim, NRRD_DIM_MAX);

  }

  /* do tacit reshaping, and possibly permuting, as needed */

                   ? ai

                      : axis));

    /* fprintf(stderr, "!%s: 1st permute[%d] = %d\n", me, ai, permute[ai]); */

  }

    /* fprintf(stderr, "!%s: ni = %d ---> diffdim = %d\n", me, ni, diffdim); */

      /* we do a tacit reshaping, which actually includes

         a tacit permuting, so we don't have to call permute

         on the parts that don't actually need it */

      /* NB: we register nrrdNix, not nrrdNuke */

      /* fprintf(stderr, "!%s: %d: tacit reshape/permute\n", me, ni); */

      }

      /* this may be done needlessly often */

        }

      }

      /* we don't have to actually call nrrdReshape(): we just nrrdWrap()

         the input data with the reshaped size array */

                 me, ni);

      }

                       (NRRD_AXIS_INFO_SIZE_BIT

                        /* HEY: this is being nixed because I can't think

                           of a sane way of keeping it consistent */

                        | NRRD_AXIS_INFO_SPACEDIRECTION_BIT));

      /* on this part, we permute (no need for a reshape) */

      }

    }

  }

  /* make sure all parts are compatible in type and shape,

     determine length of final output along axis (outlen) */

  outlen = 0;

    }

      }

    }

    }

    /* fprintf(stderr, "!%s: part %03d shape: ", me, ni); */

      /* fprintf(stderr, "%03u ", (unsigned int)ninperm[ni]->axis[ai].size);*/

      }

    }

    /*

    fprintf(stderr, "%03u\n", (unsigned int)ninperm[ni]->axis[outdim-1].size);

    */

  }

  /* fprintf(stderr, "!%s: outlen = %u\n", me, (unsigned int)outlen); */

  /* allocate temporary nrrd and concat input into it */

  outnum = 1;

    }

  }

  outnum *= size[outdim-1];

  }

    /* here is where the actual joining happens */

  }

  /* copy other axis-specific fields from nin[0] to ntmpperm */

  }

                   (NRRD_AXIS_INFO_NONE

                    /* HEY: this is being nixed because I can't think

                       of a sane way of keeping it consistent */

                    | NRRD_AXIS_INFO_SPACEDIRECTION_BIT));

  /* do the permutation required to get output in right order */

  }

  /* basic info is either already set or invalidated by joining */

  /* HEY: set content on output! */

}

/*

******** nrrdAxesSplit

**

** like reshape, but only for splitting one axis into a fast and slow part.

*/

int

nrrdAxesSplit(Nrrd *nout, const Nrrd *nin,

              unsigned int saxi, size_t sizeFast, size_t sizeSlow) {

  static const char me[]="nrrdAxesSplit", func[]="axsplit";

  unsigned int ai;

  }

             me, saxi, nin->dim-1);

  }

             me, NRRD_DIM_MAX);

  }

             "product of fast and slow sizes (%s * %s = %s)", me, saxi,

    return 1;

  }

                                 ? 0

                                 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT)))) {

    }

  }

                      NRRD_AXIS_INFO_NONE);

  }

  /* the ONLY thing we can say about the new axes are their sizes */

                        saxi, sizeFast, sizeSlow)) {

  }

  /* all basic information already copied by nrrdCopy */

}

/*

******** nrrdAxesDelete

**

** like reshape, but preserves axis information on old axes, and

** this is only for removing a "stub" axis with length 1.

*/

int

nrrdAxesDelete(Nrrd *nout, const Nrrd *nin, unsigned int daxi) {

  static const char me[]="nrrdAxesDelete", func[]="axdelete";

  unsigned int ai;

  }

  }

  }

    return 1;

  }

                                 ? 0

                                 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT)))) {

    }

  }

                      NRRD_AXIS_INFO_NONE);

  }

  }

  /* all basic information already copied by nrrdCopy */

}

/*

******** nrrdAxesMerge

**

** like reshape, but preserves axis information on old axes

** merges axis ax and ax+1 into one

*/

int

nrrdAxesMerge(Nrrd *nout, const Nrrd *nin, unsigned int maxi) {

  static const char me[]="nrrdAxesMerge", func[]="axmerge";

  unsigned int ai;

  size_t sizeFast, sizeSlow;

  }

  }

  }

                                 ? 0

                                 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT)))) {

    }

  }

                      NRRD_AXIS_INFO_NONE);

  }

  /* the ONLY thing we can say about the new axis is its size */

  }

  /* all basic information already copied by nrrdCopy */

}

/*

******** nrrdReshape_nva()

**

*/

int

nrrdReshape_nva(Nrrd *nout, const Nrrd *nin,

                unsigned int dim, const size_t *size) {

  static const char me[]="nrrdReshape_nva", func[]="reshape";

  size_t numOut;

  unsigned int ai;