1 /*
2 *  Name:
3 *     pointset.c
4 
5 *  Purpose:
6 *     Implement the PointSet class.
7 
8 *  Description:
9 *     This file implements the PointSet class. For a description of
10 *     the class and its interface, see the .h file of the same name.
11 
12 *  Inheritance:
13 *     The PointSet class inherits from the Object class.
14 
15 *  Copyright:
16 *     Copyright (C) 1997-2006 Council for the Central Laboratory of the
17 *     Research Councils
18 
19 *  Licence:
20 *     This program is free software: you can redistribute it and/or
21 *     modify it under the terms of the GNU Lesser General Public
22 *     License as published by the Free Software Foundation, either
23 *     version 3 of the License, or (at your option) any later
24 *     version.
25 *
26 *     This program is distributed in the hope that it will be useful,
27 *     but WITHOUT ANY WARRANTY; without even the implied warranty of
28 *     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
29 *     GNU Lesser General Public License for more details.
30 *
31 *     You should have received a copy of the GNU Lesser General
32 *     License along with this program.  If not, see
33 *     <http://www.gnu.org/licenses/>.
34 
35 *  Authors:
36 *     RFWS: R.F. Warren-Smith (Starlink)
37 
38 *  History:
39 *     1-FEB-1996 (RFWS):
40 *        Original version.
41 *     27-SEP-1996 (RFWS):
42 *        Added external interface and I/O facilities.
43 *     8-JAN-2003 (DSB):
44 *        Changed private InitVtab method to protected astInitPointSetVtab
45 *        method.
46 *     9-SEP-2004 (DSB):
47 *        Added astPermPoints.
48 *     5-OCT-2004 (DSB):
49 *        Bug fix in astLoadPointSet - npoint was used as size for new array
50 *        of pointers (changed to ncoord).
51 *     19-OCT-2004 (DSB):
52 *        Added astSetNpoint.
53 *     2-NOV-2004 (DSB):
54 *        - Do not write out AST__BAD axis values in the Dump function.
55 *        - Override astEqual method.
56 *        - Add protected PointAccuracy attribute.
57 *     7-JAN-2005 (DSB):
58 *        Added astAppendPoints.
59 *     14-FEB-2006 (DSB):
60 *        Override astGetObjSize.
61 *     22-FEB-2006 (DSB):
62 *        Avoid allocating memory for "acc" unless needed.
63 *     1-MAY-2009 (DSB):
64 *        Added astBndPoints.
65 *     16-JUN-2011 (DSB):
66 *        Added astReplaceNan.
67 *     2-OCT-2012 (DSB):
68 *        Check for Infs as well as NaNs.
69 */
70 
71 /* Module Macros. */
72 /* ============== */
73 /* Set the name of the class we are implementing. This indicates to
74    the header files that define class interfaces that they should make
75    "protected" symbols available. */
76 #define astCLASS PointSet
77 
78 /* Values that control the behaviour of the astReplaceNaN method. */
79 #define IGNORE_NANS  0
80 #define REPLACE_NANS 1
81 #define REPORT_NANS  2
82 
83 /*
84 *
85 *  Name:
86 *     MAKE_CLEAR
87 
88 *  Purpose:
89 *     Implement a method to clear a single value in a multi-valued attribute.
90 
91 *  Type:
92 *     Private macro.
93 
94 *  Synopsis:
95 *     #include "pointset.h"
96 *     MAKE_CLEAR(attr,component,assign)
97 
98 *  Class Membership:
99 *     Defined by the PointSet class.
100 
101 *  Description:
102 *     This macro expands to an implementation of a private member function of
103 *     the form:
104 *
105 *        static void Clear<Attribute>( AstPointSet *this, int axis )
106 *
107 *     and an external interface function of the form:
108 *
109 *        void astClear<Attribute>_( AstPointSet *this, int axis )
110 *
111 *     which implement a method for clearing a single value in a specified
112 *     multi-valued attribute for an axis of a PointSet. The "axis" value
113 *     must be in the range 0 to (ncoord-1).
114 
115 *  Parameters:
116 *     attr
117 *        The name of the attribute to be cleared, as it appears in the function
118 *        name (e.g. PointAccuracy in "astClearPointAccuracy").
119 *     component
120 *        The name of the class structure component that holds the attribute
121 *        value.
122 *     assign
123 *        An expression that evaluates to the value to assign to the component
124 *        to clear its value.
125 
126 *  Notes:
127 *     -  To avoid problems with some compilers, you should not leave any white
128 *     space around the macro arguments.
129 *
130 */
131 
132 /* Define the macro. */
133 #define MAKE_CLEAR(attr,component,assign) \
134 \
135 /* Private member function. */ \
136 /* ------------------------ */ \
137 static void Clear##attr( AstPointSet *this, int axis, int *status ) { \
138 \
139 /* Check the global error status. */ \
140    if ( !astOK ) return; \
141 \
142 /* Validate the axis index. */ \
143    if( axis < 0 || axis >= this->ncoord ){ \
144       astError( AST__AXIIN, "%s(%s): Index (%d) is invalid for attribute " \
145                 #attr " - it should be in the range 1 to %d.", status, \
146                 "astClear" #attr, astGetClass( this ), \
147                 axis + 1, this->ncoord ); \
148 \
149 /* Assign the "clear" value. */ \
150    } else if( this->component ){ \
151       this->component[ axis ] = (assign); \
152    } \
153 } \
154 \
155 /* External interface. */ \
156 /* ------------------- */ \
157 void astClear##attr##_( AstPointSet *this, int axis, int *status ) { \
158 \
159 /* Check the global error status. */ \
160    if ( !astOK ) return; \
161 \
162 /* Invoke the required method via the virtual function table. */ \
163    (**astMEMBER(this,PointSet,Clear##attr))( this, axis, status ); \
164 }
165 
166 
167 /*
168 *
169 *  Name:
170 *     MAKE_GET
171 
172 *  Purpose:
173 *     Implement a method to get a single value in a multi-valued attribute.
174 
175 *  Type:
176 *     Private macro.
177 
178 *  Synopsis:
179 *     #include "pointset.h"
180 *     MAKE_GET(attr,type,bad_value,assign)
181 
182 *  Class Membership:
183 *     Defined by the PointSet class.
184 
185 *  Description:
186 *     This macro expands to an implementation of a private member function of
187 *     the form:
188 *
189 *        static <Type> Get<Attribute>( AstPointSet *this, int axis )
190 *
191 *     and an external interface function of the form:
192 *
193 *        <Type> astGet<Attribute>_( AstPointSet *this, int axis )
194 *
195 *     which implement a method for getting a single value from a specified
196 *     multi-valued attribute for an axis of a PointSet. The "axis" value
197 *     must be in the range 0 to (ncoord-1).
198 
199 
200 *  Parameters:
201 *     attr
202 *        The name of the attribute whose value is to be obtained, as it
203 *        appears in the function name (e.g. PointAccuracy in "astGetPointAccuracy").
204 *     type
205 *        The C type of the attribute.
206 *     bad_value
207 *        A constant value to return if the global error status is set, or if
208 *        the function fails.
209 *     assign
210 *        An expression that evaluates to the value to be returned. This can
211 *        use the string "axis" to represent the zero-based value index.
212 
213 *  Notes:
214 *     -  To avoid problems with some compilers, you should not leave any white
215 *     space around the macro arguments.
216 *
217 */
218 
219 /* Define the macro. */
220 #define MAKE_GET(attr,type,bad_value,assign) \
221 \
222 /* Private member function. */ \
223 /* ------------------------ */ \
224 static type Get##attr( AstPointSet *this, int axis, int *status ) { \
225    type result;                  /* Result to be returned */ \
226 \
227 /* Initialise */ \
228    result = bad_value; \
229 \
230 /* Check the global error status. */ \
231    if ( !astOK ) return result; \
232 \
233 /* Validate the axis index. */ \
234    if( axis < 0 || axis >= this->ncoord ){ \
235       astError( AST__AXIIN, "%s(%s): Index (%d) is invalid for attribute " \
236                 #attr " - it should be in the range 1 to %d.", status, \
237                 "astGet" #attr, astGetClass( this ), \
238                 axis + 1, this->ncoord ); \
239 \
240 /* Assign the result value. */ \
241    } else { \
242       result = (assign); \
243    } \
244 \
245 /* Check for errors and clear the result if necessary. */ \
246    if ( !astOK ) result = (bad_value); \
247 \
248 /* Return the result. */ \
249    return result; \
250 } \
251 /* External interface. */ \
252 /* ------------------- */  \
253 type astGet##attr##_( AstPointSet *this, int axis, int *status ) { \
254 \
255 /* Check the global error status. */ \
256    if ( !astOK ) return (bad_value); \
257 \
258 /* Invoke the required method via the virtual function table. */ \
259    return (**astMEMBER(this,PointSet,Get##attr))( this, axis, status ); \
260 }
261 
262 /*
263 *
264 *  Name:
265 *     MAKE_SET
266 
267 *  Purpose:
268 *     Implement a method to set a single value in a multi-valued attribute
269 *     for a PointSet.
270 
271 *  Type:
272 *     Private macro.
273 
274 *  Synopsis:
275 *     #include "pointset.h"
276 *     MAKE_SET(attr,type,component,assign,null)
277 
278 *  Class Membership:
279 *     Defined by the PointSet class.
280 
281 *  Description:
282 *     This macro expands to an implementation of a private member function of
283 *     the form:
284 *
285 *        static void Set<Attribute>( AstPointSet *this, int axis, <Type> value )
286 *
287 *     and an external interface function of the form:
288 *
289 *        void astSet<Attribute>_( AstPointSet *this, int axis, <Type> value )
290 *
291 *     which implement a method for setting a single value in a specified
292 *     multi-valued attribute for a PointSet. The "axis" value must be in
293 *     the range 0 to (ncoord-1).
294 
295 *  Parameters:
296 *      attr
297 *         The name of the attribute to be set, as it appears in the function
298 *         name (e.g. PointAccuracy in "astSetPointAccuracy").
299 *      type
300 *         The C type of the attribute.
301 *      component
302 *         The name of the class structure component that holds the attribute
303 *         value.
304 *      assign
305 *         An expression that evaluates to the value to be assigned to the
306 *         component.
307 *      null
308 *         The value to initialise newly created array elements to.
309 
310 *  Notes:
311 *     -  To avoid problems with some compilers, you should not leave any white
312 *     space around the macro arguments.
313 *-
314 */
315 
316 /* Define the macro. */
317 #define MAKE_SET(attr,type,component,assign,null) \
318 \
319 /* Private member function. */ \
320 /* ------------------------ */ \
321 static void Set##attr( AstPointSet *this, int axis, type value, int *status ) { \
322    int i; \
323 \
324 /* Check the global error status. */ \
325    if ( !astOK ) return; \
326 \
327 /* Validate the axis index. */ \
328    if( axis < 0 || axis >= this->ncoord ){ \
329       astError( AST__AXIIN, "%s(%s): Index (%d) is invalid for attribute " \
330                 #attr " - it should be in the range 1 to %d.", status, \
331                 "astSet" #attr, astGetClass( this ), \
332                 axis + 1, this->ncoord ); \
333 \
334 /* Store the new value in the structure component. */ \
335    } else { \
336       if( !this->component ){ \
337          this->component = astMalloc( this->ncoord*sizeof( type ) ); \
338          for( i = 0; i < this->ncoord; i++ ) { \
339             this->component[ i ] = null; \
340          } \
341       } \
342       this->component[ axis ] = (assign); \
343    } \
344 } \
345 \
346 /* External interface. */ \
347 /* ------------------- */ \
348 void astSet##attr##_( AstPointSet *this, int axis, type value, int *status ) { \
349 \
350 /* Check the global error status. */ \
351    if ( !astOK ) return; \
352 \
353 /* Invoke the required method via the virtual function table. */ \
354    (**astMEMBER(this,PointSet,Set##attr))( this, axis, value, status ); \
355 }
356 
357 /*
358 *
359 *  Name:
360 *     MAKE_TEST
361 
362 *  Purpose:
363 *     Implement a method to test if a single value has been set in a
364 *     multi-valued attribute for a class.
365 
366 *  Type:
367 *     Private macro.
368 
369 *  Synopsis:
370 *     #include "pointset.h"
371 *     MAKE_TEST(attr,assign)
372 
373 *  Class Membership:
374 *     Defined by the PointSet class.
375 
376 *  Description:
377 *     This macro expands to an implementation of a private member function of
378 *     the form:
379 *
380 *        static int Test<Attribute>( AstPointSet *this, int axis )
381 *
382 *     and an external interface function of the form:
383 *
384 *        int astTest<Attribute>_( AstPointSet *this, int axis )
385 *
386 *     which implement a method for testing if a single value in a specified
387 *     multi-valued attribute has been set for a class. The "axis" value
388 *     must be in the range 0 to (ncoord-1).
389 
390 *  Parameters:
391 *      attr
392 *         The name of the attribute to be tested, as it appears in the function
393 *         name (e.g. PointAccuracy in "astTestPointAccuracy").
394 *      assign
395 *         An expression that evaluates to 0 or 1, to be used as the returned
396 *         value. This can use the string "axis" to represent the zero-based
397 *         index of the value within the attribute.
398 
399 *  Notes:
400 *     -  To avoid problems with some compilers, you should not leave any white
401 *     space around the macro arguments.
402 *-
403 */
404 
405 /* Define the macro. */
406 #define MAKE_TEST(attr,assign) \
407 \
408 /* Private member function. */ \
409 /* ------------------------ */ \
410 static int Test##attr( AstPointSet *this, int axis, int *status ) { \
411    int result;                   /* Value to return */ \
412 \
413    result= 0; \
414 \
415 /* Check the global error status. */ \
416    if ( !astOK ) return result; \
417 \
418 /* Validate the axis index. */ \
419    if( axis < 0 || axis >= this->ncoord ){ \
420       astError( AST__AXIIN, "%s(%s): Index (%d) is invalid for attribute " \
421                 #attr " - it should be in the range 1 to %d.", status, \
422                 "astTest" #attr, astGetClass( this ), \
423                 axis + 1, this->ncoord ); \
424 \
425 /* Assign the result value. */ \
426    } else { \
427       result = (assign); \
428    } \
429 \
430 /* Check for errors and clear the result if necessary. */ \
431    if ( !astOK ) result = 0; \
432 \
433 /* Return the result. */ \
434    return result; \
435 } \
436 /* External interface. */ \
437 /* ------------------- */ \
438 int astTest##attr##_( AstPointSet *this, int axis, int *status ) { \
439 \
440 /* Check the global error status. */ \
441    if ( !astOK ) return 0; \
442 \
443 /* Invoke the required method via the virtual function table. */ \
444    return (**astMEMBER(this,PointSet,Test##attr))( this, axis, status ); \
445 }
446 
447 /* Include files. */
448 /* ============== */
449 /* Interface definitions. */
450 /* ---------------------- */
451 
452 #include "globals.h"             /* Thread-safe global data access */
453 #include "error.h"               /* Error reporting facilities */
454 #include "memory.h"              /* Memory allocation facilities */
455 #include "globals.h"             /* Thread-safe global data access */
456 #include "object.h"              /* Base Object class */
457 #include "pointset.h"            /* Interface definition for this class */
458 
459 /* Error code definitions. */
460 /* ----------------------- */
461 #include "ast_err.h"             /* AST error codes */
462 
463 /* C header files. */
464 /* --------------- */
465 #include <stdarg.h>
466 #include <stddef.h>
467 #include <stdlib.h>
468 #include <stdio.h>
469 #include <math.h>
470 #include <string.h>
471 
472 /* Module Variables. */
473 /* ================= */
474 
475 /* Address of this static variable is used as a unique identifier for
476    member of this class. */
477 static int class_check;
478 
479 /* This static variable is used to hold an IEEE 754 quiet double precision
480    Nan value. */
481 static double ast_nan;
482 
483 /* This static variable is used to hold an IEEE 754 quiet single precision
484    Nan value. */
485 static float ast_nanf;
486 
487 /* Enable or disable the astReplaceNan method. */
488 static int replace_nan = -1;
489 
490 /* Pointers to parent class methods which are extended by this class. */
491 static const char *(* parent_getattrib)( AstObject *, const char *, int * );
492 static int (* parent_testattrib)( AstObject *, const char *, int * );
493 static void (* parent_clearattrib)( AstObject *, const char *, int * );
494 static void (* parent_setattrib)( AstObject *, const char *, int * );
495 static int (* parent_equal)( AstObject *, AstObject *, int * );
496 static int (* parent_getobjsize)( AstObject *, int * );
497 
498 /* Define macros for accessing each item of thread specific global data. */
499 #ifdef THREAD_SAFE
500 
501 /* Define how to initialise thread-specific globals. */
502 #define GLOBAL_inits \
503    globals->Class_Init = 0; \
504    globals->GetAttrib_Buff[ 0 ] = 0;
505 
506 /* Create the function that initialises global data for this module. */
507 astMAKE_INITGLOBALS(PointSet)
508 
509 /* Define macros for accessing each item of thread specific global data. */
510 #define class_init astGLOBAL(PointSet,Class_Init)
511 #define class_vtab astGLOBAL(PointSet,Class_Vtab)
512 #define getattrib_buff astGLOBAL(PointSet,GetAttrib_Buff)
513 
514 static pthread_mutex_t mutex1 = PTHREAD_MUTEX_INITIALIZER;
515 #define LOCK_MUTEX1 pthread_mutex_lock( &mutex1 );
516 #define UNLOCK_MUTEX1 pthread_mutex_unlock( &mutex1 );
517 
518 /* If thread safety is not needed, declare and initialise globals at static
519    variables. */
520 #else
521 
522 static char getattrib_buff[ 101 ];
523 
524 
525 /* Define the class virtual function table and its initialisation flag
526    as static variables. */
527 static AstPointSetVtab class_vtab;   /* Virtual function table */
528 static int class_init = 0;       /* Virtual function table initialised? */
529 
530 #define LOCK_MUTEX1
531 #define UNLOCK_MUTEX1
532 
533 #endif
534 
535 /* External Interface Function Prototypes. */
536 /* ======================================= */
537 /* The following functions have public prototypes only (i.e. no
538    protected prototypes), so we must provide local prototypes for use
539    within this module. */
540 AstPointSet *astPointSetId_( int, int, const char *, int *, ...);
541 
542 /* Prototypes for Private Member Functions. */
543 /* ======================================== */
544 static const char *GetAttrib( AstObject *, const char *, int * );
545 static double **GetPoints( AstPointSet *, int * );
546 static int Equal( AstObject *, AstObject *, int * );
547 static int GetNcoord( const AstPointSet *, int * );
548 static int GetNpoint( const AstPointSet *, int * );
549 static int GetObjSize( AstObject *, int * );
550 static int ReplaceNaN( AstPointSet *, int * );
551 static int TestAttrib( AstObject *, const char *, int * );
552 static AstPointSet *AppendPoints( AstPointSet *, AstPointSet *, int * );
553 static void BndPoints( AstPointSet *, double *, double *, int * );
554 static void CheckPerm( AstPointSet *, const int *, const char *, int * );
555 static void ClearAttrib( AstObject *, const char *, int * );
556 static void Copy( const AstObject *, AstObject *, int * );
557 static void Delete( AstObject *, int * );
558 static void Dump( AstObject *, AstChannel *, int * );
559 static void PermPoints( AstPointSet *, int, const int[], int * );
560 static void SetAttrib( AstObject *, const char *, int * );
561 static void SetPoints( AstPointSet *, double **, int * );
562 static void SetNpoint( AstPointSet *, int, int * );
563 static void SetSubPoints( AstPointSet *, int, int, AstPointSet *, int * );
564 
565 static double GetPointAccuracy( AstPointSet *, int, int * );
566 static int TestPointAccuracy( AstPointSet *, int, int * );
567 static void ClearPointAccuracy( AstPointSet *, int, int * );
568 static void SetPointAccuracy( AstPointSet *, int, double, int * );
569 
570 /* Member functions. */
571 /* ================= */
AppendPoints(AstPointSet * this,AstPointSet * that,int * status)572 static AstPointSet *AppendPoints( AstPointSet *this, AstPointSet *that, int *status ) {
573 /*
574 *+
575 *  Name:
576 *     astAppendPoints
577 
578 *  Purpose:
579 *     Append one PointSet to another.
580 
581 *  Type:
582 *     Protected virtual function.
583 
584 *  Synopsis:
585 *     #include "pointset.h"
586 *     AstPointSet *astAppendPoints( AstPointSet *this, AstPointSet *that )
587 
588 *  Class Membership:
589 *     PointSet method.
590 
591 *  Description:
592 *     This function creates a new PointSet containing all the points in
593 *     "this" followed by all the points in "that".
594 
595 *  Parameters:
596 *     this
597 *        Pointer to the first PointSet.
598 *     that
599 *        Pointer to the second PointSet.
600 
601 *  Returned Value:
602 *     Pointer to the new PointSet.
603 
604 *  Notes:
605 *     - Axis accuracies are copied from "this".
606 *     - The Ncoord attribute of the two PointSets must match.
607 *     - NULL will be returned if an error has already occurred, or if this
608 *     function should fail for any reason.
609 *-
610 */
611 
612 /* Local Variables: */
613    AstPointSet *result;
614    double **ptr;
615    double **ptr1;
616    double **ptr2;
617    int ic;
618    int n1;
619    int n2;
620    int ncoord;
621    size_t nb2;
622    size_t nb1;
623 
624 /* Initialise */
625    result = NULL;
626 
627 /* Check the global error status. */
628    if ( !astOK ) return result;
629 
630 /* Check the two PointSets have the same Ncoord value. */
631    ncoord = astGetNcoord( this );
632    if( ncoord != astGetNcoord( that ) ) {
633       astError( AST__NPTIN, "astAppendPoints(%s): Number of coordinates "
634                 "per point differ in the two supplied PointSets.", status,
635                 astGetClass( this ) );
636 
637 /* Calculate the new size for the PointSet. */
638    } else {
639       n1 = astGetNpoint( this );
640       n2 = astGetNpoint( that );
641 
642 /* Create the new PointSet and get pointers to its data. */
643       result = astPointSet( n1 + n2, ncoord, "", status );
644       ptr1 = astGetPoints( this );
645       ptr2 = astGetPoints( that );
646       ptr = astGetPoints( result );
647       if( astOK ) {
648 
649 /* Copy the axis values for each coordinate in turn. */
650          nb1 = sizeof( double )*(size_t) n1;
651          nb2 = sizeof( double )*(size_t) n2;
652          for( ic = 0; ic < ncoord; ic++ ) {
653             memcpy( ptr[ ic ], ptr1[ ic ], nb1 );
654             memcpy( ptr[ ic ] + n1, ptr2[ ic ], nb2 );
655          }
656 
657 /* Copy any axis accuracies from "this". */
658          result->acc = this->acc ?
659                   astStore( NULL, this->acc, sizeof( double )*(size_t) ncoord )
660                   : NULL;
661       }
662    }
663 
664 /* Annul the result if an error occurred. */
665    if( !astOK ) result = astAnnul( result );
666 
667 /* Return the result. */
668    return result;
669 }
670 
BndPoints(AstPointSet * this,double * lbnd,double * ubnd,int * status)671 static void BndPoints( AstPointSet *this, double *lbnd, double *ubnd, int *status ) {
672 /*
673 *+
674 *  Name:
675 *     astBndPoints
676 
677 *  Purpose:
678 *     Find the axis bounds of the points in a PointSet.
679 
680 *  Type:
681 *     Protected virtual function.
682 
683 *  Synopsis:
684 *     #include "pointset.h"
685 *     void astBndPoints( AstPointSet *this, double *lbnd, double *ubnd )
686 
687 *  Class Membership:
688 *     PointSet method.
689 
690 *  Description:
691 *     This function returns the lower and upper limits of the axis values
692 *     of the points in a PointSet.
693 
694 *  Parameters:
695 *     this
696 *        Pointer to the first PointSet.
697 *     lbnd
698 *        Pointer to an array in which to return the lowest value for
699 *        each coordinate. The length of the array should equal the number
700 *        returned by astGetNcoord.
701 *     ubnd
702 *        Pointer to an array in which to return the highest value for
703 *        each coordinate. The length of the array should equal the number
704 *        returned by astGetNcoord.
705 
706 *-
707 */
708 
709 /* Local Variables: */
710    double **ptr;
711    double *p;
712    double lb;
713    double ub;
714    int ic;
715    int ip;
716    int nc;
717    int np;
718 
719 /* Check the global error status. */
720    if ( !astOK ) return;
721 
722 /* Get pointers to the PointSet data, the number of axes adn the number
723    of points. */
724    ptr = astGetPoints( this );
725    nc = astGetNcoord( this );
726    np = astGetNpoint( this );
727 
728 /* Check the pointers can be used safely. */
729    if( astOK ) {
730 
731 /* Loop round each axis. */
732       for( ic = 0; ic < nc; ic++ ) {
733 
734 /* Initialise the bounds for this axis. */
735          lb = AST__BAD;
736          ub = AST__BAD;
737 
738 /* Search for the first good point. Use it to initialise the bounds and
739    break out of the loop. */
740          p = ptr[ ic ];
741          for( ip = 0; ip < np; ip++,p++ ) {
742             if( *p != AST__BAD ) {
743                lb = ub = *p;
744                break;
745             }
746          }
747 
748 /* Search through the remaining points. Update the bounds if the axis
749    value is good. */
750          for( ; ip < np; ip++,p++ ) {
751             if( *p != AST__BAD ) {
752                if( *p < lb ) {
753                   lb = *p;
754                } else if( *p > ub ) {
755                   ub = *p;
756                }
757             }
758          }
759 
760 /* Store the returned bounds. */
761          lbnd[ ic ] = lb;
762          ubnd[ ic ] = ub;
763       }
764    }
765 }
766 
CheckPerm(AstPointSet * this,const int * perm,const char * method,int * status)767 static void CheckPerm( AstPointSet *this, const int *perm, const char *method, int *status ) {
768 /*
769 *+
770 *  Name:
771 *     astCheckPerm
772 
773 *  Purpose:
774 *     Check that an array contains a valid permutation.
775 
776 *  Type:
777 *     Protected virtual function.
778 
779 *  Synopsis:
780 *     #include "pointset.h"
781 *     void astCheckPerm( AstPointSet *this, const int *perm, const char *method )
782 
783 *  Class Membership:
784 *     PointSet method.
785 
786 *  Description:
787 *     This function checks the validity of a permutation array that
788 *     will be used to permute the order of a PointSet's axes. If the
789 *     permutation specified by the array is not valid, an error is
790 *     reported and the global error status is set. Otherwise, the
791 *     function returns without further action.
792 
793 *  Parameters:
794 *     this
795 *        Pointer to the PointSet.
796 *     perm
797 *        Pointer to an array of integers with the same number of
798 *        elements as there are axes in the PointSet. For each axis, the
799 *        corresponding integer gives the (zero based) axis index to be
800 *        used to identify the axis values for that axis (using the
801 *        un-permuted axis numbering). To be valid, the integers in
802 *        this array should therefore all lie in the range zero to
803 *        (ncoord-1) inclusive, where "ncoord" is the number of PointSet
804 *        axes, and each value should occur exactly once.
805 *     method
806 *        Pointer to a constant null-terminated character string
807 *        containing the name of the method that invoked this function
808 *        to validate a permutation array. This method name is used
809 *        solely for constructing error messages.
810 
811 *  Notes:
812 *     - Error messages issued by this function refer to the external
813 *     (public) numbering system used for axes (which is one-based),
814 *     whereas zero-based axis indices are used internally.
815 *-
816 */
817 
818 /* Local Variables: */
819    int *there;                   /* Pointer to temporary array */
820    int coord;                     /* Loop counter for axes */
821    int ncoord;                   /* Number of PointSet axes */
822    int valid;                    /* Permutation array is valid? */
823 
824 /* Check the global error status. */
825    if ( !astOK ) return;
826 
827 /* Initialise. */
828    valid = 1;
829 
830 /* Obtain the number of PointSet axes and allocate a temporary array of
831    integers with the same number of elements. */
832    ncoord = astGetNcoord( this );
833    there = astMalloc( sizeof( int ) * (size_t) ncoord );
834    if ( astOK ) {
835 
836 /* Initialise the temporary array to zero. */
837       for ( coord = 0; coord < ncoord; coord++ ) there[ coord ] = 0;
838 
839 /* Scan the permutation array, checking that each permuted axis index it
840    contains is within the correct range. Note an error and quit checking
841    if an invalid value is found. */
842       for ( coord = 0; coord < ncoord; coord++ ) {
843          if ( ( perm[ coord ] < 0 ) || ( perm[ coord ] >= ncoord ) ) {
844             valid = 0;
845             break;
846 
847 /* Use the temporary array to count how many times each valid axis index
848    occurs. */
849 	 } else {
850             there[ perm[ coord ] ]++;
851 	 }
852       }
853 
854 /* If all the axis indices were within range, check to ensure that each value
855    occurred only once. */
856       if ( valid ) {
857          for ( coord = 0; coord < ncoord; coord++ ) {
858 
859 /* Note an error and quit checking if any value did not occur exactly once. */
860             if ( there[ coord ] != 1 ) {
861                valid = 0;
862                break;
863 	    }
864 	 }
865       }
866    }
867 
868 /* Free the temporary array. */
869    there = astFree( there );
870 
871 /* If an invalid permutation was detected and no other error has
872    occurred, then report an error (note we convert to one-based axis
873    numbering in the error message). */
874    if ( !valid && astOK ) {
875       astError( AST__PRMIN, "%s(%s): Invalid coordinate permutation array.", status,
876                 method, astGetClass( this ) );
877       astError( AST__PRMIN, "Each coordinate index should lie in the range 1 to %d "
878                 "and should occur only once.", status, ncoord );
879    }
880 }
881 
ClearAttrib(AstObject * this_object,const char * attrib,int * status)882 static void ClearAttrib( AstObject *this_object, const char *attrib, int *status ) {
883 /*
884 *  Name:
885 *     ClearAttrib
886 
887 *  Purpose:
888 *     Clear an attribute value for a PointSet.
889 
890 *  Type:
891 *     Private function.
892 
893 *  Synopsis:
894 *     #include "pointset.h"
895 *     void ClearAttrib( AstObject *this, const char *attrib, int *status )
896 
897 *  Class Membership:
898 *     PointSet member function (over-rides the astClearAttrib
899 *     protected method inherited from the Object class).
900 
901 *  Description:
902 *     This function clears the value of a specified attribute for a
903 *     PointSet, so that the default value will subsequently be used.
904 
905 *  Parameters:
906 *     this
907 *        Pointer to the PointSet.
908 *     attrib
909 *        Pointer to a null-terminated string specifying the attribute
910 *        name.  This should be in lower case with no surrounding white
911 *        space.
912 *     status
913 *        Pointer to the inherited status variable.
914 */
915 
916 /* Local Variables: */
917    AstPointSet *this;            /* Pointer to the PointSet structure */
918 
919 /* Check the global error status. */
920    if ( !astOK ) return;
921 
922 /* Obtain a pointer to the PointSet structure. */
923    this = (AstPointSet *) this_object;
924 
925 /* Check the attribute name and clear the appropriate attribute. */
926 
927 /* Test if the name matches any of the read-only attributes of this
928    class. If it does, then report an error. */
929    if ( !strcmp( attrib, "ncoord" ) ||
930         !strcmp( attrib, "npoint" ) ) {
931       astError( AST__NOWRT, "astClear: Invalid attempt to clear the \"%s\" "
932                 "value for a %s.", status, attrib, astGetClass( this ) );
933       astError( AST__NOWRT, "This is a read-only attribute." , status);
934 
935 /* If the attribute is still not recognised, pass it on to the parent
936    method for further interpretation. */
937    } else {
938       (*parent_clearattrib)( this_object, attrib, status );
939    }
940 }
941 
Equal(AstObject * this_object,AstObject * that_object,int * status)942 static int Equal( AstObject *this_object, AstObject *that_object, int *status ) {
943 /*
944 *  Name:
945 *     Equal
946 
947 *  Purpose:
948 *     Test if two PointSets are equivalent.
949 
950 *  Type:
951 *     Private function.
952 
953 *  Synopsis:
954 *     #include "pointset.h"
955 *     int Equal( AstObject *this, AstObject *that, int *status )
956 
957 *  Class Membership:
958 *     PointSet member function (over-rides the astEqual protected
959 *     method inherited from the Object class).
960 
961 *  Description:
962 *     This function returns a boolean result (0 or 1) to indicate whether
963 *     two PointSets are equivalent.
964 
965 *  Parameters:
966 *     this
967 *        Pointer to the first PointSet.
968 *     that
969 *        Pointer to the second PointSet.
970 *     status
971 *        Pointer to the inherited status variable.
972 
973 *  Returned Value:
974 *     One if the PointSets are equivalent, zero otherwise.
975 
976 *  Notes:
977 *     - The two PointSets are considered equivalent if they have the same
978 *     number of points, the same number of axis values per point, and the
979 *     same axis values to within the absolute tolerance specified by the
980 *     Accuracy attribute of the two PointSets.
981 *     - A value of zero will be returned if this function is invoked
982 *     with the global status set, or if it should fail for any reason.
983 */
984 
985 /* Local constants: */
986 #define SMALL sqrt(DBL_MIN)
987 
988 /* Local Variables: */
989    AstPointSet *that;         /* Pointer to the second PointSet structure */
990    AstPointSet *this;         /* Pointer to the first PointSet structure */
991    double **ptr_that;         /* Pointer to axis values in second PointSet */
992    double **ptr_this;         /* Pointer to axis values in first PointSet */
993    double *p_that;            /* Pointer to next axis value in second PointSet */
994    double *p_this;            /* Pointer to next axis value in first PointSet */
995    double acc1;               /* Absolute accuracy for 1st PointSet axis value */
996    double acc2;               /* Absolute accuracy for 2nd PointSet axis value */
997    double acc;                /* Combined absolute accuracy */
998    double acc_that;           /* PointAccuracy attribute for 2nd PointSet */
999    double acc_this;           /* PointAccuracy attribute for 1st PointSet */
1000    int ic;                    /* Axis index */
1001    int ip;                    /* Point index */
1002    int nc;                    /* No. of axis values per point */
1003    int np;                    /* No. of points in each PointSet */
1004    int result;                /* Result value to return */
1005 
1006 /* Initialise. */
1007    result = 0;
1008 
1009 /* Check the global error status. */
1010    if ( !astOK ) return result;
1011 
1012 /* Invoke the Equal method inherited from the parent Object class. This checks
1013    that the Objects are both of the same class (amongst other things). */
1014    if( (*parent_equal)( this_object, that_object, status ) ) {
1015 
1016 /* Obtain pointers to the two PointSet structures. */
1017       this = (AstPointSet *) this_object;
1018       that = (AstPointSet *) that_object;
1019 
1020 /* Check the number of points and the number of axis values per point are
1021    equal in the two PointSets. */
1022       np = astGetNpoint( this );
1023       nc = astGetNcoord( this );
1024       if( np == astGetNpoint( that ) && nc == astGetNcoord( that ) ) {
1025 
1026 /* Get pointers to the axis values.*/
1027          ptr_this = astGetPoints( this );
1028          ptr_that = astGetPoints( that );
1029          if( astOK ) {
1030 
1031 /* Assume the PointSets are equal until proven otherwise. */
1032             result = 1;
1033 
1034 /* Loop round each axis, until we find a difference. */
1035             for( ic = 0; ic < nc && result; ic++ ) {
1036 
1037 /* Get pointers to the next value for this axis. */
1038                p_this = ptr_this[ ic ];
1039                p_that = ptr_that[ ic ];
1040 
1041 /* Get the absolute accuracies for this axis. The default value for this
1042    attribute is AST__BAD. */
1043                acc_this = astGetPointAccuracy( this, ic );
1044                acc_that = astGetPointAccuracy( that, ic );
1045 
1046 /* If both accuracies are available, combine them in quadrature. */
1047                if( acc_this != AST__BAD && acc_that != AST__BAD ) {
1048                   acc = sqrt( acc_this*acc_this + acc_that*acc_that );
1049 
1050 /* Loop round all points on this axis */
1051                   for( ip = 0; ip < np; ip++, p_this++, p_that++ ){
1052 
1053 /* If either value is bad we do not need to compare values. */
1054                      if( *p_this == AST__BAD || *p_that == AST__BAD ) {
1055 
1056 /* If one value is bad and one is good, they differ, so break. If both
1057    values are bad they are equal so we continue. */
1058                         if( *p_this != AST__BAD || *p_that != AST__BAD ) {
1059                            result = 0;
1060                            break;
1061                         }
1062 
1063 /* Otherwise (if both axis values are good), compare axis values, and break if
1064    they differ by more than the absolute accuracy. */
1065                      } else if( fabs( *p_this - *p_that ) > acc ) {
1066                         result = 0;
1067                         break;
1068                      }
1069                   }
1070 
1071 /* If either accuracy is unavailable, we use a default relative accuracy. */
1072                } else {
1073 
1074 /* Loop round all points on this axis */
1075                   for( ip = 0; ip < np; ip++, p_this++, p_that++ ){
1076 
1077 /* If either value is bad we do not need to compare values. */
1078                      if( *p_this == AST__BAD || *p_that == AST__BAD ) {
1079 
1080 /* If one value is bad and one is good, they differ, so break. If both
1081    values are bad they are equal so we continue. */
1082                         if( *p_this != AST__BAD || *p_that != AST__BAD ) {
1083                            result = 0;
1084                            break;
1085                         }
1086 
1087 /* Otherwise (if both axis values are good), find the absolute error for
1088    both values. */
1089                      } else {
1090 
1091                         if( acc_this == AST__BAD ) {
1092                            acc1 = fabs(*p_this)*DBL_EPSILON;
1093                            if( acc1 < SMALL ) acc1 = SMALL;
1094                            acc1 *= 1.0E3;
1095                         } else {
1096                            acc1 = acc_this;
1097                         }
1098 
1099                         if( acc_that == AST__BAD ) {
1100                            acc2 = fabs(*p_that)*DBL_EPSILON;
1101                            if( acc2 < SMALL ) acc2 = SMALL;
1102                            acc2 *= 1.0E3;
1103                         } else {
1104                            acc2 = acc_that;
1105                         }
1106 
1107 /* Combine them in quadrature. */
1108                         acc = sqrt( acc1*acc1 + acc2*acc2 );
1109 
1110 /* Compare axis values, and break if they differ by more than the
1111    absolute accuracy. */
1112                         if( fabs( *p_this - *p_that ) > acc ) {
1113                            result = 0;
1114                            break;
1115                         }
1116                      }
1117                   }
1118                }
1119             }
1120          }
1121       }
1122    }
1123 
1124 /* If an error occurred, clear the result value. */
1125    if ( !astOK ) result = 0;
1126 
1127 /* Return the result, */
1128    return result;
1129 #undef SMALL
1130 }
1131 
GetAttrib(AstObject * this_object,const char * attrib,int * status)1132 static const char *GetAttrib( AstObject *this_object, const char *attrib, int *status ) {
1133 /*
1134 *  Name:
1135 *     GetAttrib
1136 
1137 *  Purpose:
1138 *     Get the value of a specified attribute for a PointSet.
1139 
1140 *  Type:
1141 *     Private function.
1142 
1143 *  Synopsis:
1144 *     #include "pointset.h"
1145 *     const char *GetAttrib( AstObject *this, const char *attrib, int *status )
1146 
1147 *  Class Membership:
1148 *     PointSet member function (over-rides the protected astGetAttrib
1149 *     method inherited from the Object class).
1150 
1151 *  Description:
1152 *     This function returns a pointer to the value of a specified
1153 *     attribute for a PointSet, formatted as a character string.
1154 
1155 *  Parameters:
1156 *     this
1157 *        Pointer to the PointSet.
1158 *     attrib
1159 *        Pointer to a null-terminated string containing the name of
1160 *        the attribute whose value is required. This name should be in
1161 *        lower case, with all white space removed.
1162 *     status
1163 *        Pointer to the inherited status variable.
1164 
1165 *  Returned Value:
1166 *     - Pointer to a null-terminated string containing the attribute
1167 *     value.
1168 
1169 *  Notes:
1170 *     - The returned string pointer may point at memory allocated
1171 *     within the PointSet, or at static memory. The contents of the
1172 *     string may be over-written or the pointer may become invalid
1173 *     following a further invocation of the same function or any
1174 *     modification of the PointSet. A copy of the string should
1175 *     therefore be made if necessary.
1176 *     - A NULL pointer will be returned if this function is invoked
1177 *     with the global error status set, or if it should fail for any
1178 *     reason.
1179 */
1180 
1181 /* Local Variables: */
1182    astDECLARE_GLOBALS           /* Pointer to thread-specific global data */
1183    AstPointSet *this;            /* Pointer to the PointSet structure */
1184    const char *result;           /* Pointer value to return */
1185    int ncoord;                   /* Ncoord attribute value */
1186    int npoint;                   /* Npoint attribute value */
1187 
1188 /* Initialise. */
1189    result = NULL;
1190 
1191 /* Check the global error status. */
1192    if ( !astOK ) return result;
1193 
1194 /* Get a pointer to the thread specific global data structure. */
1195    astGET_GLOBALS(this_object);
1196 
1197 /* Obtain a pointer to the PointSet structure. */
1198    this = (AstPointSet *) this_object;
1199 
1200 /* Compare "attrib" with each recognised attribute name in turn,
1201    obtaining the value of the required attribute. If necessary, write
1202    the value into "getattrib_buff" as a null-terminated string in an appropriate
1203    format.  Set "result" to point at the result string. */
1204 
1205 /* Ncoord. */
1206 /* ------- */
1207    if ( !strcmp( attrib, "ncoord" ) ) {
1208       ncoord = astGetNcoord( this );
1209       if ( astOK ) {
1210          (void) sprintf( getattrib_buff, "%d", ncoord );
1211          result = getattrib_buff;
1212       }
1213 
1214 /* Npoint. */
1215 /* ------- */
1216    } else if ( !strcmp( attrib, "npoint" ) ) {
1217       npoint = astGetNpoint( this );
1218       if ( astOK ) {
1219          (void) sprintf( getattrib_buff, "%d", npoint );
1220          result = getattrib_buff;
1221       }
1222 
1223 /* If the attribute name was not recognised, pass it on to the parent
1224    method for further interpretation. */
1225    } else {
1226       result = (*parent_getattrib)( this_object, attrib, status );
1227    }
1228 
1229 /* Return the result. */
1230    return result;
1231 }
1232 
GetObjSize(AstObject * this_object,int * status)1233 static int GetObjSize( AstObject *this_object, int *status ) {
1234 /*
1235 *  Name:
1236 *     GetObjSize
1237 
1238 *  Purpose:
1239 *     Return the in-memory size of an Object.
1240 
1241 *  Type:
1242 *     Private function.
1243 
1244 *  Synopsis:
1245 *     #include "pointset.h"
1246 *     int GetObjSize( AstObject *this, int *status )
1247 
1248 *  Class Membership:
1249 *     PointSet member function (over-rides the astGetObjSize protected
1250 *     method inherited from the Object class).
1251 
1252 *  Description:
1253 *     This function returns the in-memory size of the supplied PointSet,
1254 *     in bytes.
1255 
1256 *  Parameters:
1257 *     this
1258 *        Pointer to the Object.
1259 *     status
1260 *        Pointer to the inherited status variable.
1261 
1262 *  Returned Value:
1263 *     The Object size, in bytes.
1264 
1265 *  Notes:
1266 *     - A value of zero will be returned if this function is invoked
1267 *     with the global status set, or if it should fail for any reason.
1268 */
1269 
1270 /* Local Variables: */
1271    AstPointSet *this;         /* Pointer to PointSet structure */
1272    int result;                /* Result value to return */
1273 
1274 /* Initialise. */
1275    result = 0;
1276 
1277 /* Check the global error status. */
1278    if ( !astOK ) return result;
1279 
1280 /* Obtain a pointers to the PointSet structure. */
1281    this = (AstPointSet *) this_object;
1282 
1283 /* Invoke the GetObjSize method inherited from the parent class, and then
1284    add on any components of the class structure defined by thsi class
1285    which are stored in dynamically allocated memory. */
1286    result = (*parent_getobjsize)( this_object, status );
1287 
1288    result += astTSizeOf( this->ptr );
1289    result += astTSizeOf( this->values );
1290    result += astTSizeOf( this->acc );
1291 
1292 /* If an error occurred, clear the result value. */
1293    if ( !astOK ) result = 0;
1294 
1295 /* Return the result, */
1296    return result;
1297 }
1298 
GetNcoord(const AstPointSet * this,int * status)1299 static int GetNcoord( const AstPointSet *this, int *status ) {
1300 /*
1301 *+
1302 *  Name:
1303 *     astGetNcoord
1304 
1305 *  Purpose:
1306 *     Get the number of coordinate values per point from a PointSet.
1307 
1308 *  Type:
1309 *     Protected virtual function.
1310 
1311 *  Synopsis:
1312 *     #include "pointset.h"
1313 *     int astGetNcoord( const AstPointSet *this )
1314 
1315 *  Class Membership:
1316 *     PointSet method.
1317 
1318 *  Description:
1319 *     This function returns the number of coordinate values per point (1 or
1320 *     more) for a PointSet.
1321 
1322 *  Parameters:
1323 *     this
1324 *        Pointer to the PointSet.
1325 
1326 *  Returned Value:
1327 *     The number of coordinate values per point.
1328 
1329 *  Notes:
1330 *     -  A value of zero is returned if this function is invoked with the
1331 *     global error status set, or if it should fail for any reason.
1332 *-
1333 */
1334 
1335 /* Check the global error status. */
1336    if ( !astOK ) return 0;
1337 
1338 /* Return the number of coordinate values. */
1339    return this->ncoord;
1340 }
1341 
GetNpoint(const AstPointSet * this,int * status)1342 static int GetNpoint( const AstPointSet *this, int *status ) {
1343 /*
1344 *+
1345 *  Name:
1346 *     astGetNpoint
1347 
1348 *  Purpose:
1349 *     Get the number of points in a PointSet.
1350 
1351 *  Type:
1352 *     Protected virtual function.
1353 
1354 *  Synopsis:
1355 *     #include "pointset.h"
1356 *     int astGetNpoint( const AstPointSet *this )
1357 
1358 *  Class Membership:
1359 *     PointSet method.
1360 
1361 *  Description:
1362 *     This function returns the number of points (1 or more) in a PointSet.
1363 
1364 *  Parameters:
1365 *     this
1366 *        Pointer to the PointSet.
1367 
1368 *  Returned Value:
1369 *     The number of points.
1370 
1371 *  Notes:
1372 *     -  A value of zero is returned if this function is invoked with the
1373 *     global error status set, or if it should fail for any reason.
1374 *-
1375 */
1376 
1377 /* Check the global error status. */
1378    if ( !astOK ) return 0;
1379 
1380 /* Return the number of points. */
1381    return this->npoint;
1382 }
1383 
GetPoints(AstPointSet * this,int * status)1384 static double **GetPoints( AstPointSet *this, int *status ) {
1385 /*
1386 *+
1387 *  Name:
1388 *     astGetPoints
1389 
1390 *  Purpose:
1391 *     Get a pointer for the coordinate values associated with a PointSet.
1392 
1393 *  Type:
1394 *     Protected virtual function.
1395 
1396 *  Synopsis:
1397 *     #include "pointset.h"
1398 *     double **astGetPoints( AstPointSet *this )
1399 
1400 *  Class Membership:
1401 *     PointSet method.
1402 
1403 *  Description:
1404 *     This function returns a pointer which grants access to the coordinate
1405 *     values associated with a PointSet. If the PointSet has previously had
1406 *     coordinate values associated with it, this pointer will identify these
1407 *     values. Otherwise, it will point at a newly-allocated region of memory
1408 *     (associated with the PointSet) in which new coordinate values may be
1409 *     stored.
1410 
1411 *  Parameters:
1412 *     this
1413 *        Pointer to the PointSet.
1414 
1415 *  Returned Value:
1416 *     A pointer to an array of type double* with ncoord elements (where ncoord
1417 *     is the number of coordinate values per point). Each element of this array
1418 *     points at an array of double, of size npoint (where npoint is the number
1419 *     of points in the PointSet), containing the values of that coordinate for
1420 *     each point in the set. Hence, the value of the i'th coordinate for the
1421 *     j'th point (where i and j are counted from zero) is given by ptr[i][j]
1422 *     where ptr is the returned pointer value.
1423 
1424 *  Notes:
1425 *     -  The returned pointer points at an array of pointers allocated
1426 *     internally within the PointSet. The values in this array may be changed
1427 *     by the caller, who is reponsible for ensuring that they continue to
1428 *     point at valid arrays of coordinate values.
1429 *     -  No attempt should be made to de-allocate memory allocated by a
1430 *     PointSet to store coordinate values or pointers to them. This memory
1431 *     will be freed when the PointSet is deleted.
1432 *     -  No count is kept of the number of pointers issued for the PointSet
1433 *     coordinate values. The caller must keep track of these.
1434 *     -  A NULL pointer is returned if this function is invoked with the
1435 *     global error status set, or if it should fail for any reason.
1436 *-
1437 */
1438 
1439 /* Local Variables: */
1440    int i;                        /* Loop counter for coordinates */
1441    int nval;                     /* Number of values to be stored */
1442 
1443 /* Check the global error status. */
1444    if ( !astOK ) return NULL;
1445 
1446 /* If the PointSet has an existing array of pointers (which point at coordinate
1447    values), we will simply return a pointer to it. Otherwise, we must allocate
1448    space to hold new coordinate values. */
1449    if( !this->ptr ) {
1450 
1451 /* Determine the number of coordinate values to be stored and allocate memory
1452    to hold them, storing the pointer to this values array in the PointSet
1453    structure. */
1454       nval = this->npoint * this->ncoord;
1455       this->values = (double *) astMalloc( sizeof( double ) * (size_t) nval );
1456 
1457 #ifdef DEBUG
1458       for( i = 0; i < nval; i++ ) this->values[ i ] = 0.0;
1459 #endif
1460 
1461 /* If OK, also allocate memory for the array of pointers into this values
1462    array, storing a pointer to this pointer array in the PointSet structure. */
1463       if ( astOK ) {
1464          this->ptr = (double **) astMalloc( sizeof( double * )
1465                                             * (size_t) this->ncoord );
1466 
1467 /* If OK, initialise the pointer array to point into the values array. */
1468          if ( astOK ) {
1469             for ( i = 0; i < this->ncoord; i++ ) {
1470                this->ptr[ i ] = this->values + ( i * this->npoint );
1471             }
1472 
1473 /* If we failed to allocate the pointer array, then free the values array. */
1474          } else {
1475             this->values = (double *) astFree( (void *) this->values );
1476          }
1477       }
1478 
1479 #ifdef DEBUG
1480    } else {
1481 
1482 /* Check for bad values */
1483       if( this->values ) {
1484          int i, j;
1485          for( i = 0; astOK && i < this->ncoord; i++ ) {
1486             for( j = 0; j < this->npoint; j++ ) {
1487                if( !finite( (this->ptr)[ i ][ j ] ) ) {
1488                   astError( AST__INTER, "astGetPoints(PointSet): Non-finite "
1489                             "axis value returned.", status);
1490                   break;
1491                }
1492             }
1493          }
1494       }
1495 
1496 #endif
1497    }
1498 
1499 /* Return the required pointer. */
1500    return this->ptr;
1501 }
1502 
astInitPointSetVtab_(AstPointSetVtab * vtab,const char * name,int * status)1503 void astInitPointSetVtab_(  AstPointSetVtab *vtab, const char *name, int *status ) {
1504 /*
1505 *+
1506 *  Name:
1507 *     astInitPointSetVtab
1508 
1509 *  Purpose:
1510 *     Initialise a virtual function table for a PointSet.
1511 
1512 *  Type:
1513 *     Protected function.
1514 
1515 *  Synopsis:
1516 *     #include "pointset.h"
1517 *     void astInitPointSetVtab( AstPointSetVtab *vtab, const char *name )
1518 
1519 *  Class Membership:
1520 *     PointSet vtab initialiser.
1521 
1522 *  Description:
1523 *     This function initialises the component of a virtual function
1524 *     table which is used by the PointSet class.
1525 
1526 *  Parameters:
1527 *     vtab
1528 *        Pointer to the virtual function table. The components used by
1529 *        all ancestral classes will be initialised if they have not already
1530 *        been initialised.
1531 *     name
1532 *        Pointer to a constant null-terminated character string which contains
1533 *        the name of the class to which the virtual function table belongs (it
1534 *        is this pointer value that will subsequently be returned by the Object
1535 *        astClass function).
1536 *-
1537 */
1538 
1539 /* Local Variables: */
1540    AstObjectVtab *object;        /* Pointer to Object component of Vtab */
1541    astDECLARE_GLOBALS            /* Pointer to thread-specific global data */
1542    const char *envvar;           /* Pointer to environment variable value */
1543    size_t i;                     /* Index of next byte in NaN value */
1544    unsigned char *p;             /* Pointer to next byte in NaN value */
1545 
1546 /* Check the local error status. */
1547    if ( !astOK ) return;
1548 
1549 /* Get a pointer to the thread specific global data structure. */
1550    astGET_GLOBALS(NULL);
1551 
1552 /* Initialize the component of the virtual function table used by the
1553    parent class. */
1554    astInitObjectVtab( (AstObjectVtab *) vtab, name );
1555 
1556 /* Store a unique "magic" value in the virtual function table. This
1557    will be used (by astIsAPointSet) to determine if an object belongs
1558    to this class.  We can conveniently use the address of the (static)
1559    class_check variable to generate this unique value. */
1560    vtab->id.check = &class_check;
1561    vtab->id.parent = &(((AstObjectVtab *) vtab)->id);
1562 
1563 /* Initialise member function pointers. */
1564 /* ------------------------------------ */
1565 /* Store pointers to the member functions (implemented here) that
1566    provide virtual methods for this class. */
1567    vtab->AppendPoints = AppendPoints;
1568    vtab->BndPoints = BndPoints;
1569    vtab->GetNcoord = GetNcoord;
1570    vtab->GetNpoint = GetNpoint;
1571    vtab->GetPoints = GetPoints;
1572    vtab->PermPoints = PermPoints;
1573    vtab->ReplaceNaN = ReplaceNaN;
1574    vtab->SetPoints = SetPoints;
1575    vtab->SetNpoint = SetNpoint;
1576    vtab->SetSubPoints = SetSubPoints;
1577 
1578    vtab->GetPointAccuracy = GetPointAccuracy;
1579    vtab->SetPointAccuracy = SetPointAccuracy;
1580    vtab->TestPointAccuracy = TestPointAccuracy;
1581    vtab->ClearPointAccuracy = ClearPointAccuracy;
1582 
1583 /* Save the inherited pointers to methods that will be extended, and
1584    replace them with pointers to the new member functions. */
1585    object = (AstObjectVtab *) vtab;
1586 
1587    parent_clearattrib = object->ClearAttrib;
1588    object->ClearAttrib = ClearAttrib;
1589    parent_getattrib = object->GetAttrib;
1590    object->GetAttrib = GetAttrib;
1591    parent_setattrib = object->SetAttrib;
1592    object->SetAttrib = SetAttrib;
1593    parent_testattrib = object->TestAttrib;
1594    object->TestAttrib = TestAttrib;
1595    parent_equal = object->Equal;
1596    object->Equal = Equal;
1597    parent_getobjsize = object->GetObjSize;
1598    object->GetObjSize = GetObjSize;
1599 
1600 /* Declare the copy constructor, destructor and class dump function. */
1601    astSetCopy( vtab, Copy );
1602    astSetDelete( vtab, Delete );
1603    astSetDump( vtab, Dump, "PointSet", "Container for a set of points" );
1604 
1605 /* Calculate single and double precision NaN values and store in static
1606    module variables. Setting all bits to 1 produces a quiet NaN. */
1607    LOCK_MUTEX1
1608    p = (unsigned char *) &ast_nan;
1609    for( i = 0; i < sizeof( ast_nan ); i++ ) *(p++) = 255;
1610    p = (unsigned char *) &ast_nanf;
1611    for( i = 0; i < sizeof( ast_nanf ); i++ ) *(p++) = 255;
1612 
1613 /* See what  action the astReplaceNaN method should perform. This
1614    is determined by the value of the AST_REPLACE_NAN environment
1615    variable. Not set = do not check for NaNs, "1" = replace NaNs with
1616    AST__BAD silently, anything else = report an error if any NaNs are
1617    encountered. */
1618    if( replace_nan == -1 ) {
1619       envvar = getenv( "AST_REPLACE_NAN" );
1620       if( !envvar ) {
1621          replace_nan = IGNORE_NANS;
1622       } else if( !strcmp( envvar, "1" ) ) {
1623          replace_nan = REPLACE_NANS;
1624       } else {
1625          replace_nan = REPORT_NANS;
1626       }
1627    }
1628    UNLOCK_MUTEX1
1629 
1630 /* If we have just initialised the vtab for the current class, indicate
1631    that the vtab is now initialised, and store a pointer to the class
1632    identifier in the base "object" level of the vtab. */
1633    if( vtab == &class_vtab ) {
1634       class_init = 1;
1635       astSetVtabClassIdentifier( vtab, &(vtab->id) );
1636    }
1637 }
1638 
astCheckNaN_(double value)1639 double astCheckNaN_( double value ) {
1640 /*
1641 *+
1642 *  Name:
1643 *     astCheckNaN
1644 
1645 *  Purpose:
1646 *     Substitute a NaN for a supplied value if the supplied value is
1647 *     AST__NAN.
1648 
1649 *  Type:
1650 *     Protected function.
1651 
1652 *  Synopsis:
1653 *     #include "pointset.h"
1654 *     double astCheckNaN( double value )
1655 
1656 *  Class Membership:
1657 *     PointSet method.
1658 
1659 *  Description:
1660 *     If the supplied double is AST__NAN then this function returns an
1661 *     IEEE double precision NaN value. Otherwise it returns the supplied
1662 *     value.
1663 
1664 *  Parameters:
1665 *     valuethis
1666 *        The value to check.
1667 
1668 *  Returned Value:
1669 *     The suppleid value, or NaN.
1670 
1671 *  Notes:
1672 *     - This function does not check the inherited status.
1673 
1674 *-
1675 */
1676    return ( value == AST__NAN ) ? ast_nan : value;
1677 }
1678 
astCheckNaNF_(float value)1679 float astCheckNaNF_( float value ) {
1680 /*
1681 *+
1682 *  Name:
1683 *     astCheckNaNF
1684 
1685 *  Purpose:
1686 *     Substitute a NaN for a supplied value if the supplied value is
1687 *     AST__NANF.
1688 
1689 *  Type:
1690 *     Protected function.
1691 
1692 *  Synopsis:
1693 *     #include "pointset.h"
1694 *     float astCheckNaNF( float value )
1695 
1696 *  Class Membership:
1697 *     PointSet method.
1698 
1699 *  Description:
1700 *     If the supplied float is AST__NANF then this function returns an
1701 *     IEEE single precision NaN value. Otherwise it returns the supplied
1702 *     value.
1703 
1704 *  Parameters:
1705 *     valuethis
1706 *        The value to check.
1707 
1708 *  Returned Value:
1709 *     The suppleid value, or NaN.
1710 
1711 *  Notes:
1712 *     - This function does not check the inherited status.
1713 
1714 *-
1715 */
1716    return ( value == AST__NANF ) ? ast_nanf : value;
1717 }
1718 
PermPoints(AstPointSet * this,int forward,const int perm[],int * status)1719 static void PermPoints( AstPointSet *this, int forward, const int perm[], int *status ) {
1720 /*
1721 *+
1722 *  Name:
1723 *     astPermPoints
1724 
1725 *  Purpose:
1726 *     Permute the order of a PointSet's axes.
1727 
1728 *  Type:
1729 *     Protected virtual function.
1730 
1731 *  Synopsis:
1732 *     #include "pointset.h"
1733 *     void astPermPoints( AstPointSet *this, int forward, const int perm[] )
1734 
1735 *  Class Membership:
1736 *     PointSet method.
1737 
1738 *  Description:
1739 *     This function permutes the order in which a PointSet's axes occur.
1740 
1741 *  Parameters:
1742 *     this
1743 *        Pointer to the PointSet.
1744 *     forward
1745 *        The direction in which the permutation is to be applied. This
1746 *        controls the use of the "perm" arrays. If a non-zero value is
1747 *        given, then the indices into the "perm" array correspond to the
1748 *        indices of the coordinates in the returned PointSet, and the
1749 *        values stored in the "perm" array correspond to the indices of
1750 *        the coordinates in the supplied PointSet. If a zero value is
1751 *        given, then the indices into the "perm" array correspond to the
1752 *        indices of the coordinates in the supplied PointSet, and the
1753 *        values stored in the "perm" array correspond to the indices of
1754 *        the coordinates in the returnedPointSet.
1755 *     perm
1756 *        An array of int (with one element for each axis of the PointSet)
1757 *        which lists the axes in their new order. How this array is use
1758 *        depends on the value supplied for "forward".
1759 
1760 *  Notes:
1761 *     - Only genuine permutations of the axis order are permitted, so
1762 *     each axis must be referenced exactly once in the "perm" array.
1763 *     - If more than one axis permutation is applied to a PointSet, the
1764 *     effects are cumulative.
1765 *-
1766 */
1767 
1768 /* Local Variables: */
1769    double **old;                 /* Pointer to copy of old pointer array */
1770    int coord;                    /* Loop counter for axes */
1771    int ncoord;                   /* Number of axes */
1772 
1773 /* Check the global error status. Return without action if no data is
1774    associated with the PointSet. */
1775    if ( !astOK || !this->ptr ) return;
1776 
1777 /* Validate the permutation array, to check that it describes a genuine
1778    permutation. */
1779    CheckPerm( this, perm, "astPermPoints", status );
1780 
1781 /* Obtain the number of PointSet axes. */
1782    ncoord = astGetNcoord( this );
1783 
1784 /* Allocate memory and use it to store a copy of the old pointers array for
1785    the PointSet. */
1786    old = astStore( NULL, this->ptr, sizeof( double * ) * (size_t) ncoord );
1787 
1788 /* Apply the new axis permutation cumulatively to the old one and store the
1789    result in the PointSet. */
1790    if ( astOK ) {
1791       if( forward ) {
1792          for ( coord = 0; coord < ncoord; coord++ ) {
1793             this->ptr[ coord ] = old[ perm[ coord ] ];
1794          }
1795       } else {
1796          for ( coord = 0; coord < ncoord; coord++ ) {
1797             this->ptr[ perm[ coord ] ] = old[ coord ];
1798          }
1799       }
1800    }
1801 
1802 /* Free the temporary copy of the old array. */
1803    old = astFree( old );
1804 }
1805 
ReplaceNaN(AstPointSet * this,int * status)1806 static int ReplaceNaN( AstPointSet *this, int *status ) {
1807 /*
1808 *+
1809 *  Name:
1810 *     astReplaceNaN
1811 
1812 *  Purpose:
1813 *     Check for NaNs in a PointSet.
1814 
1815 *  Type:
1816 *     Protected virtual function.
1817 
1818 *  Synopsis:
1819 *     #include "pointset.h"
1820 *     int astReplaceNaNs( AstPointSet *this )
1821 
1822 *  Class Membership:
1823 *     PointSet method.
1824 
1825 *  Description:
1826 *     The behaviour of this method is determined by the AST_REPLACE_NAN
1827 *     environment variable. If AST_REPLACE_NAN is undefined, then the
1828 *     method returns without action. If AST_REPLACE_NAN is "1", then the
1829 *     method examines the supplied PointSet and replaces any NaN values
1830 *     with AST__BAD. If AST_REPLACE_NAN has any other value, any NaNs in
1831 *     the supplied PointSet are still replaced, but in addition an error
1832 *     is reported.
1833 
1834 *  Parameters:
1835 *     this
1836 *        Pointer to the PointSet.
1837 
1838 *  Returned Value:
1839 *     Non-zero if any NaN values were found in the PointSet. If AST_REPLACE_NAN
1840 *     is undefined, then zero is always returned.
1841 
1842 *  Notes:
1843 *    The value of the AST_REPLACE_NAN environment variable is obtained
1844 *    only once, when the PointSet virtual function table is first
1845 *    created. This value is saved for all subsequent invocations of this
1846 *    method.
1847 
1848 *-
1849 */
1850 
1851 /* Local Variables: */
1852    double **ptr;
1853    double *p0;
1854    double *p;
1855    int ic;
1856    int nc;
1857    int np;
1858    int result;
1859 
1860 /* Initialise */
1861    result = 0;
1862 
1863 /* Check the global error status. */
1864    if ( !astOK ) return result;
1865 
1866 /* Unless the AST_REPALCE_NAN environment variable is undefined, check
1867    for NaNs in the supplied PointSet, replacing any with AST__BAD. */
1868    if( replace_nan != IGNORE_NANS ) {
1869       ptr = astGetPoints( this );
1870       if( ptr ) {
1871          nc = astGetNcoord( this );
1872          np = astGetNpoint( this );
1873          for( ic = 0; ic < nc; ic++ ) {
1874             p = ptr[ ic ];
1875             p0 = p + np;
1876             for( ; p < p0; p++ ) {
1877                if( !astISFINITE(*p) ) {
1878                   result = 1;
1879                   *p = AST__BAD;
1880                }
1881             }
1882          }
1883 
1884 /* If any NaNs were found, and AST_REPLACE_NAN is not set to "1", report
1885    an error. */
1886          if( result && replace_nan == REPORT_NANS ) {
1887             astError( AST__ISNAN, "astReplaceNan(%s): One or more NaN values "
1888                       "were encountered within an AST PointSet.", status,
1889                       astGetClass( this ) );
1890          }
1891       }
1892    }
1893 
1894 /* Return the result. */
1895    return result;
1896 }
1897 
SetAttrib(AstObject * this_object,const char * setting,int * status)1898 static void SetAttrib( AstObject *this_object, const char *setting, int *status ) {
1899 /*
1900 *  Name:
1901 *     SetAttrib
1902 
1903 *  Purpose:
1904 *     Set an attribute value for a PointSet.
1905 
1906 *  Type:
1907 *     Private function.
1908 
1909 *  Synopsis:
1910 *     #include "pointset.h"
1911 *     void SetAttrib( AstObject *this, const char *setting )
1912 
1913 *  Class Membership:
1914 *     PointSet member function (over-rides the astSetAttrib protected
1915 *     method inherited from the Object class).
1916 
1917 *  Description:
1918 *     This function assigns an attribute value for a PointSet, the
1919 *     attribute and its value being specified by means of a string of
1920 *     the form:
1921 *
1922 *        "attribute= value "
1923 *
1924 *     Here, "attribute" specifies the attribute name and should be in
1925 *     lower case with no white space present. The value to the right
1926 *     of the "=" should be a suitable textual representation of the
1927 *     value to be assigned and this will be interpreted according to
1928 *     the attribute's data type.  White space surrounding the value is
1929 *     only significant for string attributes.
1930 
1931 *  Parameters:
1932 *     this
1933 *        Pointer to the PointSet.
1934 *     setting
1935 *        Pointer to a null-terminated string specifying the new
1936 *        attribute value.
1937 */
1938 
1939 /* Local Variables: */
1940    AstPointSet *this;            /* Pointer to the PointSet structure */
1941    int len;                      /* Length of setting string */
1942    int nc;                       /* Number of characters read by astSscanf */
1943 
1944 /* Check the global error status. */
1945    if ( !astOK ) return;
1946 
1947 /* Obtain a pointer to the PointSet structure. */
1948    this = (AstPointSet *) this_object;
1949 
1950 /* Obtain the length of the setting string. */
1951    len = (int) strlen( setting );
1952 
1953 /* Test for each recognised attribute in turn, using "astSscanf" to parse
1954    the setting string and extract the attribute value (or an offset to
1955    it in the case of string values). In each case, use the value set
1956    in "nc" to check that the entire string was matched. Once a value
1957    has been obtained, use the appropriate method to set it. */
1958 
1959 /* Define a macro to see if the setting string matches any of the
1960    read-only attributes of this class. */
1961 #define MATCH(attrib) \
1962         ( nc = 0, ( 0 == astSscanf( setting, attrib "=%*[^\n]%n", &nc ) ) && \
1963                   ( nc >= len ) )
1964 
1965 /* Use this macro to report an error if a read-only attribute has been
1966    specified. */
1967    if ( MATCH( "ncoord" ) ||
1968         MATCH( "npoint" ) ) {
1969       astError( AST__NOWRT, "astSet: The setting \"%s\" is invalid for a %s.", status,
1970                 setting, astGetClass( this ) );
1971       astError( AST__NOWRT, "This is a read-only attribute." , status);
1972 
1973 /* If the attribute is still not recognised, pass it on to the parent
1974    method for further interpretation. */
1975    } else {
1976       (*parent_setattrib)( this_object, setting, status );
1977    }
1978 
1979 /* Undefine macros local to this function. */
1980 #undef MATCH
1981 }
1982 
SetNpoint(AstPointSet * this,int npoint,int * status)1983 static void SetNpoint( AstPointSet *this, int npoint, int *status ) {
1984 /*
1985 *+
1986 *  Name:
1987 *     astSetNpoint
1988 
1989 *  Purpose:
1990 *     Reduce the number of points in a PointSet.
1991 
1992 *  Type:
1993 *     Protected virtual function.
1994 
1995 *  Synopsis:
1996 *     #include "pointset.h"
1997 *     void astSetNpoint( AstPointSet *this, int npoint )
1998 
1999 *  Class Membership:
2000 *     PointSet method.
2001 
2002 *  Description:
2003 *     This function reduces the number of points stored in a PointSet.
2004 *     Points with indices beyond the new size will be discarded.
2005 
2006 *  Parameters:
2007 *     this
2008 *        Pointer to the PointSet.
2009 *     npoint
2010 *        The new value for the number of points in the PointSet. Must be
2011 *        less than or equal to the original size of the PointSet, and
2012 *        greater than zero.
2013 *-
2014 */
2015 
2016 /* Check the global error status. */
2017    if ( !astOK ) return;
2018 
2019 /* Check the new size is valid. */
2020    if( npoint < 1 || npoint > this->npoint ) {
2021       astError( AST__NPTIN, "astSetNpoint(%s): Number of points (%d) is "
2022                 "not valid.", status, astGetClass( this ), npoint );
2023       astError( AST__NPTIN, "Should be in the range 1 to %d.", status, this->npoint );
2024 
2025 /* Store the new size. */
2026    } else {
2027       this->npoint = npoint;
2028    }
2029 }
2030 
SetPoints(AstPointSet * this,double ** ptr,int * status)2031 static void SetPoints( AstPointSet *this, double **ptr, int *status ) {
2032 /*
2033 *+
2034 *  Name:
2035 *     astSetPoints
2036 
2037 *  Purpose:
2038 *     Associate coordinate values with a PointSet.
2039 
2040 *  Type:
2041 *     Protected virtual function.
2042 
2043 *  Synopsis:
2044 *     #include "pointset.h"
2045 *     void astSetPoints( AstPointSet *this, double **ptr )
2046 
2047 *  Class Membership:
2048 *     PointSet method.
2049 
2050 *  Description:
2051 *     This function associates coordinate values with a PointSet by storing an
2052 *     array of pointers to the values within the PointSet object. A pointer to
2053 *     this pointer array will later be returned when astGetPoints is used to
2054 *     locate the coordinate values. If values are already associated with the
2055 *     PointSet, the array of pointers to them is over-written by the new values
2056 *     (any internally allocated memory holding the actual coordinate values
2057 *     first being freed).
2058 
2059 *  Parameters:
2060 *     this
2061 *        Pointer to the PointSet.
2062 *     ptr
2063 *        Pointer to an array of type double* with ncoord elements (where ncoord
2064 *        is the number of coordinate values per point in the PointSet). Each
2065 *        element of this array should point at an array of double with npoint
2066 *        elements (where npoint is the number of points in the PointSet),
2067 *        containing the values of that coordinate for each point in the set.
2068 *        Hence, the value of the i'th coordinate for the j'th point (where i
2069 *        and j are counted from zero) should be given by ptr[i][j].
2070 
2071 *  Returned Value:
2072 *     void
2073 
2074 *  Notes:
2075 *     -  It is the caller's responsibility to ensure that the pointer supplied
2076 *     points at a valid array of pointers that point at arrays of coordinate
2077 *     values. This is only superficially validated by this function, which then
2078 *     simply stores a copy of the supplied array of pointers for later use.
2079 *     The caller must also manage any allocation (and freeing) of memory for
2080 *     these coordinate values.
2081 *     -  This functon makes a copy of the array of pointers supplied, but does
2082 *     not copy the coordinate values they point at. If a PointSet containing a
2083 *     copy of the coordinate values is required, internal memory should be
2084 *     allocated within the PointSet by calling astGetPoints before storing any
2085 *     pointer, and then copying the values into this memory. Alternatively,
2086 *     using astCopy to produce a deep copy of a PointSet will also copy the
2087 *     coordinate values.
2088 *     -  A NULL pointer may be supplied as the "ptr" argument, in which case
2089 *     any previously stored array of pointers will be cancelled (and internal
2090 *     memory freed if necessary) and subsequent use of astGetPoints will then
2091 *     cause memory to be allocated internally by the PointSet to hold new
2092 *     values.
2093 *-
2094 */
2095 
2096 /* Local Variables: */
2097    int i;                        /* Loop counter for coordinates */
2098 
2099 /* Check the global error status. */
2100    if ( !astOK ) return;
2101 
2102 /* If the pointer supplied is not NULL, inspect each pointer in the array it
2103    points at to check that none of these are NULL. This validates (to some
2104    extent) the caller's data structure. Report an error and quit checking if a
2105    NULL pointer is found. */
2106    if ( ptr ) {
2107       for ( i = 0; i < this->ncoord; i++ ) {
2108          if ( !ptr[ i ] ) {
2109             astError( AST__PDSIN, "astSetPoints(%s): Invalid NULL pointer in "
2110                       "element %d of array of pointers to coordinate values.", status,
2111                       astGetClass( this ), i );
2112             break;
2113          }
2114       }
2115    }
2116 
2117 /* Do not carry on if the data structure is obviously invalid. */
2118    if ( astOK ) {
2119 
2120 /* Free any memory previously allocated to store coordinate values. */
2121       this->values = (double *) astFree( (void *) this->values );
2122 
2123 /* If a new array of pointers has been provided, (re)allocate memory and store
2124    a copy of the array in it, saving a pointer to this copy in the PointSet
2125    structure. */
2126       if ( ptr ) {
2127          this->ptr = (double **) astStore( (void *) this->ptr,
2128                                            (const void *) ptr,
2129                                            sizeof( double * )
2130                                            * (size_t) this->ncoord );
2131 
2132 /* If no pointer array was provided, free the previous one (if any). */
2133       } else {
2134          this->ptr = (double **) astFree( (void *) this->ptr );
2135       }
2136    }
2137 }
2138 
SetSubPoints(AstPointSet * point1,int point,int coord,AstPointSet * point2,int * status)2139 static void SetSubPoints( AstPointSet *point1, int point, int coord,
2140                           AstPointSet *point2, int *status ) {
2141 /*
2142 *+
2143 *  Name:
2144 *     astSetSubPoints
2145 
2146 *  Purpose:
2147 *     Associate a subset of one PointSet with another PointSet.
2148 
2149 *  Type:
2150 *     Protected virtual function.
2151 
2152 *  Synopsis:
2153 *     #include "pointset.h"
2154 *     void astSetSubPoints( AstPointSet *point1, int point, int coord,
2155 *                           AstPointSet *point2 )
2156 
2157 *  Class Membership:
2158 *     PointSet method.
2159 
2160 *  Description:
2161 *     This function selects a subset of the coordinate values associated with
2162 *     one PointSet and associates them with another PointSet. The second
2163 *     PointSet may then be used to access the subset. Any previous coordinate
2164 *     value association with the second PointSet is replaced.
2165 
2166 *  Parameters:
2167 *     point1
2168 *        Pointer to the first PointSet, from which a subset is to be selected.
2169 *     point
2170 *        The index of the first point (counting from zero) which is to appear
2171 *        in the subset (the number of points is determined by the size of the
2172 *        second PointSet).
2173 *     coord
2174 *        The index of the first coordinate (counting from zero) which is to
2175 *        appear in the subset (the number of coordinates is determined by the
2176 *        size of the second PointSet).
2177 *     point2
2178 *        Second PointSet, with which the subset of coordinate values is to be
2179 *        associated.
2180 
2181 *  Returned Value:
2182 *     void
2183 
2184 *  Notes:
2185 *     -  The range of points and coordinates selected must lie entirely within
2186 *     the first PointSet.
2187 *     -  This function does not make a copy of the coordinate values, but
2188 *     merely stores pointers to the required subset of values associated with
2189 *     the first PointSet. If a PointSet containing a copy of the subset's
2190 *     coordinate values is required, then astCopy should be used to make a
2191 *     deep copy from the second PointSet.
2192 *     -  If the first PointSet does not yet have coordinate values associated
2193 *     with it, then space will be allocated within it to hold values (so that
2194 *     the second PointSet has somewhere to point at).
2195 *-
2196 */
2197 
2198 /* Local Variables: */
2199    double ** ptr2;               /* Pointer to new pointer array */
2200    double **ptr1;                /* Pointer to original pointer array */
2201    int i;                        /* Loop counter for coordinates */
2202    int ncoord1;                  /* Number of coordinates in first PointSet */
2203    int ncoord2;                  /* Number of coordinates in second PointSet */
2204    int npoint1;                  /* Number of points in first PointSet */
2205    int npoint2;                  /* Number of points in second PointSet */
2206 
2207 /* Check the global error status. */
2208    if ( !astOK ) return;
2209 
2210 /* Obtain the sizes of both PointSets. */
2211    npoint1 = astGetNpoint( point1 );
2212    npoint2 = astGetNpoint( point2 );
2213    ncoord1 = astGetNcoord( point1 );
2214    ncoord2 = astGetNcoord( point2 );
2215 
2216 /* Check if the range of points required lies within the first PointSet and
2217    report an error if it does not. */
2218    if ( astOK ) {
2219       if ( ( point < 0 ) || ( point + npoint2 > npoint1 ) ) {
2220          astError( AST__PTRNG, "astSetSubPoints(%s): Range of points in "
2221                    "output %s (%d to %d) lies outside the input %s extent "
2222                    "(0 to %d).", status,
2223                    astGetClass( point1 ), astGetClass( point2 ), point,
2224                    point + npoint2, astGetClass( point1 ), npoint1 );
2225 
2226 /* Similarly check that the range of coordinates is valid. */
2227       } else if ( ( coord < 0 ) || ( coord + ncoord2 > ncoord1 ) ) {
2228          astError( AST__CORNG, "astSetSubPoints(%s): Range of coordinates in "
2229                    "output %s (%d to %d) lies outside the input %s extent "
2230                    "(0 to %d).", status,
2231                    astGetClass( point1 ), astGetClass( point2 ), coord,
2232                    coord + ncoord2, astGetClass( point1 ), ncoord1 );
2233 
2234 /* Obtain a pointer for the coordinate values associated with the first
2235    PointSet (this will cause internal memory to be allocated if it is not
2236    yet associated with coordinate values). */
2237       } else {
2238          ptr1 = astGetPoints( point1 );
2239 
2240 /* Allocate a temporary array to hold new pointer values. */
2241          ptr2 = (double **) astMalloc( sizeof( double * ) * (size_t) ncoord2 );
2242 
2243 /* Initialise this pointer array to point at the required subset of coordinate
2244    values. */
2245          if ( astOK ) {
2246             for ( i = 0; i < ncoord2; i++ ) {
2247                ptr2[ i ] = ptr1[ i + coord ] + point;
2248             }
2249 
2250 /* Associate the second PointSet with this new pointer array. This will free
2251    any internally allocated memory and replace any existing coordinate value
2252    association. */
2253             astSetPoints( point2, ptr2 );
2254 	 }
2255 
2256 /* Free the temporary pointer arry. */
2257          ptr2 = (double **) astFree( (void * ) ptr2 );
2258       }
2259    }
2260 }
2261 
TestAttrib(AstObject * this_object,const char * attrib,int * status)2262 static int TestAttrib( AstObject *this_object, const char *attrib, int *status ) {
2263 /*
2264 *  Name:
2265 *     TestAttrib
2266 
2267 *  Purpose:
2268 *     Test if a specified attribute value is set for a PointSet.
2269 
2270 *  Type:
2271 *     Private function.
2272 
2273 *  Synopsis:
2274 *     #include "pointset.h"
2275 *     int TestAttrib( AstObject *this, const char *attrib, int *status )
2276 
2277 *  Class Membership:
2278 *     PointSet member function (over-rides the astTestAttrib protected
2279 *     method inherited from the Object class).
2280 
2281 *  Description:
2282 *     This function returns a boolean result (0 or 1) to indicate
2283 *     whether a value has been set for one of a PointSet's attributes.
2284 
2285 *  Parameters:
2286 *     this
2287 *        Pointer to the PointSet.
2288 *     attrib
2289 *        Pointer to a null-terminated string specifying the attribute
2290 *        name.  This should be in lower case with no surrounding white
2291 *        space.
2292 *     status
2293 *        Pointer to the inherited status variable.
2294 
2295 *  Returned Value:
2296 *     One if a value has been set, otherwise zero.
2297 
2298 *  Notes:
2299 *     - A value of zero will be returned if this function is invoked
2300 *     with the global status set, or if it should fail for any reason.
2301 */
2302 
2303 /* Local Variables: */
2304    AstPointSet *this;            /* Pointer to the PointSet structure */
2305    int result;                   /* Result value to return */
2306 
2307 /* Initialise. */
2308    result = 0;
2309 
2310 /* Check the global error status. */
2311    if ( !astOK ) return result;
2312 
2313 /* Obtain a pointer to the PointSet structure. */
2314    this = (AstPointSet *) this_object;
2315 
2316 /* Check the attribute name and test the appropriate attribute. */
2317 
2318 /* Test if the name matches any of the read-only attributes of this
2319    class. If it does, then return zero. */
2320    if ( !strcmp( attrib, "ncoord" ) ||
2321         !strcmp( attrib, "npoint" ) ) {
2322       result = 0;
2323 
2324 /* If the attribute is still not recognised, pass it on to the parent
2325    method for further interpretation. */
2326    } else {
2327       result = (*parent_testattrib)( this_object, attrib, status );
2328    }
2329 
2330 /* Return the result, */
2331    return result;
2332 }
2333 
2334 /* Functions which access class attributes. */
2335 /* ---------------------------------------- */
2336 
2337 /*
2338 *att+
2339 *  Name:
2340 *     PointAccuracy
2341 
2342 *  Purpose:
2343 *     The absolute accuracies for all points in the PointSet.
2344 
2345 *  Type:
2346 *     Protected attribute.
2347 
2348 *  Synopsis:
2349 *     Floating point.
2350 
2351 *  Description:
2352 *     This attribute holds the absolute accuracy for each axis in the
2353 *     PointSet. It has a separate value for each axis. It is used when
2354 *     comparing two PointSets using the protected astEqual method inherited
2355 *     from the Object class. The default value for each axis is AST__BAD
2356 *     which causes the a default accuracy of each axis value to be calculated
2357 *     as  1.0E8*min( abs(axis value)*DBL_EPSILON, DBL_MIN ).
2358 
2359 *  Applicability:
2360 *     PointSet
2361 *        All PointSets have this attribute.
2362 *att-
2363 */
MAKE_CLEAR(PointAccuracy,acc,AST__BAD)2364 MAKE_CLEAR(PointAccuracy,acc,AST__BAD)
2365 MAKE_GET(PointAccuracy,double,AST__BAD,this->acc?this->acc[axis]:AST__BAD)
2366 MAKE_SET(PointAccuracy,double,acc,((value!=AST__BAD)?fabs(value):AST__BAD),AST__BAD)
2367 MAKE_TEST(PointAccuracy,(this->acc?this->acc[axis]!=AST__BAD:0))
2368 
2369 /* Copy constructor. */
2370 /* ----------------- */
2371 static void Copy( const AstObject *objin, AstObject *objout, int *status ) {
2372 /*
2373 *  Name:
2374 *     Copy
2375 
2376 *  Purpose:
2377 *     Copy constructor for PointSet objects.
2378 
2379 *  Type:
2380 *     Private function.
2381 
2382 *  Synopsis:
2383 *     void Copy( const AstObject *objin, AstObject *objout, int *status )
2384 
2385 *  Description:
2386 *     This function implements the copy constructor for PointSet objects.
2387 
2388 *  Parameters:
2389 *     objin
2390 *        Pointer to the object to be copied.
2391 *     objout
2392 *        Pointer to the object being constructed.
2393 *     status
2394 *        Pointer to the inherited status variable.
2395 
2396 *  Returned Value:
2397 *     void
2398 
2399 *  Notes:
2400 *     -  This constructor makes a deep copy, including a copy of the coordinate
2401 *     values (if any) associated with the input PointSet.
2402 */
2403 
2404 /* Local Variables: */
2405    AstPointSet *in;              /* Pointer to input PointSet */
2406    AstPointSet *out;             /* Pointer to output PointSet */
2407    int i;                        /* Loop counter for coordinates */
2408    int nval;                     /* Number of values to store */
2409 
2410 /* Check the global error status. */
2411    if ( !astOK ) return;
2412 
2413 /* Obtain pointers to the input and output PointSets. */
2414    in = (AstPointSet *) objin;
2415    out = (AstPointSet *) objout;
2416 
2417 /* For safety, first clear any references to the input coordinate values from
2418    the output PointSet. */
2419    out->ptr = NULL;
2420    out->values = NULL;
2421    out->acc = NULL;
2422 
2423 /* Copy axis accuracies. */
2424    if( in->acc ){
2425       out->acc = astStore( NULL, in->acc, sizeof( double )*(size_t) in->ncoord );
2426    }
2427 
2428 /* If the input PointSet is associated with coordinate values, we must
2429    allocate memory in the output PointSet to hold a copy of them. */
2430    if ( in->ptr ) {
2431 
2432 /* Determine the number of coordinate values to be stored and allocate memory
2433    to hold them, storing a pointer to this memory in the output PointSet. */
2434       nval = in->npoint * in->ncoord;
2435       out->values = (double *) astMalloc( sizeof( double ) * (size_t) nval );
2436 
2437 /* If OK, also allocate memory for the array of pointers into this values
2438    array, storing a pointer to this pointer array in the output PointSet. */
2439       if ( astOK ) {
2440          out->ptr = (double **) astMalloc( sizeof( double * )
2441                                            * (size_t) in->ncoord );
2442 
2443 /* If OK, initialise the new pointer array. */
2444          if ( astOK ) {
2445             for ( i = 0; i < in->ncoord; i++ ) {
2446                out->ptr[ i ] = out->values + ( i * in->npoint );
2447             }
2448 
2449 /* If we failed to allocate the pointer array, then free the values array. */
2450          } else {
2451             out->values = (double *) astFree( (void *) out->values );
2452          }
2453       }
2454 
2455 /* Copy the values for each coordinate from the input to the output. Use a
2456    memory copy to avoid floating point errors if the data are
2457    un-initialised. */
2458       if ( astOK ) {
2459          for ( i = 0; i < in->ncoord; i++ ) {
2460             (void) memcpy( (void *) out->ptr[ i ],
2461                            (const void *) in->ptr[ i ],
2462                            sizeof( double ) * (size_t) in->npoint );
2463          }
2464       }
2465    }
2466 }
2467 
2468 /* Destructor. */
2469 /* ----------- */
Delete(AstObject * obj,int * status)2470 static void Delete( AstObject *obj, int *status ) {
2471 /*
2472 *  Name:
2473 *     Delete
2474 
2475 *  Purpose:
2476 *     Destructor for PointSet objects.
2477 
2478 *  Type:
2479 *     Private function.
2480 
2481 *  Synopsis:
2482 *     void Delete( AstObject *obj, int *status )
2483 
2484 *  Description:
2485 *     This function implements the destructor for PointSet objects.
2486 
2487 *  Parameters:
2488 *     obj
2489 *        Pointer to the object to be deleted.
2490 *     status
2491 *        Pointer to the inherited status variable.
2492 
2493 *  Returned Value:
2494 *     void
2495 
2496 *  Notes:
2497 *     This function attempts to execute even if the global error status is
2498 *     set.
2499 */
2500 
2501 /* Local Variables: */
2502    AstPointSet *this;            /* Pointer to PointSet */
2503 
2504 /* Obtain a pointer to the PointSet structure. */
2505    this = (AstPointSet *) obj;
2506 
2507 /* Free memory holding axis accuracies. */
2508    this->acc = astFree( this->acc );
2509 
2510 /* Free any pointer array and associated coordinate values array, */
2511    this->ptr = (double **) astFree( (void *) this->ptr );
2512    this->values = (double *) astFree( (void *) this->values );
2513 
2514 /* Clear the remaining PointSet variables. */
2515    this->npoint = 0;
2516    this->ncoord = 0;
2517 }
2518 
2519 /* Dump function. */
2520 /* -------------- */
Dump(AstObject * this_object,AstChannel * channel,int * status)2521 static void Dump( AstObject *this_object, AstChannel *channel, int *status ) {
2522 /*
2523 *  Name:
2524 *     Dump
2525 
2526 *  Purpose:
2527 *     Dump function for PointSet objects.
2528 
2529 *  Type:
2530 *     Private function.
2531 
2532 *  Synopsis:
2533 *     void Dump( AstObject *this, AstChannel *channel, int *status )
2534 
2535 *  Description:
2536 *     This function implements the Dump function which writes out data
2537 *     for the PointSet class to an output Channel.
2538 
2539 *  Parameters:
2540 *     this
2541 *        Pointer to the PointSet whose data are being written.
2542 *     channel
2543 *        Pointer to the Channel to which the data are being written.
2544 *     status
2545 *        Pointer to the inherited status variable.
2546 
2547 *  Notes:
2548 *     - It is not recommended that PointSets containing large numbers
2549 *     of points be written out, as the coordinate data will be
2550 *     formatted as text and this will not be very efficient.
2551 */
2552 
2553 /* Local Constants: */
2554 #define KEY_LEN 50               /* Maximum length of a keyword */
2555 
2556 /* Local Variables: */
2557    AstPointSet *this;            /* Pointer to the PointSet structure */
2558    char key[ KEY_LEN + 1 ];      /* Buffer for keywords */
2559    int coord;                    /* Loop counter for coordinates */
2560    int i;                        /* Counter for coordinate values */
2561    int ival;                     /* Integer value */
2562    int makeComment;              /* Include a comment? */
2563    int point;                    /* Loop counter for points */
2564    int set;                      /* Attribute value set? */
2565 
2566 /* Check the global error status. */
2567    if ( !astOK ) return;
2568 
2569 /* Obtain a pointer to the PointSet structure. */
2570    this = (AstPointSet *) this_object;
2571 
2572 /* Write out values representing the instance variables for the
2573    PointSet class.  Accompany these with appropriate comment strings,
2574    possibly depending on the values being written.*/
2575 
2576 /* In the case of attributes, we first use the appropriate (private)
2577    Test...  member function to see if they are set. If so, we then use
2578    the (private) Get... function to obtain the value to be written
2579    out.
2580 
2581    For attributes which are not set, we use the astGet... method to
2582    obtain the value instead. This will supply a default value
2583    (possibly provided by a derived class which over-rides this method)
2584    which is more useful to a human reader as it corresponds to the
2585    actual default attribute value.  Since "set" will be zero, these
2586    values are for information only and will not be read back. */
2587 
2588 /* Npoint. */
2589 /* ------- */
2590    astWriteInt( channel, "Npoint", 1, 1, this->npoint,
2591                 "Number of points" );
2592 
2593 /* Ncoord. */
2594 /* ------- */
2595    astWriteInt( channel, "Ncoord", 1, 1, this->ncoord,
2596                 "Number of coordinates per point" );
2597 
2598 /* Axis Acuracies. */
2599 /* --------------- */
2600    for ( coord = 0; coord < this->ncoord; coord++ ) {
2601       if( astTestPointAccuracy( this, coord ) ) {
2602          (void) sprintf( key, "Acc%d", coord + 1 );
2603          astWriteDouble( channel, key, 1, 1, astGetPointAccuracy( this, coord ),
2604                          (coord == 0 ) ? "Axis accuracies..." : "" );
2605       }
2606    }
2607 
2608 /* Coordinate data. */
2609 /* ---------------- */
2610 /* Write an "Empty" value to indicate whether or not the PointSet
2611    contains data. */
2612    ival = ( this->ptr == NULL );
2613    set = ( ival != 0 );
2614    astWriteInt( channel, "Empty", set, 0, ival,
2615                 ival ? "PointSet is empty" :
2616                        "PointSet contains data" );
2617 
2618 /* If it contains data, create a suitable keyword for each coordinate
2619    value in turn. */
2620    if ( this->ptr ) {
2621       makeComment = 1;
2622       i = 0;
2623       for ( point = 0; point < this->npoint; point++ ) {
2624          for ( coord = 0; coord < this->ncoord; coord++ ) {
2625             i++;
2626             (void) sprintf( key, "X%d", i );
2627 
2628 /* Write the value out if good. Only supply a comment for the first good value. */
2629             if( this->ptr[ coord ][ point ] != AST__BAD ) {
2630                astWriteDouble( channel, key, 1, 1, this->ptr[ coord ][ point ],
2631                                ( makeComment ) ? "Coordinate values..." : "" );
2632                makeComment = 0;
2633             }
2634          }
2635       }
2636    }
2637 
2638 /* Undefine macros local to this function. */
2639 #undef KEY_LEN
2640 }
2641 
2642 /* Standard class functions. */
2643 /* ========================= */
2644 /* Implement the astIsAPointSet and astCheckPointSet functions using the macros
2645    defined for this purpose in the "object.h" header file. */
astMAKE_ISA(PointSet,Object)2646 astMAKE_ISA(PointSet,Object)
2647 astMAKE_CHECK(PointSet)
2648 
2649 AstPointSet *astPointSet_( int npoint, int ncoord, const char *options, int *status, ...) {
2650 /*
2651 *+
2652 *  Name:
2653 *     astPointSet
2654 
2655 *  Purpose:
2656 *     Create a PointSet.
2657 
2658 *  Type:
2659 *     Protected function.
2660 
2661 *  Synopsis:
2662 *     #include "pointset.h"
2663 *     AstPointSet *astPointSet( int npoint, int ncoord,
2664 *                               const char *options, ..., int *status )
2665 
2666 *  Class Membership:
2667 *     PointSet constructor.
2668 
2669 *  Description:
2670 *     This function creates a new PointSet and optionally initialises its
2671 *     attributes.
2672 
2673 *  Parameters:
2674 *     npoint
2675 *        The number of points to be stored in the PointSet (must be at
2676 *        least 1).
2677 *     ncoord
2678 *        The number of coordinate values associated with each point
2679 *        (must be at least 1).
2680 *     options
2681 *        Pointer to a null terminated string containing an optional
2682 *        comma-separated list of attribute assignments to be used for
2683 *        initialising the new PointSet. The syntax used is the same as
2684 *        for the astSet method and may include "printf" format
2685 *        specifiers identified by "%" symbols in the normal way.
2686 *     status
2687 *        Pointer to the inherited status variable.
2688 *     ...
2689 *        If the "options" string contains "%" format specifiers, then
2690 *        an optional list of arguments may follow it in order to
2691 *        supply values to be substituted for these specifiers. The
2692 *        rules for supplying these are identical to those for the
2693 *        astSet method (and for the C "printf" function).
2694 
2695 *  Returned Value:
2696 *     A pointer to the new PointSet.
2697 
2698 *  Notes:
2699 *     - A NULL pointer will be returned if this function is invoked
2700 *     with the global error status set, or if it should fail for any
2701 *     reason.
2702 *-
2703 */
2704 
2705 /* Local Variables: */
2706    astDECLARE_GLOBALS            /* Pointer to thread-specific global data */
2707    AstPointSet *new;             /* Pointer to new PointSet */
2708    va_list args;                 /* Variable argument list */
2709 
2710 /* Get a pointer to the thread specific global data structure. */
2711    astGET_GLOBALS(NULL);
2712 
2713 /* Check the global status. */
2714    if ( !astOK ) return NULL;
2715 
2716 /* Initialise the PointSet, allocating memory and initialising the
2717    virtual function table as well if necessary. */
2718    new = astInitPointSet( NULL, sizeof( AstPointSet ), !class_init,
2719                           &class_vtab, "PointSet", npoint, ncoord );
2720 
2721 /* If successful, note that the virtual function table has been
2722    initialised. */
2723    if ( astOK ) {
2724       class_init = 1;
2725 
2726 /* Obtain the variable argument list and pass it along with the
2727    options string to the astVSet method to initialise the new
2728    PointSet's attributes. */
2729       va_start( args, status );
2730       astVSet( new, options, NULL, args );
2731       va_end( args );
2732 
2733 /* If an error occurred, clean up by deleting the new object. */
2734       if ( !astOK ) new = astDelete( new );
2735    }
2736 
2737 /* Return a pointer to the new PointSet. */
2738    return new;
2739 }
2740 
astPointSetId_(int npoint,int ncoord,const char * options,int * status,...)2741 AstPointSet *astPointSetId_( int npoint, int ncoord,
2742                              const char *options, int *status, ...) {
2743 /*
2744 *  Name:
2745 *     astPointSetId_
2746 
2747 *  Purpose:
2748 *     Create a PointSet.
2749 
2750 *  Type:
2751 *     Private function.
2752 
2753 *  Synopsis:
2754 *     #include "pointset.h"
2755 *     AstPointSet *astPointSetId_( int npoint, int ncoord,
2756 *                                  const char *options, ... )
2757 
2758 *  Class Membership:
2759 *     PointSet constructor.
2760 
2761 *  Description:
2762 *     This function implements the external (public) interface to the
2763 *     astPointSet constructor function. It returns an ID value
2764 *     (instead of a true C pointer) to external users, and must be
2765 *     provided because astPointSet_ has a variable argument list which
2766 *     cannot be encapsulated in a macro (where this conversion would
2767 *     otherwise occur).
2768 *
2769 *     The variable argument list also prevents this function from
2770 *     invoking astPointSet_ directly, so it must be a
2771 *     re-implementation of it in all respects, except for the final
2772 *     conversion of the result to an ID value.
2773 
2774 *  Parameters:
2775 *     As for astPointSet_.
2776 
2777 *  Returned Value:
2778 *     The ID value associated with the new PointSet.
2779 */
2780 
2781 /* Local Variables: */
2782    astDECLARE_GLOBALS            /* Pointer to thread-specific global data */
2783    AstPointSet *new;             /* Pointer to new PointSet */
2784    va_list args;                 /* Variable argument list */
2785 
2786 /* Check the global status. */
2787    if ( !astOK ) return NULL;
2788 
2789 /* Get a pointer to the thread specific global data structure. */
2790    astGET_GLOBALS(NULL);
2791 
2792 /* Initialise the PointSet, allocating memory and initialising the
2793    virtual function table as well if necessary. */
2794    new = astInitPointSet( NULL, sizeof( AstPointSet ), !class_init,
2795                           &class_vtab, "PointSet", npoint, ncoord );
2796 
2797 /* If successful, note that the virtual function table has been
2798    initialised. */
2799    if ( astOK ) {
2800       class_init = 1;
2801 
2802 /* Obtain the variable argument list and pass it along with the
2803    options string to the astVSet method to initialise the new
2804    PointSet's attributes. */
2805       va_start( args, status );
2806       astVSet( new, options, NULL, args );
2807       va_end( args );
2808 
2809 /* If an error occurred, clean up by deleting the new object. */
2810       if ( !astOK ) new = astDelete( new );
2811    }
2812 
2813 /* Return an ID value for the new PointSet. */
2814    return astMakeId( new );
2815 }
2816 
astInitPointSet_(void * mem,size_t size,int init,AstPointSetVtab * vtab,const char * name,int npoint,int ncoord,int * status)2817 AstPointSet *astInitPointSet_( void *mem, size_t size, int init,
2818                                AstPointSetVtab *vtab, const char *name,
2819                                int npoint, int ncoord, int *status ) {
2820 /*
2821 *+
2822 *  Name:
2823 *     astInitPointSet
2824 
2825 *  Purpose:
2826 *     Initialise a PointSet.
2827 
2828 *  Type:
2829 *     Protected function.
2830 
2831 *  Synopsis:
2832 *     #include "pointset.h"
2833 *     AstPointSet *astInitPointSet( void *mem, size_t size, int init,
2834 *                                   AstPointSetVtab *vtab, const char *name,
2835 *                                   int npoint, int ncoord )
2836 
2837 *  Class Membership:
2838 *     PointSet initialiser.
2839 
2840 *  Description:
2841 *     This function is provided for use by class implementations to initialise
2842 *     a new PointSet object. It allocates memory (if necessary) to accommodate
2843 *     the PointSet plus any additional data associated with the derived class.
2844 *     It then initialises a PointSet structure at the start of this memory. If
2845 *     the "init" flag is set, it also initialises the contents of a virtual
2846 *     function table for a PointSet at the start of the memory passed via the
2847 *     "vtab" parameter.
2848 
2849 *  Parameters:
2850 *     mem
2851 *        A pointer to the memory in which the PointSet is to be created. This
2852 *        must be of sufficient size to accommodate the PointSet data
2853 *        (sizeof(PointSet)) plus any data used by the derived class. If a value
2854 *        of NULL is given, this function will allocate the memory itself using
2855 *        the "size" parameter to determine its size.
2856 *     size
2857 *        The amount of memory used by the PointSet (plus derived class data).
2858 *        This will be used to allocate memory if a value of NULL is given for
2859 *        the "mem" parameter. This value is also stored in the PointSet
2860 *        structure, so a valid value must be supplied even if not required for
2861 *        allocating memory.
2862 *     init
2863 *        A logical flag indicating if the PointSet's virtual function table is
2864 *        to be initialised. If this value is non-zero, the virtual function
2865 *        table will be initialised by this function.
2866 *     vtab
2867 *        Pointer to the start of the virtual function table to be associated
2868 *        with the new PointSet.
2869 *     name
2870 *        Pointer to a constant null-terminated character string which contains
2871 *        the name of the class to which the new object belongs (it is this
2872 *        pointer value that will subsequently be returned by the Object
2873 *        astClass function).
2874 *     npoint
2875 *        The number of points in the PointSet (must be at least 1).
2876 *     ncoord
2877 *        The number of coordinate values associated with each point (must be
2878 *        at least 1).
2879 
2880 *  Returned Value:
2881 *     A pointer to the new PointSet.
2882 
2883 *  Notes:
2884 *     -  A null pointer will be returned if this function is invoked with the
2885 *     global error status set, or if it should fail for any reason.
2886 *-
2887 */
2888 
2889 /* Local Variables: */
2890    AstPointSet *new;             /* Pointer to new PointSet */
2891 
2892 /* Check the global status. */
2893    if ( !astOK ) return NULL;
2894 
2895 /* If necessary, initialise the virtual function table. */
2896    if ( init ) astInitPointSetVtab( vtab, name );
2897 
2898 /* Initialise. */
2899    new = NULL;
2900 
2901 /* Check the initialisation values for validity, reporting an error if
2902    necessary. */
2903    if ( npoint < 1 ) {
2904       astError( AST__NPTIN, "astInitPointSet(%s): Number of points (%d) is "
2905                 "not valid.", status, name, npoint );
2906    } else if ( ncoord < 1 ) {
2907       astError( AST__NCOIN, "astInitPointSet(%s): Number of coordinates per "
2908                 "point (%d) is not valid.", status, name, ncoord );
2909    }
2910 
2911 /* Initialise an Object structure (the parent class) as the first component
2912    within the PointSet structure, allocating memory if necessary. */
2913    new = (AstPointSet *) astInitObject( mem, size, 0,
2914                                         (AstObjectVtab *) vtab, name );
2915 
2916    if ( astOK ) {
2917 
2918 /* Initialise the PointSet data. */
2919 /* ----------------------------- */
2920 /* Store the number of points and number of coordinate values per point. */
2921       new->npoint = npoint;
2922       new->ncoord = ncoord;
2923 
2924 /* Initialise pointers to the pointer array and associated coordinate
2925    values array. */
2926       new->ptr = NULL;
2927       new->values = NULL;
2928       new->acc = NULL;
2929 
2930 /* If an error occurred, clean up by deleting the new object. */
2931       if ( !astOK ) new = astDelete( new );
2932    }
2933 
2934 /* Return a pointer to the new object. */
2935    return new;
2936 }
2937 
astLoadPointSet_(void * mem,size_t size,AstPointSetVtab * vtab,const char * name,AstChannel * channel,int * status)2938 AstPointSet *astLoadPointSet_( void *mem, size_t size,
2939                                AstPointSetVtab *vtab, const char *name,
2940                                AstChannel *channel, int *status ) {
2941 /*
2942 *+
2943 *  Name:
2944 *     astLoadPointSet
2945 
2946 *  Purpose:
2947 *     Load a PointSet.
2948 
2949 *  Type:
2950 *     Protected function.
2951 
2952 *  Synopsis:
2953 *     #include "pointset.h"
2954 *     AstPointSet *astLoadPointSet( void *mem, size_t size,
2955 *                                   AstPointSetVtab *vtab, const char *name,
2956 *                                   AstChannel *channel )
2957 
2958 *  Class Membership:
2959 *     PointSet loader.
2960 
2961 *  Description:
2962 *     This function is provided to load a new PointSet using data read
2963 *     from a Channel. It first loads the data used by the parent class
2964 *     (which allocates memory if necessary) and then initialises a
2965 *     PointSet structure in this memory, using data read from the
2966 *     input Channel.
2967 
2968 *  Parameters:
2969 *     mem
2970 *        A pointer to the memory into which the PointSet is to be
2971 *        loaded.  This must be of sufficient size to accommodate the
2972 *        PointSet data (sizeof(PointSet)) plus any data used by
2973 *        derived classes. If a value of NULL is given, this function
2974 *        will allocate the memory itself using the "size" parameter to
2975 *        determine its size.
2976 *     size
2977 *        The amount of memory used by the PointSet (plus derived class
2978 *        data).  This will be used to allocate memory if a value of
2979 *        NULL is given for the "mem" parameter. This value is also
2980 *        stored in the PointSet structure, so a valid value must be
2981 *        supplied even if not required for allocating memory.
2982 *
2983 *        If the "vtab" parameter is NULL, the "size" value is ignored
2984 *        and sizeof(AstPointSet) is used instead.
2985 *     vtab
2986 *        Pointer to the start of the virtual function table to be
2987 *        associated with the new PointSet. If this is NULL, a pointer
2988 *        to the (static) virtual function table for the PointSet class
2989 *        is used instead.
2990 *     name
2991 *        Pointer to a constant null-terminated character string which
2992 *        contains the name of the class to which the new object
2993 *        belongs (it is this pointer value that will subsequently be
2994 *        returned by the astGetClass method).
2995 *
2996 *        If the "vtab" parameter is NULL, the "name" value is ignored
2997 *        and a pointer to the string "PointSet" is used instead.
2998 
2999 *  Returned Value:
3000 * A pointer to the new PointSet.
3001 
3002 *  Notes:
3003 *     - A null pointer will be returned if this function is invoked
3004 *     with the global error status set, or if it should fail for any
3005 *     reason.
3006 *-
3007 */
3008 
3009 /* Local Constants: */
3010    astDECLARE_GLOBALS            /* Pointer to thread-specific global data */
3011 #define KEY_LEN 50               /* Maximum length of a keyword */
3012 
3013 /* Local Variables: */
3014    AstPointSet *new;             /* Pointer to the new PointSet */
3015    char key[ KEY_LEN + 1 ];      /* Buffer for keywords */
3016    double acc;                   /* Accuracy value */
3017    int coord;                    /* Loop counter for coordinates */
3018    int empty;                    /* PointSet empty? */
3019    int i;                        /* Counter for coordinate values */
3020    int point;                    /* Loop counter for points */
3021 
3022 /* Get a pointer to the thread specific global data structure. */
3023    astGET_GLOBALS(channel);
3024 
3025 /* Initialise. */
3026    new = NULL;
3027 
3028 /* Check the global error status. */
3029    if ( !astOK ) return new;
3030 
3031 /* If a NULL virtual function table has been supplied, then this is
3032    the first loader to be invoked for this PointSet. In this case the
3033    PointSet belongs to this class, so supply appropriate values to be
3034    passed to the parent class loader (and its parent, etc.). */
3035    if ( !vtab ) {
3036       size = sizeof( AstPointSet );
3037       vtab = &class_vtab;
3038       name = "PointSet";
3039 
3040 /* If required, initialise the virtual function table for this class. */
3041       if ( !class_init ) {
3042          astInitPointSetVtab( vtab, name );
3043          class_init = 1;
3044       }
3045    }
3046 
3047 /* Invoke the parent class loader to load data for all the ancestral
3048    classes of the current one, returning a pointer to the resulting
3049    partly-built PointSet. */
3050    new = astLoadObject( mem, size, (AstObjectVtab *) vtab, name,
3051                         channel );
3052 
3053    if ( astOK ) {
3054 
3055 /* Initialise the PointSet's data pointers. */
3056       new->ptr = NULL;
3057       new->values = NULL;
3058 
3059 /* Read input data. */
3060 /* ================ */
3061 /* Request the input Channel to read all the input data appropriate to
3062    this class into the internal "values list". */
3063       astReadClassData( channel, "PointSet" );
3064 
3065 /* Now read each individual data item from this list and use it to
3066    initialise the appropriate instance variable(s) for this class. */
3067 
3068 /* In the case of attributes, we first read the "raw" input value,
3069    supplying the "unset" value as the default. If a "set" value is
3070    obtained, we then use the appropriate (private) Set... member
3071    function to validate and set the value properly. */
3072 
3073 /* Npoint. */
3074 /* ------- */
3075       new->npoint = astReadInt( channel, "npoint", 1 );
3076       if ( new->npoint < 1 ) new->npoint = 1;
3077 
3078 /* Ncoord. */
3079 /* ------- */
3080       new->ncoord = astReadInt( channel, "ncoord", 1 );
3081       if ( new->ncoord < 1 ) new->ncoord = 1;
3082 
3083 /* Axis Acuracies. */
3084 /* --------------- */
3085       new->acc = NULL;
3086       for ( coord = 0; coord < new->ncoord; coord++ ) {
3087          (void) sprintf( key, "acc%d", coord + 1 );
3088          acc = astReadDouble( channel, key, AST__BAD );
3089          if( !new->acc && acc != AST__BAD ) {
3090             new->acc = astMalloc( sizeof( double )*(size_t) new->ncoord );
3091             if( new->acc ) {
3092                for( i = 0; i < coord - 1; i++ ) new->acc[ i ] = AST__BAD;
3093             }
3094          }
3095          if( new->acc ) new->acc[ coord ] = acc;
3096       }
3097 
3098 /* Coordinate data. */
3099 /* ---------------- */
3100 /* Read a value for the "Empty" keyword to see whether the PointSet
3101    contains data. */
3102       empty = astReadInt( channel, "empty", 0 );
3103 
3104 /* If it does, allocate memory to hold the coordinate data and
3105    pointers. */
3106       if ( astOK && !empty ) {
3107          new->ptr = astMalloc( sizeof( double * ) * (size_t) new->ncoord );
3108          new->values = astMalloc( sizeof( double ) *
3109                                   (size_t) ( new->npoint * new->ncoord ) );
3110          if ( astOK ) {
3111 
3112 /* Initialise the array of pointers into the main data array. */
3113             for ( coord = 0; coord < new->ncoord; coord++ ) {
3114                new->ptr[ coord ] = new->values + ( coord * new->npoint );
3115             }
3116 
3117 /* Create a keyword for each coordinate value to be read. */
3118             i = 0;
3119             for ( point = 0; point < new->npoint; point++ ) {
3120                for ( coord = 0; coord < new->ncoord; coord++ ) {
3121                   i++;
3122                   (void) sprintf( key, "x%d", i );
3123 
3124 /* Read and assign the values. */
3125                   new->ptr[ coord ][ point ] =
3126                      astReadDouble( channel, key, AST__BAD );
3127                }
3128             }
3129          }
3130 
3131 /* If an error occurred, clean up by freeing the memory allocated
3132    above, thus emptying the PointSet. */
3133          if ( !astOK ) {
3134             new->ptr = astFree( new->ptr );
3135             new->values = astFree( new->values );
3136          }
3137       }
3138 
3139 /* If an error occurred, clean up by deleting the new PointSet. */
3140       if ( !astOK ) new = astDelete( new );
3141    }
3142 
3143 /* Return the new PointSet pointer. */
3144    return new;
3145 
3146 /* Undefine macros local to this function. */
3147 #undef KEY_LEN
3148 }
3149 
3150 /* Virtual function interfaces. */
3151 /* ============================ */
3152 /* These provide the external interface to the virtual functions defined by
3153    this class. Each simply checks the global error status and then locates and
3154    executes the appropriate member function, using the function pointer stored
3155    in the object's virtual function table (this pointer is located using the
3156    astMEMBER macro defined in "object.h").
3157 
3158    Note that the member function may not be the one defined here, as it may
3159    have been over-ridden by a derived class. However, it should still have the
3160    same interface. */
astGetNpoint_(const AstPointSet * this,int * status)3161 int astGetNpoint_( const AstPointSet *this, int *status ) {
3162    if ( !astOK ) return 0;
3163    return (**astMEMBER(this,PointSet,GetNpoint))( this, status );
3164 }
astGetNcoord_(const AstPointSet * this,int * status)3165 int astGetNcoord_( const AstPointSet *this, int *status ) {
3166    if ( !astOK ) return 0;
3167    return (**astMEMBER(this,PointSet,GetNcoord))( this, status );
3168 }
astGetPoints_(AstPointSet * this,int * status)3169 double **astGetPoints_( AstPointSet *this, int *status ) {
3170    if ( !astOK ) return NULL;
3171    return (**astMEMBER(this,PointSet,GetPoints))( this, status );
3172 }
astPermPoints_(AstPointSet * this,int forward,const int perm[],int * status)3173 void astPermPoints_( AstPointSet *this, int forward, const int perm[], int *status ) {
3174    if ( !astOK ) return;
3175    (**astMEMBER(this,PointSet,PermPoints))( this, forward, perm, status );
3176 }
astSetPoints_(AstPointSet * this,double ** ptr,int * status)3177 void astSetPoints_( AstPointSet *this, double **ptr, int *status ) {
3178    if ( !astOK ) return;
3179    (**astMEMBER(this,PointSet,SetPoints))( this, ptr, status );
3180 }
astSetNpoint_(AstPointSet * this,int npoint,int * status)3181 void astSetNpoint_( AstPointSet *this, int npoint, int *status ) {
3182    if ( !astOK ) return;
3183    (**astMEMBER(this,PointSet,SetNpoint))( this, npoint, status );
3184 }
astSetSubPoints_(AstPointSet * point1,int point,int coord,AstPointSet * point2,int * status)3185 void astSetSubPoints_( AstPointSet *point1, int point, int coord,
3186                        AstPointSet *point2, int *status ) {
3187    if ( !astOK ) return;
3188    (**astMEMBER(point1,PointSet,SetSubPoints))( point1, point, coord, point2, status );
3189 }
astAppendPoints_(AstPointSet * this,AstPointSet * that,int * status)3190 AstPointSet *astAppendPoints_( AstPointSet *this, AstPointSet *that, int *status ) {
3191    if ( !astOK ) return NULL;
3192    return (**astMEMBER(this,PointSet,AppendPoints))( this, that, status );
3193 }
astBndPoints_(AstPointSet * this,double * lbnd,double * ubnd,int * status)3194 void astBndPoints_( AstPointSet *this, double *lbnd, double *ubnd, int *status ) {
3195    if ( !astOK ) return;
3196    (**astMEMBER(this,PointSet,BndPoints))( this, lbnd, ubnd, status );
3197 }
3198 
astReplaceNaN_(AstPointSet * this,int * status)3199 int astReplaceNaN_( AstPointSet *this, int *status ) {
3200    if ( !astOK ) return 0;
3201    return (**astMEMBER(this,PointSet,ReplaceNaN))( this, status );
3202 }
3203 
3204 
3205 
3206 
3207 
3208