1 /*
2 *class++
3 * Name:
4 * LutMap
5
6 * Purpose:
7 * Transform 1-dimensional coordinates using a lookup table.
8
9 * Constructor Function:
10 c astLutMap
11 f AST_LUTMAP
12
13 * Description:
14 * A LutMap is a specialised form of Mapping which transforms
15 * 1-dimensional coordinates by using linear interpolation in a
16 * lookup table.
17 *
18 * Each input coordinate value is first scaled to give the index of
19 * an entry in the table by subtracting a starting value (the input
20 * coordinate corresponding to the first table entry) and dividing
21 * by an increment (the difference in input coordinate value
22 * between adjacent table entries).
23 *
24 * The resulting index will usually contain a fractional part, so
25 * the output coordinate value is then generated by interpolating
26 * linearly between the appropriate entries in the table. If the
27 * index lies outside the range of the table, linear extrapolation
28 * is used based on the two nearest entries (i.e. the two entries
29 * at the start or end of the table, as appropriate). If either of the
30 * entries used for the interplation has a value of AST__BAD, then the
31 * interpolated value is returned as AST__BAD.
32 *
33 * If the lookup table entries increase or decrease monotonically
34 * (ignoring any flat sections), then the inverse transformation may
35 * also be performed.
36
37 * Inheritance:
38 * The LutMap class inherits from the Mapping class.
39
40 * Attributes:
41 * In addition to those attributes common to all Mappings, every
42 * LutMap also has the following attributes:
43 *
44 * - LutInterp: The interpolation method to use between table entries.
45
46 * Functions:
47 c The LutMap class does not define any new functions beyond those
48 f The LutMap class does not define any new routines beyond those
49 * which are applicable to all Mappings.
50
51 * Copyright:
52 * Copyright (C) 1997-2006 Council for the Central Laboratory of the
53 * Research Councils
54 * Copyright (C) 2007-2011 Science & Technology Facilities Council.
55 * All Rights Reserved.
56
57 * Licence:
58 * This program is free software: you can redistribute it and/or
59 * modify it under the terms of the GNU Lesser General Public
60 * License as published by the Free Software Foundation, either
61 * version 3 of the License, or (at your option) any later
62 * version.
63 *
64 * This program is distributed in the hope that it will be useful,
65 * but WITHOUT ANY WARRANTY; without even the implied warranty of
66 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
67 * GNU Lesser General Public License for more details.
68 *
69 * You should have received a copy of the GNU Lesser General
70 * License along with this program. If not, see
71 * <http://www.gnu.org/licenses/>.
72
73 * Authors:
74 * RFWS: R.F. Warren-Smith (Starlink)
75 * DSB: David S. Berry (JAC, UCLan)
76
77 * History:
78 * 8-JUL-1997 (RFWS):
79 * Original version.
80 * 10-JUL-1997 (RFWS):
81 * Added the MapMerge function.
82 * 8-JAN-2003 (DSB):
83 * Changed private InitVtab method to protected astInitLutMapVtab
84 * method.
85 * 12-JAN-2004 (DSB):
86 * Check for AST__BAD values in the supplied lut array.
87 * 17-MAR-2006 (DSB):
88 * - MapMerge changed so that a LutMap will cancel with its own
89 * inverse.
90 * - Added attribute LutInterp
91 * 10-MAY-2006 (DSB):
92 * Override astEqual.
93 * 4-OCT-2006 (DSB):
94 * - Correct "mintick" to "lutinterp" in SetAttrib.
95 * - Do not include bad values in the dumped LUT array.
96 * 8-NOV-2007 (DSB):
97 * - Take account of the requested invert flag when comparing two
98 * neighbouring LutMaps for equality in MapMerge.
99 * 19-NOV-2010 (DSB):
100 * Added (protected) astGetLutMapInfo function.
101 * 24-JAN-2011 (DSB):
102 * Implement an inverse transformation even if the coordinate
103 * array contains sections of equal or bad values. The inverse
104 * transformation will generate bad values if used within such
105 * regions of the coordinate array.
106 * 6-JUL-2011 (DSB):
107 * Avoid indexing the lut array beyond the end when doing an
108 * inverse transform.
109 * 2-OCT-2012 (DSB):
110 * Check for Infs as well as NaNs.
111 *class--
112 */
113
114 /* Module Macros. */
115 /* ============== */
116 /* Set the name of the class we are implementing. This indicates to
117 the header files that define class interfaces that they should make
118 "protected" symbols available. */
119 #define astCLASS LutMap
120
121 #define LINEAR 0
122 #define NEAR 1
123
124 /* Include files. */
125 /* ============== */
126 /* Interface definitions. */
127 /* ---------------------- */
128
129 #include "globals.h" /* Thread-safe global data access */
130 #include "error.h" /* Error reporting facilities */
131 #include "memory.h" /* Memory management facilities */
132 #include "object.h" /* Base Object class */
133 #include "pointset.h" /* Sets of points/coordinates */
134 #include "mapping.h" /* Coordinate mappings (parent class) */
135 #include "winmap.h" /* Linear mappings between windows */
136 #include "channel.h" /* I/O channels */
137 #include "unitmap.h" /* Unit mappings */
138 #include "lutmap.h" /* Interface definition for this class */
139 #include "globals.h" /* Thread-safe global data access */
140
141 /* Error code definitions. */
142 /* ----------------------- */
143 #include "ast_err.h" /* AST error codes */
144
145 /* C header files. */
146 /* --------------- */
147 #include <float.h>
148 #include <math.h>
149 #include <limits.h>
150 #include <stdarg.h>
151 #include <stddef.h>
152 #include <stdio.h>
153 #include <string.h>
154
155 /* Module Variables. */
156 /* ================= */
157
158 /* Address of this static variable is used as a unique identifier for
159 member of this class. */
160 static int class_check;
161
162 /* Pointers to parent class methods which are extended by this class. */
163 static AstPointSet *(* parent_transform)( AstMapping *, AstPointSet *, int, AstPointSet *, int * );
164 static const char *(* parent_getattrib)( AstObject *, const char *, int * );
165 static int (* parent_testattrib)( AstObject *, const char *, int * );
166 static void (* parent_clearattrib)( AstObject *, const char *, int * );
167 static void (* parent_setattrib)( AstObject *, const char *, int * );
168
169 /* Define macros for accessing each item of thread specific global data. */
170 #ifdef THREAD_SAFE
171
172 /* Define how to initialise thread-specific globals. */
173 #define GLOBAL_inits \
174 globals->Class_Init = 0; \
175 globals->GetAttrib_Buff[ 0 ] = 0;
176
177 /* Create the function that initialises global data for this module. */
178 astMAKE_INITGLOBALS(LutMap)
179
180 /* Define macros for accessing each item of thread specific global data. */
181 #define class_init astGLOBAL(LutMap,Class_Init)
182 #define class_vtab astGLOBAL(LutMap,Class_Vtab)
183 #define getattrib_buff astGLOBAL(LutMap,GetAttrib_Buff)
184
185
186
187 /* If thread safety is not needed, declare and initialise globals at static
188 variables. */
189 #else
190
191 static char getattrib_buff[ 101 ];
192
193
194 /* Define the class virtual function table and its initialisation flag
195 as static variables. */
196 static AstLutMapVtab class_vtab; /* Virtual function table */
197 static int class_init = 0; /* Virtual function table initialised? */
198
199 #endif
200
201 /* External Interface Function Prototypes. */
202 /* ======================================= */
203 /* The following functions have public prototypes only (i.e. no
204 protected prototypes), so we must provide local prototypes for use
205 within this module. */
206 AstLutMap *astLutMapId_( int, const double [], double, double, const char *, ... );
207
208 /* Prototypes for Private Member Functions. */
209 /* ======================================== */
210 static AstPointSet *Transform( AstMapping *, AstPointSet *, int, AstPointSet *, int * );
211 static int GetLinear( AstMapping *, int * );
212 static int GetMonotonic( int, const double *, int *, double **, int **, int **, int * );
213 static int MapMerge( AstMapping *, int, int, int *, AstMapping ***, int **, int * );
214 static void Copy( const AstObject *, AstObject *, int * );
215 static void Delete( AstObject *, int * );
216 static void Dump( AstObject *, AstChannel *, int * );
217 static int Equal( AstObject *, AstObject *, int * );
218 static double *GetLutMapInfo( AstLutMap *, double *, double *, int *, int * );
219
220 static const char *GetAttrib( AstObject *, const char *, int * );
221 static int TestAttrib( AstObject *, const char *, int * );
222 static void ClearAttrib( AstObject *, const char *, int * );
223 static void SetAttrib( AstObject *, const char *, int * );
224
225 static int GetLutInterp( AstLutMap *, int * );
226 static int TestLutInterp( AstLutMap *, int * );
227 static void ClearLutInterp( AstLutMap *, int * );
228 static void SetLutInterp( AstLutMap *, int, int * );
229
230 /* Member functions. */
231 /* ================= */
ClearAttrib(AstObject * this_object,const char * attrib,int * status)232 static void ClearAttrib( AstObject *this_object, const char *attrib, int *status ) {
233 /*
234 * Name:
235 * ClearAttrib
236
237 * Purpose:
238 * Clear an attribute value for a LutMap.
239
240 * Type:
241 * Private function.
242
243 * Synopsis:
244 * #include "lutmap.h"
245 * void ClearAttrib( AstObject *this, const char *attrib, int *status )
246
247 * Class Membership:
248 * LutMap member function (over-rides the astClearAttrib protected
249 * method inherited from the Mapping class).
250
251 * Description:
252 * This function clears the value of a specified attribute for a
253 * LutMap, so that the default value will subsequently be used.
254
255 * Parameters:
256 * this
257 * Pointer to the LutMap.
258 * attrib
259 * Pointer to a null-terminated string specifying the attribute
260 * name. This should be in lower case with no surrounding white
261 * space.
262 * status
263 * Pointer to the inherited status variable.
264 */
265
266 /* Local Variables: */
267 AstLutMap *this; /* Pointer to the LutMap structure */
268
269 /* Check the global error status. */
270 if ( !astOK ) return;
271
272 /* Obtain a pointer to the LutMap structure. */
273 this = (AstLutMap *) this_object;
274
275 /* Check the attribute name and clear the appropriate attribute. */
276
277 /* LutInterp. */
278 /* ---------- */
279 if ( !strcmp( attrib, "lutinterp" ) ) {
280 astClearLutInterp( this );
281
282 /* If the attribute is still not recognised, pass it on to the parent
283 method for further interpretation. */
284 } else {
285 (*parent_clearattrib)( this_object, attrib, status );
286 }
287 }
288
Equal(AstObject * this_object,AstObject * that_object,int * status)289 static int Equal( AstObject *this_object, AstObject *that_object, int *status ) {
290 /*
291 * Name:
292 * Equal
293
294 * Purpose:
295 * Test if two LutMaps are equivalent.
296
297 * Type:
298 * Private function.
299
300 * Synopsis:
301 * #include "lutmap.h"
302 * int Equal( AstObject *this, AstObject *that, int *status )
303
304 * Class Membership:
305 * LutMap member function (over-rides the astEqual protected
306 * method inherited from the astMapping class).
307
308 * Description:
309 * This function returns a boolean result (0 or 1) to indicate whether
310 * two LutMaps are equivalent.
311
312 * Parameters:
313 * this
314 * Pointer to the first Object (a LutMap).
315 * that
316 * Pointer to the second Object.
317 * status
318 * Pointer to the inherited status variable.
319
320 * Returned Value:
321 * One if the LutMaps are equivalent, zero otherwise.
322
323 * Notes:
324 * - A value of zero will be returned if this function is invoked
325 * with the global status set, or if it should fail for any reason.
326 */
327
328 /* Local Variables: */
329 AstLutMap *that;
330 AstLutMap *this;
331 int i;
332 int nin;
333 int nout;
334 int result;
335
336 /* Initialise. */
337 result = 0;
338
339 /* Check the global error status. */
340 if ( !astOK ) return result;
341
342 /* Obtain pointers to the two LutMap structures. */
343 this = (AstLutMap *) this_object;
344 that = (AstLutMap *) that_object;
345
346 /* Check the second object is a LutMap. We know the first is a
347 LutMap since we have arrived at this implementation of the virtual
348 function. */
349 if( astIsALutMap( that ) ) {
350
351 /* Get the number of inputs and outputs and check they are the same for both. */
352 nin = astGetNin( this );
353 nout = astGetNout( this );
354 if( astGetNin( that ) == nin && astGetNout( that ) == nout ) {
355
356 /* If the Invert flags for the two LutMaps differ, it may still be possible
357 for them to be equivalent. First compare the LutMaps if their Invert
358 flags are the same. In this case all the attributes of the two LutMaps
359 must be identical. */
360 if( astGetInvert( this ) == astGetInvert( that ) ) {
361
362 if( astEQUAL( this->start, that->start ) &&
363 astEQUAL( this->inc, that->inc ) &&
364 this->nlut == that->nlut &&
365 this->lutinterp == that->lutinterp ){
366
367 result = 1;
368 for( i = 0; i < this->nlut; i++ ) {
369 if( !astEQUAL( (this->lut)[ i ], (that->lut)[ i ] ) ) {
370 result = 0;
371 break;
372 }
373 }
374 }
375
376 /* If the Invert flags for the two LutMaps differ, the attributes of the two
377 LutMaps must be inversely related to each other. */
378 } else {
379
380 /* In the specific case of a LutMap, Invert flags must be equal. */
381 result = 0;
382 }
383 }
384 }
385
386 /* If an error occurred, clear the result value. */
387 if ( !astOK ) result = 0;
388
389 /* Return the result, */
390 return result;
391 }
392
GetAttrib(AstObject * this_object,const char * attrib,int * status)393 static const char *GetAttrib( AstObject *this_object, const char *attrib, int *status ) {
394 /*
395 * Name:
396 * GetAttrib
397
398 * Purpose:
399 * Get the value of a specified attribute for a LutMap.
400
401 * Type:
402 * Private function.
403
404 * Synopsis:
405 * #include "lutmap.h"
406 * const char *GetAttrib( AstObject *this, const char *attrib, int *status )
407
408 * Class Membership:
409 * LutMap member function (over-rides the protected astGetAttrib
410 * method inherited from the Mapping class).
411
412 * Description:
413 * This function returns a pointer to the value of a specified
414 * attribute for a LutMap, formatted as a character string.
415
416 * Parameters:
417 * this
418 * Pointer to the LutMap.
419 * attrib
420 * Pointer to a null-terminated string containing the name of
421 * the attribute whose value is required. This name should be in
422 * lower case, with all white space removed.
423 * status
424 * Pointer to the inherited status variable.
425
426 * Returned Value:
427 * - Pointer to a null-terminated string containing the attribute
428 * value.
429
430 * Notes:
431 * - The returned string pointer may point at memory allocated
432 * within the LutMap, or at static memory. The contents of the
433 * string may be over-written or the pointer may become invalid
434 * following a further invocation of the same function or any
435 * modification of the LutMap. A copy of the string should
436 * therefore be made if necessary.
437 * - A NULL pointer will be returned if this function is invoked
438 * with the global error status set, or if it should fail for any
439 * reason.
440 */
441
442 /* Local Variables: */
443 astDECLARE_GLOBALS /* Pointer to thread-specific global data */
444 AstLutMap *this; /* Pointer to the LutMap structure */
445 const char *result; /* Pointer value to return */
446 int lutinterp; /* LutInterp attribute value */
447
448 /* Initialise. */
449 result = NULL;
450
451 /* Check the global error status. */
452 if ( !astOK ) return result;
453
454 /* Get a pointer to the thread specific global data structure. */
455 astGET_GLOBALS(this_object);
456
457 /* Obtain a pointer to the LutMap structure. */
458 this = (AstLutMap *) this_object;
459
460 /* Compare "attrib" with each recognised attribute name in turn,
461 obtaining the value of the required attribute. If necessary, write
462 the value into "getattrib_buff" as a null-terminated string in an appropriate
463 format. Set "result" to point at the result string. */
464
465 /* LutInterp. */
466 /* ---------- */
467 if ( !strcmp( attrib, "lutinterp" ) ) {
468 lutinterp = astGetLutInterp( this );
469 if ( astOK ) {
470 (void) sprintf( getattrib_buff, "%d", lutinterp );
471 result = getattrib_buff;
472 }
473
474 /* If the attribute name was not recognised, pass it on to the parent
475 method for further interpretation. */
476 } else {
477 result = (*parent_getattrib)( this_object, attrib, status );
478 }
479
480 /* Return the result. */
481 return result;
482
483 }
484
GetLinear(AstMapping * this_mapping,int * status)485 static int GetLinear( AstMapping *this_mapping, int *status ) {
486 /*
487 * Name:
488 * GetLinear
489
490 * Purpose:
491 * Determine if a LutMap implements a linear coordinate transformation.
492
493 * Type:
494 * Private function.
495
496 * Synopsis:
497 * #include "lutmap.h"
498 * int GetLinear( AstMapping *this, int *status )
499
500 * Class Membership:
501 * LutMap member function.
502
503 * Description:
504 * This function returns a boolean value to indicate if the LutMap
505 * supplied is equivalent to a linear coordinate
506 * transformation. This will be the case if the lookup table
507 * elements increase or decrease linearly.
508
509 * Parameters:
510 * this
511 * Pointer to the LutMap.
512 * status
513 * Pointer to the inherited status variable.
514
515 * Notes:
516 * - A value of zero will be returned if this function is invoked
517 * with the global error status set, or if it should fail for any
518 * reason.
519 */
520
521 /* Local Variables: */
522 AstLutMap *this; /* Pointer to the LutMap structure */
523 double *lut; /* Pointer to the lookup table */
524 double fract; /* Fractional position within table */
525 double hi; /* Largest value */
526 double interp; /* Interpolated value */
527 double lo; /* Smallest value */
528 double tol1; /* First tolerance estimate */
529 double tol2; /* Second tolerance estimate */
530 double tol; /* Tolerance value used */
531 int ilut; /* Loop counter for table elements */
532 int linear; /* Result to be returned */
533 int nlut; /* Number of lookup table elements */
534
535 /* Initialise. */
536 linear = 0;
537
538 /* Check the global error status. */
539 if ( !astOK ) return linear;
540
541 /* Obtain a pointer to the LutMap structure. */
542 this = (AstLutMap *) this_mapping;
543
544 /* Nearest neighbour LutMaps are not considered to be linear because of
545 the discontinuities at the start and end of the table. */
546 if( astGetLutInterp( this ) != NEAR ) {
547
548 /* Obtain the lookup table details. */
549 lut = this->lut;
550 nlut = this->nlut;
551
552 /* Loop to identify the largest and smallest values in the lookup
553 table. */
554 lo = DBL_MAX;
555 hi = -DBL_MAX;
556 for ( ilut = 0; ilut < nlut; ilut++ ) {
557 if ( lut[ ilut ] > hi ) hi = lut[ ilut ];
558 if ( lut[ ilut ] < lo ) lo = lut[ ilut ];
559 }
560
561 /* Check if the values are all the same (this makes the LutMap
562 linear, although it will have no inverse). */
563 linear = ( hi == lo );
564 if ( !linear ) {
565
566 /* Form a tolerance estimate based on the overall range of values in
567 the lookup table. */
568 tol1 = fabs( hi - lo ) * DBL_EPSILON;
569
570 /* Now loop to inspect all the lookup table elements except the first
571 and last. */
572 linear = 1;
573 for ( ilut = 1; ilut < ( nlut - 1 ); ilut++ ) {
574
575 /* Calculate the fractional position of the current element within the
576 table. */
577 fract = ( (double) ilut ) / ( (double) ( nlut - 1 ) );
578
579 /* Calculate the value it should have if the table is linear by
580 interpolating between the first and last values. */
581 interp = lut[ 0 ] * ( 1.0 - fract ) + lut[ nlut - 1 ] * fract;
582
583 /* Form a second tolerance estimate from this interpolated
584 value. Select whichever tolerance estimate is larger (this avoids
585 problems when values are near zero). */
586 tol2 = fabs( interp ) * DBL_EPSILON;
587 tol = ( tol1 > tol2 ) ? tol1 : tol2;
588
589 /* Test for linearity within a small multiple of the tolerance. */
590 linear = ( fabs( lut[ ilut ] - interp ) <= ( 2.0 * tol ) );
591 if ( !linear ) break;
592 }
593 }
594 }
595
596 /* Return the result. */
597 return linear;
598 }
599
GetLutMapInfo(AstLutMap * this,double * start,double * inc,int * nlut,int * status)600 static double *GetLutMapInfo( AstLutMap *this, double *start, double *inc,
601 int *nlut, int *status ){
602 /*
603 * Name:
604 * GetLutMapInfo
605
606 * Purpose:
607 * Return information about a LutMap.
608
609 * Type:
610 * Protected virtual function.
611
612 * Synopsis:
613 * #include "permmap.h"
614 * double *astGetLutMapInfo( AstLutMap *this, double *start, double *inc,
615 * int *nlut, int *status )
616
617 * Class Membership:
618 * LutMap method
619
620 * Description:
621 * This function returns information about the supplied LutMap.
622
623 * Parameters:
624 * this
625 * Pointer to the LutMap.
626 * start
627 * Pointer to a double in which to return the "start" value
628 * supplied when the LutMap was created.
629 * inc
630 * Pointer to a double in which to return the "inc" value
631 * supplied when the LutMap was created.
632 * nlut
633 * Pointer to a double in which to return the number of values in
634 * the LUT.
635 * status
636 * Pointer to the inherited status variable.
637
638 * Returned Value:
639 * A pointer to a dynamically allocated array holding a copy of the
640 * look-up table. This is an array of "nlut" elements, giving the
641 * output values for input values "start", "start+inc", "start+2*inc",
642 * etc. The pointer should be freed using astFree when no longer
643 * needed.
644
645 * Notes:
646 * - A value of NULL will be returned if this function is invoked
647 * with the global error status set, or if it should fail for any
648 * reason.
649 */
650
651 /* Check the global error status. */
652 if ( !astOK ) return NULL;
653
654 /* Store the scalar values. */
655 *start = this->start;
656 *inc = this->inc;
657 *nlut = this->nlut;
658
659 /* Return a copy of the look up table. */
660 return astStore( NULL, this->lut, sizeof( double )*this->nlut );
661 }
662
GetMonotonic(int nlut,const double * lut,int * nluti,double ** luti,int ** flagsi,int ** indexi,int * status)663 static int GetMonotonic( int nlut, const double *lut, int *nluti, double **luti,
664 int **flagsi, int **indexi, int *status ) {
665 /*
666 * Name:
667 * GetMonotonic
668
669 * Purpose:
670 * Determine if a array is monotonic increasing or decreasing.
671
672 * Type:
673 * Private function.
674
675 * Synopsis:
676 * #include "lutmap.h"
677 * int GetMonotonic( int nlut, const double *lut, int *nluti, double **luti,
678 * int **flagsi, int **indexi, int *status )
679
680 * Class Membership:
681 * LutMap member function.
682
683 * Description:
684 * This function returns a flag indiciating the supplied array is
685 * monotonic increasing, monotonic decreasing, or non-monotonic.
686 * Sections of equal or bad values do not invalidate an otherwise
687 * monotonic array.
688 *
689 * It also returns information needed to implement the inverse
690 * transformation of a LutMap.
691
692 * Parameters:
693 * nlut
694 * The length of the array.
695 * lut
696 * The array to check.
697 * nluti
698 * Address of an int in which to store the length of the returned
699 * "luti" and "flags" arrays.
700 * luti
701 * Address at which to return a pointer to a newly allocated array
702 * containing "*nluti" elements. This is a copy of the supplied
703 * "lut" array but with any bad or NaN values omitted. Subsequent
704 * elements are shuffled down to fill the holes left by removing
705 * these bad values. A NULL pointer is returned if there are no bad
706 * values in "lut".
707 * flagsi
708 * Address at which to return a pointer to a newly allocated array
709 * containing "*nluti" elements. Each element is non-zero if the
710 * corresponding value stored in "luti" was adjacent to a bad value
711 * in the supplied "lut" array. A NULL pointer is returned if there
712 * are no bad values in "lut".
713 * indexi
714 * Address at which to return a pointer to a newly allocated array
715 * containing "*nluti" elements. Each element is the index within
716 * "lut" of the corresponding value stored in "luti". A NULL pointer
717 * is returned if there are no bad values in "lut".
718 * status
719 * Pointer to the inherited status variable.
720
721 * Notes:
722 * - A value of zero will be returned if this function is invoked
723 * with the global error status set, or if it should fail for any
724 * reason.
725 */
726
727 /* Local Variables: */
728 const double *p3;
729 double *p1;
730 double lval;
731 int *p2;
732 int *p4;
733 int ilut;
734 int nbad;
735 int result;
736
737 /* Initialise. */
738 result = 0;
739 *nluti = 0;
740 *luti = NULL;
741 *flagsi = NULL;
742 *indexi = NULL;
743
744 /* Check the global error status. */
745 if ( !astOK ) return result;
746
747 /* As yet we have not found a good value. */
748 lval = AST__BAD;
749
750 /* As yet we have not found a bad value. */
751 nbad = 0;
752
753 /* Loop round the supplied array, ignoring bad (AST__BAD or NaN) values. */
754 for ( ilut = 0; ilut < nlut; ilut++ ) {
755 if( !astISBAD( lut[ ilut ] ) ) {
756
757 /* If this is the first good value, record it. */
758 if( lval == AST__BAD ) {
759 lval = lut[ ilut ];
760
761 /* If this is not the first good value, ignore it if it is equal to the
762 previous god value. */
763 } else if( lut[ ilut ] != lval ) {
764
765 /* Set the returned flag on the basis of the first pair of good values. */
766 if( result == 0 ) {
767 result = ( lut[ ilut ] > lval ) ? 1 : -1;
768
769 /* For subsequent pairs of good values, check the pair increases or
770 decreases in the same way as the first pair. Reset the returned value
771 to zero and break if not. */
772 } else if( result == 1 && lut[ ilut ] < lval ) {
773 result = 0;
774 break;
775
776 } else if( result == -1 && lut[ ilut ] >lval ) {
777 result = 0;
778 break;
779 }
780
781 }
782 } else {
783 nbad++;
784 }
785 }
786
787 /* If any bad values were found, we now allocate the required returned
788 arrays. */
789 if( nbad ) {
790 *nluti = nlut - nbad;
791 *luti = astMalloc( sizeof( double )*( *nluti ) );
792 *flagsi = astMalloc( sizeof( double )*( *nluti ) );
793 *indexi = astMalloc( sizeof( double )*( *nluti ) );
794
795 if( astOK ) {
796
797 /* Into "luti" copy all good values from "lut", shuffling down values to
798 fill holes left by bad values. Into "flagsi", store a flag indicating
799 if the corresponding "luti" value was adjacent to a bad value in the
800 full "lut" array. */
801 p1 = *luti;
802 p2 = *flagsi;
803 p3 = lut;
804 p4 = *indexi;
805
806 /* Do the first input point (it has no lower neighbour). */
807 if( !astISBAD( *p3 ) ) {
808 *(p1++) = *p3;
809 *(p2++) = astISBAD( p3[ +1 ] );
810 *(p4++) = 0;
811 }
812
813 /* Do all remaining input points except for the last one. */
814 for ( ilut = 1,p3++; ilut < nlut-1; ilut++,p3++ ) {
815 if( !astISBAD( *p3 ) ) {
816 *(p1++) = *p3;
817 *(p2++) = astISBAD( p3[ -1 ] ) || astISBAD( p3[ +1 ] );
818 *(p4++) = ilut;
819 }
820 }
821
822 /* Do the last input point (it has no upper neighbour). */
823 if( !astISBAD( *p3 ) ) {
824 *p1 = *p3;
825 *p2 = astISBAD( p3[ -1 ] );
826 *p4 = ilut;
827 }
828 }
829 }
830
831
832 /* Return the result. */
833 return result;
834 }
835
astInitLutMapVtab_(AstLutMapVtab * vtab,const char * name,int * status)836 void astInitLutMapVtab_( AstLutMapVtab *vtab, const char *name, int *status ) {
837 /*
838 *+
839 * Name:
840 * astInitLutMapVtab
841
842 * Purpose:
843 * Initialise a virtual function table for a LutMap.
844
845 * Type:
846 * Protected function.
847
848 * Synopsis:
849 * #include "lutmap.h"
850 * void astInitLutMapVtab( AstLutMapVtab *vtab, const char *name )
851
852 * Class Membership:
853 * LutMap vtab initialiser.
854
855 * Description:
856 * This function initialises the component of a virtual function
857 * table which is used by the LutMap class.
858
859 * Parameters:
860 * vtab
861 * Pointer to the virtual function table. The components used by
862 * all ancestral classes will be initialised if they have not already
863 * been initialised.
864 * name
865 * Pointer to a constant null-terminated character string which contains
866 * the name of the class to which the virtual function table belongs (it
867 * is this pointer value that will subsequently be returned by the Object
868 * astClass function).
869 *-
870 */
871
872 /* Local Variables: */
873 astDECLARE_GLOBALS /* Pointer to thread-specific global data */
874 AstObjectVtab *object; /* Pointer to Object component of Vtab */
875 AstMappingVtab *mapping; /* Pointer to Mapping component of Vtab */
876
877 /* Check the local error status. */
878 if ( !astOK ) return;
879
880 /* Get a pointer to the thread specific global data structure. */
881 astGET_GLOBALS(NULL);
882
883 /* Initialize the component of the virtual function table used by the
884 parent class. */
885 astInitMappingVtab( (AstMappingVtab *) vtab, name );
886
887 /* Store a unique "magic" value in the virtual function table. This
888 will be used (by astIsALutMap) to determine if an object belongs
889 to this class. We can conveniently use the address of the (static)
890 class_check variable to generate this unique value. */
891 vtab->id.check = &class_check;
892 vtab->id.parent = &(((AstMappingVtab *) vtab)->id);
893
894 /* Initialise member function pointers. */
895 /* ------------------------------------ */
896 /* Store pointers to the member functions (implemented here) that
897 provide virtual methods for this class. */
898 vtab->ClearLutInterp = ClearLutInterp;
899 vtab->GetLutInterp = GetLutInterp;
900 vtab->SetLutInterp = SetLutInterp;
901 vtab->TestLutInterp = TestLutInterp;
902 vtab->GetLutMapInfo = GetLutMapInfo;
903
904 /* Save the inherited pointers to methods that will be extended, and
905 replace them with pointers to the new member functions. */
906 object = (AstObjectVtab *) vtab;
907 mapping = (AstMappingVtab *) vtab;
908
909 parent_clearattrib = object->ClearAttrib;
910 object->ClearAttrib = ClearAttrib;
911 parent_getattrib = object->GetAttrib;
912 object->GetAttrib = GetAttrib;
913 parent_setattrib = object->SetAttrib;
914 object->SetAttrib = SetAttrib;
915 parent_testattrib = object->TestAttrib;
916 object->TestAttrib = TestAttrib;
917
918 parent_transform = mapping->Transform;
919 mapping->Transform = Transform;
920
921 /* Store replacement pointers for methods which will be over-ridden by
922 new member functions implemented here. */
923 object->Equal = Equal;
924 mapping->MapMerge = MapMerge;
925
926 /* Declare the class dump, copy and delete functions.*/
927 astSetDump( vtab, Dump, "LutMap",
928 "Map 1-d coordinates using a lookup table" );
929 astSetCopy( (AstObjectVtab *) vtab, Copy );
930 astSetDelete( (AstObjectVtab *) vtab, Delete );
931
932 /* If we have just initialised the vtab for the current class, indicate
933 that the vtab is now initialised, and store a pointer to the class
934 identifier in the base "object" level of the vtab. */
935 if( vtab == &class_vtab ) {
936 class_init = 1;
937 astSetVtabClassIdentifier( vtab, &(vtab->id) );
938 }
939 }
940
MapMerge(AstMapping * this,int where,int series,int * nmap,AstMapping *** map_list,int ** invert_list,int * status)941 static int MapMerge( AstMapping *this, int where, int series, int *nmap,
942 AstMapping ***map_list, int **invert_list, int *status ) {
943 /*
944 * Name:
945 * MapMerge
946
947 * Purpose:
948 * Simplify a sequence of Mappings containing a LutMap.
949
950 * Type:
951 * Private function.
952
953 * Synopsis:
954 * #include "mapping.h"
955 * int MapMerge( AstMapping *this, int where, int series, int *nmap,
956 * AstMapping ***map_list, int **invert_list, int *status )
957
958 * Class Membership:
959 * LutMap method (over-rides the protected astMapMerge method
960 * inherited from the Mapping class).
961
962 * Description:
963 * This function attempts to simplify a sequence of Mappings by
964 * merging a nominated LutMap in the sequence with its neighbours,
965 * so as to shorten the sequence if possible.
966 *
967 * In many cases, simplification will not be possible and the
968 * function will return -1 to indicate this, without further
969 * action.
970 *
971 * In most cases of interest, however, this function will either
972 * attempt to replace the nominated LutMap with one which it
973 * considers simpler, or to merge it with the Mappings which
974 * immediately precede it or follow it in the sequence (both will
975 * normally be considered). This is sufficient to ensure the
976 * eventual simplification of most Mapping sequences by repeated
977 * application of this function.
978 *
979 * In some cases, the function may attempt more elaborate
980 * simplification, involving any number of other Mappings in the
981 * sequence. It is not restricted in the type or scope of
982 * simplification it may perform, but will normally only attempt
983 * elaborate simplification in cases where a more straightforward
984 * approach is not adequate.
985
986 * Parameters:
987 * this
988 * Pointer to the nominated LutMap which is to be merged with
989 * its neighbours. This should be a cloned copy of the LutMap
990 * pointer contained in the array element "(*map_list)[where]"
991 * (see below). This pointer will not be annulled, and the
992 * LutMap it identifies will not be modified by this function.
993 * where
994 * Index in the "*map_list" array (below) at which the pointer
995 * to the nominated LutMap resides.
996 * series
997 * A non-zero value indicates that the sequence of Mappings to
998 * be simplified will be applied in series (i.e. one after the
999 * other), whereas a zero value indicates that they will be
1000 * applied in parallel (i.e. on successive sub-sets of the
1001 * input/output coordinates).
1002 * nmap
1003 * Address of an int which counts the number of Mappings in the
1004 * sequence. On entry this should be set to the initial number
1005 * of Mappings. On exit it will be updated to record the number
1006 * of Mappings remaining after simplification.
1007 * map_list
1008 * Address of a pointer to a dynamically allocated array of
1009 * Mapping pointers (produced, for example, by the astMapList
1010 * method) which identifies the sequence of Mappings. On entry,
1011 * the initial sequence of Mappings to be simplified should be
1012 * supplied.
1013 *
1014 * On exit, the contents of this array will be modified to
1015 * reflect any simplification carried out. Any form of
1016 * simplification may be performed. This may involve any of: (a)
1017 * removing Mappings by annulling any of the pointers supplied,
1018 * (b) replacing them with pointers to new Mappings, (c)
1019 * inserting additional Mappings and (d) changing their order.
1020 *
1021 * The intention is to reduce the number of Mappings in the
1022 * sequence, if possible, and any reduction will be reflected in
1023 * the value of "*nmap" returned. However, simplifications which
1024 * do not reduce the length of the sequence (but improve its
1025 * execution time, for example) may also be performed, and the
1026 * sequence might conceivably increase in length (but normally
1027 * only in order to split up a Mapping into pieces that can be
1028 * more easily merged with their neighbours on subsequent
1029 * invocations of this function).
1030 *
1031 * If Mappings are removed from the sequence, any gaps that
1032 * remain will be closed up, by moving subsequent Mapping
1033 * pointers along in the array, so that vacated elements occur
1034 * at the end. If the sequence increases in length, the array
1035 * will be extended (and its pointer updated) if necessary to
1036 * accommodate any new elements.
1037 *
1038 * Note that any (or all) of the Mapping pointers supplied in
1039 * this array may be annulled by this function, but the Mappings
1040 * to which they refer are not modified in any way (although
1041 * they may, of course, be deleted if the annulled pointer is
1042 * the final one).
1043 * invert_list
1044 * Address of a pointer to a dynamically allocated array which,
1045 * on entry, should contain values to be assigned to the Invert
1046 * attributes of the Mappings identified in the "*map_list"
1047 * array before they are applied (this array might have been
1048 * produced, for example, by the astMapList method). These
1049 * values will be used by this function instead of the actual
1050 * Invert attributes of the Mappings supplied, which are
1051 * ignored.
1052 *
1053 * On exit, the contents of this array will be updated to
1054 * correspond with the possibly modified contents of the
1055 * "*map_list" array. If the Mapping sequence increases in
1056 * length, the "*invert_list" array will be extended (and its
1057 * pointer updated) if necessary to accommodate any new
1058 * elements.
1059 * status
1060 * Pointer to the inherited status variable.
1061
1062 * Returned Value:
1063 * If simplification was possible, the function returns the index
1064 * in the "map_list" array of the first element which was
1065 * modified. Otherwise, it returns -1 (and makes no changes to the
1066 * arrays supplied).
1067
1068 * Notes:
1069 * - A value of -1 will be returned if this function is invoked
1070 * with the global error status set, or if it should fail for any
1071 * reason.
1072 */
1073
1074 /* Local Variables: */
1075 AstLutMap *map; /* Pointer to LutMap */
1076 AstLutMap *neb; /* Pointer to neighbouring LutMap */
1077 AstMapping *new; /* Pointer to replacement Mapping */
1078 double a1; /* First input coordinate value */
1079 double a2; /* Second input coordinate value */
1080 double b1; /* First output coordinate value */
1081 double b2; /* Second output coordinate value */
1082 int equal; /* Are LutMaps equal? */
1083 int i; /* Mapping index */
1084 int ilo; /* Index of lower LutMap */
1085 int invneb; /* Should the neigbour be used inverted? */
1086 int old_inv; /* Original Invert value for neigbour */
1087 int result; /* Result value to return */
1088 int simpler; /* Mapping simplified? */
1089
1090 /* Initialise the returned result. */
1091 result = -1;
1092
1093 /* Check the global error status. */
1094 if ( !astOK ) return result;
1095
1096 /* Obtain a pointer to the nominated LutMap. */
1097 map = (AstLutMap *) ( *map_list )[ where ];
1098
1099 /* See if the LutMap is linear. If so, it can probably be
1100 simplified. */
1101 simpler = GetLinear( (AstMapping *) map, status );
1102 if ( simpler ) {
1103
1104 /* Obtain the range of input values corresponding to the first and
1105 last lookup table elements. */
1106 a1 = map->start;
1107 a2 = map->start + map->inc * ( map->nlut - 1 );
1108
1109 /* Obtain the corresponding range of output values and check these
1110 values are not the same. */
1111 b1 = map->lut[ 0 ];
1112 b2 = map->lut[ map->nlut - 1 ];
1113 if ( b1 != b2 ) {
1114
1115 /* Create a new WinMap that implements an equivalent linear Mapping,
1116 allowing for the invert flag associated with the LutMap. */
1117 if ( !( *invert_list )[ where ] ) {
1118 new = (AstMapping *) astWinMap( 1, &a1, &a2, &b1, &b2, "", status );
1119 } else {
1120 new = (AstMapping *) astWinMap( 1, &b1, &b2, &a1, &a2, "", status );
1121 }
1122
1123 /* If OK, annul the original LutMap pointer and substitute the new
1124 one. Also clear the associated invert flag. */
1125 if ( astOK ) {
1126 (void) astAnnul( ( *map_list )[ where ] );
1127 ( *map_list )[ where ] = new;
1128 ( *invert_list )[ where ] = 0;
1129
1130 /* Assign the result value. */
1131 result = where;
1132 }
1133 }
1134
1135 /* Otherwise, see if the LutMap is in series with its own inverse. If so
1136 the pair of LutMaps can be replaced by a UnitMap. */
1137 } else if( series ) {
1138
1139 /* Is the higher neighbour a LutMap? If so get a pointer to it, and
1140 note the index of the lower of the two adjacent LutMaps. */
1141 if( where < ( *nmap - 1 ) &&
1142 astIsALutMap( ( *map_list )[ where + 1 ] ) ){
1143 neb = (AstLutMap *) ( *map_list )[ where + 1 ];
1144 invneb = ( *invert_list )[ where + 1 ];
1145 ilo = where;
1146
1147 /* If not, is the lower neighbour a LutMap? If so get a pointer to it,
1148 and note the index of the lower of the two adjacent LutMaps. */
1149 } else if( where > 0 &&
1150 astIsALutMap( ( *map_list )[ where - 1 ] ) ){
1151 neb = (AstLutMap *) ( *map_list )[ where - 1 ];
1152 invneb = ( *invert_list )[ where - 1 ];
1153 ilo = where - 1;
1154
1155 } else {
1156 neb = NULL;
1157 }
1158
1159 /* If a neighbouring LutMap was found, we can replace the pair by a
1160 UnitMap if the two LutMaps are equal but have opposite values for
1161 their Invert flags. Temporarily invert the neighbour, then compare
1162 the two LutMaps for equality, then re-invert the neighbour. */
1163 if( neb ) {
1164 old_inv = astGetInvert( neb );
1165 astSetInvert( neb, invneb );
1166 astInvert( neb );
1167 equal = astEqual( map, neb );
1168 astSetInvert( neb, old_inv );
1169
1170 /* If the two LutMaps are equal but opposite, annul the first of the two
1171 Mappings, and replace it with a UnitMap. Also set the invert flag. */
1172 if( equal ) {
1173 new = (AstMapping *) astUnitMap( 1, "", status );
1174 (void) astAnnul( ( *map_list )[ ilo ] );
1175 ( *map_list )[ ilo ] = new;
1176 ( *invert_list )[ ilo ] = 0;
1177
1178 /* Annul the second of the two Mappings, and shuffle down the rest of the
1179 list to fill the gap. */
1180 (void) astAnnul( ( *map_list )[ ilo + 1 ] );
1181 for ( i = ilo + 2; i < *nmap; i++ ) {
1182 ( *map_list )[ i - 1 ] = ( *map_list )[ i ];
1183 ( *invert_list )[ i - 1 ] = ( *invert_list )[ i ];
1184 }
1185
1186 /* Clear the vacated element at the end. */
1187 ( *map_list )[ *nmap - 1 ] = NULL;
1188 ( *invert_list )[ *nmap - 1 ] = 0;
1189
1190 /* Decrement the Mapping count and return the index of the first
1191 modified element. */
1192 ( *nmap )--;
1193 result = where;
1194 }
1195 }
1196 }
1197
1198 /* If an error occurred, clear the returned result. */
1199 if ( !astOK ) result = -1;
1200
1201 /* Return the result. */
1202 return result;
1203 }
1204
SetAttrib(AstObject * this_object,const char * setting,int * status)1205 static void SetAttrib( AstObject *this_object, const char *setting, int *status ) {
1206 /*
1207 * Name:
1208 * SetAttrib
1209
1210 * Purpose:
1211 * Set an attribute value for a LutMap.
1212
1213 * Type:
1214 * Private function.
1215
1216 * Synopsis:
1217 * #include "lutmap.h"
1218 * void SetAttrib( AstObject *this, const char *setting )
1219
1220 * Class Membership:
1221 * LutMap member function (over-rides the astSetAttrib protected
1222 * method inherited from the Mapping class).
1223
1224 * Description:
1225 * This function assigns an attribute value for a LutMap, the
1226 * attribute and its value being specified by means of a string of
1227 * the form:
1228 *
1229 * "attribute= value "
1230 *
1231 * Here, "attribute" specifies the attribute name and should be in
1232 * lower case with no white space present. The value to the right
1233 * of the "=" should be a suitable textual representation of the
1234 * value to be assigned and this will be interpreted according to
1235 * the attribute's data type. White space surrounding the value is
1236 * only significant for string attributes.
1237
1238 * Parameters:
1239 * this
1240 * Pointer to the LutMap.
1241 * setting
1242 * Pointer to a null-terminated string specifying the new attribute
1243 * value.
1244 */
1245
1246 /* Local Variables: */
1247 AstLutMap *this; /* Pointer to the LutMap structure */
1248 int lutinterp; /* LutInterp attribute value */
1249 int len; /* Length of setting string */
1250 int nc; /* Number of characters read by astSscanf */
1251
1252 /* Check the global error status. */
1253 if ( !astOK ) return;
1254
1255 /* Obtain a pointer to the LutMap structure. */
1256 this = (AstLutMap *) this_object;
1257
1258 /* Obtain the length of the setting string. */
1259 len = (int) strlen( setting );
1260
1261 /* Test for each recognised attribute in turn, using "astSscanf" to parse
1262 the setting string and extract the attribute value (or an offset to
1263 it in the case of string values). In each case, use the value set
1264 in "nc" to check that the entire string was matched. Once a value
1265 has been obtained, use the appropriate method to set it. */
1266
1267 /* LutInterp. */
1268 /* ---------- */
1269 if ( nc = 0,
1270 ( 1 == astSscanf( setting, "lutinterp= %d %n", &lutinterp, &nc ) )
1271 && ( nc >= len ) ) {
1272 astSetLutInterp( this, lutinterp );
1273
1274 /* If the attribute is still not recognised, pass it on to the parent
1275 method for further interpretation. */
1276 } else {
1277 (*parent_setattrib)( this_object, setting, status );
1278 }
1279 }
1280
TestAttrib(AstObject * this_object,const char * attrib,int * status)1281 static int TestAttrib( AstObject *this_object, const char *attrib, int *status ) {
1282 /*
1283 * Name:
1284 * TestAttrib
1285
1286 * Purpose:
1287 * Test if a specified attribute value is set for a LutMap.
1288
1289 * Type:
1290 * Private function.
1291
1292 * Synopsis:
1293 * #include "lutmap.h"
1294 * int TestAttrib( AstObject *this, const char *attrib, int *status )
1295
1296 * Class Membership:
1297 * LutMap member function (over-rides the astTestAttrib protected
1298 * method inherited from the Mapping class).
1299
1300 * Description:
1301 * This function returns a boolean result (0 or 1) to indicate whether
1302 * a value has been set for one of a LutMap's attributes.
1303
1304 * Parameters:
1305 * this
1306 * Pointer to the LutMap.
1307 * attrib
1308 * Pointer to a null-terminated string specifying the attribute
1309 * name. This should be in lower case with no surrounding white
1310 * space.
1311 * status
1312 * Pointer to the inherited status variable.
1313
1314 * Returned Value:
1315 * One if a value has been set, otherwise zero.
1316
1317 * Notes:
1318 * - A value of zero will be returned if this function is invoked
1319 * with the global status set, or if it should fail for any reason.
1320 */
1321
1322 /* Local Variables: */
1323 AstLutMap *this; /* Pointer to the LutMap structure */
1324 int result; /* Result value to return */
1325
1326 /* Initialise. */
1327 result = 0;
1328
1329 /* Check the global error status. */
1330 if ( !astOK ) return result;
1331
1332 /* Obtain a pointer to the LutMap structure. */
1333 this = (AstLutMap *) this_object;
1334
1335 /* Check the attribute name and test the appropriate attribute. */
1336
1337 /* LutInterp. */
1338 /* ---------- */
1339 if ( !strcmp( attrib, "lutinterp" ) ) {
1340 result = astTestLutInterp( this );
1341
1342 /* If the attribute is still not recognised, pass it on to the parent
1343 method for further interpretation. */
1344 } else {
1345 result = (*parent_testattrib)( this_object, attrib, status );
1346 }
1347
1348 /* Return the result, */
1349 return result;
1350 }
1351
Transform(AstMapping * this,AstPointSet * in,int forward,AstPointSet * out,int * status)1352 static AstPointSet *Transform( AstMapping *this, AstPointSet *in,
1353 int forward, AstPointSet *out, int *status ) {
1354 /*
1355 * Name:
1356 * Transform
1357
1358 * Purpose:
1359 * Apply a LutMap to transform a set of points.
1360
1361 * Type:
1362 * Private function.
1363
1364 * Synopsis:
1365 * #include "lutmap.h"
1366 * AstPointSet *Transform( AstMapping *this, AstPointSet *in,
1367 * int forward, AstPointSet *out, int *status )
1368
1369 * Class Membership:
1370 * LutMap member function (over-rides the astTransform protected
1371 * method inherited from the Mapping class).
1372
1373 * Description:
1374 * This function takes a LutMap and a set of points encapsulated
1375 * in a PointSet and transforms the points so as to apply the
1376 * lookup table transformation.
1377
1378 * Parameters:
1379 * this
1380 * Pointer to the LutMap.
1381 * in
1382 * Pointer to the PointSet holding the input coordinate data.
1383 * forward
1384 * A non-zero value indicates that the forward coordinate
1385 * transformation should be applied, while a zero value requests
1386 * the inverse transformation.
1387 * out
1388 * Pointer to a PointSet which will hold the transformed
1389 * (output) coordinate values. A NULL value may also be given,
1390 * in which case a new PointSet will be created by this
1391 * function.
1392 * status
1393 * Pointer to the inherited status variable.
1394
1395 * Returned Value:
1396 * Pointer to the output (possibly new) PointSet.
1397
1398 * Notes:
1399 * - A null pointer will be returned if this function is invoked
1400 * with the global error status set, or if it should fail for any
1401 * reason.
1402 * - The number of coordinate values per point in the input
1403 * PointSet must equal 1.
1404 * - If an output PointSet is supplied, it must have space for
1405 * sufficient number of points (with 1 coordinate value per point)
1406 * to accommodate the result. Any excess space will be ignored.
1407 */
1408
1409 /* Local Variables: */
1410 AstLutMap *map; /* Pointer to LutMap to be applied */
1411 AstPointSet *result; /* Pointer to output PointSet */
1412 double **ptr_in; /* Pointer to input coordinate data */
1413 double **ptr_out; /* Pointer to output coordinate data */
1414 double *lut; /* Pointer to LUT */
1415 double d1; /* Offset to I1 value */
1416 double d2; /* Offset to I2 value */
1417 double fract; /* Fractional interpolation distance */
1418 double scale; /* Normalising scale factor */
1419 double value_in; /* Input coordinate value */
1420 double value_out; /* Output coordinate value */
1421 double x; /* Value normalised to LUT increment */
1422 double xi; /* Integer value of "x" */
1423 int *flags; /* Flags indicating an adjacent bad value */
1424 int *index; /* Translates reduced to original indices */
1425 int i1; /* Lower adjacent LUT index */
1426 int i2; /* Upper adjacent LUT index */
1427 int i; /* New LUT index */
1428 int istart; /* Original LUT index at start of interval */
1429 int ix; /* "x" converted to an int */
1430 int near; /* Perform nearest neighbour interpolation? */
1431 int nlut; /* Number of LUT entries */
1432 int nlutm1; /* Number of LUT entries minus one */
1433 int npoint; /* Number of points */
1434 int ok; /* Lookup table is not flat */
1435 int point; /* Loop counter for points */
1436 int up; /* LUT values are increasing? */
1437
1438 /* Check the global error status. */
1439 if ( !astOK ) return NULL;
1440
1441 /* Obtain a pointer to the LutMap. */
1442 map = (AstLutMap *) this;
1443
1444 /* Apply the parent mapping using the stored pointer to the Transform
1445 member function inherited from the parent Mapping class. This
1446 function validates all arguments and generates an output PointSet
1447 if necessary, but does not actually transform any coordinate
1448 values. */
1449 result = (*parent_transform)( this, in, forward, out, status );
1450
1451 /* We will now extend the parent astTransform method by performing the
1452 calculations needed to generate the output coordinate values. */
1453
1454 /* Determine the numbers of points from the input PointSet and obtain
1455 pointers for accessing the input and output coordinate values. */
1456 npoint = astGetNpoint( in );
1457 ptr_in = astGetPoints( in );
1458 ptr_out = astGetPoints( result );
1459
1460 /* Determine whether to apply the forward or inverse mapping,
1461 according to the direction specified and whether the mapping has
1462 been inverted. */
1463 if ( astGetInvert( map ) ) forward = !forward;
1464
1465 /* Forward transformation. */
1466 /* ----------------------- */
1467 if( astOK ){
1468 if ( forward ) {
1469
1470 /* Obtain lookup table details. */
1471 lut = map->lut;
1472 nlut = map->nlut;
1473 near = ( astGetLutInterp( map ) == NEAR );
1474 nlutm1 = nlut - 1;
1475
1476 /* Calculate the scale factor required. */
1477 scale = 1.0 / map->inc;
1478
1479 /* Loop to transform each input point. */
1480 for ( point = 0; point < npoint; point++ ) {
1481
1482 /* Extract the input coordinate value. */
1483 value_in = ptr_in[ 0 ][ point ];
1484
1485 /* First check if this is the same value as we transformed last. If
1486 so, re-use the last result. */
1487 if ( value_in == map->last_fwd_in ) {
1488 value_out = map->last_fwd_out;
1489
1490 /* Check for bad input coordinates and generate a bad result if
1491 necessary. */
1492 } else if ( value_in == AST__BAD ) {
1493 value_out = AST__BAD;
1494
1495 /* For nearest-neighbour interpolation, return the value of the lookup table
1496 entry corresponding to the input coordinate. */
1497 } else if( near ){
1498 x = ( value_in - map->start ) * scale;
1499 xi = floor( x + 0.5 );
1500 ix = (int) xi;
1501 if ( ix < 0 || ix >= nlut ) {
1502 value_out = AST__BAD;
1503 } else {
1504 value_out = lut[ ix ];
1505 }
1506
1507 /* Otherwise, (for linear interpolation) identify the lookup table entry
1508 corresponding to the input coordinate. */
1509 } else {
1510 x = ( value_in - map->start ) * scale;
1511 xi = floor( x );
1512 ix = (int) xi;
1513
1514 /* If the input value lies below the first lookup table entry,
1515 extrapolate using the first two table values. */
1516 if ( ix < 0 ) {
1517 if( lut[ 0 ] != AST__BAD && lut[ 1 ] != AST__BAD ) {
1518 value_out = lut[ 0 ] + x * ( lut[ 1 ] - lut[ 0 ] );
1519 } else {
1520 value_out = AST__BAD;
1521 }
1522
1523 /* If the input value lies above the last lookup table entry (or equals
1524 it), extrapolate using the last two table values. */
1525 } else if ( ix >= nlutm1 ) {
1526 if( lut[ nlutm1 ] != AST__BAD &&
1527 lut[ nlut - 2 ] != AST__BAD ) {
1528 value_out = lut[ nlutm1 ] +
1529 ( x - (double) ( nlutm1 ) ) *
1530 ( lut[ nlutm1 ] - lut[ nlut - 2 ] );
1531 } else {
1532 value_out = AST__BAD;
1533 }
1534
1535 /* Otherwise, interpolate between the adjacent entries. */
1536 } else {
1537 if( lut[ ix ] != AST__BAD &&
1538 lut[ ix + 1 ] != AST__BAD ) {
1539 fract = x - xi;
1540 value_out = lut[ ix ] * ( 1.0 - fract ) +
1541 lut[ ix + 1 ] * fract;
1542 } else {
1543 value_out = AST__BAD;
1544 }
1545 }
1546 }
1547
1548 /* Assign the output coordinate value. */
1549 ptr_out[ 0 ][ point ] = value_out;
1550
1551 /* Retain the input and output coordinate values for possible re-use
1552 in future. */
1553 map->last_fwd_in = value_in;
1554 map->last_fwd_out = value_out;
1555 }
1556
1557 /* Inverse transformation. */
1558 /* ----------------------- */
1559 } else {
1560
1561 /* Obtain details of the lookup table to be used by the inverse
1562 transformation. This is the same as the forward transformation lookup
1563 table, except that any bad values are omitted. Also, get a pointer to a
1564 array of flags that indicate if the corresponding lookup table entries
1565 were adjacent to a bad value or not in the full lookup table. */
1566 if( map->luti ) {
1567 lut = map->luti;
1568 flags = map->flagsi;
1569 nlut = map->nluti;
1570 index = map->indexi;
1571 } else {
1572 lut = map->lut;
1573 flags = NULL;
1574 nlut = map->nlut;
1575 index = NULL;
1576 }
1577 near = ( astGetLutInterp( map ) == NEAR );
1578 nlutm1 = nlut - 1;
1579
1580 /* Loop to transform each input point. */
1581 for ( point = 0; point < npoint; point++ ) {
1582
1583 /* Extract the input coordinate value. */
1584 value_in = ptr_in[ 0 ][ point ];
1585
1586 /* First check if this is the same value as we transformed last. If
1587 so, re-use the last result. */
1588 if ( value_in == map->last_inv_in ) {
1589 value_out = map->last_inv_out;
1590
1591 /* Check for bad input coordinates and generate a bad result if
1592 necessary. */
1593 } else if ( value_in == AST__BAD ) {
1594 value_out = AST__BAD;
1595
1596 /* Otherwise, we can determine an inverse. Note the inverse transformation
1597 will not be defined, so will not be attempted, unless all the table
1598 entries are monotonically increasing or decreasing, possibly with sections
1599 of equal or bad values. */
1600 } else {
1601 up = ( lut[ nlutm1 ] > lut[ 0 ] );
1602
1603 /* Perform a binary search to identify two adjacent lookup table
1604 elements whose values bracket the input coordinate value. */
1605 i1 = -1;
1606 i2 = nlutm1;
1607 while ( i2 > ( i1 + 1 ) ) {
1608 i = ( i1 + i2 ) / 2;
1609 *( ( ( value_in >= lut[ i ] ) == up ) ? &i1 : &i2 ) = i;
1610 }
1611
1612 /* If the lower table value is equal to the required value, and either of
1613 its neighbours is also equal to the required value, then we have been
1614 asked to find the inverse in a flat region of the table, so return
1615 a bad value. Likewise, if the upper table value is equal to the required
1616 value, and either of its neighbours is also equal to the required value,
1617 then we have been asked to find the inverse in a flat region of the table,
1618 so return a bad value. */
1619 ok = 1;
1620 if( lut[ i1 ] == value_in ) {
1621 if( i1 > 0 && lut[ i1 - 1 ] == value_in ) ok = 0;
1622 if( lut[ i2 ] == value_in ) ok = 0;
1623 } else if( lut[ i2 ] == value_in ) {
1624 if( i2 < nlutm1 && lut[ i2 + 1 ] == value_in ) ok = 0;
1625 if( lut[ i1 ] == value_in ) ok = 0;
1626 }
1627
1628 if( !ok ) {
1629 value_out = AST__BAD;
1630
1631 /* If both of the two table elements were adjacent to a bad value in the
1632 full lookup table, return a bad output value. */
1633 } else if( flags && ( flags[ i1 ] && flags[ i2 ] ) ) {
1634 value_out = AST__BAD;
1635
1636 /* Nearest neighbour interpolation: return the closest of i1 or i2. Return
1637 AST__BAD if the supplied value is less than either or greater than
1638 either. */
1639 } else if( near ) {
1640 d1 = lut[ i1 ] - value_in;
1641 d2 = lut[ i2 ] - value_in;
1642 if( ( d1 > 0.0 && d2 > 0.0 ) ||
1643 ( d1 < 0.0 && d2 < 0.0 ) ) {
1644 value_out = AST__BAD;
1645
1646 } else {
1647
1648 if( fabs( d1 ) < fabs( d2 ) ){
1649 istart = index ? index[ i1 ] : i1;
1650 } else {
1651 istart = index ? index[ i2 ] : i2;
1652 }
1653 value_out = map->start + map->inc * istart;
1654
1655 }
1656
1657 /* Linear interpolation... */
1658 } else {
1659
1660 /* We are interested in the lower bracketing table element. If
1661 necessary, restrict this element's index to lie within the
1662 table. This causes extrapolation to occur (instead of
1663 interpolation) if the input value actually lies outside the range
1664 of the lookup table. */
1665 if ( i1 < 0 ) i1 = 0;
1666 if ( i1 > ( nlut - 2 ) ) i1 = nlut - 2;
1667
1668 /* Interpolate (or extrapolate) to derive the output coordinate
1669 value. */
1670 istart = index ? index[ i1 ] : i1;
1671 value_out = map->start + map->inc * ( (double) istart +
1672 ( ( value_in - lut[ i1 ] ) /
1673 ( lut[ i1 + 1 ] - lut[ i1 ] ) ) );
1674 }
1675 }
1676
1677 /* Assign the output coordinate value. */
1678 ptr_out[ 0 ][ point ] = value_out;
1679
1680 /* Retain the input and output coordinate values for possible re-use
1681 in future. */
1682 map->last_inv_in = value_in;
1683 map->last_inv_out = value_out;
1684 }
1685 }
1686 }
1687
1688 /* Return a pointer to the output PointSet. */
1689 return result;
1690 }
1691
1692 /* Functions which access class attributes. */
1693 /* ---------------------------------------- */
1694 /* Implement member functions to access the attributes associated with
1695 this class using the macros defined for this purpose in the
1696 "object.h" file. For a description of each attribute, see the class
1697 interface (in the associated .h file). */
1698
1699 /*
1700 *att++
1701 * Name:
1702 * LutInterp
1703
1704 * Purpose:
1705 * Look-up table interpolation method.
1706
1707 * Type:
1708 * Public attribute.
1709
1710 * Synopsis:
1711 * Integer.
1712
1713 * Description:
1714 * This attribute indicates the method to be used when finding the
1715 * output value of a LutMap for an input value part way between two
1716 * table entries. If it is set to 0 (the default) then linear
1717 * interpolation is used. Otherwise, nearest neighbour interpolation
1718 * is used.
1719 *
1720 * Using nearest neighbour interpolation causes AST__BAD to be returned
1721 * for any point which falls outside the bounds of the table. Linear
1722 * interpolation results in an extrapolated value being returned based
1723 * on the two end entries in the table.
1724
1725 * Applicability:
1726 * LutMap
1727 * All LutMaps have this attribute.
1728
1729 *att--
1730 */
1731 astMAKE_CLEAR(LutMap,LutInterp,lutinterp,-INT_MAX)
1732 astMAKE_GET(LutMap,LutInterp,int,LINEAR,( ( this->lutinterp == -INT_MAX ) ?
1733 LINEAR : this->lutinterp ))
1734 astMAKE_SET(LutMap,LutInterp,int,lutinterp,(( value == LINEAR ) ? LINEAR : NEAR ))
1735 astMAKE_TEST(LutMap,LutInterp,( this->lutinterp != -INT_MAX ))
1736
1737 /* Copy constructor. */
1738 /* ----------------- */
Copy(const AstObject * objin,AstObject * objout,int * status)1739 static void Copy( const AstObject *objin, AstObject *objout, int *status ) {
1740 /*
1741 * Name:
1742 * Copy
1743
1744 * Purpose:
1745 * Copy constructor for LutMap objects.
1746
1747 * Type:
1748 * Private function.
1749
1750 * Synopsis:
1751 * void Copy( const AstObject *objin, AstObject *objout, int *status )
1752
1753 * Description:
1754 * This function implements the copy constructor for LutMap objects.
1755
1756 * Parameters:
1757 * objin
1758 * Pointer to the LutMap to be copied.
1759 * objout
1760 * Pointer to the LutMap being constructed.
1761 * status
1762 * Pointer to the inherited status variable.
1763 */
1764
1765 /* Local Variables: */
1766 AstLutMap *out; /* Pointer to output LutMap */
1767 AstLutMap *in; /* Pointer to input LutMap */
1768
1769 /* Check the global error status. */
1770 if ( !astOK ) return;
1771
1772 /* Obtain a pointer to the input and output LutMaps. */
1773 in= (AstLutMap *) objin;
1774 out = (AstLutMap *) objout;
1775
1776 /* Nullify all output pointers. */
1777 out->lut = NULL;
1778 out->luti = NULL;
1779 out->flagsi = NULL;
1780 out->indexi = NULL;
1781
1782 /* Allocate memory and store a copy of the lookup table data. */
1783 out->lut = astStore( NULL, in->lut,
1784 sizeof( double ) * (size_t) in->nlut );
1785
1786 /* Do the arrays used for the inverse transformation, if they exist. */
1787 if( in->luti ) out->luti = astStore( NULL, in->luti,
1788 sizeof( double ) * (size_t) in->nluti );
1789 if( in->flagsi ) out->flagsi = astStore( NULL, in->flagsi,
1790 sizeof( double ) * (size_t) in->nluti );
1791 if( in->indexi ) out->indexi = astStore( NULL, in->indexi,
1792 sizeof( double ) * (size_t) in->nluti );
1793 }
1794
1795 /* Destructor. */
1796 /* ----------- */
Delete(AstObject * obj,int * status)1797 static void Delete( AstObject *obj, int *status ) {
1798 /*
1799 * Name:
1800 * Delete
1801
1802 * Purpose:
1803 * Destructor for LutMap objects.
1804
1805 * Type:
1806 * Private function.
1807
1808 * Synopsis:
1809 * void Delete( AstObject *obj, int *status )
1810
1811 * Description:
1812 * This function implements the destructor for LutMap objects.
1813
1814 * Parameters:
1815 * obj
1816 * Pointer to the LutMap to be deleted.
1817 * status
1818 * Pointer to the inherited status variable.
1819 */
1820
1821 /* Local Variables: */
1822 AstLutMap *this; /* Pointer to LutMap */
1823
1824 /* Obtain a pointer to the LutMap structure. */
1825 this = (AstLutMap *) obj;
1826
1827 /* Free the memory holding the lookup tables, etc. */
1828 this->lut = astFree( this->lut );
1829 this->luti = astFree( this->luti );
1830 this->flagsi = astFree( this->flagsi );
1831 this->indexi = astFree( this->indexi );
1832 }
1833
1834 /* Dump function. */
1835 /* -------------- */
Dump(AstObject * this_object,AstChannel * channel,int * status)1836 static void Dump( AstObject *this_object, AstChannel *channel, int *status ) {
1837 /*
1838 * Name:
1839 * Dump
1840
1841 * Purpose:
1842 * Dump function for LutMap objects.
1843
1844 * Type:
1845 * Private function.
1846
1847 * Synopsis:
1848 * void Dump( AstObject *this, AstChannel *channel, int *status )
1849
1850 * Description:
1851 * This function implements the Dump function which writes out data
1852 * for the LutMap class to an output Channel.
1853
1854 * Parameters:
1855 * this
1856 * Pointer to the LutMap whose data are being written.
1857 * channel
1858 * Pointer to the Channel to which the data are being written.
1859 * status
1860 * Pointer to the inherited status variable.
1861 */
1862
1863 /* Local Constants: */
1864 #define KEY_LEN 50 /* Maximum length of a keyword */
1865
1866 /* Local Variables: */
1867 AstLutMap *this; /* Pointer to the LutMap structure */
1868 char buff[ KEY_LEN + 1 ]; /* Buffer for keyword string */
1869 int ilut; /* Loop counter for table elements */
1870 int ival; /* Integer value */
1871 int set; /* Attribute value set? */
1872
1873 /* Check the global error status. */
1874 if ( !astOK ) return;
1875
1876 /* Obtain a pointer to the LutMap structure. */
1877 this = (AstLutMap *) this_object;
1878
1879 /* Write out values representing the instance variables for the LutMap
1880 class. Accompany these with appropriate comment strings, possibly
1881 depending on the values being written. */
1882
1883 /* Number of lookup table elements. */
1884 astWriteInt( channel, "Nlut", 1, 1, this->nlut,
1885 "Number of lookup table elements" );
1886
1887 /* Input coordinate at first element centre. */
1888 astWriteDouble( channel, "Start", ( this->start != 0.0 ), 1, this->start,
1889 "Input value at first element" );
1890
1891 /* Element spacing. */
1892 astWriteDouble( channel, "Incr", ( this->inc != 1.0 ), 1, this->inc,
1893 "Input value increment between elements" );
1894
1895 /* Interpolation method */
1896 set = TestLutInterp( this, status );
1897 ival = set ? GetLutInterp( this, status ) : astGetLutInterp( this );
1898 astWriteInt( channel, "LutInt", set, 1, ival, "Interpolation method" );
1899
1900 /* Lookup table contents. */
1901 for ( ilut = 0; ilut < this->nlut; ilut++ ) {
1902 if( this->lut[ ilut ] != AST__BAD ) {
1903 (void) sprintf( buff, "L%d", ilut + 1 );
1904 astWriteDouble( channel, buff, 1, 1, this->lut[ ilut ],
1905 ilut ? "" : "Lookup table elements..." );
1906 }
1907 }
1908
1909 /* Undefine macros local to this function. */
1910 #undef KEY_LEN
1911 }
1912
1913 /* Standard class functions. */
1914 /* ========================= */
1915 /* Implement the astIsALutMap and astCheckLutMap functions using the
1916 macros defined for this purpose in the "object.h" header file. */
astMAKE_ISA(LutMap,Mapping)1917 astMAKE_ISA(LutMap,Mapping)
1918 astMAKE_CHECK(LutMap)
1919
1920 AstLutMap *astLutMap_( int nlut, const double lut[],
1921 double start, double inc,
1922 const char *options, int *status, ...) {
1923 /*
1924 *++
1925 * Name:
1926 c astLutMap
1927 f AST_LUTMAP
1928
1929 * Purpose:
1930 * Create a LutMap.
1931
1932 * Type:
1933 * Public function.
1934
1935 * Synopsis:
1936 c #include "lutmap.h"
1937 c AstLutMap *astLutMap( int nlut, const double lut[],
1938 c double start, double inc,
1939 c const char *options, ... )
1940 f RESULT = AST_LUTMAP( NLUT, LUT, START, INC, OPTIONS, STATUS )
1941
1942 * Class Membership:
1943 * LutMap constructor.
1944
1945 * Description:
1946 * This function creates a new LutMap and optionally initialises
1947 * its attributes.
1948 *
1949 * A LutMap is a specialised form of Mapping which transforms
1950 * 1-dimensional coordinates by using linear interpolation in a
1951 * lookup table. Each input coordinate value is first scaled to
1952 * give the index of an entry in the table by subtracting a
1953 * starting value (the input coordinate corresponding to the first
1954 * table entry) and dividing by an increment (the difference in
1955 * input coordinate value between adjacent table entries).
1956 *
1957 * The resulting index will usually contain a fractional part, so
1958 * the output coordinate value is then generated by interpolating
1959 * linearly between the appropriate entries in the table. If the
1960 * index lies outside the range of the table, linear extrapolation
1961 * is used based on the two nearest entries (i.e. the two entries
1962 * at the start or end of the table, as appropriate).
1963 *
1964 * If the lookup table entries increase or decrease monotonically,
1965 * then the inverse transformation may also be performed.
1966
1967 * Parameters:
1968 c nlut
1969 f NLUT = INTEGER (Given)
1970 * The number of entries in the lookup table. This value must be
1971 * at least 2.
1972 c lut
1973 f LUT( NLUT ) = DOUBLE PRECISION (Given)
1974 c An array containing the "nlut"
1975 f An array containing the
1976 * lookup table entries.
1977 c start
1978 f START = DOUBLE PRECISION (Given)
1979 * The input coordinate value which corresponds to the first lookup
1980 * table entry.
1981 c inc
1982 f INC = DOUBLE PRECISION (Given)
1983 * The lookup table spacing (the increment in input coordinate
1984 * value between successive lookup table entries). This value
1985 * may be positive or negative, but must not be zero.
1986 c options
1987 f OPTIONS = CHARACTER * ( * ) (Given)
1988 c Pointer to a null-terminated string containing an optional
1989 c comma-separated list of attribute assignments to be used for
1990 c initialising the new LutMap. The syntax used is identical to
1991 c that for the astSet function and may include "printf" format
1992 c specifiers identified by "%" symbols in the normal way.
1993 f A character string containing an optional comma-separated
1994 f list of attribute assignments to be used for initialising the
1995 f new LutMap. The syntax used is identical to that for the
1996 f AST_SET routine.
1997 c ...
1998 c If the "options" string contains "%" format specifiers, then
1999 c an optional list of additional arguments may follow it in
2000 c order to supply values to be substituted for these
2001 c specifiers. The rules for supplying these are identical to
2002 c those for the astSet function (and for the C "printf"
2003 c function).
2004 f STATUS = INTEGER (Given and Returned)
2005 f The global status.
2006
2007 * Returned Value:
2008 c astLutMap()
2009 f AST_LUTMAP = INTEGER
2010 * A pointer to the new LutMap.
2011
2012 * Notes:
2013 * - If the entries in the lookup table either increase or decrease
2014 * monotonically, then the new LutMap's TranInverse attribute will
2015 * have a value of one, indicating that the inverse transformation
2016 * can be performed. Otherwise, it will have a value of zero, so
2017 * that any attempt to use the inverse transformation will result
2018 * in an error.
2019 * - A null Object pointer (AST__NULL) will be returned if this
2020 c function is invoked with the AST error status set, or if it
2021 f function is invoked with STATUS set to an error value, or if it
2022 * should fail for any reason.
2023
2024 * Status Handling:
2025 * The protected interface to this function includes an extra
2026 * parameter at the end of the parameter list descirbed above. This
2027 * parameter is a pointer to the integer inherited status
2028 * variable: "int *status".
2029
2030 *--
2031 */
2032
2033 /* Local Variables: */
2034 astDECLARE_GLOBALS /* Pointer to thread-specific global data */
2035 AstLutMap *new; /* Pointer to new LutMap */
2036 va_list args; /* Variable argument list */
2037
2038 /* Get a pointer to the thread specific global data structure. */
2039 astGET_GLOBALS(NULL);
2040
2041 /* Check the global status. */
2042 if ( !astOK ) return NULL;
2043
2044 /* Initialise the LutMap, allocating memory and initialising the
2045 virtual function table as well if necessary. */
2046 new = astInitLutMap( NULL, sizeof( AstLutMap ), !class_init, &class_vtab,
2047 "LutMap", nlut, lut, start, inc );
2048
2049 /* If successful, note that the virtual function table has been
2050 initialised. */
2051 if ( astOK ) {
2052 class_init = 1;
2053
2054 /* Obtain the variable argument list and pass it along with the
2055 options string to the astVSet method to initialise the new
2056 LutMap's attributes. */
2057 va_start( args, status );
2058 astVSet( new, options, NULL, args );
2059 va_end( args );
2060
2061 /* If an error occurred, clean up by deleting the new object. */
2062 if ( !astOK ) new = astDelete( new );
2063 }
2064
2065 /* Return a pointer to the new LutMap. */
2066 return new;
2067 }
2068
astLutMapId_(int nlut,const double lut[],double start,double inc,const char * options,...)2069 AstLutMap *astLutMapId_( int nlut, const double lut[],
2070 double start, double inc,
2071 const char *options, ... ) {
2072 /*
2073 * Name:
2074 * astLutMapId_
2075
2076 * Purpose:
2077 * Create a LutMap.
2078
2079 * Type:
2080 * Private function.
2081
2082 * Synopsis:
2083 * #include "lutmap.h"
2084 * AstLutMap *astLutMapId( int nlut, const double lut[],
2085 * double start, double inc,
2086 * const char *options, ... )
2087
2088 * Class Membership:
2089 * LutMap constructor.
2090
2091 * Description:
2092 * This function implements the external (public) interface to the
2093 * astLutMap constructor function. It returns an ID value (instead
2094 * of a true C pointer) to external users, and must be provided
2095 * because astLutMap_ has a variable argument list which cannot be
2096 * encapsulated in a macro (where this conversion would otherwise
2097 * occur).
2098 *
2099 * The variable argument list also prevents this function from
2100 * invoking astLutMap_ directly, so it must be a re-implementation
2101 * of it in all respects, except for the final conversion of the
2102 * result to an ID value.
2103
2104 * Parameters:
2105 * As for astLutMap_.
2106
2107 * Returned Value:
2108 * The ID value associated with the new LutMap.
2109 */
2110
2111 /* Local Variables: */
2112 astDECLARE_GLOBALS /* Pointer to thread-specific global data */
2113 AstLutMap *new; /* Pointer to new LutMap */
2114 va_list args; /* Variable argument list */
2115
2116 int *status; /* Pointer to inherited status value */
2117
2118 /* Get a pointer to the inherited status value. */
2119 status = astGetStatusPtr;
2120
2121 /* Get a pointer to the thread specific global data structure. */
2122 astGET_GLOBALS(NULL);
2123
2124 /* Check the global status. */
2125 if ( !astOK ) return NULL;
2126
2127 /* Initialise the LutMap, allocating memory and initialising the
2128 virtual function table as well if necessary. */
2129 new = astInitLutMap( NULL, sizeof( AstLutMap ), !class_init, &class_vtab,
2130 "LutMap", nlut, lut, start, inc );
2131
2132 /* If successful, note that the virtual function table has been
2133 initialised. */
2134 if ( astOK ) {
2135 class_init = 1;
2136
2137 /* Obtain the variable argument list and pass it along with the
2138 options string to the astVSet method to initialise the new LutMap's
2139 attributes. */
2140 va_start( args, options );
2141 astVSet( new, options, NULL, args );
2142 va_end( args );
2143
2144 /* If an error occurred, clean up by deleting the new object. */
2145 if ( !astOK ) new = astDelete( new );
2146 }
2147
2148 /* Return an ID value for the new LutMap. */
2149 return astMakeId( new );
2150 }
2151
astInitLutMap_(void * mem,size_t size,int init,AstLutMapVtab * vtab,const char * name,int nlut,const double lut[],double start,double inc,int * status)2152 AstLutMap *astInitLutMap_( void *mem, size_t size, int init,
2153 AstLutMapVtab *vtab, const char *name,
2154 int nlut, const double lut[],
2155 double start, double inc, int *status ) {
2156 /*
2157 *+
2158 * Name:
2159 * astInitLutMap
2160
2161 * Purpose:
2162 * Initialise a LutMap.
2163
2164 * Type:
2165 * Protected function.
2166
2167 * Synopsis:
2168 * #include "lutmap.h"
2169 * AstLutMap *astInitLutMap( void *mem, size_t size, int init,
2170 * AstLutMapVtab *vtab, const char *name,
2171 * int nlut, const double lut[],
2172 * double start, double inc )
2173
2174 * Class Membership:
2175 * LutMap initialiser.
2176
2177 * Description:
2178 * This function is provided for use by class implementations to
2179 * initialise a new LutMap object. It allocates memory (if
2180 * necessary) to accommodate the LutMap plus any additional data
2181 * associated with the derived class. It then initialises a LutMap
2182 * structure at the start of this memory. If the "init" flag is
2183 * set, it also initialises the contents of a virtual function
2184 * table for a LutMap at the start of the memory passed via the
2185 * "vtab" parameter.
2186
2187 * Parameters:
2188 * mem
2189 * A pointer to the memory in which the LutMap is to be
2190 * initialised. This must be of sufficient size to accommodate
2191 * the LutMap data (sizeof(LutMap)) plus any data used by the
2192 * derived class. If a value of NULL is given, this function
2193 * will allocate the memory itself using the "size" parameter to
2194 * determine its size.
2195 * size
2196 * The amount of memory used by the LutMap (plus derived class
2197 * data). This will be used to allocate memory if a value of
2198 * NULL is given for the "mem" parameter. This value is also
2199 * stored in the LutMap structure, so a valid value must be
2200 * supplied even if not required for allocating memory.
2201 * init
2202 * A logical flag indicating if the LutMap's virtual function
2203 * table is to be initialised. If this value is non-zero, the
2204 * virtual function table will be initialised by this function.
2205 * vtab
2206 * Pointer to the start of the virtual function table to be
2207 * associated with the new LutMap.
2208 * name
2209 * Pointer to a constant null-terminated character string which
2210 * contains the name of the class to which the new object
2211 * belongs (it is this pointer value that will subsequently be
2212 * returned by the astGetClass method).
2213 * nlut
2214 * The number of elements in the lookup table. This value must
2215 * be at least 2.
2216 * lut
2217 * An array containing the "nlut" lookup table elements.
2218 * start
2219 * The input coordinate value which corresponds with the first
2220 * lookup table element.
2221 * inc
2222 * The lookup table element spacing (i.e. the increment in input
2223 * coordinate value between successive lookup table elements).
2224
2225 * Returned Value:
2226 * A pointer to the new LutMap.
2227
2228 * Notes:
2229 * - A null pointer will be returned if this function is invoked
2230 * with the global error status set, or if it should fail for any
2231 * reason.
2232 *-
2233 */
2234
2235 /* Local Variables: */
2236 AstLutMap *new; /* Pointer to new LutMap */
2237 double *luti; /* Pointer to table for inverse transformation */
2238 double *p; /* Pointer to next lut element */
2239 int *flagsi; /* Pointer to flags for inverse transformation */
2240 int *indexi; /* Pointer to translation from original to reduced */
2241 int dirn; /* +1 => values increasing, -1 => values decreasing */
2242 int ilut; /* Loop counter for LUT elements */
2243 int nluti; /* Length of "luti" array */
2244
2245 /* Initialise. */
2246 new = NULL;
2247
2248 /* Check the global status. */
2249 if ( !astOK ) return new;
2250
2251 /* If necessary, initialise the virtual function table. */
2252 if ( init ) astInitLutMapVtab( vtab, name );
2253
2254 /* Check that the number of lookup table elements is valid. */
2255 if ( nlut < 2 ) {
2256 astError( AST__LUTIN, "astInitLutMap(%s): Invalid number of lookup "
2257 "table elements (%d).", status, name, nlut );
2258 astError( AST__LUTIN, "This value should be at least 2." , status);
2259
2260 /* Also check that the input value increment is not zero. */
2261 } else if ( inc == 0.0 ) {
2262 astError( AST__LUTII, "astInitLutMap(%s): An input value increment of "
2263 "zero between lookup table elements is not allowed.", status, name );
2264
2265 /* Determine if the element values increase or decrease monotonically (except
2266 that adjacent entries can be equal). We can only implement the inverse
2267 transformation if this is so. The inverse transformation will generate
2268 AST__BAD output values for sections of the table that contain equal
2269 adjacent values, or hold AST__BAD values. */
2270 } else {
2271 dirn = GetMonotonic( nlut, lut, &nluti, &luti, &flagsi, &indexi,
2272 status );
2273
2274 /* Initialise a Mapping structure (the parent class) as the first
2275 component within the LutMap structure, allocating memory if
2276 necessary. Specify that the Mapping should be defined in the
2277 forward direction, and conditionally in the inverse direction. */
2278 new = (AstLutMap *) astInitMapping( mem, size, 0,
2279 (AstMappingVtab *) vtab, name,
2280 1, 1, 1, ( dirn != 0 ) );
2281
2282 if ( astOK ) {
2283
2284 /* Initialise the LutMap data. */
2285 /* ---------------------------- */
2286 new->nlut = nlut;
2287 new->start = start;
2288 new->inc = inc;
2289 new->lutinterp = LINEAR;
2290 new->nluti = nluti;
2291 new->luti = luti;
2292 new->flagsi = flagsi;
2293 new->indexi = indexi;
2294
2295 /* Allocate memory and store the lookup table. */
2296 new->lut = astStore( NULL, lut, sizeof( double ) * (size_t) nlut );
2297
2298 /* Replace an NaN values by AST__BAD */
2299 p = new->lut;
2300 for ( ilut = 0; ilut < nlut; ilut++, p++ ) {
2301 if( !astISFINITE(*p) ) *p = AST__BAD;
2302 }
2303
2304 /* Initialise the retained input and output coordinate values. */
2305 new->last_fwd_in = AST__BAD;
2306 new->last_fwd_out = AST__BAD;
2307 new->last_inv_in = AST__BAD;
2308 new->last_inv_out = AST__BAD;
2309 }
2310
2311 /* If an error occurred, clean up by deleting the new LutMap. */
2312 if ( !astOK ) new = astDelete( new );
2313 }
2314
2315 /* Return a pointer to the new LutMap. */
2316 return new;
2317 }
2318
astLoadLutMap_(void * mem,size_t size,AstLutMapVtab * vtab,const char * name,AstChannel * channel,int * status)2319 AstLutMap *astLoadLutMap_( void *mem, size_t size,
2320 AstLutMapVtab *vtab, const char *name,
2321 AstChannel *channel, int *status ) {
2322 /*
2323 *+
2324 * Name:
2325 * astLoadLutMap
2326
2327 * Purpose:
2328 * Load a LutMap.
2329
2330 * Type:
2331 * Protected function.
2332
2333 * Synopsis:
2334 * #include "lutmap.h"
2335 * AstLutMap *astLoadLutMap( void *mem, size_t size,
2336 * AstLutMapVtab *vtab, const char *name,
2337 * AstChannel *channel )
2338
2339 * Class Membership:
2340 * LutMap loader.
2341
2342 * Description:
2343 * This function is provided to load a new LutMap using data read
2344 * from a Channel. It first loads the data used by the parent class
2345 * (which allocates memory if necessary) and then initialises a
2346 * LutMap structure in this memory, using data read from the input
2347 * Channel.
2348 *
2349 * If the "init" flag is set, it also initialises the contents of a
2350 * virtual function table for a LutMap at the start of the memory
2351 * passed via the "vtab" parameter.
2352
2353
2354 * Parameters:
2355 * mem
2356 * A pointer to the memory into which the LutMap is to be
2357 * loaded. This must be of sufficient size to accommodate the
2358 * LutMap data (sizeof(LutMap)) plus any data used by derived
2359 * classes. If a value of NULL is given, this function will
2360 * allocate the memory itself using the "size" parameter to
2361 * determine its size.
2362 * size
2363 * The amount of memory used by the LutMap (plus derived class
2364 * data). This will be used to allocate memory if a value of
2365 * NULL is given for the "mem" parameter. This value is also
2366 * stored in the LutMap structure, so a valid value must be
2367 * supplied even if not required for allocating memory.
2368 *
2369 * If the "vtab" parameter is NULL, the "size" value is ignored
2370 * and sizeof(AstLutMap) is used instead.
2371 * vtab
2372 * Pointer to the start of the virtual function table to be
2373 * associated with the new LutMap. If this is NULL, a pointer
2374 * to the (static) virtual function table for the LutMap class
2375 * is used instead.
2376 * name
2377 * Pointer to a constant null-terminated character string which
2378 * contains the name of the class to which the new object
2379 * belongs (it is this pointer value that will subsequently be
2380 * returned by the astGetClass method).
2381 *
2382 * If the "vtab" parameter is NULL, the "name" value is ignored
2383 * and a pointer to the string "LutMap" is used instead.
2384
2385 * Returned Value:
2386 * A pointer to the new LutMap.
2387
2388 * Notes:
2389 * - A null pointer will be returned if this function is invoked
2390 * with the global error status set, or if it should fail for any
2391 * reason.
2392 *-
2393 */
2394
2395 /* Local Constants: */
2396 astDECLARE_GLOBALS /* Pointer to thread-specific global data */
2397 #define KEY_LEN 50 /* Maximum length of a keyword */
2398
2399 /* Local Variables: */
2400 AstLutMap *new; /* Pointer to the new LutMap */
2401 char buff[ KEY_LEN + 1 ]; /* Buffer for keyword string */
2402 int ilut; /* Get a pointer to the thread specific global data structure. */
2403 astGET_GLOBALS(channel);
2404
2405 /* Loop counter for table elements */
2406
2407 /* Initialise. */
2408 new = NULL;
2409
2410 /* Check the global error status. */
2411 if ( !astOK ) return new;
2412
2413 /* If a NULL virtual function table has been supplied, then this is
2414 the first loader to be invoked for this LutMap. In this case the
2415 LutMap belongs to this class, so supply appropriate values to be
2416 passed to the parent class loader (and its parent, etc.). */
2417 if ( !vtab ) {
2418 size = sizeof( AstLutMap );
2419 vtab = &class_vtab;
2420 name = "LutMap";
2421
2422 /* If required, initialise the virtual function table for this class. */
2423 if ( !class_init ) {
2424 astInitLutMapVtab( vtab, name );
2425 class_init = 1;
2426 }
2427 }
2428
2429 /* Invoke the parent class loader to load data for all the ancestral
2430 classes of the current one, returning a pointer to the resulting
2431 partly-built LutMap. */
2432 new = astLoadMapping( mem, size, (AstMappingVtab *) vtab, name,
2433 channel );
2434
2435 if ( astOK ) {
2436
2437 /* Read input data. */
2438 /* ================ */
2439 /* Request the input Channel to read all the input data appropriate to
2440 this class into the internal "values list". */
2441 astReadClassData( channel, "LutMap" );
2442
2443 /* Now read each individual data item from this list and use it to
2444 initialise the appropriate instance variable(s) for this class. */
2445
2446 /* Number of lookup table elements. */
2447 new->nlut = astReadInt( channel, "nlut", 2 );
2448
2449 /* Starting input coordinate value. */
2450 new->start = astReadDouble( channel, "start", 0.0 );
2451
2452 /* Input coordinate value increment. */
2453 new->inc = astReadDouble( channel, "incr", 1.0 );
2454
2455 /* Interpolation method */
2456 new->lutinterp = astReadInt( channel, "lutint", LINEAR );
2457 if ( TestLutInterp( new, status ) ) SetLutInterp( new, new->lutinterp, status );
2458
2459 /* Allocate memory to hold the lookup table elements. */
2460 new->lut = astMalloc( sizeof( double ) * (size_t) new->nlut );
2461
2462 /* If OK, loop to read each element. */
2463 if ( astOK ) {
2464 for ( ilut = 0; ilut < new->nlut; ilut++ ) {
2465 (void) sprintf( buff, "l%d", ilut + 1 );
2466 new->lut[ ilut ] = astReadDouble( channel, buff, AST__BAD );
2467 }
2468
2469 /* Initialise the retained input and output coordinate values. */
2470 new->last_fwd_in = AST__BAD;
2471 new->last_fwd_out = AST__BAD;
2472 new->last_inv_in = AST__BAD;
2473 new->last_inv_out = AST__BAD;
2474
2475 /* See if the array is monotonic increasing or decreasing. */
2476 (void) GetMonotonic( new->nlut, new->lut, &(new->nluti),
2477 &(new->luti), &(new->flagsi), &(new->indexi),
2478 status );
2479 }
2480 }
2481
2482 /* If an error occurred, clean up by deleting the new LutMap. */
2483 if ( !astOK ) new = astDelete( new );
2484
2485 /* Return the new LutMap pointer. */
2486 return new;
2487
2488 /* Undefine macros local to this function. */
2489 #undef KEY_LEN
2490 }
2491
2492 /* Virtual function interfaces. */
2493 /* ============================ */
2494 /* These provide the external interface to the virtual functions
2495 defined by this class. Each simply checks the global error status
2496 and then locates and executes the appropriate member function,
2497 using the function pointer stored in the object's virtual function
2498 table (this pointer is located using the astMEMBER macro defined in
2499 "object.h").
2500
2501 Note that the member function may not be the one defined here, as
2502 it may have been over-ridden by a derived class. However, it should
2503 still have the same interface. */
2504
astGetLutMapInfo_(AstLutMap * this,double * start,double * inc,int * nlut,int * status)2505 double *astGetLutMapInfo_( AstLutMap *this, double *start, double *inc,
2506 int *nlut, int *status ){
2507 if( !astOK ) return NULL;
2508 return (**astMEMBER(this,LutMap,GetLutMapInfo))( this, start, inc, nlut,
2509 status );
2510 }
2511
2512
2513
2514