Head

GCC Code Coverage Report

Directory:	./		Exec	Total	Coverage
File:	src/ten/triple.c	Lines:	0	207	0.0 %
Date:	2017-05-26	Branches:	0	108	0.0 %


/*
  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"

typedef void (*tenTripleConverter)(double dst[3], const double src[3]);

#define SQRT6 2.449489742783178098197284
#define SQRT2 1.414213562373095048801689
#define SQRT3 1.732050807568877293527446

#define J1 (j[0])
#define J2 (j[1])
#define J3 (j[2])

#define K1 (k[0])
#define K2 (k[1])
#define K3 (k[2])

#define R1 (r[0])
#define R2 (r[1])
#define R3 (r[2])

#define MU1 (mu[0])
#define MU2 (mu[1])
#define MU3 (mu[2])

static void
_iden(double dst[3], const double src[3]) {
  ELL_3V_COPY(dst, src);
  return;
}

/* in the function names below, the format is _<dst>_<src>() */
static void
_mu_ev(double mu[3], const double ev[3]) {
  double mm;

  mm = mu[0] = (ev[0] + ev[1] + ev[2])/3;
  mu[1] = ((ev[0] - mm)*(ev[0] - mm) +
           (ev[1] - mm)*(ev[1] - mm) +
           (ev[2] - mm)*(ev[2] - mm))/3;
  mu[2] = ((ev[0] - mm)*(ev[0] - mm)*(ev[0] - mm) +
           (ev[1] - mm)*(ev[1] - mm)*(ev[1] - mm) +
           (ev[2] - mm)*(ev[2] - mm)*(ev[2] - mm))/3;
}

static double
_xyzmat[] = {2/SQRT6, -1/SQRT6, -1/SQRT6,
             0,        1/SQRT2, -1/SQRT2,
             1/SQRT3,  1/SQRT3,  1/SQRT3};

static void
_xyz_ev(double xyz[3], const double _ev[3]) {
  double ev[3], tmp;

  ELL_3V_COPY(ev, _ev);
  ELL_SORT3(ev[0], ev[1], ev[2], tmp);
  ELL_3MV_MUL(xyz, _xyzmat, ev);
}

static void
_ev_xyz(double ev[3], const double xyz[3]) {

  ELL_3MV_TMUL(ev, _xyzmat, xyz);
}

static void
_j_ev(double j[3], const double ev[3]) {

  J1 = ev[0] + ev[1] + ev[2];
  J2 = ev[0]*ev[1] + ev[0]*ev[2] + ev[1]*ev[2];
  J3 = ev[0]*ev[1]*ev[2];
}

static void
_k_mu(double k[3], const double mu[3]) {
  double stdv;

  K1 = 3*MU1;
  stdv = sqrt(MU2);
  K2 = SQRT3*stdv;
  K3 = stdv ? SQRT2*MU3/(stdv*stdv*stdv) : 0;
}

static void
_r_ev(double r[3], const double ev[3]) {
  double mu[3], stdv;

  _mu_ev(mu, ev);
  R1 = sqrt(ev[0]*ev[0] + ev[1]*ev[1] + ev[2]*ev[2]);
  stdv = sqrt(MU2);
  R2 = R1 ? (3/SQRT2)*stdv/R1 : 0;
  R3 = stdv ? SQRT2*MU3/(stdv*stdv*stdv) : 0;
}

static void
_r_mu(double r[3], const double mu[3]) {
  double stdv;

  R1 = sqrt(3*(MU1*MU1 + MU2));
  stdv = sqrt(MU2);
  R2 = R1 ? (3/SQRT2)*stdv/R1 : 0;
  R3 = stdv ? SQRT2*MU3/(stdv*stdv*stdv) : 0;
}

static void
_ev_wp(double ev[3], const double wp[3]) {

  ev[0] = wp[0] + wp[1]*cos(wp[2]);
  ev[1] = wp[0] + wp[1]*cos(wp[2] - 2*AIR_PI/3);
  ev[2] = wp[0] + wp[1]*cos(wp[2] + 2*AIR_PI/3);
}

static void
_wp_mu(double wp[3], const double mu[3]) {
  double stdv, mode;

  wp[0] = MU1;
  stdv = sqrt(MU2);
  wp[1] = SQRT2*stdv;
  mode = stdv ? SQRT2*MU3/(stdv*stdv*stdv) : 0;
  mode = AIR_CLAMP(-1, mode, 1);
  wp[2] = acos(AIR_CLAMP(-1, mode, 1))/3;
}

static void
_mu_j(double mu[3], const double j[3]) {

  MU1 = J1/3;
  MU2 = 2*(J1*J1 - 3*J2)/9;
  MU3 = 2*J1*J1*J1/27 - J1*J2/3 + J3;
}

static void
_r_j(double r[3], const double j[3]) {
  double mu[3], stdv;

  R1 = sqrt(J1*J1 - 2*J2);
  R2 = sqrt(J1*J1 - 3*J2)/R1;
  _mu_j(mu, j);
  stdv = sqrt(MU2);
  R3 = stdv ? SQRT2*MU3/(stdv*stdv*stdv) : 0;
}

static void
_k_r(double k[3], const double r[3]) {

  K1 = R1*sqrt(3 - 2*R2*R2);
  K2 = (SQRT2/SQRT3)*R1*R2;
  K3 = R3;
}

/*
_j_r(double j[3], const double r[3]) {
  double ss, nmu3;

  J1 = R1*sqrt(3 - 2*R2*R2);
  J2 = R1*R1*(1 - R2*R2);
  ss = R1*R2;
  nmu3 = 2*R3*ss*ss*ss;
  J3 =
}
*/

static void
_wp_k(double wp[3], const double k[3]) {

  wp[0] = K1/3;
  wp[1] = (SQRT2/SQRT3)*K2;
  wp[2] = acos(AIR_CLAMP(-1, K3, 1))/3;
}

static void
_k_wp(double k[3], const double wp[3]) {

  K1 = 3*wp[0];
  K2 = (SQRT3/SQRT2)*wp[1];
  K3 = cos(3*wp[2]);
}

static void
_rtz_xyz(double rThZ[3], const double XYZ[3]) {

  rThZ[0] = sqrt(XYZ[0]*XYZ[0] + XYZ[1]*XYZ[1]);
  rThZ[1] = atan2(XYZ[1], XYZ[0]);
  rThZ[2] = XYZ[2];
}

static void
_rtp_xyz(double RThPh[3], const double XYZ[3]) {

  RThPh[0] = sqrt(XYZ[0]*XYZ[0] + XYZ[1]*XYZ[1] + XYZ[2]*XYZ[2]);
  RThPh[1] = atan2(XYZ[1], XYZ[0]);
  RThPh[2] = atan2(sqrt(XYZ[0]*XYZ[0] + XYZ[1]*XYZ[1]), XYZ[2]);
}

static void
_xyz_rtz(double XYZ[3], const double rThZ[3]) {

  XYZ[0] = rThZ[0]*cos(rThZ[1]);
  XYZ[1] = rThZ[0]*sin(rThZ[1]);
  XYZ[2] = rThZ[2];
}

static void
_xyz_rtp(double XYZ[3], const double RThPh[3]) {

  XYZ[0] = RThPh[0]*cos(RThPh[1])*sin(RThPh[2]);
  XYZ[1] = RThPh[0]*sin(RThPh[1])*sin(RThPh[2]);
  XYZ[2] = RThPh[0]*cos(RThPh[2]);
}

static void
_rtz_k(double rThZ[3], const double k[3]) {

  rThZ[0] = K2;
  rThZ[1] = acos(AIR_CLAMP(-1, K3, 1))/3;
  rThZ[2] = K1/SQRT3;
}

static void
_k_rtz(double k[3], const double rThZ[3]) {

  K1 = SQRT3*rThZ[2];
  K2 = rThZ[0];
  K3 = cos(3*rThZ[1]);
}

static void
_rtp_r(double RThPh[3], const double r[3]) {

  RThPh[0] = R1;
  RThPh[1] = acos(AIR_CLAMP(-1, R3, 1))/3;
  RThPh[2] = asin(AIR_CLAMP(-1, (SQRT2/SQRT3)*R2, 1));
}

static void
_r_rtp(double r[3], const double RThPh[3]) {

  R1 = RThPh[0];
  R2 = sin(RThPh[2])*SQRT3/SQRT2;
  R3 = cos(3*RThPh[1]);
}

static void
_wp_rtz(double wp[3], const double rThZ[3]) {

  wp[0] = rThZ[2]/SQRT3;
  wp[1] = (SQRT2/SQRT3)*rThZ[0];
  wp[2] = rThZ[1];
}

static void
_rtz_wp(double rThZ[3], const double wp[3]) {

  rThZ[0] = (SQRT3/SQRT2)*wp[1];
  rThZ[1] = wp[2];
  rThZ[2] = SQRT3*wp[0];
}

#define CONVERT1(dst, mid, src) \
static void \
_##dst##_##src(double dst[3], const double src[3]) { \
  double mid[3]; \
  _##mid##_##src(mid, src); \
  _##dst##_##mid(dst, mid); \
}

#define CONVERT2(dst, mdB, mdA, src) \
static void \
_##dst##_##src(double dst[3], const double src[3]) { \
  double mdA[3], mdB[3];  \
  _##mdA##_##src(mdA, src); \
  _##mdB##_##mdA(mdB, mdA); \
  _##dst##_##mdB(dst, mdB); \
}

CONVERT1(ev, xyz, rtz)    /* _ev_rtz */
CONVERT1(rtz, xyz, ev)    /* _rtz_ev */
CONVERT1(ev, xyz, rtp)    /* _ev_rtp */
CONVERT1(rtp, xyz, ev)    /* _rtp_ev */

CONVERT1(k, mu, ev)       /* _k_ev */
CONVERT1(wp, mu, ev)      /* _wp_ev */

CONVERT1(ev, wp, mu)      /* _ev_mu */

CONVERT2(ev, wp, mu, j)   /* _ev_j */
CONVERT1(ev, wp, k)       /* _ev_k */

CONVERT2(ev, xyz, rtp, r) /* _ev_r */

static tenTripleConverter
_convert[TEN_TRIPLE_TYPE_MAX+1][TEN_TRIPLE_TYPE_MAX+1] = {
  {NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL},
  /* DEST:    SRC:  ev      mu      xyz       rtz       rtp       J      K       R       WP */
  /* ev */  {NULL, _iden,   _ev_mu, _ev_xyz,  _ev_rtz,  _ev_rtp,  _ev_j, _ev_k,  _ev_r,  _ev_wp},
  /* mu */  {NULL, _mu_ev,  _iden,  NULL,     NULL,     NULL,     _mu_j, NULL,   NULL,   NULL},
  /* xyz */ {NULL, _xyz_ev, NULL,   _iden,    _xyz_rtz, _xyz_rtp, NULL,  NULL,   NULL,   NULL},
  /* rtz */ {NULL, _rtz_ev, NULL,   _rtz_xyz, _iden,    NULL,     NULL,  _rtz_k, NULL,   _rtz_wp},
  /* rtp */ {NULL, _rtp_ev, NULL,   _rtp_xyz, NULL,     _iden,    NULL,  NULL,   _rtp_r, NULL},
  /* J */   {NULL, _j_ev,   NULL,   NULL,     NULL,     NULL,     _iden, NULL,   NULL,   NULL},
  /* K */   {NULL, _k_ev,   _k_mu,  NULL,     _k_rtz,   NULL,     NULL,  _iden,  _k_r,   _k_wp},
  /* R */   {NULL, _r_ev,   _r_mu,  NULL,     NULL,     _r_rtp,   _r_j,  NULL,   _iden,  NULL},
  /* WP */  {NULL, _wp_ev,  _wp_mu, NULL,     _wp_rtz,  NULL,     NULL,  _wp_k,   NULL,  _iden}};

void
tenTripleConvertSingle_d(double dst[3], int dstType,
                         const double src[3], const int srcType) {
  static const char me[]="tenTripleConvertSingle_d";
  int direct;

  if (airEnumValCheck(tenTripleType, dstType)
      || airEnumValCheck(tenTripleType, srcType)) {
    /* got invalid source or destination type */
    ELL_3V_SET(dst, AIR_NAN, AIR_NAN, AIR_NAN);
    return;
  }

  if (_convert[dstType][srcType]) {
    /* we have a direct converter */
    _convert[dstType][srcType](dst, src);
    direct = AIR_TRUE;
  } else {
    double eval[3];
    /* else, for lack of anything clever, we convert via evals */
    _convert[tenTripleTypeEigenvalue][srcType](eval, src);
    _convert[dstType][tenTripleTypeEigenvalue](dst, eval);
    direct = AIR_FALSE;
  }

  /* warn if conversion created non-existent values from
     existent input */
  if (ELL_3V_EXISTS(src) && !ELL_3V_EXISTS(dst)) {
    fprintf(stderr, "%s: problem? (%s) %g %g %g <-%s- (%s) %g %g %g\n", me,
            airEnumStr(tenTripleType, dstType),
            dst[0], dst[1], dst[2],
            direct ? "-" : "...",
            airEnumStr(tenTripleType, srcType),
            src[0], src[1], src[2]);
  }

  return;
}

void
tenTripleConvertSingle_f(float _dst[3], int dstType,
                         const float _src[3], const int srcType) {
  double dst[3], src[3];

  ELL_3V_COPY(src, _src);
  tenTripleConvertSingle_d(dst, dstType, src, srcType);
  ELL_3V_COPY_TT(_dst, float, dst);
}

void
tenTripleCalcSingle_d(double dst[3], int ttype, double ten[7]) {
  double eval[3];

  /* in time this can become more efficient ... */
  switch (ttype) {
  case tenTripleTypeEigenvalue:
    tenEigensolve_d(dst, NULL, ten);
    break;
  case tenTripleTypeMoment:
  case tenTripleTypeXYZ:
  case tenTripleTypeRThetaZ:
  case tenTripleTypeRThetaPhi:
  case tenTripleTypeK:
  case tenTripleTypeJ:
  case tenTripleTypeWheelParm:
    tenEigensolve_d(eval, NULL, ten);
    tenTripleConvertSingle_d(dst, ttype, eval, tenTripleTypeEigenvalue);
    break;
  case tenTripleTypeR:
    dst[0] = sqrt(_tenAnisoTen_d[tenAniso_S](ten));
    dst[1] = _tenAnisoTen_d[tenAniso_FA](ten);
    dst[2] = _tenAnisoTen_d[tenAniso_Mode](ten);
    break;
  default:
    /* what on earth? */
    ELL_3V_SET(dst, AIR_NAN, AIR_NAN, AIR_NAN);
  }
  return;
}

void
tenTripleCalcSingle_f(float dst[3], int ttype, float ten[7]) {
  double dst_d[3], ten_d[7];

  TEN_T_COPY(ten_d, ten);
  tenTripleCalcSingle_d(dst_d, ttype, ten_d);
  ELL_3V_COPY_TT(dst, float, dst_d);
  return;
}

int
tenTripleCalc(Nrrd *nout, int ttype, const Nrrd *nten) {
  static const char me[]="tenTripleCalc";
  size_t II, NN, size[NRRD_DIM_MAX];
  double (*ins)(void *, size_t, double), (*lup)(const void *, size_t);

  if (!( nout && nten )) {
    biffAddf(TEN, "%s: got NULL pointer", me);
    return 1;
  }
  if (airEnumValCheck(tenTripleType, ttype)) {
    biffAddf(TEN, "%s: got invalid %s (%d)", me,
             tenTripleType->name, ttype);
    return 1;
  }
  if (tenTensorCheck(nten, nrrdTypeDefault, AIR_FALSE, AIR_TRUE)) {
    biffAddf(TEN, "%s: didn't get a valid DT array", me);
    return 1;
  }
  if (!( nrrdTypeFloat == nten->type ||
         nrrdTypeDouble == nten->type )) {
    biffAddf(TEN, "%s: need input type %s or %s, not %s\n", me,
             airEnumStr(nrrdType, nrrdTypeFloat),
             airEnumStr(nrrdType, nrrdTypeFloat),
             airEnumStr(nrrdType, nten->type));
  }

  nrrdAxisInfoGet_nva(nten, nrrdAxisInfoSize, size);
  size[0] = 3;
  if (nrrdMaybeAlloc_nva(nout, nten->type, nten->dim, size)) {
    biffMovef(TEN, NRRD, "%s: couldn't alloc output", me);
    return 1;
  }

  NN = nrrdElementNumber(nten)/7;
  lup = nrrdDLookup[nten->type];
  ins = nrrdDInsert[nten->type];
  for (II=0; II<NN; II++) {
    double ten[7], trip[3];
    unsigned int vv;
    for (vv=0; vv<7; vv++) {
      ten[vv] = lup(nten->data, vv + 7*II);
    }
    tenTripleCalcSingle_d(trip, ttype, ten);
    for (vv=0; vv<3; vv++) {
      ins(nout->data, vv + 3*II, trip[vv]);
    }
  }
  if (nrrdAxisInfoCopy(nout, nten, NULL, (NRRD_AXIS_INFO_SIZE_BIT))) {
    biffMovef(TEN, NRRD, "%s: couldn't copy axis info", me);
    return 1;
  }
  nout->axis[0].kind = nrrdKindUnknown;
  if (nrrdBasicInfoCopy(nout, nten,
                        NRRD_BASIC_INFO_ALL ^ NRRD_BASIC_INFO_SPACE)) {
    biffAddf(TEN, "%s:", me);
    return 1;
  }

  return 0;
}

int
tenTripleConvert(Nrrd *nout, int dstType,
                 const Nrrd *nin, int srcType) {
  static const char me[]="tenTripleConvert";
  size_t II, NN;
  double (*ins)(void *, size_t, double), (*lup)(const void *, size_t);

  if (!( nout && nin )) {
    biffAddf(TEN, "%s: got NULL pointer", me);
    return 1;
  }
  if ( airEnumValCheck(tenTripleType, dstType) ||
       airEnumValCheck(tenTripleType, srcType) ) {
    biffAddf(TEN, "%s: got invalid %s dst (%d) or src (%d)", me,
             tenTripleType->name, dstType, srcType);
    return 1;
  }
  if (3 != nin->axis[0].size) {
    char stmp[AIR_STRLEN_SMALL];
    biffAddf(TEN, "%s: need axis[0].size 3, not %s", me,
             airSprintSize_t(stmp, nin->axis[0].size));
    return 1;
  }
  if (nrrdTypeBlock == nin->type) {
    biffAddf(TEN, "%s: input has non-scalar %s type",
             me, airEnumStr(nrrdType, nrrdTypeBlock));
    return 1;
  }

  if (nrrdCopy(nout, nin)) {
    biffMovef(TEN, NRRD, "%s: couldn't initialize output", me);
    return 1;
  }
  lup = nrrdDLookup[nin->type];
  ins = nrrdDInsert[nout->type];
  NN = nrrdElementNumber(nin)/3;
  for (II=0; II<NN; II++) {
    double src[3], dst[3];
    src[0] = lup(nin->data, 0 + 3*II);
    src[1] = lup(nin->data, 1 + 3*II);
    src[2] = lup(nin->data, 2 + 3*II);
    tenTripleConvertSingle_d(dst, dstType, src, srcType);
    ins(nout->data, 0 + 3*II, dst[0]);
    ins(nout->data, 1 + 3*II, dst[1]);
    ins(nout->data, 2 + 3*II, dst[2]);
  }

  return 0;
}


Generated by: GCOVR (Version 3.3)

Line	Branch	Exec	Source
1			/*
2			Teem: Tools to process and visualize scientific data and images .
3			Copyright (C) 2013, 2012, 2011, 2010, 2009 University of Chicago
4			Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
5			Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
6
7			This library is free software; you can redistribute it and/or
8			modify it under the terms of the GNU Lesser General Public License
9			(LGPL) as published by the Free Software Foundation; either
10			version 2.1 of the License, or (at your option) any later version.
11			The terms of redistributing and/or modifying this software also
12			include exceptions to the LGPL that facilitate static linking.
13
14			This library is distributed in the hope that it will be useful,
15			but WITHOUT ANY WARRANTY; without even the implied warranty of
16			MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17			Lesser General Public License for more details.
18
19			You should have received a copy of the GNU Lesser General Public License
20			along with this library; if not, write to Free Software Foundation, Inc.,
21			51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22			*/
23
24
25			#include "ten.h"
26			#include "privateTen.h"
27
28			typedef void (*tenTripleConverter)(double dst[3], const double src[3]);
29
30			#define SQRT6 2.449489742783178098197284
31			#define SQRT2 1.414213562373095048801689
32			#define SQRT3 1.732050807568877293527446
33
34			#define J1 (j[0])
35			#define J2 (j[1])
36			#define J3 (j[2])
37
38			#define K1 (k[0])
39			#define K2 (k[1])
40			#define K3 (k[2])
41
42			#define R1 (r[0])
43			#define R2 (r[1])
44			#define R3 (r[2])
45
46			#define MU1 (mu[0])
47			#define MU2 (mu[1])
48			#define MU3 (mu[2])
49
50			static void
51			_iden(double dst[3], const double src[3]) {
52			ELL_3V_COPY(dst, src);
53			return;
54			}
55
56			/* in the function names below, the format is _<dst>_<src>() */
57			static void
58			_mu_ev(double mu[3], const double ev[3]) {
59			double mm;
60
61			mm = mu[0] = (ev[0] + ev[1] + ev[2])/3;
62			mu[1] = ((ev[0] - mm)*(ev[0] - mm) +
63			(ev[1] - mm)*(ev[1] - mm) +
64			(ev[2] - mm)*(ev[2] - mm))/3;
65			mu[2] = ((ev[0] - mm)(ev[0] - mm)(ev[0] - mm) +
66			(ev[1] - mm)(ev[1] - mm)(ev[1] - mm) +
67			(ev[2] - mm)(ev[2] - mm)(ev[2] - mm))/3;
68			}
69
70			static double
71			_xyzmat[] = {2/SQRT6, -1/SQRT6, -1/SQRT6,
72			0, 1/SQRT2, -1/SQRT2,
73			1/SQRT3, 1/SQRT3, 1/SQRT3};
74
75			static void
76			_xyz_ev(double xyz[3], const double _ev[3]) {
77			double ev[3], tmp;
78
79			ELL_3V_COPY(ev, _ev);
80			ELL_SORT3(ev[0], ev[1], ev[2], tmp);
81			ELL_3MV_MUL(xyz, _xyzmat, ev);
82			}
83
84			static void
85			_ev_xyz(double ev[3], const double xyz[3]) {
86
87			ELL_3MV_TMUL(ev, _xyzmat, xyz);
88			}
89
90			static void
91			_j_ev(double j[3], const double ev[3]) {
92
93			J1 = ev[0] + ev[1] + ev[2];
94			J2 = ev[0]ev[1] + ev[0]ev[2] + ev[1]*ev[2];
95			J3 = ev[0]ev[1]ev[2];
96			}
97
98			static void
99			_k_mu(double k[3], const double mu[3]) {
100			double stdv;
101
102			K1 = 3*MU1;
103			stdv = sqrt(MU2);
104			K2 = SQRT3*stdv;
105			K3 = stdv ? SQRT2MU3/(stdvstdv*stdv) : 0;
106			}
107
108			static void
109			_r_ev(double r[3], const double ev[3]) {
110			double mu[3], stdv;
111
112			_mu_ev(mu, ev);
113			R1 = sqrt(ev[0]ev[0] + ev[1]ev[1] + ev[2]*ev[2]);
114			stdv = sqrt(MU2);
115			R2 = R1 ? (3/SQRT2)*stdv/R1 : 0;
116			R3 = stdv ? SQRT2MU3/(stdvstdv*stdv) : 0;
117			}
118
119			static void
120			_r_mu(double r[3], const double mu[3]) {
121			double stdv;
122
123			R1 = sqrt(3(MU1MU1 + MU2));
124			stdv = sqrt(MU2);
125			R2 = R1 ? (3/SQRT2)*stdv/R1 : 0;
126			R3 = stdv ? SQRT2MU3/(stdvstdv*stdv) : 0;
127			}
128
129			static void
130			_ev_wp(double ev[3], const double wp[3]) {
131
132			ev[0] = wp[0] + wp[1]*cos(wp[2]);
133			ev[1] = wp[0] + wp[1]cos(wp[2] - 2AIR_PI/3);
134			ev[2] = wp[0] + wp[1]cos(wp[2] + 2AIR_PI/3);
135			}
136
137			static void
138			_wp_mu(double wp[3], const double mu[3]) {
139			double stdv, mode;
140
141			wp[0] = MU1;
142			stdv = sqrt(MU2);
143			wp[1] = SQRT2*stdv;
144			mode = stdv ? SQRT2MU3/(stdvstdv*stdv) : 0;
145			mode = AIR_CLAMP(-1, mode, 1);
146			wp[2] = acos(AIR_CLAMP(-1, mode, 1))/3;
147			}
148
149			static void
150			_mu_j(double mu[3], const double j[3]) {
151
152			MU1 = J1/3;
153			MU2 = 2(J1J1 - 3*J2)/9;
154			MU3 = 2J1J1J1/27 - J1J2/3 + J3;
155			}
156
157			static void
158			_r_j(double r[3], const double j[3]) {
159			double mu[3], stdv;
160
161			R1 = sqrt(J1J1 - 2J2);
162			R2 = sqrt(J1J1 - 3J2)/R1;
163			_mu_j(mu, j);
164			stdv = sqrt(MU2);
165			R3 = stdv ? SQRT2MU3/(stdvstdv*stdv) : 0;
166			}
167
168			static void
169			_k_r(double k[3], const double r[3]) {
170
171			K1 = R1sqrt(3 - 2R2*R2);
172			K2 = (SQRT2/SQRT3)R1R2;
173			K3 = R3;
174			}
175
176			/*
177			_j_r(double j[3], const double r[3]) {
178			double ss, nmu3;
179
180			J1 = R1sqrt(3 - 2R2*R2);
181			J2 = R1R1(1 - R2*R2);
182			ss = R1*R2;
183			nmu3 = 2R3ssssss;
184			J3 =
185			}
186			*/
187
188			static void
189			_wp_k(double wp[3], const double k[3]) {
190
191			wp[0] = K1/3;
192			wp[1] = (SQRT2/SQRT3)*K2;
193			wp[2] = acos(AIR_CLAMP(-1, K3, 1))/3;
194			}
195
196			static void
197			_k_wp(double k[3], const double wp[3]) {
198
199			K1 = 3*wp[0];
200			K2 = (SQRT3/SQRT2)*wp[1];
201			K3 = cos(3*wp[2]);
202			}
203
204			static void
205			_rtz_xyz(double rThZ[3], const double XYZ[3]) {
206
207			rThZ[0] = sqrt(XYZ[0]XYZ[0] + XYZ[1]XYZ[1]);
208			rThZ[1] = atan2(XYZ[1], XYZ[0]);
209			rThZ[2] = XYZ[2];
210			}
211
212			static void
213			_rtp_xyz(double RThPh[3], const double XYZ[3]) {
214
215			RThPh[0] = sqrt(XYZ[0]XYZ[0] + XYZ[1]XYZ[1] + XYZ[2]*XYZ[2]);
216			RThPh[1] = atan2(XYZ[1], XYZ[0]);
217			RThPh[2] = atan2(sqrt(XYZ[0]XYZ[0] + XYZ[1]XYZ[1]), XYZ[2]);
218			}
219
220			static void
221			_xyz_rtz(double XYZ[3], const double rThZ[3]) {
222
223			XYZ[0] = rThZ[0]*cos(rThZ[1]);
224			XYZ[1] = rThZ[0]*sin(rThZ[1]);
225			XYZ[2] = rThZ[2];
226			}
227
228			static void
229			_xyz_rtp(double XYZ[3], const double RThPh[3]) {
230
231			XYZ[0] = RThPh[0]cos(RThPh[1])sin(RThPh[2]);
232			XYZ[1] = RThPh[0]sin(RThPh[1])sin(RThPh[2]);
233			XYZ[2] = RThPh[0]*cos(RThPh[2]);
234			}
235
236			static void
237			_rtz_k(double rThZ[3], const double k[3]) {
238
239			rThZ[0] = K2;
240			rThZ[1] = acos(AIR_CLAMP(-1, K3, 1))/3;
241			rThZ[2] = K1/SQRT3;
242			}
243
244			static void
245			_k_rtz(double k[3], const double rThZ[3]) {
246
247			K1 = SQRT3*rThZ[2];
248			K2 = rThZ[0];
249			K3 = cos(3*rThZ[1]);
250			}
251
252			static void
253			_rtp_r(double RThPh[3], const double r[3]) {
254
255			RThPh[0] = R1;
256			RThPh[1] = acos(AIR_CLAMP(-1, R3, 1))/3;
257			RThPh[2] = asin(AIR_CLAMP(-1, (SQRT2/SQRT3)*R2, 1));
258			}
259
260			static void
261			_r_rtp(double r[3], const double RThPh[3]) {
262
263			R1 = RThPh[0];
264			R2 = sin(RThPh[2])*SQRT3/SQRT2;
265			R3 = cos(3*RThPh[1]);
266			}
267
268			static void
269			_wp_rtz(double wp[3], const double rThZ[3]) {
270
271			wp[0] = rThZ[2]/SQRT3;
272			wp[1] = (SQRT2/SQRT3)*rThZ[0];
273			wp[2] = rThZ[1];
274			}
275
276			static void
277			_rtz_wp(double rThZ[3], const double wp[3]) {
278
279			rThZ[0] = (SQRT3/SQRT2)*wp[1];
280			rThZ[1] = wp[2];
281			rThZ[2] = SQRT3*wp[0];
282			}
283
284			#define CONVERT1(dst, mid, src) \
285			static void \
286			_##dst##_##src(double dst[3], const double src[3]) { \
287			double mid[3]; \
288			_##mid##_##src(mid, src); \
289			_##dst##_##mid(dst, mid); \
290			}
291
292			#define CONVERT2(dst, mdB, mdA, src) \
293			static void \
294			_##dst##_##src(double dst[3], const double src[3]) { \
295			double mdA[3], mdB[3]; \
296			_##mdA##_##src(mdA, src); \
297			_##mdB##_##mdA(mdB, mdA); \
298			_##dst##_##mdB(dst, mdB); \
299			}
300
301			CONVERT1(ev, xyz, rtz) /* _ev_rtz */
302			CONVERT1(rtz, xyz, ev) /* _rtz_ev */
303			CONVERT1(ev, xyz, rtp) /* _ev_rtp */
304			CONVERT1(rtp, xyz, ev) /* _rtp_ev */
305
306			CONVERT1(k, mu, ev) /* _k_ev */
307			CONVERT1(wp, mu, ev) /* _wp_ev */
308
309			CONVERT1(ev, wp, mu) /* _ev_mu */
310
311			CONVERT2(ev, wp, mu, j) /* _ev_j */
312			CONVERT1(ev, wp, k) /* _ev_k */
313
314			CONVERT2(ev, xyz, rtp, r) /* _ev_r */
315
316			static tenTripleConverter
317			_convert[TEN_TRIPLE_TYPE_MAX+1][TEN_TRIPLE_TYPE_MAX+1] = {
318			{NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL},
319			/* DEST: SRC: ev mu xyz rtz rtp J K R WP */
320			/* ev */ {NULL, _iden, _ev_mu, _ev_xyz, _ev_rtz, _ev_rtp, _ev_j, _ev_k, _ev_r, _ev_wp},
321			/* mu */ {NULL, _mu_ev, _iden, NULL, NULL, NULL, _mu_j, NULL, NULL, NULL},
322			/* xyz */ {NULL, _xyz_ev, NULL, _iden, _xyz_rtz, _xyz_rtp, NULL, NULL, NULL, NULL},
323			/* rtz */ {NULL, _rtz_ev, NULL, _rtz_xyz, _iden, NULL, NULL, _rtz_k, NULL, _rtz_wp},
324			/* rtp */ {NULL, _rtp_ev, NULL, _rtp_xyz, NULL, _iden, NULL, NULL, _rtp_r, NULL},
325			/* J */ {NULL, _j_ev, NULL, NULL, NULL, NULL, _iden, NULL, NULL, NULL},
326			/* K */ {NULL, _k_ev, _k_mu, NULL, _k_rtz, NULL, NULL, _iden, _k_r, _k_wp},
327			/* R */ {NULL, _r_ev, _r_mu, NULL, NULL, _r_rtp, _r_j, NULL, _iden, NULL},
328			/* WP */ {NULL, _wp_ev, _wp_mu, NULL, _wp_rtz, NULL, NULL, _wp_k, NULL, _iden}};
329
330			void
331			tenTripleConvertSingle_d(double dst[3], int dstType,
332			const double src[3], const int srcType) {
333			static const char me[]="tenTripleConvertSingle_d";
334			int direct;
335
336			if (airEnumValCheck(tenTripleType, dstType)
337			\|\| airEnumValCheck(tenTripleType, srcType)) {
338			/* got invalid source or destination type */
339			ELL_3V_SET(dst, AIR_NAN, AIR_NAN, AIR_NAN);
340			return;
341			}
342
343			if (_convert[dstType][srcType]) {
344			/* we have a direct converter */
345			_convert[dstType][srcType](dst, src);
346			direct = AIR_TRUE;
347			} else {
348			double eval[3];
349			/* else, for lack of anything clever, we convert via evals */
350			_convert[tenTripleTypeEigenvalue][srcType](eval, src);
351			_convert[dstType][tenTripleTypeEigenvalue](dst, eval);
352			direct = AIR_FALSE;
353			}
354
355			/* warn if conversion created non-existent values from
356			existent input */
357			if (ELL_3V_EXISTS(src) && !ELL_3V_EXISTS(dst)) {
358			fprintf(stderr, "%s: problem? (%s) %g %g %g <-%s- (%s) %g %g %g\n", me,
359			airEnumStr(tenTripleType, dstType),
360			dst[0], dst[1], dst[2],
361			direct ? "-" : "...",
362			airEnumStr(tenTripleType, srcType),
363			src[0], src[1], src[2]);
364			}
365
366			return;
367			}
368
369			void
370			tenTripleConvertSingle_f(float _dst[3], int dstType,
371			const float _src[3], const int srcType) {
372			double dst[3], src[3];
373
374			ELL_3V_COPY(src, _src);
375			tenTripleConvertSingle_d(dst, dstType, src, srcType);
376			ELL_3V_COPY_TT(_dst, float, dst);
377			}
378
379			void
380			tenTripleCalcSingle_d(double dst[3], int ttype, double ten[7]) {
381			double eval[3];
382
383			/* in time this can become more efficient ... */
384			switch (ttype) {
385			case tenTripleTypeEigenvalue:
386			tenEigensolve_d(dst, NULL, ten);
387			break;
388			case tenTripleTypeMoment:
389			case tenTripleTypeXYZ:
390			case tenTripleTypeRThetaZ:
391			case tenTripleTypeRThetaPhi:
392			case tenTripleTypeK:
393			case tenTripleTypeJ:
394			case tenTripleTypeWheelParm:
395			tenEigensolve_d(eval, NULL, ten);
396			tenTripleConvertSingle_d(dst, ttype, eval, tenTripleTypeEigenvalue);
397			break;
398			case tenTripleTypeR:
399			dst[0] = sqrt(_tenAnisoTen_d[tenAniso_S](ten));
400			dst[1] = _tenAnisoTen_d[tenAniso_FA](ten);
401			dst[2] = _tenAnisoTen_d[tenAniso_Mode](ten);
402			break;
403			default:
404			/* what on earth? */
405			ELL_3V_SET(dst, AIR_NAN, AIR_NAN, AIR_NAN);
406			}
407			return;
408			}
409
410			void
411			tenTripleCalcSingle_f(float dst[3], int ttype, float ten[7]) {
412			double dst_d[3], ten_d[7];
413
414			TEN_T_COPY(ten_d, ten);
415			tenTripleCalcSingle_d(dst_d, ttype, ten_d);
416			ELL_3V_COPY_TT(dst, float, dst_d);
417			return;
418			}
419
420			int
421			tenTripleCalc(Nrrd nout, int ttype, const Nrrd nten) {
422			static const char me[]="tenTripleCalc";
423			size_t II, NN, size[NRRD_DIM_MAX];
424			double (ins)(void , size_t, double), (lup)(const void , size_t);
425
426			if (!( nout && nten )) {
427			biffAddf(TEN, "%s: got NULL pointer", me);
428			return 1;
429			}
430			if (airEnumValCheck(tenTripleType, ttype)) {
431			biffAddf(TEN, "%s: got invalid %s (%d)", me,
432			tenTripleType->name, ttype);
433			return 1;
434			}
435			if (tenTensorCheck(nten, nrrdTypeDefault, AIR_FALSE, AIR_TRUE)) {
436			biffAddf(TEN, "%s: didn't get a valid DT array", me);
437			return 1;
438			}
439			if (!( nrrdTypeFloat == nten->type \|\|
440			nrrdTypeDouble == nten->type )) {
441			biffAddf(TEN, "%s: need input type %s or %s, not %s\n", me,
442			airEnumStr(nrrdType, nrrdTypeFloat),
443			airEnumStr(nrrdType, nrrdTypeFloat),
444			airEnumStr(nrrdType, nten->type));
445			}
446
447			nrrdAxisInfoGet_nva(nten, nrrdAxisInfoSize, size);
448			size[0] = 3;
449			if (nrrdMaybeAlloc_nva(nout, nten->type, nten->dim, size)) {
450			biffMovef(TEN, NRRD, "%s: couldn't alloc output", me);
451			return 1;
452			}
453
454			NN = nrrdElementNumber(nten)/7;
455			lup = nrrdDLookup[nten->type];
456			ins = nrrdDInsert[nten->type];
457			for (II=0; II<NN; II++) {
458			double ten[7], trip[3];
459			unsigned int vv;
460			for (vv=0; vv<7; vv++) {
461			ten[vv] = lup(nten->data, vv + 7*II);
462			}
463			tenTripleCalcSingle_d(trip, ttype, ten);
464			for (vv=0; vv<3; vv++) {
465			ins(nout->data, vv + 3*II, trip[vv]);
466			}
467			}
468			if (nrrdAxisInfoCopy(nout, nten, NULL, (NRRD_AXIS_INFO_SIZE_BIT))) {
469			biffMovef(TEN, NRRD, "%s: couldn't copy axis info", me);
470			return 1;
471			}
472			nout->axis[0].kind = nrrdKindUnknown;
473			if (nrrdBasicInfoCopy(nout, nten,
474			NRRD_BASIC_INFO_ALL ^ NRRD_BASIC_INFO_SPACE)) {
475			biffAddf(TEN, "%s:", me);
476			return 1;
477			}
478
479			return 0;
480			}
481
482			int
483			tenTripleConvert(Nrrd *nout, int dstType,
484			const Nrrd *nin, int srcType) {
485			static const char me[]="tenTripleConvert";
486			size_t II, NN;
487			double (ins)(void , size_t, double), (lup)(const void , size_t);
488
489			if (!( nout && nin )) {
490			biffAddf(TEN, "%s: got NULL pointer", me);
491			return 1;
492			}
493			if ( airEnumValCheck(tenTripleType, dstType) \|\|
494			airEnumValCheck(tenTripleType, srcType) ) {
495			biffAddf(TEN, "%s: got invalid %s dst (%d) or src (%d)", me,
496			tenTripleType->name, dstType, srcType);
497			return 1;
498			}
499			if (3 != nin->axis[0].size) {
500			char stmp[AIR_STRLEN_SMALL];
501			biffAddf(TEN, "%s: need axis[0].size 3, not %s", me,
502			airSprintSize_t(stmp, nin->axis[0].size));
503			return 1;
504			}
505			if (nrrdTypeBlock == nin->type) {
506			biffAddf(TEN, "%s: input has non-scalar %s type",
507			me, airEnumStr(nrrdType, nrrdTypeBlock));
508			return 1;
509			}
510
511			if (nrrdCopy(nout, nin)) {
512			biffMovef(TEN, NRRD, "%s: couldn't initialize output", me);
513			return 1;
514			}
515			lup = nrrdDLookup[nin->type];
516			ins = nrrdDInsert[nout->type];
517			NN = nrrdElementNumber(nin)/3;
518			for (II=0; II<NN; II++) {
519			double src[3], dst[3];
520			src[0] = lup(nin->data, 0 + 3*II);
521			src[1] = lup(nin->data, 1 + 3*II);
522			src[2] = lup(nin->data, 2 + 3*II);
523			tenTripleConvertSingle_d(dst, dstType, src, srcType);
524			ins(nout->data, 0 + 3*II, dst[0]);
525			ins(nout->data, 1 + 3*II, dst[1]);
526			ins(nout->data, 2 + 3*II, dst[2]);
527			}
528
529			return 0;
530			}