1 /*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 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
_limnSplineIntervalFind_Unknown(int * ii,double * ff,limnSpline * spline,double tt)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
_limnSplineIntervalFind_NonWarp(int * ii,double * ff,limnSpline * spline,double tt)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
_limnSplineIntervalFind_Warp(int * ii,double * ff,limnSpline * spline,double tt)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
_limnSplineWeightsFind_Unknown(double * wght,limnSpline * spline,double f)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
_limnSplineWeightsFind_Linear(double * wght,limnSpline * spline,double f)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
_limnSplineWeightsFind_Hermite(double * wght,limnSpline * spline,double f)109 _limnSplineWeightsFind_Hermite(double *wght, limnSpline *spline, double f) {
110 double f3, f2;
111
112 AIR_UNUSED(spline);
113 f3 = f*(f2 = f*f);
114 ELL_4V_SET(wght,
115 2*f3 - 3*f2 + 1,
116 f3 - 2*f2 + f,
117 f3 - f2,
118 -2*f3 + 3*f2);
119 return;
120 }
121
122 void
_limnSplineWeightsFind_CubicBezier(double * wght,limnSpline * spline,double f)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 g*g*g,
131 3*g*g*f,
132 3*g*f*f,
133 f*f*f);
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 + 5*C)*(x) -2*B - 8*C)*(x) + 4*B/3 + 4*C \
142 : ((2 - 3*B/2 - C)*(x) - 3 + 2*B + C)*(x)*(x) + 1 - B/3))
143
144 void
_limnSplineWeightsFind_BC(double * wght,limnSpline * spline,double f)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
_limnSplineIndexFind(int * idx,limnSpline * spline,int ii)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 + 3*ti[0], 1 + 3*ti[1], 1 + 3*ti[2], 1 + 3*ti[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 + 3*ti[0], 2 + 3*ti[1], 0 + 3*ti[2], 1 + 3*ti[3]);
209 }
210 }
211
212 void
_limnSplineFinish_Unknown(double * out,limnSpline * spline,int ii,double * wght)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
_limnSplineFinish_Scalar(double * out,limnSpline * spline,int ii,double * wght)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
_limnSplineFinish_2Vec(double * out,limnSpline * spline,int ii,double * wght)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 + 2*idx[0]] + wght[1]*cpt[0 + 2*idx[1]]
247 + wght[2]*cpt[0 + 2*idx[2]] + wght[3]*cpt[0 + 2*idx[3]]);
248 out[1] = ( wght[0]*cpt[1 + 2*idx[0]] + wght[1]*cpt[1 + 2*idx[1]]
249 + wght[2]*cpt[1 + 2*idx[2]] + wght[3]*cpt[1 + 2*idx[3]]);
250 return;
251 }
252
253 void
_limnSplineFinish_3Vec(double * out,limnSpline * spline,int ii,double * wght)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 + 3*idx[0]] + wght[1]*cpt[0 + 3*idx[1]]
262 + wght[2]*cpt[0 + 3*idx[2]] + wght[3]*cpt[0 + 3*idx[3]]);
263 out[1] = ( wght[0]*cpt[1 + 3*idx[0]] + wght[1]*cpt[1 + 3*idx[1]]
264 + wght[2]*cpt[1 + 3*idx[2]] + wght[3]*cpt[1 + 3*idx[3]]);
265 out[2] = ( wght[0]*cpt[2 + 3*idx[0]] + wght[1]*cpt[2 + 3*idx[1]]
266 + wght[2]*cpt[2 + 3*idx[2]] + wght[3]*cpt[2 + 3*idx[3]]);
267 return;
268 }
269
270 void
_limnSplineFinish_Normal(double * out,limnSpline * spline,int ii,double * wght)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
_limnSplineFinish_4Vec(double * out,limnSpline * spline,int ii,double * wght)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 + 4*idx[0]] + wght[1]*cpt[0 + 4*idx[1]]
291 + wght[2]*cpt[0 + 4*idx[2]] + wght[3]*cpt[0 + 4*idx[3]]);
292 out[1] = ( wght[0]*cpt[1 + 4*idx[0]] + wght[1]*cpt[1 + 4*idx[1]]
293 + wght[2]*cpt[1 + 4*idx[2]] + wght[3]*cpt[1 + 4*idx[3]]);
294 out[2] = ( wght[0]*cpt[2 + 4*idx[0]] + wght[1]*cpt[2 + 4*idx[1]]
295 + wght[2]*cpt[2 + 4*idx[2]] + wght[3]*cpt[2 + 4*idx[3]]);
296 out[3] = ( wght[0]*cpt[3 + 4*idx[0]] + wght[1]*cpt[3 + 4*idx[1]]
297 + wght[2]*cpt[3 + 4*idx[2]] + wght[3]*cpt[3 + 4*idx[3]]);
298 return;
299 }
300
301 /*
302 ** HEY: I have no whether Hermite splines work with this
303 */
304 void
_limnSplineFinish_Quaternion(double * out,limnSpline * spline,int ii,double * wght)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 + 4*idx[0], cpt + 4*idx[1],
314 cpt + 4*idx[2], cpt + 4*idx[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
limnSplineEvaluate(double * out,limnSpline * spline,double tt)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
limnSplineNrrdEvaluate(Nrrd * nout,limnSpline * spline,Nrrd * nin)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
limnSplineSample(Nrrd * nout,limnSpline * spline,double minT,size_t M,double maxT)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
414