path: root/ast/stcschan-demo3.c

/* Name:
      stcschan-demo3.c

   Purpose:
      A demonstration of the facilities provided by the AST library
      for reading STC metadata encoded using the STC-S linear string
      format.

   Description:
      This program reads an STC-S description from a text file, and also
      reads a set of 2-D spatial FITS-WCS headers from another (text) file.
      It then opens a specified graphics device, and displays an annotated
      coordinate grid covering the region described by the FITS headers.
      Finally, it draws the outline of the spatial extent of the STC-S
      description over the top of the annotated coordinate grid.

   Usage:
      % stcschan-demo3 <stcs-file> <header-file> <device>

      <stcs-file>: The path to a text file containing the STC-S
      description.

      <header-file>: The path to a text file containing a set of FITS-WCS
      headers.

      <device>: The name of an available PGPLOT graphics device. If not
      supplied, the available device names will be listed and the program
      will then exit.

   Example:
      % stcschan-demo3 m31.stcs andromeda.head /xserve

   To compile and link:
      Assuming your starlink distribution is in "/star":

      % gcc -o stcschan-demo3 stcschan-demo3.c -g -L/star/lib \
             -I/star/include `ast_link -pgplot` -lcpgplot

*/

/* Include system headers. */
#include <stdio.h>
#include <string.h>

/* Include the PGPLOT header. */
#include "cpgplot.h"

/* Include the AST library header. */
#include "ast.h"

/* Maximum allowed length for a single line of text from the disk file. */
#define MAX_LINE_LEN 100

/* Prototypes */
const char *source( void );
AstFrameSet *ReadFitsHeaders( const char *, int *, int * );
AstRegion *ReadStcs( const char * );


int main( int argc, char **argv ){

/* Local variables: */
   AstBox *pixbox;
   AstFrame *pixfrm;
   AstFrameSet *align_fs = NULL;
   AstFrameSet *fset = NULL;
   AstPlot *plot;
   AstRegion *reg = NULL;
   AstRegion *wcsbox = NULL;
   AstRegion *wcsreg = NULL;
   const char *dev;
   double bbox[ 4 ];
   float gbox[ 4 ];
   int overlap_flag, status, naxis1, naxis2, ibase;

/* Initialise the returned system status to indicate failure. */
   status = 1;

/* Start an AST object context. This means we do not need to annull
   each AST Object individually. Instead, all Objects created within
   this context will be annulled automatically by the corresponding
   invocation of astEnd. */
   astBegin;

/* Check there are enough command line arguments. */
   if( argc < 3 ) {
      printf( "Usage: stcschan-demo3 <stcs-file> <header-file> <device>\n" );

/* If so, attempt to read the STC-S description, creating a corresponding
   AST Region. If this is successful, attempt to read the FITS header file
   and create an equivalent AST FrameSet. */
   } else {
      reg = ReadStcs( argv[ 1 ] );
      if( reg ) fset = ReadFitsHeaders( argv[ 2 ], &naxis1, &naxis2 );
   }

/* Check we obtained a Region and a FrameSet successfully. */
   if( reg && fset ){

/* Check that we can align the FITS WCS with the spatial axes in the
   STC-S description. AST contains various built-in conversions between
   standard celestial coordinate system. The necessary conversion will be
   identified and used automatically if necessary to achieve alignment.
   The returned object "align_fs" is a FrameSet that encapsulates the
   the FITS WCS coordinate system, the STC-S spatial coordinate system and
   the Mapping between them. A NULL pointer is returned if alignment is
   not possible. Note, the astConvert method changes the base Frame in
   any supplied FrameSet to indicate which coordinate frame was used for
   alignment. So we first note the original base Frame index and then
   re-instate it afterwards. */
      ibase = astGetI( fset, "Base" );
      align_fs = astConvert( reg, fset, " " );
      astSetI( fset, "Base", ibase );

      if( !align_fs ) {
         printf( "Could not align the FITS WCS with the spatial axes "
                 "in the STC-S\n" );

/* If alignment was possible, use the alignment FrameSet to create a new
   Region representing the same region as the supplied STC-S, but
   expressed in the coordinate system of the FITS WCS. The FrameSet class
   inherits from both Frame and Mapping, and so the "align_fs" FrameSet can
   be used both as the Mapping in this call, and as the Frame. When used
   as a Mapping, a FrameSet represents the transformation between its base
   and current Frame. When used as a Frame, a FrameSet represents its
   current Frame. */
      } else {
         wcsreg = astMapRegion( reg, align_fs, align_fs );

/* It would be nice to warn the user if the STC-S AstroCoordsArea does
   not overlap the FITS grid. To do so we need a Region representing the
   FITS grid. Create one now (an AST Box). First store the bounds of the
   FITS grid, in pixel coordinates (i.e. a system in which the centre of
   the bottom left pixel is at (1.0,1.0) ). */
         bbox[ 0 ] = 0.5;
         bbox[ 1 ] = 0.5;
         bbox[ 2 ] = (double) naxis1 + 0.5;
         bbox[ 3 ] = (double) naxis2 + 0.5;

/* Get a pointer to the Frame describing the FITS pixel coordinate system
   (the base Frame in the FITS FrameSet). */
         pixfrm = astGetFrame( fset, AST__BASE );

/* Create a Box that encompasses the required range of axis values within
   the pixel coordinate Frame. */
         pixbox = astBox( pixfrm, 1, bbox, bbox + 2, NULL, " " );

/* Create another Region that represents the same area, but in the FITS
   WCS. */
         wcsbox = astMapRegion( pixbox, fset, fset );

/* If the previous step failed, it probably means the FITS header covers
   the entire sky, resulting the corners of the pixel grid having invalid
   sky positions. So cancel the error and omit the overlap test. */
         if( !wcsbox ) {
            astClearStatus;
            printf("\nContinuing, but omitting overlap test...\n\n");

/* Now see if the Region representing the FITS grid overlaps the region
   read from the STC-S description.*/
         } else {
            overlap_flag = astOverlap( wcsreg, wcsbox );

            if( overlap_flag == 1 || overlap_flag == 6 ) {
               printf( "\nThere is no overlap between the FITS grid and the "
                       "STC-S AstroCoordsArea\n\n" );

            } else if( overlap_flag == 3 || overlap_flag == 5 ) {
               printf( "\nThe FITS grid is completely contained within the "
                       "STC-S AstroCoordsArea\n\n" );
            }

         }
      }
   }

/* Check we obtained a FITS-WCS to STC-S Mapping, and that no error has
   occurred in AST. */
   if( wcsreg && astOK ){

/* Open PGPLOT using the specified device. Prompt the user for a device if
   none was supplied on the command line. */
      if( argc < 4 ) {
         dev = "?";
      } else {
         dev = argv[ 3 ];
      }
      if( cpgbeg( 0, dev, 1, 1 ) == 1 ) {

/* Clear the screen. */
         cpgpage();

/* Ensure the graphics window has equal scales on both axes. */
         cpgwnad( 0.0f, 1.0f, 0.0f, 1.0f );

/* Find the extent of the graphics window, and store in an array suitable
   for passing to the astPlot function. */
         cpgqwin( gbox, gbox + 2, gbox + 1, gbox + 3 );

/* Create an AST Plot. This is a special sort of FrameSet in which the
   base Frame corresponds to graphics coordinates. All the coordinate
   Frames and Mappings read from the FITS-WCS headers are added into the
   Plot so that graphics can be drawn in any coordinate system. The extent
   of the FITS array in pixel coordinates is mapped onto the extent of the
   graphics device as returned above by cpgqwin. The AST library comes with
   a driver module that provides primitive drawing functions by calling
   appropriate PGFPLOT functions. It is simple to write driver modules
   for other graphics systems such as Tcl/Tk, Java/Swing, etc. Set a few
   graphics attributes to show the sort of thing that can be done. */
         plot = astPlot( fset, gbox, bbox, "Colour(border)=2,Colour(ticks)=2,"
                                       "Colour(axes)=2,Grid=1,Colour(grid)=3,"
                                       "Style(grid)=4" );

/* Draw a set of annotated coordinate axes labelling the FITS WCS axes. */
         astGrid( plot );

/* If there is any overlap (or if the overlap test could not be performed),
   add the STC-S Region into the Plot. We use a UnitMap to connect it to
   the current Frame (the FITS WCS frame). */
         if( 1 || overlap_flag == 2 || overlap_flag == 4 || !wcsbox ) {
            astAddFrame( plot, AST__CURRENT, astUnitMap( 2, " " ), wcsreg );

/* Now draw the border round the STC-S Region. A Region is a sub-class of
   Mapping and so can be used to transform positions. When a Region is used
   as a Mapping, positions that are inside the Region are left unchanged
   by the transformation, and positions that are outside the Region are
   transformed into "bad" positions (i.e. every axis value has the nmagic
   value AST__BAD indicating that the axis value is undefined). The
   astBorder method is a generic function that will outline the areas
   within the current coordinate Frame of the Plot that correspond to
   valid (i.e. non-bad) positions. Set the colour and line thickness first
   to emphasise the border. */
            astSet( plot, "Colour(border)=4,Width(border)=8" );
            (void) astBorder( plot );
         }

/* Set the returned system status to indicate success. */
         status = 0;

/* Close down PGPLOT. */
         cpgend();
      }
   }

/* End the AST Object context. All Objects created since the
   corresponding invocation of astbegin will be annulled automatically. */
   astEnd;

/* If an error occurred in the AST library, set the returned system
   status non-zero. */
   if( !astOK ) status = 1;
   return status;
}







/* -------------------------------------------------------------------
 * This is a function that reads a line of text from the disk file and
 * returns it to the AST library. It is called from within the astRead
 * function.
*/

const char *source( void ){
   static char buffer[ MAX_LINE_LEN + 2 ];
   FILE *fd = astChannelData;
   return fgets( buffer, MAX_LINE_LEN + 2, fd );
}



/* -------------------------------------------------------------------
 * This function reads a set of FITS-WCS headers from a given text file,
 * and attempts to convert them into an AST FrameSet. If successful, a
 * pointer to the FrameSet is returned. A NULL pointer is returned if
 * anything goes wrong, or if the WCS is not 2-dimensional. The values of
 * the NAXIS1 and NAXIS2 headers are returned in "*naxis1" and "*naxis2".
*/

AstFrameSet *ReadFitsHeaders( const char *file, int *naxis1, int *naxis2 ){
   AstFitsChan *chan;
   AstFrameSet *result;
   AstObject *object;
   FILE *fd;

/* Initialise the returned pointer to indicate that no FrameSet has yet
   been read. */
   result = NULL;

/* Attempt to open the FITS header file */
   fd = fopen( file, "r" );
   if( !fd ) {
      printf("Failed to open FITS header file '%s'.\n", file );

/* If successful, create a FitsChan. This is the object that converts
   external FITS headers into corresponding AST Objects. Tell it to use
   the "source" function for obtaining lines of text from the disk file. */
   } else {
      chan = astFitsChan( source, NULL, " " );

/* Associate the descriptor for the input disk file with the StcsChan.
   This makes it available to the "source" function. Since this
   application is single threaded, we could instead have made "fd" a
   global variable, but the ChannelData facility is used here to illustrate
   how to pass data to a source or sink function safely in a multi-threaded
   application. */
      astPutChannelData( chan, fd );

/* Attempt to read the FITS heades and convert them into an AST FrameSet. */
      object = astRead( chan );

/* The astRead function is a generic function and so returns a generic
   AstObject pointer. Check an Object was created successfully. */
      if( !object ) {
         printf( "Failed to read an AST Object from FITS header file '%s'.\n",
                 file );

/* Now check that the object read is actually an AST FrameSet, rather than
   some other class of AST Object. */
      } else if( !astIsAFrameSet( object ) ) {
         printf( "Expected a FrameSet but read a %s from FITS header "
                 "file '%s'.\n", astGetC( object, "Class" ), file );

/* If the Object is a FrameSet, return the FrameSet pointer. */
      } else {
         result = (AstFrameSet *) object;

/* Check the WCS is 2-dimensional. If not, report an error and set the
   returned pointer to NULL. The memory used to store the FrameSet will
   be released when the current AST object context is ended (by calling
   astEnd). */
        if( astGetI( result, "Naxes" ) != 2 ) {
           printf( "The FITS WCS is not 2-dimensional.\n");
           result = NULL;

/* If it is 2-dimensional, get the NAXIS1 and NAXIS2 keyword values. */
        } else {

            if( !astGetFitsI( chan, "NAXIS1", naxis1 ) ){
               printf("Keyword 'NAXIS1' not found in header\n" );
               result = NULL;
            }

            if( !astGetFitsI( chan, "NAXIS2", naxis2 ) ){
               printf("Keyword 'NAXIS2' not found in header\n" );
               result = NULL;
            }

         }
      }

/* Close the file. */
      fclose( fd );
   }

   return result;
}


/* -------------------------------------------------------------------
 * This function reads an STC-S description from a given text file,
 * and attempts to convert them into an AST Region. If successful, a
 * pointer to the Region is returned. A NULL pointer is returned if
 * anything goes wrong.
*/

AstRegion *ReadStcs( const char *file ){
   AstStcsChan *chan;
   AstRegion *result;
   AstObject *object;
   FILE *fd;

/* Initialise the returned pointer to indicate that no Region has yet
   been read. */
   result = NULL;

/* Attempt to open the STC-S file */
   fd = fopen( file, "r" );
   if( !fd ) {
      printf("Failed to open STC-S descrption file '%s'.\n", file );

/* If successful, create an StcsChan. This is the object that converts
   external STC-S descriptions into corresponding AST Objects. Tell it to
   use the "source" function for obtaining lines of text from the disk
   file. */
   } else {
      chan = astStcsChan( source, NULL, " " );

/* Associate the descriptor for the input disk file with the StcsChan.
   This makes it available to the "source" function. Since this
   application is single threaded, we could instead have made "fd" a
   global variable, but the ChannelData facility is used here to illustrate
   how to pass data to a source or sink function safely in a multi-threaded
   application. */
      astPutChannelData( chan, fd );

/* The default behaviour of the astRead function when used on an StcsChan is
   to read and return the AstroCoordArea as an AST Region. This behaviour
   can be changed by assigning appropriate values to the StcsChan attributes
   "StcsArea", "StcsCoords" and "StcsProps". Options exist to return the
   AstroCoords as an AST PointList, and/or to return the individual
   property values read from the STC-S text in the form of an AST KeyMap
   (a sort of hashmap). For now, just take the default action of reading the
   AstroCoordsArea. */
      object = astRead( chan );

/* The astRead function is a generic function and so returns a generic
   AstObject pointer. Check an Object was created successfully. */
      if( !object ) {
         printf( "Failed to read an AST Object from STC-S description "
                 "file '%s'.\n", file );

/* Now check that the object read is actually an AST Region, rather than
   some other class of AST Object. */
      } else if( !astIsARegion( object ) ) {
         printf( "Expected a Region but read a %s from STC-S description "
                 "file '%s'.\n", astGetC( object, "Class" ), file );

/* If the Object is a Region, return the Region pointer. */
      } else {
         result = (AstRegion *) object;
      }

/* Close the file. */
      fclose( fd );
   }

   return result;
}