GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/limn/splineMethods.c Lines: 0 236 0.0 %
Date: 2017-05-26 Branches: 0 138 0.0 %

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
limnSplineTypeSpec *
28
limnSplineTypeSpecNew(int type, ...) {
29
  static const char me[]="limnSplineTypeSpecNew";
30
  limnSplineTypeSpec *spec;
31
  va_list ap;
32
33
  if (airEnumValCheck(limnSplineType, type)) {
34
    biffAddf(LIMN, "%s: given type %d not a valid limnSplineType", me, type);
35
    return NULL;
36
  }
37
  spec = AIR_CALLOC(1, limnSplineTypeSpec);
38
  spec->type = type;
39
  va_start(ap, type);
40
  if (limnSplineTypeBC == type) {
41
    spec->B = va_arg(ap, double);
42
    spec->C = va_arg(ap, double);
43
  }
44
  va_end(ap);
45
  return spec;
46
}
47
48
limnSplineTypeSpec *
49
limnSplineTypeSpecNix(limnSplineTypeSpec *spec) {
50
51
  airFree(spec);
52
  return NULL;
53
}
54
55
/*
56
** _limnSplineTimeWarpSet
57
**
58
** implements GK's patented time warping technology
59
*/
60
int
61
_limnSplineTimeWarpSet(limnSpline *spline) {
62
  static const char me[]="_limnSplineTimeWarpSet";
63
  double *cpt, *time, ss;
64
  int ii, N;
65
66
  cpt = (double*)(spline->ncpt->data);
67
  N = spline->ncpt->axis[2].size;
68
  time = spline->time;
69
70
  for (ii=0; ii<N; ii++) {
71
    if (!AIR_EXISTS(time[ii])) {
72
      biffAddf(LIMN, "%s: time[%d] doesn't exist", me, ii);
73
      return 1;
74
    }
75
    if (ii && !(time[ii-1] < time[ii])) {
76
      biffAddf(LIMN, "%s: time[%d] = %g not < time[%d] = %g", me,
77
               ii-1, time[ii-1], ii, time[ii]);
78
      return 1;
79
    }
80
    /* this will be used below */
81
    /* HEY: is there any way or reason to do any other kind of time warp? */
82
    cpt[1 + 3*ii] = ii;
83
  }
84
85
  for (ii=1; ii<N-1; ii++) {
86
    ss = (cpt[1+3*(ii+1)] - cpt[1+3*(ii-1)])/(time[ii+1] - time[ii-1]);
87
    cpt[0 + 3*ii] = ss*(time[ii] - time[ii-1]);
88
    cpt[2 + 3*ii] = ss*(time[ii+1] - time[ii]);
89
  }
90
  if (spline->loop) {
91
    ss = ((cpt[1+3*1] - cpt[1+3*0] + cpt[1+3*(N-1)] - cpt[1+3*(N-2)])
92
          / (time[1] - time[0] + time[N-1] - time[N-2]));
93
    cpt[2 + 3*0] = ss*(time[1] - time[0]);
94
    cpt[0 + 3*(N-1)] = ss*(time[N-1] - time[N-2]);
95
  } else {
96
    cpt[2 + 3*0] = ((cpt[1+3*1] - cpt[1+3*0])
97
                    * (time[1] - time[0]));
98
    cpt[0 + 3*(N-1)] = ((cpt[1+3*(N-1)] - cpt[1+3*(N-2)])
99
                        * (time[N-1] - time[N-2]));
100
  }
101
  /*
102
  fprintf(stderr, "s[0]=%g, post = %g; s[1]=%g pre = %g\n",
103
          cpt[1 + 3*0], cpt[2 + 3*0], cpt[1 + 3*1], cpt[0 + 3*1]);
104
  */
105
  return 0;
106
}
107
108
/*
109
******** limnSplineNew
110
**
111
** constructor for limnSplines.  We take all the control point information
112
** here, and copy it internally, in an effort to simplify the management of
113
** state.  The control point nrrd is 3-D, as explained in limn.h
114
**
115
** To confuse matters, the Time type of spline doesn't need the control
116
** point information in the traditional sense, but it still needs to know
117
** "when" the control points are.  So, the _ncpt nrrd is still needed, but
118
** it is only a 1-D array of times.  In this case, the internal ncpt values
119
** are set automatically.
120
**
121
** The benefit of this approach is that if this constructor returns
122
** successfully, then there is no more information or state that needs to
123
** be set-- the returned spline can be passed to evaluate or sample.
124
** LIES LIES LIES: For BC-splines, you still have to call limnSplineBCSet,
125
** but that's the only exception...
126
*/
127
limnSpline *
128
limnSplineNew(Nrrd *_ncpt, int info, limnSplineTypeSpec *spec) {
129
  static const char me[]="limnSplineNew";
130
  limnSpline *spline;
131
  size_t N;
132
  unsigned int size;
133
  airArray *mop;
134
  Nrrd *nin;
135
  char stmp[2][AIR_STRLEN_SMALL];
136
137
  if (airEnumValCheck(limnSplineInfo, info)) {
138
    biffAddf(LIMN, "%s: info %d not a valid limnSplineInfo", me, info);
139
    return NULL;
140
  }
141
  if (nrrdCheck(_ncpt)) {
142
    biffMovef(LIMN, NRRD, "%s: given nrrd has problems", me);
143
    return NULL;
144
  }
145
  if (limnSplineTypeTimeWarp == spec->type) {
146
    if (!(limnSplineInfoScalar == info)) {
147
      biffAddf(LIMN, "%s: can only time warp scalars", me);
148
      return NULL;
149
    }
150
    if (!( 1 == _ncpt->dim )) {
151
      biffAddf(LIMN, "%s: given nrrd has dimension %d, not 1", me, _ncpt->dim);
152
      return NULL;
153
    }
154
    N = _ncpt->axis[0].size;
155
  } else {
156
    if (!( 3 == _ncpt->dim )) {
157
      biffAddf(LIMN, "%s: given nrrd has dimension %d, not 3", me, _ncpt->dim);
158
      return NULL;
159
    }
160
    size = limnSplineInfoSize[info];
161
    if (!( size == _ncpt->axis[0].size && 3 == _ncpt->axis[1].size )) {
162
      biffAddf(LIMN, "%s: expected %ux3xN nrrd, not %sx%sxN", me, size,
163
               airSprintSize_t(stmp[0], _ncpt->axis[0].size),
164
               airSprintSize_t(stmp[1], _ncpt->axis[1].size));
165
      return NULL;
166
    }
167
    N = _ncpt->axis[2].size;
168
  }
169
  if (1 == N) {
170
    biffAddf(LIMN, "%s: need at least two control points", me);
171
    return NULL;
172
  }
173
174
  mop = airMopNew();
175
  if (!( spline = AIR_CALLOC(1, limnSpline) )) {
176
    biffAddf(LIMN, "%s: couldn't allocate new spline", me);
177
    airMopError(mop); return NULL;
178
  }
179
180
  airMopAdd(mop, spline, (airMopper)limnSplineNix, airMopOnError);
181
  spline->time = NULL;
182
  spline->ncpt = NULL;
183
  spline->type = spec->type;
184
  spline->info = info;
185
  spline->loop = AIR_FALSE;
186
  spline->B = spec->B;
187
  spline->C = spec->C;
188
  nin = nrrdNew();
189
  airMopAdd(mop, nin, (airMopper)nrrdNuke, airMopOnError);
190
  if (nrrdConvert(nin, _ncpt, nrrdTypeDouble)) {
191
    biffMovef(LIMN, NRRD, "%s: trouble allocating internal nrrd", me);
192
    airMopError(mop); return NULL;
193
  }
194
  if (limnSplineTypeTimeWarp == spec->type) {
195
    /* set the time array to the data of the double-converted nin,
196
       but the nin itself is scheduled to be nixed */
197
    airMopAdd(mop, nin, (airMopper)nrrdNix, airMopOnOkay);
198
    spline->time = (double*)(nin->data);
199
    /* now allocate the real control point information */
200
    spline->ncpt = nrrdNew();
201
    airMopAdd(mop, spline->ncpt, (airMopper)nrrdNuke, airMopOnError);
202
    if (nrrdMaybeAlloc_va(spline->ncpt, nrrdTypeDouble, 3,
203
                          AIR_CAST(size_t, 1),
204
                          AIR_CAST(size_t, 3),
205
                          _ncpt->axis[0].size)) {
206
      biffMovef(LIMN, NRRD,
207
                "%s: trouble allocating real control points", me);
208
      airMopError(mop); return NULL;
209
    }
210
    /* and set it all to something useful */
211
    if (_limnSplineTimeWarpSet(spline)) {
212
      biffAddf(LIMN, "%s: trouble setting time warp", me);
213
      airMopError(mop); return NULL;
214
    }
215
  } else {
216
    /* we set the control point to the double-converted nin, and we're done */
217
    spline->ncpt = nin;
218
  }
219
  airMopOkay(mop);
220
  return spline;
221
}
222
223
limnSpline *
224
limnSplineNix(limnSpline *spline) {
225
226
  if (spline) {
227
    spline->ncpt = (Nrrd *)nrrdNuke(spline->ncpt);
228
    spline->time = (double *)airFree(spline->time);
229
    airFree(spline);
230
  }
231
  return NULL;
232
}
233
234
/*
235
******** limnSplineNrrdCleverFix
236
**
237
** given that the ncpt nrrd for limnSplineNew() has to be a particular
238
** dimension and shape, and given that convenient ways of creating nrrds
239
** don't always lead to such configurations, we supply some minimal
240
** cleverness to bridge the gap.  As the name implies, you should be
241
** wary of this function.
242
**
243
** The job of this function is NOT to check the validity of a nrrd for
244
** a given spline.  That is up to limnSplineNew().
245
**
246
** Currently, this function is used by limnHestSpline, but it probably
247
** won't be used anywhere else within limn.
248
**
249
** If requested, we also take a stab at guessing limnSplineInfo.
250
*/
251
int
252
limnSplineNrrdCleverFix(Nrrd *nout, Nrrd *nin, int info, int type) {
253
  static const char me[]="limnSplineNrrdCleverFix";
254
  ptrdiff_t min[3], max[3];
255
  size_t N;
256
  unsigned int wantSize;
257
  Nrrd *ntmpA, *ntmpB;
258
  airArray *mop;
259
  char stmp[AIR_STRLEN_SMALL];
260
261
  if (!(nout && nin)) {
262
    biffAddf(LIMN, "%s: got NULL pointer", me);
263
    return 1;
264
  }
265
  if (airEnumValCheck(limnSplineInfo, info)
266
      || airEnumValCheck(limnSplineType, type)) {
267
    biffAddf(LIMN, "%s: invalid spline info (%d) or type (%d)",
268
             me, info, type);
269
    return 1;
270
  }
271
  if (nrrdCheck(nin)) {
272
    biffMovef(LIMN, NRRD, "%s: nrrd has problems", me);
273
    return 1;
274
  }
275
276
  mop = airMopNew();
277
  airMopAdd(mop, ntmpA=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
278
  airMopAdd(mop, ntmpB=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
279
  wantSize = limnSplineInfoSize[info];
280
  switch(nin->dim) {
281
  case 3:
282
    /* we assume that things are okay */
283
    if (nrrdCopy(nout, nin)) {
284
      biffMovef(LIMN, NRRD, "%s: trouble setting output", me);
285
      airMopError(mop); return 1;
286
    }
287
    break;
288
  case 2:
289
    N = nin->axis[1].size;
290
    if (wantSize != nin->axis[0].size) {
291
      biffAddf(LIMN, "%s: expected axis[0].size %d for info %s, but got %s",
292
               me, wantSize, airEnumStr(limnSplineInfo, info),
293
               airSprintSize_t(stmp, nin->axis[0].size));
294
      airMopError(mop); return 1;
295
    }
296
    if (limnSplineTypeTimeWarp == type) {
297
      /* time-warp handled differently */
298
      if (nrrdAxesDelete(nout, nin, 0)) {
299
        biffMovef(LIMN, NRRD, "%s: couldn't make data 1-D", me);
300
        airMopError(mop); return 1;
301
      }
302
    } else {
303
      if (limnSplineTypeHasImplicitTangents[type]) {
304
        ELL_3V_SET(min, 0, -1, 0);
305
        ELL_3V_SET(max, wantSize-1, 1, N-1);
306
        if (nrrdAxesInsert(ntmpA, nin, 1)
307
            || nrrdPad_va(nout, ntmpA, min, max, nrrdBoundaryPad, 0.0)) {
308
          biffMovef(LIMN, NRRD, "%s: trouble with axinsert/pad", me);
309
          airMopError(mop); return 1;
310
        }
311
      } else {
312
        /* the post- and pre- point information may be interlaced with the
313
           main control point values */
314
        if (1 != AIR_MOD((int)N, 3)) {
315
          biffAddf(LIMN,
316
                   "%s: axis[1].size must be 1+(multiple of 3) when using "
317
                   "interlaced tangent information, not %s", me,
318
                   airSprintSize_t(stmp, N));
319
          airMopError(mop); return 1;
320
        }
321
        ELL_2V_SET(min, 0, -1);
322
        ELL_2V_SET(max, wantSize-1, N);
323
        if (nrrdPad_va(ntmpA, nin, min, max, nrrdBoundaryPad, 0.0)
324
            || nrrdAxesSplit(nout, ntmpA, 1, 3, (N+2)/3)) {
325
          biffMovef(LIMN, NRRD, "%s: trouble with pad/axsplit", me);
326
          airMopError(mop); return 1;
327
        }
328
      }
329
    }
330
    break;
331
  case 1:
332
    N = nin->axis[0].size;
333
    if (limnSplineInfoScalar != info) {
334
      biffAddf(LIMN, "%s: can't have %s spline with 1-D nrrd", me,
335
               airEnumStr(limnSplineInfo, info));
336
      airMopError(mop); return 1;
337
    }
338
    if (limnSplineTypeTimeWarp == type) {
339
      /* nothing fancey needed for time-warp */
340
      if (nrrdCopy(nout, nin)) {
341
        biffMovef(LIMN, NRRD, "%s: trouble setting output", me);
342
        airMopError(mop); return 1;
343
      }
344
    } else {
345
      if (limnSplineTypeHasImplicitTangents[type]) {
346
        ELL_3V_SET(min, 0, -1, 0);
347
        ELL_3V_SET(max, 0, 1, N-1);
348
        if (nrrdAxesInsert(ntmpA, nin, 0)
349
            || nrrdAxesInsert(ntmpB, ntmpA, 0)
350
            || nrrdPad_va(nout, ntmpB, min, max, nrrdBoundaryPad, 0.0)) {
351
          biffMovef(LIMN, NRRD, "%s: trouble with axinsert/axinsert/pad", me);
352
          airMopError(mop); return 1;
353
        }
354
      } else {
355
        /* the post- and pre- point information may be interlaced with the
356
           main control point values */
357
        if (1 != AIR_MOD((int)N, 3)) {
358
          biffAddf(LIMN,
359
                   "%s: axis[1].size must be 1+(multiple of 3) when using "
360
                   "interlaced tangent information, not %s", me,
361
                   airSprintSize_t(stmp, N));
362
          airMopError(mop); return 1;
363
        }
364
        ELL_2V_SET(min, 0, -1);
365
        ELL_2V_SET(max, 0, N+1);
366
        if (nrrdAxesInsert(ntmpA, nin, 0)
367
            || nrrdPad_va(ntmpB, ntmpA, min, max, nrrdBoundaryPad, 0.0)
368
            || nrrdAxesSplit(nout, ntmpB, 1, 3, (N+2)/3)) {
369
          biffMovef(LIMN, NRRD, "%s: trouble with axinsert/pad/axsplit", me);
370
          airMopError(mop); return 1;
371
        }
372
      }
373
    }
374
    break;
375
  default:
376
    biffAddf(LIMN, "%s: input nrrd dim %d baffling", me, nin->dim);
377
    return 1;
378
  }
379
  if (nrrdCheck(nout)) {
380
    biffMovef(LIMN, NRRD, "%s: oops: didn't create valid output", me);
381
    airMopError(mop); return 1;
382
  }
383
  airMopOkay(mop);
384
  return 0;
385
}
386
387
limnSpline *
388
limnSplineCleverNew(Nrrd *ncpt, int info, limnSplineTypeSpec *spec) {
389
  static const char me[]="limnSplineCleverNew";
390
  limnSpline *spline;
391
  Nrrd *ntmp;
392
  airArray *mop;
393
394
  if (!( ncpt && spec )) {
395
    biffAddf(LIMN, "%s: got NULL pointer", me);
396
    return NULL;
397
  }
398
  mop = airMopNew();
399
  airMopAdd(mop, ntmp = nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
400
  if (limnSplineNrrdCleverFix(ntmp, ncpt, info, spec->type)) {
401
    biffAddf(LIMN, "%s: couldn't fix up given control point nrrd", me);
402
    airMopError(mop); return NULL;
403
  }
404
  if (!( spline = limnSplineNew(ntmp, info, spec) )) {
405
    biffAddf(LIMN, "%s: couldn't create spline", me);
406
    airMopError(mop); return NULL;
407
  }
408
409
  airMopOkay(mop);
410
  return spline;
411
}
412
413
int
414
limnSplineUpdate(limnSpline *spline, Nrrd *_ncpt) {
415
  static const char me[]="limnSplineUpdate";
416
  Nrrd *ntmp;
417
  char stmp[2][AIR_STRLEN_SMALL];
418
419
  if (!(spline && _ncpt)) {
420
    biffAddf(LIMN, "%s: got NULL pointer", me);
421
    return 1;
422
  }
423
  if (nrrdCheck(_ncpt)) {
424
    biffMovef(LIMN, NRRD, "%s: given nrrd has problems", me);
425
    return 1;
426
  }
427
  if (limnSplineTypeTimeWarp == spline->type) {
428
    if (!( 1 == _ncpt->dim )) {
429
      biffAddf(LIMN, "%s: given nrrd has dimension %d, not 1", me, _ncpt->dim);
430
      return 1;
431
    }
432
    if (!( spline->ncpt->axis[2].size == _ncpt->axis[0].size )) {
433
      biffAddf(LIMN, "%s: have %s time points, but got %s", me,
434
               airSprintSize_t(stmp[0], spline->ncpt->axis[2].size),
435
               airSprintSize_t(stmp[1], _ncpt->axis[0].size));
436
      return 1;
437
    }
438
  } else {
439
    if (!( nrrdSameSize(spline->ncpt, _ncpt, AIR_TRUE) )) {
440
      biffMovef(LIMN, NRRD, "%s: given ncpt doesn't match original one", me);
441
      return 1;
442
    }
443
  }
444
445
  if (limnSplineTypeTimeWarp == spline->type) {
446
    ntmp = nrrdNew();
447
    if (nrrdWrap_va(ntmp, spline->time, nrrdTypeDouble, 1,
448
                    _ncpt->axis[0].size)
449
        || nrrdConvert(ntmp, _ncpt, nrrdTypeDouble)) {
450
      biffMovef(LIMN, NRRD, "%s: trouble copying info", me);
451
      nrrdNix(ntmp); return 1;
452
    }
453
    if (_limnSplineTimeWarpSet(spline)) {
454
      biffAddf(LIMN, "%s: trouble setting time warp", me);
455
      nrrdNix(ntmp); return 1;
456
    }
457
    nrrdNix(ntmp);
458
  } else {
459
    if (nrrdConvert(spline->ncpt, _ncpt, nrrdTypeDouble)) {
460
      biffMovef(LIMN, NRRD, "%s: trouble converting to internal nrrd", me);
461
      return 1;
462
    }
463
  }
464
465
  return 0;
466
}