Head

GCC Code Coverage Report

Directory:	./		Exec	Total	Coverage
File:	src/limn/splineEval.c	Lines:	0	126	0.0 %
Date:	2017-05-26	Branches:	0	98	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 "limn.h"

void
_limnSplineIntervalFind_Unknown(int *ii, double *ff,
                                limnSpline *spline, double tt) {
  static const char me[]="_limnSplineIntervalFind_Unknown";

  AIR_UNUSED(ii);
  AIR_UNUSED(ff);
  AIR_UNUSED(spline);
  AIR_UNUSED(tt);
  fprintf(stderr, "%s: WARNING: spline type unset somewhere\n", me);
  return;
}

void
_limnSplineIntervalFind_NonWarp(int *ii, double *ff,
                                limnSpline *spline, double tt) {
  int N;

  N = spline->ncpt->axis[2].size + (spline->loop ? 1 : 0);
  tt = AIR_CLAMP(0, tt, N-1);
  *ii = (int)tt;
  *ff = tt - *ii;
  return;
}

void
_limnSplineIntervalFind_Warp(int *ii, double *ff,
                             limnSpline *spline, double tt) {
  int N;

  N = spline->ncpt->axis[2].size;
  tt = AIR_CLAMP(spline->time[0], tt, spline->time[N-1]);
  *ii = AIR_CLAMP(0, *ii, N-2);
  /* the last value of ii may be the right one */
  if (!AIR_IN_CL(spline->time[*ii], tt, spline->time[*ii+1])) {
    /* HEY: make this a binary search */
    for (*ii=0; *ii<N-2; (*ii)++) {
      if (AIR_IN_CL(spline->time[*ii], tt, spline->time[*ii+1])) {
        break;
      }
    }
  }
  *ff = (tt - spline->time[*ii])/(spline->time[*ii+1] - spline->time[*ii]);
  return;
}

typedef void (*_limnSplineIntervalFind_t)(int *, double *,
                                          limnSpline *, double);
_limnSplineIntervalFind_t
_limnSplineIntervalFind[LIMN_SPLINE_TYPE_MAX+1] = {
  _limnSplineIntervalFind_Unknown,
  _limnSplineIntervalFind_NonWarp,
  _limnSplineIntervalFind_Warp,
  _limnSplineIntervalFind_NonWarp,
  _limnSplineIntervalFind_NonWarp,
  _limnSplineIntervalFind_NonWarp
};

void
_limnSplineWeightsFind_Unknown(double *wght, limnSpline *spline, double f) {
  static const char me[]="_limnSplineWeights_Unknown";

  AIR_UNUSED(wght);
  AIR_UNUSED(spline);
  AIR_UNUSED(f);
  fprintf(stderr, "%s: WARNING: spline type unset somewhere\n", me);
  return;
}

void
_limnSplineWeightsFind_Linear(double *wght, limnSpline *spline, double f) {

  AIR_UNUSED(spline);
  ELL_4V_SET(wght, 0, 1-f, f, 0);
  /*
  fprintf(stderr, "%g ----> %g %g %g %g\n", f,
          wght[0], wght[1], wght[2], wght[3]);
  */
  return;
}

void
_limnSplineWeightsFind_Hermite(double *wght, limnSpline *spline, double f) {
  double f3, f2;

  AIR_UNUSED(spline);
  f3 = f*(f2 = f*f);
  ELL_4V_SET(wght,
             2*f3 - 3*f2 + 1,
             f3 - 2*f2 + f,
             f3 - f2,
             -2*f3 + 3*f2);
  return;
}

void
_limnSplineWeightsFind_CubicBezier(double *wght,
                                   limnSpline *spline, double f) {
  double g;

  AIR_UNUSED(spline);
  g = 1 - f;
  ELL_4V_SET(wght,
             g*g*g,
             3*g*g*f,
             3*g*f*f,
             f*f*f);
  return;
}

/* lifted from nrrd/kernel.c */
#define _BCCUBIC(x, B, C)                                           \
  ((x) >= 2.0 ? 0 :                                                 \
  ((x) >= 1.0                                                       \
   ? (((-B/6 - C)*(x) + B + 5*C)*(x) -2*B - 8*C)*(x) + 4*B/3 + 4*C  \
   : ((2 - 3*B/2 - C)*(x) - 3 + 2*B + C)*(x)*(x) + 1 - B/3))

void
_limnSplineWeightsFind_BC(double *wght, limnSpline *spline, double f) {
  double B, C, f0, f1, f2, f3;

  B = spline->B;
  C = spline->C;
  f0 = f+1;
  f1 = f;
  f2 = AIR_ABS(f-1);
  f3 = AIR_ABS(f-2);
  ELL_4V_SET(wght,
             _BCCUBIC(f0, B, C),
             _BCCUBIC(f1, B, C),
             _BCCUBIC(f2, B, C),
             _BCCUBIC(f3, B, C));
  return;
}

typedef void (*_limnSplineWeightsFind_t)(double *, limnSpline *, double);

_limnSplineWeightsFind_t
_limnSplineWeightsFind[LIMN_SPLINE_TYPE_MAX+1] = {
  _limnSplineWeightsFind_Unknown,
  _limnSplineWeightsFind_Linear,
  _limnSplineWeightsFind_Hermite, /* TimeWarp */
  _limnSplineWeightsFind_Hermite,
  _limnSplineWeightsFind_CubicBezier,
  _limnSplineWeightsFind_BC
};

void
_limnSplineIndexFind(int *idx, limnSpline *spline, int ii) {
  int N, ti[4];

  N = spline->ncpt->axis[2].size;
  if (limnSplineTypeHasImplicitTangents[spline->type]) {
    if (spline->loop) {
      ELL_4V_SET(ti,
                 AIR_MOD(ii-1, N),
                 AIR_MOD(ii+0, N),
                 AIR_MOD(ii+1, N),
                 AIR_MOD(ii+2, N));
    } else {
      ELL_4V_SET(ti,
                 AIR_CLAMP(0, ii-1, N-1),
                 AIR_CLAMP(0, ii+0, N-1),
                 AIR_CLAMP(0, ii+1, N-1),
                 AIR_CLAMP(0, ii+2, N-1));
    }
    ELL_4V_SET(idx, 1 + 3*ti[0], 1 + 3*ti[1], 1 + 3*ti[2], 1 + 3*ti[3]);
  } else {
    if (spline->loop) {
      ELL_4V_SET(ti,
                 AIR_MOD(ii+0, N),
                 AIR_MOD(ii+0, N),
                 AIR_MOD(ii+1, N),
                 AIR_MOD(ii+1, N));
    } else {
      ELL_4V_SET(ti,
                 AIR_CLAMP(0, ii+0, N-1),
                 AIR_CLAMP(0, ii+0, N-1),
                 AIR_CLAMP(0, ii+1, N-1),
                 AIR_CLAMP(0, ii+1, N-1));
    }
    ELL_4V_SET(idx, 1 + 3*ti[0], 2 + 3*ti[1], 0 + 3*ti[2], 1 + 3*ti[3]);
  }
}

void
_limnSplineFinish_Unknown(double *out, limnSpline *spline,
                          int ii, double *wght) {
  static const char me[]="_limnSplineFinish_Unknown";

  AIR_UNUSED(out);
  AIR_UNUSED(spline);
  AIR_UNUSED(ii);
  AIR_UNUSED(wght);
  fprintf(stderr, "%s: WARNING: spline info unset somewhere\n", me);
  return;
}

void
_limnSplineFinish_Scalar(double *out, limnSpline *spline,
                         int ii, double *wght) {
  int idx[4];
  double *cpt;

  cpt = (double*)(spline->ncpt->data);
  _limnSplineIndexFind(idx, spline, ii);
  *out = (  wght[0]*cpt[idx[0]] + wght[1]*cpt[idx[1]]
          + wght[2]*cpt[idx[2]] + wght[3]*cpt[idx[3]]);
  return;
}

void
_limnSplineFinish_2Vec(double *out, limnSpline *spline,
                       int ii, double *wght) {
  int idx[4];
  double *cpt;

  cpt = (double*)(spline->ncpt->data);
  _limnSplineIndexFind(idx, spline, ii);
  out[0] = (  wght[0]*cpt[0 + 2*idx[0]] + wght[1]*cpt[0 + 2*idx[1]]
            + wght[2]*cpt[0 + 2*idx[2]] + wght[3]*cpt[0 + 2*idx[3]]);
  out[1] = (  wght[0]*cpt[1 + 2*idx[0]] + wght[1]*cpt[1 + 2*idx[1]]
            + wght[2]*cpt[1 + 2*idx[2]] + wght[3]*cpt[1 + 2*idx[3]]);
  return;
}

void
_limnSplineFinish_3Vec(double *out, limnSpline *spline,
                       int ii, double *wght) {
  int idx[4];
  double *cpt;

  cpt = (double*)(spline->ncpt->data);
  _limnSplineIndexFind(idx, spline, ii);
  out[0] = (  wght[0]*cpt[0 + 3*idx[0]] + wght[1]*cpt[0 + 3*idx[1]]
            + wght[2]*cpt[0 + 3*idx[2]] + wght[3]*cpt[0 + 3*idx[3]]);
  out[1] = (  wght[0]*cpt[1 + 3*idx[0]] + wght[1]*cpt[1 + 3*idx[1]]
            + wght[2]*cpt[1 + 3*idx[2]] + wght[3]*cpt[1 + 3*idx[3]]);
  out[2] = (  wght[0]*cpt[2 + 3*idx[0]] + wght[1]*cpt[2 + 3*idx[1]]
            + wght[2]*cpt[2 + 3*idx[2]] + wght[3]*cpt[2 + 3*idx[3]]);
  return;
}

void
_limnSplineFinish_Normal(double *out, limnSpline *spline,
                         int ii, double *wght) {

  AIR_UNUSED(out);
  AIR_UNUSED(spline);
  AIR_UNUSED(ii);
  AIR_UNUSED(wght);
  fprintf(stderr, "%s: NOT IMPLEMENTED\n", "_limnSplineFinish_Normal");
  return;
}

void
_limnSplineFinish_4Vec(double *out, limnSpline *spline,
                       int ii, double *wght) {
  int idx[4];
  double *cpt;

  cpt = (double*)(spline->ncpt->data);
  _limnSplineIndexFind(idx, spline, ii);
  out[0] = (  wght[0]*cpt[0 + 4*idx[0]] + wght[1]*cpt[0 + 4*idx[1]]
            + wght[2]*cpt[0 + 4*idx[2]] + wght[3]*cpt[0 + 4*idx[3]]);
  out[1] = (  wght[0]*cpt[1 + 4*idx[0]] + wght[1]*cpt[1 + 4*idx[1]]
            + wght[2]*cpt[1 + 4*idx[2]] + wght[3]*cpt[1 + 4*idx[3]]);
  out[2] = (  wght[0]*cpt[2 + 4*idx[0]] + wght[1]*cpt[2 + 4*idx[1]]
            + wght[2]*cpt[2 + 4*idx[2]] + wght[3]*cpt[2 + 4*idx[3]]);
  out[3] = (  wght[0]*cpt[3 + 4*idx[0]] + wght[1]*cpt[3 + 4*idx[1]]
            + wght[2]*cpt[3 + 4*idx[2]] + wght[3]*cpt[3 + 4*idx[3]]);
  return;
}

/*
** HEY: I have no whether Hermite splines work with this
*/
void
_limnSplineFinish_Quaternion(double *out, limnSpline *spline,
                             int ii, double *wght) {
  int idx[4];
  double *cpt;

  cpt = (double*)(spline->ncpt->data);
  _limnSplineIndexFind(idx, spline, ii);
  ell_q_avg4_d(out, NULL,
               cpt + 4*idx[0], cpt + 4*idx[1],
               cpt + 4*idx[2], cpt + 4*idx[3],
               wght, LIMN_SPLINE_Q_AVG_EPS, 30 /* maxIter */);
  return;
}

typedef void (*_limnSplineFinish_t)(double *, limnSpline *, int, double *);
_limnSplineFinish_t
_limnSplineFinish[LIMN_SPLINE_INFO_MAX+1] = {
  _limnSplineFinish_Unknown,
  _limnSplineFinish_Scalar,
  _limnSplineFinish_2Vec,
  _limnSplineFinish_3Vec,
  _limnSplineFinish_Normal,
  _limnSplineFinish_4Vec,
  _limnSplineFinish_Quaternion
};

void
limnSplineEvaluate(double *out, limnSpline *spline, double tt) {
  int ii=0;
  double ff, wght[4];

  if (out && spline) {
    _limnSplineIntervalFind[spline->type](&ii, &ff, spline, tt);
    _limnSplineWeightsFind[spline->type](wght, spline, ff);
    _limnSplineFinish[spline->info](out, spline, ii, wght);
  }
  return;
}

int
limnSplineNrrdEvaluate(Nrrd *nout, limnSpline *spline, Nrrd *nin) {
  static const char me[]="limnSplineNrrdEvaluate";
  double tt, *out, (*lup)(const void *, size_t);
  int odim, infoSize;
  size_t I, M, size[NRRD_DIM_MAX+1];

  if (!(nout && spline && nin)) {
    biffAddf(LIMN, "%s: got NULL pointer", me);
    return 1;
  }
  if (limnSplineInfoScalar == spline->info) {
    nrrdAxisInfoGet_va(nin, nrrdAxisInfoSize, size);
    infoSize = 1;
    odim = nin->dim;
  } else {
    nrrdAxisInfoGet_va(nin, nrrdAxisInfoSize, size+1);
    infoSize = size[0] = limnSplineInfoSize[spline->info];
    odim = 1 + nin->dim;
  }
  if (nrrdMaybeAlloc_nva(nout, nrrdTypeDouble, odim, size)) {
    biffMovef(LIMN, NRRD, "%s: output allocation failed", me);
    return 1;
  }
  lup = nrrdDLookup[nin->type];
  out = (double*)(nout->data);
  M = nrrdElementNumber(nin);
  for (I=0; I<M; I++) {
    tt = lup(nin->data, I);
    limnSplineEvaluate(out, spline, tt);
    out += infoSize;
  }

  /* HEY: peripheral info copying? */

  return 0;
}

int
limnSplineSample(Nrrd *nout, limnSpline *spline,
                 double minT, size_t M, double maxT) {
  static const char me[]="limnSplineSample";
  airArray *mop;
  Nrrd *ntt;
  double *tt;
  size_t I;

  if (!(nout && spline)) {
    biffAddf(LIMN, "%s: got NULL pointer", me);
    return 1;
  }
  mop = airMopNew();
  airMopAdd(mop, ntt=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
  if (nrrdMaybeAlloc_va(ntt, nrrdTypeDouble, 1,
                        M)) {
    biffMovef(LIMN, NRRD, "%s: trouble allocating tmp nrrd", me);
    airMopError(mop); return 1;
  }
  tt = (double*)(ntt->data);
  for (I=0; I<M; I++) {
    tt[I] = AIR_AFFINE(0, I, M-1, minT, maxT);
  }
  if (limnSplineNrrdEvaluate(nout, spline, ntt)) {
    biffAddf(LIMN, "%s: trouble", me);
    airMopError(mop); return 1;
  }
  airMopOkay(mop);
  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 "limn.h"
26
27			void
28			_limnSplineIntervalFind_Unknown(int ii, double ff,
29			limnSpline *spline, double tt) {
30			static const char me[]="_limnSplineIntervalFind_Unknown";
31
32			AIR_UNUSED(ii);
33			AIR_UNUSED(ff);
34			AIR_UNUSED(spline);
35			AIR_UNUSED(tt);
36			fprintf(stderr, "%s: WARNING: spline type unset somewhere\n", me);
37			return;
38			}
39
40			void
41			_limnSplineIntervalFind_NonWarp(int ii, double ff,
42			limnSpline *spline, double tt) {
43			int N;
44
45			N = spline->ncpt->axis[2].size + (spline->loop ? 1 : 0);
46			tt = AIR_CLAMP(0, tt, N-1);
47			*ii = (int)tt;
48			ff = tt - ii;
49			return;
50			}
51
52			void
53			_limnSplineIntervalFind_Warp(int ii, double ff,
54			limnSpline *spline, double tt) {
55			int N;
56
57			N = spline->ncpt->axis[2].size;
58			tt = AIR_CLAMP(spline->time[0], tt, spline->time[N-1]);
59			ii = AIR_CLAMP(0, ii, N-2);
60			/* the last value of ii may be the right one */
61			if (!AIR_IN_CL(spline->time[ii], tt, spline->time[ii+1])) {
62			/* HEY: make this a binary search */
63			for (ii=0; ii<N-2; (*ii)++) {
64			if (AIR_IN_CL(spline->time[ii], tt, spline->time[ii+1])) {
65			break;
66			}
67			}
68			}
69			ff = (tt - spline->time[ii])/(spline->time[ii+1] - spline->time[ii]);
70			return;
71			}
72
73			typedef void (_limnSplineIntervalFind_t)(int , double *,
74			limnSpline *, double);
75			_limnSplineIntervalFind_t
76			_limnSplineIntervalFind[LIMN_SPLINE_TYPE_MAX+1] = {
77			_limnSplineIntervalFind_Unknown,
78			_limnSplineIntervalFind_NonWarp,
79			_limnSplineIntervalFind_Warp,
80			_limnSplineIntervalFind_NonWarp,
81			_limnSplineIntervalFind_NonWarp,
82			_limnSplineIntervalFind_NonWarp
83			};
84
85			void
86			_limnSplineWeightsFind_Unknown(double wght, limnSpline spline, double f) {
87			static const char me[]="_limnSplineWeights_Unknown";
88
89			AIR_UNUSED(wght);
90			AIR_UNUSED(spline);
91			AIR_UNUSED(f);
92			fprintf(stderr, "%s: WARNING: spline type unset somewhere\n", me);
93			return;
94			}
95
96			void
97			_limnSplineWeightsFind_Linear(double wght, limnSpline spline, double f) {
98
99			AIR_UNUSED(spline);
100			ELL_4V_SET(wght, 0, 1-f, f, 0);
101			/*
102			fprintf(stderr, "%g ----> %g %g %g %g\n", f,
103			wght[0], wght[1], wght[2], wght[3]);
104			*/
105			return;
106			}
107
108			void
109			_limnSplineWeightsFind_Hermite(double wght, limnSpline spline, double f) {
110			double f3, f2;
111
112			AIR_UNUSED(spline);
113			f3 = f(f2 = ff);
114			ELL_4V_SET(wght,
115			2f3 - 3f2 + 1,
116			f3 - 2*f2 + f,
117			f3 - f2,
118			-2f3 + 3f2);
119			return;
120			}
121
122			void
123			_limnSplineWeightsFind_CubicBezier(double *wght,
124			limnSpline *spline, double f) {
125			double g;
126
127			AIR_UNUSED(spline);
128			g = 1 - f;
129			ELL_4V_SET(wght,
130			ggg,
131			3gg*f,
132			3gf*f,
133			fff);
134			return;
135			}
136
137			/* lifted from nrrd/kernel.c */
138			#define _BCCUBIC(x, B, C) \
139			((x) >= 2.0 ? 0 : \
140			((x) >= 1.0 \
141			? (((-B/6 - C)(x) + B + 5C)(x) -2B - 8C)(x) + 4B/3 + 4C \
142			: ((2 - 3B/2 - C)(x) - 3 + 2B + C)(x)*(x) + 1 - B/3))
143
144			void
145			_limnSplineWeightsFind_BC(double wght, limnSpline spline, double f) {
146			double B, C, f0, f1, f2, f3;
147
148			B = spline->B;
149			C = spline->C;
150			f0 = f+1;
151			f1 = f;
152			f2 = AIR_ABS(f-1);
153			f3 = AIR_ABS(f-2);
154			ELL_4V_SET(wght,
155			_BCCUBIC(f0, B, C),
156			_BCCUBIC(f1, B, C),
157			_BCCUBIC(f2, B, C),
158			_BCCUBIC(f3, B, C));
159			return;
160			}
161
162			typedef void (_limnSplineWeightsFind_t)(double , limnSpline *, double);
163
164			_limnSplineWeightsFind_t
165			_limnSplineWeightsFind[LIMN_SPLINE_TYPE_MAX+1] = {
166			_limnSplineWeightsFind_Unknown,
167			_limnSplineWeightsFind_Linear,
168			_limnSplineWeightsFind_Hermite, /* TimeWarp */
169			_limnSplineWeightsFind_Hermite,
170			_limnSplineWeightsFind_CubicBezier,
171			_limnSplineWeightsFind_BC
172			};
173
174			void
175			_limnSplineIndexFind(int idx, limnSpline spline, int ii) {
176			int N, ti[4];
177
178			N = spline->ncpt->axis[2].size;
179			if (limnSplineTypeHasImplicitTangents[spline->type]) {
180			if (spline->loop) {
181			ELL_4V_SET(ti,
182			AIR_MOD(ii-1, N),
183			AIR_MOD(ii+0, N),
184			AIR_MOD(ii+1, N),
185			AIR_MOD(ii+2, N));
186			} else {
187			ELL_4V_SET(ti,
188			AIR_CLAMP(0, ii-1, N-1),
189			AIR_CLAMP(0, ii+0, N-1),
190			AIR_CLAMP(0, ii+1, N-1),
191			AIR_CLAMP(0, ii+2, N-1));
192			}
193			ELL_4V_SET(idx, 1 + 3ti[0], 1 + 3ti[1], 1 + 3ti[2], 1 + 3ti[3]);
194			} else {
195			if (spline->loop) {
196			ELL_4V_SET(ti,
197			AIR_MOD(ii+0, N),
198			AIR_MOD(ii+0, N),
199			AIR_MOD(ii+1, N),
200			AIR_MOD(ii+1, N));
201			} else {
202			ELL_4V_SET(ti,
203			AIR_CLAMP(0, ii+0, N-1),
204			AIR_CLAMP(0, ii+0, N-1),
205			AIR_CLAMP(0, ii+1, N-1),
206			AIR_CLAMP(0, ii+1, N-1));
207			}
208			ELL_4V_SET(idx, 1 + 3ti[0], 2 + 3ti[1], 0 + 3ti[2], 1 + 3ti[3]);
209			}
210			}
211
212			void
213			_limnSplineFinish_Unknown(double out, limnSpline spline,
214			int ii, double *wght) {
215			static const char me[]="_limnSplineFinish_Unknown";
216
217			AIR_UNUSED(out);
218			AIR_UNUSED(spline);
219			AIR_UNUSED(ii);
220			AIR_UNUSED(wght);
221			fprintf(stderr, "%s: WARNING: spline info unset somewhere\n", me);
222			return;
223			}
224
225			void
226			_limnSplineFinish_Scalar(double out, limnSpline spline,
227			int ii, double *wght) {
228			int idx[4];
229			double *cpt;
230
231			cpt = (double*)(spline->ncpt->data);
232			_limnSplineIndexFind(idx, spline, ii);
233			out = ( wght[0]cpt[idx[0]] + wght[1]*cpt[idx[1]]
234			+ wght[2]cpt[idx[2]] + wght[3]cpt[idx[3]]);
235			return;
236			}
237
238			void
239			_limnSplineFinish_2Vec(double out, limnSpline spline,
240			int ii, double *wght) {
241			int idx[4];
242			double *cpt;
243
244			cpt = (double*)(spline->ncpt->data);
245			_limnSplineIndexFind(idx, spline, ii);
246			out[0] = ( wght[0]cpt[0 + 2idx[0]] + wght[1]cpt[0 + 2idx[1]]
247			+ wght[2]cpt[0 + 2idx[2]] + wght[3]cpt[0 + 2idx[3]]);
248			out[1] = ( wght[0]cpt[1 + 2idx[0]] + wght[1]cpt[1 + 2idx[1]]
249			+ wght[2]cpt[1 + 2idx[2]] + wght[3]cpt[1 + 2idx[3]]);
250			return;
251			}
252
253			void
254			_limnSplineFinish_3Vec(double out, limnSpline spline,
255			int ii, double *wght) {
256			int idx[4];
257			double *cpt;
258
259			cpt = (double*)(spline->ncpt->data);
260			_limnSplineIndexFind(idx, spline, ii);
261			out[0] = ( wght[0]cpt[0 + 3idx[0]] + wght[1]cpt[0 + 3idx[1]]
262			+ wght[2]cpt[0 + 3idx[2]] + wght[3]cpt[0 + 3idx[3]]);
263			out[1] = ( wght[0]cpt[1 + 3idx[0]] + wght[1]cpt[1 + 3idx[1]]
264			+ wght[2]cpt[1 + 3idx[2]] + wght[3]cpt[1 + 3idx[3]]);
265			out[2] = ( wght[0]cpt[2 + 3idx[0]] + wght[1]cpt[2 + 3idx[1]]
266			+ wght[2]cpt[2 + 3idx[2]] + wght[3]cpt[2 + 3idx[3]]);
267			return;
268			}
269
270			void
271			_limnSplineFinish_Normal(double out, limnSpline spline,
272			int ii, double *wght) {
273
274			AIR_UNUSED(out);
275			AIR_UNUSED(spline);
276			AIR_UNUSED(ii);
277			AIR_UNUSED(wght);
278			fprintf(stderr, "%s: NOT IMPLEMENTED\n", "_limnSplineFinish_Normal");
279			return;
280			}
281
282			void
283			_limnSplineFinish_4Vec(double out, limnSpline spline,
284			int ii, double *wght) {
285			int idx[4];
286			double *cpt;
287
288			cpt = (double*)(spline->ncpt->data);
289			_limnSplineIndexFind(idx, spline, ii);
290			out[0] = ( wght[0]cpt[0 + 4idx[0]] + wght[1]cpt[0 + 4idx[1]]
291			+ wght[2]cpt[0 + 4idx[2]] + wght[3]cpt[0 + 4idx[3]]);
292			out[1] = ( wght[0]cpt[1 + 4idx[0]] + wght[1]cpt[1 + 4idx[1]]
293			+ wght[2]cpt[1 + 4idx[2]] + wght[3]cpt[1 + 4idx[3]]);
294			out[2] = ( wght[0]cpt[2 + 4idx[0]] + wght[1]cpt[2 + 4idx[1]]
295			+ wght[2]cpt[2 + 4idx[2]] + wght[3]cpt[2 + 4idx[3]]);
296			out[3] = ( wght[0]cpt[3 + 4idx[0]] + wght[1]cpt[3 + 4idx[1]]
297			+ wght[2]cpt[3 + 4idx[2]] + wght[3]cpt[3 + 4idx[3]]);
298			return;
299			}
300
301			/*
302			** HEY: I have no whether Hermite splines work with this
303			*/
304			void
305			_limnSplineFinish_Quaternion(double out, limnSpline spline,
306			int ii, double *wght) {
307			int idx[4];
308			double *cpt;
309
310			cpt = (double*)(spline->ncpt->data);
311			_limnSplineIndexFind(idx, spline, ii);
312			ell_q_avg4_d(out, NULL,
313			cpt + 4idx[0], cpt + 4idx[1],
314			cpt + 4idx[2], cpt + 4idx[3],
315			wght, LIMN_SPLINE_Q_AVG_EPS, 30 /* maxIter */);
316			return;
317			}
318
319			typedef void (_limnSplineFinish_t)(double , limnSpline , int, double );
320			_limnSplineFinish_t
321			_limnSplineFinish[LIMN_SPLINE_INFO_MAX+1] = {
322			_limnSplineFinish_Unknown,
323			_limnSplineFinish_Scalar,
324			_limnSplineFinish_2Vec,
325			_limnSplineFinish_3Vec,
326			_limnSplineFinish_Normal,
327			_limnSplineFinish_4Vec,
328			_limnSplineFinish_Quaternion
329			};
330
331			void
332			limnSplineEvaluate(double out, limnSpline spline, double tt) {
333			int ii=0;
334			double ff, wght[4];
335
336			if (out && spline) {
337			_limnSplineIntervalFind[spline->type](&ii, &ff, spline, tt);
338			_limnSplineWeightsFind[spline->type](wght, spline, ff);
339			_limnSplineFinish[spline->info](out, spline, ii, wght);
340			}
341			return;
342			}
343
344			int
345			limnSplineNrrdEvaluate(Nrrd nout, limnSpline spline, Nrrd *nin) {
346			static const char me[]="limnSplineNrrdEvaluate";
347			double tt, out, (lup)(const void *, size_t);
348			int odim, infoSize;
349			size_t I, M, size[NRRD_DIM_MAX+1];
350
351			if (!(nout && spline && nin)) {
352			biffAddf(LIMN, "%s: got NULL pointer", me);
353			return 1;
354			}
355			if (limnSplineInfoScalar == spline->info) {
356			nrrdAxisInfoGet_va(nin, nrrdAxisInfoSize, size);
357			infoSize = 1;
358			odim = nin->dim;
359			} else {
360			nrrdAxisInfoGet_va(nin, nrrdAxisInfoSize, size+1);
361			infoSize = size[0] = limnSplineInfoSize[spline->info];
362			odim = 1 + nin->dim;
363			}
364			if (nrrdMaybeAlloc_nva(nout, nrrdTypeDouble, odim, size)) {
365			biffMovef(LIMN, NRRD, "%s: output allocation failed", me);
366			return 1;
367			}
368			lup = nrrdDLookup[nin->type];
369			out = (double*)(nout->data);
370			M = nrrdElementNumber(nin);
371			for (I=0; I<M; I++) {
372			tt = lup(nin->data, I);
373			limnSplineEvaluate(out, spline, tt);
374			out += infoSize;
375			}
376
377			/* HEY: peripheral info copying? */
378
379			return 0;
380			}
381
382			int
383			limnSplineSample(Nrrd nout, limnSpline spline,
384			double minT, size_t M, double maxT) {
385			static const char me[]="limnSplineSample";
386			airArray *mop;
387			Nrrd *ntt;
388			double *tt;
389			size_t I;
390
391			if (!(nout && spline)) {
392			biffAddf(LIMN, "%s: got NULL pointer", me);
393			return 1;
394			}
395			mop = airMopNew();
396			airMopAdd(mop, ntt=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
397			if (nrrdMaybeAlloc_va(ntt, nrrdTypeDouble, 1,
398			M)) {
399			biffMovef(LIMN, NRRD, "%s: trouble allocating tmp nrrd", me);
400			airMopError(mop); return 1;
401			}
402			tt = (double*)(ntt->data);
403			for (I=0; I<M; I++) {
404			tt[I] = AIR_AFFINE(0, I, M-1, minT, maxT);
405			}
406			if (limnSplineNrrdEvaluate(nout, spline, ntt)) {
407			biffAddf(LIMN, "%s: trouble", me);
408			airMopError(mop); return 1;
409			}
410			airMopOkay(mop);
411			return 0;
412			}
413