1 /*
2 *class++
3 * Name:
4 * Polygon
5
6 * Purpose:
7 * A polygonal region within a 2-dimensional Frame.
8
9 * Constructor Function:
10 c astPolygon
11 f AST_POLYGON
12
13 * Description:
14 * The Polygon class implements a polygonal area, defined by a
15 * collection of vertices, within a 2-dimensional Frame. The vertices
16 * are connected together by geodesic curves within the encapsulated Frame.
17 * For instance, if the encapsulated Frame is a simple Frame then the
18 * geodesics will be straight lines, but if the Frame is a SkyFrame then
19 * the geodesics will be great circles. Note, the vertices must be
20 * supplied in an order such that the inside of the polygon is to the
21 * left of the boundary as the vertices are traversed. Supplying them
22 * in the reverse order will effectively negate the polygon.
23 *
24 * Within a SkyFrame, neighbouring vertices are always joined using the
25 * shortest path. Thus if an edge of 180 degrees or more in length is
26 * required, it should be split into section each of which is less
27 * than 180 degrees. The closed path joining all the vertices in order
28 * will divide the celestial sphere into two disjoint regions. The
29 * inside of the polygon is the region which is circled in an
30 * anti-clockwise manner (when viewed from the inside of the celestial
31 * sphere) when moving through the list of vertices in the order in
32 * which they were supplied when the Polygon was created (i.e. the
33 * inside is to the left of the boundary when moving through the
34 * vertices in the order supplied).
35
36 * Inheritance:
37 * The Polygon class inherits from the Region class.
38
39 * Attributes:
40 * In addition to those attributes common to all Regions, every
41 * Polygon also has the following attributes:
42 * - SimpVertices: Simplify by transforming the vertices?
43
44 * Functions:
45 c In addition to those functions applicable to all Regions, the
46 c following functions may also be applied to all Polygons:
47 f In addition to those routines applicable to all Regions, the
48 f following routines may also be applied to all Polygons:
49 *
50 c - astDownsize: Reduce the number of vertices in a Polygon.
51 f - AST_DOWNSIZE: Reduce the number of vertices in a Polygon.
52 c - astConvex<X>: Create a Polygon giving the convex hull of a pixel array
53 f - AST_CONVEX<X>: Create a Polygon giving the convex hull of a pixel array
54 c - astOutline<X>: Create a Polygon outlining values in a pixel array
55 f - AST_OUTLINE<X>: Create a Polygon outlining values in a pixel array
56
57 * Copyright:
58 * Copyright (C) 1997-2006 Council for the Central Laboratory of the
59 * Research Councils
60 * Copyright (C) 2009 Science & Technology Facilities Council.
61 * All Rights Reserved.
62
63 * Licence:
64 * This program is free software: you can redistribute it and/or
65 * modify it under the terms of the GNU Lesser General Public
66 * License as published by the Free Software Foundation, either
67 * version 3 of the License, or (at your option) any later
68 * version.
69 *
70 * This program is distributed in the hope that it will be useful,
71 * but WITHOUT ANY WARRANTY; without even the implied warranty of
72 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
73 * GNU Lesser General Public License for more details.
74 *
75 * You should have received a copy of the GNU Lesser General
76 * License along with this program. If not, see
77 * <http://www.gnu.org/licenses/>.
78
79 * Authors:
80 * DSB: David S. Berry (Starlink)
81
82 * History:
83 * 26-OCT-2004 (DSB):
84 * Original version.
85 * 28-MAY-2009 (DSB):
86 * Added astDownsize.
87 * 29-MAY-2009 (DSB):
88 * Added astOutline<X>.
89 * 30-JUN-2009 (DSB):
90 * Override astGetBounded.
91 * 4-NOV-2013 (DSB):
92 * Modify RegPins so that it can handle uncertainty regions that straddle
93 * a discontinuity. Previously, such uncertainty Regions could have a huge
94 * bounding box resulting in matching region being far too big.
95 * 6-DEC-2013 (DSB):
96 * Reverse the order of the vertices when the Polygon is created,
97 * if necessary, to ensure that the unnegated Polygon is bounded.
98 * The parent Region class assumes that unnegated regions are
99 * bounded.
100 * 6-JAN-2014 (DSB):
101 * Free edges when clearing the cache, not when establishing a new
102 * cache, as the number of edges may have changed.
103 * 10-JAN-2014 (DSB):
104 * - Remove unused parameter description in prologue of for astOutline<X>
105 * 24-FEB-2014 (DSB):
106 * Added astConvex<X>.
107 * 25-FEB-2014 (DSB):
108 * Added attribute SimpVertices.
109 *class--
110 */
111
112 /* Module Macros. */
113 /* ============== */
114 /* Set the name of the class we are implementing. This indicates to
115 the header files that define class interfaces that they should make
116 "protected" symbols available. */
117 #define astCLASS Polygon
118
119 /* Macros which return the maximum and minimum of two values. */
120 #define MAX(aa,bb) ((aa)>(bb)?(aa):(bb))
121 #define MIN(aa,bb) ((aa)<(bb)?(aa):(bb))
122
123 /* Macro to check for equality of floating point values. We cannot
124 compare bad values directory because of the danger of floating point
125 exceptions, so bad values are dealt with explicitly. */
126 #define EQUAL(aa,bb) (((aa)==AST__BAD)?(((bb)==AST__BAD)?1:0):(((bb)==AST__BAD)?0:(fabs((aa)-(bb))<=1.0E9*MAX((fabs(aa)+fabs(bb))*DBL_EPSILON,DBL_MIN))))
127
128 /* Define a macro for testing if a pixel value <V> satisfies the requirements
129 specified by <Oper> and <Value>. Compiler optimisation should remove
130 all the "if" testing from this expression. */
131 #define ISVALID(V,OperI,Value) ( \
132 ( OperI == AST__LT ) ? ( (V) < Value ) : ( \
133 ( OperI == AST__LE ) ? ( (V) <= Value ) : ( \
134 ( OperI == AST__EQ ) ? ( (V) == Value ) : ( \
135 ( OperI == AST__GE ) ? ( (V) >= Value ) : ( \
136 ( OperI == AST__NE ) ? ( (V) != Value ) : ( \
137 (V) > Value \
138 ) \
139 ) \
140 ) \
141 ) \
142 ) \
143 )
144
145 /* Macros specifying whether a point is inside, outside or on the
146 boundary of the polygon. */
147 #define UNKNOWN 0
148 #define IN 1
149 #define OUT 2
150 #define ON 3
151
152 /* Include files. */
153 /* ============== */
154 /* Interface definitions. */
155 /* ---------------------- */
156
157 #include "globals.h" /* Thread-safe global data access */
158 #include "error.h" /* Error reporting facilities */
159 #include "memory.h" /* Memory allocation facilities */
160 #include "object.h" /* Base Object class */
161 #include "pointset.h" /* Sets of points/coordinates */
162 #include "region.h" /* Coordinate regions (parent class) */
163 #include "channel.h" /* I/O channels */
164 #include "box.h" /* Box Regions */
165 #include "wcsmap.h" /* Definitons of AST__DPI etc */
166 #include "polygon.h" /* Interface definition for this class */
167 #include "mapping.h" /* Position mappings */
168 #include "unitmap.h" /* Unit Mapping */
169 #include "pal.h" /* SLALIB library interface */
170 #include "frame.h" /* Coordinate system description */
171
172 /* Error code definitions. */
173 /* ----------------------- */
174 #include "ast_err.h" /* AST error codes */
175
176 /* C header files. */
177 /* --------------- */
178 #include <float.h>
179 #include <math.h>
180 #include <stdarg.h>
181 #include <stdlib.h>
182 #include <stddef.h>
183 #include <stdio.h>
184 #include <string.h>
185 #include <limits.h>
186
187 /* Type definitions. */
188 /* ================= */
189
190 /* A structure that holds information about an edge of the new Polygon
191 being created by astDownsize. The edge is a line betweeen two of the
192 vertices of the original Polygon. */
193 typedef struct Segment {
194 int i1; /* Index of starting vertex within old Polygon */
195 int i2; /* Index of ending vertex within old Polygon */
196 double error; /* Max geodesic distance from any old vertex to the line */
197 int imax; /* Index of the old vertex at which max error is reached */
198 struct Segment *next;/* Pointer to next Segment in a double link list */
199 struct Segment *prev;/* Pointer to previous Segment in a double link list */
200 } Segment;
201
202
203 /* Module Variables. */
204 /* ================= */
205
206 /* Address of this static variable is used as a unique identifier for
207 member of this class. */
208 static int class_check;
209
210 /* Pointers to parent class methods which are extended by this class. */
211 static AstPointSet *(* parent_transform)( AstMapping *, AstPointSet *, int, AstPointSet *, int * );
212 static AstMapping *(* parent_simplify)( AstMapping *, int * );
213 static void (* parent_setregfs)( AstRegion *, AstFrame *, int * );
214 static void (* parent_resetcache)( AstRegion *, int * );
215 static const char *(* parent_getattrib)( AstObject *, const char *, int * );
216 static int (* parent_testattrib)( AstObject *, const char *, int * );
217 static void (* parent_clearattrib)( AstObject *, const char *, int * );
218 static void (* parent_setattrib)( AstObject *, const char *, int * );
219
220
221 #ifdef THREAD_SAFE
222 /* Define how to initialise thread-specific globals. */
223 #define GLOBAL_inits \
224 globals->Class_Init = 0; \
225 globals->GetAttrib_Buff[ 0 ] = 0;
226
227 /* Create the function that initialises global data for this module. */
228 astMAKE_INITGLOBALS(Polygon)
229
230 /* Define macros for accessing each item of thread specific global data. */
231 #define class_init astGLOBAL(Polygon,Class_Init)
232 #define class_vtab astGLOBAL(Polygon,Class_Vtab)
233 #define getattrib_buff astGLOBAL(Polygon,GetAttrib_Buff)
234
235
236 #include <pthread.h>
237
238
239 #else
240
241 static char getattrib_buff[ 51 ];
242
243 /* Define the class virtual function table and its initialisation flag
244 as static variables. */
245 static AstPolygonVtab class_vtab; /* Virtual function table */
246 static int class_init = 0; /* Virtual function table initialised? */
247
248 #endif
249
250 /* External Interface Function Prototypes. */
251 /* ======================================= */
252 /* The following functions have public prototypes only (i.e. no
253 protected prototypes), so we must provide local prototypes for use
254 within this module. */
255 AstPolygon *astPolygonId_( void *, int, int, const double *, void *, const char *, ... );
256
257 /* Prototypes for Private Member Functions. */
258 /* ======================================== */
259
260 /* Define a macro that expands to a single prototype for function
261 FindInsidePoint for a given data type and operation. */
262 #define FINDINSIDEPOINT_PROTO0(X,Xtype,Oper) \
263 static void FindInsidePoint##Oper##X( Xtype, const Xtype *, const int[2], const int[2], int *, int *, int *, int * );
264
265 /* Define a macro that expands to a set of prototypes for all operations
266 for function FindInsidePoint for a given data type. */
267 #define FINDINSIDEPOINT_PROTO(X,Xtype) \
268 FINDINSIDEPOINT_PROTO0(X,Xtype,LT) \
269 FINDINSIDEPOINT_PROTO0(X,Xtype,LE) \
270 FINDINSIDEPOINT_PROTO0(X,Xtype,EQ) \
271 FINDINSIDEPOINT_PROTO0(X,Xtype,GE) \
272 FINDINSIDEPOINT_PROTO0(X,Xtype,GT) \
273 FINDINSIDEPOINT_PROTO0(X,Xtype,NE)
274
275 /* Use the above macros to define all FindInsidePoint prototypes for all
276 data types and operations. */
277 #if HAVE_LONG_DOUBLE /* Not normally implemented */
278 FINDINSIDEPOINT_PROTO(LD,long double)
279 #endif
280 FINDINSIDEPOINT_PROTO(D,double)
281 FINDINSIDEPOINT_PROTO(L,long int)
282 FINDINSIDEPOINT_PROTO(UL,unsigned long int)
283 FINDINSIDEPOINT_PROTO(I,int)
284 FINDINSIDEPOINT_PROTO(UI,unsigned int)
285 FINDINSIDEPOINT_PROTO(S,short int)
286 FINDINSIDEPOINT_PROTO(US,unsigned short int)
287 FINDINSIDEPOINT_PROTO(B,signed char)
288 FINDINSIDEPOINT_PROTO(UB,unsigned char)
289 FINDINSIDEPOINT_PROTO(F,float)
290
291 /* Define a macro that expands to a single prototype for function
292 TraceEdge for a given data type and operation. */
293 #define TRACEEDGE_PROTO0(X,Xtype,Oper) \
294 static AstPointSet *TraceEdge##Oper##X( Xtype, const Xtype *, const int[2], const int[2], int, int, int, int, int, int * );
295
296 /* Define a macro that expands to a set of prototypes for all operations
297 for function TraceEdge for a given data type. */
298 #define TRACEEDGE_PROTO(X,Xtype) \
299 TRACEEDGE_PROTO0(X,Xtype,LT) \
300 TRACEEDGE_PROTO0(X,Xtype,LE) \
301 TRACEEDGE_PROTO0(X,Xtype,EQ) \
302 TRACEEDGE_PROTO0(X,Xtype,GE) \
303 TRACEEDGE_PROTO0(X,Xtype,GT) \
304 TRACEEDGE_PROTO0(X,Xtype,NE)
305
306 /* Use the above macros to define all TraceEdge prototypes for all
307 data types and operations. */
308 #if HAVE_LONG_DOUBLE /* Not normally implemented */
309 TRACEEDGE_PROTO(LD,long double)
310 #endif
311 TRACEEDGE_PROTO(D,double)
312 TRACEEDGE_PROTO(L,long int)
313 TRACEEDGE_PROTO(UL,unsigned long int)
314 TRACEEDGE_PROTO(I,int)
315 TRACEEDGE_PROTO(UI,unsigned int)
316 TRACEEDGE_PROTO(S,short int)
317 TRACEEDGE_PROTO(US,unsigned short int)
318 TRACEEDGE_PROTO(B,signed char)
319 TRACEEDGE_PROTO(UB,unsigned char)
320 TRACEEDGE_PROTO(F,float)
321
322 /* Define a macro that expands to a single prototype for function
323 PartHull for a given data type and operation. */
324 #define PARTHULL_PROTO0(X,Xtype,Oper) \
325 static void PartHull##Oper##X( Xtype, const Xtype[], int, int, int, int, int, int, int, const int[2], double **, double **, int *, int * );
326
327 /* Define a macro that expands to a set of prototypes for all operations
328 for function PartHull for a given data type. */
329 #define PARTHULL_PROTO(X,Xtype) \
330 PARTHULL_PROTO0(X,Xtype,LT) \
331 PARTHULL_PROTO0(X,Xtype,LE) \
332 PARTHULL_PROTO0(X,Xtype,EQ) \
333 PARTHULL_PROTO0(X,Xtype,GE) \
334 PARTHULL_PROTO0(X,Xtype,GT) \
335 PARTHULL_PROTO0(X,Xtype,NE)
336
337 /* Use the above macros to define all PartHull prototypes for all
338 data types and operations. */
339 #if HAVE_LONG_DOUBLE /* Not normally implemented */
340 PARTHULL_PROTO(LD,long double)
341 #endif
342 PARTHULL_PROTO(D,double)
343 PARTHULL_PROTO(L,long int)
344 PARTHULL_PROTO(UL,unsigned long int)
345 PARTHULL_PROTO(I,int)
346 PARTHULL_PROTO(UI,unsigned int)
347 PARTHULL_PROTO(S,short int)
348 PARTHULL_PROTO(US,unsigned short int)
349 PARTHULL_PROTO(B,signed char)
350 PARTHULL_PROTO(UB,unsigned char)
351 PARTHULL_PROTO(F,float)
352
353 /* Define a macro that expands to a single prototype for function
354 ConvexHull for a given data type and operation. */
355 #define CONVEXHULL_PROTO0(X,Xtype,Oper) \
356 static AstPointSet *ConvexHull##Oper##X( Xtype, const Xtype[], const int[2], int, int, int, int * );
357
358 /* Define a macro that expands to a set of prototypes for all operations
359 for function ConvexHull for a given data type. */
360 #define CONVEXHULL_PROTO(X,Xtype) \
361 CONVEXHULL_PROTO0(X,Xtype,LT) \
362 CONVEXHULL_PROTO0(X,Xtype,LE) \
363 CONVEXHULL_PROTO0(X,Xtype,EQ) \
364 CONVEXHULL_PROTO0(X,Xtype,GE) \
365 CONVEXHULL_PROTO0(X,Xtype,GT) \
366 CONVEXHULL_PROTO0(X,Xtype,NE)
367
368 /* Use the above macros to define all ConvexHull prototypes for all
369 data types and operations. */
370 #if HAVE_LONG_DOUBLE /* Not normally implemented */
371 CONVEXHULL_PROTO(LD,long double)
372 #endif
373 CONVEXHULL_PROTO(D,double)
374 CONVEXHULL_PROTO(L,long int)
375 CONVEXHULL_PROTO(UL,unsigned long int)
376 CONVEXHULL_PROTO(I,int)
377 CONVEXHULL_PROTO(UI,unsigned int)
378 CONVEXHULL_PROTO(S,short int)
379 CONVEXHULL_PROTO(US,unsigned short int)
380 CONVEXHULL_PROTO(B,signed char)
381 CONVEXHULL_PROTO(UB,unsigned char)
382 CONVEXHULL_PROTO(F,float)
383
384 /* Define a macro that expands to a single prototype for function
385 FindBoxEdge for a given data type and operation. */
386 #define FINDBOXEDGE_PROTO0(X,Xtype,Oper) \
387 static void FindBoxEdge##Oper##X( Xtype, const Xtype[], int, int, int, int, int *, int *, int *, int * );
388
389 /* Define a macro that expands to a set of prototypes for all operations
390 for function FindBoxEdge for a given data type. */
391 #define FINDBOXEDGE_PROTO(X,Xtype) \
392 FINDBOXEDGE_PROTO0(X,Xtype,LT) \
393 FINDBOXEDGE_PROTO0(X,Xtype,LE) \
394 FINDBOXEDGE_PROTO0(X,Xtype,EQ) \
395 FINDBOXEDGE_PROTO0(X,Xtype,GE) \
396 FINDBOXEDGE_PROTO0(X,Xtype,GT) \
397 FINDBOXEDGE_PROTO0(X,Xtype,NE)
398
399 /* Use the above macros to define all FindBoxEdge prototypes for all
400 data types and operations. */
401 #if HAVE_LONG_DOUBLE /* Not normally implemented */
402 FINDBOXEDGE_PROTO(LD,long double)
403 #endif
404 FINDBOXEDGE_PROTO(D,double)
405 FINDBOXEDGE_PROTO(L,long int)
406 FINDBOXEDGE_PROTO(UL,unsigned long int)
407 FINDBOXEDGE_PROTO(I,int)
408 FINDBOXEDGE_PROTO(UI,unsigned int)
409 FINDBOXEDGE_PROTO(S,short int)
410 FINDBOXEDGE_PROTO(US,unsigned short int)
411 FINDBOXEDGE_PROTO(B,signed char)
412 FINDBOXEDGE_PROTO(UB,unsigned char)
413 FINDBOXEDGE_PROTO(F,float)
414
415
416
417
418
419
420
421
422
423 /* Non-generic function prototypes. */
424 static AstMapping *Simplify( AstMapping *, int * );
425 static AstPointSet *DownsizePoly( AstPointSet *, double, int, AstFrame *, int * );
426 static AstPointSet *RegBaseMesh( AstRegion *, int * );
427 static AstPointSet *Transform( AstMapping *, AstPointSet *, int, AstPointSet *, int * );
428 static AstPolygon *Downsize( AstPolygon *, double, int, int * );
429 static Segment *AddToChain( Segment *, Segment *, int * );
430 static Segment *NewSegment( Segment *, int, int, int, int * );
431 static Segment *RemoveFromChain( Segment *, Segment *, int * );
432 static double Polywidth( AstFrame *, AstLineDef **, int, int, double[ 2 ], int * );
433 static int GetBounded( AstRegion *, int * );
434 static int IntCmp( const void *, const void * );
435 static int RegPins( AstRegion *, AstPointSet *, AstRegion *, int **, int * );
436 static int RegTrace( AstRegion *, int, double *, double **, int * );
437 static void Cache( AstPolygon *, int * );
438 static void Copy( const AstObject *, AstObject *, int * );
439 static void Delete( AstObject *, int * );
440 static void Dump( AstObject *, AstChannel *, int * );
441 static void EnsureInside( AstPolygon *, int * );
442 static void FindMax( Segment *, AstFrame *, double *, double *, int, int, int * );
443 static void RegBaseBox( AstRegion *this, double *, double *, int * );
444 static void ResetCache( AstRegion *this, int * );
445 static void SetPointSet( AstPolygon *, AstPointSet *, int * );
446 static void SetRegFS( AstRegion *, AstFrame *, int * );
447 static void SmoothPoly( AstPointSet *, int, double, int * );
448
449 static const char *GetAttrib( AstObject *, const char *, int * );
450 static void ClearAttrib( AstObject *, const char *, int * );
451 static void SetAttrib( AstObject *, const char *, int * );
452 static int TestAttrib( AstObject *, const char *, int * );
453
454 static int GetSimpVertices( AstPolygon *, int * );
455 static int TestSimpVertices( AstPolygon *, int * );
456 static void ClearSimpVertices( AstPolygon *, int * );
457 static void SetSimpVertices( AstPolygon *, int, int * );
458
459
460 /* Member functions. */
461 /* ================= */
AddToChain(Segment * head,Segment * seg,int * status)462 static Segment *AddToChain( Segment *head, Segment *seg, int *status ){
463 /*
464 * Name:
465 * AddToChain
466
467 * Purpose:
468 * Add a Segment into the linked list of Segments, maintaining the
469 * required order in the list.
470
471 * Type:
472 * Private function.
473
474 * Synopsis:
475 * #include "polygon.h"
476 * Segment *AddToChain( Segment *head, Segment *seg, int *status )
477
478 * Class Membership:
479 * Polygon member function
480
481 * Description:
482 * The linked list of Segments maintained by astDownsize is searched
483 * from the high error end (the head), until a Segment is foound which
484 * has a lower error than the supplied segment. The supplied Segment
485 * is then inserted into the list at that point.
486
487 * Parameters:
488 * head
489 * The Segment structure at the head of the list (i.e. the segment
490 * with maximum error).
491 * seg
492 * The Segment to be added into the list.
493 * status
494 * Pointer to the inherited status variable.
495
496 * Returned Value:
497 * Pointer to the link head (which will have changed if "seg" has a
498 * higher error than the original head).
499
500 */
501
502 /* Local Variables: */
503 Segment *tseg;
504
505 /* Check the global error status. */
506 if ( !astOK ) return head;
507
508 /* If the list is empty, return the supplied segment as the new list
509 head. */
510 if( !head ) {
511 head = seg;
512
513 /* If the supplied segment has a higher error than the original head,
514 insert the new segment in front of the original head. */
515 } else if( seg->error > head->error ){
516 seg->next = head;
517 head->prev = seg;
518 head = seg;
519
520 /* Otherwise, move down the list from the head until a segment is found
521 which has a lower error than the supplied Segment. Then insert the
522 supplied segment into the list in front of it. */
523 } else {
524 tseg = head;
525 seg->next = NULL;
526
527 while( tseg->next ) {
528 if( seg->error > tseg->next->error ) {
529 seg->next = tseg->next;
530 seg->prev = tseg;
531 tseg->next->prev = seg;
532 tseg->next = seg;
533 break;
534 }
535 tseg = tseg->next;
536 }
537
538 if( !seg->next ) {
539 tseg->next = seg;
540 seg->prev = tseg;
541 }
542 }
543
544 /* Return the new head. */
545 return head;
546 }
547
Cache(AstPolygon * this,int * status)548 static void Cache( AstPolygon *this, int *status ){
549 /*
550 * Name:
551 * Cache
552
553 * Purpose:
554 * Calculate intermediate values and cache them in the Polygon structure.
555
556 * Type:
557 * Private function.
558
559 * Synopsis:
560 * #include "polygon.h"
561 * void Cache( AstPolygon *this, int *status )
562
563 * Class Membership:
564 * Polygon member function
565
566 * Description:
567 * This function uses the PointSet stored in the parent Region to calculate
568 * some intermediate values which are useful in other methods. These
569 * values are stored within the Polygon structure.
570
571 * Parameters:
572 * this
573 * Pointer to the Polygon.
574 * status
575 * Pointer to the inherited status variable.
576
577 */
578
579 /* Local Variables: */
580 AstFrame *frm; /* Pointer to base Frame in Polygon */
581 double **ptr; /* Pointer to data in the encapsulated PointSet */
582 double end[ 2 ]; /* Start position for edge */
583 double maxwid; /* Maximum polygon width found so far */
584 double polcen[ 2 ]; /* Polygon centre perpendicular to current edge */
585 double polwid; /* Polygon width perpendicular to current edge */
586 double start[ 2 ]; /* Start position for edge */
587 int i; /* Axis index */
588 int nv; /* Number of vertices in Polygon */
589
590 /* Check the global error status. */
591 if ( !astOK ) return;
592
593 /* Do nothing if the cached information is up to date. */
594 if( this->stale ) {
595
596 /* Get a pointer to the base Frame. */
597 frm = astGetFrame( ((AstRegion *) this)->frameset, AST__BASE );
598
599 /* Get the number of vertices. */
600 nv = astGetNpoint( ((AstRegion *) this)->points );
601
602 /* Get pointers to the coordinate data in the parent Region structure. */
603 ptr = astGetPoints( ((AstRegion *) this)->points );
604
605 /* Free any existing edge information in the Polygon structure. */
606 if( this->edges ) {
607 for( i = 0; i < nv; i++ ) {
608 this->edges[ i ] = astFree( this->edges[ i ] );
609 }
610
611 /* Allocate memory to store new edge information if necessary. */
612 } else {
613 this->edges = astMalloc( sizeof( AstLineDef *)*(size_t) nv );
614 this->startsat = astMalloc( sizeof( double )*(size_t) nv );
615 }
616
617 /* Check pointers can be used safely. */
618 if( this->edges ) {
619
620 /* Create and store a description of each edge. Also form the total
621 distance round the polygon, and the distance from the first vertex
622 at which each edge starts. */
623 this->totlen = 0.0;
624 start[ 0 ] = ptr[ 0 ][ nv - 1 ];
625 start[ 1 ] = ptr[ 1 ][ nv - 1 ];
626
627 for( i = 0; i < nv; i++ ) {
628 end[ 0 ] = ptr[ 0 ][ i ];
629 end[ 1 ] = ptr[ 1 ][ i ];
630 this->edges[ i ] = astLineDef( frm, start, end );
631 start[ 0 ] = end[ 0 ];
632 start[ 1 ] = end[ 1 ];
633
634 this->startsat[ i ] = this->totlen;
635 this->totlen += this->edges[ i ]->length;
636 }
637
638 /* We now look for a point that is inside the polygon. We want a point
639 that is well within the polygon, since points that are only just inside
640 the polygon can give numerical problems. Loop round each edge with
641 non-zero length. */
642 maxwid = -1.0;
643 for( i = 0; i < nv; i++ ) {
644 if( this->edges[ i ]->length > 0.0 ) {
645
646 /* We define another line perpendicular to the current edge, passing
647 through the mid point of the edge, extending towards the inside of the
648 polygon. The following function returns the distance we can travel
649 along this line before we hit any of the other polygon edges. It also
650 puts the position corresponding to half that distance into "polcen". */
651 polwid = Polywidth( frm, this->edges, i, nv, polcen, status );
652
653 /* If the width of the polygon perpendicular to the current edge is
654 greater than the width perpdeicular to any other edge, record the
655 width and also store the current polygon centre. */
656 if( polwid > maxwid && polwid != AST__BAD ) {
657 maxwid = polwid;
658 (this->in)[ 0 ] = polcen[ 0 ];
659 (this->in)[ 1 ] = polcen[ 1 ];
660 }
661 }
662 }
663
664 /* If no width was found it probably means that the polygon vertices were
665 given in clockwise order, resulting in the above process probing the
666 infinite extent outside the polygonal hole. In this case any point
667 outside the hole will do, so we use the current contents of the
668 "polcen" array. Set a flag indicating if the vertices are stored in
669 anti-clockwise order. */
670 if( maxwid < 0.0 ) {
671 (this->in)[ 0 ] = polcen[ 0 ];
672 (this->in)[ 1 ] = polcen[ 1 ];
673 this->acw = 0;
674 } else {
675 this->acw = 1;
676 }
677 }
678
679 /* Free resources */
680 frm = astAnnul( frm );
681
682 /* Indicate cached information is up to date. */
683 this->stale = 0;
684 }
685 }
686
ClearAttrib(AstObject * this_object,const char * attrib,int * status)687 static void ClearAttrib( AstObject *this_object, const char *attrib, int *status ) {
688 /*
689 * Name:
690 * ClearAttrib
691
692 * Purpose:
693 * Clear an attribute value for a Polygon.
694
695 * Type:
696 * Private function.
697
698 * Synopsis:
699 * #include "polygon.h"
700 * void ClearAttrib( AstObject *this, const char *attrib, int *status )
701
702 * Class Membership:
703 * Polygon member function (over-rides the astClearAttrib protected
704 * method inherited from the Region class).
705
706 * Description:
707 * This function clears the value of a specified attribute for a
708 * Polygon, so that the default value will subsequently be used.
709
710 * Parameters:
711 * this
712 * Pointer to the Polygon.
713 * attrib
714 * Pointer to a null terminated string specifying the attribute
715 * name. This should be in lower case with no surrounding white
716 * space.
717 * status
718 * Pointer to the inherited status variable.
719 */
720
721 /* Local Variables: */
722 AstPolygon *this; /* Pointer to the Polygon structure */
723
724 /* Check the global error status. */
725 if ( !astOK ) return;
726
727 /* Obtain a pointer to the Polygon structure. */
728 this = (AstPolygon *) this_object;
729
730 /* Check the attribute name and clear the appropriate attribute. */
731
732 /* SimpVertices. */
733 /* ------------- */
734 if ( !strcmp( attrib, "simpvertices" ) ) {
735 astClearSimpVertices( this );
736
737 /* If the attribute is not recognised, pass it on to the parent method
738 for further interpretation. */
739 } else {
740 (*parent_clearattrib)( this_object, attrib, status );
741 }
742 }
743
744 /*
745 *++
746 * Name:
747 c astConvex<X>
748 f AST_CONVEX<X>
749
750 * Purpose:
751 * Create a new Polygon representing the convex hull of a 2D data grid.
752
753 * Type:
754 * Public function.
755
756 * Synopsis:
757 c #include "polygon.h"
758 c AstPolygon *astConvex<X>( <Xtype> value, int oper, const <Xtype> array[],
759 c const int lbnd[2], const int ubnd[2], int starpix )
760 f RESULT = AST_CONVEX<X>( VALUE, OPER, ARRAY, LBND, UBND, STARPIX, STATUS )
761
762 * Class Membership:
763 * Polygon method.
764
765 * Description:
766 * This is a set of functions that create the shortest Polygon that
767 * encloses all pixels with a specified value within a gridded
768 * 2-dimensional data array (e.g. an image).
769 *
770 * A basic 2-dimensional Frame is used to represent the pixel coordinate
771 * system in the returned Polygon. The Domain attribute is set to
772 * "PIXEL", the Title attribute is set to "Pixel coordinates", and the
773 * Unit attribute for each axis is set to "pixel". All other
774 * attributes are left unset. The nature of the pixel coordinate system
775 * is determined by parameter
776 c "starpix".
777 f STARPIX.
778 *
779 * You should use a function which matches the numerical type of the
780 * data you are processing by replacing <X> in the generic function
781 * name
782 c astConvex<X>
783 f AST_CONVEX<X>
784 c by an appropriate 1- or 2-character type code. For example, if you
785 * are procesing data with type
786 c "float", you should use the function astConvexF
787 f REAL, you should use the function AST_CONVEXR
788 * (see the "Data Type Codes" section below for the codes appropriate to
789 * other numerical types).
790
791 * Parameters:
792 c value
793 f VALUE = <Xtype> (Given)
794 * A data value that specifies the pixels to be included within the
795 * convex hull.
796 c oper
797 f OPER = INTEGER (Given)
798 * Indicates how the
799 c "value"
800 f VALUE
801 * parameter is used to select the required pixels. It can
802 * have any of the following values:
803 c - AST__LT: include pixels with value less than "value".
804 c - AST__LE: include pixels with value less than or equal to "value".
805 c - AST__EQ: include pixels with value equal to "value".
806 c - AST__NE: include pixels with value not equal to "value".
807 c - AST__GE: include pixels with value greater than or equal to "value".
808 c - AST__GT: include pixels with value greater than "value".
809 f - AST__LT: include pixels with value less than VALUE.
810 f - AST__LE: include pixels with value less than or equal to VALUE.
811 f - AST__EQ: include pixels with value equal to VALUE.
812 f - AST__NE: include pixels with value not equal to VALUE.
813 f - AST__GE: include pixels with value greater than or equal to VALUE.
814 f - AST__GT: include pixels with value greater than VALUE.
815 c array
816 f ARRAY( * ) = <Xtype> (Given)
817 c Pointer to a
818 f A
819 * 2-dimensional array containing the data to be processed. The
820 * numerical type of this array should match the 1- or
821 * 2-character type code appended to the function name (e.g. if
822 c you are using astConvexF, the type of each array element
823 c should be "float").
824 f you are using AST_CONVEXR, the type of each array element
825 f should be REAL).
826 *
827 * The storage order of data within this array should be such
828 * that the index of the first grid dimension varies most
829 * rapidly and that of the second dimension least rapidly
830 c (i.e. Fortran array indexing is used).
831 f (i.e. normal Fortran array storage order).
832 c lbnd
833 f LBND( 2 ) = INTEGER (Given)
834 c Pointer to an array of two integers
835 f An array
836 * containing the coordinates of the centre of the first pixel
837 * in the input grid along each dimension.
838 c ubnd
839 f UBND( 2) = INTEGER (Given)
840 c Pointer to an array of two integers
841 f An array
842 * containing the coordinates of the centre of the last pixel in
843 * the input grid along each dimension.
844 *
845 c Note that "lbnd" and "ubnd" together define the shape
846 f Note that LBND and UBND together define the shape
847 * and size of the input grid, its extent along a particular
848 c (j'th) dimension being ubnd[j]-lbnd[j]+1 (assuming the
849 c index "j" to be zero-based). They also define
850 f (J'th) dimension being UBND(J)-LBND(J)+1. They also define
851 * the input grid's coordinate system, each pixel having unit
852 * extent along each dimension with integral coordinate values
853 * at its centre or upper corner, as selected by parameter
854 c "starpix".
855 f STARPIX.
856 c starpix
857 f STARPIX = LOGICAL (Given)
858 * A flag indicating the nature of the pixel coordinate system used
859 * to describe the vertex positions in the returned Polygon. If
860 c non-zero,
861 f .TRUE.,
862 * the standard Starlink definition of pixel coordinate is used in
863 * which a pixel with integer index I spans a range of pixel coordinate
864 * from (I-1) to I (i.e. pixel corners have integral pixel coordinates).
865 c If zero,
866 f If .FALSE.,
867 * the definition of pixel coordinate used by other AST functions
868 c such as astResample, astMask,
869 f such as AST_RESAMPLE, AST_MASK,
870 * etc., is used. In this definition, a pixel with integer index I
871 * spans a range of pixel coordinate from (I-0.5) to (I+0.5) (i.e.
872 * pixel centres have integral pixel coordinates).
873 f STATUS = INTEGER (Given and Returned)
874 f The global status.
875
876 * Returned Value:
877 c astConvex<X>()
878 f AST_CONVEX<X> = INTEGER
879 * A pointer to the required Polygon.
880 c NULL
881 f AST__NULL
882 * is returned without error if the array contains no pixels that
883 * satisfy the criterion specified by
884 c "value" and "oper".
885 f VALUE and OPER.
886
887 * Notes:
888 c - NULL
889 f - AST__NULL
890 * will be returned if this function is invoked with the global
891 * error status set, or if it should fail for any reason.
892
893 * Data Type Codes:
894 * To select the appropriate masking function, you should
895 c replace <X> in the generic function name astConvex<X> with a
896 f replace <X> in the generic function name AST_CONVEX<X> with a
897 * 1- or 2-character data type code, so as to match the numerical
898 * type <Xtype> of the data you are processing, as follows:
899 c - D: double
900 c - F: float
901 c - L: long int
902 c - UL: unsigned long int
903 c - I: int
904 c - UI: unsigned int
905 c - S: short int
906 c - US: unsigned short int
907 c - B: byte (signed char)
908 c - UB: unsigned byte (unsigned char)
909 f - D: DOUBLE PRECISION
910 f - R: REAL
911 f - I: INTEGER
912 f - UI: INTEGER (treated as unsigned)
913 f - S: INTEGER*2 (short integer)
914 f - US: INTEGER*2 (short integer, treated as unsigned)
915 f - B: BYTE (treated as signed)
916 f - UB: BYTE (treated as unsigned)
917 *
918 c For example, astConvexD would be used to process "double"
919 c data, while astConvexS would be used to process "short int"
920 c data, etc.
921 f For example, AST_CONVEXD would be used to process DOUBLE
922 f PRECISION data, while AST_CONVEXS would be used to process
923 f short integer data (stored in an INTEGER*2 array), etc.
924 f
925 f For compatibility with other Starlink facilities, the codes W
926 f and UW are provided as synonyms for S and US respectively (but
927 f only in the Fortran interface to AST).
928
929 *--
930 */
931 /* Define a macro to implement the function for a specific data
932 type. Note, this function cannot be a virtual function since the
933 argument list does not include a Polygon, and so no virtual function
934 table is available. */
935 #define MAKE_CONVEX(X,Xtype) \
936 AstPolygon *astConvex##X##_( Xtype value, int oper, const Xtype array[], \
937 const int lbnd[2], const int ubnd[2], \
938 int starpix, int *status ) { \
939 \
940 /* Local Variables: */ \
941 AstFrame *frm; /* Frame in which to define the Polygon */ \
942 AstPointSet *candidate; /* Candidate polygon vertices */ \
943 AstPolygon *result; /* Result value to return */ \
944 int xdim; /* Number of pixels per row */ \
945 int ydim; /* Number of rows */ \
946 static double junk[ 6 ] = {0.0, 0.0, 1.0, 1.0, 0.0, 1.0 }; /* Junk poly */ \
947 \
948 /* Initialise. */ \
949 result = NULL; \
950 \
951 /* Check the global error status. */ \
952 if ( !astOK ) return result; \
953 \
954 /* Get the array dimensions. */ \
955 xdim = ubnd[ 0 ] - lbnd[ 0 ] + 1; \
956 ydim = ubnd[ 1 ] - lbnd[ 1 ] + 1; \
957 \
958 /* Get the basic concvex hull. */ \
959 if( oper == AST__LT ) { \
960 candidate = ConvexHullLT##X( value, array, lbnd, starpix, xdim, ydim, status ); \
961 \
962 } else if( oper == AST__LE ) { \
963 candidate = ConvexHullLE##X( value, array, lbnd, starpix, xdim, ydim, status ); \
964 \
965 } else if( oper == AST__EQ ) { \
966 candidate = ConvexHullEQ##X( value, array, lbnd, starpix, xdim, ydim, status ); \
967 \
968 } else if( oper == AST__NE ) { \
969 candidate = ConvexHullNE##X( value, array, lbnd, starpix, xdim, ydim, status ); \
970 \
971 } else if( oper == AST__GE ) { \
972 candidate = ConvexHullGE##X( value, array, lbnd, starpix, xdim, ydim, status ); \
973 \
974 } else if( oper == AST__GT ) { \
975 candidate = ConvexHullGT##X( value, array, lbnd, starpix, xdim, ydim, status ); \
976 \
977 } else if( astOK ){ \
978 astError( AST__OPRIN, "astConvex"#X": Invalid operation code " \
979 "(%d) supplied (programming error).", status, oper ); \
980 } \
981 \
982 /* Check some good selected values were found. */ \
983 if( candidate ) { \
984 \
985 /* Create a default Polygon with 3 junk vertices. */ \
986 frm = astFrame( 2, "Domain=PIXEL,Unit(1)=pixel,Unit(2)=pixel," \
987 "Title=Pixel coordinates", status ); \
988 result = astPolygon( frm, 3, 3, junk, NULL, " ", status ); \
989 \
990 /* Change the PointSet within the Polygon to the one created above. */ \
991 SetPointSet( result, candidate, status ); \
992 \
993 /* Free resources. */ \
994 frm = astAnnul( frm ); \
995 candidate = astAnnul( candidate ); \
996 } \
997 \
998 /* If an error occurred, clear the returned result. */ \
999 if ( !astOK ) result = astAnnul( result ); \
1000 \
1001 /* Return the result. */ \
1002 return result; \
1003 }
1004
1005
1006 /* Expand the above macro to generate a function for each required
1007 data type. */
1008 #if HAVE_LONG_DOUBLE /* Not normally implemented */
MAKE_CONVEX(LD,long double)1009 MAKE_CONVEX(LD,long double)
1010 #endif
1011 MAKE_CONVEX(D,double)
1012 MAKE_CONVEX(L,long int)
1013 MAKE_CONVEX(UL,unsigned long int)
1014 MAKE_CONVEX(I,int)
1015 MAKE_CONVEX(UI,unsigned int)
1016 MAKE_CONVEX(S,short int)
1017 MAKE_CONVEX(US,unsigned short int)
1018 MAKE_CONVEX(B,signed char)
1019 MAKE_CONVEX(UB,unsigned char)
1020 MAKE_CONVEX(F,float)
1021
1022 /* Undefine the macros. */
1023 #undef MAKE_CONVEX
1024
1025 /*
1026 * Name:
1027 * ConvexHull
1028
1029 * Purpose:
1030 * Find the convex hull enclosing selected pixels in a 2D array.
1031
1032 * Type:
1033 * Private function.
1034
1035 * Synopsis:
1036 * #include "polygon.h"
1037 * AstPointSet *ConvexHull<Oper><X>( <Xtype> value, const <Xtype> array[],
1038 * const int lbnd[2], int starpix,
1039 * int xdim, int ydim, int *status )
1040
1041 * Class Membership:
1042 * Polygon member function
1043
1044 * Description:
1045 * This function uses an algorithm similar to "Andrew's Monotone Chain
1046 * Algorithm" to create a list of vertices describing the convex hull
1047 * enclosing the selected pixels in the supplied array. The vertices
1048 * are returned in a PointSet.
1049
1050 * Parameters:
1051 * value
1052 * A data value that specifies the pixels to be selected.
1053 * array
1054 * Pointer to a 2-dimensional array containing the data to be
1055 * processed. The numerical type of this array should match the 1-
1056 * or 2-character type code appended to the function name.
1057 * lbnd
1058 * The lower pixel index bounds of the array.
1059 * starpix
1060 * If non-zero, the usual Starlink definition of pixel coordinate
1061 * is used (integral values at pixel corners). Otherwise, the
1062 * system used by other AST functions such as astResample is used
1063 * (integral values at pixel centres).
1064 * xdim
1065 * The number of pixels along each row of the array.
1066 * ydim
1067 * The number of rows in the array.
1068 * status
1069 * Pointer to the inherited status variable.
1070
1071 * Returned Value:
1072 * A PointSet holding the vertices of the convex hull, or NULL if an
1073 * error occurs.
1074
1075 * Notes:
1076 * - The following code has been designed with a view to it being
1077 * multi-threaded at some point.
1078
1079 */
1080
1081 /* Define a macro to implement the function for a specific data
1082 type and operation. */
1083 #define MAKE_CONVEXHULL(X,Xtype,Oper,OperI) \
1084 static AstPointSet *ConvexHull##Oper##X( Xtype value, const Xtype array[], \
1085 const int lbnd[2], int starpix, \
1086 int xdim, int ydim, int *status ) { \
1087 \
1088 /* Local Variables: */ \
1089 AstPointSet *result; \
1090 double **ptr; \
1091 double *xv1; \
1092 double *xv2; \
1093 double *xv3; \
1094 double *xv4; \
1095 double *xvert; \
1096 double *yv1; \
1097 double *yv2; \
1098 double *yv3; \
1099 double *yv4; \
1100 double *yvert; \
1101 int nv1; \
1102 int nv2; \
1103 int nv3; \
1104 int nv4; \
1105 int nv; \
1106 int xhi; \
1107 int xhiymax; \
1108 int xhiymin; \
1109 int xlo; \
1110 int xloymax; \
1111 int xloymin; \
1112 int yhi; \
1113 int yhixmax; \
1114 int yhixmin; \
1115 int ylo; \
1116 int yloxmax; \
1117 int yloxmin; \
1118 \
1119 /* Initialise */ \
1120 result = NULL; \
1121 \
1122 /* Check the global error status. */ \
1123 if ( !astOK ) return result; \
1124 \
1125 /* Find the lowest Y value at any selected pixel, and find the max and \
1126 min X value of the selected pixels at that lowest Y value. */ \
1127 FindBoxEdge##Oper##X( value, array, xdim, ydim, 1, 1, &ylo, &yloxmax, \
1128 &yloxmin, status ); \
1129 \
1130 /* Skip if there are no selected values in the array. */ \
1131 if( ylo > 0 ) { \
1132 \
1133 /* Find the highest Y value at any selected pixel, and find the max and \
1134 min X value of the selected pixels at that highest Y value. */ \
1135 FindBoxEdge##Oper##X( value, array, xdim, ydim, 1, 0, &yhi, &yhixmax, \
1136 &yhixmin, status ); \
1137 \
1138 /* Find the lowest X value at any selected pixel, and find the max and \
1139 min Y value of the selected pixels at that lowest X value. */ \
1140 FindBoxEdge##Oper##X( value, array, xdim, ydim, 0, 1, &xlo, &xloymax, \
1141 &xloymin, status ); \
1142 \
1143 /* Find the highest X value at any selected pixel, and find the max and \
1144 min Y value of the selected pixels at that highest X value. */ \
1145 FindBoxEdge##Oper##X( value, array, xdim, ydim, 0, 0, &xhi, &xhiymax, \
1146 &xhiymin, status ); \
1147 \
1148 /* Create a list of vertices for the bottom right corner of the bounding \
1149 box of the selected pixels. */ \
1150 PartHull##Oper##X( value, array, xdim, ydim, yloxmax, ylo, xhi, xhiymin, \
1151 starpix, lbnd, &xv1, &yv1, &nv1, status ); \
1152 \
1153 /* Create a list of vertices for the top right corner of the bounding \
1154 box of the selected pixels. */ \
1155 PartHull##Oper##X( value, array, xdim, ydim, xhi, xhiymax, yhixmax, yhi, \
1156 starpix, lbnd, &xv2, &yv2, &nv2, status ); \
1157 \
1158 /* Create a list of vertices for the top left corner of the bounding \
1159 box of the selected pixels. */ \
1160 PartHull##Oper##X( value, array, xdim, ydim, yhixmin, yhi, xlo, xloymax, \
1161 starpix, lbnd, &xv3, &yv3, &nv3, status ); \
1162 \
1163 /* Create a list of vertices for the bottom left corner of the bounding \
1164 box of the selected pixels. */ \
1165 PartHull##Oper##X( value, array, xdim, ydim, xlo, xloymin, yloxmin, ylo, \
1166 starpix, lbnd, &xv4, &yv4, &nv4, status ); \
1167 \
1168 /* Concatenate the four vertex lists and store them in the returned \
1169 PointSet. */ \
1170 nv = nv1 + nv2 + nv3 + nv4; \
1171 result = astPointSet( nv, 2, " ", status ); \
1172 ptr = astGetPoints( result ); \
1173 if( astOK ) { \
1174 xvert = ptr[ 0 ]; \
1175 yvert = ptr[ 1 ]; \
1176 \
1177 memcpy( xvert, xv1, nv1*sizeof( double ) ); \
1178 memcpy( yvert, yv1, nv1*sizeof( double ) ); \
1179 xvert += nv1; \
1180 yvert += nv1; \
1181 \
1182 memcpy( xvert, xv2, nv2*sizeof( double ) ); \
1183 memcpy( yvert, yv2, nv2*sizeof( double ) ); \
1184 xvert += nv2; \
1185 yvert += nv2; \
1186 \
1187 memcpy( xvert, xv3, nv3*sizeof( double ) ); \
1188 memcpy( yvert, yv3, nv3*sizeof( double ) ); \
1189 xvert += nv3; \
1190 yvert += nv3; \
1191 \
1192 memcpy( xvert, xv4, nv4*sizeof( double ) ); \
1193 memcpy( yvert, yv4, nv4*sizeof( double ) ); \
1194 } \
1195 \
1196 /* Free resources. */ \
1197 xv1 = astFree( xv1 ); \
1198 xv2 = astFree( xv2 ); \
1199 xv3 = astFree( xv3 ); \
1200 xv4 = astFree( xv4 ); \
1201 yv1 = astFree( yv1 ); \
1202 yv2 = astFree( yv2 ); \
1203 yv3 = astFree( yv3 ); \
1204 yv4 = astFree( yv4 ); \
1205 } \
1206 \
1207 /* Free the returned PointSet if an error occurred. */ \
1208 if( result && !astOK ) result = astAnnul( result ); \
1209 \
1210 /* Return the result. */ \
1211 return result; \
1212 }
1213
1214 /* Define a macro that uses the above macro to to create implementations
1215 of ConvexHull for all operations. */
1216 #define MAKEALL_CONVEXHULL(X,Xtype) \
1217 MAKE_CONVEXHULL(X,Xtype,LT,AST__LT) \
1218 MAKE_CONVEXHULL(X,Xtype,LE,AST__LE) \
1219 MAKE_CONVEXHULL(X,Xtype,EQ,AST__EQ) \
1220 MAKE_CONVEXHULL(X,Xtype,NE,AST__NE) \
1221 MAKE_CONVEXHULL(X,Xtype,GE,AST__GE) \
1222 MAKE_CONVEXHULL(X,Xtype,GT,AST__GT)
1223
1224 /* Expand the above macro to generate a function for each required
1225 data type and operation. */
1226 #if HAVE_LONG_DOUBLE /* Not normally implemented */
1227 MAKEALL_CONVEXHULL(LD,long double)
1228 #endif
1229 MAKEALL_CONVEXHULL(D,double)
1230 MAKEALL_CONVEXHULL(L,long int)
1231 MAKEALL_CONVEXHULL(UL,unsigned long int)
1232 MAKEALL_CONVEXHULL(I,int)
1233 MAKEALL_CONVEXHULL(UI,unsigned int)
1234 MAKEALL_CONVEXHULL(S,short int)
1235 MAKEALL_CONVEXHULL(US,unsigned short int)
1236 MAKEALL_CONVEXHULL(B,signed char)
1237 MAKEALL_CONVEXHULL(UB,unsigned char)
1238 MAKEALL_CONVEXHULL(F,float)
1239
1240 /* Undefine the macros. */
1241 #undef MAKE_CONVEXHULL
1242 #undef MAKEALL_CONVEXHULL
1243
1244 static AstPolygon *Downsize( AstPolygon *this, double maxerr, int maxvert,
1245 int *status ) {
1246 /*
1247 *++
1248 * Name:
1249 c astDownsize
1250 f AST_DOWNSIZE
1251
1252 * Purpose:
1253 * Reduce the number of vertices in a Polygon.
1254
1255 * Type:
1256 * Public virtual function.
1257
1258 * Synopsis:
1259 c #include "polygon.h"
1260 c AstPolygon *astDownsize( AstPolygon *this, double maxerr, int maxvert )
1261 f RESULT = AST_DOWNSIZE( THIS, MAXERR, MAXVERT, STATUS )
1262
1263 * Class Membership:
1264 * Polygon method.
1265
1266 * Description:
1267 * This function returns a pointer to a new Polygon that contains a
1268 * subset of the vertices in the supplied Polygon. The subset is
1269 * chosen so that the returned Polygon is a good approximation to
1270 * the supplied Polygon, within the limits specified by the supplied
1271 * parameter values. That is, the density of points in the returned
1272 * Polygon is greater at points where the curvature of the boundary of
1273 * the supplied Polygon is greater.
1274
1275 * Parameters:
1276 c this
1277 f THIS = INTEGER (Given)
1278 * Pointer to the Polygon.
1279 c maxerr
1280 f MAXERR = DOUBLE PRECISION (Given)
1281 * The maximum allowed discrepancy between the supplied and
1282 * returned Polygons, expressed as a geodesic distance within the
1283 * Polygon's coordinate frame. If this is zero or less, the
1284 * returned Polygon will have the number of vertices specified by
1285 c maxvert.
1286 f MAXVERT.
1287 c maxvert
1288 f MAXVERT = INTEGER (Given)
1289 * The maximum allowed number of vertices in the returned Polygon.
1290 * If this is less than 3, the number of vertices in the returned
1291 * Polygon will be the minimum needed to achieve the maximum
1292 * discrepancy specified by
1293 c maxerr.
1294 f MAXERR.
1295 f STATUS = INTEGER (Given and Returned)
1296 f The global status.
1297
1298 * Returned Value:
1299 c astDownsize()
1300 f AST_DOWNSIZE = INTEGER
1301 * Pointer to the new Polygon.
1302
1303 * Notes:
1304 * - A null Object pointer (AST__NULL) will be returned if this
1305 c function is invoked with the AST error status set, or if it
1306 f function is invoked with STATUS set to an error value, or if it
1307 * should fail for any reason.
1308 *--
1309 */
1310
1311 /* Local Variables: */
1312 AstFrame *frm; /* Base Frame from the Polygon */
1313 AstPointSet *pset; /* PointSet holding vertices of downsized polygon */
1314 AstPolygon *result; /* Returned pointer to new Polygon */
1315
1316 /* Initialise. */
1317 result = NULL;
1318
1319 /* Check the global error status. */
1320 if ( !astOK ) return result;
1321
1322 /* Get a pointer to the base Frame of the Polygon. */
1323 frm = astGetFrame( ((AstRegion *) this)->frameset, AST__BASE );
1324
1325 /* Create a PointSet holding the vertices of the downsized polygon. */
1326 pset = DownsizePoly( ((AstRegion *) this)->points, maxerr, maxvert,
1327 frm, status );
1328
1329 /* Take a deep copy of the supplied Polygon. */
1330 result = astCopy( this );
1331
1332 /* Change the PointSet within the result Polygon to the one created above. */ \
1333 SetPointSet( result, pset, status ); \
1334
1335 /* Free resources. */
1336 frm = astAnnul( frm );
1337 pset = astAnnul( pset );
1338
1339 /* If an error occurred, annul the returned Polygon. */
1340 if ( !astOK ) result = astAnnul( result );
1341
1342 /* Return the result. */
1343 return result;
1344 }
1345
DownsizePoly(AstPointSet * pset,double maxerr,int maxvert,AstFrame * frm,int * status)1346 static AstPointSet *DownsizePoly( AstPointSet *pset, double maxerr,
1347 int maxvert, AstFrame *frm, int *status ) {
1348 /*
1349 * Name:
1350 * DownsizePoly
1351
1352 * Purpose:
1353 * Reduce the number of vertices in a Polygon.
1354
1355 * Type:
1356 * Private function.
1357
1358 * Synopsis:
1359 * #include "polygon.h"
1360 * AstPointSet *DownsizePoly( AstPointSet *pset, double maxerr, int maxvert,
1361 * AstFrame *frm, int *status )
1362
1363 * Class Membership:
1364 * Polygon member function.
1365
1366 * Description:
1367 * This function returns a pointer to a new PointSet that contains a
1368 * subset of the vertices in the supplied PointSet. The subset is
1369 * chosen so that the returned polygon is a good approximation to
1370 * the supplied polygon, within the limits specified by the supplied
1371 * parameter values. That is, the density of points in the returned
1372 * polygon is greater at points where the curvature of the boundary of
1373 * the supplied polygon is greater.
1374
1375 * Parameters:
1376 * pset
1377 * Pointer to the PointSet holding the polygon vertices.
1378 * maxerr
1379 * The maximum allowed discrepancy between the supplied and
1380 * returned Polygons, expressed as a geodesic distance within the
1381 * Polygon's coordinate frame. If this is zero or less, the
1382 * returned Polygon will have the number of vertices specified by
1383 * maxvert.
1384 * maxvert
1385 * The maximum allowed number of vertices in the returned Polygon.
1386 * If this is less than 3, the number of vertices in the returned
1387 * Polygon will be the minimum needed to achieve the maximum
1388 * discrepancy specified by
1389 * maxerr.
1390 * frm
1391 * Pointer to the Frame in which the polygon is defined.
1392 * status
1393 * Pointer to the inherited status variable.
1394
1395 * Returned Value:
1396 * Pointer to the new PointSet.
1397
1398 * Notes:
1399 * - A null Object pointer (AST__NULL) will be returned if this
1400 * function is invoked with the AST error status set, or if it
1401 * should fail for any reason.
1402 */
1403
1404 /* Local Variables: */
1405 AstPointSet *result; /* Returned pointer to new PointSet */
1406 Segment *head; /* Pointer to new polygon edge with highest error */
1407 Segment *seg1; /* Pointer to new polygon edge */
1408 Segment *seg2; /* Pointer to new polygon edge */
1409 Segment *seg3; /* Pointer to new polygon edge */
1410 double **ptr; /* Pointer to arrays of axis values */
1411 double *x; /* Pointer to array of X values for old Polygon */
1412 double *xnew; /* Pointer to array of X values for new Polygon */
1413 double *y; /* Pointer to array of Y values for old Polygon */
1414 double *ynew; /* Pointer to array of Y values for new Polygon */
1415 double x1; /* Lowest X value at any vertex */
1416 double x2; /* Highest X value at any vertex */
1417 double y1; /* Lowest Y value at any vertex */
1418 double y2; /* Highest Y value at any vertex */
1419 int *newpoly; /* Holds indices of retained input vertices */
1420 int i1; /* Index of first vertex added to output polygon */
1421 int i1x; /* Index of vertex with lowest X */
1422 int i1y; /* Index of vertex with lowest Y */
1423 int i2; /* Index of second vertex added to output polygon */
1424 int i2x; /* Index of vertex with highest X */
1425 int i2y; /* Index of vertex with highest Y */
1426 int i3; /* Index of third vertex added to output polygon */
1427 int iadd; /* Normalised vertex index */
1428 int iat; /* Index at which to store new vertex index */
1429 int newlen; /* Number of vertices currently in new Polygon */
1430 int nv; /* Number of vertices in old Polygon */
1431
1432 /* Initialise. */
1433 result = NULL;
1434
1435 /* Check the global error status. */
1436 if ( !astOK ) return result;
1437
1438 /* Get the number of vertices in the supplied polygon. */
1439 nv = astGetNpoint( pset );
1440
1441 /* If the maximum error allowed is zero, and the maximum number of
1442 vertices is equal to or greater than the number in the supplied
1443 polygon, just return a deep copy of the supplied PointSet. */
1444 if( maxerr <= 0.0 && ( maxvert < 3 || maxvert >= nv ) ) {
1445 result = astCopy( pset );
1446
1447 /* Otherwise... */
1448 } else {
1449
1450 /* Get pointers to the X and Y arrays holding the vertex coordinates in
1451 the supplied polygon. */
1452 ptr = astGetPoints( pset );
1453 x = ptr[ 0 ];
1454 y = ptr[ 1 ];
1455
1456 /* Allocate memory for an array to hold the original indices of the vertices
1457 to be used to create the returned Polygon. This array is expanded as
1458 needed. */
1459 newpoly = astMalloc( 10*sizeof( int ) );
1460
1461 /* Check the pointers can be used safely. */
1462 if( astOK ) {
1463
1464 /* We first need to decide on three widely spaced vertices which form a
1465 reasonable triangular approximation to the whole polygon. First find
1466 the vertices with the highest and lowest Y values, and the highest and
1467 lowest X values. */
1468 x1 = DBL_MAX;
1469 x2 = -DBL_MAX;
1470 y1 = DBL_MAX;
1471 y2 = -DBL_MAX;
1472
1473 i1y = i1x = 0;
1474 i2y = i2x = nv/2;
1475
1476 for( i3 = 0; i3 < nv; i3++ ) {
1477 if( y[ i3 ] < y1 ) {
1478 y1 = y[ i3 ];
1479 i1y = i3;
1480 } else if( y[ i3 ] > y2 ) {
1481 y2 = y[ i3 ];
1482 i2y = i3;
1483 }
1484
1485 if( x[ i3 ] < x1 ) {
1486 x1 = x[ i3 ];
1487 i1x = i3;
1488 } else if( x[ i3 ] > x2 ) {
1489 x2 = x[ i3 ];
1490 i2x = i3;
1491 }
1492 }
1493
1494 /* Use the axis which spans the greater range. */
1495 if( y2 - y1 > x2 - x1 ) {
1496 i1 = i1y;
1497 i2 = i2y;
1498 } else {
1499 i1 = i1x;
1500 i2 = i2x;
1501 }
1502
1503 /* The index with vertex i1 is definitely going to be one of our three
1504 vertices. We are going to use the line from i1 to i2 to choose the two
1505 other vertices to use. Create a structure describing the segment of the
1506 Polygon from the lowest value on the selected axis (X or Y) to the
1507 highest value. As always, the polygons is traversed in an anti-clockwise
1508 direction. */
1509 seg1 = NewSegment( NULL, i1, i2, nv, status );
1510
1511 /* Find the vertex within this segment which is furthest away from the
1512 line on the right hand side (as moving from vertex i1 to vertex i2). */
1513 FindMax( seg1, frm, x, y, nv, 0, status );
1514
1515 /* Likewise, create a structure describing the remained of the Polygon
1516 (i.e. the segment from the highest value on the selected axis to the
1517 lowest value). Then find the vertex within this segment which is
1518 furthest away from the line on the right hand side. */
1519 seg2 = NewSegment( NULL, i2, i1, nv, status );
1520 FindMax( seg2, frm, x, y, nv, 0, status );
1521
1522 /* Select the segment for which the found vertex is furthest from the
1523 line. */
1524 if( seg2->error > seg1->error ) {
1525
1526 /* If the second segment, we will use the vertex that is farthest from
1527 the line as one of our threee vertices. To ensure that movement from
1528 vertex i1 to i2 to i3 is anti-clockwise, we must use this new vertex
1529 as vertex i3, not i2. */
1530 i3 = seg2->imax;
1531
1532 /* Create a description of the polygon segment from vertex i1 to i3, and
1533 find the vertex which is furthest to the right of the line joining the
1534 two vertices. We use this as the middle vertex (i2). */
1535 seg1 = NewSegment( seg1, i1, i3, nv, status );
1536 FindMax( seg1, frm, x, y, nv, 0, status );
1537 i2 = seg1->imax;
1538
1539 /* Do the same if we are choosing the other segment, ordering the
1540 vertices to retain anti-clockwise movement from i1 to i2 to i3. */
1541 } else {
1542 i2 = seg1->imax;
1543 seg1 = NewSegment( seg1, i2, i1, nv, status );
1544 FindMax( seg1, frm, x, y, nv, 0, status );
1545 i3 = seg1->imax;
1546 }
1547
1548 /* Ensure the vertex indices are in the first cycle. */
1549 if( i2 >= nv ) i2 -= nv;
1550 if( i3 >= nv ) i3 -= nv;
1551
1552 /* Create Segment structures to describe each of these three edges. */
1553 seg1 = NewSegment( seg1, i1, i2, nv, status );
1554 seg2 = NewSegment( seg2, i2, i3, nv, status );
1555 seg3 = NewSegment( NULL, i3, i1, nv, status );
1556
1557 /* Record these 3 vertices in an array holding the original indices of
1558 the vertices to be used to create the returned Polygon. */
1559 newpoly[ 0 ] = i1;
1560 newpoly[ 1 ] = i2;
1561 newpoly[ 2 ] = i3;
1562
1563 /* Indicate the new polygon currently has 3 vertices. */
1564 newlen = 3;
1565
1566 /* Search the old vertices between the start and end of segment 3, looking
1567 for the vertex which lies furthest from the line of segment 3. The
1568 residual between this point and the line is stored in the Segment
1569 structure, as is the index of the vertex at which this maximum residual
1570 occurred. */
1571 FindMax( seg3, frm, x, y, nv, 1, status );
1572
1573 /* The "head" variable points to the head of a double linked list of
1574 Segment structures. This list is ordered by residual, so that the
1575 Segment with the maximum residual is at the head, and the Segment
1576 with the minimum residual is at the tail. Initially "seg3" is at the
1577 head. */
1578 head = seg3;
1579
1580 /* Search the old vertices between the start and end of segment 1, looking
1581 for the vertex which lies furthest from the line of segment 1. The
1582 residual between this point and the line is stored in the Segment
1583 structure, as is the index of the vertex at which this maximum residual
1584 occurred. */
1585 FindMax( seg1, frm, x, y, nv, 1, status );
1586
1587 /* Insert segment 1 into the linked list of Segments, at a position that
1588 maintains the ordering of the segments by error. Thus the head of the
1589 list will still have the max error. */
1590 head = AddToChain( head, seg1, status );
1591
1592 /* Do the same for segment 2. */
1593 FindMax( seg2, frm, x, y, nv, 1, status );
1594 head = AddToChain( head, seg2, status );
1595
1596 /* If the maximum allowed number of vertices in the output Polygon is
1597 less than 3, allow any number of vertices up to the number in the
1598 input Polygon (termination will then be determined just by "maxerr"). */
1599 if( maxvert < 3 ) maxvert = nv;
1600
1601 /* Loop round adding an extra vertex to the returned Polygon until the
1602 maximum residual between the new and old polygons is no more than
1603 "maxerr". Abort early if the specified maximum number of vertices is
1604 reached. */
1605 while( head->error > maxerr && newlen < maxvert ) {
1606
1607 /* The segment at the head of the list has the max error (that is, it is
1608 the segment that departs most from the supplied Polygon). To make the
1609 new polygon a better fit to the old polygon, we add the vertex that is
1610 furthest away from this segment to the new polygon. Remember that a
1611 polygon is cyclic so if the vertex has an index that is greater than the
1612 number of vertices in the old polygon, reduce the index by the number
1613 of vertices in the old polygon. */
1614 iadd = head->imax;
1615 if( iadd >= nv ) iadd -= nv;
1616 iat = newlen++;
1617 newpoly = astGrow( newpoly, newlen, sizeof( int ) );
1618 if( !astOK ) break;
1619 newpoly[ iat ] = iadd;
1620
1621 /* We now split the segment that had the highest error into two segments.
1622 The split occurs at the vertex that had the highest error. */
1623 seg1 = NewSegment( NULL, head->imax, head->i2, nv, status );
1624 seg2 = head;
1625 seg2->i2 = head->imax;
1626
1627 /* We do not know where these two new segments should be in the ordered
1628 linked list, so remove them from the list. */
1629 head = RemoveFromChain( head, seg1, status );
1630 head = RemoveFromChain( head, seg2, status );
1631
1632 /* Find the vertex that deviates most from the first of these two new
1633 segments, and then add the segment into the list of vertices, using
1634 the maximum deviation to determine the position of the segment within
1635 the list. */
1636 FindMax( seg1, frm, x, y, nv, 1, status );
1637 head = AddToChain( head, seg1, status );
1638
1639 /* Do the same for the second new segment. */
1640 FindMax( seg2, frm, x, y, nv, 1, status );
1641 head = AddToChain( head, seg2, status );
1642 }
1643
1644 /* Now we have reached the required accuracy, free resources. */
1645 while( head ) {
1646 seg1 = head;
1647 head = head->next;
1648 seg1 = astFree( seg1 );
1649 }
1650
1651 /* If no vertices have been left out, return a deep copy of the supplied
1652 PointSet. */
1653 if( newlen == nv ) {
1654 result = astCopy( pset );
1655
1656 /* Otherwise, sort the indices of the vertices to be retained so that they
1657 are in the same order as they were in the supplied Polygon. */
1658 } else if( astOK ){
1659 qsort( newpoly, newlen, sizeof( int ), IntCmp );
1660
1661 /* Create a new PointSet and get pointers to its axis values. */
1662 result = astPointSet( newlen, 2, " ", status );
1663 ptr = astGetPoints( result );
1664 xnew = ptr[ 0 ];
1665 ynew = ptr[ 1 ];
1666
1667 /* Copy the axis values for the retained vertices from the old to the new
1668 PointSet. */
1669 if( astOK ) {
1670 for( iat = 0; iat < newlen; iat++ ) {
1671 *(xnew++) = x[ newpoly[ iat ] ];
1672 *(ynew++) = y[ newpoly[ iat ] ];
1673 }
1674 }
1675 }
1676 }
1677
1678 /* Free resources. */
1679 newpoly = astFree( newpoly );
1680 }
1681
1682 /* If an error occurred, annul the returned PointSet. */
1683 if ( !astOK ) result = astAnnul( result );
1684
1685 /* Return the result. */
1686 return result;
1687 }
1688
EnsureInside(AstPolygon * this,int * status)1689 static void EnsureInside( AstPolygon *this, int *status ){
1690 /*
1691 * Name:
1692 * EnsureInside
1693
1694 * Purpose:
1695 * Ensure the unnegated Polygon represents the inside of the polygon.
1696
1697 * Type:
1698 * Private function.
1699
1700 * Synopsis:
1701 * #include "polygon.h"
1702 * void EnsureInside( AstPolygon *this, int *status )
1703
1704 * Class Membership:
1705 * Polygon member function
1706
1707 * Description:
1708 * Reversing the order of the vertices of a Polygon is like negating
1709 * the Polygon. But the parent Region class assumes that an unnegated
1710 * region bounded by closed curves (e.g. boxes, circles, ellipses, etc)
1711 * is bounded. So we need to have a way to ensure that a Polygon also
1712 * follows this convention. So this function reverses the order of the
1713 * vertices in the Polygon, if necessary, to ensure that the unnegated
1714 * Polygon is bounded.
1715
1716 * Parameters:
1717 * this
1718 * The Polygon.
1719 * status
1720 * Pointer to the inherited status variable.
1721
1722 */
1723
1724 /* Local Variables: */
1725 AstRegion *this_region;
1726 double **ptr;
1727 double *p;
1728 double *q;
1729 double tmp;
1730 int bounded;
1731 int i;
1732 int j;
1733 int jmid;
1734 int negated;
1735 int np;
1736
1737 /* Check the global error status. */
1738 if ( !astOK ) return;
1739
1740 /* Is the unnegated Polygon unbounded? If so, we need to reverse the
1741 vertices. */
1742 bounded = astGetBounded( this );
1743 negated = astGetNegated( this );
1744 if( ( bounded && negated ) || ( !bounded && !negated ) ) {
1745 this_region = (AstRegion *) this;
1746
1747 /* Get a pointer to the arrays holding the coordinates at the Polygon
1748 vertices. */
1749 ptr = astGetPoints( this_region->points );
1750
1751 /* Get the number of vertices. */
1752 np = astGetNpoint( this_region->points );
1753
1754 /* Store the index of the last vertex to swap. For odd "np" the central
1755 vertex does not need to be swapped. */
1756 jmid = np/2;
1757
1758 /* Loop round the two axes spanned by the Polygon. */
1759 for( i = 0; i < 2; i++ ) {
1760
1761 /* Get pointers to the first pair of axis values to be swapped - i.e. the
1762 first and last axis values. */
1763 p = ptr[ i ];
1764 q = p + np - 1;
1765
1766 /* Loop round all pairs of axis values. */
1767 for( j = 0; j < jmid; j++ ) {
1768
1769 /* Swap the pair. */
1770 tmp = *p;
1771 *(p++) = *q;
1772 *(q--) = tmp;
1773 }
1774 }
1775
1776 /* Invert the value of the "Negated" attribute to cancel out the effect
1777 of the above vertex reversal. */
1778 astNegate( this );
1779
1780 /* Indicate the cached information in the Polygon structure is stale. */
1781 this->stale = 1;
1782 }
1783 }
1784
1785 /*
1786 * Name:
1787 * FindBoxEdge
1788
1789 * Purpose:
1790 * Find an edge of the bounding box containing the selected pixels.
1791
1792 * Type:
1793 * Private function.
1794
1795 * Synopsis:
1796 * #include "polygon.h"
1797 * void FindBoxEdge<Oper><X>( <Xtype> value, const <Xtype> array[],
1798 * int xdim, int ydim, int axis, int low,
1799 * int *val, int *valmax, int *valmin,
1800 * int *status )
1801
1802 * Class Membership:
1803 * Polygon member function
1804
1805 * Description:
1806 * This function search for an edge of the bounding box containing the
1807 * selected pixels in the supplied array.
1808
1809 * Parameters:
1810 * value
1811 * A data value that specifies the pixels to be selected.
1812 * array
1813 * Pointer to a 2-dimensional array containing the data to be
1814 * processed. The numerical type of this array should match the 1-
1815 * or 2-character type code appended to the function name.
1816 * xdim
1817 * The number of pixels along each row of the array.
1818 * ydim
1819 * The number of rows in the array.
1820 * axis
1821 * The index (0 or 1) of the pixel axis perpendicular to the edge of the
1822 * bounding box being found.
1823 * low
1824 * If non-zero, the lower edge of the bounding box on the axis
1825 * specified by "axis" is found and returned. Otherwise, the higher
1826 * edge of the bounding box on the axis specified by "axis" is
1827 * found and returned.
1828 * val
1829 * Address of an int in which to return the value on axis "axis" of
1830 * the higher or lower (as specified by "low") edge of the bounding
1831 * box.
1832 * valmax
1833 * Address of an int in which to return the highest value on axis
1834 * "1-axis" of the selected pixels on the returned edge.
1835 * valmin
1836 * Address of an int in which to return the lowest value on axis
1837 * "1-axis" of the selected pixels on the returned edge.
1838 * status
1839 * Pointer to the inherited status variable.
1840
1841 * Notes;
1842 * - Zero is returned for "*val" if no good values are found, or if
1843 * an error occurs.
1844
1845 */
1846
1847 /* Define a macro to implement the function for a specific data
1848 type and operation. */
1849 #define MAKE_FINDBOXEDGE(X,Xtype,Oper,OperI) \
1850 static void FindBoxEdge##Oper##X( Xtype value, const Xtype array[], int xdim, \
1851 int ydim, int axis, int low, int *val, \
1852 int *valmax, int *valmin, int *status ) { \
1853 \
1854 /* Local Variables: */ \
1855 int astep; \
1856 int bstep; \
1857 int a0; \
1858 int a1; \
1859 int b0; \
1860 int b1; \
1861 int inc; \
1862 int a; \
1863 int b; \
1864 const Xtype *pc; \
1865 \
1866 /* Initialise. */ \
1867 *val = 0; \
1868 *valmin = 0; \
1869 *valmax = 0; \
1870 \
1871 /* Check the global error status. */ \
1872 if ( !astOK ) return; \
1873 \
1874 /* If we are finding an edge that is parallel to the X axis... */ \
1875 if( axis ) { \
1876 \
1877 /* Get the vector step between adjacent pixels on the selected axis, and \
1878 on the other axis. */ \
1879 astep = xdim; \
1880 bstep = 1; \
1881 \
1882 /* Get the first and last value to check on the selected axis, and the \
1883 increment between checks. */ \
1884 if( low ) { \
1885 a0 = 1; \
1886 a1 = ydim; \
1887 inc = 1; \
1888 } else { \
1889 a0 = ydim; \
1890 a1 = 1; \
1891 inc = -1; \
1892 } \
1893 \
1894 /* The first and last value to check on the other axis. */ \
1895 b0 = 1; \
1896 b1 = xdim; \
1897 \
1898 /* Do the same if we are finding an edge that is parallel to the Y axis. */ \
1899 } else { \
1900 astep = 1; \
1901 bstep = xdim; \
1902 \
1903 if( low ) { \
1904 a0 = 1; \
1905 a1 = xdim; \
1906 inc = 1; \
1907 } else { \
1908 a0 = xdim; \
1909 a1 = 1; \
1910 inc = -1; \
1911 } \
1912 \
1913 b0 = 1; \
1914 b1 = ydim; \
1915 } \
1916 \
1917 /* Loop round the axis values. */ \
1918 a = a0; \
1919 while( 1 ) { \
1920 \
1921 /* Get a pointer to the first value to be checked at this axis value. */ \
1922 pc = array + (a - 1)*astep + (b0 - 1)*bstep; \
1923 \
1924 /* Scan the other axis to find the first and last selected pixel. */ \
1925 for( b = b0; b <= b1; b++ ) { \
1926 if( ISVALID(*pc,OperI,value) ) { \
1927 if( *valmin == 0 ) *valmin = b; \
1928 *valmax = b; \
1929 } \
1930 pc += bstep; \
1931 } \
1932 \
1933 /* If any selected pixels were found, store the axis value and exit. */ \
1934 if( *valmax ) { \
1935 *val = a; \
1936 break; \
1937 } \
1938 \
1939 /* Move on to the next axis value. */ \
1940 if( a != a1 ) { \
1941 a += inc; \
1942 } else { \
1943 break; \
1944 } \
1945 \
1946 } \
1947 }
1948
1949 /* Define a macro that uses the above macro to to create implementations
1950 of FindBoxEdge for all operations. */
1951 #define MAKEALL_FINDBOXEDGE(X,Xtype) \
1952 MAKE_FINDBOXEDGE(X,Xtype,LT,AST__LT) \
1953 MAKE_FINDBOXEDGE(X,Xtype,LE,AST__LE) \
1954 MAKE_FINDBOXEDGE(X,Xtype,EQ,AST__EQ) \
1955 MAKE_FINDBOXEDGE(X,Xtype,NE,AST__NE) \
1956 MAKE_FINDBOXEDGE(X,Xtype,GE,AST__GE) \
1957 MAKE_FINDBOXEDGE(X,Xtype,GT,AST__GT)
1958
1959 /* Expand the above macro to generate a function for each required
1960 data type and operation. */
1961 #if HAVE_LONG_DOUBLE /* Not normally implemented */
MAKEALL_FINDBOXEDGE(LD,long double)1962 MAKEALL_FINDBOXEDGE(LD,long double)
1963 #endif
1964 MAKEALL_FINDBOXEDGE(D,double)
1965 MAKEALL_FINDBOXEDGE(L,long int)
1966 MAKEALL_FINDBOXEDGE(UL,unsigned long int)
1967 MAKEALL_FINDBOXEDGE(I,int)
1968 MAKEALL_FINDBOXEDGE(UI,unsigned int)
1969 MAKEALL_FINDBOXEDGE(S,short int)
1970 MAKEALL_FINDBOXEDGE(US,unsigned short int)
1971 MAKEALL_FINDBOXEDGE(B,signed char)
1972 MAKEALL_FINDBOXEDGE(UB,unsigned char)
1973 MAKEALL_FINDBOXEDGE(F,float)
1974
1975 /* Undefine the macros. */
1976 #undef MAKE_FINDBOXEDGE
1977 #undef MAKEALL_FINDBOXEDGE
1978
1979 /*
1980 * Name:
1981 * FindInsidePoint
1982
1983 * Purpose:
1984 * Find a point that is inside the required outline.
1985
1986 * Type:
1987 * Private function.
1988
1989 * Synopsis:
1990 * #include "polygon.h"
1991 * void FindInsidePoint<Oper><X>( <Xtype> value, const <Xtype> array[],
1992 * const int lbnd[ 2 ], const int ubnd[ 2 ],
1993 * int *inx, int *iny, int *iv,
1994 * int *status );
1995
1996 * Class Membership:
1997 * Polygon member function
1998
1999 * Description:
2000 * The central pixel in the array is checked to see if its value meets
2001 * the requirements implied by <Oper> and "value". If so, its pixel
2002 * indices and vector index are returned> if not, a spiral search is
2003 * made until such a pixel value is found.
2004
2005 * Parameters:
2006 * value
2007 * The data value defining valid pixels.
2008 * array
2009 * The data array.
2010 * lbnd
2011 * The lower pixel index bounds of the array.
2012 * ubnd
2013 * The upper pixel index bounds of the array.
2014 * inx
2015 * Pointer to an int in which to return the X pixel index of the
2016 * first point that meets the requirements implied by <oper> and
2017 * "value".
2018 * iny
2019 * Pointer to an int in which to return the Y pixel index of the
2020 * first point that meets the requirements implied by <oper> and
2021 * "value".
2022 * iv
2023 * Pointer to an int in which to return the vector index of the
2024 * first point that meets the requirements implied by <oper> and
2025 * "value".
2026 * status
2027 * Pointer to the inherited status variable.
2028
2029 * Notes:
2030 * - <Oper> must be one of LT, LE, EQ, NE, GE, GT.
2031
2032
2033 */
2034
2035 /* Define a macro to implement the function for a specific data
2036 type and operation. */
2037 #define MAKE_FINDINSIDEPOINT(X,Xtype,Oper,OperI) \
2038 static void FindInsidePoint##Oper##X( Xtype value, const Xtype array[], \
2039 const int lbnd[ 2 ], const int ubnd[ 2 ], \
2040 int *inx, int *iny, int *iv, \
2041 int *status ){ \
2042 \
2043 /* Local Variables: */ \
2044 const Xtype *pv; /* Pointer to next data value to test */ \
2045 const char *text; /* Pointer to text describing oper */ \
2046 int cy; /* Central row index */ \
2047 int iskin; /* Index of spiral layer being searched */ \
2048 int nskin; /* Number of spiral layers to search */ \
2049 int nx; /* Pixel per row */ \
2050 int tmp; /* Temporary storage */ \
2051 int xhi; /* High X pixel index bound of current skin */ \
2052 int xlo; /* Low X pixel index bound of current skin */ \
2053 int yhi; /* High X pixel index bound of current skin */ \
2054 int ylo; /* Low X pixel index bound of current skin */ \
2055 \
2056 /* Check the global error status. */ \
2057 if ( !astOK ) return; \
2058 \
2059 /* Find number of pixels in one row of the array. */ \
2060 nx = ( ubnd[ 0 ] - lbnd[ 0 ] + 1 ); \
2061 \
2062 /* Find the pixel indices of the central pixel */ \
2063 *inx = ( ubnd[ 0 ] + lbnd[ 0 ] )/2; \
2064 *iny = ( ubnd[ 1 ] + lbnd[ 1 ] )/2; \
2065 \
2066 /* Initialise the vector index and pointer to refer to the central pixel. */ \
2067 *iv = ( *inx - lbnd[ 0 ] ) + nx*( *iny - lbnd[ 1 ] ) ; \
2068 pv = array + (*iv); \
2069 \
2070 /* Test the pixel value, returning if it is valid. This \
2071 relies on the compiler optimisation to remove the "if" statements \
2072 for all but the operation appropriate to the function being defined. */ \
2073 if( OperI == AST__LT ) { \
2074 if( *pv < value ) return; \
2075 \
2076 } else if( OperI == AST__LE ) { \
2077 if( *pv <= value ) return; \
2078 \
2079 } else if( OperI == AST__EQ ) { \
2080 if( *pv == value ) return; \
2081 \
2082 } else if( OperI == AST__NE ) { \
2083 if( *pv != value ) return; \
2084 \
2085 } else if( OperI == AST__GE ) { \
2086 if( *pv >= value ) return; \
2087 \
2088 } else { \
2089 if( *pv > value ) return; \
2090 } \
2091 \
2092 /* The central pixel is invalid if we arrive here. So we need to do a \
2093 spiral search out from the centre looking for a valid pixel. Record \
2094 the central row index (this is the row at which we jump to the next \
2095 outer skin when doing the spiral search below). */ \
2096 cy = *iny; \
2097 \
2098 /* Find how many skins can be searched as part of the spiral search \
2099 before the edge of the array is encountered. */ \
2100 nskin = ubnd[ 0 ] - *inx; \
2101 tmp = *inx - lbnd[ 0 ]; \
2102 if( tmp < nskin ) nskin = tmp; \
2103 tmp = ubnd[ 1 ] - *iny; \
2104 if( tmp < nskin ) nskin = tmp; \
2105 tmp = *iny - lbnd[ 1 ]; \
2106 if( tmp < nskin ) nskin = tmp; \
2107 \
2108 /* Initialise the skin box bounds to be just the central pixel. */ \
2109 xlo = xhi = *inx; \
2110 ylo = yhi = *iny; \
2111 \
2112 /* Loop round each skin looking for a valid test pixel. */ \
2113 for( iskin = 0; iskin < nskin; iskin++ ) { \
2114 \
2115 /* Increment the upper and lower bounds of the box forming the next \
2116 skin. */ \
2117 xhi++; \
2118 xlo--; \
2119 yhi++; \
2120 ylo--; \
2121 \
2122 /* Initialise things for the first pixel in the new skin by moving one \
2123 pixel to the right. */ \
2124 pv++; \
2125 (*iv)++; \
2126 (*inx)++; \
2127 \
2128 /* Move up the right hand edge of the box corresponding to the current \
2129 skin. We start at the middle of the right hand edge. */ \
2130 for( *iny = cy; *iny <= yhi; (*iny)++ ) { \
2131 \
2132 /* Test the pixel value, returning if it is valid. This relies on the \
2133 compiler optimisation to remove the "if" statements for all but the \
2134 operation appropriate to the function being defined. */ \
2135 if( OperI == AST__LT ) { \
2136 if( *pv < value ) return; \
2137 \
2138 } else if( OperI == AST__LE ) { \
2139 if( *pv <= value ) return; \
2140 \
2141 } else if( OperI == AST__EQ ) { \
2142 if( *pv == value ) return; \
2143 \
2144 } else if( OperI == AST__NE ) { \
2145 if( *pv != value ) return; \
2146 \
2147 } else if( OperI == AST__GE ) { \
2148 if( *pv >= value ) return; \
2149 \
2150 } else { \
2151 if( *pv > value ) return; \
2152 } \
2153 \
2154 /* Move up a pixel. */ \
2155 pv += nx; \
2156 *iv += nx; \
2157 } \
2158 \
2159 /* Move down a pixel so that *iny == yhi. */ \
2160 pv -= nx; \
2161 *iv -= nx; \
2162 (*iny)--; \
2163 \
2164 /* Move left along the top edge of the box corresponding to the current \
2165 skin. */ \
2166 for( *inx = xhi; *inx >= xlo; (*inx)-- ) { \
2167 \
2168 /* Test the pixel value, returning if it is valid. */ \
2169 if( OperI == AST__LT ) { \
2170 if( *pv < value ) return; \
2171 \
2172 } else if( OperI == AST__LE ) { \
2173 if( *pv <= value ) return; \
2174 \
2175 } else if( OperI == AST__EQ ) { \
2176 if( *pv == value ) return; \
2177 \
2178 } else if( OperI == AST__NE ) { \
2179 if( *pv != value ) return; \
2180 \
2181 } else if( OperI == AST__GE ) { \
2182 if( *pv >= value ) return; \
2183 \
2184 } else { \
2185 if( *pv > value ) return; \
2186 } \
2187 \
2188 /* Move left a pixel. */ \
2189 pv--; \
2190 (*iv)--; \
2191 } \
2192 \
2193 /* Move right a pixel so that *inx == xlo. */ \
2194 pv++; \
2195 (*iv)++; \
2196 (*inx)++; \
2197 \
2198 /* Move down along the left hand edge of the box corresponding to the current \
2199 skin. */ \
2200 for( *iny = yhi; *iny >= ylo; (*iny)-- ) { \
2201 \
2202 /* Test the pixel value, returning if it is valid. */ \
2203 if( OperI == AST__LT ) { \
2204 if( *pv < value ) return; \
2205 \
2206 } else if( OperI == AST__LE ) { \
2207 if( *pv <= value ) return; \
2208 \
2209 } else if( OperI == AST__EQ ) { \
2210 if( *pv == value ) return; \
2211 \
2212 } else if( OperI == AST__NE ) { \
2213 if( *pv != value ) return; \
2214 \
2215 } else if( OperI == AST__GE ) { \
2216 if( *pv >= value ) return; \
2217 \
2218 } else { \
2219 if( *pv > value ) return; \
2220 } \
2221 \
2222 /* Move down a pixel. */ \
2223 pv -= nx; \
2224 *iv -= nx; \
2225 } \
2226 \
2227 /* Move up a pixel so that *iny == ylo. */ \
2228 pv += nx; \
2229 *iv += nx; \
2230 (*iny)++; \
2231 \
2232 /* Move right along the bottom edge of the box corresponding to the current \
2233 skin. */ \
2234 for( *inx = xlo; *inx <= xhi; (*inx)++ ) { \
2235 \
2236 /* Test the pixel value, returning if it is valid. */ \
2237 if( OperI == AST__LT ) { \
2238 if( *pv < value ) return; \
2239 \
2240 } else if( OperI == AST__LE ) { \
2241 if( *pv <= value ) return; \
2242 \
2243 } else if( OperI == AST__EQ ) { \
2244 if( *pv == value ) return; \
2245 \
2246 } else if( OperI == AST__NE ) { \
2247 if( *pv != value ) return; \
2248 \
2249 } else if( OperI == AST__GE ) { \
2250 if( *pv >= value ) return; \
2251 \
2252 } else { \
2253 if( *pv > value ) return; \
2254 } \
2255 \
2256 /* Move right a pixel. */ \
2257 pv++; \
2258 (*iv)++; \
2259 } \
2260 \
2261 /* Move left a pixel so that *inx == xhi. */ \
2262 pv--; \
2263 (*iv)--; \
2264 (*inx)--; \
2265 \
2266 /* Move up the right hand edge of the box correspionding to the current \
2267 skin. We stop just below the middle of the right hand edge. */ \
2268 for( *iny = ylo; *iny < cy; (*iny)++ ) { \
2269 \
2270 /* Test the pixel value, returning if it is valid. This relies on the \
2271 compiler optimisation to remove the "if" statements for all but the \
2272 operation appropriate to the function being defined. */ \
2273 if( OperI == AST__LT ) { \
2274 if( *pv < value ) return; \
2275 \
2276 } else if( OperI == AST__LE ) { \
2277 if( *pv <= value ) return; \
2278 \
2279 } else if( OperI == AST__EQ ) { \
2280 if( *pv == value ) return; \
2281 \
2282 } else if( OperI == AST__NE ) { \
2283 if( *pv != value ) return; \
2284 \
2285 } else if( OperI == AST__GE ) { \
2286 if( *pv >= value ) return; \
2287 \
2288 } else { \
2289 if( *pv > value ) return; \
2290 } \
2291 \
2292 /* Move up a pixel. */ \
2293 pv += nx; \
2294 *iv += nx; \
2295 } \
2296 } \
2297 \
2298 /* Report an error if no inside pooint could be found. */ \
2299 if( OperI == AST__LT ) { \
2300 text = "less than"; \
2301 } else if( OperI == AST__LE ) { \
2302 text = "less than or equal to"; \
2303 } else if( OperI == AST__EQ ) { \
2304 text = "equal to"; \
2305 } else if( OperI == AST__NE ) { \
2306 text = "not equal to"; \
2307 } else if( OperI == AST__GE ) { \
2308 text = "greater than or equal to"; \
2309 } else { \
2310 text = "greater than"; \
2311 } \
2312 astError( AST__NONIN, "astOutline"#X": Could not find a pixel value %s " \
2313 "%g in the supplied array.", status, text, (double) value ); \
2314 }
2315
2316 /* Define a macro that uses the above macro to to create implementations
2317 of FindInsidePoint for all operations. */
2318 #define MAKEALL_FINDINSIDEPOINT(X,Xtype) \
2319 MAKE_FINDINSIDEPOINT(X,Xtype,LT,AST__LT) \
2320 MAKE_FINDINSIDEPOINT(X,Xtype,LE,AST__LE) \
2321 MAKE_FINDINSIDEPOINT(X,Xtype,EQ,AST__EQ) \
2322 MAKE_FINDINSIDEPOINT(X,Xtype,GE,AST__GE) \
2323 MAKE_FINDINSIDEPOINT(X,Xtype,GT,AST__GT) \
2324 MAKE_FINDINSIDEPOINT(X,Xtype,NE,AST__NE)
2325
2326 /* Expand the above macro to generate a function for each required
2327 data type and operation. */
2328 #if HAVE_LONG_DOUBLE /* Not normally implemented */
2329 MAKEALL_FINDINSIDEPOINT(LD,long double)
2330 #endif
2331 MAKEALL_FINDINSIDEPOINT(D,double)
2332 MAKEALL_FINDINSIDEPOINT(L,long int)
2333 MAKEALL_FINDINSIDEPOINT(UL,unsigned long int)
2334 MAKEALL_FINDINSIDEPOINT(I,int)
2335 MAKEALL_FINDINSIDEPOINT(UI,unsigned int)
2336 MAKEALL_FINDINSIDEPOINT(S,short int)
2337 MAKEALL_FINDINSIDEPOINT(US,unsigned short int)
2338 MAKEALL_FINDINSIDEPOINT(B,signed char)
2339 MAKEALL_FINDINSIDEPOINT(UB,unsigned char)
2340 MAKEALL_FINDINSIDEPOINT(F,float)
2341
2342 /* Undefine the macros. */
2343 #undef MAKE_FINDINSIDEPOINT
2344 #undef MAKEALL_FINDINSIDEPOINT
2345
2346 static void FindMax( Segment *seg, AstFrame *frm, double *x, double *y,
2347 int nv, int abs, int *status ){
2348 /*
2349 * Name:
2350 * FindMax
2351
2352 * Purpose:
2353 * Find the maximum discrepancy between a given line segment and the
2354 * Polygon being downsized.
2355
2356 * Type:
2357 * Private function.
2358
2359 * Synopsis:
2360 * #include "polygon.h"
2361 * void FindMax( Segment *seg, AstFrame *frm, double *x, double *y,
2362 * int nv, int abs, int *status )
2363
2364 * Class Membership:
2365 * Polygon member function
2366
2367 * Description:
2368 * The supplied Segment structure describes a range of vertices in
2369 * the polygon being downsized. This function checks each of these
2370 * vertices to find the one that lies furthest from the line joining the
2371 * first and last vertices in the segment. The maximum error, and the
2372 * vertex index at which this maximum error is found, is stored in the
2373 * Segment structure.
2374
2375 * Parameters:
2376 * seg
2377 * The structure describing the range of vertices to check. It
2378 * corresponds to a candidate edge in the downsized polygon.
2379 * frm
2380 * The Frame in which the positions are defined.
2381 * x
2382 * Pointer to the X axis values in the original Polygon.
2383 * y
2384 * Pointer to the Y axis values in the original Polygon.
2385 * nv
2386 * Total number of vertics in the old Polygon..
2387 * abs
2388 * If non-zero, then the stored maximum is the position with
2389 * maximum absolute error. Otherwise, the stored maximum is the
2390 * position with maximum positive error (positive errors are to the
2391 * right when travelling from start to end of the segment).
2392 * status
2393 * Pointer to the inherited status variable.
2394
2395 */
2396
2397 /* Local Variables: */
2398 AstPointSet *pset1; /* PointSet holding vertex positions */
2399 AstPointSet *pset2; /* PointSet holding offset par/perp components */
2400 double **ptr1; /* Pointers to "pset1" data arrays */
2401 double **ptr2; /* Pointers to "pset2" data arrays */
2402 double *px; /* Pointer to next X value */
2403 double *py; /* Pointer to next Y value */
2404 double ax; /* X value at start */
2405 double ay; /* Y value at start */
2406 double ba; /* Distance between a and b */
2407 double bax; /* X increment from a to b */
2408 double bay; /* Y increment from a to b */
2409 double cax; /* X increment from a to c */
2410 double cay; /* Y increment from a to c */
2411 double end[ 2 ]; /* Position of starting vertex */
2412 double error; /* Error value for current vertex */
2413 double start[ 2 ]; /* Position of starting vertex */
2414 int i1; /* Starting index (always in first cycle) */
2415 int i2; /* Ending index converted to first cycle */
2416 int i2b; /* Upper vertex limit in first cycle */
2417 int i; /* Loop count */
2418 int n; /* Number of vertices to check */
2419
2420 /* Check the global error status. */
2421 if ( !astOK ) return;
2422
2423 /* Stuff needed for handling cyclic redundancy of vertex indices. */
2424 i1 = seg->i1;
2425 i2 = seg->i2;
2426 n = i2 - i1 - 1;
2427 i2b = i2;
2428 if( i2 >= nv ) {
2429 i2 -= nv;
2430 i2b = nv;
2431 }
2432
2433 /* If the segment has no intermediate vertices, set the segment error to
2434 zero and return. */
2435 if( n < 1 ) {
2436 seg->error = 0.0;
2437 seg->imax = i1;
2438
2439 /* For speed, we use simple plane geometry if the Polygon is defined in a
2440 simple Frame. */
2441 } else if( !strcmp( astGetClass( frm ), "Frame" ) ) {
2442
2443 /* Point "a" is the vertex that marks the start of the segment. Point "b"
2444 is the vertex that marks the end of the segment. */
2445 ax = x[ i1 ];
2446 ay = y[ i1 ];
2447 bax = x[ i2 ] - ax;
2448 bay = y[ i2 ] - ay;
2449 ba = sqrt( bax*bax + bay*bay );
2450
2451 /* Initialise the largest error found so far. */
2452 seg->error = -1.0;
2453
2454 /* Check the vertices from the start (plus one) up to the end (minus one)
2455 or the last vertex (which ever comes first). */
2456 for( i = i1 + 1; i < i2b; i++ ) {
2457
2458 /* Position "c" is the vertex being tested. Find the squared distance from
2459 "c" to the line joining "a" and "b". */
2460 cax = x[ i ] - ax;
2461 cay = y[ i ] - ay;
2462 error = ( bay*cax - cay*bax )/ba;
2463 if( abs ) error = fabs( error );
2464
2465 /* If this is the largest value found so far, record it. Note the error
2466 here is a squared distance. */
2467 if( error > seg->error ) {
2468 seg->error = error;
2469 seg->imax = i;
2470 }
2471 }
2472
2473 /* If the end vertex is in the next cycle, check the remaining vertex
2474 posI would have thought a telentitions in the same way. */
2475 if( i2b != i2 ) {
2476
2477 for( i = 0; i < i2; i++ ) {
2478
2479 cax = x[ i ] - ax;
2480 cay = y[ i ] - ay;
2481 error = ( bay*cax - cay*bax )/ba;
2482 if( abs ) error = fabs( error );
2483
2484 if( error > seg->error ) {
2485 seg->error = error;
2486 seg->imax = i + i2b;
2487 }
2488
2489 }
2490 }
2491
2492 /* If the polygon is not defined in a simple Frame, we use the overloaded
2493 Frame methods to do the geometry. */
2494 } else {
2495
2496 /* Create a PointSet to hold the positions of the vertices to test. We do
2497 not need to test the start or end vertex. */
2498 pset1 = astPointSet( n, 2, " ", status );
2499 ptr1 = astGetPoints( pset1 );
2500 if( astOK ) {
2501
2502 /* Copy the vertex axis values form the start (plus one) up to the end
2503 (minus one) vertex or the last vertex (which ever comes first). */
2504 px = ptr1[ 0 ];
2505 py = ptr1[ 1 ];
2506
2507 for( i = i1 + 1; i < i2b; i++ ){
2508 *(px++) = x[ i ];
2509 *(py++) = y[ i ];
2510 }
2511
2512 /* If the end vertex is in the next cycle, copy the remaining vertex
2513 positions into the PointSet. */
2514 if( i2b != i2 ) {
2515 for( i = 0; i < i2; i++ ) {
2516 *(px++) = x[ i ];
2517 *(py++) = y[ i ];
2518 }
2519 }
2520
2521 /* Record the start and end vertex positions. */
2522 start[ 0 ] = x[ i1 ];
2523 start[ 1 ] = y[ i1 ];
2524 end[ 0 ] = x[ i2 ];
2525 end[ 1 ] = y[ i2 ];
2526
2527 /* Resolve the vector from the start to each vertex into two components,
2528 parallel and perpendicular to the start->end vector. */
2529 pset2 = astResolvePoints( frm, start, end, pset1, NULL );
2530 ptr2 = astGetPoints( pset2 );
2531 if( astOK ) {
2532
2533 /* Find the vertex with largest perpendicular component. */
2534 seg->error = -1.0;
2535 py = ptr2[ 1 ];
2536 for( i = 1; i <= n; i++ ) {
2537
2538 error = *(py++);
2539 if( abs ) error = fabs( error );
2540
2541 if( error > seg->error ) {
2542 seg->error = error;
2543 seg->imax = i + i1;
2544 }
2545
2546 }
2547 }
2548
2549 /* Free resources. */
2550 pset2 = astAnnul( pset2 );
2551 }
2552 pset1 = astAnnul( pset1 );
2553 }
2554 }
2555
GetAttrib(AstObject * this_object,const char * attrib,int * status)2556 static const char *GetAttrib( AstObject *this_object, const char *attrib, int *status ) {
2557 /*
2558 * Name:
2559 * GetAttrib
2560
2561 * Purpose:
2562 * Get the value of a specified attribute for a Polygon.
2563
2564 * Type:
2565 * Private function.
2566
2567 * Synopsis:
2568 * #include "polygon.h"
2569 * const char *GetAttrib( AstObject *this, const char *attrib, int *status )
2570
2571 * Class Membership:
2572 * Polygon member function (over-rides the protected astGetAttrib
2573 * method inherited from the Region class).
2574
2575 * Description:
2576 * This function returns a pointer to the value of a specified
2577 * attribute for a Polygon, formatted as a character string.
2578
2579 * Parameters:
2580 * this
2581 * Pointer to the Polygon.
2582 * attrib
2583 * Pointer to a null-terminated string containing the name of
2584 * the attribute whose value is required. This name should be in
2585 * lower case, with all white space removed.
2586 * status
2587 * Pointer to the inherited status variable.
2588
2589 * Returned Value:
2590 * - Pointer to a null-terminated string containing the attribute
2591 * value.
2592
2593 * Notes:
2594 * - The returned string pointer may point at memory allocated
2595 * within the Polygon, or at static memory. The contents of the
2596 * string may be over-written or the pointer may become invalid
2597 * following a further invocation of the same function or any
2598 * modification of the Polygon. A copy of the string should
2599 * therefore be made if necessary.
2600 * - A NULL pointer will be returned if this function is invoked
2601 * with the global error status set, or if it should fail for any
2602 * reason.
2603 */
2604
2605 /* Local Variables: */
2606 astDECLARE_GLOBALS /* Pointer to thread-specific global data */
2607 AstPolygon *this; /* Pointer to the Polygon structure */
2608 const char *result; /* Pointer value to return */
2609 int ival; /* Integer attribute value */
2610
2611 /* Initialise. */
2612 result = NULL;
2613
2614 /* Check the global error status. */
2615 if ( !astOK ) return result;
2616
2617 /* Get a pointer to the thread specific global data structure. */
2618 astGET_GLOBALS(this_object);
2619
2620 /* Obtain a pointer to the Polygon structure. */
2621 this = (AstPolygon *) this_object;
2622
2623 /* Compare "attrib" with each recognised attribute name in turn,
2624 obtaining the value of the required attribute. If necessary, write
2625 the value into "getattrib_buff" as a null-terminated string in an appropriate
2626 format. Set "result" to point at the result string. */
2627
2628 /* SimpVertices. */
2629 /* ------------- */
2630 if ( !strcmp( attrib, "simpvertices" ) ) {
2631 ival = astGetSimpVertices( this );
2632 if ( astOK ) {
2633 (void) sprintf( getattrib_buff, "%d", ival );
2634 result = getattrib_buff;
2635 }
2636
2637 /* If the attribute name was not recognised, pass it on to the parent
2638 method for further interpretation. */
2639 } else {
2640 result = (*parent_getattrib)( this_object, attrib, status );
2641 }
2642
2643 /* Return the result. */
2644 return result;
2645
2646 }
2647
GetBounded(AstRegion * this,int * status)2648 static int GetBounded( AstRegion *this, int *status ) {
2649 /*
2650 * Name:
2651 * GetBounded
2652
2653 * Purpose:
2654 * Is the Region bounded?
2655
2656 * Type:
2657 * Private function.
2658
2659 * Synopsis:
2660 * #include "polygon.h"
2661 * int GetBounded( AstRegion *this, int *status )
2662
2663 * Class Membership:
2664 * Polygon method (over-rides the astGetBounded method inherited from
2665 * the Region class).
2666
2667 * Description:
2668 * This function returns a flag indicating if the Region is bounded.
2669 * The implementation provided by the base Region class is suitable
2670 * for Region sub-classes representing the inside of a single closed
2671 * curve (e.g. Circle, Interval, Box, etc). Other sub-classes (such as
2672 * CmpRegion, PointList, etc ) may need to provide their own
2673 * implementations.
2674
2675 * Parameters:
2676 * this
2677 * Pointer to the Region.
2678 * status
2679 * Pointer to the inherited status variable.
2680
2681 * Returned Value:
2682 * Non-zero if the Region is bounded. Zero otherwise.
2683
2684 */
2685
2686 /* Local Variables: */
2687 int neg; /* Has the Polygon been negated? */
2688 int result; /* Returned result */
2689
2690 /* Initialise */
2691 result = 0;
2692
2693 /* Check the global error status. */
2694 if ( !astOK ) return result;
2695
2696 /* Ensure cached information is available. */
2697 Cache( (AstPolygon *) this, status );
2698
2699 /* See if the Polygon has been negated. */
2700 neg = astGetNegated( this );
2701
2702 /* If the polygon vertices are stored in anti-clockwise order, then the
2703 polygon is bounded if it has not been negated. */
2704 if( ( (AstPolygon *) this)->acw ) {
2705 result = (! neg );
2706
2707 /* If the polygon vertices are stored in clockwise order, then the
2708 polygon is bounded if it has been negated. */
2709 } else {
2710 result = neg;
2711 }
2712
2713 /* Return the result. */
2714 return result;
2715 }
2716
astInitPolygonVtab_(AstPolygonVtab * vtab,const char * name,int * status)2717 void astInitPolygonVtab_( AstPolygonVtab *vtab, const char *name, int *status ) {
2718 /*
2719 *+
2720 * Name:
2721 * astInitPolygonVtab
2722
2723 * Purpose:
2724 * Initialise a virtual function table for a Polygon.
2725
2726 * Type:
2727 * Protected function.
2728
2729 * Synopsis:
2730 * #include "polygon.h"
2731 * void astInitPolygonVtab( AstPolygonVtab *vtab, const char *name )
2732
2733 * Class Membership:
2734 * Polygon vtab initialiser.
2735
2736 * Description:
2737 * This function initialises the component of a virtual function
2738 * table which is used by the Polygon class.
2739
2740 * Parameters:
2741 * vtab
2742 * Pointer to the virtual function table. The components used by
2743 * all ancestral classes will be initialised if they have not already
2744 * been initialised.
2745 * name
2746 * Pointer to a constant null-terminated character string which contains
2747 * the name of the class to which the virtual function table belongs (it
2748 * is this pointer value that will subsequently be returned by the Object
2749 * astClass function).
2750 *-
2751 */
2752
2753 /* Local Variables: */
2754 astDECLARE_GLOBALS /* Pointer to thread-specific global data */
2755 AstMappingVtab *mapping; /* Pointer to Mapping component of Vtab */
2756 AstRegionVtab *region; /* Pointer to Region component of Vtab */
2757 AstObjectVtab *object; /* Pointer to Object component of Vtab */
2758
2759 /* Check the local error status. */
2760 if ( !astOK ) return;
2761
2762 /* Get a pointer to the thread specific global data structure. */
2763 astGET_GLOBALS(NULL);
2764
2765 /* Initialize the component of the virtual function table used by the
2766 parent class. */
2767 astInitRegionVtab( (AstRegionVtab *) vtab, name );
2768
2769 /* Store a unique "magic" value in the virtual function table. This
2770 will be used (by astIsAPolygon) to determine if an object belongs
2771 to this class. We can conveniently use the address of the (static)
2772 class_check variable to generate this unique value. */
2773 vtab->id.check = &class_check;
2774 vtab->id.parent = &(((AstRegionVtab *) vtab)->id);
2775
2776 /* Initialise member function pointers. */
2777 /* ------------------------------------ */
2778 /* Store pointers to the member functions (implemented here) that provide
2779 virtual methods for this class. */
2780 vtab->Downsize = Downsize;
2781
2782 /* Save the inherited pointers to methods that will be extended, and
2783 replace them with pointers to the new member functions. */
2784 object = (AstObjectVtab *) vtab;
2785 mapping = (AstMappingVtab *) vtab;
2786 region = (AstRegionVtab *) vtab;
2787
2788 parent_transform = mapping->Transform;
2789 mapping->Transform = Transform;
2790
2791 parent_simplify = mapping->Simplify;
2792 mapping->Simplify = Simplify;
2793
2794 parent_setregfs = region->SetRegFS;
2795 region->SetRegFS = SetRegFS;
2796
2797 parent_resetcache = region->ResetCache;
2798 region->ResetCache = ResetCache;
2799
2800 parent_clearattrib = object->ClearAttrib;
2801 object->ClearAttrib = ClearAttrib;
2802 parent_getattrib = object->GetAttrib;
2803 object->GetAttrib = GetAttrib;
2804 parent_setattrib = object->SetAttrib;
2805 object->SetAttrib = SetAttrib;
2806 parent_testattrib = object->TestAttrib;
2807 object->TestAttrib = TestAttrib;
2808
2809 region->RegPins = RegPins;
2810 region->RegBaseMesh = RegBaseMesh;
2811 region->RegBaseBox = RegBaseBox;
2812 region->RegTrace = RegTrace;
2813 region->GetBounded = GetBounded;
2814
2815 vtab->ClearSimpVertices = ClearSimpVertices;
2816 vtab->GetSimpVertices = GetSimpVertices;
2817 vtab->SetSimpVertices = SetSimpVertices;
2818 vtab->TestSimpVertices = TestSimpVertices;
2819
2820 /* Store replacement pointers for methods which will be over-ridden by
2821 new member functions implemented here. */
2822
2823 /* Declare the copy constructor, destructor and class dump
2824 functions. */
2825 astSetDump( vtab, Dump, "Polygon", "Polygonal region" );
2826 astSetDelete( vtab, Delete );
2827 astSetCopy( vtab, Copy );
2828
2829 /* If we have just initialised the vtab for the current class, indicate
2830 that the vtab is now initialised, and store a pointer to the class
2831 identifier in the base "object" level of the vtab. */
2832 if( vtab == &class_vtab ) {
2833 class_init = 1;
2834 astSetVtabClassIdentifier( vtab, &(vtab->id) );
2835 }
2836 }
2837
IntCmp(const void * a,const void * b)2838 static int IntCmp( const void *a, const void *b ){
2839 /*
2840 * Name:
2841 * IntCmp
2842
2843 * Purpose:
2844 * An integer comparison function for the "qsort" function.
2845
2846 * Type:
2847 * Private function.
2848
2849 * Synopsis:
2850 * #include "polygon.h"
2851 * int IntCmp( const void *a, const void *b )
2852
2853 * Class Membership:
2854 * Polygon member function
2855
2856 * Description:
2857 * See the docs for "qsort".
2858
2859 * Parameters:
2860 * a
2861 * Pointer to the first int
2862 * b
2863 * Pointer to the second int
2864
2865 * Returnd Value:
2866 * Positive, negative or zero, depending on whether a is larger than,
2867 * equal to, or less than b.
2868
2869 */
2870
2871 return *((int*)a) - *((int*)b);
2872 }
2873
NewSegment(Segment * seg,int i1,int i2,int nvert,int * status)2874 static Segment *NewSegment( Segment *seg, int i1, int i2, int nvert,
2875 int *status ){
2876 /*
2877 * Name:
2878 * NewSegment
2879
2880 * Purpose:
2881 * Initialise a structure describing a segment of the new Polygon
2882 * created by astDownsize.
2883
2884 * Type:
2885 * Private function.
2886
2887 * Synopsis:
2888 * #include "polygon.h"
2889 * Segment *NewSegment( Segment *seg, int i1, int i2, int nvert,
2890 * int *status )
2891
2892 * Class Membership:
2893 * Polygon member function
2894
2895 * Description:
2896 * This function initialises the contents of a structure describing
2897 * the specified range of vertices within a Polygon. The cyclic nature
2898 * of vertex indices is taken into account.
2899 *
2900 * If no structure is supplied, memory is allocated to hold a new
2901 * structure.
2902
2903 * Parameters:
2904 * seg
2905 * Pointer to a structure to initialise, or NULL if a new structure
2906 * is to be allocated.
2907 * i1
2908 * The index of a vertex within the old Polygon (supplied to
2909 * astDownsize) that marks the start of the new line segment in
2910 * the downsized polygon.
2911 * i2
2912 * The index of a vertex within the old Polygon (supplied to
2913 * astDownsize) that marks the end of the new line segment in
2914 * the downsized polygon.
2915 * nvert
2916 * Total number of vertics in the old Polygon..
2917 * status
2918 * Pointer to the inherited status variable.
2919
2920 * Returnd Value:
2921 * Pointer to the initialised Segment structure. It should be freed using
2922 * astFree when no longer needed.
2923
2924 */
2925
2926 /* Local Variables: */
2927 Segment *result;
2928
2929 /* Check the global error status. */
2930 if ( !astOK ) return NULL;
2931
2932 /* Get a pointer to the structure to be initialised, allocating memory
2933 for a new structure if none was supplied. */
2934 result = seg ? seg : astMalloc( sizeof( Segment ) );
2935
2936 /* Check the pointer can be used safely. */
2937 if( result ) {
2938
2939 /* If the supplied ending index is less than the starting index, the
2940 ending index must have gone all the way round the polygon and started
2941 round again. Increase the ending index by the number of vertices to
2942 put it in the same cycle as the starting index. */
2943 if( i2 < i1 ) i2 += nvert;
2944
2945 /* If the supplied starting index is within the first cycle (i.e. zero ->
2946 nvert-1 ), use the indices as they are (which may mean that the ending
2947 index is greater than nvert, but this is handled correctly by other
2948 functions). */
2949 if( i1 < nvert ) {
2950 result->i1 = i1;
2951 result->i2 = i2;
2952
2953 /* If the supplied starting index is within the second cycle (i.e. nvert
2954 or greater) the ending index will be even greater, so we can reduce
2955 both by "nvert" to put them both in the first cycle. The goal is that
2956 the starting index should always be in the first cycle, but the ending
2957 index may possibly be in the second cycle. */
2958 } else {
2959 result->i1 = i1 - nvert;
2960 result->i2 = i2 - nvert;
2961 }
2962
2963 /* Nullify the links to other Segments */
2964 result->next = NULL;
2965 result->prev = NULL;
2966 }
2967
2968 /* Return the pointer to the new Segment structure. */
2969 return result;
2970 }
2971
2972 /*
2973 *++
2974 * Name:
2975 c astOutline<X>
2976 f AST_OUTLINE<X>
2977
2978 * Purpose:
2979 * Create a new Polygon outling values in a 2D data grid.
2980
2981 * Type:
2982 * Public function.
2983
2984 * Synopsis:
2985 c #include "polygon.h"
2986 c AstPolygon *astOutline<X>( <Xtype> value, int oper, const <Xtype> array[],
2987 c const int lbnd[2], const int ubnd[2], double maxerr,
2988 c int maxvert, const int inside[2], int starpix )
2989 f RESULT = AST_OUTLINE<X>( VALUE, OPER, ARRAY, LBND, UBND, MAXERR,
2990 f MAXVERT, INSIDE, STARPIX, STATUS )
2991
2992 * Class Membership:
2993 * Polygon method.
2994
2995 * Description:
2996 * This is a set of functions that create a Polygon enclosing a single
2997 * contiguous set of pixels that have a specified value within a gridded
2998 * 2-dimensional data array (e.g. an image).
2999 *
3000 * A basic 2-dimensional Frame is used to represent the pixel coordinate
3001 * system in the returned Polygon. The Domain attribute is set to
3002 * "PIXEL", the Title attribute is set to "Pixel coordinates", and the
3003 * Unit attribute for each axis is set to "pixel". All other
3004 * attributes are left unset. The nature of the pixel coordinate system
3005 * is determined by parameter
3006 c "starpix".
3007 f STARPIX.
3008 *
3009 * The
3010 c "maxerr" and "maxvert"
3011 f MAXERR and MAXVERT
3012 * parameters can be used to control how accurately the returned
3013 * Polygon represents the required region in the data array. The
3014 * number of vertices in the returned Polygon will be the minimum
3015 * needed to achieve the required accuracy.
3016 *
3017 * You should use a function which matches the numerical type of the
3018 * data you are processing by replacing <X> in the generic function
3019 * name
3020 c astOutline<X>
3021 f AST_OUTLINE<X>
3022 c by an appropriate 1- or 2-character type code. For example, if you
3023 * are procesing data with type
3024 c "float", you should use the function astOutlineF
3025 f REAL, you should use the function AST_OUTLINER
3026 * (see the "Data Type Codes" section below for the codes appropriate to
3027 * other numerical types).
3028
3029 * Parameters:
3030 c value
3031 f VALUE = <Xtype> (Given)
3032 * A data value that specifies the pixels to be outlined.
3033 c oper
3034 f OPER = INTEGER (Given)
3035 * Indicates how the
3036 c "value"
3037 f VALUE
3038 * parameter is used to select the outlined pixels. It can
3039 * have any of the following values:
3040 c - AST__LT: outline pixels with value less than "value".
3041 c - AST__LE: outline pixels with value less than or equal to "value".
3042 c - AST__EQ: outline pixels with value equal to "value".
3043 c - AST__NE: outline pixels with value not equal to "value".
3044 c - AST__GE: outline pixels with value greater than or equal to "value".
3045 c - AST__GT: outline pixels with value greater than "value".
3046 f - AST__LT: outline pixels with value less than VALUE.
3047 f - AST__LE: outline pixels with value less than or equal to VALUE.
3048 f - AST__EQ: outline pixels with value equal to VALUE.
3049 f - AST__NE: outline pixels with value not equal to VALUE.
3050 f - AST__GE: outline pixels with value greater than or equal to VALUE.
3051 f - AST__GT: outline pixels with value greater than VALUE.
3052 c array
3053 f ARRAY( * ) = <Xtype> (Given)
3054 c Pointer to a
3055 f A
3056 * 2-dimensional array containing the data to be processed. The
3057 * numerical type of this array should match the 1- or
3058 * 2-character type code appended to the function name (e.g. if
3059 c you are using astOutlineF, the type of each array element
3060 c should be "float").
3061 f you are using AST_OUTLINER, the type of each array element
3062 f should be REAL).
3063 *
3064 * The storage order of data within this array should be such
3065 * that the index of the first grid dimension varies most
3066 * rapidly and that of the second dimension least rapidly
3067 c (i.e. Fortran array indexing is used).
3068 f (i.e. normal Fortran array storage order).
3069 c lbnd
3070 f LBND( 2 ) = INTEGER (Given)
3071 c Pointer to an array of two integers
3072 f An array
3073 * containing the coordinates of the centre of the first pixel
3074 * in the input grid along each dimension.
3075 c ubnd
3076 f UBND( 2) = INTEGER (Given)
3077 c Pointer to an array of two integers
3078 f An array
3079 * containing the coordinates of the centre of the last pixel in
3080 * the input grid along each dimension.
3081 *
3082 c Note that "lbnd" and "ubnd" together define the shape
3083 f Note that LBND and UBND together define the shape
3084 * and size of the input grid, its extent along a particular
3085 c (j'th) dimension being ubnd[j]-lbnd[j]+1 (assuming the
3086 c index "j" to be zero-based). They also define
3087 f (J'th) dimension being UBND(J)-LBND(J)+1. They also define
3088 * the input grid's coordinate system, each pixel having unit
3089 * extent along each dimension with integral coordinate values
3090 * at its centre or upper corner, as selected by parameter
3091 c "starpix".
3092 f STARPIX.
3093 c maxerr
3094 f MAXERR = DOUBLE PRECISION (Given)
3095 * Together with
3096 c "maxvert",
3097 f MAXVERT,
3098 * this determines how accurately the returned Polygon represents
3099 * the required region of the data array. It gives the target
3100 * discrepancy between the returned Polygon and the accurate outline
3101 * in the data array, expressed as a number of pixels. Insignificant
3102 * vertices are removed from the accurate outline, one by one, until
3103 * the number of vertices remaining in the returned Polygon equals
3104 c "maxvert",
3105 f MAXVERT,
3106 * or the largest discrepancy between the accurate outline and the
3107 * returned Polygon is greater than
3108 c "maxerr". If "maxerr"
3109 f MAXERR. If MAXERR
3110 * is zero or less, its value is ignored and the returned Polygon will
3111 * have the number of vertices specified by
3112 c "maxvert".
3113 f MAXVERT.
3114 c maxvert
3115 f MAXVERT = INTEGER (Given)
3116 * Together with
3117 c "maxerr",
3118 f MAXERR,
3119 * this determines how accurately the returned Polygon represents
3120 * the required region of the data array. It gives the maximum
3121 * allowed number of vertices in the returned Polygon. Insignificant
3122 * vertices are removed from the accurate outline, one by one, until
3123 * the number of vertices remaining in the returned Polygon equals
3124 c "maxvert",
3125 f MAXVERT,
3126 * or the largest discrepancy between the accurate outline and the
3127 * returned Polygon is greater than
3128 c "maxerr". If "maxvert"
3129 f MAXERR. If MAXVERT
3130 * is less than 3, its value is ignored and the number of vertices in
3131 * the returned Polygon will be the minimum needed to ensure that the
3132 * discrepancy between the accurate outline and the returned
3133 * Polygon is less than
3134 c "maxerr".
3135 f MAXERR.
3136 c inside
3137 f INSIDE( 2 ) = INTEGER (Given)
3138 c Pointer to an array of two integers
3139 f An array
3140 * containing the indices of a pixel known to be inside the
3141 * required region. This is needed because the supplied data
3142 * array may contain several disjoint areas of pixels that satisfy
3143 * the criterion specified by
3144 c "value" and "oper".
3145 f VALUE and OPER.
3146 * In such cases, the area described by the returned Polygon will
3147 * be the one that contains the pixel specified by
3148 c "inside".
3149 f INSIDE.
3150 * If the specified pixel is outside the bounds given by
3151 c "lbnd" and "ubnd",
3152 f LBND and UBND,
3153 * or has a value that does not meet the criterion specified by
3154 c "value" and "oper",
3155 f VALUE and OPER,
3156 * then this function will search for a suitable pixel. The search
3157 * starts at the central pixel and proceeds in a spiral manner until
3158 * a pixel is found that meets the specified crierion.
3159 c starpix
3160 f STARPIX = LOGICAL (Given)
3161 * A flag indicating the nature of the pixel coordinate system used
3162 * to describe the vertex positions in the returned Polygon. If
3163 c non-zero,
3164 f .TRUE.,
3165 * the standard Starlink definition of pixel coordinate is used in
3166 * which a pixel with integer index I spans a range of pixel coordinate
3167 * from (I-1) to I (i.e. pixel corners have integral pixel coordinates).
3168 c If zero,
3169 f If .FALSE.,
3170 * the definition of pixel coordinate used by other AST functions
3171 c such as astResample, astMask,
3172 f such as AST_RESAMPLE, AST_MASK,
3173 * etc., is used. In this definition, a pixel with integer index I
3174 * spans a range of pixel coordinate from (I-0.5) to (I+0.5) (i.e.
3175 * pixel centres have integral pixel coordinates).
3176 f STATUS = INTEGER (Given and Returned)
3177 f The global status.
3178
3179 * Returned Value:
3180 c astOutline<X>()
3181 f AST_OUTLINE<X> = INTEGER
3182 * A pointer to the required Polygon.
3183
3184 * Notes:
3185 * - This function proceeds by first finding a very accurate polygon,
3186 * and then removing insignificant vertices from this fine polygon
3187 * using
3188 c astDownsize.
3189 f AST_DOWNSIZE.
3190 * - The returned Polygon is the outer boundary of the contiguous set
3191 * of pixels that includes ths specified "inside" point, and satisfy
3192 * the specified value requirement. This set of pixels may potentially
3193 * include "holes" where the pixel values fail to meet the specified
3194 * value requirement. Such holes will be ignored by this function.
3195 c - NULL
3196 f - AST__NULL
3197 * will be returned if this function is invoked with the global
3198 * error status set, or if it should fail for any reason.
3199
3200 * Data Type Codes:
3201 * To select the appropriate masking function, you should
3202 c replace <X> in the generic function name astOutline<X> with a
3203 f replace <X> in the generic function name AST_OUTLINE<X> with a
3204 * 1- or 2-character data type code, so as to match the numerical
3205 * type <Xtype> of the data you are processing, as follows:
3206 c - D: double
3207 c - F: float
3208 c - L: long int
3209 c - UL: unsigned long int
3210 c - I: int
3211 c - UI: unsigned int
3212 c - S: short int
3213 c - US: unsigned short int
3214 c - B: byte (signed char)
3215 c - UB: unsigned byte (unsigned char)
3216 f - D: DOUBLE PRECISION
3217 f - R: REAL
3218 f - I: INTEGER
3219 f - UI: INTEGER (treated as unsigned)
3220 f - S: INTEGER*2 (short integer)
3221 f - US: INTEGER*2 (short integer, treated as unsigned)
3222 f - B: BYTE (treated as signed)
3223 f - UB: BYTE (treated as unsigned)
3224 *
3225 c For example, astOutlineD would be used to process "double"
3226 c data, while astOutlineS would be used to process "short int"
3227 c data, etc.
3228 f For example, AST_OUTLINED would be used to process DOUBLE
3229 f PRECISION data, while AST_OUTLINES would be used to process
3230 f short integer data (stored in an INTEGER*2 array), etc.
3231 f
3232 f For compatibility with other Starlink facilities, the codes W
3233 f and UW are provided as synonyms for S and US respectively (but
3234 f only in the Fortran interface to AST).
3235
3236 *--
3237 */
3238 /* Define a macro to implement the function for a specific data
3239 type. Note, this function cannot be a virtual function since the
3240 argument list does not include a Polygon, and so no virtual function
3241 table is available. */
3242 #define MAKE_OUTLINE(X,Xtype) \
3243 AstPolygon *astOutline##X##_( Xtype value, int oper, const Xtype array[], \
3244 const int lbnd[2], const int ubnd[2], double maxerr, \
3245 int maxvert, const int inside[2], int starpix, \
3246 int *status ) { \
3247 \
3248 /* Local Variables: */ \
3249 AstFrame *frm; /* Frame in which to define the Polygon */ \
3250 AstPointSet *candidate; /* Candidate polygon vertices */ \
3251 AstPointSet *pset; /* PointSet holding downsized polygon vertices */ \
3252 AstPolygon *result; /* Result value to return */ \
3253 const Xtype *pv; /* Pointer to next test point */ \
3254 Xtype v; /* Value of current pixel */ \
3255 double **ptr; /* Pointers to PointSet arrays */ \
3256 int boxsize; /* Half width of smoothign box in vertices */ \
3257 int inx; /* X index of inside point */ \
3258 int iny; /* Y index of inside point */ \
3259 int iv; /* Vector index of next test pixel */ \
3260 int ixv; /* X pixel index of next test point */ \
3261 int nv0; /* Number of vertices in accurate outline */ \
3262 int nx; /* Length of array x axis */ \
3263 int smooth; /* Do we need to smooth the polygon? */ \
3264 int stop_at_invalid; /* Indicates when to stop rightwards search */ \
3265 int tmp; /* Alternative boxsize */ \
3266 int valid; /* Does current pixel meet requirements? */ \
3267 static double junk[ 6 ] = {0.0, 0.0, 1.0, 1.0, 0.0, 1.0 }; /* Junk poly */ \
3268 \
3269 /* Initialise. */ \
3270 result = NULL; \
3271 \
3272 /* Check the global error status. */ \
3273 if ( !astOK ) return result; \
3274 \
3275 /* Avoid compiler warnings. */ \
3276 iv = 0; \
3277 \
3278 /* If we are going to be smoothing the polygon before downsizing it, we \
3279 need to ensure that the full polygon is retained within TraceEdge. if \
3280 this is not the case, TraceEdge can remove all vertices from straight \
3281 lines, except for the vertices that mark the beinning and end of the \
3282 straight line. */ \
3283 smooth = ( maxerr > 0.0 || maxvert > 2 ); \
3284 \
3285 /* Store the X dimension of the array. */ \
3286 nx = ubnd[ 0 ] - lbnd[ 0 ] + 1; \
3287 \
3288 /* See if a valid inside point was supplied. It must be inside the bounds \
3289 of the array, and must have a pixel value that meets the specified \
3290 requirement. */ \
3291 inx = inside[ 0 ]; \
3292 iny = inside[ 1 ]; \
3293 valid = ( inx >= lbnd[ 0 ] && inx <= ubnd[ 0 ] && \
3294 iny >= lbnd[ 1 ] && iny <= ubnd[ 1 ] ); \
3295 \
3296 if( valid ) { \
3297 iv = ( inx - lbnd[ 0 ] ) + (iny - lbnd[ 1 ] )*nx; \
3298 v = array[ iv ]; \
3299 \
3300 if( oper == AST__LT ) { \
3301 valid = ( v < value ); \
3302 \
3303 } else if( oper == AST__LE ) { \
3304 valid = ( v <= value ); \
3305 \
3306 } else if( oper == AST__EQ ) { \
3307 valid = ( v == value ); \
3308 \
3309 } else if( oper == AST__NE ) { \
3310 valid = ( v != value ); \
3311 \
3312 } else if( oper == AST__GE ) { \
3313 valid = ( v >= value ); \
3314 \
3315 } else if( oper == AST__GT ) { \
3316 valid = ( v > value ); \
3317 \
3318 } else if( astOK ){ \
3319 astError( AST__OPRIN, "astOutline"#X": Invalid operation code " \
3320 "(%d) supplied (programming error).", status, oper ); \
3321 } \
3322 } \
3323 \
3324 /* If no valid inside point was supplied, find one now. */ \
3325 if( !valid ) { \
3326 \
3327 if( oper == AST__LT ) { \
3328 FindInsidePointLT##X( value, array, lbnd, ubnd, &inx, &iny, &iv, status ); \
3329 \
3330 } else if( oper == AST__LE ) { \
3331 FindInsidePointLE##X( value, array, lbnd, ubnd, &inx, &iny, &iv, status ); \
3332 \
3333 } else if( oper == AST__EQ ) { \
3334 FindInsidePointEQ##X( value, array, lbnd, ubnd, &inx, &iny, &iv, status ); \
3335 \
3336 } else if( oper == AST__NE ) { \
3337 FindInsidePointNE##X( value, array, lbnd, ubnd, &inx, &iny, &iv, status ); \
3338 \
3339 } else if( oper == AST__GE ) { \
3340 FindInsidePointGE##X( value, array, lbnd, ubnd, &inx, &iny, &iv, status ); \
3341 \
3342 } else if( oper == AST__GT ) { \
3343 FindInsidePointGT##X( value, array, lbnd, ubnd, &inx, &iny, &iv, status ); \
3344 \
3345 } else if( astOK ){ \
3346 astError( AST__OPRIN, "astOutline"#X": Invalid operation code " \
3347 "(%d) supplied (programming error).", status, oper ); \
3348 } \
3349 } \
3350 \
3351 /* We now need to find a point on the boundary of the region containing \
3352 the inside point. Starting at the inside point, move to the right \
3353 through the array until a pixel is found which fails to meet the value \
3354 requirement or the edge of the array is reached. */ \
3355 \
3356 candidate = NULL; \
3357 pv = array + iv; \
3358 ixv = inx; \
3359 stop_at_invalid = 1; \
3360 \
3361 while( ++ixv <= ubnd[ 0 ] ) { \
3362 \
3363 /* Get the next pixel value. */ \
3364 v = *(++pv); \
3365 \
3366 /* See if it meets the value requirement. */ \
3367 if( oper == AST__LT ) { \
3368 valid = ( v < value ); \
3369 \
3370 } else if( oper == AST__LE ) { \
3371 valid = ( v <= value ); \
3372 \
3373 } else if( oper == AST__EQ ) { \
3374 valid = ( v == value ); \
3375 \
3376 } else if( oper == AST__NE ) { \
3377 valid = ( v != value ); \
3378 \
3379 } else if( oper == AST__GE ) { \
3380 valid = ( v >= value ); \
3381 \
3382 } else if( oper == AST__GT ) { \
3383 valid = ( v > value ); \
3384 \
3385 } else if( astOK ){ \
3386 astError( AST__OPRIN, "astOutline"#X": Invalid operation code " \
3387 "(%d) supplied (programming error).", status, oper ); \
3388 break; \
3389 } \
3390 \
3391 /* If we are currently looking for the next invalid pixel, and this pixel \
3392 is invalid... */ \
3393 if( stop_at_invalid ) { \
3394 if( ! valid ) { \
3395 \
3396 /* The current pixel may be on the required polygon, or it may be on the \
3397 edge of a hole contained within the region being outlined. We would \
3398 like to jump over such holes so that we can continue to look for the \
3399 real edge of the region being outlined. In order to determine if we \
3400 have reached a hole, we trace the edge that passes through the current \
3401 pixel, forming a candidate polygon in the process. In the process, We \
3402 see if the inside point falls within this candidate polygon. If it does \
3403 then the polygon is accepted as the required polygon. Otherwise, it is \
3404 rejected as a mere hole, and we continue moving away from the inside \
3405 point, looking for a new edge. */ \
3406 if( oper == AST__LT ) { \
3407 candidate = TraceEdgeLT##X( value, array, lbnd, ubnd, iv - 1, \
3408 ixv - 1, iny, starpix, smooth, status ); \
3409 \
3410 } else if( oper == AST__LE ) { \
3411 candidate = TraceEdgeLE##X( value, array, lbnd, ubnd, iv - 1, \
3412 ixv - 1, iny, starpix, smooth, status ); \
3413 \
3414 } else if( oper == AST__EQ ) { \
3415 candidate = TraceEdgeEQ##X( value, array, lbnd, ubnd, iv - 1, \
3416 ixv - 1, iny, starpix, smooth, status ); \
3417 \
3418 } else if( oper == AST__NE ) { \
3419 candidate = TraceEdgeNE##X( value, array, lbnd, ubnd, iv - 1, \
3420 ixv - 1, iny, starpix, smooth, status ); \
3421 \
3422 } else if( oper == AST__GE ) { \
3423 candidate = TraceEdgeGE##X( value, array, lbnd, ubnd, iv - 1, \
3424 ixv - 1, iny, starpix, smooth, status ); \
3425 \
3426 } else if( oper == AST__GT ) { \
3427 candidate = TraceEdgeGT##X( value, array, lbnd, ubnd, iv - 1, \
3428 ixv - 1, iny, starpix, smooth, status ); \
3429 \
3430 } else if( astOK ){ \
3431 astError( AST__OPRIN, "astOutline"#X": Invalid operation code " \
3432 "(%d) supplied (programming error).", status, oper ); \
3433 } \
3434 \
3435 /* If the candidate polygon is the required polygon, break out of the \
3436 loop. Otherwise, indicate that we want to continue moving right, \
3437 across the hole, until we reach the far side of the hole (i.e. find \
3438 the next valid pixel). */ \
3439 if( candidate ) { \
3440 break; \
3441 } else { \
3442 stop_at_invalid = 0; \
3443 } \
3444 } \
3445 \
3446 /* If we are currently looking for the next valid pixel, and the current \
3447 pixel is valid... */ \
3448 } else if( valid ) { \
3449 \
3450 /* We have reached the far side of a hole. Continue moving right, looking \
3451 now for the next invalid pixel (which may be on the required polygon). */ \
3452 stop_at_invalid = 1; \
3453 } \
3454 } \
3455 \
3456 /* If we have not yet found the required polygon, we must have reached \
3457 the right hand edge of the array. So we now follow the edge of the \
3458 array round until we meet the boundary of the required region. */ \
3459 if( !candidate ) { \
3460 if( oper == AST__LT ) { \
3461 candidate = TraceEdgeLT##X( value, array, lbnd, ubnd, iv - 1, \
3462 ixv - 1, iny, starpix, smooth, status ); \
3463 \
3464 } else if( oper == AST__LE ) { \
3465 candidate = TraceEdgeLE##X( value, array, lbnd, ubnd, iv - 1, \
3466 ixv - 1, iny, starpix, smooth, status ); \
3467 \
3468 } else if( oper == AST__EQ ) { \
3469 candidate = TraceEdgeEQ##X( value, array, lbnd, ubnd, iv - 1, \
3470 ixv - 1, iny, starpix, smooth, status ); \
3471 \
3472 } else if( oper == AST__NE ) { \
3473 candidate = TraceEdgeNE##X( value, array, lbnd, ubnd, iv - 1, \
3474 ixv - 1, iny, starpix, smooth, status ); \
3475 \
3476 } else if( oper == AST__GE ) { \
3477 candidate = TraceEdgeGE##X( value, array, lbnd, ubnd, iv - 1, \
3478 ixv - 1, iny, starpix, smooth, status ); \
3479 \
3480 } else if( oper == AST__GT ) { \
3481 candidate = TraceEdgeGT##X( value, array, lbnd, ubnd, iv - 1, \
3482 ixv - 1, iny, starpix, smooth, status ); \
3483 \
3484 } else if( astOK ){ \
3485 astError( AST__OPRIN, "astOutline"#X": Invalid operation code " \
3486 "(%d) supplied (programming error).", status, oper ); \
3487 } \
3488 } \
3489 \
3490 /* If required smooth the full resolution polygon before downsizing it. */ \
3491 if( smooth ) { \
3492 \
3493 /* Initially, set the boxsize to be equal to the required accouracy. */ \
3494 if( maxerr > 0 ) { \
3495 boxsize = (int) maxerr; \
3496 } else { \
3497 boxsize = INT_MAX; \
3498 } \
3499 \
3500 /* Determine a second box size equal to the average number of vertices in \
3501 the accurate outline, per vertex in the returned Polygon. */ \
3502 nv0 = astGetNpoint( candidate ); \
3503 if( maxvert > 2 ) { \
3504 tmp = nv0/(2*maxvert); \
3505 } else { \
3506 tmp = INT_MAX; \
3507 } \
3508 \
3509 /* Use the minimum of the two box sizes. */ \
3510 if( tmp < boxsize ) boxsize = tmp; \
3511 \
3512 /* Ensure the box is sufficiently small to allow at least 10 full boxes \
3513 (=20 half boxes) around the polygon. */ \
3514 tmp = nv0/20; \
3515 if( tmp < boxsize ) boxsize = tmp; \
3516 if( boxsize == 0 ) boxsize = 1; \
3517 \
3518 /* Smooth the polygon. */ \
3519 SmoothPoly( candidate, boxsize, 1.0, status ); \
3520 } \
3521 \
3522 /* Reduce the number of vertices in the outline. */ \
3523 frm = astFrame( 2, "Domain=PIXEL,Unit(1)=pixel,Unit(2)=pixel," \
3524 "Title=Pixel coordinates", status ); \
3525 pset = DownsizePoly( candidate, maxerr, maxvert, frm, status ); \
3526 \
3527 /* Create a default Polygon with 3 junk vertices. */ \
3528 result = astPolygon( frm, 3, 3, junk, NULL, " ", status ); \
3529 \
3530 /* Change the PointSet within the Polygon to the one created above. */ \
3531 SetPointSet( result, pset, status ); \
3532 \
3533 /* Free resources. Note, we need to free the arrays within the candidate \
3534 PointSet explicitly, since they were not created as part of the \
3535 construction of the PointSet (see TraceEdge). */ \
3536 pset = astAnnul( pset ); \
3537 frm = astAnnul( frm ); \
3538 ptr = astGetPoints( candidate ); \
3539 if( astOK ) { \
3540 astFree( ptr[ 0 ] ); \
3541 astFree( ptr[ 1 ] ); \
3542 } \
3543 candidate = astAnnul( candidate ); \
3544 \
3545 /* If an error occurred, clear the returned result. */ \
3546 if ( !astOK ) result = astAnnul( result ); \
3547 \
3548 /* Return the result. */ \
3549 return result; \
3550 }
3551
3552
3553 /* Expand the above macro to generate a function for each required
3554 data type. */
3555 #if HAVE_LONG_DOUBLE /* Not normally implemented */
MAKE_OUTLINE(LD,long double)3556 MAKE_OUTLINE(LD,long double)
3557 #endif
3558 MAKE_OUTLINE(D,double)
3559 MAKE_OUTLINE(L,long int)
3560 MAKE_OUTLINE(UL,unsigned long int)
3561 MAKE_OUTLINE(I,int)
3562 MAKE_OUTLINE(UI,unsigned int)
3563 MAKE_OUTLINE(S,short int)
3564 MAKE_OUTLINE(US,unsigned short int)
3565 MAKE_OUTLINE(B,signed char)
3566 MAKE_OUTLINE(UB,unsigned char)
3567 MAKE_OUTLINE(F,float)
3568
3569 /* Undefine the macros. */
3570 #undef MAKE_OUTLINE
3571
3572 /*
3573 * Name:
3574 * PartHull
3575
3576 * Purpose:
3577 * Find the convex hull enclosing selected pixels in one corner of a 2D array.
3578
3579 * Type:
3580 * Private function.
3581
3582 * Synopsis:
3583 * #include "polygon.h"
3584 * void PartHull<Oper><X>( <Xtype> value, const <Xtype> array[], int xdim,
3585 * int ydim, int xs, int ys, int xe, int ye,
3586 * int starpix, const int lbnd[2], double **xvert,
3587 * double **yvert, int *nvert, int *status )
3588
3589 * Class Membership:
3590 * Polygon member function
3591
3592 * Description:
3593 * This function uses an algorithm similar to "Andrew's Monotone Chain
3594 * Algorithm" to create a list of vertices describing one corner of the
3595 * convex hull enclosing the selected pixels in the supplied array.
3596 * The corner is defined to be the area of the array to the right of
3597 * the line from (xs,ys) to (xe,ye).
3598
3599 * Parameters:
3600 * value
3601 * A data value that specifies the pixels to be selected.
3602 * array
3603 * Pointer to a 2-dimensional array containing the data to be
3604 * processed. The numerical type of this array should match the 1-
3605 * or 2-character type code appended to the function name.
3606 * xdim
3607 * The number of pixels along each row of the array.
3608 * ydim
3609 * The number of rows in the array.
3610 * xs
3611 * The X GRID index of the first pixel on the line to be checked.
3612 * ys
3613 * The Y GRID index of the first pixel on the line to be checked.
3614 * xe
3615 * The X GRID index of the last pixel on the line to be checked.
3616 * ye
3617 * The Y GRID index of the last pixel on the line to be checked.
3618 * starpix
3619 * If non-zero, the usual Starlink definition of pixel coordinate
3620 * is used (integral values at pixel corners). Otherwise, the
3621 * system used by other AST functions such as astResample is used
3622 * (integral values at pixel centres).
3623 * lbnd
3624 * The lower pixel index bounds of the array.
3625 * xvert
3626 * Address of a pointer in which to return a pointer to the list
3627 * of GRID x values on the hull.
3628 * yvert
3629 * Address of a pointer in which to return a pointer to the list
3630 * of GRID y values on the hull.
3631 * nvert
3632 * Address of a pointer in which to return the number of points in
3633 * the returned xvert and yvert arrays.
3634 * status
3635 * Pointer to the inherited status variable.
3636
3637 */
3638
3639 /* Define a macro to implement the function for a specific data
3640 type and operation. */
3641 #define MAKE_PARTHULL(X,Xtype,Oper,OperI) \
3642 static void PartHull##Oper##X( Xtype value, const Xtype array[], int xdim, \
3643 int ydim, int xs, int ys, int xe, int ye, \
3644 int starpix, const int lbnd[2], \
3645 double **xvert, double **yvert, int *nvert, \
3646 int *status ) { \
3647 \
3648 /* Local Variables: */ \
3649 const Xtype *pc; \
3650 double *pxy; \
3651 double dx2; \
3652 double dx1; \
3653 double dy1; \
3654 double dy2; \
3655 double off; \
3656 double xdelta; \
3657 int ivert; \
3658 int ix; \
3659 int iy; \
3660 int x0; \
3661 int x1; \
3662 int xl; \
3663 int xlim; \
3664 int xr; \
3665 int yinc; \
3666 \
3667 /* Initialise */ \
3668 *yvert = NULL; \
3669 *xvert = NULL; \
3670 *nvert = 0; \
3671 \
3672 /* Check the global error status. */ \
3673 if ( !astOK ) return; \
3674 \
3675 /* If the line has zero length. just return a single vertex. */ \
3676 if( xs == xe && ys == ye ) { \
3677 *xvert = astMalloc( sizeof( double ) ); \
3678 *yvert = astMalloc( sizeof( double ) ); \
3679 if( astOK ) { \
3680 if( starpix ) { \
3681 (*xvert)[ 0 ] = xs + lbnd[ 0 ] - 1.5; \
3682 (*yvert)[ 0 ] = ys + lbnd[ 1 ] - 1.5; \
3683 } else { \
3684 (*xvert)[ 0 ] = xs + lbnd[ 0 ] - 1.0; \
3685 (*yvert)[ 0 ] = ys + lbnd[ 1 ] - 1.0; \
3686 } \
3687 *nvert = 1; \
3688 } \
3689 return; \
3690 } \
3691 \
3692 /* Otherwise check the line is sloping. */ \
3693 if( xs == xe ) { \
3694 astError( AST__INTER, "astOutline(Polygon): Bounding box " \
3695 "has zero width (internal AST programming error).", \
3696 status ); \
3697 return; \
3698 } else if( ys == ye ) { \
3699 astError( AST__INTER, "astOutline(Polygon): Bounding box " \
3700 "has zero height (internal AST programming error).", \
3701 status ); \
3702 return; \
3703 } \
3704 \
3705 /* Calculate the difference in length between adjacent rows of the area \
3706 to be tested. */ \
3707 xdelta = ((double)( xe - xs ))/((double)( ye - ys )); \
3708 \
3709 /* The left and right X limits */ \
3710 if( xe > xs ) { \
3711 xl = xs; \
3712 xr = xe; \
3713 } else { \
3714 xl = xe; \
3715 xr = xs; \
3716 } \
3717 \
3718 /* Get the increment in row number as we move from the start to the end \
3719 of the line. */ \
3720 yinc = ( ye > ys ) ? 1 : -1; \
3721 \
3722 /* Loop round all rows that cross the region to be tested, from start to \
3723 end of the supplied line. */ \
3724 iy = ys; \
3725 while( astOK ) { \
3726 \
3727 /* Get the GRID X coord where the line crosses this row. */ \
3728 xlim = (int)( 0.5 + xs + xdelta*( iy - ys ) ); \
3729 \
3730 /* Get the index of the first and last columns to be tested on this row. */ \
3731 if( yinc < 0 ) { \
3732 x0 = xl; \
3733 x1 = xlim; \
3734 } else { \
3735 x0 = xlim; \
3736 x1 = xr; \
3737 } \
3738 \
3739 /* Get a pointer to the first pixel to be tested at this row. */ \
3740 pc = array + ( iy - 1 )*xdim + x0 - 1; \
3741 \
3742 /* Loop round all columns to be tested in this row. */ \
3743 for( ix = x0; ix <= x1 && astOK; ix++,pc++ ) { \
3744 \
3745 /* Ignore pixels that are not selected. */ \
3746 if( ISVALID(*pc,OperI,value) ) { \
3747 \
3748 /* If this is the very first pixel, initialise the hull to contain just \
3749 the first pixel. */ \
3750 if( *nvert == 0 ){ \
3751 *xvert = astMalloc( 200*sizeof( double ) ); \
3752 *yvert = astMalloc( 200*sizeof( double ) ); \
3753 if( astOK ) { \
3754 (*xvert)[ 0 ] = ix; \
3755 (*yvert)[ 0 ] = iy; \
3756 *nvert = 1; \
3757 } else { \
3758 break; \
3759 } \
3760 \
3761 /* Otherwise.... */ \
3762 } else { \
3763 \
3764 /* Loop until the hull has been corrected to include the current pixel. */ \
3765 while( 1 ) { \
3766 \
3767 /* If the hull currently contains only one pixel, add the current pixel to \
3768 the end of the hull. */ \
3769 if( *nvert == 1 ){ \
3770 (*xvert)[ 1 ] = ix; \
3771 (*yvert)[ 1 ] = iy; \
3772 *nvert = 2; \
3773 break; \
3774 \
3775 /* Otherwise... */ \
3776 } else { \
3777 \
3778 /* Extend the line from the last-but-one pixel on the hull to the last \
3779 pixel on the hull, and see if the current pixel is to the left of \
3780 this line. If it is, it too is on the hull and so push it onto the end \
3781 of the list of vertices. */ \
3782 dx1 = (*xvert)[ *nvert - 1 ] - (*xvert)[ *nvert - 2 ]; \
3783 dy1 = (*yvert)[ *nvert - 1 ] - (*yvert)[ *nvert - 2 ]; \
3784 dx2 = ix - (*xvert)[ *nvert - 2 ]; \
3785 dy2 = iy - (*yvert)[ *nvert - 2 ]; \
3786 \
3787 if( dx1*dy2 > dx2*dy1 ) { \
3788 ivert = (*nvert)++; \
3789 *xvert = astGrow( *xvert, *nvert, sizeof( double ) ); \
3790 *yvert = astGrow( *yvert, *nvert, sizeof( double ) ); \
3791 if( astOK ) { \
3792 (*xvert)[ ivert ] = ix; \
3793 (*yvert)[ ivert ] = iy; \
3794 } \
3795 \
3796 /* Leave the loop now that the new point is on the hull. */ \
3797 break; \
3798 \
3799 /* If the new point is to the left of the line, then the last point \
3800 previously thought to be on hull is in fact not on the hull, so remove \
3801 it. We then loop again to compare the new pixel with modified hull. */ \
3802 } else { \
3803 (*nvert)--; \
3804 } \
3805 } \
3806 } \
3807 } \
3808 } \
3809 } \
3810 \
3811 if( iy == ye ) { \
3812 break; \
3813 } else { \
3814 iy += yinc; \
3815 } \
3816 \
3817 } \
3818 \
3819 /* Convert GRID coords to PIXEL coords. */ \
3820 if( astOK ) { \
3821 pxy = *xvert; \
3822 off = starpix ? lbnd[ 0 ] - 1.5 : lbnd[ 0 ] - 1.0; \
3823 for( ivert = 0; ivert < *nvert; ivert++ ) *(pxy++) += off; \
3824 \
3825 pxy = *yvert; \
3826 off = starpix ? lbnd[ 1 ] - 1.5 : lbnd[ 1 ] - 1.0; \
3827 for( ivert = 0; ivert < *nvert; ivert++ ) *(pxy++) += off; \
3828 \
3829 /* Free lists if an error has occurred. */ \
3830 } else { \
3831 *xvert = astFree( *xvert ); \
3832 *yvert = astFree( *yvert ); \
3833 *nvert = 0; \
3834 } \
3835 }
3836
3837 /* Define a macro that uses the above macro to to create implementations
3838 of PartHull for all operations. */
3839 #define MAKEALL_PARTHULL(X,Xtype) \
3840 MAKE_PARTHULL(X,Xtype,LT,AST__LT) \
3841 MAKE_PARTHULL(X,Xtype,LE,AST__LE) \
3842 MAKE_PARTHULL(X,Xtype,EQ,AST__EQ) \
3843 MAKE_PARTHULL(X,Xtype,NE,AST__NE) \
3844 MAKE_PARTHULL(X,Xtype,GE,AST__GE) \
3845 MAKE_PARTHULL(X,Xtype,GT,AST__GT)
3846
3847 /* Expand the above macro to generate a function for each required
3848 data type and operation. */
3849 #if HAVE_LONG_DOUBLE /* Not normally implemented */
3850 MAKEALL_PARTHULL(LD,long double)
3851 #endif
3852 MAKEALL_PARTHULL(D,double)
3853 MAKEALL_PARTHULL(L,long int)
3854 MAKEALL_PARTHULL(UL,unsigned long int)
3855 MAKEALL_PARTHULL(I,int)
3856 MAKEALL_PARTHULL(UI,unsigned int)
3857 MAKEALL_PARTHULL(S,short int)
3858 MAKEALL_PARTHULL(US,unsigned short int)
3859 MAKEALL_PARTHULL(B,signed char)
3860 MAKEALL_PARTHULL(UB,unsigned char)
3861 MAKEALL_PARTHULL(F,float)
3862
3863 /* Undefine the macros. */
3864 #undef MAKE_PARTHULL
3865 #undef MAKEALL_PARTHULL
3866
3867 static double Polywidth( AstFrame *frm, AstLineDef **edges, int i, int nv,
3868 double cen[ 2 ], int *status ){
3869 /*
3870 * Name:
3871 * Polywidth
3872
3873 * Purpose:
3874 * Find the width of a polygon perpendicular to a given edge.
3875
3876 * Type:
3877 * Private function.
3878
3879 * Synopsis:
3880 * #include "polygon.h"
3881 * double Polywidth( AstFrame *frm, AstLineDef **edges, int i, int nv,
3882 * double cen[ 2 ], int *status )
3883
3884 * Class Membership:
3885 * Polygon member function
3886
3887 * Description:
3888 * This function defines a line perpendicular to a given polygon edge,
3889 * passing through the mid point of the edge, extending towards the
3890 * inside of the polygon. It returns the distance that can be travelled
3891 * along this line before any of the other polygon edges are hit (the
3892 * "width" of the polygon perpendicular to the given edge). It also
3893 * puts the position corresponding to half that distance into "cen".
3894
3895 * Parameters:
3896 * frm
3897 * The Frame in which the lines are defined.
3898 * edges
3899 * Array of "nv" pointers to AstLineDef structures, each defining an
3900 * edge of the polygon.
3901 * i
3902 * The index of the edge that is to define the polygon width.
3903 * nv
3904 * Total number of edges.
3905 * cen
3906 * An array into which are put the coords of the point half way
3907 * along the polygon width line.
3908 * status
3909 * Pointer to the inherited status variable.
3910
3911 * Returnd Value:
3912 * The width of the polygon perpendicular to the given edge, or
3913 * AST__BAD if the width cannot be determined (usually because the
3914 * vertices been supplied in a clockwise direction, effectively
3915 * negating the Polygon).
3916
3917 */
3918
3919 /* Local Variables: */
3920 AstLineDef *line;
3921 double *cross;
3922 double d;
3923 double end[ 2 ];
3924 double l1;
3925 double l2;
3926 double result;
3927 double start[ 2 ];
3928 int j;
3929
3930 /* Check the global error status. */
3931 result = AST__BAD;
3932 if ( !astOK ) return result;
3933
3934 /* Create a Line description for a line perpendicular to the specified
3935 edge, passing through the mid point of the edge, and extending towards
3936 the inside of the polygon. First move away from the start along the
3937 line to the mid point. This gives the start of the new line. */
3938 l1 = 0.5*( edges[ i ]->length );
3939 astLineOffset( frm, edges[ i ], l1, 0.0, start );
3940
3941 /* We now move away from this position at right angles to the line. We
3942 start off by moving 5 times the length of the specified edge. For
3943 some Frames (e.g. SkyFrames) this may result in a position that is
3944 much too close (i.e. if it goes all the way round the great circle
3945 and comes back to the beginning). Therefore, we check that the end
3946 point is the requested distance from the start point, and if not, we
3947 halve the length of the line and try again. */
3948 l2 = 10.0*l1;
3949 while( 1 ) {
3950 astLineOffset( frm, edges[ i ], l1, l2, end );
3951 d = astDistance( frm, start, end );
3952 if( d != AST__BAD && fabs( d - l2 ) < 1.0E-6*l2 ) break;
3953 l2 *= 0.5;
3954 }
3955
3956 /* Create a description of the required line. */
3957 line = astLineDef( frm, start, end );
3958
3959 /* Loop round every edge, except for the supplied edge. */
3960 for( j = 0; j < nv; j++ ) {
3961 if( j != i ) {
3962
3963 /* Find the position at which the line created above crosses the current
3964 edge. Skip to the next edge if the line does not intersect the edge
3965 within the length of the edge. */
3966 if( astLineCrossing( frm, line, edges[ j ], &cross ) ) {
3967
3968 /* Find the distance between the crossing point and the line start. */
3969 d = astDistance( frm, start, cross );
3970
3971 /* If this is less than the smallest found so far, record it. */
3972 if( d != AST__BAD && ( d < result || result == AST__BAD ) ) {
3973 result = d;
3974 }
3975 }
3976
3977 /* Free resources */
3978 cross = astFree( cross );
3979 }
3980 }
3981 line = astFree( line );
3982
3983 /* If a width was found, return the point half way across the polygon. */
3984 if( result != AST__BAD ) {
3985 astOffset( frm, start, end, 0.5*result, cen );
3986
3987 /* The usual reason for not finding a width is if the polygon vertices
3988 are supplied in clockwise order, effectively negating the polygon, and
3989 resulting in the "inside" of the polygon being the infinite region
3990 outside a polygonal hole. In this case, the end point of the line
3991 perpendicular to the initial edge can be returned as a representative
3992 "inside" point. */
3993 } else {
3994 cen[ 0 ] = end[ 0 ];
3995 cen[ 1 ] = end[ 1 ];
3996 }
3997
3998 /* Return the width. */
3999 return result;
4000
4001 }
4002
RegBaseBox(AstRegion * this_region,double * lbnd,double * ubnd,int * status)4003 static void RegBaseBox( AstRegion *this_region, double *lbnd, double *ubnd, int *status ){
4004 /*
4005 * Name:
4006 * RegBaseBox
4007
4008 * Purpose:
4009 * Returns the bounding box of an un-negated Region in the base Frame of
4010 * the encapsulated FrameSet.
4011
4012 * Type:
4013 * Private function.
4014
4015 * Synopsis:
4016 * #include "polygon.h"
4017 * void RegBaseBox( AstRegion *this, double *lbnd, double *ubnd, int *status )
4018
4019 * Class Membership:
4020 * Polygon member function (over-rides the astRegBaseBox protected
4021 * method inherited from the Region class).
4022
4023 * Description:
4024 * This function returns the upper and lower axis bounds of a Region in
4025 * the base Frame of the encapsulated FrameSet, assuming the Region
4026 * has not been negated. That is, the value of the Negated attribute
4027 * is ignored.
4028
4029 * Parameters:
4030 * this
4031 * Pointer to the Region.
4032 * lbnd
4033 * Pointer to an array in which to return the lower axis bounds
4034 * covered by the Region in the base Frame of the encapsulated
4035 * FrameSet. It should have at least as many elements as there are
4036 * axes in the base Frame.
4037 * ubnd
4038 * Pointer to an array in which to return the upper axis bounds
4039 * covered by the Region in the base Frame of the encapsulated
4040 * FrameSet. It should have at least as many elements as there are
4041 * axes in the base Frame.
4042 * status
4043 * Pointer to the inherited status variable.
4044
4045 */
4046
4047 /* Local Variables: */
4048 AstFrame *frm; /* Pointer to encapsulated Frame */
4049 AstPointSet *pset; /* Pointer to PointSet defining the Region */
4050 AstPolygon *this; /* Pointer to Polygon structure */
4051 AstRegion *reg; /* Base Frame equivalent of supplied Polygon */
4052 double **ptr; /* Pointer to PointSet data */
4053 double *x; /* Pointer to next X axis value */
4054 double *y; /* Pointer to next Y axis value */
4055 double dist; /* Offset along an axis */
4056 double x0; /* The first X axis value */
4057 double y0; /* The first Y axis value */
4058 int ip; /* Point index */
4059 int np; /* No. of points in PointSet */
4060
4061 /* Check the global error status. */
4062 if ( !astOK ) return;
4063
4064 /* Get a pointer to the Polygon structure. */
4065 this = (AstPolygon *) this_region;
4066
4067 /* If the base Frame bounding box has already been found, return the
4068 values stored in the Polygon structure. */
4069 if( this->lbnd[ 0 ] != AST__BAD ) {
4070 lbnd[ 0 ] = this->lbnd[ 0 ];
4071 lbnd[ 1 ] = this->lbnd[ 1 ];
4072 ubnd[ 0 ] = this->ubnd[ 0 ];
4073 ubnd[ 1 ] = this->ubnd[ 1 ];
4074
4075 /* If the base Frame bounding box has not yet been found, find it now and
4076 store it in the Polygon structure. */
4077 } else {
4078
4079 /* Get the axis values for the PointSet which defines the location and
4080 extent of the region in the base Frame of the encapsulated FrameSet. */
4081 pset = this_region->points;
4082 ptr = astGetPoints( pset );
4083 np = astGetNpoint( pset );
4084
4085 /* Get a pointer to the base Frame in the frameset encapsulated by the
4086 parent Region structure. */
4087 frm = astGetFrame( this_region->frameset, AST__BASE );
4088
4089 /* Find the upper and lower bounds of the box enclosing all the vertices.
4090 The box is expressed in terms of axis offsets from the first vertex, in
4091 order to avoid problems with boxes that cross RA=0 or RA=12h */
4092 lbnd[ 0 ] = 0.0;
4093 lbnd[ 1 ] = 0.0;
4094 ubnd[ 0 ] = 0.0;
4095 ubnd[ 1 ] = 0.0;
4096
4097 x = ptr[ 0 ];
4098 y = ptr[ 1 ];
4099
4100 x0 = *x;
4101 y0 = *y;
4102
4103 for( ip = 0; ip < np; ip++, x++, y++ ) {
4104
4105 dist = astAxDistance( frm, 1, x0, *x );
4106 if( dist < lbnd[ 0 ] ) {
4107 lbnd[ 0 ] = dist;
4108 } else if( dist > ubnd[ 0 ] ) {
4109 ubnd[ 0 ] = dist;
4110 }
4111
4112 dist = astAxDistance( frm, 2, y0, *y );
4113 if( dist < lbnd[ 1 ] ) {
4114 lbnd[ 1 ] = dist;
4115 } else if( dist > ubnd[ 1 ] ) {
4116 ubnd[ 1 ] = dist;
4117 }
4118
4119 }
4120
4121 /* Convert the box bounds to absolute values, rather than values relative
4122 to the first vertex. */
4123 lbnd[ 0 ] += x0;
4124 lbnd[ 1 ] += y0;
4125 ubnd[ 0 ] += x0;
4126 ubnd[ 1 ] += y0;
4127
4128 /* The astNormBox requires a Mapping which can be used to test points in
4129 this base Frame. Create a copy of the Polygon and then set its
4130 FrameSet so that the current Frame in the copy is the same as the base
4131 Frame in the original. */
4132 reg = astCopy( this );
4133 astSetRegFS( reg, frm );
4134 astSetNegated( reg, 0 );
4135
4136 /* Normalise this box. */
4137 astNormBox( frm, lbnd, ubnd, reg );
4138
4139 /* Free resources */
4140 reg = astAnnul( reg );
4141 frm = astAnnul( frm );
4142
4143 /* Store it in the olygon structure for future use. */
4144 this->lbnd[ 0 ] = lbnd[ 0 ];
4145 this->lbnd[ 1 ] = lbnd[ 1 ];
4146 this->ubnd[ 0 ] = ubnd[ 0 ];
4147 this->ubnd[ 1 ] = ubnd[ 1 ];
4148
4149 }
4150 }
4151
RegBaseMesh(AstRegion * this_region,int * status)4152 static AstPointSet *RegBaseMesh( AstRegion *this_region, int *status ){
4153 /*
4154 * Name:
4155 * RegBaseMesh
4156
4157 * Purpose:
4158 * Return a PointSet containing a mesh of points on the boundary of a
4159 * Region in its base Frame.
4160
4161 * Type:
4162 * Private function.
4163
4164 * Synopsis:
4165 * #include "polygon.h"
4166 * AstPointSet *astRegBaseMesh( AstRegion *this, int *status )
4167
4168 * Class Membership:
4169 * Polygon member function (over-rides the astRegBaseMesh protected
4170 * method inherited from the Region class).
4171
4172 * Description:
4173 * This function returns a PointSet containing a mesh of points on the
4174 * boundary of the Region. The points refer to the base Frame of
4175 * the encapsulated FrameSet.
4176
4177 * Parameters:
4178 * this
4179 * Pointer to the Region.
4180 * status
4181 * Pointer to the inherited status variable.
4182
4183 * Returned Value:
4184 * Pointer to the PointSet. Annul the pointer using astAnnul when it
4185 * is no longer needed.
4186
4187 * Notes:
4188 * - A NULL pointer is returned if an error has already occurred, or if
4189 * this function should fail for any reason.
4190
4191 */
4192
4193 /* Local Variables: */
4194 AstFrame *frm; /* Base Frame in encapsulated FrameSet */
4195 AstPointSet *result; /* Returned pointer */
4196 AstPolygon *this; /* The Polygon structure */
4197 double **rptr; /* Pointers to returned mesh data */
4198 double **vptr; /* Pointers to vertex data */
4199 double *lens; /* Pointer to work space holding edge lengths */
4200 double d; /* Length of this edge */
4201 double delta; /* Angular separation of points */
4202 double end[ 2 ]; /* End position */
4203 double mp; /* No. of mesh points per unit distance */
4204 double p[ 2 ]; /* Position in 2D Frame */
4205 double start[ 2 ]; /* Start position */
4206 double total; /* Total length of polygon */
4207 int ip; /* Point index */
4208 int iv; /* Vertex index */
4209 int n; /* No. of points on this edge */
4210 int next; /* Index of next point in returned PointSet */
4211 int np; /* No. of points in returned PointSet */
4212 int nv; /* No. of polygon vertices */
4213
4214 /* Initialise */
4215 result= NULL;
4216
4217 /* Check the global error status. */
4218 if ( !astOK ) return result;
4219
4220 /* If the Region structure contains a pointer to a PointSet holding
4221 a previously created mesh, return it. */
4222 if( this_region->basemesh ) {
4223 result = astClone( this_region->basemesh );
4224
4225 /* Otherwise, create a new mesh. */
4226 } else {
4227
4228 /* Get a pointer to the Polygon structure. */
4229 this = (AstPolygon *) this_region;
4230
4231 /* Get a pointer to the base Frame in the encapsulated FrameSet. */
4232 frm = astGetFrame( this_region->frameset, AST__BASE );
4233
4234 /* Get the number of vertices and pointers to the vertex axis values. */
4235 nv = astGetNpoint( this_region->points );
4236 vptr = astGetPoints( this_region->points );
4237
4238 /* Allocate memory to hold the geodesic length of each edge. */
4239 lens = astMalloc( sizeof( double )*(size_t) nv );
4240
4241 if( astOK ) {
4242
4243 /* Find the total geodesic distance around the boundary. */
4244 total = 0.0;
4245
4246 start[ 0 ] = vptr[ 0 ][ 0 ];
4247 start[ 1 ] = vptr[ 1 ][ 0 ];
4248
4249 for( iv = 1; iv < nv; iv++ ) {
4250 end[ 0 ] = vptr[ 0 ][ iv ];
4251 end[ 1 ] = vptr[ 1 ][ iv ];
4252
4253 d = astDistance( frm, start, end );
4254 if( d != AST__BAD ) total += fabs( d );
4255 lens[ iv ] = d;
4256 start[ 0 ] = end[ 0 ];
4257 start[ 1 ] = end[ 1 ];
4258 }
4259
4260 end[ 0 ] = vptr[ 0 ][ 0 ];
4261 end[ 1 ] = vptr[ 1 ][ 0 ];
4262
4263 d = astDistance( frm, start, end );
4264 if( d != AST__BAD ) total += fabs( d );
4265 lens[ 0 ] = d;
4266
4267 /* Find the number of mesh points per unit geodesic distance. */
4268 if( total > 0.0 ){
4269 mp = astGetMeshSize( this )/total;
4270
4271 /* Find the total number of mesh points required. */
4272 np = 0;
4273 for( iv = 0; iv < nv; iv++ ) {
4274 if( lens[ iv ] != AST__BAD ) np += 1 + (int)( lens[ iv ]*mp );
4275 }
4276
4277 /* Create a suitable PointSet to hold the returned positions. */
4278 result = astPointSet( np, 2, "", status );
4279 rptr = astGetPoints( result );
4280 if( astOK ) {
4281
4282 /* Initialise the index of the next point to be added to the returned
4283 PointSet. */
4284 next = 0;
4285
4286 /* Loop round each good edge of the polygon. The edge ends at vertex "iv". */
4287 start[ 0 ] = vptr[ 0 ][ 0 ];
4288 start[ 1 ] = vptr[ 1 ][ 0 ];
4289
4290 for( iv = 1; iv < nv; iv++ ) {
4291 end[ 0 ] = vptr[ 0 ][ iv ];
4292 end[ 1 ] = vptr[ 1 ][ iv ];
4293 if( lens[ iv ] != AST__BAD ) {
4294
4295 /* Add the position of the starting vertex to the returned mesh. */
4296 rptr[ 0 ][ next ] = start[ 0 ];
4297 rptr[ 1 ][ next ] = start[ 1 ];
4298 next++;
4299
4300 /* Find the number of points on this edge, and the geodesic distance
4301 between them. */
4302 n = 1 + (int) ( lens[ iv ]*mp );
4303
4304 /* If more than one point, find the distance between points. */
4305 if( n > 1 ) {
4306 delta = lens[ iv ]/n;
4307
4308 /* Loop round the extra points. */
4309 for( ip = 1; ip < n; ip++ ) {
4310
4311 /* Find the position of the next mesh point. */
4312 astOffset( frm, start, end, delta*ip, p );
4313
4314 /* Add it to the mesh. */
4315 rptr[ 0 ][ next ] = p[ 0 ];
4316 rptr[ 1 ][ next ] = p[ 1 ];
4317 next++;
4318
4319 }
4320 }
4321 }
4322
4323 /* The end of this edge becomes the start of the next. */
4324 start[ 0 ] = end[ 0 ];
4325 start[ 1 ] = end[ 1 ];
4326 }
4327
4328 /* Now do the edge which ends at the first vertex. */
4329 end[ 0 ] = vptr[ 0 ][ 0 ];
4330 end[ 1 ] = vptr[ 1 ][ 0 ];
4331 if( lens[ 0 ] != AST__BAD ) {
4332 rptr[ 0 ][ next ] = start[ 0 ];
4333 rptr[ 1 ][ next ] = start[ 1 ];
4334 next++;
4335
4336 n = 1 + (int)( lens[ 0 ]*mp );
4337 if( n > 1 ) {
4338 delta = lens[ 0 ]/n;
4339 for( ip = 1; ip < n; ip++ ) {
4340 astOffset( frm, start, end, delta*ip, p );
4341 rptr[ 0 ][ next ] = p[ 0 ];
4342 rptr[ 1 ][ next ] = p[ 1 ];
4343 next++;
4344 }
4345 }
4346 }
4347
4348 /* Check the PointSet size was correct. */
4349 if( next != np && astOK ) {
4350 astError( AST__INTER, "astRegBaseMesh(%s): Error in the "
4351 "allocated PointSet size (%d) - should have "
4352 "been %d (internal AST programming error).", status,
4353 astGetClass( this ), np, next );
4354 }
4355
4356 /* Save the returned pointer in the Region structure so that it does not
4357 need to be created again next time this function is called. */
4358 if( astOK ) this_region->basemesh = astClone( result );
4359
4360 }
4361
4362 } else if( astOK ) {
4363 astError( AST__BADIN, "astRegBaseMesh(%s): The boundary of "
4364 "the supplied %s has an undefined length.", status,
4365 astGetClass( this ), astGetClass( this ) );
4366 }
4367
4368 }
4369
4370 /* Free resources. */
4371 frm = astAnnul( frm );
4372 lens = astFree( lens );
4373 }
4374
4375 /* Annul the result if an error has occurred. */
4376 if( !astOK ) result = astAnnul( result );
4377
4378 /* Return a pointer to the output PointSet. */
4379 return result;
4380 }
4381
RegPins(AstRegion * this_region,AstPointSet * pset,AstRegion * unc,int ** mask,int * status)4382 static int RegPins( AstRegion *this_region, AstPointSet *pset, AstRegion *unc,
4383 int **mask, int *status ){
4384 /*
4385 * Name:
4386 * RegPins
4387
4388 * Purpose:
4389 * Check if a set of points fall on the boundary of a given Polygon.
4390
4391 * Type:
4392 * Private function.
4393
4394 * Synopsis:
4395 * #include "polygon.h"
4396 * int RegPins( AstRegion *this, AstPointSet *pset, AstRegion *unc,
4397 * int **mask, int *status )
4398
4399 * Class Membership:
4400 * Polygon member function (over-rides the astRegPins protected
4401 * method inherited from the Region class).
4402
4403 * Description:
4404 * This function returns a flag indicating if the supplied set of
4405 * points all fall on the boundary of the given Polygon.
4406 *
4407 * Some tolerance is allowed, as specified by the uncertainty Region
4408 * stored in the supplied Polygon "this", and the supplied uncertainty
4409 * Region "unc" which describes the uncertainty of the supplied points.
4410
4411 * Parameters:
4412 * this
4413 * Pointer to the Polygon.
4414 * pset
4415 * Pointer to the PointSet. The points are assumed to refer to the
4416 * base Frame of the FrameSet encapsulated by "this".
4417 * unc
4418 * Pointer to a Region representing the uncertainties in the points
4419 * given by "pset". The Region is assumed to represent the base Frame
4420 * of the FrameSet encapsulated by "this". Zero uncertainity is assumed
4421 * if NULL is supplied.
4422 * mask
4423 * Pointer to location at which to return a pointer to a newly
4424 * allocated dynamic array of ints. The number of elements in this
4425 * array is equal to the value of the Npoint attribute of "pset".
4426 * Each element in the returned array is set to 1 if the
4427 * corresponding position in "pset" is on the boundary of the Region
4428 * and is set to zero otherwise. A NULL value may be supplied
4429 * in which case no array is created. If created, the array should
4430 * be freed using astFree when no longer needed.
4431 * status
4432 * Pointer to the inherited status variable.
4433
4434 * Returned Value:
4435 * Non-zero if the points all fall on the boundary of the given
4436 * Region, to within the tolerance specified. Zero otherwise.
4437
4438 */
4439
4440 /* Local variables: */
4441 AstFrame *frm; /* Base Frame in supplied Polygon */
4442 AstPointSet *pset1; /* Pointer to copy of supplied PointSet */
4443 AstPointSet *pset2; /* Pointer to PointSet holding resolved components */
4444 AstPolygon *this; /* Pointer to the Polygon structure. */
4445 AstRegion *tunc; /* Uncertainity Region from "this" */
4446 double **ptr1; /* Pointer to axis values in "pset1" */
4447 double **ptr2; /* Pointer to axis values in "pset2" */
4448 double **vptr; /* Pointer to axis values at vertices */
4449 double *safe; /* An interior point in "this" */
4450 double edge_len; /* Length of current edge */
4451 double end[2]; /* Position of end of current edge */
4452 double l1; /* Length of bounding box diagonal */
4453 double l2; /* Length of bounding box diagonal */
4454 double lbnd_tunc[2]; /* Lower bounds of "this" uncertainty Region */
4455 double lbnd_unc[2]; /* Lower bounds of supplied uncertainty Region */
4456 double par; /* Parallel component */
4457 double parmax; /* Max acceptable parallel component */
4458 double prp; /* Perpendicular component */
4459 double start[2]; /* Position of start of current edge */
4460 double ubnd_tunc[2]; /* Upper bounds of "this" uncertainty Region */
4461 double ubnd_unc[2]; /* Upper bounds of supplied uncertainty Region */
4462 double wid; /* Width of acceptable margin around polygon */
4463 int *m; /* Pointer to next mask value */
4464 int i; /* Edge index */
4465 int ip; /* Point index */
4466 int np; /* No. of supplied points */
4467 int nv; /* No. of vertices */
4468 int result; /* Returned flag */
4469
4470 /* Initialise */
4471 result = 0;
4472 if( mask ) *mask = NULL;
4473
4474 /* Check the inherited status. */
4475 if( !astOK ) return result;
4476
4477 /* Get a pointer to the Polygon structure. */
4478 this = (AstPolygon *) this_region;
4479
4480 /* Check the supplied PointSet has 2 axis values per point. */
4481 if( astGetNcoord( pset ) != 2 && astOK ) {
4482 astError( AST__INTER, "astRegPins(%s): Illegal number of axis "
4483 "values per point (%d) in the supplied PointSet - should be "
4484 "2 (internal AST programming error).", status, astGetClass( this ),
4485 astGetNcoord( pset ) );
4486 }
4487
4488 /* Get the number of axes in the uncertainty Region and check it is also 2. */
4489 if( unc && astGetNaxes( unc ) != 2 && astOK ) {
4490 astError( AST__INTER, "astRegPins(%s): Illegal number of axes (%d) "
4491 "in the supplied uncertainty Region - should be 2 "
4492 "(internal AST programming error).", status, astGetClass( this ),
4493 astGetNaxes( unc ) );
4494 }
4495
4496 /* Get pointers to the axis values at the polygon vertices. */
4497 vptr = astGetPoints( this_region->points );
4498
4499 /* Get the number of vertices. */
4500 nv = astGetNpoint( this_region->points );
4501
4502 /* Take a copy of the supplied PointSet and get pointers to its axis
4503 values,and its size */
4504 pset1 = astCopy( pset );
4505 ptr1 = astGetPoints( pset1 );
4506 np = astGetNpoint( pset1 );
4507
4508 /* Create a PointSet to hold the resolved components and get pointers to its
4509 axis data. */
4510 pset2 = astPointSet( np, 2, "", status );
4511 ptr2 = astGetPoints( pset2 );
4512
4513 /* Create a mask array if required. */
4514 if( mask ) *mask = astMalloc( sizeof(int)*(size_t) np );
4515
4516 /* Get the centre of the region in the base Frame. We use this as a "safe"
4517 interior point within the region. */
4518 safe = astRegCentre( this, NULL, NULL, 0, AST__BASE );
4519
4520 /* We now find the maximum distance on each axis that a point can be from the
4521 boundary of the Polygon for it still to be considered to be on the boundary.
4522 First get the Region which defines the uncertainty within the Polygon
4523 being checked (in its base Frame), re-centre it on the interior point
4524 found above (to avoid problems if the uncertainty region straddles a
4525 discontinuity), and get its bounding box. The current Frame of the
4526 uncertainty Region is the same as the base Frame of the Polygon. */
4527 tunc = astGetUncFrm( this, AST__BASE );
4528 if( safe ) astRegCentre( tunc, safe, NULL, 0, AST__CURRENT );
4529 astGetRegionBounds( tunc, lbnd_tunc, ubnd_tunc );
4530
4531 /* Find the geodesic length within the base Frame of "this" of the diagonal of
4532 the bounding box. */
4533 frm = astGetFrame( this_region->frameset, AST__BASE );
4534 l1 = astDistance( frm, lbnd_tunc, ubnd_tunc );
4535
4536 /* Also get the Region which defines the uncertainty of the supplied
4537 points and get its bounding box. First re-centre the uncertainty at the
4538 interior position to avoid problems from uncertainties that straddle a
4539 discontinuity. */
4540 if( unc ) {
4541 if( safe ) astRegCentre( unc, safe, NULL, 0, AST__CURRENT );
4542 astGetRegionBounds( unc, lbnd_unc, ubnd_unc );
4543
4544 /* Find the geodesic length of the diagonal of this bounding box. */
4545 l2 = astDistance( frm, lbnd_unc, ubnd_unc );
4546
4547 /* Assume zero uncertainty if no "unc" Region was supplied. */
4548 } else {
4549 l2 = 0.0;
4550 }
4551
4552 /* The required border width is half of the total diagonal of the two bounding
4553 boxes. */
4554 if( astOK ) {
4555 wid = 0.5*( l1 + l2 );
4556
4557 /* Loop round each edge of the polygon. Edge "i" starts at vertex "i-1"
4558 and ends at vertex "i". Edge zero starts at vertex "nv-1" and ends at
4559 vertex zero. */
4560 start[ 0 ] = vptr[ 0 ][ nv - 1 ];
4561 start[ 1 ] = vptr[ 1 ][ nv - 1 ];
4562 for( i = 0; i < nv; i++ ) {
4563 end[ 0 ] = vptr[ 0 ][ i ];
4564 end[ 1 ] = vptr[ 1 ][ i ];
4565
4566 /* Find the length of this edge. */
4567 edge_len = astDistance( frm, start, end );
4568
4569 /* Resolve all the supplied mesh points into components parallel and
4570 perpendicular to this edge. */
4571 (void) astResolvePoints( frm, start, end, pset1, pset2 );
4572
4573 /* A point is effectively on this edge if the parallel component is
4574 greater than (-wid) and less than (edge_len+wid) AND the perpendicular
4575 component has an absolute value less than wid. Identify such positions
4576 and set them bad in pset1. */
4577 parmax = edge_len + wid;
4578 for( ip = 0; ip < np; ip++ ) {
4579 par = ptr2[ 0 ][ ip ];
4580 prp = ptr2[ 1 ][ ip ];
4581
4582 if( par != AST__BAD && prp != AST__BAD ) {
4583 if( par > -wid && par < parmax && prp > -wid && prp < wid ) {
4584 ptr1[ 0 ][ ip ] = AST__BAD;
4585 ptr1[ 1 ][ ip ] = AST__BAD;
4586 }
4587 }
4588 }
4589
4590 /* The end of the current edge becomes the start of the next. */
4591 start[ 0 ] = end[ 0 ];
4592 start[ 1 ] = end[ 1 ];
4593 }
4594
4595 /* See if any good points are left in pset1. If so, it means that those
4596 points were not on any edge of hte Polygon. We use two alogorithms
4597 here depending on whether we are creating a mask array, since we can
4598 abort the check upon finding the first good point if we are not
4599 producing a mask. */
4600 result = 1;
4601 if( mask ) {
4602 m = *mask;
4603 for( ip = 0; ip < np; ip++, m++ ) {
4604 if( ptr1[ 0 ][ ip ] != AST__BAD &&
4605 ptr1[ 1 ][ ip ] != AST__BAD ) {
4606 *m = 0;
4607 result = 0;
4608 } else {
4609 *m = 1;
4610 }
4611 }
4612 } else {
4613 for( ip = 0; ip < np; ip++, m++ ) {
4614 if( ptr1[ 0 ][ ip ] != AST__BAD &&
4615 ptr1[ 1 ][ ip ] != AST__BAD ) {
4616 result = 0;
4617 break;
4618 }
4619 }
4620 }
4621 }
4622
4623 /* Free resources. */
4624 tunc = astAnnul( tunc );
4625 frm = astAnnul( frm );
4626 safe = astFree( safe );
4627 pset1 = astAnnul( pset1 );
4628 pset2 = astAnnul( pset2 );
4629
4630 /* If an error has occurred, return zero. */
4631 if( !astOK ) {
4632 result = 0;
4633 if( mask ) *mask = astAnnul( *mask );
4634 }
4635
4636 /* Return the result. */
4637 return result;
4638 }
4639
RegTrace(AstRegion * this_region,int n,double * dist,double ** ptr,int * status)4640 static int RegTrace( AstRegion *this_region, int n, double *dist, double **ptr,
4641 int *status ){
4642 /*
4643 *+
4644 * Name:
4645 * RegTrace
4646
4647 * Purpose:
4648 * Return requested positions on the boundary of a 2D Region.
4649
4650 * Type:
4651 * Private function.
4652
4653 * Synopsis:
4654 * #include "polygon.h"
4655 * int astTraceRegion( AstRegion *this, int n, double *dist, double **ptr );
4656
4657 * Class Membership:
4658 * Polygon member function (overrides the astTraceRegion method
4659 * inherited from the parent Region class).
4660
4661 * Description:
4662 * This function returns positions on the boundary of the supplied
4663 * Region, if possible. The required positions are indicated by a
4664 * supplied list of scalar parameter values in the range zero to one.
4665 * Zero corresponds to some arbitrary starting point on the boundary,
4666 * and one corresponds to the end (which for a closed region will be
4667 * the same place as the start).
4668
4669 * Parameters:
4670 * this
4671 * Pointer to the Region.
4672 * n
4673 * The number of positions to return. If this is zero, the function
4674 * returns without action (but the returned function value still
4675 * indicates if the method is supported or not).
4676 * dist
4677 * Pointer to an array of "n" scalar parameter values in the range
4678 * 0 to 1.0.
4679 * ptr
4680 * A pointer to an array of pointers. The number of elements in
4681 * this array should equal tthe number of axes in the Frame spanned
4682 * by the Region. Each element of the array should be a pointer to
4683 * an array of "n" doubles, in which to return the "n" values for
4684 * the corresponding axis. The contents of the arrays are unchanged
4685 * if the supplied Region belongs to a class that does not
4686 * implement this method.
4687
4688 * Returned Value:
4689 * Non-zero if the astTraceRegion method is implemented by the class
4690 * of Region supplied, and zero if not.
4691
4692 *-
4693 */
4694
4695 /* Local Variables; */
4696 AstFrame *frm;
4697 AstMapping *map;
4698 AstPointSet *bpset;
4699 AstPointSet *cpset;
4700 AstPolygon *this;
4701 double **bptr;
4702 double d;
4703 double p[ 2 ];
4704 int i;
4705 int j0;
4706 int j;
4707 int ncur;
4708 int nv;
4709 int monotonic;
4710
4711 /* Check inherited status, and the number of points to return, returning
4712 a non-zero value to indicate that this class supports the astRegTrace
4713 method. */
4714 if( ! astOK || n == 0 ) return 1;
4715
4716 /* Get a pointer to the Polygon structure. */
4717 this = (AstPolygon *) this_region;
4718
4719 /* Ensure cached information is available. */
4720 Cache( this, status );
4721
4722 /* Get a pointer to the base Frame in the encapsulated FrameSet. */
4723 frm = astGetFrame( this_region->frameset, AST__BASE );
4724
4725 /* We first determine the required positions in the base Frame of the
4726 Region, and then transform them into the current Frame. Get the
4727 base->current Mapping, and the number of current Frame axes. */
4728 map = astGetMapping( this_region->frameset, AST__BASE, AST__CURRENT );
4729
4730 /* If it's a UnitMap we do not need to do the transformation, so put the
4731 base Frame positions directly into the supplied arrays. */
4732 if( astIsAUnitMap( map ) ) {
4733 bpset = NULL;
4734 bptr = ptr;
4735 ncur = 2;
4736
4737 /* Otherwise, create a PointSet to hold the base Frame positions (known
4738 to be 2D since this is an polygon). */
4739 } else {
4740 bpset = astPointSet( n, 2, " ", status );
4741 bptr = astGetPoints( bpset );
4742 ncur = astGetNout( map );
4743 }
4744
4745 /* Check the pointers can be used safely. */
4746 if( astOK ) {
4747
4748 /* Get the number of vertices. */
4749 nv = astGetNpoint( this_region->points );
4750
4751 /* If we have a reasonable number of pointsand there are a reasonable
4752 number of vertices, we can be quicker if we know if the parameteric
4753 distances are monotonic increasing. Find out now. */
4754 if( n > 5 && nv > 5 ) {
4755
4756 monotonic = 1;
4757 for( i = 1; i < n; i++ ) {
4758 if( dist[ i ] < dist[ i - 1 ] ) {
4759 monotonic = 0;
4760 break;
4761 }
4762 }
4763
4764 } else {
4765 monotonic = 0;
4766 }
4767
4768 /* Loop round each point. */
4769 j0 = 1;
4770 for( i = 0; i < n; i++ ) {
4771
4772 /* Get the required round the polygon, starting from vertex zero. */
4773 d = dist[ i ]*this->totlen;
4774
4775 /* Loop round each vertex until we find one which is beyond the required
4776 point. If the supplied distances are monotonic increasing, we can
4777 start the checks at the same vertex that was used for the previous
4778 since we know there will never be a backward step. */
4779 for( j = j0; j < nv; j++ ) {
4780 if( this->startsat[ j ] > d ) break;
4781 }
4782
4783 /* If the distances are monotonic increasing, record the vertex that we
4784 have reached so far. */
4785 if( monotonic ) j0 = j;
4786
4787 /* Find the distance to travel beyond the previous vertex. */
4788 d -= this->startsat[ j - 1 ];
4789
4790 /* Find the position, that is the required distance from the previous
4791 vertex towards the next vertex. */
4792 astLineOffset( frm, this->edges[ j - 1 ], d, 0.0, p );
4793
4794 /* Store the resulting axis values. */
4795 bptr[ 0 ][ i ] = p[ 0 ];
4796 bptr[ 1 ][ i ] = p[ 1 ];
4797 }
4798 }
4799
4800 /* If required, transform the base frame positions into the current
4801 Frame, storing them in the supplied array. Then free resources. */
4802 if( bpset ) {
4803 cpset = astPointSet( n, ncur, " ", status );
4804 astSetPoints( cpset, ptr );
4805
4806 (void) astTransform( map, bpset, 1, cpset );
4807
4808 cpset = astAnnul( cpset );
4809 bpset = astAnnul( bpset );
4810 }
4811
4812 /* Free remaining resources. */
4813 map = astAnnul( map );
4814 frm = astAnnul( frm );
4815
4816 /* Return a non-zero value to indicate that this class supports the
4817 astRegTrace method. */
4818 return 1;
4819 }
4820
RemoveFromChain(Segment * head,Segment * seg,int * status)4821 static Segment *RemoveFromChain( Segment *head, Segment *seg, int *status ){
4822 /*
4823 * Name:
4824 * RemoveFromChain
4825
4826 * Purpose:
4827 * Remove a Segment from the link list maintained by astDownsize.
4828
4829 * Type:
4830 * Private function.
4831
4832 * Synopsis:
4833 * #include "polygon.h"
4834 * Segment *RemoveFromChain( Segment *head, Segment *seg, int *status )
4835
4836 * Class Membership:
4837 * Polygon member function
4838
4839 * Description:
4840 * The supplied Segment is removed form the list, and the gap is
4841 * closed up.
4842
4843 * Parameters:
4844 * head
4845 * The Segment structure at the head of the list.
4846 * seg
4847 * The Segment to be removed from the list.
4848 * status
4849 * Pointer to the inherited status variable.
4850
4851 * Returned Value:
4852 * Pointer to the link head (which will have changed if "seg" was the
4853 * original head).
4854
4855 */
4856
4857 /* Check the global error status. */
4858 if ( !astOK ) return head;
4859
4860 /* If the Segment was the original head, make the next segment the new
4861 head. */
4862 if( head == seg ) head = seg->next;
4863
4864 /* Close up the links between the Segments on either side of the segment
4865 being removed. */
4866 if( seg->prev ) seg->prev->next = seg->next;
4867 if( seg->next ) seg->next->prev = seg->prev;
4868
4869 /* Nullify the links in the segment being removed. */
4870 seg->next = NULL;
4871 seg->prev = NULL;
4872
4873 /* Return the new head. */
4874 return head;
4875 }
4876
ResetCache(AstRegion * this_region,int * status)4877 static void ResetCache( AstRegion *this_region, int *status ){
4878 /*
4879 * Name:
4880 * ResetCache
4881
4882 * Purpose:
4883 * Clear cached information within the supplied Region.
4884
4885 * Type:
4886 * Private function.
4887
4888 * Synopsis:
4889 * #include "polygon.h"
4890 * void ResetCache( AstRegion *this, int *status )
4891
4892 * Class Membership:
4893 * Region member function (overrides the astResetCache method
4894 * inherited from the parent Region class).
4895
4896 * Description:
4897 * This function clears cached information from the supplied Region
4898 * structure.
4899
4900 * Parameters:
4901 * this
4902 * Pointer to the Region.
4903 * status
4904 * Pointer to the inherited status variable.
4905 */
4906
4907 /* Local Variables: */
4908 AstPolygon *this;
4909 int i;
4910 int nv;
4911
4912 /* Get a pointer to the Polygon structure. */
4913 this = (AstPolygon *) this_region;
4914
4915 /* If a Polygon was supplied, indicate cached information needs to be
4916 recalculated. */
4917 if( this ) {
4918 this->stale = 1;
4919 this->lbnd[ 0 ] = AST__BAD;
4920
4921 /* Free any edge structures (number of vertices may be about to change so
4922 this cannot be left until the next call to "Cache()". */
4923 if( this->edges ) {
4924 nv = astGetNpoint( this_region->points );
4925 for( i = 0; i < nv; i++ ) {
4926 this->edges[ i ] = astFree( this->edges[ i ] );
4927 }
4928 this->edges = astFree( this->edges );
4929 }
4930
4931 /* Clear the cache of the parent class. */
4932 (*parent_resetcache)( this_region, status );
4933 }
4934 }
4935
SetAttrib(AstObject * this_object,const char * setting,int * status)4936 static void SetAttrib( AstObject *this_object, const char *setting, int *status ) {
4937 /*
4938 * Name:
4939 * SetAttrib
4940
4941 * Purpose:
4942 * Set an attribute value for a Polygon.
4943
4944 * Type:
4945 * Private function.
4946
4947 * Synopsis:
4948 * #include "polygon.h"
4949 * void SetAttrib( AstObject *this, const char *setting, int *status )
4950
4951 * Class Membership:
4952 * Polygon member function (extends the astSetAttrib method inherited from
4953 * the Region class).
4954
4955 * Description:
4956 * This function assigns an attribute value for a Polygon, the attribute
4957 * and its value being specified by means of a string of the form:
4958 *
4959 * "attribute= value "
4960 *
4961 * Here, "attribute" specifies the attribute name and should be in lower
4962 * case with no white space present. The value to the right of the "="
4963 * should be a suitable textual representation of the value to be assigned
4964 * and this will be interpreted according to the attribute's data type.
4965 * White space surrounding the value is only significant for string
4966 * attributes.
4967
4968 * Parameters:
4969 * this
4970 * Pointer to the Polygon.
4971 * setting
4972 * Pointer to a null terminated string specifying the new attribute
4973 * value.
4974 * status
4975 * Pointer to the inherited status variable.
4976
4977 * Returned Value:
4978 * void
4979 */
4980
4981 /* Local Vaiables: */
4982 AstPolygon *this; /* Pointer to the Polygon structure */
4983 int ival; /* Integer attribute value */
4984 int len; /* Length of setting string */
4985 int nc; /* Number of characters read by astSscanf */
4986
4987 /* Check the global error status. */
4988 if ( !astOK ) return;
4989
4990 /* Obtain a pointer to the Polygon structure. */
4991 this = (AstPolygon *) this_object;
4992
4993 /* Obtain the length of the setting string. */
4994 len = strlen( setting );
4995
4996 /* Test for each recognised attribute in turn, using "astSscanf" to parse the
4997 setting string and extract the attribute value (or an offset to it in the
4998 case of string values). In each case, use the value set in "nc" to check
4999 that the entire string was matched. Once a value has been obtained, use the
5000 appropriate method to set it. */
5001
5002 /* SimpVertices. */
5003 /* ------------- */
5004 if ( nc = 0,
5005 ( 1 == astSscanf( setting, "simpvertices= %d %n", &ival, &nc ) )
5006 && ( nc >= len ) ) {
5007 astSetSimpVertices( this, ival );
5008
5009 /* Pass any unrecognised setting to the parent method for further
5010 interpretation. */
5011 } else {
5012 (*parent_setattrib)( this_object, setting, status );
5013 }
5014 }
5015
SetPointSet(AstPolygon * this,AstPointSet * pset,int * status)5016 static void SetPointSet( AstPolygon *this, AstPointSet *pset, int *status ){
5017 /*
5018 * Name:
5019 * SetPointSet
5020
5021 * Purpose:
5022 * Store a new set of vertices in an existing Polygon.
5023
5024 * Type:
5025 * Private function.
5026
5027 * Synopsis:
5028 * #include "polygon.h"
5029 * void SetPointSet( AstPolygon *this, AstPointSet *pset, int *status )
5030
5031 * Class Membership:
5032 * Polygon member function
5033
5034 * Description:
5035 * The PointSet in the supplied Polygon is annulled, and replaced by a
5036 * clone of the supplied PointSet pointer.
5037
5038 * Parameters:
5039 * this
5040 * Pointer to the Polygon to be changed.
5041 * pset
5042 * The PointSet containing the new vertex information.
5043 * status
5044 * Pointer to the inherited status variable.
5045
5046 */
5047
5048
5049 /* Check the global error status. */
5050 if ( !astOK ) return;
5051
5052 /* Indicate the cached information in the polygon will need to be
5053 re-calculated when needed. */
5054 astResetCache( this );
5055
5056 /* Annul the pointer to the PointSet already in the supplied Polygon. */
5057 (void) astAnnul( ((AstRegion *) this)->points );
5058
5059 /* Store a clone of the supplied new PointSet pointer. */
5060 ((AstRegion *) this)->points = astClone( pset );
5061
5062 }
5063
SetRegFS(AstRegion * this_region,AstFrame * frm,int * status)5064 static void SetRegFS( AstRegion *this_region, AstFrame *frm, int *status ) {
5065 /*
5066 * Name:
5067 * SetRegFS
5068
5069 * Purpose:
5070 * Stores a new FrameSet in a Region
5071
5072 * Type:
5073 * Private function.
5074
5075 * Synopsis:
5076 * #include "polygon.h"
5077 * void SetRegFS( AstRegion *this_region, AstFrame *frm, int *status )
5078
5079 * Class Membership:
5080 * Polygon method (over-rides the astSetRegFS method inherited from
5081 * the Region class).
5082
5083 * Description:
5084 * This function creates a new FrameSet and stores it in the supplied
5085 * Region. The new FrameSet contains two copies of the supplied
5086 * Frame, connected by a UnitMap.
5087
5088 * Parameters:
5089 * this
5090 * Pointer to the Region.
5091 * frm
5092 * The Frame to use.
5093 * status
5094 * Pointer to the inherited status variable.
5095
5096 */
5097
5098
5099 /* Check the global error status. */
5100 if ( !astOK ) return;
5101
5102 /* Indicate cached information eeds re-calculating. */
5103 astResetCache( this_region );
5104
5105 /* Invoke the parent method to store the FrameSet in the parent Region
5106 structure. */
5107 (* parent_setregfs)( this_region, frm, status );
5108
5109 }
5110
Simplify(AstMapping * this_mapping,int * status)5111 static AstMapping *Simplify( AstMapping *this_mapping, int *status ) {
5112 /*
5113 * Name:
5114 * Simplify
5115
5116 * Purpose:
5117 * Simplify the Mapping represented by a Region.
5118
5119 * Type:
5120 * Private function.
5121
5122 * Synopsis:
5123 * #include "polygon.h"
5124 * AstMapping *Simplify( AstMapping *this, int *status )
5125
5126 * Class Membership:
5127 * Polygon method (over-rides the astSimplify method inherited
5128 * from the Region class).
5129
5130 * Description:
5131 * This function invokes the parent Region Simplify method, and then
5132 * performs any further region-specific simplification.
5133 *
5134 * If the Mapping from base to current Frame is not a UnitMap, this
5135 * will include attempting to fit a new Region to the boundary defined
5136 * in the current Frame.
5137
5138 * Parameters:
5139 * this
5140 * Pointer to the original Region.
5141 * status
5142 * Pointer to the inherited status variable.
5143
5144 * Returned Value:
5145 * A pointer to the simplified Region. A cloned pointer to the
5146 * supplied Region will be returned if no simplication could be
5147 * performed.
5148
5149 * Notes:
5150 * - A NULL pointer value will be returned if this function is
5151 * invoked with the AST error status set, or if it should fail for
5152 * any reason.
5153 */
5154
5155 /* Local Variables: */
5156 AstFrame *frm; /* Current Frame */
5157 AstMapping *map; /* Base -> current Mapping */
5158 AstMapping *result; /* Result pointer to return */
5159 AstPointSet *mesh; /* Mesh of current Frame positions */
5160 AstPointSet *ps2; /* Polygon PointSet in current Frame */
5161 AstPolygon *newpol; /* New Polygon */
5162 AstRegion *new; /* Pointer to simplified Region */
5163 AstRegion *this; /* Pointer to supplied Region structure */
5164 AstRegion *unc; /* Pointer to uncertainty Region */
5165 double **ptr2; /* Pointer axis values in "ps2" */
5166 double *mem; /* Pointer to array holding new vertex coords */
5167 double *p; /* Pointer to next vertex coords */
5168 double *q; /* Pointer to next vertex coords */
5169 int iv; /* Vertex index */
5170 int nv; /* Number of vertices in polygon */
5171 int ok; /* Are the new polygon vertices good? */
5172 int simpler; /* Has some simplication taken place? */
5173
5174 /* Initialise. */
5175 result = NULL;
5176
5177 /* Check the global error status. */
5178 if ( !astOK ) return result;
5179
5180 /* Get a pointer to the supplied Region structure. */
5181 this = (AstRegion *) this_mapping;
5182
5183 /* Invoke the parent Simplify method inherited from the Region class. This
5184 will simplify the encapsulated FrameSet and uncertainty Region. */
5185 new = (AstRegion *) (*parent_simplify)( this_mapping, status );
5186
5187 /* Note if any simplification took place. This is assumed to be the case
5188 if the pointer returned by the above call is different to the supplied
5189 pointer. */
5190 simpler = ( new != this );
5191
5192 /* We attempt to simplify the Polygon by re-defining it within its current
5193 Frame. Transforming the Polygon from its base to its current Frame may
5194 result in the region no having the same edges. If required, we test this
5195 by transforming a set of bounds on the Polygon boundary. This can only
5196 be done if the current Frame is 2-dimensional. Also, there is only any
5197 point in doing it if the Mapping from base to current Frame in the
5198 Polygon is not a UnitMap. */
5199 map = astGetMapping( new->frameset, AST__BASE, AST__CURRENT );
5200 if( !astIsAUnitMap( map ) && astGetNout( map ) == 2 ) {
5201
5202 /* Get a pointer to the Frame represented by the Polgon. */
5203 frm = astGetFrame( new->frameset, AST__CURRENT );
5204
5205 /* Get the Region describing the positional uncertainty in this Frame. */
5206 unc = astGetUncFrm( new, AST__CURRENT );
5207
5208 /* Get the positions of the vertices within this Frame. */
5209 ps2 = astRegTransform( this, this->points, 1, NULL, NULL );
5210 ptr2 = astGetPoints( ps2 );
5211
5212 /* Get the number of vertices. */
5213 nv = astGetNpoint( ps2 );
5214
5215 /* Re-organise the vertex axis values into the form required by the
5216 Polygon constructor function. */
5217 mem = astMalloc( sizeof( double)*(size_t)( 2*nv ) );
5218 if( astOK ) {
5219 ok = 1;
5220 p = mem;
5221 q = ptr2[ 0 ];
5222 for( iv = 0; iv < nv; iv++ ) {
5223 if( ( *(p++) = *(q++) ) == AST__BAD ) ok = 0;
5224 }
5225 q = ptr2[ 1 ];
5226 for( iv = 0; iv < nv; iv++ ) *(p++) = *(q++);
5227
5228 /* Create a new Polygon using these transformed vertices. */
5229 if( ok ) {
5230 newpol = astPolygon( frm, nv, nv, mem, unc, "", status );
5231
5232 /* If the SimpVertices attribute is zero, we now check that the
5233 transformation has not bent the edges of the polygon significantly.
5234 If it has, we annul the new Polygon. */
5235 if( !astGetSimpVertices( this ) ) {
5236
5237 /* Get a mesh of points covering the Polygon in this Frame. */
5238 mesh = astRegMesh( new );
5239
5240 /* See if all points within the mesh created from the original Polygon fall
5241 on the boundary of the new Polygon, to within the uncertainty of the
5242 Region. If not, annul the new Polgon. */
5243 if( !astRegPins( newpol, mesh, NULL, NULL ) ) {
5244 newpol = astAnnul( newpol );
5245 }
5246
5247 /* Free the mesh. */
5248 mesh = astAnnul( mesh );
5249 }
5250
5251 /* If we still have a new polygon, use the new Polygon in place of the
5252 original Region. */
5253 if( newpol ) {
5254 (void) astAnnul( new );
5255 new = (AstRegion *) newpol;
5256 simpler = 1;
5257 }
5258 }
5259 }
5260
5261 /* Free other resources. */
5262 frm = astAnnul( frm );
5263 unc = astAnnul( unc );
5264 ps2 = astAnnul( ps2 );
5265 mem = astFree( mem );
5266 }
5267 map = astAnnul( map );
5268
5269 /* If any simplification could be performed, copy Region attributes from
5270 the supplied Region to the returned Region, and return a pointer to it.
5271 If the supplied Region had no uncertainty, ensure the returned Region
5272 has no uncertainty. Otherwise, return a clone of the supplied pointer. */
5273 if( simpler ){
5274 astRegOverlay( new, this, 1 );
5275 result = (AstMapping *) new;
5276
5277 } else {
5278 new = astAnnul( new );
5279 result = astClone( this );
5280 }
5281
5282 /* If an error occurred, annul the returned pointer. */
5283 if ( !astOK ) result = astAnnul( result );
5284
5285 /* Return the result. */
5286 return result;
5287 }
5288
SmoothPoly(AstPointSet * pset,int boxsize,double strength,int * status)5289 static void SmoothPoly( AstPointSet *pset, int boxsize, double strength,
5290 int *status ) {
5291 /*
5292 * Name:
5293 * SmoothPoly
5294
5295 * Purpose:
5296 * Smoooth a polygon assuming plane geometry.
5297
5298 * Type:
5299 * Private function.
5300
5301 * Synopsis:
5302 * #include "polygon.h"
5303 * void SmoothPoly( AstPointSet *pset, int boxsize, double strength,
5304 * int *status )
5305
5306 * Class Membership:
5307 * Polygon member function
5308
5309 * Description:
5310 * This function smooths a polygon, without changing the number of
5311 * vertices. It assumes plane geometry, so should not be used to
5312 * smooth polygons defined within a SkyFrame or CmpFrame.
5313 *
5314 * Each vertex is replaced by a new vertex determined as follows: the
5315 * mean X and Y axis value of the vertices in a section of the polygon
5316 * centred on the vertex being replaced are found (the length of the
5317 * section is given by parameter "boxsize"). The new vertex position
5318 * is then the weighted mean of this mean (X,Y) position and the old
5319 * vertex position. The weight for the mean (X,Y) position is given
5320 * by parameter "strength", and the weight for the old vertex
5321 * position is (1.0 - strength)
5322
5323 * Parameters:
5324 * pset
5325 * A PointSet holding the polygon vertices.
5326 * boxsize
5327 * Half width of the box filter, given as a number of vertices.
5328 * strength
5329 * The weight to use for the mean (X,Y) position when finding each
5330 * new vertex position. Should be in the range 0.0 to 1.0. A value
5331 * of zero results in no change being made to the polygon. A value
5332 * of 1.0 results in the returned polygon being fully smoothed.
5333 * status
5334 * Pointer to the inherited status variable.
5335
5336 */
5337
5338 /* Local Variables: */
5339 double **ptr;
5340 double *newx;
5341 double *newy;
5342 double *nx;
5343 double *ny;
5344 double *oldx;
5345 double *oldy;
5346 double *ox;
5347 double *oy;
5348 double *px;
5349 double *py;
5350 double *qx;
5351 double *qy;
5352 double a;
5353 double b;
5354 double sx;
5355 double sy;
5356 int half_width;
5357 int i;
5358 int nv;
5359 int top;
5360
5361 /* Check the global error status. */
5362 if ( !astOK ) return;
5363
5364 /* Get the number of vertices. */
5365 nv = astGetNpoint( pset );
5366
5367 /* Get pointers to arrays holding the supplied vertex positions. */
5368 ptr = astGetPoints( pset );
5369 oldx = ptr[ 0 ];
5370 oldy = ptr[ 1 ];
5371
5372 /* Allocate arrays to hold the returned vertex positions. */
5373 newx = astMalloc( nv*sizeof( double ) );
5374 newy = astMalloc( nv*sizeof( double ) );
5375
5376 /* Check these pointers can be used safely. */
5377 if( astOK ) {
5378
5379 /* Get weighting factors for the fully smoothed and unsmoothed positions. */
5380 a = strength;
5381 b = 1.0 - a;
5382
5383 /* Ensure the box size is sufficiently small for there to be room for
5384 two boxes along the polygon. */
5385 half_width = nv/4 - 1;
5386 if( boxsize < half_width ) half_width = boxsize;
5387 if( half_width < 1 ) half_width = 1;
5388
5389 /* Modify the weight for the fully smoothed position to include the
5390 normalisation factor needed to account for the box width. */
5391 a /= 2*half_width + 1;
5392
5393 /* Find the sum of the x and y coordinates within a box centred on the
5394 first vertex. This includes vertices from the end of the polygon. */
5395 px = oldx + 1;
5396 qx = oldx + nv;
5397 sx = (oldx)[ 0 ];
5398
5399 py = oldy + 1;
5400 qy = oldy + nv;
5401 sy = (oldy)[ 0 ];
5402
5403 for( i = 0; i < half_width; i++ ) {
5404 sx += *(px++) + *(--qx);
5405 sy += *(py++) + *(--qy);
5406 }
5407
5408 /* Replacing vertices within the first half box will include vertices at
5409 both ends of the polygon. Set up the pointers accordingly, and then
5410 find replacements for each vertex in the first half box.*/
5411 ox = oldx;
5412 oy = oldy;
5413 nx = newx;
5414 ny = newy;
5415 for( i = 0; i < half_width; i++ ) {
5416
5417 /* Form the new vertex (x,y) values as the weighted mean of the mean
5418 (x,y) values in the box, and the old (x,y) values. */
5419 *(nx++) = a*sx + b*( *(ox++) );
5420 *(ny++) = a*sy + b*( *(oy++) );
5421
5422 /* Add in the next vertex X and Y axis values to the running sums, and
5423 remove the position that has now passed out of the box. */
5424 sx += *(px++) - *(qx++);
5425 sy += *(py++) - *(qy++);
5426 }
5427
5428 /* Adjust the pointer for the rest of the polygon, up to one half box away
5429 from the end. In this section, the smoothing box does not touch either
5430 end of the polygon. */
5431 top = nv - half_width - 1;
5432 qx = oldx;
5433 qy = oldy;
5434 for( ; i < top; i++ ){
5435
5436 /* Form the new vertex (x,y) values as the weighted mean of the mean
5437 (x,y) values in the box, and the old (x,y) values. */
5438 *(nx++) = a*sx + b*( *(ox++) );
5439 *(ny++) = a*sy + b*( *(oy++) );
5440
5441 /* Add in the next vertex X and Y axis values to the running sums, and
5442 remove the position that has now passed out of the box. */
5443 sx += *(px++) - *(qx++);
5444 sy += *(py++) - *(qy++);
5445 }
5446
5447 /* Now do the last half box (which includes vertices from the start of
5448 the polygon). */
5449 top = nv;
5450 px = oldx;
5451 py = oldy;
5452 for( ; i < top; i++ ){
5453 *(nx++) = a*sx + b*( *(ox++) );
5454 *(ny++) = a*sy + b*( *(oy++) );
5455 sx += *(px++) - *(qx++);
5456 sy += *(py++) - *(qy++);
5457 }
5458
5459 /* Replace the data points in the PointSet. */
5460 ptr[ 0 ] = newx;
5461 ptr[ 1 ] = newy;
5462 oldx = astFree( oldx );
5463 oldy = astFree( oldy );
5464 }
5465 }
5466
TestAttrib(AstObject * this_object,const char * attrib,int * status)5467 static int TestAttrib( AstObject *this_object, const char *attrib, int *status ) {
5468 /*
5469 * Name:
5470 * TestAttrib
5471
5472 * Purpose:
5473 * Test if a specified attribute value is set for a Polygon.
5474
5475 * Type:
5476 * Private function.
5477
5478 * Synopsis:
5479 * #include "polygon.h"
5480 * int TestAttrib( AstObject *this, const char *attrib, int *status )
5481
5482 * Class Membership:
5483 * Polygon member function (over-rides the astTestAttrib protected
5484 * method inherited from the Region class).
5485
5486 * Description:
5487 * This function returns a boolean result (0 or 1) to indicate whether
5488 * a value has been set for one of a Polygon's attributes.
5489
5490 * Parameters:
5491 * this
5492 * Pointer to the Polygon.
5493 * attrib
5494 * Pointer to a null terminated string specifying the attribute
5495 * name. This should be in lower case with no surrounding white
5496 * space.
5497 * status
5498 * Pointer to the inherited status variable.
5499
5500 * Returned Value:
5501 * One if a value has been set, otherwise zero.
5502
5503 * Notes:
5504 * - A value of zero will be returned if this function is invoked
5505 * with the global status set, or if it should fail for any reason.
5506 */
5507
5508 /* Local Variables: */
5509 AstPolygon *this; /* Pointer to the Polygon structure */
5510 int result; /* Result value to return */
5511
5512 /* Initialise. */
5513 result = 0;
5514
5515 /* Check the global error status. */
5516 if ( !astOK ) return result;
5517
5518 /* Obtain a pointer to the Polygon structure. */
5519 this = (AstPolygon *) this_object;
5520
5521 /* Check the attribute name and test the appropriate attribute. */
5522
5523 /* SimpVertices. */
5524 /* ------------- */
5525 if ( !strcmp( attrib, "simpvertices" ) ) {
5526 result = astTestSimpVertices( this );
5527
5528 /* If the attribute is not recognised, pass it on to the parent method
5529 for further interpretation. */
5530 } else {
5531 result = (*parent_testattrib)( this_object, attrib, status );
5532 }
5533
5534 /* Return the result, */
5535 return result;
5536 }
5537
5538 /*
5539 * Name:
5540 * TraceEdge
5541
5542 * Purpose:
5543 * Find a point that is inside the required outline.
5544
5545 * Type:
5546 * Private function.
5547
5548 * Synopsis:
5549 * #include "polygon.h"
5550 * void TraceEdge<Oper><X>( <Xtype> value, const <Xtype> array[],
5551 * const int lbnd[ 2 ], const int ubnd[ 2 ],
5552 * int iv0, int ix0, int iy0, int starpix,
5553 * int full, int *status );
5554
5555 * Class Membership:
5556 * Polygon member function
5557
5558 * Description:
5559 * This function forms a polygon enclosing the region of the data
5560 * array specified by <Oper> and "value". If this polygon contains
5561 * the point "(inx,iny)", then a PointSet is returned holding the
5562 * pixel coordinates of the Polygon vertices. If the polygon
5563 * does not contain "(inx,iny)", a NULL pointer is returned.
5564 *
5565 * Each vertex in the polygon corresponds to a corner of a pixel in
5566 * the data array.
5567
5568 * Parameters:
5569 * value
5570 * The data value defining valid pixels.
5571 * array
5572 * The data array.
5573 * lbnd
5574 * The lower pixel index bounds of the array.
5575 * ubnd
5576 * The upper pixel index bounds of the array.
5577 * iv0
5578 * The vector index of a pixel inside the region such that the
5579 * pixel to the right is NOT inside the region. This defines the
5580 * start of the polygon.
5581 * ix0
5582 * The X pixel index of the pixel specified by "iv0".
5583 * inx
5584 * The X pixel index of a point which must be inside the polygon
5585 * for the polygon to be acceptable.
5586 * iny
5587 * The Y pixel index of a point which must be inside the polygon
5588 * for the polygon to be acceptable.
5589 * starpix
5590 * If non-zero, the usual Starlink definition of pixel coordinate
5591 * is used (integral values at pixel corners). Otherwise, the
5592 * system used by other AST functions such as astResample is used
5593 * (integral values at pixel centres).
5594 * full
5595 * If non-zero, the full polygon is stored. If zero, vertices in the
5596 * middle of straight sections of the Polygon are omitted.
5597 * status
5598 * Pointer to the inherited status variable.
5599
5600 * Returned Value:
5601 * A pointer to a PointSet holding the vertices of the polygon, or
5602 * NULL if the polygon did not contain "(inx,iny)".
5603
5604 * Notes:
5605 * - <Oper> must be one of LT, LE, EQ, GE, GT, NE.
5606
5607
5608 */
5609
5610 /* Define a macro to implement the function for a specific data
5611 type and operation. */
5612 #define MAKE_TRACEEDGE(X,Xtype,Oper,OperI) \
5613 static AstPointSet *TraceEdge##Oper##X( Xtype value, const Xtype array[], \
5614 const int lbnd[ 2 ], const int ubnd[ 2 ], \
5615 int iv0, int ix0, int iy0, \
5616 int starpix, int full, \
5617 int *status ){ \
5618 \
5619 /* Local Variables: */ \
5620 AstPointSet *result; /* Pointer to text describing oper */ \
5621 const Xtype *pa; /* Pointer to current valid pixel value */ \
5622 const Xtype *pb; /* Pointer to neigbouring valid pixel value */ \
5623 const Xtype *pc; /* Pointer to neigbouring valid pixel value */ \
5624 double *ptr[ 2 ]; /* PointSet data pointers */ \
5625 double *xvert; /* Pointer to array holding vertex X axis values */ \
5626 double *yvert; /* Pointer to array holding vertex Y axis values */ \
5627 double xx; /* Pixel X coord at corner */ \
5628 double yy; /* Pixel Y coord at corner */ \
5629 int at; /* The pixel corner to draw to */ \
5630 int done; /* Have we arrived back at the start of the polygon? */ \
5631 int ii; /* Index of new vertex */ \
5632 int ix; /* X pixel index of current valid pixel */ \
5633 int iy; /* Y pixel index of current valid pixel */ \
5634 int nright; /* Overall number of right hand turns along polygon */ \
5635 int nvert; /* Number of vertices */ \
5636 int nx; /* Pixels per row */ \
5637 \
5638 /* Initialise */ \
5639 result = NULL; \
5640 \
5641 /* Check the global error status. */ \
5642 if ( !astOK ) return result; \
5643 \
5644 \
5645 /* Initialise pointers to arrays holding the X and Y pixel coordinates at \
5646 the vertices of the polygon. */ \
5647 xvert = NULL; \
5648 yvert = NULL; \
5649 nvert = 0; \
5650 \
5651 /* Find number of pixels in one row of the array. */ \
5652 nx = ( ubnd[ 0 ] - lbnd[ 0 ] + 1 ); \
5653 \
5654 /* The four corners of a pixel are numbered as follows: 0=bottom left, \
5655 1=top left, 2=top right, 3=bottom right. The following algorithm moves \
5656 along pixel edges, from corner to corner, using the above numbering \
5657 scheme to identify the corners. We start the polygon by moving from the \
5658 bottom right to the top right corner of pixel "(ix0,iy0)". */ \
5659 ix = ix0; \
5660 iy = iy0; \
5661 at = 2; \
5662 \
5663 /* Store a pointer to the first good pixel value. */ \
5664 pa = array + ( ix - lbnd[ 0 ] ) + nx*( iy - lbnd[ 1 ] ) ; \
5665 \
5666 /* We count the number of times the polygon turns to the right compared \
5667 to the left. Initialise it to zero. */ \
5668 nright = 0; \
5669 \
5670 /* Loop round tracing out the polygon until we arrive back at the \
5671 beginning. The Polygon class requires that the inside of the polygon \
5672 is to the left as we move round the polygon in an anti-clockwise \
5673 direction. So at each corner, we attempt to move to the next \
5674 anti-clockwise good pixel corner. */ \
5675 done = 0; \
5676 while( !done ) { \
5677 \
5678 /* If we have arrived at the bottom left corner of the good pixel, we must \
5679 have come from the top left corner since all movements around a pixel \
5680 must be anti-clockwise. */ \
5681 if( at == 0 ) { \
5682 \
5683 /* Note the pixel coordinates at the bottom left corner of the current \
5684 pixel. */ \
5685 if( starpix ) { \
5686 xx = ix - 1.0; \
5687 yy = iy - 1.0; \
5688 } else { \
5689 xx = ix - 0.5; \
5690 yy = iy - 0.5; \
5691 } \
5692 \
5693 /* Get a pointer to lower left pixel value */ \
5694 pb = pa - nx - 1; \
5695 \
5696 /* Get a pointer to lower mid pixel value. */ \
5697 pc = pb + 1; \
5698 \
5699 /* If the lower left pixel is within the array and meets the validity \
5700 requirements, move to the left along its top edge. */ \
5701 if( iy > lbnd[ 1 ] && ix > lbnd[ 0 ] && ISVALID(*pb,OperI,value) ) { \
5702 nright++; \
5703 pa = pb; \
5704 at = 1; \
5705 ix--; \
5706 iy--; \
5707 \
5708 /* Otherwise, if lower mid pixel is good, move down its left edge. */ \
5709 } else if( iy > lbnd[ 1 ] && ISVALID(*pc,OperI,value) ) { \
5710 pa = pc; \
5711 at = 0; \
5712 iy--; \
5713 \
5714 /* Otherwise, move to the right along the bottom edge of the current pixel. */ \
5715 } else { \
5716 nright--; \
5717 at = 3; \
5718 } \
5719 \
5720 /* If the polygon bends at this point, or if we will be smoothing the \
5721 polygon, append the pixel coordinates at this pixel corner to the \
5722 polygon. */ \
5723 if( full || pa != pc ) ADD( xx, yy ); \
5724 \
5725 /* If we have arrived at the top left corner of the good pixel, we must \
5726 have come from the top right corner. */ \
5727 } else if( at == 1 ) { \
5728 \
5729 /* Note the pixel coordinates at the top left corner of the current \
5730 pixel. */ \
5731 if( starpix ) { \
5732 xx = ix - 1.0; \
5733 yy = iy; \
5734 } else { \
5735 xx = ix - 0.5; \
5736 yy = iy + 0.5; \
5737 } \
5738 \
5739 /* Get a pointer to upper left pixel value */ \
5740 pb = pa + nx - 1; \
5741 \
5742 /* Get a pointer to mid left pixel value. */ \
5743 pc = pa - 1; \
5744 \
5745 /* If upper left pixel is good, move up its left edge. */ \
5746 if( iy < ubnd[ 1 ] && ix > lbnd[ 0 ] && ISVALID(*pb,OperI,value) ) { \
5747 nright++; \
5748 pa = pb; \
5749 at = 2; \
5750 ix--; \
5751 iy++; \
5752 \
5753 /* Otherwise, if left mid pixel is good, move left along its top edge. */ \
5754 } else if( ix > lbnd[ 0 ] && ISVALID(*pc,OperI,value) ) { \
5755 pa = pc; \
5756 at = 1; \
5757 ix--; \
5758 \
5759 /* Otherwise, move down the left edge of the current pixel. */ \
5760 } else { \
5761 nright--; \
5762 at = 0; \
5763 } \
5764 \
5765 /* If the polygon bends at this point, or if we will be smoothing the \
5766 polygon, append the pixel coordinates at this pixel corner to the \
5767 polygon. */ \
5768 if( full || pa != pc ) ADD( xx, yy ); \
5769 \
5770 /* If we have arrived at the top right corner of the good pixel, we must \
5771 have come from the bottom right corner. */ \
5772 } else if( at == 2 ) { \
5773 \
5774 /* Note the pixel coordinates at the top right corner of the current \
5775 pixel. */ \
5776 if( starpix ) { \
5777 xx = ix; \
5778 yy = iy; \
5779 } else { \
5780 xx = ix + 0.5; \
5781 yy = iy + 0.5; \
5782 } \
5783 \
5784 /* Pointer to upper right pixel value */ \
5785 pb = pa + nx + 1; \
5786 \
5787 /* Pointer to top mid pixel value. */ \
5788 pc = pa + nx; \
5789 \
5790 /* If upper right pixel is good, move right along its bottom edge. */ \
5791 if( iy < ubnd[ 1 ] && ix < ubnd[ 0 ] && ISVALID(*pb,OperI,value) ){ \
5792 nright++; \
5793 pa = pb; \
5794 at = 3; \
5795 ix++; \
5796 iy++; \
5797 \
5798 /* Otherwise, if upper mid pixel is good, move up its right edge. */ \
5799 } else if( iy < ubnd[ 1 ] && ISVALID(*pc,OperI,value) ) { \
5800 pa = pc; \
5801 at = 2; \
5802 iy++; \
5803 \
5804 /* Otherwise, move left along the top edge of the current pixel. */ \
5805 } else { \
5806 nright--; \
5807 at = 1; \
5808 } \
5809 \
5810 /* If the polygon bends at this point, or if we will be smoothing the \
5811 polygon, append the pixel coordinates at this pixel corner to the \
5812 polygon. */ \
5813 if( full || pa != pc ) ADD( xx, yy ); \
5814 \
5815 /* Arrived at bottom right corner of good pixel from lower left. */ \
5816 } else { \
5817 \
5818 /* Note the pixel coordinates at the bottom right corner of the current \
5819 pixel. */ \
5820 if( starpix ) { \
5821 xx = ix; \
5822 yy = iy - 1.0; \
5823 } else { \
5824 xx = ix + 0.5; \
5825 yy = iy - 0.5; \
5826 } \
5827 \
5828 /* Pointer to lower right pixel value */ \
5829 pb = pa - ( nx - 1 ); \
5830 \
5831 /* Pointer to mid right pixel value. */ \
5832 pc = pa + 1; \
5833 \
5834 /* If lower right pixel is good, move down its left edge. */ \
5835 if( iy > lbnd[ 1 ] && ix < ubnd[ 0 ] && ISVALID(*pb,OperI,value) ) { \
5836 nright++; \
5837 pa = pb; \
5838 at = 0; \
5839 ix++; \
5840 iy--; \
5841 \
5842 /* Otherwise, if right mid pixel is good, move right along its lower edge. */ \
5843 } else if( ix < ubnd[ 0 ] && ISVALID(*pc,OperI,value) ) { \
5844 pa = pc; \
5845 at = 3; \
5846 ix++; \
5847 \
5848 /* Otherwise, move up the left edge of the current pixel. */ \
5849 } else { \
5850 nright--; \
5851 at = 2; \
5852 } \
5853 \
5854 /* If the polygon bends at this point, or if we will be smoothing the \
5855 polygon, append the pixel coordinates at this pixel corner to the \
5856 polygon. */ \
5857 if( full || pa != pc ) ADD( xx, yy ); \
5858 } \
5859 \
5860 /* If we have arrived back at the start, break out of the loop. */ \
5861 if( ix == ix0 && iy == iy0 && at == 2 ) done = 1; \
5862 } \
5863 \
5864 /* If we have circled round to the right, the polygon will not enclosed \
5865 the specified position, so free resources and return a NULL pointer. */ \
5866 if( nright > 0 ) { \
5867 xvert = astFree( xvert ); \
5868 yvert = astFree( yvert ); \
5869 \
5870 /* If we have circled round to the left, the polygon will enclose \
5871 the specified position, so create and return a PointSet. */ \
5872 } else { \
5873 result = astPointSet( nvert, 2, " ", status ); \
5874 ptr[ 0 ] = xvert; \
5875 ptr[ 1 ] = yvert; \
5876 astSetPoints( result, ptr ); \
5877 } \
5878 \
5879 /* Annul the returned PointSet if anythign went wrong. */ \
5880 if( !astOK && result ) result = astAnnul( result ); \
5881 \
5882 /* Return the PointSet pointer. */ \
5883 return result; \
5884 }
5885
5886 /* Define a macro to add a vertex position to dynamically allocated
5887 arrays of X and Y positions. */
5888 #define ADD(X,Y) {\
5889 ii = nvert++; \
5890 xvert = (double *) astGrow( xvert, nvert, sizeof( double ) ); \
5891 yvert = (double *) astGrow( yvert, nvert, sizeof( double ) ); \
5892 if( astOK ) { \
5893 xvert[ ii ] = (X); \
5894 yvert[ ii ] = (Y); \
5895 } \
5896 }
5897
5898 /* Define a macro that uses the above macro to to create implementations
5899 of TraceEdge for all operations. */
5900 #define MAKEALL_TRACEEDGE(X,Xtype) \
5901 MAKE_TRACEEDGE(X,Xtype,LT,AST__LT) \
5902 MAKE_TRACEEDGE(X,Xtype,LE,AST__LE) \
5903 MAKE_TRACEEDGE(X,Xtype,EQ,AST__EQ) \
5904 MAKE_TRACEEDGE(X,Xtype,NE,AST__NE) \
5905 MAKE_TRACEEDGE(X,Xtype,GE,AST__GE) \
5906 MAKE_TRACEEDGE(X,Xtype,GT,AST__GT)
5907
5908 /* Expand the above macro to generate a function for each required
5909 data type and operation. */
5910 #if HAVE_LONG_DOUBLE /* Not normally implemented */
MAKEALL_TRACEEDGE(LD,long double)5911 MAKEALL_TRACEEDGE(LD,long double)
5912 #endif
5913 MAKEALL_TRACEEDGE(D,double)
5914 MAKEALL_TRACEEDGE(L,long int)
5915 MAKEALL_TRACEEDGE(UL,unsigned long int)
5916 MAKEALL_TRACEEDGE(I,int)
5917 MAKEALL_TRACEEDGE(UI,unsigned int)
5918 MAKEALL_TRACEEDGE(S,short int)
5919 MAKEALL_TRACEEDGE(US,unsigned short int)
5920 MAKEALL_TRACEEDGE(B,signed char)
5921 MAKEALL_TRACEEDGE(UB,unsigned char)
5922 MAKEALL_TRACEEDGE(F,float)
5923
5924 /* Undefine the macros. */
5925 #undef MAKE_TRACEEDGE
5926 #undef MAKEALL_TRACEEDGE
5927
5928 static AstPointSet *Transform( AstMapping *this_mapping, AstPointSet *in,
5929 int forward, AstPointSet *out, int *status ) {
5930 /*
5931 * Name:
5932 * Transform
5933
5934 * Purpose:
5935 * Apply a Polygon to transform a set of points.
5936
5937 * Type:
5938 * Private function.
5939
5940 * Synopsis:
5941 * #include "polygon.h"
5942 * AstPointSet *Transform( AstMapping *this, AstPointSet *in,
5943 * int forward, AstPointSet *out, int *status )
5944
5945 * Class Membership:
5946 * Polygon member function (over-rides the astTransform protected
5947 * method inherited from the Region class).
5948
5949 * Description:
5950 * This function takes a Polygon and a set of points encapsulated in a
5951 * PointSet and transforms the points by setting axis values to
5952 * AST__BAD for all points which are outside the region. Points inside
5953 * the region are copied unchanged from input to output.
5954
5955 * Parameters:
5956 * this
5957 * Pointer to the Polygon.
5958 * in
5959 * Pointer to the PointSet holding the input coordinate data.
5960 * forward
5961 * A non-zero value indicates that the forward coordinate transformation
5962 * should be applied, while a zero value requests the inverse
5963 * transformation.
5964 * out
5965 * Pointer to a PointSet which will hold the transformed (output)
5966 * coordinate values. A NULL value may also be given, in which case a
5967 * new PointSet will be created by this function.
5968 * status
5969 * Pointer to the inherited status variable.
5970
5971 * Returned Value:
5972 * Pointer to the output (possibly new) PointSet.
5973
5974 * Notes:
5975 * - The forward and inverse transformations are identical for a
5976 * Region.
5977 * - A null pointer will be returned if this function is invoked with the
5978 * global error status set, or if it should fail for any reason.
5979 * - The number of coordinate values per point in the input PointSet must
5980 * match the number of axes in the Frame represented by the Polygon.
5981 * - If an output PointSet is supplied, it must have space for sufficient
5982 * number of points and coordinate values per point to accommodate the
5983 * result. Any excess space will be ignored.
5984 */
5985
5986 /* Local Variables: */
5987 AstFrame *frm; /* Pointer to base Frame in FrameSet */
5988 AstLineDef *a; /* Line from inside point to test point */
5989 AstLineDef *b; /* Polygon edge */
5990 AstPointSet *in_base; /* PointSet holding base Frame input positions*/
5991 AstPointSet *result; /* Pointer to output PointSet */
5992 AstPolygon *this; /* Pointer to Polygon */
5993 double **ptr_in; /* Pointer to input base Frame coordinate data */
5994 double **ptr_out; /* Pointer to output current Frame coordinate data */
5995 double *px; /* Pointer to array of first axis values */
5996 double *py; /* Pointer to array of second axis values */
5997 double p[ 2 ]; /* Current test position */
5998 int closed; /* Is the boundary part of the Region? */
5999 int i; /* Edge index */
6000 int icoord; /* Coordinate index */
6001 int in_region; /* Is the point inside the Region? */
6002 int ncoord_out; /* No. of current Frame axes */
6003 int ncross; /* Number of crossings */
6004 int neg; /* Has the Region been negated? */
6005 int npoint; /* No. of input points */
6006 int nv; /* No. of vertices */
6007 int point; /* Loop counter for input points */
6008 int pos; /* Is test position in, on, or outside boundary? */
6009
6010 /* Check the global error status. */
6011 if ( !astOK ) return NULL;
6012
6013 /* Obtain a pointer to the Polygon structure. */
6014 this = (AstPolygon *) this_mapping;
6015
6016 /* Apply the parent mapping using the stored pointer to the Transform member
6017 function inherited from the parent Region class. This function validates
6018 all arguments and generates an output PointSet if necessary,
6019 containing a copy of the input PointSet. */
6020 result = (*parent_transform)( this_mapping, in, forward, out, status );
6021
6022 /* Get the number of points to be transformed. */
6023 npoint = astGetNpoint( result );
6024
6025 /* Get a pointer to the output axis values. */
6026 ptr_out = astGetPoints( result );
6027
6028 /* Find the number of axes in the current Frame. This need not be 2 (the
6029 number of axes in the *base* Frame must be 2 however). */
6030 ncoord_out = astGetNcoord( result );
6031
6032 /* We will now extend the parent astTransform method by performing the
6033 calculations needed to generate the output coordinate values. */
6034
6035 /* First use the encapsulated FrameSet to transform the supplied positions
6036 from the current Frame in the encapsulated FrameSet (the Frame
6037 represented by the Region), to the base Frame (the Frame in which the
6038 Region is defined). This call also returns a pointer to the base Frame
6039 of the encapsulated FrameSet. Note, the returned pointer may be a
6040 clone of the "in" pointer, and so we must be carefull not to modify the
6041 contents of the returned PointSet. */
6042 in_base = astRegTransform( this, in, 0, NULL, &frm );
6043 ptr_in = astGetPoints( in_base );
6044
6045 /* Get the number of vertices in the polygon. */
6046 nv = astGetNpoint( ((AstRegion *) this)->points );
6047
6048 /* See if the boundary is part of the Region. */
6049 closed = astGetClosed( this );
6050
6051 /* See if the Region has been negated. */
6052 neg = astGetNegated( this );
6053
6054 /* Perform coordinate arithmetic. */
6055 /* ------------------------------ */
6056 if ( astOK ) {
6057 px = ptr_in[ 0 ];
6058 py = ptr_in[ 1 ];
6059
6060 /* Loop round each supplied point in the base Frame of the polygon. */
6061 for ( point = 0; point < npoint; point++, px++, py++ ) {
6062
6063 /* If the input point is bad, indicate that bad output values should be
6064 returned. */
6065 if( *px == AST__BAD || *py == AST__BAD ) {
6066 in_region = 0;
6067
6068 /* Otherwise, we first determine if the point is inside, outside, or on,
6069 the Polygon boundary. Initialially it is unknown. */
6070 } else {
6071
6072 /* Ensure cached information is available.*/
6073 Cache( this, status );
6074
6075 /* Create a definition of the line from a point which is inside the
6076 polygon to the supplied point. This is a structure which includes
6077 cached intermediate information which can be used to speed up
6078 subsequent calculations. */
6079 p[ 0 ] = *px;
6080 p[ 1 ] = *py;
6081 a = astLineDef( frm, this->in, p );
6082
6083 /* We now determine the number of times this line crosses the polygon
6084 boundary. Initialise the number of crossings to zero. */
6085 ncross = 0;
6086 pos = UNKNOWN;
6087
6088 /* Loop rouind all edges of the polygon. */
6089 for( i = 0; i < nv; i++ ) {
6090 b = this->edges[ i ];
6091
6092 /* If this point is on the current edge, then we need do no more checks
6093 since we know it is either inside or outside the polygon (depending on
6094 whether the polygon is closed or not). */
6095 if( astLineContains( frm, b, 0, p ) ) {
6096 pos = ON;
6097 break;
6098
6099 /* Otherwise, see if the two lines cross within their extent. If so,
6100 increment the number of crossings. */
6101 } else if( astLineCrossing( frm, b, a, NULL ) ) {
6102 ncross++;
6103 }
6104 }
6105
6106 /* Free resources */
6107 a = astFree( a );
6108
6109 /* If the position is not on the boundary, it is inside the boundary if
6110 the number of crossings is even, and outside otherwise. */
6111 if( pos == UNKNOWN ) pos = ( ncross % 2 == 0 )? IN : OUT;
6112
6113 /* Whether the point is in the Region depends on whether the point is
6114 inside the polygon boundary, whether the Polygon has been negated, and
6115 whether the polygon is closed. */
6116 if( neg ) {
6117 if( pos == IN ) {
6118 in_region = 0;
6119 } else if( pos == OUT ) {
6120 in_region = 1;
6121 } else if( closed ) {
6122 in_region = 1;
6123 } else {
6124 in_region = 0;
6125 }
6126
6127 } else {
6128 if( pos == IN ) {
6129 in_region = 1;
6130 } else if( pos == OUT ) {
6131 in_region = 0;
6132 } else if( closed ) {
6133 in_region = 1;
6134 } else {
6135 in_region = 0;
6136 }
6137 }
6138 }
6139
6140 /* If the point is not inside the Region, store bad output values. */
6141 if( !in_region ) {
6142 for ( icoord = 0; icoord < ncoord_out; icoord++ ) {
6143 ptr_out[ icoord ][ point ] = AST__BAD;
6144 }
6145 }
6146 }
6147 }
6148
6149 /* Free resources */
6150 in_base = astAnnul( in_base );
6151 frm = astAnnul( frm );
6152
6153 /* Annul the result if an error has occurred. */
6154 if( !astOK ) result = astAnnul( result );
6155
6156 /* Return a pointer to the output PointSet. */
6157 return result;
6158 }
6159
6160 /* Functions which access class attributes. */
6161 /* ---------------------------------------- */
6162 /* Implement member functions to access the attributes associated with
6163 this class using the macros defined for this purpose in the
6164 "object.h" file. For a description of each attribute, see the class
6165 interface (in the associated .h file). */
6166
6167 /*
6168 *att++
6169 * Name:
6170 * SimpVertices
6171
6172 * Purpose:
6173 * Simplify a Polygon by transforming its vertices?
6174
6175 * Type:
6176 * Public attribute.
6177
6178 * Synopsis:
6179 * Integer (boolean).
6180
6181 * Description:
6182 * This attribute controls the behaviour of the
6183 c astSimplify
6184 f AST_SIMPLIFY
6185 * method when applied to a Polygon. The simplified Polygon is created
6186 * by transforming the vertices from the Frame in which the Polygon
6187 * was originally defined into the Frame currently represented by the
6188 * Polygon. If SimpVertices is non-zero (the default) then this
6189 * simplified Polygon is returned without further checks. If SimpVertices
6190 * is zero, a check is made that the edges of the new Polygon do not
6191 * depart significantly from the edges of the original Polygon (as
6192 * determined by the uncertainty associated with the Polygon). This
6193 * could occur, for instance, if the Mapping frrm the original to the
6194 * current Frame is highly non-linear. If this check fails, the
6195 * original unsimplified Polygon is returned without change.
6196
6197 * Applicability:
6198 * Polygon
6199 * All Polygons have this attribute.
6200
6201 *att--
6202 */
6203 astMAKE_CLEAR(Polygon,SimpVertices,simp_vertices,-INT_MAX)
6204 astMAKE_GET(Polygon,SimpVertices,int,0,( ( this->simp_vertices != -INT_MAX ) ?
6205 this->simp_vertices : 1 ))
6206 astMAKE_SET(Polygon,SimpVertices,int,simp_vertices,( value != 0 ))
6207 astMAKE_TEST(Polygon,SimpVertices,( this->simp_vertices != -INT_MAX ))
6208
6209 /* Copy constructor. */
6210 /* ----------------- */
Copy(const AstObject * objin,AstObject * objout,int * status)6211 static void Copy( const AstObject *objin, AstObject *objout, int *status ) {
6212 /*
6213 * Name:
6214 * Copy
6215
6216 * Purpose:
6217 * Copy constructor for Polygon objects.
6218
6219 * Type:
6220 * Private function.
6221
6222 * Synopsis:
6223 * void Copy( const AstObject *objin, AstObject *objout, int *status )
6224
6225 * Description:
6226 * This function implements the copy constructor for Polygon objects.
6227
6228 * Parameters:
6229 * objin
6230 * Pointer to the object to be copied.
6231 * objout
6232 * Pointer to the object being constructed.
6233 * status
6234 * Pointer to the inherited status variable.
6235
6236 * Notes:
6237 * - This constructor makes a deep copy.
6238 */
6239
6240 /* Local Variables: */
6241 AstPolygon *out; /* Pointer to output Polygon */
6242
6243 /* Check the global error status. */
6244 if ( !astOK ) return;
6245
6246 /* Obtain pointers to the output Polygon. */
6247 out = (AstPolygon *) objout;
6248
6249 /* For safety, first clear any references to the input memory from
6250 the output Polygon. */
6251 out->edges = NULL;
6252 out->startsat = NULL;
6253
6254 /* Indicate cached information needs nre-calculating. */
6255 astResetCache( (AstPolygon *) out );
6256 }
6257
6258
6259 /* Destructor. */
6260 /* ----------- */
Delete(AstObject * obj,int * status)6261 static void Delete( AstObject *obj, int *status ) {
6262 /*
6263 * Name:
6264 * Delete
6265
6266 * Purpose:
6267 * Destructor for Polygon objects.
6268
6269 * Type:
6270 * Private function.
6271
6272 * Synopsis:
6273 * void Delete( AstObject *obj, int *status )
6274
6275 * Description:
6276 * This function implements the destructor for Polygon objects.
6277
6278 * Parameters:
6279 * obj
6280 * Pointer to the object to be deleted.
6281 * status
6282 * Pointer to the inherited status variable.
6283
6284 * Notes:
6285 * This function attempts to execute even if the global error status is
6286 * set.
6287 */
6288
6289 /* Local Variables: */
6290 AstPointSet *ps; /* Pointer to PointSet inside Region */
6291 AstPolygon *this; /* Pointer to Polygon */
6292 int i; /* Index of vertex */
6293 int istat; /* Original AST error status */
6294 int nv; /* Number of vertices */
6295 int rep; /* Original error reporting state */
6296
6297 /* Obtain a pointer to the Polygon structure. */
6298 this = (AstPolygon *) obj;
6299
6300 /* Annul all resources. */
6301 ps = ((AstRegion *) this)->points;
6302 if( this->edges && ps ) {
6303
6304 /* Ensure we get a value for "nv" even if an error has occurred. */
6305 istat = astStatus;
6306 astClearStatus;
6307 rep = astReporting( 0 );
6308
6309 nv = astGetNpoint( ps );
6310
6311 astSetStatus( istat );
6312 astReporting( rep );
6313
6314 /* Free the structures holding the edge information. */
6315 for( i = 0; i < nv; i++ ) {
6316 this->edges[ i ] = astFree( this->edges[ i ] );
6317 }
6318 this->edges = astFree( this->edges );
6319 this->startsat = astFree( this->startsat );
6320
6321 }
6322 }
6323
6324 /* Dump function. */
6325 /* -------------- */
Dump(AstObject * this_object,AstChannel * channel,int * status)6326 static void Dump( AstObject *this_object, AstChannel *channel, int *status ) {
6327 /*
6328 * Name:
6329 * Dump
6330
6331 * Purpose:
6332 * Dump function for Polygon objects.
6333
6334 * Type:
6335 * Private function.
6336
6337 * Synopsis:
6338 * void Dump( AstObject *this, AstChannel *channel, int *status )
6339
6340 * Description:
6341 * This function implements the Dump function which writes out data
6342 * for the Polygon class to an output Channel.
6343
6344 * Parameters:
6345 * this
6346 * Pointer to the Polygon whose data are being written.
6347 * channel
6348 * Pointer to the Channel to which the data are being written.
6349 * status
6350 * Pointer to the inherited status variable.
6351 */
6352
6353 /* Local Variables: */
6354 AstPolygon *this; /* Pointer to the Polygon structure */
6355 int ival; /* Integer attribute value */
6356 int set; /* Attribute value set? */
6357
6358 /* Check the global error status. */
6359 if ( !astOK ) return;
6360
6361 /* Obtain a pointer to the Polygon structure. */
6362 this = (AstPolygon *) this_object;
6363
6364 /* Write out values representing the instance variables for the
6365 Polygon class. Accompany these with appropriate comment strings,
6366 possibly depending on the values being written.*/
6367
6368 /* In the case of attributes, we first use the appropriate (private)
6369 Test... member function to see if they are set. If so, we then use
6370 the (private) Get... function to obtain the value to be written
6371 out.
6372
6373 For attributes which are not set, we use the astGet... method to
6374 obtain the value instead. This will supply a default value
6375 (possibly provided by a derived class which over-rides this method)
6376 which is more useful to a human reader as it corresponds to the
6377 actual default attribute value. Since "set" will be zero, these
6378 values are for information only and will not be read back. */
6379
6380 /* SimpVertices. */
6381 /* ------------ */
6382 /* Write out the forward-inverse simplification flag. */
6383 set = TestSimpVertices( this, status );
6384 ival = set ? GetSimpVertices( this, status ) : astGetSimpVertices( this );
6385 astWriteInt( channel, "SimpVT", set, 0, ival, "Simplify by transforming vertices?" );
6386
6387 /* A flag indicating the convention used for determining the interior of
6388 the polygon. A zero value indicates that the old AST system is in
6389 use (inside to the left when moving anti-clockwise round the vertices
6390 as viewed from the outside of the celestial sphere). A non-zero value
6391 indicates the STC system is in use (inside to the left when moving
6392 anti-clockwise round the vertices as viewed from the inside of the
6393 celestial sphere). AST currently uses the STC system. */
6394 astWriteInt( channel, "Order", 1, 0, 1, "Polygon uses STC vertex order convention" );
6395 }
6396
6397 /* Standard class functions. */
6398 /* ========================= */
6399 /* Implement the astIsAPolygon and astCheckPolygon functions using the macros
6400 defined for this purpose in the "object.h" header file. */
astMAKE_ISA(Polygon,Region)6401 astMAKE_ISA(Polygon,Region)
6402 astMAKE_CHECK(Polygon)
6403
6404 AstPolygon *astPolygon_( void *frame_void, int npnt, int dim,
6405 const double *points, AstRegion *unc,
6406 const char *options, int *status, ...) {
6407 /*
6408 *++
6409 * Name:
6410 c astPolygon
6411 f AST_POLYGON
6412
6413 * Purpose:
6414 * Create a Polygon.
6415
6416 * Type:
6417 * Public function.
6418
6419 * Synopsis:
6420 c #include "polygon.h"
6421 c AstPolygon *astPolygon( AstFrame *frame, int npnt, int dim,
6422 c const double *points, AstRegion *unc,
6423 c const char *options, ... )
6424 f RESULT = AST_POLYGON( FRAME, NPNT, DIM, POINTS, UNC, OPTIONS, STATUS )
6425
6426 * Class Membership:
6427 * Polygon constructor.
6428
6429 * Description:
6430 * This function creates a new Polygon object and optionally initialises
6431 * its attributes.
6432 *
6433 * The Polygon class implements a polygonal area, defined by a
6434 * collection of vertices, within a 2-dimensional Frame. The vertices
6435 * are connected together by geodesic curves within the encapsulated Frame.
6436 * For instance, if the encapsulated Frame is a simple Frame then the
6437 * geodesics will be straight lines, but if the Frame is a SkyFrame then
6438 * the geodesics will be great circles. Note, the vertices must be
6439 * supplied in an order such that the inside of the polygon is to the
6440 * left of the boundary as the vertices are traversed. Supplying them
6441 * in the reverse order will effectively negate the polygon.
6442 *
6443 * Within a SkyFrame, neighbouring vertices are always joined using the
6444 * shortest path. Thus if an edge of 180 degrees or more in length is
6445 * required, it should be split into section each of which is less
6446 * than 180 degrees. The closed path joining all the vertices in order
6447 * will divide the celestial sphere into two disjoint regions. The
6448 * inside of the polygon is the region which is circled in an
6449 * anti-clockwise manner (when viewed from the inside of the celestial
6450 * sphere) when moving through the list of vertices in the order in
6451 * which they were supplied when the Polygon was created (i.e. the
6452 * inside is to the left of the boundary when moving through the
6453 * vertices in the order supplied).
6454
6455 * Parameters:
6456 c frame
6457 f FRAME = INTEGER (Given)
6458 * A pointer to the Frame in which the region is defined. It must
6459 * have exactly 2 axes. A deep copy is taken of the supplied Frame.
6460 * This means that any subsequent changes made to the Frame using the
6461 * supplied pointer will have no effect the Region.
6462 c npnt
6463 f NPNT = INTEGER (Given)
6464 * The number of points in the Region.
6465 c dim
6466 f DIM = INTEGER (Given)
6467 c The number of elements along the second dimension of the "points"
6468 f The number of elements along the first dimension of the POINTS
6469 * array (which contains the point coordinates). This value is
6470 * required so that the coordinate values can be correctly
6471 * located if they do not entirely fill this array. The value
6472 c given should not be less than "npnt".
6473 f given should not be less than NPNT.
6474 c points
6475 f POINTS( DIM, 2 ) = DOUBLE PRECISION (Given)
6476 c The address of the first element of a 2-dimensional array of
6477 c shape "[2][dim]" giving the physical coordinates of the vertices.
6478 c These should be stored such that the value of coordinate
6479 c number "coord" for point number "pnt" is found in element
6480 c "in[coord][pnt]".
6481 f A 2-dimensional array giving the physical coordinates of the
6482 f vertices. These should be stored such that the value of coordinate
6483 f number COORD for point number PNT is found in element IN(PNT,COORD).
6484 c unc
6485 f UNC = INTEGER (Given)
6486 * An optional pointer to an existing Region which specifies the
6487 * uncertainties associated with the boundary of the Box being created.
6488 * The uncertainty in any point on the boundary of the Box is found by
6489 * shifting the supplied "uncertainty" Region so that it is centred at
6490 * the boundary point being considered. The area covered by the
6491 * shifted uncertainty Region then represents the uncertainty in the
6492 * boundary position. The uncertainty is assumed to be the same for
6493 * all points.
6494 *
6495 * If supplied, the uncertainty Region must be of a class for which
6496 * all instances are centro-symetric (e.g. Box, Circle, Ellipse, etc.)
6497 * or be a Prism containing centro-symetric component Regions. A deep
6498 * copy of the supplied Region will be taken, so subsequent changes to
6499 * the uncertainty Region using the supplied pointer will have no
6500 * effect on the created Box. Alternatively,
6501 f a null Object pointer (AST__NULL)
6502 c a NULL Object pointer
6503 * may be supplied, in which case a default uncertainty is used
6504 * equivalent to a box 1.0E-6 of the size of the Box being created.
6505 *
6506 * The uncertainty Region has two uses: 1) when the
6507 c astOverlap
6508 f AST_OVERLAP
6509 * function compares two Regions for equality the uncertainty
6510 * Region is used to determine the tolerance on the comparison, and 2)
6511 * when a Region is mapped into a different coordinate system and
6512 * subsequently simplified (using
6513 c astSimplify),
6514 f AST_SIMPLIFY),
6515 * the uncertainties are used to determine if the transformed boundary
6516 * can be accurately represented by a specific shape of Region.
6517 c options
6518 f OPTIONS = CHARACTER * ( * ) (Given)
6519 c Pointer to a null-terminated string containing an optional
6520 c comma-separated list of attribute assignments to be used for
6521 c initialising the new Polygon. The syntax used is identical to
6522 c that for the astSet function and may include "printf" format
6523 c specifiers identified by "%" symbols in the normal way.
6524 f A character string containing an optional comma-separated
6525 f list of attribute assignments to be used for initialising the
6526 f new Polygon. The syntax used is identical to that for the
6527 f AST_SET routine.
6528 c ...
6529 c If the "options" string contains "%" format specifiers, then
6530 c an optional list of additional arguments may follow it in
6531 c order to supply values to be substituted for these
6532 c specifiers. The rules for supplying these are identical to
6533 c those for the astSet function (and for the C "printf"
6534 c function).
6535 f STATUS = INTEGER (Given and Returned)
6536 f The global status.
6537
6538 * Returned Value:
6539 c astPolygon()
6540 f AST_POLYGON = INTEGER
6541 * A pointer to the new Polygon.
6542
6543 * Notes:
6544 * - A null Object pointer (AST__NULL) will be returned if this
6545 c function is invoked with the AST error status set, or if it
6546 f function is invoked with STATUS set to an error value, or if it
6547 * should fail for any reason.
6548
6549 * Status Handling:
6550 * The protected interface to this function includes an extra
6551 * parameter at the end of the parameter list descirbed above. This
6552 * parameter is a pointer to the integer inherited status
6553 * variable: "int *status".
6554
6555 *--
6556 */
6557
6558 /* Local Variables: */
6559 astDECLARE_GLOBALS /* Pointer to thread-specific global data */
6560 AstFrame *frame; /* Pointer to Frame structure */
6561 AstPolygon *new; /* Pointer to new Polygon */
6562 va_list args; /* Variable argument list */
6563
6564 /* Get a pointer to the thread specific global data structure. */
6565 astGET_GLOBALS(NULL);
6566
6567 /* Check the global status. */
6568 if ( !astOK ) return NULL;
6569
6570 /* Obtain and validate a pointer to the supplied Frame structure. */
6571 frame = astCheckFrame( frame_void );
6572
6573 /* Initialise the Polygon, allocating memory and initialising the
6574 virtual function table as well if necessary. */
6575 new = astInitPolygon( NULL, sizeof( AstPolygon ), !class_init,
6576 &class_vtab, "Polygon", frame, npnt,
6577 dim, points, unc );
6578
6579 /* If successful, note that the virtual function table has been
6580 initialised. */
6581 if ( astOK ) {
6582 class_init = 1;
6583
6584 /* Obtain the variable argument list and pass it along with the options string
6585 to the astVSet method to initialise the new Polygon's attributes. */
6586 va_start( args, status );
6587 astVSet( new, options, NULL, args );
6588 va_end( args );
6589
6590 /* If an error occurred, clean up by deleting the new object. */
6591 if ( !astOK ) new = astDelete( new );
6592 }
6593
6594 /* Return a pointer to the new Polygon. */
6595 return new;
6596 }
6597
astPolygonId_(void * frame_void,int npnt,int dim,const double * points,void * unc_void,const char * options,...)6598 AstPolygon *astPolygonId_( void *frame_void, int npnt, int dim,
6599 const double *points, void *unc_void,
6600 const char *options, ... ) {
6601 /*
6602 * Name:
6603 * astPolygonId_
6604
6605 * Purpose:
6606 * Create a Polygon.
6607
6608 * Type:
6609 * Private function.
6610
6611 * Synopsis:
6612 * #include "polygon.h"
6613 * AstPolygon *astPolygonId_( void *frame_void, int npnt,
6614 * int dim, const double *points, void *unc_void,
6615 * const char *options, ... )
6616
6617 * Class Membership:
6618 * Polygon constructor.
6619
6620 * Description:
6621 * This function implements the external (public) interface to the
6622 * astPolygon constructor function. It returns an ID value (instead
6623 * of a true C pointer) to external users, and must be provided
6624 * because astPolygon_ has a variable argument list which cannot be
6625 * encapsulated in a macro (where this conversion would otherwise
6626 * occur).
6627 *
6628 * The variable argument list also prevents this function from
6629 * invoking astPolygon_ directly, so it must be a re-implementation
6630 * of it in all respects, except for the final conversion of the
6631 * result to an ID value.
6632
6633 * Parameters:
6634 * As for astPolygon_.
6635
6636 * Returned Value:
6637 * The ID value associated with the new Polygon.
6638 */
6639
6640 /* Local Variables: */
6641 astDECLARE_GLOBALS /* Pointer to thread-specific global data */
6642 AstFrame *frame; /* Pointer to Frame structure */
6643 AstPolygon *new; /* Pointer to new Polygon */
6644 AstRegion *unc; /* Pointer to Region structure */
6645 va_list args; /* Variable argument list */
6646
6647 int *status; /* Get a pointer to the thread specific global data structure. */
6648 astGET_GLOBALS(NULL);
6649
6650 /* Get a pointer to the inherited status value. */
6651 status = astGetStatusPtr;
6652
6653 /* Check the global status. */
6654 if ( !astOK ) return NULL;
6655
6656 /* Obtain a Frame pointer from the supplied ID and validate the
6657 pointer to ensure it identifies a valid Frame. */
6658 frame = astVerifyFrame( astMakePointer( frame_void ) );
6659
6660 /* Obtain a Region pointer from the supplied "unc" ID and validate the
6661 pointer to ensure it identifies a valid Region . */
6662 unc = unc_void ? astCheckRegion( astMakePointer( unc_void ) ) : NULL;
6663
6664 /* Initialise the Polygon, allocating memory and initialising the
6665 virtual function table as well if necessary. */
6666 new = astInitPolygon( NULL, sizeof( AstPolygon ), !class_init,
6667 &class_vtab, "Polygon", frame, npnt, dim,
6668 points, unc );
6669
6670 /* If successful, note that the virtual function table has been
6671 initialised. */
6672 if ( astOK ) {
6673 class_init = 1;
6674
6675 /* Obtain the variable argument list and pass it along with the options string
6676 to the astVSet method to initialise the new Polygon's attributes. */
6677 va_start( args, options );
6678 astVSet( new, options, NULL, args );
6679 va_end( args );
6680
6681 /* If an error occurred, clean up by deleting the new object. */
6682 if ( !astOK ) new = astDelete( new );
6683 }
6684
6685 /* Return an ID value for the new Polygon. */
6686 return astMakeId( new );
6687 }
6688
6689
astInitPolygon_(void * mem,size_t size,int init,AstPolygonVtab * vtab,const char * name,AstFrame * frame,int npnt,int dim,const double * points,AstRegion * unc,int * status)6690 AstPolygon *astInitPolygon_( void *mem, size_t size, int init, AstPolygonVtab *vtab,
6691 const char *name, AstFrame *frame, int npnt,
6692 int dim, const double *points, AstRegion *unc, int *status ) {
6693 /*
6694 *+
6695 * Name:
6696 * astInitPolygon
6697
6698 * Purpose:
6699 * Initialise a Polygon.
6700
6701 * Type:
6702 * Protected function.
6703
6704 * Synopsis:
6705 * #include "polygon.h"
6706 * AstPolygon *astInitPolygon( void *mem, size_t size, int init, AstPolygonVtab *vtab,
6707 * const char *name, AstFrame *frame, int npnt,
6708 * int dim, const double *points, AstRegion *unc )
6709
6710 * Class Membership:
6711 * Polygon initialiser.
6712
6713 * Description:
6714 * This function is provided for use by class implementations to initialise
6715 * a new Polygon object. It allocates memory (if necessary) to accommodate
6716 * the Polygon plus any additional data associated with the derived class.
6717 * It then initialises a Polygon structure at the start of this memory. If
6718 * the "init" flag is set, it also initialises the contents of a virtual
6719 * function table for a Polygon at the start of the memory passed via the
6720 * "vtab" parameter.
6721
6722 * Parameters:
6723 * mem
6724 * A pointer to the memory in which the Polygon is to be initialised.
6725 * This must be of sufficient size to accommodate the Polygon data
6726 * (sizeof(Polygon)) plus any data used by the derived class. If a value
6727 * of NULL is given, this function will allocate the memory itself using
6728 * the "size" parameter to determine its size.
6729 * size
6730 * The amount of memory used by the Polygon (plus derived class data).
6731 * This will be used to allocate memory if a value of NULL is given for
6732 * the "mem" parameter. This value is also stored in the Polygon
6733 * structure, so a valid value must be supplied even if not required for
6734 * allocating memory.
6735 * init
6736 * A logical flag indicating if the Polygon's virtual function table is
6737 * to be initialised. If this value is non-zero, the virtual function
6738 * table will be initialised by this function.
6739 * vtab
6740 * Pointer to the start of the virtual function table to be associated
6741 * with the new Polygon.
6742 * name
6743 * Pointer to a constant null-terminated character string which contains
6744 * the name of the class to which the new object belongs (it is this
6745 * pointer value that will subsequently be returned by the astGetClass
6746 * method).
6747 * frame
6748 * A pointer to the Frame in which the region is defined.
6749 * npnt
6750 * The number of points in the Region.
6751 * dim
6752 * The number of elements along the second dimension of the "points"
6753 * array (which contains the point coordinates). This value is
6754 * required so that the coordinate values can be correctly
6755 * located if they do not entirely fill this array. The value
6756 * given should not be less than "npnt".
6757 * points
6758 * The address of the first element of a 2-dimensional array of
6759 * shape "[2][dim]" giving the physical coordinates of the
6760 * points. These should be stored such that the value of coordinate
6761 * number "coord" for point number "pnt" is found in element
6762 * "in[coord][pnt]".
6763 * unc
6764 * A pointer to a Region which specifies the uncertainty in the
6765 * supplied positions (all points in the new Polygon being
6766 * initialised are assumed to have the same uncertainty). A NULL
6767 * pointer can be supplied, in which case default uncertainties equal
6768 * to 1.0E-6 of the dimensions of the new Polygon's bounding box are
6769 * used. If an uncertainty Region is supplied, it must be either a Box,
6770 * a Circle or an Ellipse, and its encapsulated Frame must be related
6771 * to the Frame supplied for parameter "frame" (i.e. astConvert
6772 * should be able to find a Mapping between them). Two positions
6773 * the "frame" Frame are considered to be co-incident if their
6774 * uncertainty Regions overlap. The centre of the supplied
6775 * uncertainty Region is immaterial since it will be re-centred on the
6776 * point being tested before use. A deep copy is taken of the supplied
6777 * Region.
6778
6779 * Returned Value:
6780 * A pointer to the new Polygon.
6781
6782 * Notes:
6783 * - A null pointer will be returned if this function is invoked with the
6784 * global error status set, or if it should fail for any reason.
6785 *-
6786 */
6787
6788 /* Local Variables: */
6789 AstPolygon *new; /* Pointer to new Polygon */
6790 AstPointSet *pset; /* Pointer to PointSet holding points */
6791 const double *q; /* Pointer to next supplied axis value */
6792 double **ptr; /* Pointer to data in pset */
6793 double *p; /* Pointer to next PointSet axis value */
6794 int i; /* Axis index */
6795 int j; /* Point index */
6796 int nin; /* No. of axes */
6797
6798 /* Check the global status. */
6799 if ( !astOK ) return NULL;
6800
6801 /* If necessary, initialise the virtual function table. */
6802 if ( init ) astInitPolygonVtab( vtab, name );
6803
6804 /* Initialise. */
6805 new = NULL;
6806
6807 /* Check the number of axis values per position is correct. */
6808 nin = astGetNaxes( frame );
6809 if( nin != 2 ) {
6810 astError( AST__BADIN, "astInitPolygon(%s): The supplied %s has %d "
6811 "axes - polygons must have exactly 2 axes.", status, name,
6812 astGetClass( frame ), nin );
6813
6814 /* If so create a PointSet and store the supplied points in it. Check
6815 none are bad. */
6816 } else {
6817 pset = astPointSet( npnt, 2, "", status );
6818 ptr = astGetPoints( pset );
6819 for( i = 0; i < 2 && astOK; i++ ) {
6820 p = ptr[ i ];
6821 q = points + i*dim;
6822 for( j = 0; j < npnt; j++ ) {
6823 if( (*(p++) = *(q++)) == AST__BAD ) {
6824 astError( AST__BADIN, "astInitPolygon(%s): One or more "
6825 "bad axis values supplied for the vertex "
6826 "number %d.", status, name, j + 1 );
6827 break;
6828 }
6829 }
6830 }
6831
6832 /* Initialise a Region structure (the parent class) as the first component
6833 within the Polygon structure, allocating memory if necessary. */
6834 new = (AstPolygon *) astInitRegion( mem, size, 0, (AstRegionVtab *) vtab,
6835 name, frame, pset, unc );
6836 if ( astOK ) {
6837
6838 /* Initialise the Polygon data. */
6839 /* ------------------------------ */
6840 new->lbnd[ 0 ] = AST__BAD;
6841 new->ubnd[ 0 ] = AST__BAD;
6842 new->lbnd[ 1 ] = AST__BAD;
6843 new->ubnd[ 1 ] = AST__BAD;
6844 new->simp_vertices = -INT_MAX;
6845 new->edges = NULL;
6846 new->startsat = NULL;
6847 new->totlen = 0.0;
6848 new->acw = 1;
6849 new->stale = 1;
6850
6851 /* Ensure the vertices are stored such that the unnegated Polygon
6852 represents the inside of the polygon. */
6853 EnsureInside( new, status );
6854
6855 /* If an error occurred, clean up by deleting the new Polygon. */
6856 if ( !astOK ) new = astDelete( new );
6857 }
6858
6859 /* Free resources. */
6860 pset = astAnnul( pset );
6861
6862 }
6863
6864 /* Return a pointer to the new Polygon. */
6865 return new;
6866 }
6867
astLoadPolygon_(void * mem,size_t size,AstPolygonVtab * vtab,const char * name,AstChannel * channel,int * status)6868 AstPolygon *astLoadPolygon_( void *mem, size_t size, AstPolygonVtab *vtab,
6869 const char *name, AstChannel *channel, int *status ) {
6870 /*
6871 *+
6872 * Name:
6873 * astLoadPolygon
6874
6875 * Purpose:
6876 * Load a Polygon.
6877
6878 * Type:
6879 * Protected function.
6880
6881 * Synopsis:
6882 * #include "polygon.h"
6883 * AstPolygon *astLoadPolygon( void *mem, size_t size, AstPolygonVtab *vtab,
6884 * const char *name, AstChannel *channel )
6885
6886 * Class Membership:
6887 * Polygon loader.
6888
6889 * Description:
6890 * This function is provided to load a new Polygon using data read
6891 * from a Channel. It first loads the data used by the parent class
6892 * (which allocates memory if necessary) and then initialises a
6893 * Polygon structure in this memory, using data read from the input
6894 * Channel.
6895 *
6896 * If the "init" flag is set, it also initialises the contents of a
6897 * virtual function table for a Polygon at the start of the memory
6898 * passed via the "vtab" parameter.
6899
6900 * Parameters:
6901 * mem
6902 * A pointer to the memory into which the Polygon is to be
6903 * loaded. This must be of sufficient size to accommodate the
6904 * Polygon data (sizeof(Polygon)) plus any data used by derived
6905 * classes. If a value of NULL is given, this function will
6906 * allocate the memory itself using the "size" parameter to
6907 * determine its size.
6908 * size
6909 * The amount of memory used by the Polygon (plus derived class
6910 * data). This will be used to allocate memory if a value of
6911 * NULL is given for the "mem" parameter. This value is also
6912 * stored in the Polygon structure, so a valid value must be
6913 * supplied even if not required for allocating memory.
6914 *
6915 * If the "vtab" parameter is NULL, the "size" value is ignored
6916 * and sizeof(AstPolygon) is used instead.
6917 * vtab
6918 * Pointer to the start of the virtual function table to be
6919 * associated with the new Polygon. If this is NULL, a pointer
6920 * to the (static) virtual function table for the Polygon class
6921 * is used instead.
6922 * name
6923 * Pointer to a constant null-terminated character string which
6924 * contains the name of the class to which the new object
6925 * belongs (it is this pointer value that will subsequently be
6926 * returned by the astGetClass method).
6927 *
6928 * If the "vtab" parameter is NULL, the "name" value is ignored
6929 * and a pointer to the string "Polygon" is used instead.
6930
6931 * Returned Value:
6932 * A pointer to the new Polygon.
6933
6934 * Notes:
6935 * - A null pointer will be returned if this function is invoked
6936 * with the global error status set, or if it should fail for any
6937 * reason.
6938 *-
6939 */
6940
6941 /* Local Variables: */
6942 astDECLARE_GLOBALS /* Pointer to thread-specific global data */
6943 AstPolygon *new; /* Pointer to the new Polygon */
6944 int order; /* Is the new (STC) order convention used? */
6945
6946 /* Initialise. */
6947 new = NULL;
6948
6949 /* Check the global error status. */
6950 if ( !astOK ) return new;
6951
6952 /* Get a pointer to the thread specific global data structure. */
6953 astGET_GLOBALS(channel);
6954
6955 /* If a NULL virtual function table has been supplied, then this is
6956 the first loader to be invoked for this Polygon. In this case the
6957 Polygon belongs to this class, so supply appropriate values to be
6958 passed to the parent class loader (and its parent, etc.). */
6959 if ( !vtab ) {
6960 size = sizeof( AstPolygon );
6961 vtab = &class_vtab;
6962 name = "Polygon";
6963
6964 /* If required, initialise the virtual function table for this class. */
6965 if ( !class_init ) {
6966 astInitPolygonVtab( vtab, name );
6967 class_init = 1;
6968 }
6969 }
6970
6971 /* Invoke the parent class loader to load data for all the ancestral
6972 classes of the current one, returning a pointer to the resulting
6973 partly-built Polygon. */
6974 new = astLoadRegion( mem, size, (AstRegionVtab *) vtab, name,
6975 channel );
6976
6977 if ( astOK ) {
6978
6979 /* Read input data. */
6980 /* ================ */
6981 /* Request the input Channel to read all the input data appropriate to
6982 this class into the internal "values list". */
6983 astReadClassData( channel, "Polygon" );
6984
6985 /* Now read each individual data item from this list and use it to
6986 initialise the appropriate instance variable(s) for this class. */
6987
6988 /* In the case of attributes, we first read the "raw" input value,
6989 supplying the "unset" value as the default. If a "set" value is
6990 obtained, we then use the appropriate (private) Set... member
6991 function to validate and set the value properly. */
6992
6993 new->simp_vertices = astReadInt( channel, "simpvt", -INT_MAX );
6994 if ( TestSimpVertices( new, status ) ) SetSimpVertices( new, new->simp_vertices, status );
6995
6996 /* A flag indicating what order the vertices are stored in. See the Dump
6997 function. */
6998 order = astReadInt( channel, "order", 0 );
6999
7000 /* Initialise other class properties. */
7001 new->lbnd[ 0 ] = AST__BAD;
7002 new->ubnd[ 0 ] = AST__BAD;
7003 new->lbnd[ 1 ] = AST__BAD;
7004 new->ubnd[ 1 ] = AST__BAD;
7005 new->edges = NULL;
7006 new->startsat = NULL;
7007 new->totlen = 0.0;
7008 new->acw = 1;
7009 new->stale = 1;
7010
7011 /* If the order in which the vertices were written used the old AST
7012 convention, negate the Polygon so that it is consistent with the
7013 current conevtion (based on STC). */
7014 if( ! order ) astNegate( new );
7015
7016 /* Ensure the vertices are stored such that the unnegated Polygon
7017 represents the inside of the polygon. */
7018 EnsureInside( new, status );
7019
7020 /* If an error occurred, clean up by deleting the new Polygon. */
7021 if ( !astOK ) new = astDelete( new );
7022 }
7023
7024 /* Return the new Polygon pointer. */
7025 return new;
7026 }
7027
7028 /* Virtual function interfaces. */
7029 /* ============================ */
7030 /* These provide the external interface to the virtual functions defined by
7031 this class. Each simply checks the global error status and then locates and
7032 executes the appropriate member function, using the function pointer stored
7033 in the object's virtual function table (this pointer is located using the
7034 astMEMBER macro defined in "object.h").
7035
7036 Note that the member function may not be the one defined here, as it may
7037 have been over-ridden by a derived class. However, it should still have the
7038 same interface. */
7039
7040
astDownsize_(AstPolygon * this,double maxerr,int maxvert,int * status)7041 AstPolygon *astDownsize_( AstPolygon *this, double maxerr, int maxvert,
7042 int *status ) {
7043 if ( !astOK ) return NULL;
7044 return (**astMEMBER(this,Polygon,Downsize))( this, maxerr, maxvert, status );
7045 }
7046
7047
7048