diff options
| author | William Joye <wjoye@cfa.harvard.edu> | 2019-05-10 18:41:55 (GMT) |
|---|---|---|
| committer | William Joye <wjoye@cfa.harvard.edu> | 2019-05-10 18:41:55 (GMT) |
| commit | 75fab9d80911f74aaab738fa9ab8a4f9b0f57a6b (patch) | |
| tree | 8d0d428ac62291f834ea8927bfa82c115ff1689d /ast/matrixmap.c | |
| parent | 0609bafa2688ce94fbfdc1ea496b20c3fac732f1 (diff) | |
| download | blt-75fab9d80911f74aaab738fa9ab8a4f9b0f57a6b.zip blt-75fab9d80911f74aaab738fa9ab8a4f9b0f57a6b.tar.gz blt-75fab9d80911f74aaab738fa9ab8a4f9b0f57a6b.tar.bz2 | |
upgrade ast 8.7.1
Diffstat (limited to 'ast/matrixmap.c')
| -rw-r--r-- | ast/matrixmap.c | 5767 |
1 files changed, 0 insertions, 5767 deletions
diff --git a/ast/matrixmap.c b/ast/matrixmap.c deleted file mode 100644 index ed9600b..0000000 --- a/ast/matrixmap.c +++ /dev/null @@ -1,5767 +0,0 @@ -/* -*class++ -* Name: -* MatrixMap - -* Purpose: -* Map coordinates by multiplying by a matrix. - -* Constructor Function: -c astMatrixMap -f AST_MATRIXMAP - -* Description: -* A MatrixMap is form of Mapping which performs a general linear -* transformation. Each set of input coordinates, regarded as a -* column-vector, are pre-multiplied by a matrix (whose elements -* are specified when the MatrixMap is created) to give a new -* column-vector containing the output coordinates. If appropriate, -* the inverse transformation may also be performed. - -* Inheritance: -* The MatrixMap class inherits from the Mapping class. - -* Attributes: -* The MatrixMap class does not define any new attributes beyond -* those which are applicable to all Mappings. - -* Functions: -c The MatrixMap class does not define any new functions beyond those -f The MatrixMap class does not define any new routines beyond those -* which are applicable to all Mappings. - -* Copyright: -* Copyright (C) 1997-2006 Council for the Central Laboratory of the -* Research Councils -* Copyright (C) 2009 Science & Technology Facilities Council. -* All Rights Reserved. - -* Licence: -* This program is free software: you can redistribute it and/or -* modify it under the terms of the GNU Lesser General Public -* License as published by the Free Software Foundation, either -* version 3 of the License, or (at your option) any later -* version. -* -* This program is distributed in the hope that it will be useful, -* but WITHOUT ANY WARRANTY; without even the implied warranty of -* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -* GNU Lesser General Public License for more details. -* -* You should have received a copy of the GNU Lesser General -* License along with this program. If not, see -* <http://www.gnu.org/licenses/>. - -* Authors: -* DSB: D.S. Berry (Starlink) -* RFWS: R.F. Warren-Smith (Starlink) - -* History: -* 9-FEB-1996 (DSB): -* Original version. -* 13-NOV-1996 (DSB): -* Updated to support attributes, I/O and an external interface. -* 3-JUN-1997 (DSB): -* astMtrMult and astMtrRot made protected instead of public. -* 16-JUN-1997 (RFWS): -* Tidied public prologues. -* 24-JUN-1997 (DSB): -* Zero returned for coordinates which are indeterminate as a -* result of using an inverted, non-square, diagonal matrix. -* 10-OCT-1997 (DSB): -* o The inverse matrix is no longer dumped by the Dump function. -* Instead, it is re-calculated by the Load function. -* o The description of argument "form" in astMatrixMap corrected -* to indicate that a value of 2 produces a unit matrix. -* o String values used to represent choices externally, instead -* of integers. -* 24-NOV-1997 (DSB): -* Use of error code AST__OPT replaced by AST__RDERR. -* 28-JAN-1998 (DSB): -* Bug fix in astMtrMult: the matrix (forward or inverse) used for -* the "a" MatrixMap was determined by the Invert flag of the other -* ("this") MatrixMap. -* 14-APR-1998 (DSB): -* Bug fix in Dump. Previously, matrix elements with value AST__BAD -* were explicitly written out. Now they are not written out, since -* AST__BAD can have different values on different machines. Missing -* elements default to AST__BAD when read back in using astLoadMatrixMap. -* 20-APR-1998 (DSB): -* Bug fix in astLoadMatrixMap: initialise the pointer to the inverse -* matrix array to NULL if no inverse matrix is needed. -* 25-AUG-1998 (DSB): -* - Transform changed so that bad input axis values are not -* propagated to output axes which are independant of the input axis. -* - CompressMatrix changed to allow a tolerance of DBL_EPSILON when -* determining if a matrix is a unit matrix, or a diagonal matrix. -* - MapMerge changed to allow MatrixMaps to swap with PermMaps -* in order to move the MatrixMap closer to a Mapping with which it -* could merge. -* 22-FEB-1999 (DSB): -* Changed logic of MapMerge to avoid infinite looping. -* 5-MAY-1999 (DSB): -* More corrections to MapMerge: Cleared up errors in the use of the -* supplied invert flags, and corrected logic for deciding which -* neighbouring Mapping to swap with. -* 16-JUL-1999 (DSB): -* Fixed memory leaks in MatWin and MapMerge. -* 8-JAN-2003 (DSB): -* Changed private InitVtab method to protected astInitatrixMapVtab -* method. -* 11-SEP-2003 (DSB): -* Increased tolerance on checks for unit matrices within -* CompressMatrix. Now uses sqrt(DBL_EPSILON)*diag (previously was -* DBL_EPSILON*DIAG ). -* 10-NOV-2003 (DSB): -* Modified functions which swap a MatrixMap with another Mapping -* (e.g. MatSwapPerm, etc), to simplify the returned Mappings. -* 13-JAN-2003 (DSB): -* Modified the tolerance used by CompressMatrix when checking for -* zero matrix elements. Old system compared each element to thre -* size of the diagonal, but different scalings on different axes could -* cause this to trat as zero values which should nto be treated as -* zero. -* 23-APR-2004 (DSB): -* Changes to simplification algorithm. -* 8-JUL-2004 (DSB): -* astMtrMult - Report an error if either MatrixMap does not have a -* defined forward transformation. -* 1-SEP-2004 (DSB): -* Ensure do1 and do2 are initialised before use in MapMerge. -* 7-SEP-2005 (DSB): -* Take account of the Invert flag when using the zoom factor from -* a ZoomMap. -* 14-FEB-2006 (DSB): -* Correct row/col confusion in CompressMatrix. -* 15-MAR-2006 (DSB): -* Override astEqual. -* 15-MAR-2009 (DSB): -* MapSplit: Only create the returned Mapping if it would have some -* outputs. Also, do not create the returned Mapping if any output -* depends on a mixture of selected and unselected inputs. -* 16-JUL-2009 (DSB): -* MatPerm: Fix memory leak (mm2 was not being annulled). -* 2-OCT-2012 (DSB): -* - Check for Infs as well as NaNs. -* - In MapSplit do not split the MatrixMap if the resulting -* matrix would contain only bad elements. -* - Report an error if an attempt is made to create a MatrixMap -* containing only bad elements. -* 4-NOV-2013 (DSB): -* Allow a full form MatrixMap to be simplified to a diagonal form -* MatrixMap if all the off-diagonal values are zero. -* 23-APR-2015 (DSB): -* Improve MapMerge. If a MatrixMap can merge with its next-but-one -* neighbour, then swap the MatrixMap with its neighbour, so that -* it is then next its next-but-one neighbour, and then merge the -* two Mappings into a single Mapping. Previously, only the swap -* was performed - not the merger. And the swap was only performed -* if the intervening neighbour could not itself merge. This could -* result in an infinite simplification loop, which was detected by -* CmpMap and and aborted, resulting in no useful simplification. -* 15-JUN-2017 (DSB): -* A diagonal MatrixMap in which the diagonal elements are all zero -* cannot be simplified to a ZoomMap, since ZoomMaps cannot have -* zero zoom factor. -* 16-JUN-2017 (DSB): -* Fix error checking bug in MtrMult - it was checking for the -* inverse transformation of "this" instead of the forward -* transformation of "a". -* 7-NOW-2017 (DSB): -* Allow a diagonal MatrixMap to merge with a WinMap. -* 5-JUN-2018 (DSB): -* Include the inverse matrix in the dump of a MatrixMap. Previously, -* the inverse matrix was calculated afresh using function InvertMatrix -* when a MatrixMap was read from a dump. However this could introduce -* small round-trip errors if the inverse matrix in the original -* MatrixMap was created by astMtrRot etc, rather than the InvertMatrix -* function. -*class-- -*/ - -/* Module Macros. */ -/* ============== */ -/* Set the name of the class we are implementing. This indicates to - the header files that define class interfaces that they should make - "protected" symbols available. */ -#define astCLASS MatrixMap - -/* Define identifiers for the different forms of matrix storage. */ -#define FULL 0 -#define DIAGONAL 1 -#define UNIT 2 - -/* Include files. */ -/* ============== */ -/* Interface definitions. */ -/* ---------------------- */ - -#include "globals.h" /* Thread-safe global data access */ -#include "error.h" /* Error reporting facilities */ -#include "memory.h" /* Memory allocation facilities */ -#include "object.h" /* Base Object class */ -#include "pointset.h" /* Sets of points/coordinates */ -#include "mapping.h" /* Coordinate mappings (parent class) */ -#include "matrixmap.h" /* Interface definition for this class */ -#include "pal.h" /* SLALIB function definitions */ -#include "permmap.h" -#include "zoommap.h" -#include "unitmap.h" -#include "winmap.h" - -/* Error code definitions. */ -/* ----------------------- */ -#include "ast_err.h" /* AST error codes */ - -/* C header files. */ -/* --------------- */ -#include <ctype.h> -#include <math.h> -#include <stdio.h> -#include <stdlib.h> -#include <string.h> - -/* Module Variables. */ -/* ================= */ - -/* Address of this static variable is used as a unique identifier for - member of this class. */ -static int class_check; -static const char *Form[3] = { "Full", "Diagonal", "Unit" }; /* Text values - used to represent storage form externally */ - -/* Pointers to parent class methods which are extended by this class. */ -static AstPointSet *(* parent_transform)( AstMapping *, AstPointSet *, int, AstPointSet *, int * ); -static int *(* parent_mapsplit)( AstMapping *, int, const int *, AstMapping **, int * ); - - -#ifdef THREAD_SAFE -/* Define how to initialise thread-specific globals. */ -#define GLOBAL_inits \ - globals->Class_Init = 0; - -/* Create the function that initialises global data for this module. */ -astMAKE_INITGLOBALS(MatrixMap) - -/* Define macros for accessing each item of thread specific global data. */ -#define class_init astGLOBAL(MatrixMap,Class_Init) -#define class_vtab astGLOBAL(MatrixMap,Class_Vtab) - - -#include <pthread.h> - - -#else - - -/* Define the class virtual function table and its initialisation flag - as static variables. */ -static AstMatrixMapVtab class_vtab; /* Virtual function table */ -static int class_init = 0; /* Virtual function table initialised? */ - -#endif - -/* External Interface Function Prototypes. */ -/* ======================================= */ -/* The following functions have public prototypes only (i.e. no - protected prototypes), so we must provide local prototypes for use - within this module. */ -AstMatrixMap *astMatrixMapId_( int, int, int, const double [], const char *, ... ); - -/* Prototypes for Private Member Functions. */ -/* ======================================== */ -static AstMatrixMap *MatMat( AstMapping *, AstMapping *, int, int, int * ); -static AstMatrixMap *MatPerm( AstMatrixMap *, AstPermMap *, int, int, int, int * ); -static AstMatrixMap *MatZoom( AstMatrixMap *, AstZoomMap *, int, int, int * ); -static AstMatrixMap *MtrMult( AstMatrixMap *, AstMatrixMap *, int * ); -static AstMatrixMap *MtrRot( AstMatrixMap *, double, const double[], int * ); -static AstPointSet *Transform( AstMapping *, AstPointSet *, int, AstPointSet *, int * ); -static AstWinMap *MatWin2( AstMatrixMap *, AstWinMap *, int, int, int, int * ); -static double *InvertMatrix( int, int, int, double *, int * ); -static double Rate( AstMapping *, double *, int, int, int * ); -static int Equal( AstObject *, AstObject *, int * ); -static int FindString( int, const char *[], const char *, const char *, const char *, const char *, int * ); -static int Ustrcmp( const char *, const char *, int * ); -static int GetTranForward( AstMapping *, int * ); -static int GetIsLinear( AstMapping *, int * ); -static int GetTranInverse( AstMapping *, int * ); -static int CanSwap( AstMapping *, AstMapping *, int, int, int *, int * ); -static int MapMerge( AstMapping *, int, int, int *, AstMapping ***, int **, int * ); -static int PermOK( AstMapping *, int * ); -static int ScalingRowCol( AstMatrixMap *, int, int * ); -static void CompressMatrix( AstMatrixMap *, int * ); -static void Copy( const AstObject *, AstObject *, int * ); -static void Delete( AstObject *obj, int * ); -static void Dump( AstObject *, AstChannel *, int * ); -static void ExpandMatrix( AstMatrixMap *, int * ); -static void MatWin( AstMapping **, int *, int, int * ); -static void MatPermSwap( AstMapping **, int *, int, int * ); -static void PermGet( AstPermMap *, int **, int **, double **, int * ); -static void SMtrMult( int, int, int, const double *, double *, double*, int * ); -static int *MapSplit( AstMapping *, int, const int *, AstMapping **, int * ); - -/* Member functions. */ -/* ================= */ -static int CanSwap( AstMapping *map1, AstMapping *map2, int inv1, int inv2, - int *simpler, int *status ){ -/* -* Name: -* CanSwap - -* Purpose: -* Determine if two Mappings could be swapped. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* int CanSwap( AstMapping *map1, AstMapping *map2, int inv1, int inv2, -* int *simpler, int *status ) - -* Class Membership: -* MatrixMap member function - -* Description: -* This function returns a flag indicating if the pair of supplied -* Mappings could be replaced by an equivalent pair of Mappings from the -* same classes as the supplied pair, but in reversed order. Each pair -* of Mappings is considered to be compunded in series. The supplied -* Mapings are not changed in any way. - -* Parameters: -* map1 -* The Mapping to be applied first. -* map2 -* The Mapping to be applied second. -* inv1 -* The invert flag to use with map1. A value of zero causes the forward -* mapping to be used, and a non-zero value causes the inverse -* mapping to be used. -* inv2 -* The invert flag to use with map2. -* simpler -* Addresss of a location at which to return a flag indicating if -* the swapped Mappings would be intrinsically simpler than the -* original Mappings. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* 1 if the Mappings could be swapped, 0 otherwise. - -* Notes: -* - One of the supplied pair of Mappings must be a MatrixMap. -* - A value of 0 is returned if an error has already occurred, or if -* this function should fail for any reason. -*/ - -/* Local Variables: */ - AstMatrixMap *mat; /* Pointer to MatrixMap Mapping */ - AstMapping *nomat; /* Pointer to non-MatrixMap Mapping */ - const char *class1; /* Pointer to map1 class string */ - const char *class2; /* Pointer to map2 class string */ - const char *nomat_class; /* Pointer to non-MatrixMap class string */ - double *consts; /* Pointer to constants array */ - int *inperm; /* Pointer to input axis permutation array */ - int *outperm; /* Pointer to output axis permutation array */ - int i; /* Loop count */ - int invert[ 2 ]; /* Original invert flags */ - int nax; /* No. of in/out coordinates for the MatrixMap */ - int nin; /* No. of input coordinates for the PermMap */ - int nout; /* No. of output coordinates for the PermMap */ - int ret; /* Returned flag */ - -/* Check the global error status. */ - if ( !astOK ) return 0; - -/* Initialise */ - ret = 0; - *simpler = 0; - -/* Temporarily set the Invert attributes of both Mappings to the supplied - values. */ - invert[ 0 ] = astGetInvert( map1 ); - astSetInvert( map1, inv1 ); - - invert[ 1 ] = astGetInvert( map2 ); - astSetInvert( map2, inv2 ); - -/* Get the classes of the two mappings. */ - class1 = astGetClass( map1 ); - class2 = astGetClass( map2 ); - if( astOK ){ - -/* Get a pointer to the MatrixMap and non-MatrixMap Mappings. */ - if( !strcmp( class1, "MatrixMap" ) ){ - mat = (AstMatrixMap *) map1; - nomat = map2; - nomat_class = class2; - } else { - nomat = map1; - mat = (AstMatrixMap *) map2; - nomat_class = class1; - } - -/* Get the number of input axes for the MatrixMap. */ - nax = astGetNin( mat ); - -/* If it is a WinMap, the Mappings can be swapped. */ - if( !strcmp( nomat_class, "WinMap" ) ){ - ret = 1; - -/* If it is a PermMap, the Mappings can be swapped so long as: - 1) all links between input and output axes in the PermMap are - bi-directional. This does not preclude the existence of unconnected - axes, which do not have links (bi-directional or otherwise). - 2) The MatrixMap is square, and invertable. - 3) If the permMap is applied first, then each output of the PermMap - which is assigned a constant value must correspond to a "scaling" row - and column in the MatrixMap. I.e. if PermMap output axis "i" is - assigned a constant value, then row i and column i of the following - MatrixMap must contain only zeros, EXCEPT for the diagonal term (row - i, column i) which must be non-zero. If the Mappings are in the other - order, then the same applies to PermMap input axes assigned a constant - value. */ - -/* Check the other Mapping is a PermMap, and that the MatrixMap is square - and has an inverse. */ - } else if( !strcmp( nomat_class, "PermMap" ) && - nax == astGetNout( mat ) && ( mat->form == UNIT || - ( mat->i_matrix != NULL && - mat->f_matrix != NULL ) ) ) { - -/* Get the number of input and output coordinates for the PermMap. */ - nin = astGetNin( nomat ); - nout = astGetNout( nomat ); - -/* We need to know the axis permutation arrays and constants array for - the PermMap. */ - PermGet( (AstPermMap *) nomat, &outperm, &inperm, &consts, status ); - if( astOK ) { - -/* Indicate we can swap with the PermMap. */ - ret = 1; - -/* Check each output axis. If any links between axes are found which are - not bi-directional, indicate that we cannot swap with the PermMap. */ - for( i = 0; i < nout; i++ ){ - if( outperm[ i ] >= 0 && outperm[ i ] < nin ) { - if( inperm[ outperm[ i ] ] != i ) { - ret = 0; - break; - } - } - } - -/* Check each input axis. If any links between axes are found which are - not bi-directional, indicate that we cannot swap with the PermMap. */ - for( i = 0; i < nin; i++ ){ - if( inperm[ i ] >= 0 && inperm[ i ] < nout ) { - if( outperm[ inperm[ i ] ] != i ) { - ret = 0; - break; - } - } - } - -/* If the PermMap is suitable, check that any constant values fed from the - PermMap into the MatrixMap (in either forward or inverse direction) - are not changed by the MatrixMap. This requires the row and column for - each constant axis to be zeros, ecept for a value of 1.0 on the - diagonal. First deal with the cases where the PermMap is applied - first, so the outputs of the PermMap are fed into the MatrixMap in the - forward direction. */ - if( ret && ( nomat == map1 ) ) { - - if( nout != nax ){ - astError( AST__RDERR, "PermMap produces %d outputs, but the following" - "MatrixMap has %d inputs\n", status, nout, nax ); - ret = 0; - } - -/* Consider each output axis of the PermMap. */ - for( i = 0; i < nout && astOK ; i++ ) { - -/* If this PermMap output is assigned a constant... */ - if( outperm[ i ] < 0 || outperm[ i ] >= nin ) { - -/* Check the i'th row of the MatrixMap is all zero except for the i'th - column which must be non-zero. If not indicate that the MatrixMap cannot - swap with the PermMap and leave the loop. */ - if( !ScalingRowCol( mat, i, status ) ) { - ret = 0; - break; - } - } - } - } - -/* Now deal with the cases where the PermMap is applied second, so the inputs - of the PermMap are fed into the MatrixMap in the inverse direction. */ - if( ret && ( nomat == map2 ) ) { - - if( nin != nax ){ - astError( AST__RDERR, "Inverse PermMap produces %d inputs, but the " - "preceding MatrixMap has %d outputs\n", status, nin, nax ); - ret = 0; - } - -/* Consider each input axis of the PermMap. */ - for( i = 0; i < nin && astOK; i++ ){ - -/* If this PermMap input is assigned a constant (by the inverse Mapping)... */ - if( inperm[ i ] < 0 || inperm[ i ] >= nout ) { - -/* Check the i'th row of the MatrixMap is all zero except for the i'th - column which must be non-zero. If not indicate that the MatrixMap cannot - swap with the PermMap and leave the loop. */ - if( !ScalingRowCol( mat, i, status ) ) { - ret = 0; - break; - } - } - } - } - -/* If we can swap with the PermMap, the swapped Mappings may be - intrinsically simpler than the original mappings. */ - if( ret ) { - -/* If the PermMap precedes the WinMap, this will be the case if the PermMap - has more outputs than inputs. If the WinMap precedes the PermMap, this - will be the case if the PermMap has more inputs than outputs. */ - *simpler = ( nomat == map1 ) ? nout > nin : nin > nout; - } - -/* Free the axis permutation and constants arrays. */ - outperm = (int *) astFree( (void *) outperm ); - inperm = (int *) astFree( (void *) inperm ); - consts = (double *) astFree( (void *) consts ); - } - } - } - -/* Re-instate the original settings of the Invert attributes for the - supplied MatrixMaps. */ - astSetInvert( map1, invert[ 0 ] ); - astSetInvert( map2, invert[ 1 ] ); - -/* Return the answer. */ - return astOK ? ret : 0; -} - -static void CompressMatrix( AstMatrixMap *this, int *status ){ -/* -* Name: -* CompressMatrix - -* Purpose: -* If possible, reduce the amount of storage needed to store a MatrixMap. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* void CompressMatrix( AstMatrixMap *this, int *status ) - -* Class Membership: -* MatrixMap member function. - -* Description: -* The supplid MatrixMap is converted to its most compressed form -* (i.e no element values if it is a unit matrix, diagonal elements only -* if it is a diagonal matrix, or all elements otherwise). - -* Parameters: -* this -* A pointer to the MatrixMap to be compressed. -* status -* Pointer to the inherited status variable. - -*/ - -/* Local Variables: */ - double *a; /* Pointer to next element */ - double *colmax; /* Pointer to array holding column max values */ - double *fmat; /* Pointer to compressed forward matrix */ - double *rowmax; /* Pointer to array holding row max values */ - double mval; /* Matrix element value */ - int i; /* Loop count */ - int j; /* Loop count */ - int k; /* Loop count */ - int ncol; /* No. of columns in forward matrix */ - int ndiag; /* No. of diagonal elements in matrix */ - int new_form; /* Compressed storage form */ - int new_inv; /* New inverse requied? */ - int next_diag; /* Index of next diagonal element */ - int nrow; /* No. of rows in forward matrix */ - -/* Check the global error status. */ - if ( !astOK || !this ) return; - -/* Initialise variables to avoid "used of uninitialised variable" - messages from dumb compilers. */ - new_inv = 0; - -/* Get the dimensions of the forward matrix. */ - if( astGetInvert( this ) ){ - nrow = astGetNin( this ); - ncol = astGetNout( this ); - } else { - ncol = astGetNin( this ); - nrow = astGetNout( this ); - } - -/* Store the number of diagonal elements in the matrix. This is the - minimum of the number of rows and columns. */ - if( ncol < nrow ){ - ndiag = ncol; - } else { - ndiag = nrow; - } - -/* If the MatrixMap is already stored in UNIT form, it cannot be compressed - any further. */ - if( this->form == UNIT){ - return; - -/* Otherwise, if the MatrixMap is stored in DIAGONAL form, it could be - compressed into a UNIT MatrixMap if all the supplied element values are - one. */ - } else if( this->form == DIAGONAL ){ - new_form = UNIT; - for( i = 0; i < ndiag; i++ ){ - if( !astEQUAL( (this->f_matrix)[ i ], 1.0 ) ){ - new_form = DIAGONAL; - break; - } - } - -/* If it can be compressed, change the storage form and free the arrays - holding the diagonal element values. */ - if( new_form == UNIT ) { - this->f_matrix = (double *) astFree( (void *)( this->f_matrix ) ); - this->i_matrix = (double *) astFree( (void *)( this->i_matrix ) ); - this->form = UNIT; - } - -/* Otherwise, a full MatrixMap has been supplied, but this could be stored - in a unit or diagonal MatrixMap if the element values are appropriate. */ - } else { - new_form = FULL; - -/* Find the maximum absolute value in each column. Scale by - sqrt(DBL_EPSILON) to be come a lower limit for non-zero values. */ - colmax = astMalloc( ncol*sizeof( double ) ); - if( colmax ) { - for( j = 0; j < ncol; j++ ) { - colmax[ j ] = 0.0; - i = j; - for( k = 0; k < nrow; k++ ) { - mval = (this->f_matrix)[ i ]; - if( mval != AST__BAD ) { - mval = fabs( mval ); - if( mval > colmax[ j ] ) colmax[ j ] = mval; - } - i += ncol; - } - colmax[ j ] *= sqrt( DBL_EPSILON ); - } - } - -/* Find the maximum absolute value in each row. Scale by - sqrt(DBL_EPSILON) to be come a lower limit for non-zero values. */ - rowmax = astMalloc( nrow*sizeof( double ) ); - if( rowmax ) { - for( k = 0; k < nrow; k++ ) { - rowmax[ k ] = 0.0; - i = k*ncol; - for( j = 0; j < ncol; j++ ) { - mval = (this->f_matrix)[ i ]; - if( mval != AST__BAD ) { - mval = fabs( mval ); - if( mval > rowmax[ k ] ) rowmax[ k ] = mval; - } - i++; - } - rowmax[ k ] *= sqrt( DBL_EPSILON ); - } - } - -/* Check memory can be used */ - if( astOK ) { - -/* Initialise a flag indicating that the inverse matrix does not need to - be re-calculated. */ - new_inv = 0; - -/* Initially assume that the forward matrix is a unit matrix. */ - new_form = UNIT; - -/* Store a pointer to the next matrix element. */ - a = this->f_matrix; - -/* Loop through all the rows in the forward matrix array. */ - for( k = 0; k < nrow; k++ ) { - -/* Loop through all the elements in this column. */ - for( j = 0; j < ncol; j++, a++ ) { - -/* If this element is bad, use full form. */ - if( *a == AST__BAD ) { - new_form = FULL; - -/* Otherwise, if this is a diagonal term, check its value. If it is not one, - then the matrix cannot be a unit matrix, but it could still be a diagonal - matrix. */ - } else { - if( j == k ) { - if( *a != 1.0 && new_form == UNIT ) new_form = DIAGONAL; - -/* If this is not a diagonal element, and the element value is not zero, - then the matrix is not a diagonal matrix. Allow a tolerance of - SQRT(DBL_EPSILON) times the largest value in the same row or column as - the current matrix element. That is, an element must be insignificant - to both its row and its column to be considered as effectively zero. - Replace values less than this limit with zero. */ - } else { - mval = fabs( *a ); - if( mval <= rowmax[ k ] && - mval <= colmax[ j ] ) { - -/* If the element will change value, set a flag indicating that the inverse - matrix needs to be re-calculated. */ - if( *a != 0.0 ) new_inv = 1; - -/* Ensure this element value is zero. */ - *a = 0.0; - - } else { - new_form = FULL; - } - } - } - } - } - } - -/* Free memory. */ - colmax = astFree( colmax ); - rowmax = astFree( rowmax ); - -/* If it can be compressed into a UNIT MatrixMap, change the storage form and - free the arrays holding the element values. */ - if( new_form == UNIT ) { - this->f_matrix = (double *) astFree( (void *)( this->f_matrix ) ); - this->i_matrix = (double *) astFree( (void *)( this->i_matrix ) ); - this->form = UNIT; - -/* Otherwise, if it can be compressed into a DIAGONAL MatrixMap, copy the - diagonal elements from the full forward matrix into a newly allocated - array, use this array to replace the forward matrix array in the MatrixMap, - create a new inverse matrix, and change the storage form. */ - } else if( new_form == DIAGONAL ) { - fmat = astMalloc( sizeof(double)*(size_t)ndiag ); - if( fmat ){ - - next_diag = 0; - for( i = 0; i < ndiag; i++ ){ - fmat[ i ] = (this->f_matrix)[ next_diag ]; - next_diag += ncol + 1; - } - - (void) astFree( (void *) this->f_matrix ); - (void) astFree( (void *) this->i_matrix ); - - this->f_matrix = fmat; - this->i_matrix = InvertMatrix( DIAGONAL, nrow, ncol, fmat, status ); - this->form = DIAGONAL; - - } - -/* Calculate a new inverse matrix if necessary. */ - } else if( new_inv ) { - (void) astFree( (void *) this->i_matrix ); - this->i_matrix = InvertMatrix( FULL, nrow, ncol, this->f_matrix, status ); - } - } - - return; - -} - -static int Equal( AstObject *this_object, AstObject *that_object, int *status ) { -/* -* Name: -* Equal - -* Purpose: -* Test if two MatrixMaps are equivalent. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* int Equal( AstObject *this, AstObject *that, int *status ) - -* Class Membership: -* MatrixMap member function (over-rides the astEqual protected -* method inherited from the astMapping class). - -* Description: -* This function returns a boolean result (0 or 1) to indicate whether -* two MatrixMaps are equivalent. - -* Parameters: -* this -* Pointer to the first Object (a MatrixMap). -* that -* Pointer to the second Object. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* One if the MatrixMaps are equivalent, zero otherwise. - -* Notes: -* - A value of zero will be returned if this function is invoked -* with the global status set, or if it should fail for any reason. -*/ - -/* Local Variables: */ - AstMatrixMap *that; - AstMatrixMap *this; - double *that_matrix; - double *this_matrix; - int i; - int nin; - int nout; - int result; - -/* Initialise. */ - result = 0; - -/* Check the global error status. */ - if ( !astOK ) return result; - -/* Obtain pointers to the two MatrixMap structures. */ - this = (AstMatrixMap *) this_object; - that = (AstMatrixMap *) that_object; - -/* Check the second object is a MatrixMap. We know the first is a - MatrixMap since we have arrived at this implementation of the virtual - function. */ - if( astIsAMatrixMap( that ) ) { - -/* Get the number of inputs and outputs and check they are the same for both. */ - nin = astGetNin( this ); - nout = astGetNout( this ); - if( astGetNout( that ) == nout && astGetNin( that ) == nin ) { - -/* Assume the MatrixMaps are equivalent. */ - result = 1; - -/* Ensure both MatrixMaps are stored in full form. */ - ExpandMatrix( this, status ); - ExpandMatrix( that, status ); - -/* Get pointers to the arrays holding the elements of the forward matrix - for both MatrixMaps. */ - if( astGetInvert( this ) ) { - this_matrix = this->i_matrix; - } else { - this_matrix = this->f_matrix; - } - - if( astGetInvert( that ) ) { - that_matrix = that->i_matrix; - } else { - that_matrix = that->f_matrix; - } - -/* If either of the above arrays is not available, try to get the inverse - matrix arrays. */ - if( !this_matrix || !that_matrix ) { - if( astGetInvert( this ) ) { - this_matrix = this->f_matrix; - } else { - this_matrix = this->i_matrix; - } - - if( astGetInvert( that ) ) { - that_matrix = that->f_matrix; - } else { - that_matrix = that->i_matrix; - } - } - -/* If both arrays are now available compare their elements. */ - if( this_matrix && that_matrix ) { - result = 1; - for( i = 0; i < nin*nout; i++ ) { - if( !astEQUAL( this_matrix[ i ], that_matrix[ i ] ) ){ - result = 0; - break; - } - } - } - -/* Ensure the supplied MatrixMaps are stored back in compressed form. */ - CompressMatrix( this, status ); - CompressMatrix( that, status ); - } - } - -/* If an error occurred, clear the result value. */ - if ( !astOK ) result = 0; - -/* Return the result, */ - return result; -} - -static void ExpandMatrix( AstMatrixMap *this, int *status ){ -/* -* Name: -* ExpandMatrix - -* Purpose: -* Ensure the MatrixMap is stored in full (non-compressed) form. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* void ExpandMatrix( AstMatrixMap *this, int *status ) - -* Class Membership: -* MatrixMap member function. - -* Description: -* If the supplid MatrixMap is stored in a compressed form (i.e no -* element values if it is a unit matrix, diagonal elements only -* if it is a diagonal matrix), it is expanded into a full MatrixMap -* in which all elements are stored. - -* Parameters: -* this -* A pointer to the MatrixMap to be expanded. -* status -* Pointer to the inherited status variable. - -*/ - -/* Local Variables: */ - double *fmat; /* Pointer to full forward matrix */ - double *imat; /* Pointer to full inverse matrix */ - int i; /* Loop count */ - int ncol; /* No. of columns in forward matrix */ - int ndiag; /* No. of diagonal elements in matrix */ - int nrow; /* No. of rows in forward matrix */ - -/* Check the global error status. Also return if the MatrixMap - pointer is null. */ - if ( !astOK || !this ) return; - -/* Return without action if the MatrixMap is already in full form. */ - if( this->form == FULL ) return; - -/* Get the dimensions of the forward matrix. */ - if( astGetInvert( this ) ){ - nrow = astGetNin( this ); - ncol = astGetNout( this ); - } else { - ncol = astGetNin( this ); - nrow = astGetNout( this ); - } - -/* Store the number of diagonal elements. */ - if( nrow > ncol ){ - ndiag = ncol; - } else { - ndiag = nrow; - } - -/* Allocate arrays to hold the full forward and inverse matrices. */ - fmat = (double *) astMalloc( sizeof( double )*(size_t)( nrow*ncol ) ); - imat = (double *) astMalloc( sizeof( double )*(size_t)( nrow*ncol ) ); - if( imat && fmat ){ - -/* Fill them both with zeros. */ - for( i = 0; i < nrow*ncol; i++ ) { - fmat[ i ] = 0.0; - imat[ i ] = 0.0; - } - -/* If a unit MatrixMap was supplied, put ones on the diagonals. */ - if( this->form == UNIT ){ - for( i = 0; i < ndiag; i++ ) { - fmat[ i*( ncol + 1 ) ] = 1.0; - imat[ i*( nrow + 1 ) ] = 1.0; - } - -/* If a diagonal MatrixMap was supplied, copy the diagonal terms from - the supplied MatrixMap. */ - } else if( this->form == DIAGONAL ){ - for( i = 0; i < ndiag; i++ ) { - fmat[ i*( ncol + 1 ) ] = (this->f_matrix)[ i ]; - imat[ i*( nrow + 1 ) ] = (this->i_matrix)[ i ]; - } - } - -/* Free any existing arrays in the MatrixMap and store the new ones. */ - (void) astFree( (void *) this->f_matrix ); - (void) astFree( (void *) this->i_matrix ); - - this->f_matrix = fmat; - this->i_matrix = imat; - -/* Update the storage form. */ - this->form = FULL; - -/* If either of the new matrices could not be allocated, ensure that - both have been freed. */ - } else { - fmat = (double *) astFree( (void *) fmat ); - imat = (double *) astFree( (void *) imat ); - } - - return; - -} - -static int FindString( int n, const char *list[], const char *test, - const char *text, const char *method, - const char *class, int *status ){ -/* -* Name: -* FindString - -* Purpose: -* Find a given string within an array of character strings. - -* Type: -* Private function. - -* Synopsis: -* #include "matrix.h" -* int FindString( int n, const char *list[], const char *test, -* const char *text, const char *method, const char *class, int *status ) - -* Class Membership: -* MatrixMap method. - -* Description: -* This function identifies a supplied string within a supplied -* array of valid strings, and returns the index of the string within -* the array. The test option may not be abbreviated, but case is -* insignificant. - -* Parameters: -* n -* The number of strings in the array pointed to be "list". -* list -* A pointer to an array of legal character strings. -* test -* A candidate string. -* text -* A string giving a description of the object, parameter, -* attribute, etc, to which the test value refers. -* This is only for use in constructing error messages. It should -* start with a lower case letter. -* method -* Pointer to a string holding the name of the calling method. -* This is only for use in constructing error messages. -* class -* Pointer to a string holding the name of the supplied object class. -* This is only for use in constructing error messages. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* The index of the identified string within the supplied array, starting -* at zero. - -* Notes: -* - A value of -1 is returned if an error has already occurred, or -* if this function should fail for any reason (for instance if the -* supplied option is not specified in the supplied list). - -*/ - -/* Local Variables: */ - int ret; /* The returned index */ - -/* Check global status. */ - if( !astOK ) return -1; - -/* Compare the test string with each element of the supplied list. Leave - the loop when a match is found. */ - for( ret = 0; ret < n; ret++ ) { - if( !Ustrcmp( test, list[ ret ], status ) ) break; - } - -/* Report an error if the supplied test string does not match any element - in the supplied list. */ - if( ret >= n ) { - astError( AST__RDERR, "%s(%s): Illegal value '%s' supplied for %s.", status, - method, class, test, text ); - ret = -1; - } - -/* Return the answer. */ - return ret; -} - -static int GetIsLinear( AstMapping *this_mapping, int *status ){ -/* -* Name: -* GetIsLinear - -* Purpose: -* Return the value of the IsLinear attribute for a MatrixMap. - -* Type: -* Private function. - -* Synopsis: -* #include "mapping.h" -* void GetIsLinear( AstMapping *this, int *status ) - -* Class Membership: -* MatrixMap member function (over-rides the protected astGetIsLinear -* method inherited from the Mapping class). - -* Description: -* This function returns the value of the IsLinear attribute for a -* Frame, which is always one. - -* Parameters: -* this -* Pointer to the MatrixMap. -* status -* Pointer to the inherited status variable. -*/ - return 1; -} - -static int Ustrcmp( const char *a, const char *b, int *status ){ -/* -* Name: -* Ustrncmp - -* Purpose: -* A case blind version of strcmp. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* int Ustrcmp( const char *a, const char *b ) - -* Class Membership: -* MatrixMap member function. - -* Description: -* Returns 0 if there are no differences between the two strings, and 1 -* otherwise. Comparisons are case blind. - -* Parameters: -* a -* Pointer to first string. -* b -* Pointer to second string. - -* Returned Value: -* Zero if the strings match, otherwise one. - -* Notes: -* - This function does not consider the sign of the difference between -* the two strings, whereas "strcmp" does. -* - This function attempts to execute even if an error has occurred. - -*/ - -/* Local Variables: */ - const char *aa; /* Pointer to next "a" character */ - const char *bb; /* Pointer to next "b" character */ - int ret; /* Returned value */ - -/* Initialise the returned value to indicate that the strings match. */ - ret = 0; - -/* Initialise pointers to the start of each string. */ - aa = a; - bb = b; - -/* Loop round each character. */ - while( 1 ){ - -/* We leave the loop if either of the strings has been exhausted. */ - if( !(*aa ) || !(*bb) ){ - -/* If one of the strings has not been exhausted, indicate that the - strings are different. */ - if( *aa || *bb ) ret = 1; - -/* Break out of the loop. */ - break; - -/* If neither string has been exhausted, convert the next characters to - upper case and compare them, incrementing the pointers to the next - characters at the same time. If they are different, break out of the - loop. */ - } else { - - if( toupper( (int) *(aa++) ) != toupper( (int) *(bb++) ) ){ - ret = 1; - break; - } - - } - - } - -/* Return the result. */ - return ret; - -} - -void astInitMatrixMapVtab_( AstMatrixMapVtab *vtab, const char *name, int *status ) { -/* -*+ -* Name: -* astInitMatrixMapVtab - -* Purpose: -* Initialise a virtual function table for a MatrixMap. - -* Type: -* Protected function. - -* Synopsis: -* #include "matrixmap.h" -* void astInitMatrixMapVtab( AstMatrixMapVtab *vtab, const char *name ) - -* Class Membership: -* MatrixMap vtab initialiser. - -* Description: -* This function initialises the component of a virtual function -* table which is used by the MatrixMap class. - -* Parameters: -* vtab -* Pointer to the virtual function table. The components used by -* all ancestral classes will be initialised if they have not already -* been initialised. -* name -* Pointer to a constant null-terminated character string which contains -* the name of the class to which the virtual function table belongs (it -* is this pointer value that will subsequently be returned by the Object -* astClass function). -*- -*/ - -/* Local Variables: */ - astDECLARE_GLOBALS /* Pointer to thread-specific global data */ - AstObjectVtab *object; /* Pointer to Object component of Vtab */ - AstMappingVtab *mapping; /* Pointer to Mapping component of Vtab */ - -/* Check the local error status. */ - if ( !astOK ) return; - -/* Get a pointer to the thread specific global data structure. */ - astGET_GLOBALS(NULL); - -/* Initialize the component of the virtual function table used by the - parent class. */ - astInitMappingVtab( (AstMappingVtab *) vtab, name ); - -/* Store a unique "magic" value in the virtual function table. This - will be used (by astIsAMatrixMap) to determine if an object belongs - to this class. We can conveniently use the address of the (static) - class_check variable to generate this unique value. */ - vtab->id.check = &class_check; - vtab->id.parent = &(((AstMappingVtab *) vtab)->id); - -/* Initialise member function pointers. */ -/* ------------------------------------ */ -/* Store pointers to the member functions (implemented here) that provide - virtual methods for this class. */ - vtab->MtrRot = MtrRot; - vtab->MtrMult = MtrMult; - -/* Save the inherited pointers to methods that will be extended, and - replace them with pointers to the new member functions. */ - object = (AstObjectVtab *) vtab; - mapping = (AstMappingVtab *) vtab; - - parent_transform = mapping->Transform; - mapping->Transform = Transform; - - parent_mapsplit = mapping->MapSplit; - mapping->MapSplit = MapSplit; - -/* Store replacement pointers for methods which will be over-ridden by - new member functions implemented here. */ - object->Equal = Equal; - mapping->GetIsLinear = GetIsLinear; - mapping->GetTranForward = GetTranForward; - mapping->GetTranInverse = GetTranInverse; - mapping->MapMerge = MapMerge; - mapping->Rate = Rate; - -/* Declare the destructor and copy constructor. */ - astSetDelete( (AstObjectVtab *) vtab, Delete ); - astSetCopy( (AstObjectVtab *) vtab, Copy ); - -/* Declare the class dump function. */ - astSetDump( vtab, Dump, "MatrixMap", "Matrix transformation" ); - -/* If we have just initialised the vtab for the current class, indicate - that the vtab is now initialised, and store a pointer to the class - identifier in the base "object" level of the vtab. */ - if( vtab == &class_vtab ) { - class_init = 1; - astSetVtabClassIdentifier( vtab, &(vtab->id) ); - } -} - - -static double *InvertMatrix( int form, int nrow, int ncol, double *matrix, int *status ){ -/* -* Name: -* InvertMatrix - -* Purpose: -* Invert a suplied matrix. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* double *InvertMatrix( int form, int nrow, int ncol, double *matrix, int *status ) - -* Class Membership: -* MatrixMap member function. - -* Description: -* This function returns a pointer to a matrix holding the inverse of -* the supplied matrix, or a NULL pointer if the inverse is not defined. -* The memory to store the inverse matrix is allocated internally, and -* should be freed using astFree when no longer required. -* -* The correspondence between a full matrix and its inverse is only -* unique if the matrix is square, and so a NULL pointer is returned if -* the supplied matrix is not square. - -* Parameters: -* form -* The form of the MatrixMap; UNIT, DIAGONAL or FULL. -* nrow -* Number of rows in the supplied matrix. -* ncol -* Number of columns in the supplied matrix. -* matrix -* A pointer to the input matrix. Elements should be stored in row -* order (i.e. (row 1,column 1 ), (row 1,column 2 )... (row 2,column 1), -* etc). -* status -* Pointer to the inherited status variable. - -* Returned Value: -* Pointer to the output matrix. - -* Notes: -* - A NULL pointer is returned if a unit matrix is supplied. -* - A NULL pointer will be returned if this function is invoked with the -* global error status set, or if it should fail for any reason. -* - No error is reported if the inverse is not defined. -*/ - -/* Local Variables: */ - double det; /* Determinant of supplied matrix */ - double mval; /* Matrix element value */ - double *out; /* Pointer to returned inverse matrix */ - double *vector; /* Pointer to vector used by palDmat */ - int i; /* Matrix element number */ - int *iw; /* Pointer to workspace used by palDmat */ - int nel; /* No. of elements in square matrix */ - int ndiag; /* No. of diagonal elements */ - int ok; /* Zero if any bad matrix values found */ - int sing; /* Zero if matrix is not singular */ - -/* Check the global error status. */ - if ( !astOK ) return NULL; - -/* Return a NULL pointer if the input matrix is NULL. */ - if( !matrix ) return NULL; - -/* If a unit matrix map has been supplied, return NULL. */ - if( form == UNIT ){ - return NULL; - -/* If a diagonal matrix has been supplied, allocate an array to hold - the diagonal terms of the inverse matrix. Store the reciprocal - of the input matrix diagonal terms in it. If any of the input diagonal - terms are zero or BAD, set the associated elements of the inverse matrix - BAD. */ - } else if( form == DIAGONAL ){ - if( nrow > ncol ) { - ndiag = ncol; - } else { - ndiag = nrow; - } - - out = (double *) astMalloc( sizeof( double )*(size_t)ndiag ); - - if( out ) { - for( i = 0; i < ndiag; i++ ) { - mval = matrix[ i ]; - if( mval != 0.0 && mval != AST__BAD ){ - out[ i ] = 1.0/mval; - } else { - out[ i ] = AST__BAD; - } - } - } - -/* If a full matrix has been supplied, initialise the returned pointer. */ - } else { - out = NULL; - -/* Check that the matrix is square. */ - if( nrow == ncol ){ - -/* Find the number of elements in the matrix. */ - nel = nrow*ncol; - -/* See if there are any bad values in the matrix. */ - ok = 1; - for ( i=0; i<nel; i++ ) { - if ( matrix[i] == AST__BAD ) { - ok = 0; - break; - } - } - -/* Only continue if there are no bad matrix values. */ - if( ok ) { - -/* Take a copy of the supplied matrix */ - out = (double *) astStore( NULL, (void *) matrix, - astSizeOf( (void *) matrix ) ); - -/* The SLALIB function which inverts the matrix also applies the inverse - matrix to a vector. We are not interested in the vector in this - instance, but we still have to provide one for SLALIB to use. Allocate - memory for the vector. */ - vector = (double *) astMalloc( sizeof(double)*(size_t) nrow ); - -/* If it was allocated succesfully, fill it with zeros. */ - if( astOK ){ - for ( i=0; i<nrow; i++ ) vector[i] = 0.0; - -/* Obtain work space and attempt to invert the matrix using SLALIB, then - free the work space. */ - iw = (int *) astMalloc( sizeof(int)*(size_t) nrow ); - if( astOK ) palDmat( nrow, out, vector, &det, &sing, iw ); - iw = (int *) astFree( (void *) iw ); - - } - -/* If the matrix could not be inverted, free the memory used to hold the - square matrix, and return the NULL pointer. */ - if ( !astOK || sing != 0 ){ - out = (double *) astFree( (void *) out ); - } - -/* Free the memory used to hold the vector. */ - vector = (double *) astFree( (void *) vector ); - } - } - } - -/* Return the pointer. */ - - return out; - -} - -static int MapMerge( AstMapping *this, int where, int series, int *nmap, - AstMapping ***map_list, int **invert_list, int *status ) { -/* -* Name: -* MapMerge - -* Purpose: -* Simplify a sequence of Mappings containing a MatrixMap. - -* Type: -* Private function. - -* Synopsis: -* #include "mapping.h" -* int MapMerge( AstMapping *this, int where, int series, int *nmap, -* AstMapping ***map_list, int **invert_list, int *status ) - -* Class Membership: -* MatrixMap method (over-rides the protected astMapMerge method -* inherited from the Mapping class). - -* Description: -* This function attempts to simplify a sequence of Mappings by -* merging a nominated MatrixMap in the sequence with its neighbours, -* so as to shorten the sequence if possible. -* -* In many cases, simplification will not be possible and the -* function will return -1 to indicate this, without further -* action. -* -* In most cases of interest, however, this function will either -* attempt to replace the nominated MatrixMap with a Mapping which it -* considers simpler, or to merge it with the Mappings which -* immediately precede it or follow it in the sequence (both will -* normally be considered). This is sufficient to ensure the -* eventual simplification of most Mapping sequences by repeated -* application of this function. -* -* In some cases, the function may attempt more elaborate -* simplification, involving any number of other Mappings in the -* sequence. It is not restricted in the type or scope of -* simplification it may perform, but will normally only attempt -* elaborate simplification in cases where a more straightforward -* approach is not adequate. - -* Parameters: -* this -* Pointer to the nominated MatrixMap which is to be merged with -* its neighbours. This should be a cloned copy of the MatrixMap -* pointer contained in the array element "(*map_list)[where]" -* (see below). This pointer will not be annulled, and the -* MatrixMap it identifies will not be modified by this function. -* where -* Index in the "*map_list" array (below) at which the pointer -* to the nominated MatrixMap resides. -* series -* A non-zero value indicates that the sequence of Mappings to -* be simplified will be applied in series (i.e. one after the -* other), whereas a zero value indicates that they will be -* applied in parallel (i.e. on successive sub-sets of the -* input/output coordinates). -* nmap -* Address of an int which counts the number of Mappings in the -* sequence. On entry this should be set to the initial number -* of Mappings. On exit it will be updated to record the number -* of Mappings remaining after simplification. -* map_list -* Address of a pointer to a dynamically allocated array of -* Mapping pointers (produced, for example, by the astMapList -* method) which identifies the sequence of Mappings. On entry, -* the initial sequence of Mappings to be simplified should be -* supplied. -* -* On exit, the contents of this array will be modified to -* reflect any simplification carried out. Any form of -* simplification may be performed. This may involve any of: (a) -* removing Mappings by annulling any of the pointers supplied, -* (b) replacing them with pointers to new Mappings, (c) -* inserting additional Mappings and (d) changing their order. -* -* The intention is to reduce the number of Mappings in the -* sequence, if possible, and any reduction will be reflected in -* the value of "*nmap" returned. However, simplifications which -* do not reduce the length of the sequence (but improve its -* execution time, for example) may also be performed, and the -* sequence might conceivably increase in length (but normally -* only in order to split up a Mapping into pieces that can be -* more easily merged with their neighbours on subsequent -* invocations of this function). -* -* If Mappings are removed from the sequence, any gaps that -* remain will be closed up, by moving subsequent Mapping -* pointers along in the array, so that vacated elements occur -* at the end. If the sequence increases in length, the array -* will be extended (and its pointer updated) if necessary to -* accommodate any new elements. -* -* Note that any (or all) of the Mapping pointers supplied in -* this array may be annulled by this function, but the Mappings -* to which they refer are not modified in any way (although -* they may, of course, be deleted if the annulled pointer is -* the final one). -* invert_list -* Address of a pointer to a dynamically allocated array which, -* on entry, should contain values to be assigned to the Invert -* attributes of the Mappings identified in the "*map_list" -* array before they are applied (this array might have been -* produced, for example, by the astMapList method). These -* values will be used by this function instead of the actual -* Invert attributes of the Mappings supplied, which are -* ignored. -* -* On exit, the contents of this array will be updated to -* correspond with the possibly modified contents of the -* "*map_list" array. If the Mapping sequence increases in -* length, the "*invert_list" array will be extended (and its -* pointer updated) if necessary to accommodate any new -* elements. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* If simplification was possible, the function returns the index -* in the "map_list" array of the first element which was -* modified. Otherwise, it returns -1 (and makes no changes to the -* arrays supplied). - -* Notes: -* - A value of -1 will be returned if this function is invoked -* with the global error status set, or if it should fail for any -* reason. -*/ - -/* Local Variables: */ - AstMapping **maplt; /* New mappings list pointer */ - AstMapping *map2; /* Pointer to replacement Mapping */ - AstMapping *mc[2]; /* Copies of supplied Mappings to swap */ - AstMapping *newmap; /* Pointer to replacement MatrixMap */ - AstMapping *smc0; /* Simplied Mapping */ - AstMapping *smc1; /* Simplied Mapping */ - AstMatrixMap *mm; /* Pointer to supplied MatrixMap */ - const char *class1; /* Pointer to first Mapping class string */ - const char *class2; /* Pointer to second Mapping class string */ - const char *nclass; /* Pointer to neighbouring Mapping class */ - double *b; /* Pointer to scale terms */ - double *new_mat; /* Pointer to elements of new MatrixMap */ - double factor; /* Zoom factor for new ZoomMap */ - int *invlt; /* New invert flags list pointer */ - int do1; /* Would a backward swap make a simplification? */ - int do2; /* Would a forward swap make a simplification? */ - int i1; /* Index of first MatrixMap to merge */ - int i2; /* Index of last MatrixMap to merge */ - int i; /* Loop counter */ - int ic[2]; /* Copies of supplied invert flags to swap */ - int invert; /* Should the inverted Mapping be used? */ - int j; /* Loop counter */ - int nin; /* Number of input coordinates for MatrixMap */ - int nmapt; /* No. of Mappings in list */ - int nout; /* Number of output coordinates for MatrixMap */ - int nstep1; /* No. of Mappings backwards to next mergable Mapping */ - int nstep2; /* No. of Mappings forward to next mergable Mapping */ - int result; /* Result value to return */ - int swaphi; /* Can MatrixMap be swapped with higher neighbour? */ - int swaplo; /* Can MatrixMap be swapped with lower neighbour? */ - int zoom; /* Can MatrixMap be replaced by a ZoomMap? */ - -/* Initialise. */ - result = -1; - -/* Check the global error status. */ - if ( !astOK ) return result; - -/* Initialise variables to avoid "used of uninitialised variable" - messages from dumb compilers. */ - i1 = 0; - i2 = 0; - -/* Get the Invert attribute for the specified mapping. */ - invert = astGetInvert( ( *map_list )[ where ] ); - -/* Get the number of input and output axes for the MatrixMap. Swap these - if the supplied invert flag is not the same as the Invert attribute of - the Mapping. */ - if( ( invert && !( *invert_list )[ where ] ) || - ( !invert && ( *invert_list )[ where ] ) ) { - nout = astGetNin( ( *map_list )[ where ] ); - nin = astGetNout( ( *map_list )[ where ] ); - - } else { - nin = astGetNin( ( *map_list )[ where ] ); - nout = astGetNout( ( *map_list )[ where ] ); - } - -/* First of all, see if the MatrixMap can be replaced by a simpler Mapping, - without reference to the neighbouring Mappings in the list. */ -/* ======================================================================*/ - map2 = NULL; - mm = (AstMatrixMap *) ( *map_list )[ where ]; - -/* If the MatrixMap is a square unit matrix, it can be replaced by a - UnitMap. */ - if( mm->form == UNIT && nin == nout ){ - map2 = (AstMapping *) astUnitMap( nin, "", status ); - -/* If the MatrixMap is a square diagonal matrix with equal diagonal - terms, then it can be replaced by a ZoomMap, so long as the - diagonal elements are not all zero. */ - } else if( mm->form == DIAGONAL && nin == nout && - mm->f_matrix && mm->i_matrix && - (mm->f_matrix)[ 0 ] != AST__BAD ){ - zoom = 1; - b = mm->f_matrix + 1; - for( i = 1; i < nin; i++ ){ - if( !astEQUAL( *b, *( b - 1 ) ) ){ - zoom = 0; - break; - } - b++; - } - - if( zoom ){ - if( ( *invert_list )[ where ] ){ - factor = (mm->i_matrix)[ 0 ]; - } else { - factor = (mm->f_matrix)[ 0 ]; - } - - if( factor != 0.0 ){ - map2 = (AstMapping *) astZoomMap( nin, factor, "", status ); - } - } - -/* If the MatrixMap is a full matrix but all off-diagonal elements are - zero, it can be replaced by a diagonal MatrixMap. */ - } else if( mm->form == FULL && nin == nout && mm->f_matrix ){ - new_mat = astMalloc( sizeof( double )*nin ); - b = mm->f_matrix; - for( i = 0; i < nin && new_mat; i++ ){ - for( j = 0; j < nout; j++,b++ ){ - if( i == j ) { - new_mat[ i ] = *b; - } else if( *b != 0.0 ) { - new_mat = astFree( new_mat ); - break; - } - } - } - - if( new_mat ) { - map2 = (AstMapping *) astMatrixMap( nin, nout, 1, new_mat, "", - status ); - new_mat = astFree( new_mat ); - } - } - -/* If the MatrixMap can be replaced, annul the MatrixMap pointer in the - list and replace it with the new Mapping pointer, and indicate that the - forward transformation of the returned Mapping should be used. */ - if( map2 ){ - (void) astAnnul( ( *map_list )[ where ] ); - ( *map_list )[ where ] = map2; - ( *invert_list )[ where ] = 0; - -/* Return the index of the first modified element. */ - result = where; - -/* If the MatrixMap itself could not be simplified, see if it can be merged - with the Mappings on either side of it in the list. */ -/*==========================================================================*/ - } else { - -/* Store the classes of the neighbouring Mappings in the list. */ - class1 = ( where > 0 ) ? astGetClass( ( *map_list )[ where - 1 ] ) : NULL; - class2 = ( where < *nmap - 1 ) ? astGetClass( ( *map_list )[ where + 1 ] ) : NULL; - -/* In series. */ -/* ========== */ - if ( series ) { - -/* We first look to see if the MatrixMap can be merged with one of its - neighbours, resulting in a reduction of one in the number of Mappings - in the list. MatrixMaps can merge directly with another MatrixMap, a - ZoomMap, an invertable PermMap, or a UnitMap. */ - if( class1 && ( !strcmp( class1, "MatrixMap" ) || - !strcmp( class1, "ZoomMap" ) || - !strcmp( class1, "PermMap" ) || - !strcmp( class1, "UnitMap" ) ) ){ - nclass = class1; - i1 = where - 1; - i2 = where; - - } else if( class2 && ( !strcmp( class2, "MatrixMap" ) || - !strcmp( class2, "ZoomMap" ) || - !strcmp( class2, "PermMap" ) || - !strcmp( class2, "UnitMap" ) ) ){ - nclass = class2; - i1 = where; - i2 = where + 1; - - } else { - nclass = NULL; - } - -/* Only some PermMaps can be merged with (those which have consistent - forward and inverse mappings). If this is not one of them, set nclass - NULL to indicate this. */ - if( nclass && !strcmp( nclass, "PermMap" ) && - !PermOK( ( *map_list )[ (i1==where)?i2:i1 ], status ) ) nclass = NULL; - -/* If the MatrixMap is diagonal it can also merge with a WinMap. */ - if( !nclass && mm->form == DIAGONAL) { - if( class1 && ( !strcmp( class1, "WinMap" ) ) ){ - nclass = class1; - i1 = where - 1; - i2 = where; - - } else if( class2 && ( !strcmp( class2, "WinMap" ) ) ){ - nclass = class2; - i1 = where; - i2 = where + 1; - - } - } - -/* If the MatrixMap can merge with one of its neighbours, create the merged - Mapping. */ - if( nclass ){ - - if( !strcmp( nclass, "MatrixMap" ) ){ - newmap = (AstMapping *) MatMat( ( *map_list )[ i1 ], ( *map_list )[ i2 ], - ( *invert_list )[ i1 ], ( *invert_list )[ i2 ], status ); - invert = 0; - - } else if( !strcmp( nclass, "ZoomMap" ) ){ - if( i1 == where ){ - newmap = (AstMapping *) MatZoom( (AstMatrixMap *)( *map_list )[ i1 ], - (AstZoomMap *)( *map_list )[ i2 ], - ( *invert_list )[ i1 ], ( *invert_list )[ i2 ], status ); - } else { - newmap = (AstMapping *) MatZoom( (AstMatrixMap *)( *map_list )[ i2 ], - (AstZoomMap *)( *map_list )[ i1 ], - ( *invert_list )[ i2 ], ( *invert_list )[ i1 ], status ); - } - invert = 0; - - } else if( !strcmp( nclass, "PermMap" ) ){ - if( i1 == where ){ - newmap = (AstMapping *) MatPerm( (AstMatrixMap *)( *map_list )[ i1 ], - (AstPermMap *)( *map_list )[ i2 ], - ( *invert_list )[ i1 ], ( *invert_list )[ i2 ], 1, status ); - } else { - newmap = (AstMapping *) MatPerm( (AstMatrixMap *)( *map_list )[ i2 ], - (AstPermMap *)( *map_list )[ i1 ], - ( *invert_list )[ i2 ], ( *invert_list )[ i1 ], 0, status ); - } - invert = 0; - - } else if( !strcmp( nclass, "WinMap" ) ){ - if( i1 == where ){ - newmap = (AstMapping *) MatWin2( (AstMatrixMap *)( *map_list )[ i1 ], - (AstWinMap *)( *map_list )[ i2 ], - ( *invert_list )[ i1 ], ( *invert_list )[ i2 ], 1, status ); - } else { - newmap = (AstMapping *) MatWin2( (AstMatrixMap *)( *map_list )[ i2 ], - (AstWinMap *)( *map_list )[ i1 ], - ( *invert_list )[ i2 ], ( *invert_list )[ i1 ], 0, status ); - } - invert = 0; - - } else { - newmap = astClone( ( *map_list )[ where ] ); - invert = ( *invert_list )[ where ]; - } - -/* If succesfull... */ - if( astOK ){ - -/* Annul the first of the two Mappings, and replace it with the merged - MatrixMap. Also set the invert flag. */ - (void) astAnnul( ( *map_list )[ i1 ] ); - ( *map_list )[ i1 ] = newmap; - ( *invert_list )[ i1 ] = invert; - -/* Annul the second of the two Mappings, and shuffle down the rest of the - list to fill the gap. */ - (void) astAnnul( ( *map_list )[ i2 ] ); - for ( i = i2 + 1; i < *nmap; i++ ) { - ( *map_list )[ i - 1 ] = ( *map_list )[ i ]; - ( *invert_list )[ i - 1 ] = ( *invert_list )[ i ]; - } - -/* Clear the vacated element at the end. */ - ( *map_list )[ *nmap - 1 ] = NULL; - ( *invert_list )[ *nmap - 1 ] = 0; - -/* Decrement the Mapping count and return the index of the first - modified element. */ - ( *nmap )--; - result = i1; - - } - -/* If the MatrixMap could not merge directly with either of its neighbours, - we consider whether it would be worthwhile to swap the MatrixMap with - either of its neighbours. This can only be done for certain classes - of Mapping (WinMaps and some PermMaps), and will usually require both - Mappings to be modified (unless they are commutative). The advantage of - swapping the order of the Mappings is that it may result in the MatrixMap - being adjacent to a Mapping with which it can merge directly on the next - invocation of this function, thus reducing the number of Mappings - in the list. */ - } else { - -/* Set a flag if we could swap the MatrixMap with its higher neighbour. "do2" - is returned if swapping the Mappings would simplify either of the Mappings. */ - if( where + 1 < *nmap ){ - swaphi = CanSwap( ( *map_list )[ where ], - ( *map_list )[ where + 1 ], - ( *invert_list )[ where ], - ( *invert_list )[ where + 1 ], &do2, status ); - } else { - swaphi = 0; - do2 = 0; - } - -/* If so, step through each of the Mappings which follow the MatrixMap, - looking for a Mapping with which the MatrixMap could merge directly. Stop - when such a Mapping is found, or if a Mapping is found with which the - MatrixMap could definitely not swap. Note the number of Mappings which - separate the MatrixMap from the Mapping with which it could merge (if - any). */ - nstep2 = -1; - if( swaphi ){ - for( i2 = where + 1; i2 < *nmap; i2++ ){ - -/* See if we can merge with this Mapping. If so, note the number of steps - between the two Mappings and leave the loop. */ - nclass = astGetClass( ( *map_list )[ i2 ] ); - if( !strcmp( nclass, "MatrixMap" ) || - !strcmp( nclass, "ZoomMap" ) || - ( !strcmp( nclass, "PermMap" ) && PermOK( ( *map_list )[ i2 ], status ) ) || - !strcmp( nclass, "UnitMap" ) ) { - nstep2 = i2 - where - 1; - break; - } - -/* If there is no chance that we can swap with this Mapping, leave the loop - with -1 for the number of steps to indicate that no merging is possible. - MatrixMaps can swap with WinMaps and some permmaps. */ - if( strcmp( nclass, "WinMap" ) && - strcmp( nclass, "PermMap" ) ) { - break; - } - - } - - } - -/* Do the same working forward from the MatrixMap towards the start of the map - list. */ - if( where > 0 ){ - swaplo = CanSwap( ( *map_list )[ where - 1 ], - ( *map_list )[ where ], - ( *invert_list )[ where - 1 ], - ( *invert_list )[ where ], &do1, status ); - } else { - swaplo = 0; - do1 = 0; - } - - nstep1 = -1; - if( swaplo ){ - for( i1 = where - 1; i1 >= 0; i1-- ){ - - nclass = astGetClass( ( *map_list )[ i1 ] ); - if( !strcmp( nclass, "MatrixMap" ) || - ( !strcmp( nclass, "PermMap" ) && PermOK( ( *map_list )[ i1 ], status ) ) || - !strcmp( nclass, "ZoomMap" ) || - !strcmp( nclass, "UnitMap" ) ) { - nstep1 = where - 1 - i1; - break; - } - - if( strcmp( nclass, "WinMap" ) && - strcmp( nclass, "PermMap" ) ) { - break; - } - - } - - } - -/* Choose which neighbour to swap with so that the MatrixMap moves towards the - nearest Mapping with which it can merge. */ - if( do1 || ( - nstep1 != -1 && ( nstep2 == -1 || nstep2 > nstep1 ) ) ){ - nclass = class1; - i1 = where - 1; - i2 = where; - } else if( do2 || nstep2 != -1 ){ - nclass = class2; - i1 = where; - i2 = where + 1; - } else { - nclass = NULL; - } - -/* If there is a target Mapping in the list with which the MatrixMap could - merge, consider replacing the supplied Mappings with swapped Mappings to - bring the MatrixMap closer to the target Mapping. */ - if( nclass ){ - -/* Swap the Mappings. */ - if (!strcmp( nclass, "WinMap" ) ){ - MatWin( (*map_list) + i1, (*invert_list) + i1, where - i1, status ); - - } else if( !strcmp( nclass, "PermMap" ) ){ - MatPermSwap( (*map_list) + i1, (*invert_list) + i1, where - i1, status ); - } - -/* And then merge them. */ - if( where == i1 && where + 1 < *nmap ) { /* Merging upwards */ - map2 = astClone( (*map_list)[ where + 1 ] ); - nmapt = *nmap - where - 1; - maplt = *map_list + where + 1; - invlt = *invert_list + where + 1; - - (void) astMapMerge( map2, 0, series, &nmapt, &maplt, &invlt ); - map2 = astAnnul( map2 ); - *nmap = where + 1 + nmapt; - - } else if( where - 2 >= 0 ) { /* Merging downwards */ - map2 = astClone( (*map_list)[ where - 2 ] ); - nmapt = *nmap - where + 2; - maplt = *map_list + where - 2 ; - invlt = *invert_list + where - 2; - - (void) astMapMerge( map2, 0, series, &nmapt, &maplt, &invlt ); - map2 = astAnnul( map2 ); - *nmap = where - 2 + nmapt; - } - - result = i1; - -/* If there is no Mapping available for merging, it may still be - advantageous to swap with a neighbour because the swapped Mapping may - be simpler than the original Mappings. For instance, a PermMap may - strip rows of the MatrixMap leaving only a UnitMap. */ - } else if( swaphi || swaplo ) { - -/* Try swapping with each possible neighbour in turn. */ - for( i = 0; i < 2; i++ ) { - -/* Set up the class and pointers for the mappings to be swapped, first - the lower neighbour, then the upper neighbour. */ - if( i == 0 && swaplo ){ - nclass = class1; - i1 = where - 1; - i2 = where; - - } else if( i == 1 && swaphi ){ - nclass = class2; - i1 = where; - i2 = where + 1; - - } else { - nclass = NULL; - } - -/* If we have a Mapping to swap with... */ - if( nclass ) { - -/* Take copies of the Mapping and Invert flag arrays so we do not change - the supplied values. */ - mc[ 0 ] = (AstMapping *) astCopy( ( (*map_list) + i1 )[0] ); - mc[ 1 ] = (AstMapping *) astCopy( ( (*map_list) + i1 )[1] ); - ic[ 0 ] = ( (*invert_list) + i1 )[0]; - ic[ 1 ] = ( (*invert_list) + i1 )[1]; - -/* Swap these Mappings. */ - if( !strcmp( nclass, "WinMap" ) ){ - MatWin( mc, ic, where - i1, status ); - } else if( !strcmp( nclass, "PermMap" ) ){ - MatPermSwap( mc, ic, where - i1, status ); - } - -/* If neither of the swapped Mappings can be simplified further, then there - is no point in swapping the Mappings, so just annul the map copies. */ - smc0 = astSimplify( mc[0] ); - smc1 = astSimplify( mc[1] ); - - if( astGetClass( smc0 ) == astGetClass( mc[0] ) && - astGetClass( smc1 ) == astGetClass( mc[1] ) ) { - - mc[ 0 ] = (AstMapping *) astAnnul( mc[ 0 ] ); - mc[ 1 ] = (AstMapping *) astAnnul( mc[ 1 ] ); - -/* If one or both of the swapped Mappings could be simplified, then annul - the supplied Mappings and return the swapped mappings, storing the index - of the first modified Mapping. */ - } else { - (void ) astAnnul( ( (*map_list) + i1 )[0] ); - (void ) astAnnul( ( (*map_list) + i1 )[1] ); - - ( (*map_list) + i1 )[0] = mc[ 0 ]; - ( (*map_list) + i1 )[1] = mc[ 1 ]; - - ( (*invert_list) + i1 )[0] = ic[ 0 ]; - ( (*invert_list) + i1 )[1] = ic[ 1 ]; - - result = i1; - break; - } - -/* Annul the simplied Mappings */ - smc0 = astAnnul( smc0 ); - smc1 = astAnnul( smc1 ); - - } - } - } - } - } - } - -/* Return the result. */ - return result; -} - -static int *MapSplit( AstMapping *this_map, int nin, const int *in, AstMapping **map, int *status ){ -/* -* Name: -* MapSplit - -* Purpose: -* Create a Mapping representing a subset of the inputs of an existing -* MatrixMap. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* int *MapSplit( AstMapping *this, int nin, const int *in, AstMapping **map, int *status ) - -* Class Membership: -* MatrixMap method (over-rides the protected astMapSplit method -* inherited from the Mapping class). - -* Description: -* This function creates a new Mapping by picking specified inputs from -* an existing MatrixMap. This is only possible if the specified inputs -* correspond to some subset of the MatrixMap outputs. That is, there -* must exist a subset of the MatrixMap outputs for which each output -* depends only on the selected MatrixMap inputs, and not on any of the -* inputs which have not been selected. In addition, outputs that are -* not in this subset must not depend on any selected inputs. If these -* conditions are not met by the supplied MatrixMap, then a NULL Mapping -* is returned. - -* Parameters: -* this -* Pointer to the MatrixMap to be split (the MatrixMap is not actually -* modified by this function). -* nin -* The number of inputs to pick from "this". -* in -* Pointer to an array of indices (zero based) for the inputs which -* are to be picked. This array should have "nin" elements. If "Nin" -* is the number of inputs of the supplied MatrixMap, then each element -* should have a value in the range zero to Nin-1. -* map -* Address of a location at which to return a pointer to the new -* Mapping. This Mapping will have "nin" inputs (the number of -* outputs may be different to "nin"). A NULL pointer will be -* returned if the supplied MatrixMap has no subset of outputs which -* depend only on the selected inputs. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* A pointer to a dynamically allocated array of ints. The number of -* elements in this array will equal the number of outputs for the -* returned Mapping. Each element will hold the index of the -* corresponding output in the supplied MatrixMap. The array should be -* freed using astFree when no longer needed. A NULL pointer will -* be returned if no output Mapping can be created. - -* Notes: -* - If this function is invoked with the global error status set, -* or if it should fail for any reason, then NULL values will be -* returned as the function value and for the "map" pointer. -*/ - -/* Local Variables: */ - AstMatrixMap *this; /* Pointer to MatrixMap structure */ - double *mat; /* Pointer to matrix for supplied MatrixMap */ - double *pmat; /* Pointer to row start in returned matrix */ - double *prow; /* Pointer to row start in supplied matrix */ - double *rmat; /* Pointer to matrix for returned MatrixMap */ - double el; /* Next element value in supplied matrix */ - int *result; /* Pointer to returned array */ - int good; /* Would new matrix contain any good values/ */ - int i; /* Loop count */ - int icol; /* Column index within supplied MatrixMap */ - int iel; /* Index of next element from the input matrix */ - int irow; /* Row index within supplied MatrixMap */ - int isel; /* Does output depend on any selected inputs? */ - int ncol; /* Number of columns (inputs) in supplied MatrixMap */ - int nout; /* Number of outputs in returned MatrixMap */ - int nrow; /* Number of rows (outputs) in supplied MatrixMap */ - int ok; /* Are input indices OK? */ - int sel; /* Does any output depend on selected inputs? */ - int unsel; /* Does any output depend on unselected inputs? */ - -/* Initialise */ - result = NULL; - *map = NULL; - -/* Check the global error status. */ - if ( !astOK ) return result; - -/* Invoke the parent astMapSplit method to see if it can do the job. */ - result = (*parent_mapsplit)( this_map, nin, in, map, status ); - -/* If not, we provide a special implementation here. */ - if( !result ) { - -/* Get a pointer to the MatrixMap structure. */ - this = (AstMatrixMap *) this_map; - -/* Get the number of inputs and outputs. */ - ncol = astGetNin( this ); - nrow = astGetNout( this ); - -/* Check the supplied input indices are usable. */ - ok = 1; - for( i = 0; i < nin; i++ ) { - if( in[ i ] < 0 || in[ i ] >= ncol ) { - ok = 0; - break; - } - } - - if( ok ) { - -/* Ensure the MatrixMap is stored in full form. */ - ExpandMatrix( this, status ); - -/* Allocate the largest array that could be necessary to hold the - returned array of Mapping outputs. */ - result = astMalloc( sizeof(int)*(size_t) nrow ); - -/* Allocate the largest array that could be necessary to hold the - matrix representing the returned MatrixMap. */ - rmat = astMalloc( sizeof(double)*(size_t) (nrow*ncol) ); - -/* Get the matrix which defines the current forward transformation. This - takes into account whether the MatrixMap has been inverted or not. */ - if( astGetInvert( this ) ) { - mat = this->i_matrix; - } else { - mat = this->f_matrix; - } - -/* We cannot create the require Mapping if the matrix is undefined. */ - if( !mat || !astOK ) { - ok = 0; - nout = 0; - good = 0; - -/* Otherwise, loop round all the rows in the matrix. */ - } else { - nout = 0; - good = 0; - pmat = rmat; - iel = 0; - for( irow = 0; irow < nrow; irow++ ) { - -/* Indicate that this output (i.e. row of the matrix) depends on neither - selected nor unselected inputs as yet. */ - sel = 0; - unsel = 0; - -/* Save a pointer to the first element of this row in the MatrixMap - matrix. */ - prow = mat + iel; - -/* Loop round all the elements in the current row of the matrix. */ - for( icol = 0; icol < ncol; icol++ ) { - -/* If this element is non-zero and non-bad, then output "irow" depends on - input "icol". */ - el = mat[ iel++ ]; - if( el != 0.0 && el != AST__BAD ) { - -/* Is input "icol" one of the selected inputs? */ - isel = 0; - for( i = 0; i < nin; i++ ) { - if( in[ i ] == icol ) { - isel = 1; - break; - } - } - -/* If so, note that this output depends on selected inputs. Otherwise note - it depends on unselected inputs. */ - if( isel ) { - sel = 1; - } else { - unsel = 1; - } - } - } - -/* If this output depends only on selected inputs, we can include it in - the returned Mapping.*/ - if( sel && !unsel ) { - -/* Store the index of the output within the original MatrixMap. */ - result[ nout ] = irow; - -/* Increment the number of outputs in the returned Mapping. */ - nout++; - -/* Copy the elements of the current matrix row which correspond to the - selected inputs into the new matrix. */ - for( i = 0; i < nin; i++ ) { - if( astISGOOD( prow[ in[ i ] ] ) ) { - *(pmat++) = prow[ in[ i ] ]; - good = 1; - } - } - } - -/* If this output depends on a selected input, but also depends on an - unselected input, we cannot split the MatrixMap. */ - if( sel && unsel ) { - ok = 0; - break; - } - } - } - - -/* If the returned Mapping can be created, create it. */ - if( ok && nout > 0 && good ) { - *map = (AstMapping *) astMatrixMap( nin, nout, 0, rmat, "", status ); - -/* Otherwise, free the returned array. */ - } else { - result = astFree( result ); - } - -/* Free resources. */ - rmat = astFree( rmat ); - -/* Re-compress the supplied MatrixMap. */ - CompressMatrix( this, status ); - } - } - -/* Free returned resources if an error has occurred. */ - if( !astOK ) { - result = astFree( result ); - *map = astAnnul( *map ); - } - -/* Return the list of output indices. */ - return result; -} - -static AstMatrixMap *MatMat( AstMapping *map1, AstMapping *map2, int inv1, - int inv2, int *status ){ -/* -* Name: -* MatMat - -* Purpose: -* Create a merged MatrixMap from two supplied MatrixMaps. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* AstMatrixMap *MatMat( AstMapping *map1, AstMapping *map2, int inv1, -* int inv2, int *status ) - -* Class Membership: -* MatrixMap member function - -* Description: -* This function creates a new MatrixMap which performs a mapping -* equivalent to applying the two supplied MatrixMaps in series, in the -* directions specified by the "invert" flags (the Invert attributes of -* the supplied MatrixMaps are ignored). - -* Parameters: -* map1 -* A pointer to the MatrixMap to apply first. -* map2 -* A pointer to the MatrixMap to apply second. -* inv1 -* The invert flag to use with map1. A value of zero causes the forward -* mapping to be used, and a non-zero value causes the inverse -* mapping to be used. -* inv2 -* The invert flag to use with map2. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* Pointer to the new MatrixMap. - -* Notes: -* - The forward direction of the returned MatrixMap is equivalent to the -* combined effect of the two supplied MatrixMap, operating in the -* directions specified by "inv1" and "inv2". -* - A null pointer will be returned if this function is invoked with the -* global error status set, or if it should fail for any reason. -*/ - -/* Local Variables: */ - AstMatrixMap *result; /* Pointer to output MatrixMap */ - int invert[ 2 ]; /* Original invert flags */ - -/* Check the global error status. */ - if ( !astOK ) return NULL; - -/* Initialise the returned pointer. */ - result = NULL; - -/* Temporarily set their Invert attributes to the supplied values. */ - invert[ 0 ] = astGetInvert( map1 ); - astSetInvert( map1, inv1 ); - - invert[ 1 ] = astGetInvert( map2 ); - astSetInvert( map2, inv2 ); - -/* Create a new MatrixMap by multiplying them together. */ - result = astMtrMult( (AstMatrixMap *) map1, (AstMatrixMap *) map2 ); - -/* Re-instate the original settings of the Invert attributes for the - supplied MatrixMaps. */ - astSetInvert( map1, invert[ 0 ] ); - astSetInvert( map2, invert[ 1 ] ); - -/* If an error has occurred, annull the returned MatrixMap. */ - if( !astOK ) result = astAnnul( result ); - -/* Return a pointer to the output MatrixMap. */ - return result; -} - -static AstMatrixMap *MatPerm( AstMatrixMap *mm, AstPermMap *pm, int minv, - int pinv, int mat1, int *status ){ -/* -* Name: -* MatPerm - -* Purpose: -* Create a MatrixMap by merging a MatrixMap and a PermMap. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* AstMatrixMap *MatPerm( AstMatrixMap *mm, AstPermMap *pm, int minv, -* int pinv, int mat1, int *status ) - -* Class Membership: -* MatrixMap member function - -* Description: -* This function creates a new MatrixMap which performs a mapping -* equivalent to applying the two supplied Mappings in series in the -* directions specified by the "invert" flags (the Invert attributes of -* the supplied MatrixMaps are ignored). - -* Parameters: -* mm -* A pointer to the MatrixMap. -* pm -* A pointer to the PermMap. -* minv -* The invert flag to use with mm. A value of zero causes the forward -* mapping to be used, and a non-zero value causes the inverse -* mapping to be used. -* pinv -* The invert flag to use with pm. -* mat1 -* If non-zero, then "mm" is applied first followed by "pm". Otherwise, -* "pm" is applied first followed by "mm". -* status -* Pointer to the inherited status variable. - -* Returned Value: -* Pointer to the new MatrixMap. - -* Notes: -* - The forward direction of the returned MatrixMap is equivalent to the -* combined effect of the two supplied Mappings, operating in the -* directions specified by "pinv" and "minv". -* - A null pointer will be returned if this function is invoked with the -* global error status set, or if it should fail for any reason. -*/ - -/* Local Variables: */ - AstMatrixMap *mm2; /* Pointer to an intermediate MatrixMap */ - AstMatrixMap *result; /* Pointer to output MatrixMap */ - AstPointSet *pset1; /* Pointer to a PointSet holding unpermuted unit vectors */ - AstPointSet *pset2; /* Pointer to a PointSet holding permuted unit vectors */ - double *matrix; /* Pointer to a matrix representing the PermMap */ - double *p; /* Pointer to next matrix element */ - double **ptr1; /* Pointer to the data in pset1 */ - double **ptr2; /* Pointer to the data in pset2 */ - int i; /* Axis index */ - int j; /* Point index */ - int nax; /* No. of axes in the PermMap */ - int old_minv; /* Original setting of MatrixMap Invert attribute */ - int old_pinv; /* Original setting of PermMap Invert attribute */ - -/* Check the global error status. */ - if ( !astOK ) return NULL; - -/* Initialise the returned pointer. */ - result = NULL; - -/* Temporarily set the Invert attributes of both Mappings to the supplied - values. */ - old_minv = astGetInvert( mm ); - astSetInvert( mm, minv ); - - old_pinv = astGetInvert( pm ); - astSetInvert( pm, pinv ); - -/* Get the number of axes in the PermMap. The PermMap will have the same - number of input and output axes because a check has already been made on - it to ensure that this is so (in function PermOK). */ - nax = astGetNin( pm ); - -/* We first represent the PermMap as a MatrixMap containing elements with - values zero or one. Each row of this matrix is obtained by transforming a - unit vector along each axis using the inverse PermMap. Allocate memory - to hold the matrix array, and create a PointSet holding the unit - vectors. */ - matrix = (double *) astMalloc( sizeof( double )*(size_t)( nax*nax ) ); - - pset1 = astPointSet( nax, nax, "", status ); - ptr1 = astGetPoints( pset1 ); - - pset2 = astPointSet( nax, nax, "", status ); - ptr2 = astGetPoints( pset2 ); - - if( astOK ){ - for( i = 0; i < nax; i++ ){ - for( j = 0; j < nax; j++ ) ptr1[ i ][ j ] = 0.0; - ptr1[ i ][ i ] = 1.0; - } - -/* Transform these unit vectors using the inverse PermMap. */ - (void) astTransform( pm, pset1, 0, pset2 ); - -/* Copy the transformed vectors into the matrix array. */ - p = matrix; - for( j = 0; j < nax; j++ ){ - for( i = 0; i < nax; i++ ) *(p++) = ptr2[ i ][ j ]; - } - -/* Create a MatrixMap holding this array. */ - mm2 = astMatrixMap( nax, nax, 0, matrix, "", status ); - -/* Create a new MatrixMap equal to the product of the supplied MatrixMap - and the MatrixMap just created from the PermMap. */ - if( mat1 ){ - result = astMtrMult( mm, mm2 ); - } else { - result = astMtrMult( mm2, mm ); - } - -/* Free everything. */ - mm2 = astAnnul( mm2 ) ; - } - - pset2 = astAnnul( pset2 ); - pset1 = astAnnul( pset1 ); - matrix = (double *) astFree( (void *) matrix ); - -/* Re-instate the original settings of the Invert attribute for the - supplied Mappings. */ - astSetInvert( mm, old_minv ); - astSetInvert( pm, old_pinv ); - -/* If an error has occurred, annull the returned MatrixMap. */ - if( !astOK ) result = astAnnul( result ); - -/* Return a pointer to the output MatrixMap. */ - return result; -} - -static void MatPermSwap( AstMapping **maps, int *inverts, int imm, int *status ){ -/* -* Name: -* MatPermSwap - -* Purpose: -* Swap a PermMap and a MatrixMap. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* void MatPermSwap( AstMapping **maps, int *inverts, int imm ) - -* Class Membership: -* MatrixMap member function - -* Description: -* A list of two Mappings is supplied containing a PermMap and a -* MatrixMap. These Mappings are annulled, and replaced with -* another pair of Mappings consisting of a PermMap and a MatrixMap -* in the opposite order. These Mappings are chosen so that their -* combined effect is the same as the original pair of Mappings. - -* Parameters: -* maps -* A pointer to an array of two Mapping pointers. -* inverts -* A pointer to an array of two invert flags. -* imm -* The index within "maps" of the MatrixMap. - -* Notes: -* - There are restictions on the sorts of PermMaps which can be -* swapped with a MatrixMap -- see function CanSwap. It is assumed -* that the supplied MatrixMap and PermMap satisfy these requirements. - -*/ - -/* Local Variables: */ - AstMatrixMap *mm; /* Pointer to the supplied MatrixMap */ - AstMatrixMap *mmnew; /* Pointer to new MatrixMap */ - AstMatrixMap *smmnew; /* Pointer to new simplified MatrixMap */ - AstPermMap *pm; /* Pointer to the supplied PermMap */ - AstPermMap *pmnew; /* Pointer to new PermMap */ - AstPermMap *spmnew; /* Pointer to new simplified PermMap */ - double *consts; /* Pointer to constants array */ - double *matrix; /* Supplied array of matrix elements */ - double *out_el; /* Pointer to next element of new MatrixMap */ - double *out_mat; /* Matrix elements for new MatrixMap */ - double c; /* Constant */ - double matel; /* Matrix element */ - int *inperm; /* Pointer to input axis permutation array */ - int *outperm; /* Pointer to output axis permutation array */ - int col; /* Index of matrix column */ - int i; /* Axis count */ - int k; /* Axis count */ - int nin; /* No. of axes in supplied PermMap */ - int nout; /* No. of axes in returned PermMap */ - int old_pinv; /* Invert value for the supplied PermMap */ - int row; /* Index of matrix row */ - -/* Check the global error status. */ - if ( !astOK ) return; - -/* Initialise variables to avoid "used of uninitialised variable" - messages from dumb compilers. */ - mmnew = NULL; - pmnew = NULL; - -/* Store pointers to the supplied PermMap and the MatrixMap. */ - pm = (AstPermMap *) maps[ 1 - imm ]; - mm = (AstMatrixMap *) maps[ imm ]; - -/* Temporarily set the Invert attribute of the supplied PermMap to the - supplied value. */ - old_pinv = astGetInvert( pm ); - astSetInvert( pm, inverts[ 1 - imm ] ); - -/* Ensure the MatrixMap is stored in full form. */ - ExpandMatrix( mm, status ); - -/* Store a pointer to the required array of matrix elements. */ - if( inverts[ imm ] ) { - matrix = mm->i_matrix; - } else { - matrix = mm->f_matrix; - } - -/* Get the number of input and output axes of the PermMap. */ - nin = astGetNin( pm ); - nout = astGetNout( pm ); - -/* Allocate memory to hold the matrix elements for the swapped MatrixMap. - The number of rows and olumns in the new matrix must equal the number of - input or output axes for the PermMap, depending on whether the PermMap - or MatrixMap is applied first. */ - if( imm == 0 ) { - out_mat = (double *) astMalloc( sizeof( double )*(size_t)( nout*nout ) ); - } else { - out_mat = (double *) astMalloc( sizeof( double )*(size_t)( nin*nin ) ); - } - -/* We need to know the axis permutation arrays and constants array for - the PermMap. */ - PermGet( pm, &outperm, &inperm, &consts, status ); - if( astOK ) { - -/* First deal with cases where the MatrixMap is applied first. */ - if( imm == 0 ) { - -/* Consider each output axis of the PermMap. */ - for( i = 0; i < nout; i++ ) { - -/* If this output is connected to one of the input axes... */ - row = outperm[ i ]; - if( row >= 0 && row < nin ) { - -/* Permute the row of the supplied matrix which feeds the corresponding - PermMap input axis (i.e. axis outperm[k] ) using the forward PermMap. - Store zeros for any output axes which are assigned constants. This forms - row i of the new MatrixMap. */ - out_el = out_mat + nout*i; - for( k = 0; k < nout; k++ ){ - col = outperm[ k ]; - if( col >= 0 && col < nin ) { - *(out_el++) = *( matrix + nin*row + col ); - } else { - *(out_el++) = 0.0; - } - } - -/* If this output is asigned a constant value, use a "diagonal" vector for - row i of the new MatrixMap (i.e. all zeros except for a 1.0 in column - i ). */ - } else { - out_el = out_mat + nout*i; - for( k = 0; k < nout; k++ ) { - if( k != i ) { - *(out_el++) = 0.0; - } else { - *(out_el++) = 1.0; - } - } - } - } - -/* Create the new MatrixMap. */ - mmnew = astMatrixMap( nout, nout, 0, out_mat, "", status ); - -/* Any PermMap inputs which are assigned a constant value need to be - changed now, since they will no longer be scaled by the inverse - MatrixMap. CanSwap ensures that the inverse MatrixMap produces a - simple scaling for constant axes, so we change the PermMap constant - to be the constant AFTER scaling by the inverse MatrixMap. - - Consider each input axis of the PermMap. */ - for( i = 0; i < nin; i++ ) { - -/* If this input is assigned a constant value... */ - if( inperm[ i ] < 0 ) { - -/* Divide the supplied constant value by the corresponding diagonal term - in the supplied MatrixMap. */ - c = consts[ -inperm[ i ] - 1 ]; - if( c != AST__BAD ) { - matel = matrix[ i*( nin + 1 ) ]; - if( matel != 0.0 && matel != AST__BAD ) { - consts[ -inperm[ i ] - 1 ] /= matel; - } else { - consts[ -inperm[ i ] - 1 ] = AST__BAD; - } - } - } - } - -/* Now deal with cases where the PermMap is applied first. */ - } else { - -/* Consider each input axis of the PermMap. */ - for( i = 0; i < nin; i++ ) { - -/* If this input is connected to one of the output axes... */ - row = inperm[ i ]; - if( row >= 0 && row < nout ) { - -/* Permute the row of the supplied matrix which feeds the corresponding - PermMap output axis (i.e. axis inperm[k] ) using the inverse PermMap. - Store zeros for any input axes which are assigned constants. This forms - row i of the new MatrixMap. */ - out_el = out_mat + nin*i; - for( k = 0; k < nin; k++ ){ - col = inperm[ k ]; - if( col >= 0 && col < nout ) { - *(out_el++) = *( matrix + nout*row + col ); - } else { - *(out_el++) = 0.0; - } - } - -/* If this input is asigned a constant value, use a "diagonal" vector for - row i of the new MatrixMap (i.e. all zeros except for a 1.0 in column - i ). */ - } else { - out_el = out_mat + nin*i; - for( k = 0; k < nin; k++ ) { - if( k != i ) { - *(out_el++) = 0.0; - } else { - *(out_el++) = 1.0; - } - } - } - } - -/* Create the new MatrixMap. */ - mmnew = astMatrixMap( nin, nin, 0, out_mat, "", status ); - -/* Any PermMap outputs which are assigned a constant value need to be - changed now, since they will no longer be scaled by the forward - MatrixMap. CanSwap ensures that the forward MatrixMap produces a - simple scaling for constant axes, so we change the PermMap constant - to be the constant AFTER scaling by the forward MatrixMap. - - Consider each output axis of the PermMap. */ - for( i = 0; i < nout; i++ ) { - -/* If this output is assigned a constant value... */ - if( outperm[ i ] < 0 ) { - -/* Multiple the supplied constant value by the corresponding diagonal term in - the supplied MatrixMap. */ - c = consts[ -outperm[ i ] - 1 ]; - if( c != AST__BAD ) { - matel = matrix[ i*( nout + 1 ) ]; - if( matel != AST__BAD ) { - consts[ -outperm[ i ] - 1 ] *= matel; - } else { - consts[ -outperm[ i ] - 1 ] = AST__BAD; - } - } - } - } - } - -/* Create a new PermMap (since the constants may have changed). */ - pmnew = astPermMap( nin, inperm, nout, outperm, consts, "", status ); - -/* Free the axis permutation and constants arrays. */ - outperm = (int *) astFree( (void *) outperm ); - inperm = (int *) astFree( (void *) inperm ); - consts = (double *) astFree( (void *) consts ); - } - -/* Free the memory used to hold the new matrix elements. */ - out_mat = (double *) astFree( (void *) out_mat ); - -/* Ensure the supplied MatrixMap is stored back in compressed form. */ - CompressMatrix( mm, status ); - -/* Re-instate the original value of the Invert attribute of the supplied - PermMap. */ - astSetInvert( pm, old_pinv ); - - if( astOK ) { - -/* Annul the supplied PermMap. */ - (void) astAnnul( pm ); - -/* Simplify the returned Mappings. */ - spmnew = astSimplify( pmnew ); - pmnew = astAnnul( pmnew ); - - smmnew = astSimplify( mmnew ); - mmnew = astAnnul( mmnew ); - -/* Store a pointer to the new PermMap in place of the supplied MatrixMap. This - PermMap should be used in its forward direction. */ - maps[ imm ] = (AstMapping *) spmnew; - inverts[ imm ] = astGetInvert( spmnew ); - -/* Annul the supplied matrixMap. */ - (void) astAnnul( mm ); - -/* Store a pointer to the new MatrixMap. This MatrixMap should be used in - its forward direction. */ - maps[ 1 - imm ] = (AstMapping *) smmnew; - inverts[ 1 - imm ] = astGetInvert( smmnew ); - } - -/* Return. */ - return; -} - -static void MatWin( AstMapping **maps, int *inverts, int imm, int *status ){ -/* -* Name: -* MatWin - -* Purpose: -* Swap a WinMap and a MatrixMap. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* void MatWin( AstMapping **maps, int *inverts, int imm, int *status ) - -* Class Membership: -* WinMap member function - -* Description: -* A list of two Mappings is supplied containing a WinMap and a -* MatrixMap. These Mappings are annulled, and replaced with -* another pair of Mappings consisting of a WinMap and a MatrixMap -* in the opposite order. These Mappings are chosen so that their -* combined effect is the same as the original pair of Mappings. -* The scale factors in the returned WinMap are always unity (i.e. -* the differences in scaling get absorbed into the returned -* MatrixMap). - -* Parameters: -* maps -* A pointer to an array of two Mapping pointers. -* inverts -* A pointer to an array of two invert flags. -* imm -* The index within "maps" of the MatrixMap. -* status -* Pointer to the inherited status variable. - -*/ - -/* Local Variables: */ - AstMatrixMap *m1; /* Pointer to Diagonal scale factor MatrixMap */ - AstMatrixMap *m2; /* Pointer to returned MatrixMap */ - AstMatrixMap *sm2; /* Pointer to simplified returned MatrixMap */ - AstMatrixMap *mm; /* Pointer to the supplied MatrixMap */ - AstPointSet *pset1; /* Shift terms from supplied WinMap */ - AstPointSet *pset2; /* Shift terms for returned WinMap */ - AstWinMap *w1; /* Pointer to the returned WinMap */ - AstWinMap *sw1; /* Pointer to the simplified returned WinMap */ - AstWinMap *wm; /* Pointer to the supplied WinMap */ - double **ptr1; /* Pointer to pset1 data */ - double **ptr2; /* Pointer to pset2 data */ - double *a; /* Array of shift terms from supplied WinMap */ - double *aa; /* Pointer to next shift term */ - double *b; /* Array of scale terms from supplied WinMap */ - double *bb; /* Pointer to next scale term */ - int i; /* Axis count */ - int nin; /* No. of axes in supplied WinMap */ - int nout; /* No. of axes in returned WinMap */ - int old_minv; /* Invert value for the supplied MatrixMap */ - int old_winv; /* Invert value for the supplied WinMap */ - -/* Check the global error status. */ - if ( !astOK ) return; - -/* Store pointers to the supplied WinMap and the MatrixMap. */ - wm = (AstWinMap *) maps[ 1 - imm ]; - mm = (AstMatrixMap *) maps[ imm ]; - -/* Temporarily set the Invert attribute of the supplied Mappings to the - supplied values. */ - old_winv = astGetInvert( wm ); - astSetInvert( wm, inverts[ 1 - imm ] ); - - old_minv = astGetInvert( mm ); - astSetInvert( mm, inverts[ imm ] ); - -/* Get copies of the shift and scale terms used by the WinMap. This - also returns the number of axes in the WinMap. */ - nin = astWinTerms( wm, &a, &b ); - -/* Create a diagonal MatrixMap holding the scale factors from the - supplied WinMap. */ - m1 = astMatrixMap( nin, nin, 1, b, "", status ); - -/* Create a PointSet holding a single position given by the shift terms - in the supplied WinMap. */ - pset1 = astPointSet( 1, nin, "", status ); - ptr1 = astGetPoints( pset1 ); - if( astOK ){ - aa = a; - for( i = 0; i < nin; i++ ) ptr1[ i ][ 0 ] = *(aa++); - } - -/* First deal with cases when the WinMap is applied first, followed by - the MatrixMap. */ - if( imm == 1 ){ - -/* Multiply the diagonal matrix holding the WinMap scale factors by the - supplied matrix. The resulting MatrixMap is the one to return in the - map list. */ - m2 = astMtrMult( m1, mm ); - -/* Transform the position given by the shift terms from the supplied - WinMap using the supplied MatrixMap to get the shift terms for - the returned WinMap. */ - pset2 = astTransform( mm, pset1, 1, NULL ); - -/* Now deal with cases when the MatrixMap is applied first, followed by - the WinMap. */ - } else { - -/* Multiply the supplied MatrixMap by the diagonal matrix holding scale - factors from the supplied WinMap. The resulting MatrixMap is the one to - return in the map list. */ - m2 = astMtrMult( mm, m1 ); - -/* Transform the position given by the shift terms from the supplied - WinMap using the inverse of the returned MatrixMap to get the shift - terms for the returned WinMap. */ - pset2 = astTransform( m2, pset1, 0, NULL ); - - } - -/* Re-instate the original value of the Invert attributes of the supplied - Mappings. */ - astSetInvert( wm, old_winv ); - astSetInvert( mm, old_minv ); - -/* Get pointers to the shift terms for the returned WinMap. */ - ptr2 = astGetPoints( pset2 ); - -/* Create the returned WinMap, initially with undefined corners. The number of - axes in the WinMap must equal the number of shift terms. */ - nout = astGetNcoord( pset2 ); - w1 = astWinMap( nout, NULL, NULL, NULL, NULL, "", status ); - -/* If succesful, store the scale and shift terms in the WinMap. The scale - terms are always unity. */ - if( astOK ){ - bb = w1->b; - aa = w1->a; - for( i = 0; i < nout; i++ ) { - *(bb++) = 1.0; - *(aa++) = ptr2[ i ][ 0 ]; - } - -/* Replace the supplied Mappings and invert flags with the ones found - above. Remember that the order of the Mappings is now swapped */ - (void) astAnnul( maps[ 0 ] ); - (void) astAnnul( maps[ 1 ] ); - - sw1 = astSimplify( w1 ); - w1 = astAnnul( w1 ); - - maps[ imm ] = (AstMapping *) sw1; - inverts[ imm ] = astGetInvert( sw1 ); - - sm2 = astSimplify( m2 ); - m2 = astAnnul( m2 ); - - maps[ 1 - imm ] = (AstMapping *) sm2; - inverts[ 1 - imm ] = astGetInvert( sm2 ); - - } - -/* Annul the MatrixMap and PointSet holding the scale and shift terms from the - supplied WinMap. */ - m1 = astAnnul( m1 ); - pset1 = astAnnul( pset1 ); - pset2 = astAnnul( pset2 ); - -/* Free the copies of the scale and shift terms from the supplied WinMap. */ - b = (double *) astFree( (void *) b ); - a = (double *) astFree( (void *) a ); - -/* Return. */ - return; -} - -static AstWinMap *MatWin2( AstMatrixMap *mm, AstWinMap *wm, int minv, - int winv, int mat1, int *status ){ -/* -* Name: -* MatWin2 - -* Purpose: -* Create a WinMap by merging a diagonal MatrixMap and a WinMap. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* AstWinMap *MatWin2( AstMatrixMap *mm, AstWinMap *wm, int minv, -* int winv, int mat1, int *status ) - -* Class Membership: -* MatrixMap member function - -* Description: -* This function creates a new WinMap which performs a mapping -* equivalent to applying the two supplied Mappings in series in the -* directions specified by the "invert" flags (the Invert attributes of -* the supplied MatrixMaps are ignored), in the order specified by -* "mat1". - -* Parameters: -* mm -* A pointer to the MatrixMap. Assumed to be diagonal. -* wm -* A pointer to the WinMap. -* minv -* The invert flag to use with mm. A value of zero causes the forward -* mapping to be used, and a non-zero value causes the inverse -* mapping to be used. -* winv -* The invert flag to use with wm. -* mat1 -* If non-zero, then "mm" is applied first followed by "wm". Otherwise, -* "wm" is applied first followed by "mm". -* status -* Pointer to the inherited status variable. - -* Returned Value: -* Pointer to the new MatrixMap. - -* Notes: -* - The forward direction of the returned MatrixMap is equivalent to the -* combined effect of the two supplied Mappings, operating in the -* directions specified by "winv" and "minv". -* - A null pointer will be returned if this function is invoked with the -* global error status set, or if it should fail for any reason. -*/ - -/* Local Variables: */ - AstWinMap *result; /* Pointer to output WinMap */ - double *ina; /* Input corner A in new WinMap */ - double *inb; /* Input corner B in new WinMap */ - double *newscales; /* Scales for new WinMap */ - double *newshifts; /* Shifts for new WinMap */ - double *outa; /* Output corner A in new WinMap */ - double *outb; /* Output corner B in new WinMap */ - double *scales2; /* Pointer to extended WinMap scales array */ - double *scales; /* Pointer to WinMap scales array */ - double *shifts; /* Pointer to WinMap shifts array */ - int i; /* Axis index */ - int ncol; /* No. of columns in the MatrixMap */ - int nrow; /* No. of rows in the MatrixMap */ - int nt; /* Number of axes in WinMap */ - int old_minv; /* Original setting of MatrixMap Invert attribute */ - int old_winv; /* Original setting of WinMap Invert attribute */ - -/* Check the global error status. */ - if ( !astOK ) return NULL; - -/* Initialise the returned pointer. */ - result = NULL; - -/* Temporarily set the Invert attributes of both Mappings to the supplied - values. */ - old_minv = astGetInvert( mm ); - astSetInvert( mm, minv ); - - old_winv = astGetInvert( wm ); - astSetInvert( wm, winv ); - -/* Get the number of inputs (columns) and outputs (rows) for the MatrixMap. */ - ncol = astGetNin( mm ); - nrow = astGetNout( mm ); - -/* Get the scales and shifts implemented by the WinMap. These take into - account the current Invert attribute of the WinMap. */ - nt = astWinTerms( wm, &shifts, &scales ); - -/* First deal with cases where the MatrixMap is applied first. */ - if( mat1 ){ - -/* Sanity check. */ - if( nt != nrow ) { - if( astOK ) astError( AST__INTER, "astMapMerge(%s): WinMap has %d axes, " - "but MatrixMap has %d rows (internal AST programming " - "error).", status, astGetClass(mm), nt, nrow ); - - } else { - -/* Allocate the array to hold the scale terms for the new WinMap. */ - newscales = astMalloc( nrow*sizeof(double) ); - -/* Ensure that the original scales array is padded with sufficient zeros - to allow it to be transformed using the matrixmap. */ - scales2 = astCalloc( ncol, sizeof(double) ); - if( astOK ) memcpy( scales2, scales, - (ncol<nrow?ncol:nrow)*sizeof(double) ); - -/* Use the MatrixMap to transform the scale terms from the WinMap. */ - astTranN( mm, 1, ncol, 1, scales2, 1, nrow, 1, newscales ); - -/* Free resources. */ - scales2 = astFree( scales2 ); - -/* The shifts are unchanged. */ - newshifts = shifts; - } - -/* Now deal with cases where the WinMap is applied first. */ - } else { - -/* Sanity check. */ - if( nt != ncol ) { - if( astOK ) astError( AST__INTER, "astMapMerge(%s): WinMap has %d axes, " - "but MatrixMap has %d columns (internal AST programming " - "error).", status, astGetClass(mm), nt, ncol ); - - } else { - -/* Allocate the array to hold the scale and shift terms for the new WinMap. */ - newscales = astMalloc( nrow*sizeof(double) ); - newshifts = astMalloc( nrow*sizeof(double) ); - -/* Use the MatrixMap to transform the scale terms from the WinMap. */ - astTranN( mm, 1, ncol, 1, scales, 1, nrow, 1, newscales ); - -/* Use the MatrixMap to transform the shift terms from the WinMap. */ - astTranN( mm, 1, ncol, 1, shifts, 1, nrow, 1, newshifts ); - - } - } - -/* Create the new WinMap. */ - ina = astMalloc( nt*sizeof(double) ); - inb = astMalloc( nt*sizeof(double) ); - outa = astMalloc( nt*sizeof(double) ); - outb = astMalloc( nt*sizeof(double) ); - if( astOK ) { - for( i = 0; i < nt; i++ ) { - ina[ i ] = 0.0; - inb[ i ] = 1.0; - outa[ i ] = newshifts[ i ]; - outb[ i ] = newscales[ i ] + newshifts[ i ]; - } - result = astWinMap( nt, ina, inb, outa, outb, "", status ); - } - -/* Re-instate the original settings of the Invert attribute for the - supplied Mappings. */ - astSetInvert( mm, old_minv ); - astSetInvert( wm, old_winv ); - -/* Free resources. */ - ina = astFree( ina ); - inb = astFree( inb ); - outa = astFree( outa ); - outb = astFree( outb ); - if( newscales != scales ) newscales = astFree( newscales ); - if( newshifts != shifts ) newshifts = astFree( newshifts ); - scales = astFree( scales ); - shifts = astFree( shifts ); - -/* If an error has occurred, annull the returned MatrixMap. */ - if( !astOK ) result = astAnnul( result ); - -/* Return a pointer to the output MatrixMap. */ - return result; -} - -static AstMatrixMap *MatZoom( AstMatrixMap *mm, AstZoomMap *zm, int minv, - int zinv, int *status ){ -/* -* Name: -* MatZoom - -* Purpose: -* Create a MatrixMap by merging a MatrixMap and a ZoomMap. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* AstMatrixMap *MatZoom( AstMatrixMap *mm, AstZoomMap *zm, int minv, -* int zinv, int *status ) - -* Class Membership: -* MatrixMap member function - -* Description: -* This function creates a new MatrixMap which performs a mapping -* equivalent to applying the two supplied Mappings in series in the -* directions specified by the "invert" flags (the Invert attributes of -* the supplied MatrixMaps are ignored). - -* Parameters: -* mm -* A pointer to the MatrixMap. -* zm -* A pointer to the ZoomMap. -* minv -* The invert flag to use with mm. A value of zero causes the forward -* mapping to be used, and a non-zero value causes the inverse -* mapping to be used. -* zinv -* The invert flag to use with zm. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* Pointer to the new MatrixMap. - -* Notes: -* - The forward direction of the returned MatrixMap is equivalent to the -* combined effect of the two supplied Mappings, operating in the -* directions specified by "zinv" and "minv". -* - A null pointer will be returned if this function is invoked with the -* global error status set, or if it should fail for any reason. -*/ - -/* Local Variables: */ - AstMatrixMap *mm2; /* Pointer to intermediate MatrixMap */ - AstMatrixMap *result; /* Pointer to output MatrixMap */ - double *matrix; /* Pointer to diagonal matrix elements array */ - double zfac; /* Zoom factor */ - int i; /* Axis index */ - int nrow; /* No. of rows in the MatrixMap */ - int old_minv; /* Original setting of MatrixMap Invert attribute */ - int old_zinv; /* Original setting of ZoomMap Invert attribute */ - -/* Check the global error status. */ - if ( !astOK ) return NULL; - -/* Initialise the returned pointer. */ - result = NULL; - -/* Temporarily set the Invert attributes of both Mappings to the supplied - values. */ - old_minv = astGetInvert( mm ); - astSetInvert( mm, minv ); - - old_zinv = astGetInvert( zm ); - astSetInvert( zm, zinv ); - -/* Get the number of rows in the MatrixMap (i.e. the number of output - axes). */ - nrow = astGetNout( mm ); - -/* Get the zoom factor implemented by the ZoomMap. Invert it if necessary - since astGetZoom does not take account of the Invert setting. */ - zfac = astGetZoom( zm ); - if( zinv ) zfac = 1.0 / zfac; - -/* Create a diagonal matrix map in which each diagonal element is equal - to the zoom factor. */ - matrix = (double *) astMalloc( sizeof( double )*(size_t) nrow ); - if( astOK ) { - for( i = 0; i < nrow; i++ ) matrix[ i ] = zfac; - } - mm2 = astMatrixMap( nrow, nrow, 1, matrix, "", status ); - matrix = (double *) astFree( (void *) matrix ); - -/* Create a new MatrixMap holding the product of the supplied MatrixMap - and the diagonal MatrixMap just created. */ - result = astMtrMult( mm, mm2 ); - mm2 = astAnnul( mm2 ); - -/* Re-instate the original settings of the Invert attribute for the - supplied Mappings. */ - astSetInvert( mm, old_minv ); - astSetInvert( zm, old_zinv ); - -/* If an error has occurred, annull the returned MatrixMap. */ - if( !astOK ) result = astAnnul( result ); - -/* Return a pointer to the output MatrixMap. */ - return result; -} - -/* Functions which access class attributes. */ -/* ---------------------------------------- */ -/* Implement member functions to access the attributes associated with - this class using the macros defined for this purpose in the - "object.h" file. For a description of each attribute, see the class - interface (in the associated .h file). */ - -static AstMatrixMap *MtrMult( AstMatrixMap *this, AstMatrixMap *a, int *status ){ -/* -*+ -* Name: -* astMtrMult - -* Purpose: -* Multiply a MatrixMap by another MatrixMap. - -* Type: -* Protected virtual function. - -* Synopsis: -* #include "matrixmap.h" -* AstMatrixMap *MtrMult( astMatrixMap *this, astMatrixMap *a ) - -* Class Membership: -* MatrixMap method - -* Description: -* This function multiples the matrices given by "this" and "a", returning -* a pointer to a new MatrixMap holding the product "a x this". -* -* The number of columns in the "a" matrix must match the number of -* rows in the "this" matrix. The number of rows in the returned -* MatrixMap is equal to the number of rows in "a", and the number of -* columns is the same as the number of rows in "this". - -* Parameters: -* this -* Pointer to the first MatrixMap. -* a -* Pointer to a second MatrixMap. - -* Returned Value: -* A pointer to the product MatrixMap. - -* Notes: -* - An error is reported if the two MatrixMaps have incompatible -* shapes, or if either MatrixMap does not have a defined forward -* transformation. -* - A null Object pointer will also be returned if this function -* is invoked with the AST error status set, or if it should fail -* for any reason. -*- -*/ - -/* Local variables. */ - astDECLARE_GLOBALS /* Pointer to thread-specific global data */ - AstMatrixMap *new; /* New MatrixMap holding the product matrix */ - double *a_matrix; /* Pointer to the forward "a" matrix */ - double *a_row; /* Pointer to start of current row in "a" */ - double a_val; /* Current element value from "a" */ - double factor; /* Diagonal matrix term */ - double *new_matrix; /* Pointer to the new forward "this" matrix */ - double *new_val; /* Pointer to current output element value */ - double sum; /* Dot product value */ - double *this_col; /* Pointer to start of current column in "this" */ - double *this_matrix; /* Pointer to the forward "this" matrix */ - double this_val; /* Current element value from "this" */ - int col; /* Current output column number */ - int i; /* Loop count */ - int minrow; /* Min. number of rows in "a" or "this" */ - int ncol_a; /* No. of columns in the "a" MatrixMap */ - int ncol_this; /* No. of columns in the "this" MatrixMap */ - int nrow_a; /* No. of rows in the "a" MatrixMap */ - int nrow_this; /* No. of rows in the "this" MatrixMap */ - int row; /* Current output row number */ - -/* Return a NULL pointer if an error has already occurred. */ - if ( !astOK ) return NULL; - -/* Get a pointer to the thread specific global data structure. */ - astGET_GLOBALS(NULL); - -/* Initialise */ - new = NULL; - -/* Report an error if eitherof the MatrixMaps doe snot have a defined - forward transformation.*/ - if( !astGetTranForward( this ) ){ - astError( AST__MTRML, "astMtrMult(%s): Cannot find the product of 2 " - "MatrixMaps- the first MatrixMap has no forward transformation.", status, - astClass(this) ); - return NULL; - } - - if( !astGetTranForward( a ) ){ - astError( AST__MTRML, "astMtrMult(%s): Cannot find the product of 2 " - "MatrixMaps- the second MatrixMap has no forward transformation.", status, - astClass(this) ); - return NULL; - } - -/* Report an error if the shapes of the two matrices are incompatible. */ - nrow_a = astGetNout( a ); - ncol_a = astGetNin( a ); - nrow_this = astGetNout( this ); - ncol_this = astGetNin( this ); - - if( ncol_a != nrow_this && astOK ){ - astError( AST__MTRML, "astMtrMult(%s): Number of rows in the first " - "MatrixMap (%d) does not equal number of columns in the " - "second MatrixMap (%d).", status, astClass(this), nrow_this, ncol_a ); - return NULL; - } - -/* Store the minimum number of rows in either matrix for later use. */ - if( nrow_a < nrow_this ){ - minrow = nrow_a; - } else { - minrow = nrow_this; - } - -/* Ensure that "this" is stored in FULL form (i.e. with all elements - stored explicitly, even if the matrix is a unit or diagonal matrix). */ - ExpandMatrix( this, status ); - -/* Store pointers to the current forward matrices (taking into - account the current states of the Mapping inversion flags ). */ - this_matrix = astGetInvert( this ) ? this->i_matrix : this->f_matrix; - a_matrix = astGetInvert( a ) ? a->i_matrix : a->f_matrix; - -/* Get memory to hold the product matrix in full form. */ - new_matrix = (double *) astMalloc( sizeof( double )* - (size_t)( nrow_a*ncol_this ) ); - if( astOK ){ - -/* First deal with cases where the "a" MatrixMap represents a unit - matrix. */ - if( a->form == UNIT ){ - -/* Copy the required number of rows from "this" to "new". */ - (void) memcpy( (void *) new_matrix, (const void *) this_matrix, - sizeof(double)*(size_t)( minrow*ncol_this ) ); - -/* If there are insufficient rows in "this", append some zero-filled rows. */ - if( minrow < nrow_a ){ - for( i = minrow*ncol_this; i < nrow_a*ncol_this; i++ ){ - new_matrix[ i ] = 0.0; - } - } - -/* Now deal with cases where the "a" MatrixMap represents a diagonal - matrix. */ - } else if( a->form == DIAGONAL ){ - -/* Scale the required number of rows from "this" storing them in "new", - and checking for bad values. */ - i = 0; - - for( row = 0; row < minrow; row++ ){ - factor = a_matrix[ row ]; - - if( factor != AST__BAD ){ - - for( col = 0; col < ncol_this; col++ ){ - this_val = this_matrix[ i ]; - if( this_val != AST__BAD ){ - new_matrix[ i ] = this_val*factor; - } else { - new_matrix[ i ] = AST__BAD; - } - i++; - } - - } else { - - for( col = 0; col < ncol_this; col++ ){ - new_matrix[ i++ ] = AST__BAD; - } - - } - } - -/* If there are insufficient rows in "this", append some zero-filled rows. */ - if( minrow < nrow_a ){ - for( i = minrow*ncol_this; i < nrow_a*ncol_this; i++ ){ - new_matrix[ i ] = 0.0; - } - } - - -/* Now deal with cases where the "a" MatrixMap represents a full, non-diagonal - matrix. */ - } else { - -/* Initialise a pointer to the next element in the product matrix. */ - new_val = new_matrix; - -/* Get a pointer to the start of each row of the "a" matrix. */ - for( row = 0; row < nrow_a; row++ ){ - a_row = a_matrix + ncol_a*row; - -/* Get a pointer to the start of each column of the "this" matrix. */ - for( col = 0; col < ncol_this; col++ ){ - this_col = this_matrix + col; - -/* Form the dot product of the current row from "a", and the current - column from "this", checking for bad values. */ - sum = 0.0; - for( i = 0; i < ncol_a; i++ ){ - a_val = a_row[ i ]; - this_val = this_col[ i*ncol_this ]; - if( a_val != AST__BAD && this_val != AST__BAD ){ - sum += a_val*this_val; - } else { - sum = AST__BAD; - break; - } - } - -/* Store the output matrix element value. */ - *(new_val++) = sum; - - } - } - } - -/* Create the new MatrixMap. */ - new = astInitMatrixMap( NULL, sizeof( AstMatrixMap ), !class_init, - &class_vtab, "MatrixMap", ncol_this, nrow_a, - FULL, new_matrix ); - -/* If possible, compress the new MatrixMap by removing off-diagonal zero - elements. */ - CompressMatrix( new, status ); - -/* Re-compress the original "this" MatrixMap. */ - CompressMatrix( this, status ); - - } - -/* Free the memory used to hold the product matrix in full form. */ - new_matrix = (double *) astFree( (void *) new_matrix ); - - return new; - -} - -static AstMatrixMap *MtrRot( AstMatrixMap *this, double theta, - const double axis[], int *status ){ -/* -*+ -* Name: -* astMtrRot - -* Purpose: -* Multiply a MatrixMap by a rotation matrix. - -* Type: -* Protected virtual function. - -* Synopsis: -* #include "matrixmap.h" -* AstMatrixMap *astMtrRot( astMatrixMap *this, double theta, -* const double axis[] ) - -* Class Membership: -* MatrixMap method. - -* Description: -* This function creates a new MatrixMap which is a copy of "this", -* rotated by a specified angle. It can only be used on MatrixMaps which -* have either 2 or 3 output coordinates. In the 3-D case, the rotation -* is about an arbitrary axis passing through the origin. - -* Parameters: -* this -* Pointer to the MatrixMap. -* theta -* The angle by which to rotate the matrix, in radians. If the matrix -* is applied to a 2-D vector position, the resulting vector is -* rotated clockwise about the origin (i.e. from the positive direction -* of the second axis to the positive direction of the first axis). If -* the vector positions are three dimensional, the rotation is clockwise -* when looking along the vector given by "axis". Note, "theta" measures -f when looking along the vector given by AXIS. Note, THETA measures -* the movemement of the vectors relative to a fixed reference frame. -* Alternatively, the reference frame can be thought of as rotating by -* (-theta) relative to the fixed vectors. -* axis -* A 3-D vector specifying the axis of rotation. This parameter is -* ignored if the output from MatrixMap is 2-dimensional. - -* Returned Value: -* A pointer to the rotated MatrixMap. - -* Notes: -* - A null Object pointer will also be returned if this function -* is invoked with the AST error status set, or if it should fail -* for any reason. -*- -*/ - -/* Local variables. */ - AstMatrixMap *new; /* New MatrixMap holding the rotated matrix */ - double as,a,b,c,d,e,f,g; /* Intermediate quantities */ - double axlen; /* Length of axis vector */ - double axlen2; /* Squared length of axis vector */ - double costh; /* Cos(rotation angle) */ - double sinth; /* Sin(rotation angle) */ - double rotmat[9]; /* Rotation matrix */ - double work[3]; /* Work space for matrix multiplication */ - int ncol; /* No. of columns in the MatrixMap */ - int nrow; /* No. of rows in the MatrixMap */ - -/* Return a NULL pointer if an error has already occurred. */ - if ( !astOK ) return NULL; - -/* Initialise the returned MarixMap to be a copy of the supplied MatrixMap. */ - new = astCopy( this ); - -/* Save the cos and sin of the rotation angle for future use. */ - costh = cos( theta ); - sinth = sin( theta ); - -/* Return without changing the MatrixMap if the rotation angle is a - multiple of 360 degrees. */ - if ( costh == 1.0 ) return new; - -/* Get the dimensions of the MatrixMap. */ - nrow = astGetNout( new ); - ncol = astGetNin( new ); - -/* First do rotation of a plane about the origin. */ - if( nrow == 2 ){ - -/* Ensure that the MatrixMap is stored in full form rather than - compressed form. */ - ExpandMatrix( new, status ); - -/* Form the 2x2 forward rotation matrix. Theta is the clockwise angle - of rotation. */ - rotmat[0] = costh; - rotmat[1] = sinth; - rotmat[2] = -sinth; - rotmat[3] = costh; - -/* Post-multiply the current forward matrix (depending on whether or not - the MatrixMap has been inverted) by the forward rotation matrix. */ - if( !astGetInvert( new ) ){ - SMtrMult( 1, 2, ncol, rotmat, new->f_matrix, work, status ); - } else { - SMtrMult( 1, 2, ncol, rotmat, new->i_matrix, work, status ); - } - -/* Now form the 2x2 inverse rotation matrix (the diagonal elements - don't change). */ - rotmat[1] = -sinth; - rotmat[2] = sinth; - -/* Pre-multiply the current inverse matrix (depending on whether or - not the MatrixMap has been inverted) by the inverse rotation matrix. */ - if( !astGetInvert( new ) ){ - SMtrMult( 0, ncol, 2, rotmat, new->i_matrix, work, status ); - } else { - SMtrMult( 0, ncol, 2, rotmat, new->f_matrix, work, status ); - } - -/* See if the matrix can be stored as a UNIT or DIAGONAL matrix. */ - CompressMatrix( new, status ); - -/* Now do rotation of a volume about an axis passing through the origin. */ - } else if( nrow == 3 ){ - -/* Find the length of the axis vector. Report an error if it has zero - length or has not been supplied. */ - if( axis ) { - axlen2 = axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2]; - } else { - axlen2 = 0.0; - } - if( axlen2 <= 0.0 ) { - astError( AST__MTRAX, "astMtrRot(%s): NULL or zero length " - "axis vector supplied.", status, astClass(new) ); - } - axlen = sqrt( axlen2 ); - -/* Ensure that the MatrixMap is stored in full form rather than - compressed form. */ - ExpandMatrix( new, status ); - -/* Form commonly used terms in the rotation matrix. */ - as = sinth/axlen; - a = (1.0 - costh)/axlen2; - b = a*axis[0]*axis[1]; - c = as*axis[2]; - d = a*axis[0]*axis[2]; - e = as*axis[1]; - f = a*axis[1]*axis[2]; - g = as*axis[0]; - -/* Form the 3x3 forward rotation matrix. Theta is the clockwise angle - of rotation looking in the direction of the axis vector. */ - rotmat[0] = a*axis[0]*axis[0] + costh; - rotmat[1] = b - c; - rotmat[2] = d + e; - rotmat[3] = b + c; - rotmat[4] = a*axis[1]*axis[1] + costh; - rotmat[5] = f - g; - rotmat[6] = d - e; - rotmat[7] = f + g; - rotmat[8] = a*axis[2]*axis[2] + costh; - -/* Post-multiply the current forward matrix (depending on whether or not - the MatrixMap has been inverted) by the forward rotation matrix. */ - if( !astGetInvert( new ) ){ - SMtrMult( 1, 3, ncol, rotmat, new->f_matrix, work, status ); - } else { - SMtrMult( 1, 3, ncol, rotmat, new->i_matrix, work, status ); - } - -/* Now form the 3x3 inverse rotation matrix (the diagonal elements - don't change). */ - rotmat[1] = b + c; - rotmat[2] = d - e; - rotmat[3] = b - c; - rotmat[5] = f + g; - rotmat[6] = d + e; - rotmat[7] = f - g; - -/* Pre-multiply the current inverse matrix (depending on whether or - not the MatrixMap has been inverted) by the inverse rotation matrix. */ - if( !astGetInvert( new ) ){ - SMtrMult( 0, ncol, 3, rotmat, new->i_matrix, work, status ); - } else { - SMtrMult( 0, ncol, 3, rotmat, new->f_matrix, work, status ); - } - -/* See if the matrix can be stored as a UNIT or DIAGONAL matrix. */ - CompressMatrix( new, status ); - -/* Report an error if the matrix is not suitable for rotation. */ - } else { - astError( AST__MTR23, "astMtrRot(%s): Cannot rotate a %dx%d" - " MatrixMap.", status, astClass(new), nrow, ncol ); - } - -/* Delete the new MatrixMap if an error has occurred. */ - if( !astOK ) new = astDelete( new ); - - return new; - -} - -static void PermGet( AstPermMap *map, int **outperm, int **inperm, - double **consts, int *status ){ -/* -* Name: -* PermGet - -* Purpose: -* Get the axis permutation and constants array for a PermMap. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* void PermGet( AstPermMap *map, int **outperm, int **inperm, -* double **const, int *status ) - -* Class Membership: -* MatrixMap member function - -* Description: -* This function returns axis permutation and constants arrays which can -* be used to create a PermMap which is equivalent to the supplied PermMap. - -* Parameters: -* map -* The PermMap. -* outperm -* An address at which to return a popinter to an array of ints -* holding the output axis permutation array. The array should be -* released using astFree when no longer needed. -* inperm -* An address at which to return a popinter to an array of ints -* holding the input axis permutation array. The array should be -* released using astFree when no longer needed. -* consts -* An address at which to return a popinter to an array of doubles -* holding the constants array. The array should be released using -* astFree when no longer needed. -* status -* Pointer to the inherited status variable. - -* Notes: -* - NULL pointers are returned if an error has already occurred, or if -* this function should fail for any reason. -*/ - -/* Local Variables: */ - AstPointSet *pset1; /* PointSet holding input positions for PermMap */ - AstPointSet *pset2; /* PointSet holding output positions for PermMap */ - double **ptr1; /* Pointer to pset1 data */ - double **ptr2; /* Pointer to pset2 data */ - double *cnst; /* Pointer to constants array */ - double cn; /* Potential new constant value */ - double ip; /* Potential output axis index */ - double op; /* Potential input axis index */ - int *inprm; /* Pointer to input axis permutation array */ - int *outprm; /* Pointer to output axis permutation array */ - int i; /* Axis count */ - int nc; /* Number of constants stored so far */ - int nin; /* No. of input coordinates for the PermMap */ - int nout; /* No. of output coordinates for the PermMap */ - -/* Initialise. */ - if( outperm ) *outperm = NULL; - if( inperm ) *inperm = NULL; - if( consts ) *consts = NULL; - -/* Check the global error status and the supplied pointers. */ - if ( !astOK || !outperm || !inperm || !consts ) return; - -/* Get the number of input and output axes for the supplied PermMap. */ - nin = astGetNin( map ); - nout = astGetNout( map ); - -/* Allocate the memory for the returned arrays. */ - outprm = (int *) astMalloc( sizeof( int )* (size_t) nout ); - inprm = (int *) astMalloc( sizeof( int )* (size_t) nin ); - cnst = (double *) astMalloc( sizeof( double )* (size_t) ( nout + nin ) ); - -/* Returned the pointers to these arrays.*/ - *outperm = outprm; - *inperm = inprm; - *consts = cnst; - -/* Create two PointSets, each holding two points, which can be used for - input and output positions with the PermMap. */ - pset1 = astPointSet( 2, nin, "", status ); - pset2 = astPointSet( 2, nout, "", status ); - -/* Set up the two input positions to be [0,1,2...] and [-1,-1,-1,...]. The - first position is used to enumerate the axes, and the second is used to - check for constant axis values. */ - ptr1 = astGetPoints( pset1 ); - if( astOK ){ - for( i = 0; i < nin; i++ ){ - ptr1[ i ][ 0 ] = ( double ) i; - ptr1[ i ][ 1 ] = -1.0; - } - } - -/* Use the PermMap to transform these positions in the forward direction. */ - (void) astTransform( map, pset1, 1, pset2 ); - -/* No constant axis valeus found yet. */ - nc = 0; - -/* Look at the mapped positions to determine the output axis permutation - array. */ - ptr2 = astGetPoints( pset2 ); - if( astOK ){ - -/* Do each output axis. */ - for( i = 0; i < nout; i++ ){ - -/* If the output axis value is copied from an input axis value, the index - of the appropriate input axis will be in the mapped first position. */ - op = ptr2[ i ][ 0 ]; - -/* If the output axis value is assigned a constant value, the result of - mapping the two different input axis values will be the same. */ - cn = ptr2[ i ][ 1 ]; - if( op == cn ) { - -/* We have found another constant. Store it in the constants array, and - store the index of the constant in the output axis permutation array. */ - cnst[ nc ] = cn; - outprm[ i ] = -( nc + 1 ); - nc++; - -/* If the output axis values are different, then the output axis value - must be copied from the input axis value. */ - } else { - outprm[ i ] = (int) ( op + 0.5 ); - } - } - } - -/* Now do the same thing to determine the input permutation array. */ - if( astOK ){ - for( i = 0; i < nout; i++ ){ - ptr2[ i ][ 0 ] = ( double ) i; - ptr2[ i ][ 1 ] = -1.0; - } - } - - (void) astTransform( map, pset2, 0, pset1 ); - - if( astOK ){ - - for( i = 0; i < nin; i++ ){ - - ip = ptr1[ i ][ 0 ]; - cn = ptr1[ i ][ 1 ]; - if( ip == cn ) { - - cnst[ nc ] = cn; - inprm[ i ] = -( nc + 1 ); - nc++; - - } else { - inprm[ i ] = (int) ( ip + 0.5 ); - } - } - } - -/* Annul the PointSets. */ - pset1 = astAnnul( pset1 ); - pset2 = astAnnul( pset2 ); - -/* If an error has occurred, attempt to free the returned arrays. */ - if( !astOK ) { - *outperm = (int *) astFree( (void *) *outperm ); - *inperm = (int *) astFree( (void *) *inperm ); - *consts = (double *) astFree( (void *) *consts ); - } - -/* Return. */ - return; -} - -static int PermOK( AstMapping *pm, int *status ){ -/* -* Name: -* PermOK - -* Purpose: -* Determine if a PermMap can be merged with a MatrixMap. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* int PermOK( AstMapping *pm, int *status ) - -* Class Membership: -* PermMap member function - -* Description: -* This function returns a flag indicating if the supplied PermMap -* could be merged with a MatrixMap. For thios to be possible, the -* PermMap must have the same number of input and output axes, and the -* inverse and forward mappings must be consistent. - -* Parameters: -* pm -* The PermMap. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* 1 if the PermMap can be merged, 0 otherwise. - -* Notes: -* - A value of 0 is returned if an error has already occurred, or if -* this function should fail for any reason. -*/ - -/* Local Variables: */ - AstPointSet *pset1; /* PointSet holding input positions for PermMap */ - AstPointSet *pset2; /* PointSet holding output positions for PermMap */ - double **ptr1; /* Pointer to pset1 data */ - int i; /* Loop count */ - int nin; /* No. of input coordinates for the PermMap */ - int nout; /* No. of output coordinates for the PermMap */ - int ret; /* Returned flag */ - -/* Check the global error status. */ - if ( !astOK ) return 0; - -/* Initialise */ - ret = 0; - -/* The PermMap must have the same number of input and output coordinates. */ - nin = astGetNin( pm ); - nout = astGetNout( pm ); - if( nin == nout ){ - -/* Create two PointSets, each holding two points, which can be used for - the input and output positions with the PermMap. */ - pset1 = astPointSet( 2, nin, "", status ); - pset2 = astPointSet( 2, nout, "", status ); - -/* Set up the two input positions to be [1,2,3...] and [0,-1,-2,...] */ - ptr1 = astGetPoints( pset1 ); - if( astOK ){ - for( i = 0; i < nin; i++ ){ - ptr1[ i ][ 0 ] = ( double )( i + 1 ); - ptr1[ i ][ 1 ] = ( double )( -i ); - } - - } - -/* Use the PermMap to transform these positions in the forward direction. */ - (void) astTransform( pm, pset1, 1, pset2 ); - -/* Now transform the results back again using the inverse PermMap. */ - (void) astTransform( pm, pset2, 0, pset1 ); - -/* See if the input positions have changed. If they have, then the PermMap - does not have a consistent pair of transformations. If they have not, - then the transformations must be consistent because we used two - different input positions and only one could come out unchanged by - chance. */ - if( astOK ){ - ret = 1; - for( i = 0; i < nin; i++ ){ - if( ptr1[ i ][ 0 ] != ( double )( i + 1 ) || - ptr1[ i ][ 1 ] != ( double )( -i ) ){ - ret = 0; - break; - } - } - } - -/* Annul the PointSets. */ - pset1 = astAnnul( pset1 ); - pset2 = astAnnul( pset2 ); - } - -/* Return the answer. */ - return astOK ? ret : 0; -} - -static double Rate( AstMapping *this, double *at, int ax1, int ax2, int *status ){ -/* -* Name: -* Rate - -* Purpose: -* Calculate the rate of change of a Mapping output. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* result = Rate( AstMapping *this, double *at, int ax1, int ax2, int *status ) - -* Class Membership: -* MatrixMap member function (overrides the astRate method inherited -* from the Mapping class ). - -* Description: -* This function returns the rate of change of a specified output of -* the supplied Mapping with respect to a specified input, at a -* specified input position. - -* Parameters: -* this -* Pointer to the Mapping to be applied. -* at -* The address of an array holding the axis values at the position -* at which the rate of change is to be evaluated. The number of -* elements in this array should equal the number of inputs to the -* Mapping. -* ax1 -* The index of the Mapping output for which the rate of change is to -* be found (output numbering starts at 0 for the first output). -* ax2 -* The index of the Mapping input which is to be varied in order to -* find the rate of change (input numbering starts at 0 for the first -* input). -* status -* Pointer to the inherited status variable. - -* Returned Value: -* The rate of change of Mapping output "ax1" with respect to input -* "ax2", evaluated at "at", or AST__BAD if the value cannot be -* calculated. - -*/ - -/* Local Variables: */ - AstMatrixMap *map; - double *matrix; - double result; - -/* Check inherited status */ - if( !astOK ) return AST__BAD; - -/* Get a pointer to the MatrixMap structure. */ - map = (AstMatrixMap *) this; - -/* Get a pointer to the array holding the required matrix elements, according - to whether the MatrixMap has been inverted. */ - if( !astGetInvert( this ) ) { - matrix = map->f_matrix; - } else { - matrix = map->i_matrix; - } - -/* First deal with full MatrixMaps in which all matrix elements are stored. */ - if( map->form == FULL ){ - result = matrix[ ax1*astGetNin( this ) + ax2 ]; - -/* For unit matrices, the rate is unity if the input and output axes are - equal, and zero otherwise. */ - } else if( map->form == UNIT ){ - result = (ax1 == ax2 ) ? 1.0 : 0.0; - -/* For diagonal matrices, the rate is zero for off diagonal elements and - the matrix array stored the on-diagonal rates. */ - } else if( ax1 == ax2 ) { - result = matrix[ ax1 ]; - - } else { - result = 0.0; - } - -/* Return the result. */ - return result; -} - -static void SMtrMult( int post, int m, int n, const double *mat1, - double *mat2, double *work, int *status ){ -/* -* Name: -* SMtrMult - -* Purpose: -* Multiply a square matrix and a non-square matrix. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* void SMtrMult( int post, int m, int n, const double *mat1, -* double *mat2, double *work, int *status ) - -* Class Membership: -* MatrixMap member function. - -* Description: -* The matrix pointed to by "mat2" is modified by multiplying it by -* the square matrix pointed to by "mat1". If "post" is 1, then: -* -* mat2 -> mat1*mat2 (mat1 is mxm and mat2 is mxn) -* -* If "post" is zero, then: -* -* mat2 -> mat2*mat1 (mat1 is nxn and mat2 is mxn) -* -* The restriction that "mat1" must be square is imposed so that the -* returned matrix will have the same shape as the supplied matrix (mat1). - -* Parameters: -* post -* Specifies whether to post- or pre- multiply mat2 by mat1. -* m -* The number of rows in mat2. -* n -* The number of columns in mat2. -* mat1 -* The multiplier matrix. It must be square of size m or n, depending -* on "post". -* mat2 -* The multiplicand matrix. -* work -* Pointer to work space containing room for m doubles (if "post" -* is 1), or n doubles (if "post" is 0). -* status -* Pointer to the inherited status variable. - -* Notes: -* - No error is reported if "mat2" is supplied NULL. In this case -* it will also be returned NULL. -*/ - -/* Local Variables */ - double dot; /* Output matrix element value */ - const double *mat1_col; /* Pointer to start of current column of mat1 */ - const double *mat1_row; /* Pointer to start of current row of mat1 */ - double *mat2_col; /* Pointer to start of current column of mat2 */ - double *mat2_row; /* Pointer to start of current row of mat2 */ - double cel; /* Column element value */ - double rel; /* Row element value */ - int i; /* Index of current output row */ - int j; /* Index of current output column */ - int k; /* Dot product index */ - -/* Do nothing if mat2 is NULL */ - if ( mat2 ){ - -/* First deal with cases where the supplied matrix is post-multiplied - (i.e. the returned matrix is mat1*mat2). */ - if( post ){ - -/* Loop round each column of the output matrix, storing a pointer to - the start of the corresponding column of mat2. */ - for( j=0; j<n; j++ ){ - mat2_col = mat2 + j; - -/* Loop round each row of the output matrix, storing a pointer to - the start of the corresponding row of mat1. */ - for( i=0; i<m; i++ ){ - mat1_row = mat1 + i*m; - -/* Get the dot product of the corresponding row from mat1 and the - corresponding column from mat2 and store it in the work array. */ - dot = 0.0; - for( k=0; k<m; k++ ) { - rel = mat1_row[ k ]; - cel = mat2_col[ k*n ]; - if( rel != AST__BAD && cel != AST__BAD ){ - dot += rel*cel; - } else { - dot = AST__BAD; - break; - } - } - work[ i ] = dot; - } - -/* Copy the values stored in the work array to the current column of - the output matrix. */ - for( i=0; i<m; i++ ) mat2_col[ i*n ] = work[ i ]; - } - -/* Now deal with cases where the supplied matrix is pre-multiplied - (i.e. the returned matrix is mat2*mat1). */ - } else { - -/* Loop round each row of the output matrix, storing a pointer to - the start of the corresponding row of mat2. */ - for( i=0; i<m; i++ ){ - mat2_row = mat2 + i*n; - -/* Loop round each column of the output matrix, storing a pointer to - the start of the corresponding column of mat1. */ - for( j=0; j<n; j++ ){ - mat1_col = mat1 + j; - -/* Get the dot product of the corresponding row from mat2 and the - corresponding column from mat1 and store it in the work array. */ - dot = 0.0; - for( k=0; k<n; k++ ) { - rel = mat2_row[ k ]; - cel = mat1_col[ k*n ]; - if( rel != AST__BAD && cel != AST__BAD ){ - dot += rel*cel; - } else { - dot = AST__BAD; - break; - } - } - work[ j ] = dot; - } - -/* Copy the values stored in the work array to the current row of - the output matrix. */ - for( j=0; j<n; j++ ) mat2_row[ j ] = work[ j ]; - } - } - } - - return; - -} - -static int GetTranForward( AstMapping *this, int *status ) { -/* -* -* Name: -* GetTranForward - -* Purpose: -* Determine if a MatrixMap defines a forward coordinate transformation. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* int GetTranForward( AstMapping *this, int *status ) - -* Class Membership: -* MatrixMap member function (over-rides the astGetTranForward method -* inherited from the Mapping class). - -* Description: -* This function returns a value indicating if the MatrixMap is able -* to perform a forward coordinate transformation. - -* Parameters: -* this -* Pointer to the MatrixMap. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* Zero if the forward coordinate transformation is not defined, or 1 if it -* is. - -* Notes: -* - A value of zero will be returned if this function is invoked with the -* global error status set, or if it should fail for any reason. -*/ - -/* Local Variables: */ - AstMatrixMap *map; /* Pointer to MatrixMap to be queried */ - int invert; /* Has the mapping been inverted? */ - int result; /* The returned value */ - -/* Initialise. */ - result = 0; - -/* Check the global error status. */ - if ( !astOK ) return result; - -/* Obtain a pointer to the MatrixMap. */ - map = (AstMatrixMap *) this; - -/* All unit MatrixMaps are defined in both directions. */ - if( map->form == UNIT ) { - result = 1; - -/* Otherwise, check that the appropriate array is defined in the - MatrixMap. */ - } else { - -/* Determine if the Mapping has been inverted. */ - invert = astGetInvert( this ); - -/* If OK, obtain the result. */ - if ( astOK ) { - - if( invert ){ - result = ( map->i_matrix != NULL ); - } else { - result = ( map->f_matrix != NULL ); - } - - } - - } - -/* Return the result. */ - return result; - -} - -static int GetTranInverse( AstMapping *this, int *status ) { -/* -* -* Name: -* GetTranInverse - -* Purpose: -* Determine if a MatrixMap defines an inverse coordinate transformation. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* int GetTranInverse( AstMapping *this, int *status ) - -* Class Membership: -* MatrixMap member function (over-rides the astGetTranInverse method -* inherited from the Mapping class). - -* Description: -* This function returns a value indicating if the MatrixMap is able -* to perform an inverse coordinate transformation. - -* Parameters: -* this -* Pointer to the MatrixMap. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* Zero if the inverse coordinate transformation is not defined, or 1 if it -* is. - -* Notes: -* - A value of zero will be returned if this function is invoked with the -* global error status set, or if it should fail for any reason. -*/ - -/* Local Variables: */ - AstMatrixMap *map; /* Pointer to MatrixMap to be queried */ - int invert; /* Has the mapping been inverted? */ - int result; /* The returned value */ - -/* Initialise. */ - result = 0; - -/* Check the global error status. */ - if ( !astOK ) return result; - -/* Obtain a pointer to the MatrixMap. */ - map = (AstMatrixMap *) this; - -/* All unit MatrixMaps are defined in both directions. */ - if( map->form == UNIT ) { - result = 1; - -/* Otherwise, check that the appropriate array is defined in the - MatrixMap. */ - } else { - -/* Determine if the Mapping has been inverted. */ - invert = astGetInvert( this ); - -/* If OK, obtain the result. */ - if ( astOK ) { - - if( invert ){ - result = ( map->f_matrix != NULL ); - } else { - result = ( map->i_matrix != NULL ); - } - - } - - } - -/* Return the result. */ - return result; - -} - -static AstPointSet *Transform( AstMapping *this, AstPointSet *in, - int forward, AstPointSet *out, int *status ) { -/* -* Name: -* Transform - -* Purpose: -* Apply a MatrixMap to transform a set of points. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* AstPointSet *Transform( AstMapping *this, AstPointSet *in, -* int forward, AstPointSet *out, int *status ) - -* Class Membership: -* MatrixMap member function (over-rides the astTransform protected -* method inherited from the Mapping class). - -* Description: -* This function takes a MatrixMap and a set of points encapsulated in a -* PointSet and transforms the points by multiplying them by the matrix. - -* Parameters: -* this -* Pointer to the MatrixMap. -* in -* Pointer to the PointSet holding the input coordinate data. -* forward -* A non-zero value indicates that the forward coordinate transformation -* should be applied, while a zero value requests the inverse -* transformation. -* out -* Pointer to a PointSet which will hold the transformed (output) -* coordinate values. A NULL value may also be given, in which case a -* new PointSet will be created by this function. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* Pointer to the output (possibly new) PointSet. - -* Notes: -* - A null pointer will be returned if this function is invoked with the -* global error status set, or if it should fail for any reason. -* - The number of coordinate values per point in the input PointSet must -* match the number of columns in the MatrixMap being applied. -* - The number of coordinate values per point in the output PointSet will -* equal the number of rows in the MatrixMap being applied. -* - If an output PointSet is supplied, it must have space for sufficient -* number of points and coordinate values per point to accommodate the -* result. Any excess space will be ignored. -*/ - -/* Local Variables: */ - AstPointSet *result; /* Pointer to output PointSet */ - AstMatrixMap *map; /* Pointer to MatrixMap to be applied */ - double diag_term; /* Current diagonal element value */ - double *indata; /* Pointer to next input data value */ - double *matrix; /* Pointer to start of matrix element array */ - double *matrix_element; /* Pointer to current matrix element value */ - double *outdata; /* Pointer to next output data value */ - double **ptr_in; /* Pointer to input coordinate data */ - double **ptr_out; /* Pointer to output coordinate data */ - double sum; /* Partial output value */ - double val; /* Data value */ - int in_coord; /* Index of output coordinate */ - int nax; /* Output axes for which input axes exist */ - int ncoord_in; /* Number of coordinates per input point */ - int ncoord_out; /* Number of coordinates per output point */ - int npoint; /* Number of points */ - int out_coord; /* Index of output coordinate */ - int point; /* Loop counter for points */ - -/* Check the global error status. */ - if ( !astOK ) return NULL; - -/* Obtain a pointer to the MatrixMap. */ - map = (AstMatrixMap *) this; - -/* Apply the parent mapping using the stored pointer to the Transform member - function inherited from the parent Mapping class. This function validates - all arguments and generates an output PointSet if necessary, but does not - actually transform any coordinate values. */ - result = (*parent_transform)( this, in, forward, out, status ); - -/* We will now extend the parent astTransform method by performing the - calculations needed to generate the output coordinate values. */ - -/* Determine the numbers of points and coordinates per point from the input - and output PointSets and obtain pointers for accessing the input and - output coordinate values. */ - ncoord_in = astGetNcoord( in ); - ncoord_out = astGetNcoord( result ); - npoint = astGetNpoint( in ); - ptr_in = astGetPoints( in ); - ptr_out = astGetPoints( result ); - -/* Determine whether to apply the forward or inverse mapping, according to the - direction specified and whether the mapping has been inverted. */ - if ( astGetInvert( map ) ) forward = !forward; - -/* Get a pointer to the array holding the required matrix elements, according - to the direction of mapping required. */ - if ( forward ) { - matrix = map->f_matrix; - } else { - matrix = map->i_matrix; - } - -/* Perform coordinate arithmetic. */ -/* ------------------------------ */ - if ( astOK ) { - -/* First deal with full MatrixMaps in which all matrix elements are stored. */ - if( map->form == FULL ){ - -/* Loop to apply the matrix to each point in turn, checking for - (and propagating) bad values in the process. The matrix elements are - accessed sequentially in row order. The next matrix element to be - used is identified by a pointer which is initialised to point to the - first element of the matrix prior to processing each point. */ - for ( point = 0; point < npoint; point++ ) { - matrix_element = matrix; - -/* Each output co-ordinate value is created by summing the product of the - corresponding input co-ordinates and the elements of one row of the - matrix. */ - for ( out_coord = 0; out_coord < ncoord_out; out_coord++ ) { - sum = 0.0; - - for ( in_coord = 0; in_coord < ncoord_in; in_coord++ ) { - -/* Check the current input coordinate value and the current matrix element. - If the coordinate value is bad, then the output value will also be - bad unless the matrix element is zero. That is, a zero matrix element - results in the input coordinate value being ignored, even if it is bad. - This prevents bad input values being propagated to output axes which - are independant of the bad input axis. A bad matrix element always results - in the output value being bad. In either of these cases, break out of the - loop, remembering to advance the pointer to the next matrix element so - that it points to the start of the next row ready for doing the next - output coordinate. */ - if ( ( ptr_in[ in_coord ][ point ] == AST__BAD && - (*matrix_element) != 0.0 ) || - (*matrix_element) == AST__BAD ) { - sum = AST__BAD; - matrix_element += ncoord_in - in_coord; - break; - -/* If the input coordinate and the current matrix element are both - valid, increment the sum by their product, and step to the next matrix - element pointer If we arrive here with a bad input value, then the - matrix element must be zero, in which case the running sum is left - unchanged. */ - } else { - if ( ptr_in[ in_coord ][ point ] != AST__BAD ) { - sum += ptr_in[ in_coord ][ point ] * (*matrix_element); - } - matrix_element++; - } - } - -/* Store the output coordinate value. */ - ptr_out[ out_coord ][ point ] = sum; - - } - - } - -/* Now deal with unit and diagonal MatrixMaps. */ - } else { - -/* Find the number of output axes for which input data is available. */ - if( ncoord_in < ncoord_out ){ - nax = ncoord_in; - } else { - nax = ncoord_out; - } - -/* For unit matrices, copy the input axes to the corresponding output axes. */ - if( map->form == UNIT ){ - for( out_coord = 0; out_coord < nax; out_coord++ ) { - (void) memcpy( ptr_out[ out_coord ], - (const void *) ptr_in[ out_coord ], - sizeof( double )*(size_t)npoint ); - } - -/* For diagonal matrices, scale each input axis using the appropriate - diagonal element from the matrix, and store in the output. */ - } else { - for( out_coord = 0; out_coord < nax; out_coord++ ){ - diag_term = matrix[ out_coord ]; - outdata = ptr_out[ out_coord ]; - indata = ptr_in[ out_coord ]; - - if( diag_term != AST__BAD ){ - for( point = 0; point < npoint; point++ ){ - val = *(indata++); - if( val != AST__BAD ){ - *(outdata++) = diag_term*val; - } else { - *(outdata++) = AST__BAD; - } - } - - } else { - for( point = 0; point < npoint; point++ ){ - *(outdata++) = AST__BAD; - } - } - } - } - -/* If there are any remaining output axes, fill the first one with zeros. */ - if( nax < ncoord_out ){ - outdata = ptr_out[ nax ]; - for( point = 0; point < npoint; point++ ) *(outdata++) = 0.0; - -/* Copy this axis to any remaining output axes. */ - outdata = ptr_out[ nax ]; - for( out_coord = nax + 1; out_coord < ncoord_out; out_coord++ ) { - (void) memcpy( ptr_out[ out_coord ], (const void *) outdata, - sizeof( double )*(size_t)npoint ); - } - } - } - } - -/* Return a pointer to the output PointSet. */ - return result; -} - -static int ScalingRowCol( AstMatrixMap *map, int axis, int *status ){ -/* -* Name: -* ScalingRowCol - -* Purpose: -* Determine if a given row and column of a MatrixMap are zeros -* with a non-zero diagonal term. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* int ScalingRowCol( AstMatrixMap *map, int axis, int *status ) - -* Class Membership: -* MatrixMap member function - -* Description: -* This function returns a flag indicating if a MatrixMap presents a -* simple scaling for a given axis in both directions. The MatrixMap -* must be square. A value of one is returned if every element of the -* row and column corresponding to the given axis is zero, except for -* the diagonal term which must be non-zero. - -* Parameters: -* map -* The MatrixMap. -* axis -* The zero-based index of the axis to check. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* 1 if the row/column produces a simple scaling, 0 otherwise. - -*/ - -/* Local Variables: */ - double *el; /* Pointer to matrix element */ - int i; /* Element count */ - int ncol; /* No. of input coordinates */ - int ret; /* Returned flag */ - -/* Initialise */ - ret = 0; - -/* Check the global error status. */ - if ( !astOK ) return ret; - -/* If a unit or diagonal MatrixMap has been supplied, return 1. */ - if( map->form != FULL ){ - ret = 1; - -/* If a full matrix has been supplied... */ - } else { - -/* Assume the row/column gives a unit mapping. */ - ret = 1; - -/* Get the number of input axes for the MatrixMap. */ - ncol = astGetNin( map ); - -/* Check that all elements of the "axis"th row are effectively zero, except - for the "axis"th element which must be non-zero. */ - el = map->f_matrix + axis*ncol; - for( i = 0; i < ncol; i++ ) { - if( i == axis ) { - if( fabs( *el ) <= DBL_EPSILON ) { - ret = 0; - break; - } - } else if( fabs( *el ) > DBL_EPSILON ) { - ret = 0; - break; - } - el++; - } - -/* Check that all elements of the "axis"th column are effectively zero, except - for the "axis"th element which must be non-zero. */ - if( ret ) { - el = map->f_matrix + axis; - for( i = 0; i < ncol; i++ ) { - if( i == axis ) { - if( fabs( *el ) <= DBL_EPSILON ) { - ret = 0; - break; - } - } else if( fabs( *el ) > DBL_EPSILON ) { - ret = 0; - break; - } - el += ncol; - } - } - } - -/* Return the answer. */ - return astOK ? ret : 0; -} - -/* Functions which access class attributes. */ -/* ---------------------------------------- */ -/* Implement member functions to access the attributes associated with - this class using the macros defined for this purpose in the - "object.h" file. For a description of each attribute, see the class - interface (in the associated .h file). */ - -/* Copy constructor. */ -/* ----------------- */ -static void Copy( const AstObject *objin, AstObject *objout, int *status ) { -/* -* Name: -* Copy - -* Purpose: -* Copy constructor for MatrixMap objects. - -* Type: -* Private function. - -* Synopsis: -* void Copy( const AstObject *objin, AstObject *objout, int *status ) - -* Description: -* This function implements the copy constructor for MatrixMap objects. - -* Parameters: -* objin -* Pointer to the object to be copied. -* objout -* Pointer to the object being constructed. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* void - -* Notes: -* - This constructor makes a deep copy, including a copy of the matrix -* element values associated with the input MatrixMap. -*/ - - -/* Local Variables: */ - AstMatrixMap *in; /* Pointer to input MatrixMap */ - AstMatrixMap *out; /* Pointer to output MatrixMap */ - int nel; /* No. of elements in the matrix */ - int nin; /* No. of input coordinates */ - int nout; /* No. of output coordinates */ - -/* Check the global error status. */ - if ( !astOK ) return; - -/* Obtain pointers to the input and output MatrixMaps. */ - in = (AstMatrixMap *) objin; - out = (AstMatrixMap *) objout; - -/* Nullify the pointers stored in the output object since these will - currently be pointing at the input data (since the output is a simple - byte-for-byte copy of the input). Otherwise, the input data could be - freed by accidient if the output object is deleted due to an error - occuring in this function. */ - out->f_matrix = NULL; - out->i_matrix = NULL; - -/* If the input MatrixMap is a unit mapping, then no matrix elements are - stored with it, so do nothing in this case. */ - if( out->form != UNIT ){ - -/* Obtain the number of stored values in the MatrixMap. This is independant of - whether the Mapping has been inverted or not. If the MatrixMap is diagonal, - only the diagonal terms are stored. */ - nin = astGetNin( in ); - nout = astGetNout( in ); - - if( out->form == DIAGONAL ){ - if( nin < nout ){ - nel = nin; - } else { - nel = nout; - } - - } else { - nel = nin*nout; - } - -/* Store the forward matrix elements in the output MatrixMap. */ - out->f_matrix = (double *) astStore( NULL, (void *) in->f_matrix, - sizeof( double )*(size_t) nel ); - -/* Store the inverse matrix elements (if defined) in the output - MatrixMap. */ - if( in->i_matrix ){ - out->i_matrix = (double *) astStore( NULL, (void *) in->i_matrix, - sizeof( double )*(size_t) nel ); - } - -/* If an error has occurred, free the output MatrixMap arrays. */ - if( !astOK ) { - out->f_matrix = (double *) astFree( (void *) out->f_matrix ); - out->i_matrix = (double *) astFree( (void *) out->i_matrix ); - } - } - - return; - -} - -/* Destructor. */ -/* ----------- */ -static void Delete( AstObject *obj, int *status ) { -/* -* Name: -* Delete - -* Purpose: -* Destructor for MatrixMap objects. - -* Type: -* Private function. - -* Synopsis: -* void Delete( AstObject *obj, int *status ) - -* Description: -* This function implements the destructor for MatrixMap objects. - -* Parameters: -* obj -* Pointer to the object to be deleted. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* void - -* Notes: -* This function attempts to execute even if the global error status is -* set. -*/ - -/* Local Variables: */ - AstMatrixMap *this; /* Pointer to MatrixMap */ - -/* Obtain a pointer to the MatrixMap structure. */ - this = (AstMatrixMap *) obj; - -/* Free the arrays used to store element values for forward and inverse - matrices. */ - this->f_matrix = (double *) astFree( (void *) this->f_matrix ); - this->i_matrix = (double *) astFree( (void *) this->i_matrix ); -} - -/* Dump function. */ -/* -------------- */ -static void Dump( AstObject *this_object, AstChannel *channel, int *status ) { -/* -* Name: -* Dump - -* Purpose: -* Dump function for MatrixMap objects. - -* Type: -* Private function. - -* Synopsis: -* void Dump( AstObject *this, AstChannel *channel, int *status ) - -* Description: -* This function implements the Dump function which writes out data -* for the MatrixMap class to an output Channel. - -* Parameters: -* this -* Pointer to the MatrixMap whose data are being written. -* channel -* Pointer to the Channel to which the data are being written. -* status -* Pointer to the inherited status variable. -*/ - -#define KEY_LEN 50 /* Maximum length of a keyword */ - -/* Local Variables: */ - AstMatrixMap *this; /* Pointer to the MatrixMap structure */ - char buff[ KEY_LEN + 1 ]; /* Buffer for keyword string */ - int el; /* Element index */ - int nel; /* No. of elements in the matrix */ - int nin; /* No. of input coords */ - int nout; /* No. of output coords */ - -/* Check the global error status. */ - if ( !astOK ) return; - -/* Obtain a pointer to the MatrixMap structure. */ - this = (AstMatrixMap *) this_object; - -/* Find the number of elements stored for each matrix. */ - nin = astGetNin( this ); - nout = astGetNout( this ); - - if( this->form == FULL ){ - nel = nin*nout; - - } else if( this->form == DIAGONAL ){ - nel = astMIN( nin, nout ); - - } else { - nel = 0; - } - -/* Write out values representing the instance variables for the - MatrixMap class. */ - -/* The forward matrix. Note BAD values are not written out as the - AST__BAD value may differ on different machines. If a matrix element - is not found when reading the matrix back in again (in astLoadMatrixMap), - then it is assigned a default value of AST__BAD. */ - if( this->f_matrix ){ - for( el = 0; el < nel; el++ ){ - if( (this->f_matrix)[ el ] != AST__BAD ) { - (void) sprintf( buff, "M%d", el ); - astWriteDouble( channel, buff, 1, 1, (this->f_matrix)[ el ], - "Forward matrix value" ); - } - } - } - -/* The inverse matrix. */ - if( this->i_matrix ){ - for( el = 0; el < nel; el++ ){ - if( (this->i_matrix)[ el ] != AST__BAD ) { - (void) sprintf( buff, "IM%d", el ); - astWriteDouble( channel, buff, 1, 1, (this->i_matrix)[ el ], - "Inverse matrix value" ); - } - } - } - -/* The matrix storage form. */ - astWriteString( channel, "Form", 1, 1, Form[ this->form ], - "Matrix storage form" ); - -/* Undefine macros local to this function. */ -#undef KEY_LEN -} - -/* Standard class functions. */ -/* ========================= */ -/* Implement the astIsAMatrixMap and astCheckMatrixMap functions using the macros - defined for this purpose in the "object.h" header file. */ -astMAKE_ISA(MatrixMap,Mapping) -astMAKE_CHECK(MatrixMap) - -AstMatrixMap *astMatrixMap_( int nin, int nout, int form, - const double matrix[], const char *options, int *status, ...){ -/* -*++ -* Name: -c astMatrixMap -f AST_MATRIXMAP - -* Purpose: -* Create a MatrixMap. - -* Type: -* Public function. - -* Synopsis: -c #include "matrixmap.h" -c AstMatrixMap *astMatrixMap( int nin, int nout, int form, -c const double matrix[], -c const char *options, ... ) -f RESULT = AST_MATRIXMAP( NIN, NOUT, FORM, MATRIX, OPTIONS, STATUS ) - -* Class Membership: -* MatrixMap constructor. - -* Description: -* This function creates a new MatrixMap and optionally initialises -* its attributes. -* -* A MatrixMap is a form of Mapping which performs a general linear -* transformation. Each set of input coordinates, regarded as a -* column-vector, are pre-multiplied by a matrix (whose elements -* are specified when the MatrixMap is created) to give a new -* column-vector containing the output coordinates. If appropriate, -* the inverse transformation may also be performed. - -* Parameters: -c nin -f NIN = INTEGER (Given) -* The number of input coordinates, which determines the number -* of columns in the matrix. -c nout -f NOUT = INTEGER (Given) -* The number of output coordinates, which determines the number -* of rows in the matrix. -c form -f FORM = INTEGER (Given) -* An integer which indicates the form in which the matrix -* elements will be supplied. -* -c A value of zero indicates that a full "nout" x "nin" matrix -f A value of zero indicates that a full NOUT x NIN matrix -c of values will be supplied via the "matrix" parameter -f of values will be supplied via the MATRIX argument -* (below). In this case, the elements should be given in row -* order (the elements of the first row, followed by the -* elements of the second row, etc.). -* -* A value of 1 indicates that only the diagonal elements of the -* matrix will be supplied, and that all others should be -c zero. In this case, the elements of "matrix" should contain -f zero. In this case, the elements of MATRIX should contain -* only the diagonal elements, stored consecutively. -* -* A value of 2 indicates that a "unit" matrix is required, -* whose diagonal elements are set to unity (with all other -c elements zero). In this case, the "matrix" parameter is -c ignored and a NULL pointer may be supplied. -f elements zero). In this case, the MATRIX argument is not used. -c matrix -f MATRIX( * ) = DOUBLE PRECISION (Given) -* The array of matrix elements to be used, stored according to -c the value of "form". -f the value of FORM. -c options -f OPTIONS = CHARACTER * ( * ) (Given) -c Pointer to a null-terminated string containing an optional -c comma-separated list of attribute assignments to be used for -c initialising the new MatrixMap. The syntax used is identical to -c that for the astSet function and may include "printf" format -c specifiers identified by "%" symbols in the normal way. -f A character string containing an optional comma-separated -f list of attribute assignments to be used for initialising the -f new MatrixMap. The syntax used is identical to that for the -f AST_SET routine. -c ... -c If the "options" string contains "%" format specifiers, then -c an optional list of additional arguments may follow it in -c order to supply values to be substituted for these -c specifiers. The rules for supplying these are identical to -c those for the astSet function (and for the C "printf" -c function). -f STATUS = INTEGER (Given and Returned) -f The global status. - -* Returned Value: -c astMatrixMap() -f AST_MATRIXMAP = INTEGER -* A pointer to the new MatrixMap. - -* Notes: -* - In general, a MatrixMap's forward transformation will always -* be available (as indicated by its TranForward attribute), but -* its inverse transformation (TranInverse attribute) will only be -* available if the associated matrix is square and non-singular. -* - As an exception to this, the inverse transformation is always -* available if a unit or diagonal matrix is specified. In this -* case, if the matrix is not square, one or more of the input -* coordinate values may not be recoverable from a set of output -* coordinates. Any coordinates affected in this way will simply be -* set to the value zero. -* - A null Object pointer (AST__NULL) will be returned if this -c function is invoked with the AST error status set, or if it -f function is invoked with STATUS set to an error value, or if it -* should fail for any reason. - -* Status Handling: -* The protected interface to this function includes an extra -* parameter at the end of the parameter list descirbed above. This -* parameter is a pointer to the integer inherited status -* variable: "int *status". - -*-- -*/ - -/* Local Variables: */ - astDECLARE_GLOBALS /* Pointer to thread-specific global data */ - AstMatrixMap *new; /* Pointer to new MatrixMap */ - va_list args; /* Variable argument list */ - -/* Check the global status. */ - if ( !astOK ) return NULL; - -/* Get a pointer to the thread specific global data structure. */ - astGET_GLOBALS(NULL); - -/* Initialise the MatrixMap, allocating memory and initialising the - virtual function table as well if necessary. */ - new = astInitMatrixMap( NULL, sizeof( AstMatrixMap ), !class_init, - &class_vtab, "MatrixMap", nin, nout, form, matrix); - -/* If successful, note that the virtual function table has been - initialised. */ - if ( astOK ) { - class_init = 1; - -/* Obtain the variable argument list and pass it along with the options string - to the astVSet method to initialise the new MatrixMap's attributes. */ - va_start( args, status ); - astVSet( new, options, NULL, args ); - va_end( args ); - -/* If an error occurred, clean up by deleting the new object. */ - if ( !astOK ) new = astDelete( new ); - } - -/* Return a pointer to the new MatrixMap. */ - return new; -} - -AstMatrixMap *astMatrixMapId_( int nin, int nout, int form, const double matrix[], - const char *options, ... ) { -/* -* Name: -* astMatrixMapId_ - -* Purpose: -* Create a MatrixMap. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* AstMatrixMap *astMatrixMapId_( int nin, int nout, int form, -* const double matrix[], const char *options, -* ... ) - -* Class Membership: -* MatrixMap constructor. - -* Description: -* This function implements the external (public) interface to the -* astMatrixMap constructor function. It returns an ID value (instead -* of a true C pointer) to external users, and must be provided -* because astMatrixMap_ has a variable argument list which cannot be -* encapsulated in a macro (where this conversion would otherwise -* occur). -* -* The variable argument list also prevents this function from -* invoking astMatrixMap_ directly, so it must be a re-implementation -* of it in all respects, except for the final conversion of the -* result to an ID value. - -* Parameters: -* As for astMatrixMap_. - -* Returned Value: -* The ID value associated with the new MatrixMap. -*/ - -/* Local Variables: */ - astDECLARE_GLOBALS /* Pointer to thread-specific global data */ - AstMatrixMap *new; /* Pointer to new MatrixMap */ - va_list args; /* Variable argument list */ - int *status; /* Pointer to inherited status value */ - -/* Get a pointer to the inherited status value. */ - status = astGetStatusPtr; - -/* Get a pointer to the thread specific global data structure. */ - astGET_GLOBALS(NULL); - -/* Check the global status. */ - if ( !astOK ) return NULL; - -/* Initialise the MatrixMap, allocating memory and initialising the - virtual function table as well if necessary. */ - new = astInitMatrixMap( NULL, sizeof( AstMatrixMap ), !class_init, &class_vtab, - "MatrixMap", nin, nout, form, matrix ); - -/* If successful, note that the virtual function table has been - initialised. */ - if ( astOK ) { - class_init = 1; - -/* Obtain the variable argument list and pass it along with the options string - to the astVSet method to initialise the new MatrixMap's attributes. */ - va_start( args, options ); - astVSet( new, options, NULL, args ); - va_end( args ); - -/* If an error occurred, clean up by deleting the new object. */ - if ( !astOK ) new = astDelete( new ); - } - -/* Return an ID value for the new MatrixMap. */ - return astMakeId( new ); -} - -AstMatrixMap *astInitMatrixMap_( void *mem, size_t size, int init, - AstMatrixMapVtab *vtab, const char *name, - int nin, int nout, int form, - const double *matrix, int *status ) { -/* -*+ -* Name: -* astInitMatrixMap - -* Purpose: -* Initialise a MatrixMap. - -* Type: -* Protected function. - -* Synopsis: -* #include "matrixmap.h" -* AstMatrixMap *astInitMatrixMap( void *mem, size_t size, int init, -* AstMatrixMapVtab *vtab, const char *name, -* int nin, int nout, int form, -* const double *matrix ) - -* Class Membership: -* MatrixMap initialiser. - -* Description: -* This function is provided for use by class implementations to initialise -* a new MatrixMap object. It allocates memory (if necessary) to accommodate -* the MatrixMap plus any additional data associated with the derived class. -* It then initialises a MatrixMap structure at the start of this memory. If -* the "init" flag is set, it also initialises the contents of a virtual -* function table for a MatrixMap at the start of the memory passed via the -* "vtab" parameter. - -* Parameters: -* mem -* A pointer to the memory in which the MatrixMap is to be initialised. -* This must be of sufficient size to accommodate the MatrixMap data -* (sizeof(MatrixMap)) plus any data used by the derived class. If a value -* of NULL is given, this function will allocate the memory itself using -* the "size" parameter to determine its size. -* size -* The amount of memory used by the MatrixMap (plus derived class data). -* This will be used to allocate memory if a value of NULL is given for -* the "mem" parameter. This value is also stored in the MatrixMap -* structure, so a valid value must be supplied even if not required for -* allocating memory. -* init -* A logical flag indicating if the MatrixMap's virtual function table is -* to be initialised. If this value is non-zero, the virtual function -* table will be initialised by this function. -* vtab -* Pointer to the start of the virtual function table to be associated -* with the new MatrixMap. -* name -* Pointer to a constant null-terminated character string which contains -* the name of the class to which the new object belongs (it is this -* pointer value that will subsequently be returned by the astGetClass -* method). -* nin -* The number of input coordinate values per point. This is the -* same as the number of columns in the matrix. -* nout -* The number of output coordinate values per point. This is the -* same as the number of rows in the matrix. -* form -* If "form" is 2 or larger, then a unit MatrixMap is created. In this -* case "matrix" is ignored and can be supplied as NULL. If "form" is -* 1, then a diagonal MatrixMap is created. In this case, the number of -* values in "matrix" should be equal to the minimum of nin and nout, -* and "matrix" should contain the corresponding diagonal terms, in row -* order. If "form" is 0 or less, then a full MatrixMap is created, and -* "matrix" should contain all nin*nout element values. -* matrix -* A pointer to an array of matrix element values. The values should be -* supplied in row order. The content of this array is determined by -* "form". If a full MatrixMap is to be created then the array starts -* with (row 1, column 1), then comes (row 1, column 2), up to (row 1, -* column nin), then (row 2, column 1), (row 2, column 2), and so on, -* finishing with (row nout, column nin) ). An error is reported if a -* NULL value is supplied unless "form" is 2 or more. - -* Returned Value: -* A pointer to the new MatrixMap. - -* Notes: -* - A null pointer will be returned if this function is invoked with the -* global error status set, or if it should fail for any reason. -*- -*/ - -/* Local Variables: */ - AstMatrixMap *new; /* Pointer to new MatrixMap */ - double *fmat; /* Pointer to the forward matrix */ - double *imat; /* Pointer to the inverse matrix */ - int i; /* Loop count */ - int nel; /* No. of elements in matrix array */ - int nuse; /* Number of usable matrix elements */ - int used_form; /* Form limited to 0, 1 or 2 */ - -/* Check the global status. */ - if ( !astOK ) return NULL; - -/* If necessary, initialise the virtual function table. */ - if ( init ) astInitMatrixMapVtab( vtab, name ); - -/* Initialise. */ - new = NULL; - -/* Report an error if a NULL matrix was supplied, unless a unit MatrixMap - has been requested. */ - if( form < 2 && !matrix ){ - astError( AST__MTRMT, "astInitMatrixMap(%s): NULL matrix supplied.", status, - name ); - - } else { - -/* Initialise a Mapping structure (the parent class) as the first component - within the MatrixMap structure, allocating memory if necessary. Specify that - the Mapping should be defined in both the forward and inverse directions. */ - new = (AstMatrixMap *) astInitMapping( mem, size, 0, - (AstMappingVtab *) vtab, name, - nin, nout, 1, 1 ); - - if ( astOK ) { - -/* Initialise the MatrixMap data. */ -/* ---------------------------- */ -/* If a unit MatrixMap is being created, then no additional storage is - required. */ - if( form > 1 ){ - nel = 0; - used_form = UNIT; - -/* If a diagonal MatrixMap is being created, then memory is needed to hold - the diagonal terms. */ - } else if( form == 1 ){ - if( nin < nout ){ - nel = nin; - } else { - nel = nout; - } - used_form = DIAGONAL; - -/* If a full MatrixMap is being created, then memory is needed to hold - all the terms. */ - } else { - nel = nin*nout ; - used_form = FULL; - } - -/* Allocate memory for the forward matrix, storing the supplied matrix - values in it. */ - fmat = (double *) astStore( NULL, (void *) matrix, - sizeof(double)*(size_t)nel ); - -/* Replace any NaNs by AST__BAD and count the number of usable values. */ - if( nel > 0 ) { - nuse = 0; - for( i = 0; i < nel; i++ ) { - if( !astISFINITE(fmat[ i ]) ) { - fmat[ i ] = AST__BAD; - } else if( fmat[ i ] != AST__BAD ) { - nuse++; - } - } - -/* Report an error if there are no usable values. */ - if( nuse == 0 && astOK ) { - astError( AST__MTRMT, "astInitMatrixMap(%s): Supplied matrix " - "contains only bad values.", status, name ); - } - } - -/* Create an inverse matrix if possible. */ - imat = InvertMatrix( used_form, nout, nin, fmat, status ); - -/* Store the matrix arrays. */ - new->form = used_form; - new->f_matrix = fmat; - new->i_matrix = imat; - -/* Attempt to compress the MatrixMap into DIAGONAL or UNIT form. */ - CompressMatrix( new, status ); - -/* If an error occurred, clean up by deleting the new MatrixMap. */ - if ( !astOK ) new = astDelete( new ); - } - } - -/* Return a pointer to the new MatrixMap. */ - return new; -} - -AstMatrixMap *astLoadMatrixMap_( void *mem, size_t size, - AstMatrixMapVtab *vtab, const char *name, - AstChannel *channel, int *status ) { -/* -*+ -* Name: -* astLoadMatrixMap - -* Purpose: -* Load a MatrixMap. - -* Type: -* Protected function. - -* Synopsis: -* #include "matrixmap.h" -* AstMatrixMap *astLoadMatrixMap( void *mem, size_t size, -* AstMatrixMapVtab *vtab, const char *name, -* AstChannel *channel ) - -* Class Membership: -* MatrixMap loader. - -* Description: -* This function is provided to load a new MatrixMap using data read -* from a Channel. It first loads the data used by the parent class -* (which allocates memory if necessary) and then initialises a -* MatrixMap structure in this memory, using data read from the input -* Channel. -* -* If the "init" flag is set, it also initialises the contents of a -* virtual function table for a MatrixMap at the start of the memory -* passed via the "vtab" parameter. - - -* Parameters: -* mem -* A pointer to the memory into which the MatrixMap is to be -* loaded. This must be of sufficient size to accommodate the -* MatrixMap data (sizeof(MatrixMap)) plus any data used by derived -* classes. If a value of NULL is given, this function will -* allocate the memory itself using the "size" parameter to -* determine its size. -* size -* The amount of memory used by the MatrixMap (plus derived class -* data). This will be used to allocate memory if a value of -* NULL is given for the "mem" parameter. This value is also -* stored in the MatrixMap structure, so a valid value must be -* supplied even if not required for allocating memory. -* -* If the "vtab" parameter is NULL, the "size" value is ignored -* and sizeof(AstMatrixMap) is used instead. -* vtab -* Pointer to the start of the virtual function table to be -* associated with the new MatrixMap. If this is NULL, a pointer -* to the (static) virtual function table for the MatrixMap class -* is used instead. -* name -* Pointer to a constant null-terminated character string which -* contains the name of the class to which the new object -* belongs (it is this pointer value that will subsequently be -* returned by the astGetClass method). -* -* If the "vtab" parameter is NULL, the "name" value is ignored -* and a pointer to the string "MatrixMap" is used instead. - -* Returned Value: -* A pointer to the new MatrixMap. - -* Notes: -* - A null pointer will be returned if this function is invoked -* with the global error status set, or if it should fail for any -* reason. -*- -*/ - -#define KEY_LEN 50 /* Maximum length of a keyword */ - - astDECLARE_GLOBALS /* Pointer to thread-specific global data */ -/* Local Variables: */ - AstMatrixMap *new; /* Pointer to the new MatrixMap */ - char buff[ KEY_LEN + 1 ]; /* Buffer for keyword string */ - const char *form; /* String form */ - int def; /* Is the matrix defined? */ - int el; /* Element index */ - int nel; /* No. of elements in the matrix */ - int nin; /* No. of input coords */ - int nout; /* No. of output coords */ - -/* Get a pointer to the thread specific global data structure. */ - astGET_GLOBALS(channel); - -/* Initialise. */ - new = NULL; - -/* Check the global error status. */ - if ( !astOK ) return new; - -/* If a NULL virtual function table has been supplied, then this is - the first loader to be invoked for this MatrixMap. In this case the - MatrixMap belongs to this class, so supply appropriate values to be - passed to the parent class loader (and its parent, etc.). */ - if ( !vtab ) { - size = sizeof( AstMatrixMap ); - vtab = &class_vtab; - name = "MatrixMap"; - -/* If required, initialise the virtual function table for this class. */ - if ( !class_init ) { - astInitMatrixMapVtab( vtab, name ); - class_init = 1; - } - } - -/* Invoke the parent class loader to load data for all the ancestral - classes of the current one, returning a pointer to the resulting - partly-built MatrixMap. */ - new = astLoadMapping( mem, size, (AstMappingVtab *) vtab, name, - channel ); - - if ( astOK ) { - -/* Read input data. */ -/* ================ */ -/* Request the input Channel to read all the input data appropriate to - this class into the internal "values list". */ - astReadClassData( channel, "MatrixMap" ); - -/* Now obtain the Matrix storage form from this list. */ - form = astReadString( channel, "form", Form[FULL] ); - new->form = FindString( 3, Form, form, "the MatrixMap component 'Form'", - "astRead", astGetClass( channel ), status ); - form = astFree( (void *) form ); - -/* Find the number of elements stored for each matrix. */ - nin = astGetNin( (AstMapping *) new ); - nout = astGetNout( (AstMapping *) new ); - - if( new->form == FULL ){ - nel = nin*nout; - - } else if( new->form == DIAGONAL ){ - nel = astMIN( nin, nout ); - - } else { - nel = 0; - } - -/* Allocate memory to hold the forward matrix. */ - new->f_matrix = (double *) astMalloc( sizeof(double)*(size_t)nel ); - -/* Now read the other data items from the list and use them to - initialise the appropriate instance variable(s) for this class. */ - -/* The forward matrix. */ - if( new->f_matrix ){ - def = 0; - - for( el = 0; el < nel; el++ ){ - (void) sprintf( buff, "m%d", el ); - (new->f_matrix)[ el ] = astReadDouble( channel, buff, AST__BAD ); - if( (new->f_matrix)[ el ] != AST__BAD ) def = 1; - } - -/* Store a NULL pointer if no elements of the matrix were found. */ - if( !def ) new->f_matrix = (double *) astFree( (void *) new->f_matrix ); - - } - -/* The inverse matrix. */ - new->i_matrix = (double *) astMalloc( sizeof(double)*(size_t)nel ); - if( new->i_matrix ){ - def = 0; - - for( el = 0; el < nel; el++ ){ - (void) sprintf( buff, "im%d", el ); - (new->i_matrix)[ el ] = astReadDouble( channel, buff, AST__BAD ); - if( (new->i_matrix)[ el ] != AST__BAD ) def = 1; - } - -/* If no elements of the matrix were found, create an inverse matrix if - possible, otherwise store a NULL pointer. Note, prior to AST 8.6.3, the - inverse matrix was not included in the dump - it was always recalculated - using InvertMatrix, but this led to small round-trip errors in cases, - where the original inverse matrix was not created using InvertMatrix - (e.g. was created by astMtrRot). */ - if( !def ) { - new->i_matrix = (double *) astFree( (void *) new->i_matrix ); - if( new->f_matrix ){ - new->i_matrix = InvertMatrix( new->form, nout, nin, new->f_matrix, status ); - } else { - new->i_matrix = NULL; - } - } - } - -/* If an error occurred, clean up by deleting the new MatrixMap. */ - if ( !astOK ) new = astDelete( new ); - } - -/* Return the new MatrixMap pointer. */ - return new; - -/* Undefine macros local to this function. */ -#undef KEY_LEN -} - -/* Virtual function interfaces. */ -/* ============================ */ -/* These provide the external interface to the virtual functions defined by - this class. Each simply checks the global error status and then locates and - executes the appropriate member function, using the function pointer stored - in the object's virtual function table (this pointer is located using the - astMEMBER macro defined in "object.h"). - - Note that the member function may not be the one defined here, as it may - have been over-ridden by a derived class. However, it should still have the - same interface. */ - -AstMatrixMap *astMtrRot_( AstMatrixMap *this, double theta, - const double axis[], int *status ){ - if( !astOK ) return NULL; - return (**astMEMBER(this,MatrixMap,MtrRot))( this, theta, axis, status ); -} - -AstMatrixMap *astMtrMult_( AstMatrixMap *this, AstMatrixMap *a, int *status ){ - if( !astOK ) return NULL; - return (**astMEMBER(this,MatrixMap,MtrMult))( this, a, status ); -} - - - - |