New Region method astPointInRegion to test a single point

A convenience wrapper that just uses astTranN to transform the supplied
point and then tests the resulting axis values.

ast.news

@@ -1,6 +1,6 @@
 AST Library
 -----------
-   A new release (V9.0.3) of the Starlink AST (astrometry) library is
+   A new release (V9.1.0) of the Starlink AST (astrometry) library is
 now available.
 
    AST provides a comprehensive range of facilities for attaching
@@ -14,7 +14,6 @@ celestial coordinate systems. AST is portable and
 environment-independent.
 
 
-
 Main Changes in this Version
 ----------------------------
 
@@ -24,6 +23,10 @@ source files into a subdirectory named "src".
 - A bug has been fixed in the TimeFrame class that caused the time
 returned by astCurrentTime to be wrong by 37 seconds.
 
+- The Region class has a new convenience method (astPointInRegion) to
+test a if a single point is inside a Region.
+
+
 Main Changes in V9.0.2
 ----------------------
 
@@ -32,6 +35,7 @@ Main Changes in V9.0.2
 in conjunction with the TabOK attribute to force the use of the
 "-TAB" algorithm when writing a FrameSet out using the FITS-WCS encoding.
 
+
 Main Changes in V9.0.0
 ----------------------
 

ast_par.source

@@ -504,6 +504,7 @@ c      PARAMETER ( AST__UKERNN = 2 )
       INTEGER*8 AST_MASK8US
       INTEGER*8 AST_MASK8UW
       INTEGER*8 AST_MASK8W
+      LOGICAL AST_POINTINREGION
 
 *  Box class.
       INTEGER AST_BOX

ast_tester/testregions.f

@@ -193,7 +193,7 @@ subroutine checkInterval( status )
       integer status, frm1, frm2, frm3, unc, int1, int2, int3, int4,
      :        int5, frm4, map, outperm(6), inperm(6), pm, reg
       double precision lbnd(3), ubnd(3), p(5,3), q(5,3),in(4,3),out(4,3)
-      double precision xin(9), yin(9), xout(9), yout(9)
+      double precision xin(9), yin(9), xout(9), yout(9), point(3)
 
       logical hasframeset
 
@@ -261,6 +261,20 @@ subroutine checkInterval( status )
       if( q(5,2) .ne. AST__BAD ) call stopit( status, 'Interval 5b' )
       if( q(5,3) .ne. AST__BAD ) call stopit( status, 'Interval 5c' )
 
+      point(1) = p(1,1)
+      point(2) = p(1,2)
+      point(3) = p(1,3)
+      if( .not. ast_pointinregion( int1, point, status ) ) then
+         call stopit( status, 'Interval 5d' )
+      end if
+
+      point(1) = p(5,1)
+      point(2) = p(5,2)
+      point(3) = p(5,3)
+      if( ast_pointinregion( int1, point, status ) ) then
+         call stopit( status, 'Interval 5e' )
+      end if
+
       call ast_negate( int1, status )
       call ast_trann( int1, 5, 3, 5, p, .true., 3, 5, q, status )
       if( q(1,1) .ne. p(1,1)) call stopit( status, 'Interval 6' )

component.xml

@@ -2,7 +2,7 @@
 <!DOCTYPE component SYSTEM "componentinfo.dtd">
 
 <component id="ast" support="S">
-  <version>9.0.3</version>
+  <version>9.1.0</version>
   <path>libext/ast</path>
   <description>WCS library</description>
   <abstract>

configure.ac

@@ -8,8 +8,8 @@ dnl   automake 1.6 or better.
 
 dnl   Initialisation: package name and version number
 m4_define([v_maj], [9])
-m4_define([v_min], [0])
-m4_define([v_mic], [3])
+m4_define([v_min], [1])
+m4_define([v_mic], [0])
 m4_define([project_version], [v_maj.v_min.v_mic])
 AC_INIT([ast],[project_version],[[email protected]])
 AC_CONFIG_AUX_DIR([build-aux])

src/fregion.c

@@ -21,6 +21,7 @@
 *     AST_GETREGIONFRAME
 *     AST_GETREGIONFRAMESET
 *     AST_OVERLAP
+*     AST_POINTINREGION
 *     AST_SETUNC
 *     AST_GETUNC
 *     AST_SHOWMESH
@@ -175,6 +176,20 @@ F77_INTEGER_FUNCTION(ast_overlap)( INTEGER(THIS),
    return RESULT;
 }
 
+F77_LOGICAL_FUNCTION(ast_pointinregion)( INTEGER(THIS),
+                                         DOUBLE_ARRAY(POINT),
+                                         INTEGER(STATUS) ) {
+   GENPTR_INTEGER(THIS)
+   GENPTR_DOUBLE(POINT)
+   F77_LOGICAL_TYPE(RESULT);
+
+   astAt( "AST_POINTINREGION", NULL, 0 );
+   astWatchSTATUS(
+      RESULT = astPointInRegion( astI2P( *THIS ), POINT ) ? F77_TRUE : F77_FALSE;
+   )
+   return RESULT;
+}
+
 /* AST_MASK<X> requires a function for each possible data type, so
    define it via a macro. */
 #define MAKE_AST_MASK(f,F,Ftype,X,Xtype) \

src/region.c

@@ -37,7 +37,10 @@ f     AST_FORMAT routine
 c     astTran<X> functions
 f     AST_TRAN<X> routines
 *     is the way to determine if a given position is inside or outside the
-*     Region. When used as a Mapping, most classes of Frame are equivalent to
+*     Region (but see also the convenience function
+c     astPointInRegion).
+f     AST_POINTINREGION).
+*     When used as a Mapping, most classes of Frame are equivalent to
 *     a UnitMap. However, the Region class modifies this behaviour so that a
 *     Region acts like a UnitMap only for input positions which are within the
 *     area represented by the Region. Input positions which are outside the
@@ -127,6 +130,8 @@ c     - astNegate: Toggle the value of the Negated attribute
 f     - AST_NEGATE: Toggle the value of the Negated attribute
 c     - astOverlap: Determines the nature of the overlap between two Regions
 f     - AST_OVERLAP: Determines the nature of the overlap between two Regions
+c     - astPointInRegion: Test if a single point is inside a Region
+f     - AST_POINTINREGION: Test if a single point is inside a Region
 c     - astMask<X>: Mask a region of a data grid
 f     - AST_MASK<X>: Mask a region of a data grid
 c     - astSetUnc: Associate a new uncertainty with a Region
@@ -138,6 +143,7 @@ f     - AST_SHOWMESH: Display a mesh of points on the surface of a Region
 *     Copyright (C) 1997-2006 Council for the Central Laboratory of the
 *     Research Councils
 *     Copyright (C) 2009 Science & Technology Facilities Council.
+*     Copyright (C) 2020 East Asian Observatory
 *     All Rights Reserved.
 
 *  Licence:
@@ -254,6 +260,8 @@ f     - AST_SHOWMESH: Display a mesh of points on the surface of a Region
 *        regions defiend in cartesian spaces).
 *     24-SEP-2019 (DSB):
 *        Added 8-byte interface for astMask<X>.
+*     23-JAN-2020 (DSB):
+*        Added astPointInRegion.
 *class--
 
 *  Implementation Notes:
@@ -964,6 +972,7 @@ static int LineCrossing( AstFrame *, AstLineDef *, AstLineDef *, double[5], int
 static int Match( AstFrame *, AstFrame *, int, int **, int **, AstMapping **, AstFrame **, int * );
 static int Overlap( AstRegion *, AstRegion *, int * );
 static int OverlapX( AstRegion *, AstRegion *, int * );
+static int PointInRegion( AstRegion *, const double *, int * );
 static int RegDummyFS( AstRegion *, int * );
 static int RegPins( AstRegion *, AstPointSet *, AstRegion *, int **, int * );
 static int SubFrame( AstFrame *, AstFrame *, int, const int *, const int *, AstMapping **, AstFrame **, int * );
@@ -4624,6 +4633,7 @@ void astInitRegionVtab_(  AstRegionVtab *vtab, const char *name, int *status ) {
    vtab->GetRegionMesh = GetRegionMesh;
    vtab->GetRegionPoints = GetRegionPoints;
    vtab->GetRegionDisc = GetRegionDisc;
+   vtab->PointInRegion = PointInRegion;
    vtab->RegOverlay = RegOverlay;
    vtab->RegFrame = RegFrame;
    vtab->RegDummyFS = RegDummyFS;
@@ -7632,6 +7642,111 @@ static AstFrame *PickAxes( AstFrame *this_frame, int naxes, const int axes[],
    return frame;
 }
 
+static int PointInRegion( AstRegion *this, const double *point, int *status ){
+/*
+*++
+*  Name:
+c     astPointInRegion
+f     AST_POINTINREGION
+
+*  Purpose:
+*     Tests if a single point is inside a Region.
+
+*  Type:
+*     Public virtual function.
+
+*  Synopsis:
+c     #include "region.h"
+c     int astPointInRegion( AstRegion *this, const double point[] )
+f     RESULT = AST_POINTINREGION( THIS, POINT, STATUS )
+
+*  Class Membership:
+*     Region method.
+
+*  Description:
+c     This function
+f     This routine
+*     returns a logical value indicating if a supplied point is inside a
+*     supplied Region.
+
+*  Parameters:
+c     this
+f     THIS = INTEGER (Given)
+*        Pointer to the Region.
+c     point
+f     POINT( * ) = DOUBLE PRECISION (Given)
+c        Pointer to an
+f        An
+*        array holding the axis values of the point to be tested. The
+*        number of values in this array should match the number of axes in
+*        the supplied Region.
+f     STATUS = INTEGER (Given and Returned)
+f        The global status.
+
+*  Returned Value:
+c     astPointInRegion()
+f     AST_POINTINREGION = LOGICAL
+c        Zero if the point is outside the Region, and non-zero otherwise.
+f        .FALSE. if the point is outside the Region, and .TRUE. otherwise.
+
+*  Notes:
+*    - If the supplied value is AST__BAD on every axis, then
+c    0
+f    .FALSE.
+*    is always returned.
+*    - If many points need to be tested, then it is usually more
+*    efficient to use the Region as a Mapping to transform all the
+*    points. This can be done using one of the Mapping transformation
+*    methods
+c    (astTran<X>).
+f    (AST_TRAN<X>).
+*    If a transformed axis value is AST__BAD then the corresponding
+*    input point is outside the Region.
+c    - Zero will be returned if an error occurs.
+f    - .FALSE. will be returned if an error occurs.
+
+*--
+*/
+
+/* Local Variables: */
+   double *out;
+   int i;
+   int naxes;
+   int result;
+
+/* Initialise the radius. */
+   result = 0;
+
+/* Check the inherited status. */
+   if( !astOK ) return result;
+
+/* Get the number of axes in the Region. */
+   naxes = astGetNaxes( this );
+
+/* Allocate an array in which to store the transformed axis values. */
+   out = astMalloc( naxes*sizeof(*out) );
+
+/* Use the Region to transform the supplied position. */
+   astTranN( this, 1, naxes, 1, point, 1, naxes, 1, out );
+
+/* If any good output axis value has been set to AST__BAD, then the point is
+   outside the Region. */
+   if( astOK ) {
+      for( i = 0; i < naxes; i++ ) {
+         if( point[ i ] != AST__BAD ) {
+            if( out[ i ] != AST__BAD ) result = 1;
+            break;
+         }
+      }
+   }
+
+/* Free resources. */
+   out = astFree( out );
+
+/* Return the result */
+   return result;
+}
+
 static void RegBaseBox( AstRegion *this, double *lbnd, double *ubnd, int *status ){
 /*
 *+
@@ -9418,8 +9533,8 @@ f     AST_GETREGIONDisc
 
 *  Synopsis:
 c     #include "region.h"
-c     void GetRegionDisc( AstRegion *this, double centre[3],
-c                         double *radius )
+c     void astGetRegionDisc( AstRegion *this, double centre[3],
+c                            double *radius )
 f     CALL AST_GETREGIONDISC( THIS, CENTRE, RADIUS, STATUS )
 
 *  Class Membership:
@@ -13545,6 +13660,10 @@ double *astRegCentre_( AstRegion *this, double *cen, double **ptr, int index,
    if ( !astOK ) return NULL;
    return (**astMEMBER(this,Region,RegCentre))( this, cen, ptr, index, ifrm, status );
 }
+int astPointInRegion_( AstRegion *this, const double *point, int *status ){
+   if ( !astOK ) return 0;
+   return (**astMEMBER(this,Region,PointInRegion))( this, point, status );
+}
 AstRegion *astGetNegation_( AstRegion *this, int *status ){
    if ( !astOK ) return NULL;
    return (**astMEMBER(this,Region,GetNegation))( this, status );

src/region.h

@@ -120,6 +120,7 @@ typedef struct AstRegionVtab {
 /* Properties (e.g. methods) specific to this class. */
    int (* Overlap)( AstRegion *, AstRegion *, int * );
    int (* OverlapX)( AstRegion *, AstRegion *, int * );
+   int (* PointInRegion)( AstRegion *, const double *, int * );
    AstRegion *(* MapRegion)( AstRegion *, AstMapping *, AstFrame *, int * );
    AstFrame *(* GetRegionFrame)( AstRegion *, int * );
    AstFrameSet *(* GetRegionFrameSet)( AstRegion *, int * );
@@ -255,6 +256,7 @@ AstFrame *astGetRegionFrame_( AstRegion *, int * );
 AstFrameSet *astGetRegionFrameSet_( AstRegion *, int * );
 int astOverlap_( AstRegion *, AstRegion *, int * );
 void astNegate_( AstRegion *, int * );
+int astPointInRegion_( AstRegion *, const double *, int * );
 
 #if HAVE_LONG_DOUBLE     /* Not normally implemented */
 int astMask4LD_( AstRegion *, AstMapping *, int, int, const int[], const int[], long double [], long double, int * );
@@ -476,6 +478,8 @@ astINVOKE(O,astGetRegionFrameSet_(astCheckRegion(this),STATUS_PTR))
 astINVOKE(V,astNegate_(astCheckRegion(this),STATUS_PTR))
 #define astOverlap(this,that) \
 astINVOKE(V,astOverlap_(astCheckRegion(this),astCheckRegion(that),STATUS_PTR))
+#define astPointInRegion(this,point) \
+astINVOKE(V,astPointInRegion_(astCheckRegion(this),point,STATUS_PTR))
 #define astSetUnc(this,unc) astINVOKE(V,astSetUnc_(astCheckRegion(this),unc?astCheckRegion(unc):NULL,STATUS_PTR))
 #define astGetUnc(this,def) astINVOKE(O,astGetUnc_(astCheckRegion(this),def,STATUS_PTR))
 #define astGetRegionBounds(this,lbnd,ubnd) astINVOKE(V,astGetRegionBounds_(astCheckRegion(this),lbnd,ubnd,STATUS_PTR))

sun_master.tex

@@ -18,7 +18,7 @@
 \stardocnumber      {210.29}
 f-
 \stardocauthors     {R.F. Warren-Smith \& D.S. Berry}
-\stardocdate        {25th Novemeber 2019}
+\stardocdate        {23rd January 2020}
 \stardoctitle     {AST\linebreak%
                                 A Library for Handling\linebreak%
                                 World Coordinate Systems\linebreak%
@@ -22013,9 +22013,9 @@ \subsection{Changes Introduced in V9.0.2}
 
 
 \subsection{\xlabel{changes}\xlabel{list_of_most_recent_changes}Changes
-Introduced in V9.0.3}
+Introduced in V9.1.0}
 The following describes the most significant changes which have
-occurred in the AST library between versions V9.0.2 and V9.0.3 (the
+occurred in the AST library between versions V9.0.2 and V9.1.0 (the
 current version):
 
 \begin{enumerate}
@@ -22025,6 +22025,9 @@ \subsection{\xlabel{changes}\xlabel{list_of_most_recent_changes}Changes
 \item A bug has been fixed in the TimeFrame class that caused the time
 returned by astCurrentTime to be wrong by 37 seconds.
 
+\item The Region class has a new convenience method (astPointInRegion) to
+test a if a single point is inside a Region.
+
 \end{enumerate}
 
 Programs which are statically linked will need to be re-linked in