Attachment 'stcschan-demo3.c'

Download

   1 /* Name:
   2       stcschan-demo3.c
   3 
   4    Purpose:
   5       A demonstration of the facilities provided by the AST library 
   6       for reading STC metadata encoded using the STC-S linear string 
   7       format.
   8 
   9    Description:
  10       This program reads an STC-S description from a text file, and also
  11       reads a set of 2-D spatial FITS-WCS headers from another (text) file. 
  12       It then opens a specified graphics device, and displays an annotated 
  13       coordinate grid covering the region described by the FITS headers.
  14       Finally, it draws the outline of the spatial extent of the STC-S 
  15       description over the top of the annotated coordinate grid.
  16 
  17    Usage:
  18       % stcschan-demo3 <stcs-file> <header-file> <device>
  19 
  20       <stcs-file>: The path to a text file containing the STC-S
  21       description.
  22 
  23       <header-file>: The path to a text file containing a set of FITS-WCS
  24       headers.
  25 
  26       <device>: The name of an available PGPLOT graphics device. If not
  27       supplied, the available device names will be listed and the program
  28       will then exit.
  29 
  30    Example:
  31       % stcschan-demo3 m31.stcs andromeda.head /xserve
  32 
  33    To compile and link:
  34       Assuming your starlink distribution is in "/star":
  35 
  36       % gcc -o stcschan-demo3 stcschan-demo3.c -g -L/star/lib \
  37              -I/star/include `ast_link -pgplot` -lcpgplot
  38 
  39 */
  40 
  41 /* Include system headers. */
  42 #include <stdio.h>
  43 #include <string.h>
  44 
  45 /* Include the PGPLOT header. */
  46 #include "cpgplot.h"
  47 
  48 /* Include the AST library header. */
  49 #include "ast.h"
  50 
  51 /* Maximum allowed length for a single line of text from the disk file. */
  52 #define MAX_LINE_LEN 100
  53 
  54 /* Prototypes */
  55 const char *source( void );
  56 AstFrameSet *ReadFitsHeaders( const char *, int *, int * );
  57 AstRegion *ReadStcs( const char * );
  58 
  59 
  60 int main( int argc, char **argv ){
  61 
  62 /* Local variables: */
  63    AstBox *pixbox;   
  64    AstFrame *pixfrm;
  65    AstFrameSet *align_fs = NULL;
  66    AstFrameSet *fset = NULL;
  67    AstPlot *plot;
  68    AstRegion *reg = NULL;
  69    AstRegion *wcsbox = NULL;
  70    AstRegion *wcsreg = NULL;
  71    const char *dev;
  72    double bbox[ 4 ];
  73    float gbox[ 4 ];
  74    int overlap_flag, status, naxis1, naxis2, ibase;
  75    
  76 /* Initialise the returned system status to indicate failure. */
  77    status = 1;
  78 
  79 /* Start an AST object context. This means we do not need to annull 
  80    each AST Object individually. Instead, all Objects created within 
  81    this context will be annulled automatically by the corresponding
  82    invocation of astEnd. */
  83    astBegin;
  84 
  85 /* Check there are enough command line arguments. */
  86    if( argc < 3 ) {
  87       printf( "Usage: stcschan-demo3 <stcs-file> <header-file> <device>\n" );
  88 
  89 /* If so, attempt to read the STC-S description, creating a corresponding 
  90    AST Region. If this is successful, attempt to read the FITS header file 
  91    and create an equivalent AST FrameSet. */
  92    } else {
  93       reg = ReadStcs( argv[ 1 ] );
  94       if( reg ) fset = ReadFitsHeaders( argv[ 2 ], &naxis1, &naxis2 );
  95    }
  96 
  97 /* Check we obtained a Region and a FrameSet successfully. */
  98    if( reg && fset ){
  99 
 100 /* Check that we can align the FITS WCS with the spatial axes in the
 101    STC-S description. AST contains various built-in conversions between 
 102    standard celestial coordinate system. The necessary conversion will be
 103    identified and used automatically if necessary to achieve alignment. 
 104    The returned object "align_fs" is a FrameSet that encapsulates the
 105    the FITS WCS coordinate system, the STC-S spatial coordinate system and
 106    the Mapping between them. A NULL pointer is returned if alignment is
 107    not possible. Note, the astConvert method changes the base Frame in
 108    any supplied FrameSet to indicate which coordinate frame was used for
 109    alignment. So we first note the original base Frame index and then
 110    re-instate it afterwards. */
 111       ibase = astGetI( fset, "Base" );
 112       align_fs = astConvert( reg, fset, " " ); 
 113       astSetI( fset, "Base", ibase );
 114 
 115       if( !align_fs ) {
 116          printf( "Could not align the FITS WCS with the spatial axes "
 117                  "in the STC-S\n" );
 118 
 119 /* If alignment was possible, use the alignment FrameSet to create a new 
 120    Region representing the same region as the supplied STC-S, but 
 121    expressed in the coordinate system of the FITS WCS. The FrameSet class
 122    inherits from both Frame and Mapping, and so the "align_fs" FrameSet can 
 123    be used both as the Mapping in this call, and as the Frame. When used
 124    as a Mapping, a FrameSet represents the transformation between its base
 125    and current Frame. When used as a Frame, a FrameSet represents its
 126    current Frame. */
 127       } else {
 128          wcsreg = astMapRegion( reg, align_fs, align_fs );
 129 
 130 /* It would be nice to warn the user if the STC-S AstroCoordsArea does
 131    not overlap the FITS grid. To do so we need a Region representing the
 132    FITS grid. Create one now (an AST Box). First store the bounds of the 
 133    FITS grid, in pixel coordinates (i.e. a system in which the centre of 
 134    the bottom left pixel is at (1.0,1.0) ). */
 135          bbox[ 0 ] = 0.5;
 136          bbox[ 1 ] = 0.5;
 137          bbox[ 2 ] = (double) naxis1 + 0.5;
 138          bbox[ 3 ] = (double) naxis2 + 0.5;
 139 
 140 /* Get a pointer to the Frame describing the FITS pixel coordinate system
 141    (the base Frame in the FITS FrameSet). */
 142          pixfrm = astGetFrame( fset, AST__BASE );
 143 
 144 /* Create a Box that encompasses the required range of axis values within
 145    the pixel coordinate Frame. */
 146          pixbox = astBox( pixfrm, 1, bbox, bbox + 2, NULL, " " );
 147 
 148 /* Create another Region that represents the same area, but in the FITS
 149    WCS. */
 150          wcsbox = astMapRegion( pixbox, fset, fset );
 151 
 152 /* If the previous step failed, it probably means the FITS header covers
 153    the entire sky, resulting the corners of the pixel grid having invalid
 154    sky positions. So cancel the error and omit the overlap test. */
 155          if( !wcsbox ) {
 156             astClearStatus;
 157             printf("\nContinuing, but omitting overlap test...\n\n");
 158                 
 159 /* Now see if the Region representing the FITS grid overlaps the region 
 160    read from the STC-S description.*/
 161          } else {
 162             overlap_flag = astOverlap( wcsreg, wcsbox );
 163 
 164             if( overlap_flag == 1 || overlap_flag == 6 ) {
 165                printf( "\nThere is no overlap between the FITS grid and the "
 166                        "STC-S AstroCoordsArea\n\n" );
 167    
 168             } else if( overlap_flag == 3 || overlap_flag == 5 ) {
 169                printf( "\nThe FITS grid is completely contained within the "
 170                        "STC-S AstroCoordsArea\n\n" );
 171             }
 172 
 173          }
 174       }
 175    }
 176 
 177 /* Check we obtained a FITS-WCS to STC-S Mapping, and that no error has
 178    occurred in AST. */
 179    if( wcsreg && astOK ){
 180 
 181 /* Open PGPLOT using the specified device. Prompt the user for a device if
 182    none was supplied on the command line. */
 183       if( argc < 4 ) {
 184          dev = "?";
 185       } else {
 186          dev = argv[ 3 ];
 187       }
 188       if( cpgbeg( 0, dev, 1, 1 ) == 1 ) {
 189 
 190 /* Clear the screen. */
 191          cpgpage();
 192 
 193 /* Ensure the graphics window has equal scales on both axes. */
 194          cpgwnad( 0.0f, 1.0f, 0.0f, 1.0f );
 195 
 196 /* Find the extent of the graphics window, and store in an array suitable
 197    for passing to the astPlot function. */
 198          cpgqwin( gbox, gbox + 2, gbox + 1, gbox + 3 );
 199 
 200 /* Create an AST Plot. This is a special sort of FrameSet in which the
 201    base Frame corresponds to graphics coordinates. All the coordinate
 202    Frames and Mappings read from the FITS-WCS headers are added into the
 203    Plot so that graphics can be drawn in any coordinate system. The extent 
 204    of the FITS array in pixel coordinates is mapped onto the extent of the 
 205    graphics device as returned above by cpgqwin. The AST library comes with
 206    a driver module that provides primitive drawing functions by calling
 207    appropriate PGFPLOT functions. It is simple to write driver modules
 208    for other graphics systems such as Tcl/Tk, Java/Swing, etc. Set a few
 209    graphics attributes to show the sort of thing that can be done. */
 210          plot = astPlot( fset, gbox, bbox, "Colour(border)=2,Colour(ticks)=2,"
 211                                        "Colour(axes)=2,Grid=1,Colour(grid)=3,"
 212                                        "Style(grid)=4" );
 213 
 214 /* Draw a set of annotated coordinate axes labelling the FITS WCS axes. */
 215          astGrid( plot );
 216 
 217 /* If there is any overlap (or if the overlap test could not be performed), 
 218    add the STC-S Region into the Plot. We use a UnitMap to connect it to 
 219    the current Frame (the FITS WCS frame). */
 220          if( 1 || overlap_flag == 2 || overlap_flag == 4 || !wcsbox ) {
 221             astAddFrame( plot, AST__CURRENT, astUnitMap( 2, " " ), wcsreg );
 222 
 223 /* Now draw the border round the STC-S Region. A Region is a sub-class of
 224    Mapping and so can be used to transform positions. When a Region is used 
 225    as a Mapping, positions that are inside the Region are left unchanged
 226    by the transformation, and positions that are outside the Region are
 227    transformed into "bad" positions (i.e. every axis value has the nmagic
 228    value AST__BAD indicating that the axis value is undefined). The
 229    astBorder method is a generic function that will outline the areas
 230    within the current coordinate Frame of the Plot that correspond to 
 231    valid (i.e. non-bad) positions. Set the colour and line thickness first
 232    to emphasise the border. */
 233             astSet( plot, "Colour(border)=4,Width(border)=8" );
 234             (void) astBorder( plot );
 235          }
 236 
 237 /* Set the returned system status to indicate success. */
 238          status = 0;
 239 
 240 /* Close down PGPLOT. */
 241          cpgend();
 242       }
 243    }
 244 
 245 /* End the AST Object context. All Objects created since the
 246    corresponding invocation of astbegin will be annulled automatically. */
 247    astEnd;
 248 
 249 /* If an error occurred in the AST library, set the returned system
 250    status non-zero. */
 251    if( !astOK ) status = 1;
 252    return status;
 253 }
 254 
 255 
 256 
 257 
 258 
 259 
 260 
 261 /* -------------------------------------------------------------------
 262  * This is a function that reads a line of text from the disk file and
 263  * returns it to the AST library. It is called from within the astRead
 264  * function. 
 265 */
 266 
 267 const char *source( void ){
 268    static char buffer[ MAX_LINE_LEN + 2 ];
 269    FILE *fd = astChannelData;
 270    return fgets( buffer, MAX_LINE_LEN + 2, fd );
 271 }
 272 
 273 
 274 
 275 /* -------------------------------------------------------------------
 276  * This function reads a set of FITS-WCS headers from a given text file,
 277  * and attempts to convert them into an AST FrameSet. If successful, a
 278  * pointer to the FrameSet is returned. A NULL pointer is returned if
 279  * anything goes wrong, or if the WCS is not 2-dimensional. The values of 
 280  * the NAXIS1 and NAXIS2 headers are returned in "*naxis1" and "*naxis2". 
 281 */
 282 
 283 AstFrameSet *ReadFitsHeaders( const char *file, int *naxis1, int *naxis2 ){
 284    AstFitsChan *chan;
 285    AstFrameSet *result;
 286    AstObject *object;
 287    FILE *fd;
 288 
 289 /* Initialise the returned pointer to indicate that no FrameSet has yet
 290    been read. */
 291    result = NULL;
 292 
 293 /* Attempt to open the FITS header file */
 294    fd = fopen( file, "r" );
 295    if( !fd ) {
 296       printf("Failed to open FITS header file '%s'.\n", file );
 297 
 298 /* If successful, create a FitsChan. This is the object that converts 
 299    external FITS headers into corresponding AST Objects. Tell it to use 
 300    the "source" function for obtaining lines of text from the disk file. */
 301    } else {
 302       chan = astFitsChan( source, NULL, " " );
 303 
 304 /* Associate the descriptor for the input disk file with the StcsChan.
 305    This makes it available to the "source" function. Since this
 306    application is single threaded, we could instead have made "fd" a 
 307    global variable, but the ChannelData facility is used here to illustrate 
 308    how to pass data to a source or sink function safely in a multi-threaded 
 309    application. */
 310       astPutChannelData( chan, fd );
 311 
 312 /* Attempt to read the FITS heades and convert them into an AST FrameSet. */
 313       object = astRead( chan );
 314       
 315 /* The astRead function is a generic function and so returns a generic
 316    AstObject pointer. Check an Object was created successfully. */
 317       if( !object ) {
 318          printf( "Failed to read an AST Object from FITS header file '%s'.\n", 
 319                  file );
 320 
 321 /* Now check that the object read is actually an AST FrameSet, rather than
 322    some other class of AST Object. */
 323       } else if( !astIsAFrameSet( object ) ) {      
 324          printf( "Expected a FrameSet but read a %s from FITS header "
 325                  "file '%s'.\n", astGetC( object, "Class" ), file );
 326 
 327 /* If the Object is a FrameSet, return the FrameSet pointer. */
 328       } else {
 329          result = (AstFrameSet *) object;
 330 
 331 /* Check the WCS is 2-dimensional. If not, report an error and set the
 332    returned pointer to NULL. The memory used to store the FrameSet will
 333    be released when the current AST object context is ended (by calling
 334    astEnd). */
 335         if( astGetI( result, "Naxes" ) != 2 ) {
 336            printf( "The FITS WCS is not 2-dimensional.\n");
 337            result = NULL;
 338 
 339 /* If it is 2-dimensional, get the NAXIS1 and NAXIS2 keyword values. */
 340         } else {
 341 
 342             if( !astGetFitsI( chan, "NAXIS1", naxis1 ) ){
 343                printf("Keyword 'NAXIS1' not found in header\n" );
 344                result = NULL;
 345             } 
 346 
 347             if( !astGetFitsI( chan, "NAXIS2", naxis2 ) ){
 348                printf("Keyword 'NAXIS2' not found in header\n" );
 349                result = NULL;
 350             } 
 351        
 352          }
 353       }
 354 
 355 /* Close the file. */
 356       fclose( fd );
 357    }
 358 
 359    return result;
 360 }
 361 
 362 
 363 /* -------------------------------------------------------------------
 364  * This function reads an STC-S description from a given text file,
 365  * and attempts to convert them into an AST Region. If successful, a
 366  * pointer to the Region is returned. A NULL pointer is returned if
 367  * anything goes wrong. 
 368 */
 369 
 370 AstRegion *ReadStcs( const char *file ){
 371    AstStcsChan *chan;
 372    AstRegion *result;
 373    AstObject *object;
 374    FILE *fd;
 375 
 376 /* Initialise the returned pointer to indicate that no Region has yet
 377    been read. */
 378    result = NULL;
 379 
 380 /* Attempt to open the STC-S file */
 381    fd = fopen( file, "r" );
 382    if( !fd ) {
 383       printf("Failed to open STC-S descrption file '%s'.\n", file );
 384 
 385 /* If successful, create an StcsChan. This is the object that converts 
 386    external STC-S descriptions into corresponding AST Objects. Tell it to 
 387    use the "source" function for obtaining lines of text from the disk 
 388    file. */
 389    } else {
 390       chan = astStcsChan( source, NULL, " " );
 391 
 392 /* Associate the descriptor for the input disk file with the StcsChan.
 393    This makes it available to the "source" function. Since this
 394    application is single threaded, we could instead have made "fd" a 
 395    global variable, but the ChannelData facility is used here to illustrate 
 396    how to pass data to a source or sink function safely in a multi-threaded 
 397    application. */
 398       astPutChannelData( chan, fd );
 399 
 400 /* The default behaviour of the astRead function when used on an StcsChan is 
 401    to read and return the AstroCoordArea as an AST Region. This behaviour
 402    can be changed by assigning appropriate values to the StcsChan attributes
 403    "StcsArea", "StcsCoords" and "StcsProps". Options exist to return the 
 404    AstroCoords as an AST PointList, and/or to return the individual
 405    property values read from the STC-S text in the form of an AST KeyMap
 406    (a sort of hashmap). For now, just take the default action of reading the 
 407    AstroCoordsArea. */
 408       object = astRead( chan );
 409       
 410 /* The astRead function is a generic function and so returns a generic
 411    AstObject pointer. Check an Object was created successfully. */
 412       if( !object ) {
 413          printf( "Failed to read an AST Object from STC-S description "
 414                  "file '%s'.\n", file );
 415 
 416 /* Now check that the object read is actually an AST Region, rather than
 417    some other class of AST Object. */
 418       } else if( !astIsARegion( object ) ) {      
 419          printf( "Expected a Region but read a %s from STC-S description "
 420                  "file '%s'.\n", astGetC( object, "Class" ), file );
 421 
 422 /* If the Object is a Region, return the Region pointer. */
 423       } else {
 424          result = (AstRegion *) object;
 425       }
 426 
 427 /* Close the file. */
 428       fclose( fd );
 429    }
 430 
 431    return result;
 432 }
 433 
 434 
 435 
 436 

Attached Files

To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.

You are not allowed to attach a file to this page.