Skip to content

Commit

Permalink
astMapRegion now removes any unused base Frame axes
Browse files Browse the repository at this point in the history 
If the astMapRegion method is used to map a Region into a new Frame that
has fewer axes than the original Region, and if the inverse
transformation of the supplied Mapping does not specify a value for the
missing axes, then those axes are removed entirely from the Region.
Previously they were retained, but supplied automatically with bad
values. This affects the number of mesh points per axes for such Regions,
and so affects the accuracy of overlap determination.
  • Loading branch information
David Berry committed Dec 1, 2016
1 parent b2a56c3 commit 21fc804
Show file tree
Hide file tree
Showing 4 changed files with 219 additions and 15 deletions.
8 changes: 8 additions & 0 deletions ast.news
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,14 @@ extra argument that gives the number of values supplied in the arguments
array. Note, any existing code that uses this method will need to be
changed.

- If the astMapRegion (AST_MAPREGION) method is used to map a Region into
a new Frame that has fewer axes than the original Region, and if the
inverse transformation of the supplied Mapping does not specify a value
for the missing axes, then those axes are removed entirely from the
Region. Previously they were retained, but supplied with bad values. This
affects the number of mesh points per axes for such Regions, and so
affects the accuracy of overlap determination.


Main Changes in V8.3.0
----------------------
Expand Down
5 changes: 3 additions & 2 deletions ast_tester/teststc.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -544,15 +544,16 @@ subroutine Example1( status )

if( abs( lbnd(2) + 2.42406841E-06 ) .gt. 0.0000001E-06 )
: call stopit( status, 'Error 11' )
if( abs( 86400.0D0*lbnd(3) + 5.0D-5 ) .gt. 0.1E-10 )

if( abs( 86400.0D0*lbnd(3) + 4.9662776D-5 ) .gt. 0.1E-10 )
: call stopit( status, 'Error 12' )
if( abs( lbnd(4) - 0.07 ) .gt. 0.0001 )
: call stopit( status, 'Error 13' )
if( abs( ubnd(1) - 2.42406841E-06 ) .gt. 0.0000001E-06 )
: call stopit( status, 'Error 14' )
if( abs( ubnd(2) - 2.42406841E-06 ) .gt. 0.0000001E-06 )
: call stopit( status, 'Error 15' )
if( abs( 86400.0D0*ubnd(3) - 5.0D-5 ) .gt. 0.1E-10 )
if( abs( 86400.0D0*ubnd(3) - 4.9662776D-5 ) .gt. 0.1E-10 )
: call stopit( status, 'Error 16' )
if( abs( ubnd(4) - 0.17 ) .gt. 0.0001 )
: call stopit( status, 'Error 17' )
Expand Down
205 changes: 193 additions & 12 deletions region.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -226,9 +226,12 @@ f - AST_SHOWMESH: Display a mesh of points on the surface of a Region
* - Use astAxNorm to fix a bug in astGetRegionBounds for cases where
* the Region cross a longitude=0 singularity.
* 11-NOV-2016 (DSB):
* When loading a Region, use the dimensionality of the base Frame
* of the FrameSet (rather than the current Frame as it used to be)
* When loading a Region, use the dimensionality of the base Frame
* of the FrameSet (rather than the current Frame as it used to be)
* to define the number of axes required in the PointSet.
* 1-DEC-2016 (DSB):
* Changed MapRegion to remove any unnecessary base frame axes in
* the returned Region.
*class--
* Implementation Notes:
Expand Down Expand Up @@ -5359,27 +5362,40 @@ static AstRegion *MapRegion( AstRegion *this, AstMapping *map0,
/* Local Variables: */
AstFrame *frame;
AstFrameSet *fs;
AstMapping *tmap;
AstMapping *usemap;
AstMapping *map;
AstPointSet *ps2;
AstPointSet *ps1;
AstPointSet *pst;
AstPointSet *ps2;
AstRegion *usethis;
AstRegion *result;
double **ptr1;
double **ptr2;
int *axflags;
int *inax;
int *keep;
int *outax;
int i;
int icurr;
int j;
int nax1;
int nax2;
int nkept;
int nnew;
int nold;
int np;
int ntotal;
int ok;

/* Initialise. */
/* Initialise returned value. */
result = NULL;

/* Check the global error status. */
if ( !astOK ) return result;

/* Initialise local variables */
axflags = NULL;

/* If a FrameSet was supplied for the Mapping, use the base->current
Mapping */
if( astIsAFrameSet( map0 ) ) {
Expand Down Expand Up @@ -5407,14 +5423,24 @@ static AstRegion *MapRegion( AstRegion *this, AstMapping *map0,
astGetClass( map ) );
}

/* It must not introduce any bad axis values. We can only perform this
test reliably if the supplied Region has not bad axis values. */

/* Get the number of axes in the supplied Region. */
nold = astGetNaxes( this );

/* Get the number of axes in the returned Region. */
nnew = astGetNaxes( frame );

/* The forward transformation must not introduce any bad axis values. We
can only perform this test reliably if the supplied Region has no bad
axis values. */
ps1 = this->points;
if( ps1 ) {
nax1 = astGetNcoord( ps1 );
np = astGetNpoint( ps1 );
ptr1 = astGetPoints( ps1 );
if( ptr1 ) {

/* Check the axis values defining the Region are good. */
ok = 1;
for( i = 0; i < nax1 && ok; i++ ){
for( j = 0; j < np; j++ ) {
Expand All @@ -5425,12 +5451,18 @@ static AstRegion *MapRegion( AstRegion *this, AstMapping *map0,
}
}
if( ok ) {

/* Transform the points defining the supplied Region into the current Frame
of the Region. */
pst = astRegTransform( this, ps1, 1, NULL, NULL );

/* The use the supplied Mapping to transfom them into the new Frame. */
ps2 = astTransform( map, pst, 1, NULL );
nax2 = astGetNcoord( ps2 );

/* Test if any of these axis values are bad. */
ptr2 = astGetPoints( ps2 );
if( ptr2 ) {
for( i = 0; i < nax2 && ok; i++ ){
for( i = 0; i < nnew && ok; i++ ){
for( j = 0; j < np; j++ ) {
if( ptr2[ i ][ j ] == AST__BAD ){
ok = 0;
Expand All @@ -5443,6 +5475,29 @@ static AstRegion *MapRegion( AstRegion *this, AstMapping *map0,
"results from using the supplied %s to transform "
"the supplied %s is undefined.", status, astGetClass( this ),
astGetClass( map ), astGetClass( this ) );

/* If all axis values are good, use the inverse transformation of the
supplied Mapping to transform them back to the Frame of the supplied
Region. */
} else {
pst = astAnnul( pst );
pst = astTransform( map, ps2, 0, NULL );

/* Get a flag for each input of the supplied Mapping (i.e. each axis of
the supplied Region) which is non-zero if the inverse transformation
of the Mapping has produced any good axis values. */
ptr1 = astGetPoints( pst );
axflags = astCalloc( nold, sizeof(int) );
if( astOK ) {
for( i = 0; i < nold; i++ ){
for( j = 0; j < np; j++ ) {
if( ptr1[ i ][ j ] != AST__BAD ){
axflags[ i ] = 1;
break;
}
}
}
}
}
}
ps2 = astAnnul( ps2 );
Expand All @@ -5451,8 +5506,131 @@ static AstRegion *MapRegion( AstRegion *this, AstMapping *map0,
}
}

/* Assume we will be using the supplied Region (this) and Mapping (map). */
usethis = astClone( this );
usemap = astClone( map );

/* If the new Frame for the Region has fewer axes than the old Frame, see
if we can discard some base-frame axes. We only do this if the inverse
transformation would otherwise supply bad values for the unused axes.
Using bad axis values is not a good idea as some operations cannot be
performed if any bad values are supplied. Also having more axes than
needed is bad as it results in fewer mesh points per axis. */
if( nnew < nold ) {

/* First invert the Mapping since astMapSplit only allows selection of
inputs, and we want to select outputs. */
astInvert( map );

/* Create an array holding the indices of the required inputs. */
inax = astMalloc( nnew*sizeof(int) );
if( astOK ) for( i = 0; i < nnew; i++ ) inax[i] = i;

/* Attempt to split the mapping to extract a new mapping that omits any
unnecessary outputs (i.e. outputs that are indepenent of the selected
inputs). Check the Mapping was split successfully. */
outax = astMapSplit( map, nnew, inax, &tmap );
if( outax ) {

/* Get the number of old axes that have been retained in the Mapping. */
nkept = astGetNout( tmap );

/* Now we need to ensure that any unused axes for which the Mapping
creates non-bad values are retained (such values are significant
since they may determine whether the new Region is null or not).
We only need to do this if any of the outputs that were split off
above have been found to generate good values, as indicated by the
"axflags" array. Count the number of extra axes that need to be kept,
over and above those that are kept by the Mapping returned by
astMapSplit above. */
ntotal = 0;
keep = NULL;
if( axflags ) {
keep = astMalloc( nold*sizeof(int) );
if( astOK ) {

/* Loop round each axis in the supplied Region. */
for( i = 0; i < nold; i++ ) {

/* Search the "outax" array to see if this axis was retained by astMapSplit.
If it was, append its index to the total list of axes to keep. */
ok = 0;
for( j = 0; j < nkept; j++ ) {
if( outax[ j ] == i ) {
keep[ ntotal++ ] = i;
ok = 1;
break;
}
}

/* If the axis was not retained by astMapSplit, but generates good axis
values, also append its index to the total list of axes to keep. */
if( !ok && axflags[ i ] ) {
keep[ ntotal++ ] = i;
}
}
}
}

/* If there are no extra axes to keep, then the Mapping returned by
astMapSplit above can be used as it is. */
if( ntotal == nkept ) {

/* The values in the "outax" array will hold the zero-based indices of
the original old axes that are retained by the new Mapping. We need to
create a copy of the supplied Region that includes only these axes. */
usethis = astAnnul( usethis );
usethis = astPickAxes( this, nkept, outax, NULL );

/* Use the temportary Mapping returned by astMapSplit in place of the
supplied Mapping (remember to invert to undo the inversion performed
above). */
usemap = astAnnul( usemap );
usemap = astClone( tmap );
astInvert( usemap );

/* Free temporary resources. */
tmap = astAnnul( tmap );
outax = astFree( outax );

/* If we need to retain some extra axes because they generate good values
(even though they are independent of the new Frame axes)... */
} else if( ntotal > nkept ){

/* We know the old Frame axes that we want to keep, so use astMapSplit
in the opposite direction - i.e. use it on the Mapping form old to
new. */
astInvert( map );

tmap = astAnnul( tmap );
outax = astFree( outax );

outax = astMapSplit( map, ntotal, keep, &tmap );
if( outax ) {

usethis = astAnnul( usethis );
usethis = astPickAxes( this, ntotal, keep, NULL );

usemap = astAnnul( usemap );
usemap = astClone( tmap );

tmap = astAnnul( tmap );
outax = astFree( outax );

}

astInvert( map );
}
keep = astFree( keep );
}
inax = astFree( inax );

/* Invert the Mapping again to bring it back to its original state. */
astInvert( map );
}

/* Take a deep copy of the supplied Region. */
result = astCopy( this );
result = astCopy( usethis );

/* Get a pointer to the encapsulated FrameSet. */
if( astOK ) {
Expand All @@ -5461,7 +5639,7 @@ static AstRegion *MapRegion( AstRegion *this, AstMapping *map0,
/* Add in the new Frame and Mapping. First note the index of the original
current Frame. */
icurr = astGetCurrent( fs );
astAddFrame( fs, AST__CURRENT, map, frame );
astAddFrame( fs, AST__CURRENT, usemap, frame );

/* Remove the original current Frame. */
astRemoveFrame( fs, icurr );
Expand All @@ -5474,11 +5652,14 @@ static AstRegion *MapRegion( AstRegion *this, AstMapping *map0,

/* Since the Mapping has been changed, any cached information calculated
on the basis of the Mapping properties may no longer be up to date. */
astResetCache( this );
astResetCache( result );

/* Free resources */
usemap = astAnnul( usemap );
usethis = astAnnul( usethis );
map = astAnnul( map );
frame = astAnnul( frame );
axflags = astFree( axflags );

/* If not OK, annul the returned pointer. */
if( !astOK ) result = astAnnul( result );
Expand Down
16 changes: 15 additions & 1 deletion sun_master.tex
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21740,8 +21740,22 @@ \subsection{\xlabel{changes}\xlabel{list_of_most_recent_changes}Changes
\item A NormMap that encapsulates a basic Frame will now be simplified to a
UnitMap.

\end{enumerate}
\item If the
c+
astMapRegion
c-
f+
AST\_MAPREGION
f-
method is used to map a Region into a new Frame that has fewer axes than
the original Region, and if the inverse transformation of the supplied
Mapping does not specify a value for the missing axes, then those axes
are removed entirely from the Region. Previously they were retained, but
automatically supplied with bad values. This affects the number of mesh
points per axes for such Regions, and so affects the accuracy of overlap
determination.

\end{enumerate}

Programs which are statically linked will need to be re-linked in
order to take advantage of these new facilities.
Expand Down

0 comments on commit 21fc804

Please sign in to comment.