/*

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

#include "privateTen.h"

#define INFO "Generate twisting helical tensor field"

static const char *_tend_helixInfoL =

  (INFO

   ". The main utility of such a field is to debug handling of coordinate "

   "systems in tensor field visualization.  The \"space directions\" and "

   "\"space origin\" fields of the NRRD header determines the mapping from "

   "coordinates in the index space of the image to coordinates in the "

   "world space in which the image is "

   "sampled.  The \"measurement frame\" field determines the mapping from "

   "the coordinates of the tensor itself, to coordinates of the world space. "

   "When these are correctly handled, the "

   "region of high anisotropy is a right-handed helix (same as DNA). "

   "Using differing axes sizes (via \"-s\") helps make sure that the "

   "raster ordering of axes is correct.  In addition, the tensors twist "

   "relative to the helix, which exposes handling of the measurement frame. "

   "If you trace paths guided by the principal eigenvector of the tensors, "

   "along the surface of the helical cylinder, you get another "

   "right-handed helix, as if the the tensor field is modeling the result "

   "if twisting a set of fibers into single-stranded helical bundle. ");

void

tend_helixDoit(Nrrd *nout, double bnd,

               double orig[3], double i2w[9], double mf[9],

               double r, double R, double S, double angle, int incrtwist,

               double ev[3], double bgEval, int verbose) {

  int sx, sy, sz, xi, yi, zi;

  double th, t0, t1, t2, t3, v1, v2,

    wpos[3], vpos[3], mfT[9],

    W2H[9], H2W[9], H2C[9], C2H[9], fv[3], rv[3], uv[3], mA[9], mB[9], inside,

    tmp[3], len;

  float *out;

    }

#define WPOS(pos, th) ELL_3V_SET((pos),R*cos(th), R*sin(th), S*(th)/(2*AIR_PI))

#define VAL(th) (WPOS(wpos, th), ELL_3V_DIST(wpos, vpos))

#define RR 0.61803399

#define CC (1.0-RR)

#define SHIFT3(a,b,c,d) (a)=(b); (b)=(c); (c)=(d)

#define SHIFT2(a,b,c)   (a)=(b); (b)=(c)

        }

          }

        }

        /* t1 (and t2) are now the th for which the point on the helix

           (R*cos(th), R*sin(th), S*(th)/(2*AIR_PI)) is closest to vpos */

        ELL_3V_COPY(rv, wpos);

                                      pointing to vpos */

        ELL_3MV_ROW0_SET(W2H, uv); /* as is (uv,rv,fv) */

        ELL_3MV_ROW1_SET(W2H, rv);

        ELL_3MV_ROW2_SET(W2H, fv);

        ELL_3M_TRANSPOSE(H2W, W2H);

          th = angle;

        }

        ELL_3M_TRANSPOSE(C2H, H2C);

                         AIR_LERP(inside, bgEval, ev[1]),

                         AIR_LERP(inside, bgEval, ev[2]),

                         AIR_LERP(inside, bgEval, ev[0]));

      }

    }

  }

  return;

int

tend_helixMain(int argc, const char **argv, const char *me,

               hestParm *hparm) {

  int pret;

  airArray *mop;

  Nrrd *nout;

    orig[3], i2w[9], rot[9], mf[9], spd[4][3], bge;

             "sizes along fast, medium, and slow axes of the sampled volume, "

             "often called \"X\", \"Y\", and \"Z\".  It is best to use "

             "slightly different sizes here, to expose errors in interpreting "

             "axis ordering (e.g. \"-s 39 40 41\")");

             "0 0 0",

             "quaternion quotient space orientation of image");

             "0 0 0",

             "quaternion quotient space orientation of measurement frame");

             "parameter governing how fuzzy the boundary between high and "

             "low anisotropy is. Use \"-b 0\" for no fuzziness");

             "(minor) radius of cylinder tracing helix");

             "(major) radius of helical turns");

             "spacing between turns of helix (along its axis)");

             "maximal angle of twist of tensors along path.  There is no "

             "twist at helical core of path, and twist increases linearly "

             "with radius around this path.  Positive twist angle with "

             "positive spacing resulting in a right-handed twist around a "

             "right-handed helix. ");

             "changes behavior of twist angle as function of distance from "

             "center of helical core: instead of increasing linearly as "

             "describe above, be at a constant angle");

             "0.006 0.002 0.001",

             "eigenvalues of tensors (in order) along direction of coil, "

             "circumferential around coil, and radial around coil. ");

             "eigenvalue of isotropic background");

             "verbose output");

             "output file");

                        AIR_CAST(size_t, 7),

  }

             -2*R + 2*R/size[0],

             -2*R + 2*R/size[1],

             -2*R + 2*R/size[2]);

                     nrrdCenterUnknown, nrrdCenterCell,

                     nrrdCenterCell, nrrdCenterCell);

                     nrrdKind3DMaskedSymMatrix, nrrdKindSpace,

                     nrrdKindSpace, nrrdKindSpace);

  }

TEND_CMD(helix, INFO);