diff options
author | William Joye <wjoye@cfa.harvard.edu> | 2018-01-09 19:26:44 (GMT) |
---|---|---|
committer | William Joye <wjoye@cfa.harvard.edu> | 2018-01-09 19:26:44 (GMT) |
commit | 1332d38f2805d986ea130e43218c0d2e870b4dc1 (patch) | |
tree | aa72853cb8d0d8fcd53a6f5eddf196a374226706 /ast/polymap.c | |
parent | 5e545ec8058cc5238dc870468b34b5d4617f307f (diff) | |
download | blt-1332d38f2805d986ea130e43218c0d2e870b4dc1.zip blt-1332d38f2805d986ea130e43218c0d2e870b4dc1.tar.gz blt-1332d38f2805d986ea130e43218c0d2e870b4dc1.tar.bz2 |
update ast 8.6.2
Diffstat (limited to 'ast/polymap.c')
-rw-r--r-- | ast/polymap.c | 6107 |
1 files changed, 0 insertions, 6107 deletions
diff --git a/ast/polymap.c b/ast/polymap.c deleted file mode 100644 index 498000f..0000000 --- a/ast/polymap.c +++ /dev/null @@ -1,6107 +0,0 @@ -/* -*class++ -* Name: -* PolyMap - -* Purpose: -* Map coordinates using polynomial functions. - -* Constructor Function: -c astPolyMap -f AST_POLYMAP - -* Description: -* A PolyMap is a form of Mapping which performs a general polynomial -* transformation. Each output coordinate is a polynomial function of -* all the input coordinates. The coefficients are specified separately -* for each output coordinate. The forward and inverse transformations -* are defined independantly by separate sets of coefficients. If no -* inverse transformation is supplied, the default behaviour is to use -* an iterative method to evaluate the inverse based only on the forward -* transformation (see attribute IterInverse). - -* Inheritance: -* The PolyMap class inherits from the Mapping class. - -* Attributes: -* In addition to those attributes common to all Mappings, every -* PolyMap also has the following attributes: -* -* - IterInverse: Provide an iterative inverse transformation? -* - NiterInverse: Maximum number of iterations for iterative inverse -* - TolInverse: Target relative error for iterative inverse - -* Functions: -c In addition to those functions applicable to all Objects, the -c following functions may also be applied to all Mappings: -f In addition to those routines applicable to all Objects, the -f following routines may also be applied to all Mappings: -* -c - astPolyCoeffs: Retrieve the coefficients of a PolyMap transformation -c - astPolyTran: Fit a PolyMap inverse or forward transformation -f - AST_POLYCOEFFS: Retrieve the coefficients of a PolyMap transformation -f - AST_POLYTRAN: Fit a PolyMap inverse or forward transformation - -* Copyright: -* Copyright (C) 2017 East Asian Observatory. -* Copyright (C) 1997-2006 Council for the Central Laboratory of the -* Research Councils -* Copyright (C) 2011 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) - -* History: -* 27-SEP-2003 (DSB): -* Original version. -* 13-APR-2005 (DSB): -* Changed the keys used by the Dump/astLoadPolyMap functions. They -* used to exceed 8 characters and consequently caused problems for -* FitsChans. -* 20-MAY-2005 (DSB): -* Correct the indexing of keywords produced in the Dump function. -* 20-APR-2006 (DSB): -* Guard against undefined transformations in Copy. -* 10-MAY-2006 (DSB): -* Override astEqual. -* 4-JUL-2008 (DSB): -* Fixed loop indexing problems in Equal function. -* 27-MAY-2011 (DSB): -* Added public method astPolyTran. -* 18-JUL-2011 (DSB): -* - Added attributes IterInverse, NiterInverse and TolInverse. -* - Do not report an error if astPolyTran fails to fit an inverse. -* 15-OCT-2011 (DSB): -* Improve argument checking and error reporting in PolyTran -* 8-MAY-2014 (DSB): -* Move to using CMinPack for minimisations. -* 11-NOV-2016 (DSB): -* - Fix bug in MapMerge that could cause a seg fault. It did not -* check that the PolyMap had a defined transformation before accessing -* the transformation's coefficient array. -* - Fix similar bugs in Equal that could cause seg faults. -* 15-MAR-2017 (DSB): -* - Change the GetTranForward and GetTranInverse functions so that they -* take into account the state of the Invert attribute. -* - Improve docs for the IterInverse attribute to explain that the -* inverse transformation replaced is always the original inverse -* transformation, as defined by the arguments supplied to the PolyMap -* constructor, regardless of the state of the Invert attribute. -* 17-MAR-2017 (DSB): -* - Add the astPolyCoeffs method. -* 30-MAR-2017 (DSB): -* Modify the astPolyTran method so that it can be used by the -* ChebyMap class to determine new transformations implemented as -* Chebyshev polynomials. -* 27-JUN-2017 (DSB): -* In SamplePoly1D/2D ensure the final sample on each axis does not -* go above the supplied upper bound. This can happen due to rounding -* error. This is important for ChebyMaps since points outside the -* bounds are set bad when transformed using a ChebyMap, causing NaNs -* to be generated in lmder1 (cminpack minimisation function). -* 3-JUL-2017 (DSB): -* Within FitPoly1D and FitPoly2D, use an initial guess that represents -* a unit mapping between normalised input and output values, rather -* than unnormalised PolyMap values. This is in case the PolyMap being -* fitted includes a change of scale (e.g. the PolyMap input is in "mm" -* but the output is in "rads" and includes some large scaling factor -* to do the conversion). -* 9-JAN-2018 (DSB): -* Correct Transform to take account of AST__BAD coeffs correctly. -* Previously bad coeffs could generate NaN output values. -*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 PolyMap - -/* 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 "cmpmap.h" /* Compound mappings */ -#include "polymap.h" /* Interface definition for this class */ -#include "unitmap.h" /* Unit mappings */ -#include "cminpack/cminpack.h" /* Levenberg - Marquardt minimization */ -#include "pal.h" /* SLALIB function definitions */ - -/* 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> -#include <limits.h> -#include <float.h> - -/* Module Variables. */ -/* ================= */ - -/* Address of this static variable is used as a unique identifier for - member of this class. */ -static int class_check; - -/* Pointers to parent class methods which are extended by this class. */ -static AstPointSet *(* parent_transform)( AstMapping *, AstPointSet *, int, AstPointSet *, int * ); -static const char *(* parent_getattrib)( AstObject *, const char *, int * ); -static int (* parent_testattrib)( AstObject *, const char *, int * ); -static void (* parent_clearattrib)( AstObject *, const char *, int * ); -static void (* parent_setattrib)( AstObject *, const char *, int * ); -static int (* parent_getobjsize)( AstObject *, int * ); - -#if defined(THREAD_SAFE) -static int (* parent_managelock)( AstObject *, int, int, AstObject **, int * ); -#endif - - -#ifdef THREAD_SAFE -/* Define how to initialise thread-specific globals. */ -#define GLOBAL_inits \ - globals->Class_Init = 0; \ - globals->GetAttrib_Buff[ 0 ] = 0; - -/* Create the function that initialises global data for this module. */ -astMAKE_INITGLOBALS(PolyMap) - -/* Define macros for accessing each item of thread specific global data. */ -#define class_init astGLOBAL(PolyMap,Class_Init) -#define class_vtab astGLOBAL(PolyMap,Class_Vtab) -#define getattrib_buff astGLOBAL(LutMap,GetAttrib_Buff) - -#include <pthread.h> - - -#else - -static char getattrib_buff[ 101 ]; - -/* Define the class virtual function table and its initialisation flag - as static variables. */ -static AstPolyMapVtab 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. */ -AstPolyMap *astPolyMapId_( int, int, int, const double[], int, const double[], const char *, ... ); - - -/* Prototypes for Private Member Functions. */ -/* ======================================== */ -static AstPointSet *Transform( AstMapping *, AstPointSet *, int, AstPointSet *, int * ); -static AstPolyMap **GetJacobian( AstPolyMap *, int * ); -static AstPolyMap *PolyTran( AstPolyMap *, int, double, double, int, const double *, const double *, int * ); -static double **SamplePoly1D( AstPolyMap *, int, double **, double, double, int, int *, double[2], int * ); -static double **SamplePoly2D( AstPolyMap *, int, double **, const double *, const double *, int, int *, double[4], int * ); -static double *FitPoly1D( AstPolyMap *, int, int, double, int, double **, double[2], int *, double *, int * ); -static double *FitPoly2D( AstPolyMap *, int, int, double, int, double **, double[4], int *, double *, int * ); -static int Equal( AstObject *, AstObject *, int * ); -static int GetObjSize( AstObject *, int * ); -static int GetTranForward( AstMapping *, int * ); -static int GetTranInverse( AstMapping *, int * ); -static int MPFunc1D( void *, int, int, const double *, double *, double *, int, int ); -static int MPFunc2D( void *, int, int, const double *, double *, double *, int, int ); -static int MapMerge( AstMapping *, int, int, int *, AstMapping ***, int **, int * ); -static int ReplaceTransformation( AstPolyMap *, int, double, double, int, const double *, const double *, int * ); -static void Copy( const AstObject *, AstObject *, int * ); -static void Delete( AstObject *obj, int * ); -static void Dump( AstObject *, AstChannel *, int * ); -static void FreeArrays( AstPolyMap *, int, int * ); -static void IterInverse( AstPolyMap *, AstPointSet *, AstPointSet *, int * ); -static void LMFunc1D( const double *, double *, int, int, void * ); -static void LMFunc2D( const double *, double *, int, int, void * ); -static void LMJacob1D( const double *, double *, int, int, void * ); -static void LMJacob2D( const double *, double *, int, int, void * ); -static void PolyPowers( AstPolyMap *, double **, int, const int *, double **, int, int, int * ); -static void StoreArrays( AstPolyMap *, int, int, const double *, int * ); -static void PolyCoeffs( AstPolyMap *, int, int, double *, int *, int * ); -static void FitPoly1DInit( AstPolyMap *, int, double **, AstMinPackData *, double *, int *); -static void FitPoly2DInit( AstPolyMap *, int, double **, AstMinPackData *, double *, int *); - -#if defined(THREAD_SAFE) -static int ManageLock( AstObject *, int, int, AstObject **, int * ); -#endif - -static const char *GetAttrib( AstObject *, const char *, int * ); -static int TestAttrib( AstObject *, const char *, int * ); -static void ClearAttrib( AstObject *, const char *, int * ); -static void SetAttrib( AstObject *, const char *, int * ); - -static int GetIterInverse( AstPolyMap *, int * ); -static int TestIterInverse( AstPolyMap *, int * ); -static void ClearIterInverse( AstPolyMap *, int * ); -static void SetIterInverse( AstPolyMap *, int, int * ); - -static int GetNiterInverse( AstPolyMap *, int * ); -static int TestNiterInverse( AstPolyMap *, int * ); -static void ClearNiterInverse( AstPolyMap *, int * ); -static void SetNiterInverse( AstPolyMap *, int, int * ); - -static double GetTolInverse( AstPolyMap *, int * ); -static int TestTolInverse( AstPolyMap *, int * ); -static void ClearTolInverse( AstPolyMap *, int * ); -static void SetTolInverse( AstPolyMap *, double, int * ); - - -/* Member functions. */ -/* ================= */ - -static void ClearAttrib( AstObject *this_object, const char *attrib, int *status ) { -/* -* Name: -* ClearAttrib - -* Purpose: -* Clear an attribute value for a PolyMap. - -* Type: -* Private function. - -* Synopsis: -* #include "polymap.h" -* void ClearAttrib( AstObject *this, const char *attrib, int *status ) - -* Class Membership: -* PolyMap member function (over-rides the astClearAttrib protected -* method inherited from the Mapping class). - -* Description: -* This function clears the value of a specified attribute for a -* PolyMap, so that the default value will subsequently be used. - -* Parameters: -* this -* Pointer to the PolyMap. -* attrib -* Pointer to a null-terminated string specifying the attribute -* name. This should be in lower case with no surrounding white -* space. -* status -* Pointer to the inherited status variable. -*/ - -/* Local Variables: */ - AstPolyMap *this; /* Pointer to the PolyMap structure */ - -/* Check the global error status. */ - if ( !astOK ) return; - -/* Obtain a pointer to the PolyMap structure. */ - this = (AstPolyMap *) this_object; - -/* Check the attribute name and clear the appropriate attribute. */ - -/* IterInverse. */ -/* ------------ */ - if ( !strcmp( attrib, "iterinverse" ) ) { - astClearIterInverse( this ); - -/* NiterInverse. */ -/* ------------- */ - } else if ( !strcmp( attrib, "niterinverse" ) ) { - astClearNiterInverse( this ); - -/* TolInverse. */ -/* ----------- */ - } else if ( !strcmp( attrib, "tolinverse" ) ) { - astClearTolInverse( this ); - -/* If the attribute is still not recognised, pass it on to the parent - method for further interpretation. */ - } else { - (*parent_clearattrib)( this_object, attrib, status ); - } -} - -static int Equal( AstObject *this_object, AstObject *that_object, int *status ) { -/* -* Name: -* Equal - -* Purpose: -* Test if two PolyMaps are equivalent. - -* Type: -* Private function. - -* Synopsis: -* #include "polymap.h" -* int Equal( AstObject *this, AstObject *that, int *status ) - -* Class Membership: -* PolyMap 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 PolyMaps are equivalent. - -* Parameters: -* this -* Pointer to the first Object (a PolyMap). -* that -* Pointer to the second Object. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* One if the PolyMaps 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: */ - AstPolyMap *that; - AstPolyMap *this; - int i, j, k; - int nin; - int nout; - int result; - int tmp; - -/* Initialise. */ - result = 0; - -/* Check the global error status. */ - if ( !astOK ) return result; - -/* Obtain pointers to the two PolyMap structures. */ - this = (AstPolyMap *) this_object; - that = (AstPolyMap *) that_object; - -/* Check the second object is a PolyMap. We know the first is a - PolyMap since we have arrived at this implementation of the virtual - function. */ - if( astIsAPolyMap( that ) ) { - -/* Get the number of inputs and outputs and check they are the same for both. */ - nin = astGetNin( this ); - nout = astGetNout( this ); - if( astGetNin( that ) == nin && astGetNout( that ) == nout ) { - -/* If the Invert flags for the two PolyMaps differ, it may still be possible - for them to be equivalent. First compare the PolyMaps if their Invert - flags are the same. In this case all the attributes of the two PolyMaps - must be identical. */ - if( astGetInvert( this ) == astGetInvert( that ) ) { - -/* Assume they are the same until we find a difference. */ - result = 1; - -/* The "_f" and "_i" suffixes in the PolyMap structure array names refer - to the forward and inverse transformations of the original uninverted - PolyMap. So we need to ensure the "nout" and "nin" values also refer - to the uninverted values. So if the PolyMaps are inverted, swap nout - and nin. */ - if( astGetInvert( this ) ) { - tmp = nin; - nin = nout; - nout = tmp; - } - -/* Check properties of the forward transformation. */ - if( this->ncoeff_f && that->ncoeff_f ) { - for( i = 0; i < nout && result; i++ ) { - if( this->ncoeff_f[ i ] != that->ncoeff_f[ i ] ){ - result = 0; - } - } - } else if( this->ncoeff_f || that->ncoeff_f ) { - result = 0; - } - - if( this->mxpow_f && that->mxpow_f ) { - for( i = 0; i < nout && result; i++ ) { - if( this->mxpow_f[ i ] != that->mxpow_f[ i ] ){ - result = 0; - } - } - } else if( this->mxpow_f || that->mxpow_f ) { - result = 0; - } - - if( this->coeff_f && that->coeff_f ) { - for( i = 0; i < nout && result; i++ ) { - for( j = 0; j < this->ncoeff_f[ i ] && result; j++ ) { - if( !astEQUAL( this->coeff_f[ i ][ j ], - that->coeff_f[ i ][ j ] ) ) { - result = 0; - } - } - } - } else if( this->coeff_f || that->coeff_f ) { - result = 0; - } - - if( this->power_f && that->power_f ) { - for( i = 0; i < nout && result; i++ ) { - for( j = 0; j < this->ncoeff_f[ i ] && result; j++ ) { - for( k = 0; k < nin && result; k++ ) { - if( this->power_f[ i ][ j ][ k ] != - that->power_f[ i ][ j ][ k ] ) { - result = 0; - } - } - } - } - } else if( this->power_f || that->power_f ) { - result = 0; - } - -/* Check properties of the inverse transformation. */ - if( this->ncoeff_i && that->ncoeff_i ) { - for( i = 0; i < nout && result; i++ ) { - if( this->ncoeff_i[ i ] != that->ncoeff_i[ i ] ){ - result = 0; - } - } - } else if( this->ncoeff_i || that->ncoeff_i ) { - result = 0; - } - - if( this->mxpow_i && that->mxpow_i ) { - for( i = 0; i < nout && result; i++ ) { - if( this->mxpow_i[ i ] != that->mxpow_i[ i ] ){ - result = 0; - } - } - } else if( this->mxpow_i || that->mxpow_i ) { - result = 0; - } - - if( this->coeff_f && that->coeff_f ) { - for( i = 0; i < nout && result; i++ ) { - for( j = 0; j < this->ncoeff_f[ i ] && result; j++ ) { - if( !astEQUAL( this->coeff_f[ i ][ j ], - that->coeff_f[ i ][ j ] ) ) { - result = 0; - } - } - } - } else if( this->coeff_f || that->coeff_f ) { - result = 0; - } - - if( this->power_f && that->power_f ) { - for( i = 0; i < nout && result; i++ ) { - for( j = 0; j < this->ncoeff_f[ i ] && result; j++ ) { - for( k = 0; k < nin && result; k++ ) { - if( this->power_f[ i ][ j ][ k ] != - that->power_f[ i ][ j ][ k ] ) { - result = 0; - } - } - } - } - } else if( this->power_f || that->power_f ) { - result = 0; - } - -/* If the Invert flags for the two PolyMaps differ, the attributes of the two - PolyMaps must be inversely related to each other. */ - } else { - -/* In the specific case of a PolyMap, Invert flags must be equal. */ - result = 0; - - } - } - } - -/* If an error occurred, clear the result value. */ - if ( !astOK ) result = 0; - -/* Return the result, */ - return result; -} - -static double *FitPoly1D( AstPolyMap *this, int forward, int nsamp, double acc, - int order, double **table, double scales[2], int *ncoeff, - double *racc, int *status ){ -/* -* Name: -* FitPoly1D - -* Purpose: -* Fit a (1-in,1-out) polynomial to a supplied set of data. - -* Type: -* Private function. - -* Synopsis: -* double *FitPoly1D( AstPolyMap *this, int forward, int nsamp, double acc, -* int order, double **table, double scales[2], int *ncoeff, -* double *racc, int *status ) - -* Description: -* This function fits a least squares 1D polynomial curve to the -* positions in a supplied table, and returns the coefficients of a -* PolyMap to describe the fit. For the purposes of this function, -* the input to the fit is refered to as x1 and the output as y1. So -* the returned coefficients describe a PolyMap with forward -* transformation: -* -* y1 = P1( x1 ) - -* Parameters: -* this -* Pointer to the PolyMap. -* forward -* Non-zero if the forward transformation of "this" is being -* replaced. Zero if the inverse transformation of "this" is being -* replaced. -* nsamp -* The number of (x1,y1) positions in the supplied table. -* acc -* The required accuracy, expressed as a geodesic distance within -* the polynomials output space (not the normalised tabular space). -* order -* The maximum power (plus one) of x1 within P1. So for instance, a -* value of 3 refers to a quadratic polynomial. -* table -* Pointer to an array of 2 pointers. Each of these pointers points -* to an array of "nsamp" doubles, being the normalised and sampled -* values for x1 and y1 in that order. -* scales -* Array holding the scaling factors used to produced the normalised -* values in the two columns of the table. Multiplying the normalised -* table values by the scale factor produces input or output axis -* values for the returned PolyMap -* ncoeff -* Pointer to an int in which to return the number of coefficients -* described by the returned array. -* racc -* Pointer to a double in which to return the achieved accuracy -* (which may be greater than "acc"). -* status -* Pointer to the inherited status variable. - -* Returned Value: -* A pointer to an array of doubles defining the polynomial in the -* form required by the PolyMap contructor. The number of coefficients -* is returned via "ncoeff". If the polynomial could not be found, -* then NULL is returned. The returned pointer should be freed using -* astFree when no longer needed. - -*/ - -/* Local Variables: */ - AstMinPackData data; - double *coeffs; - double *pc; - double *pr; - double *pxp1; - double *result; - double *work1; - double *work2; - double *work4; - double f1; - double f2; - double maxterm; - double term; - double tv; - int *work3; - int info; - int k; - int ncof; - int w1; - -/* Initialise returned value */ - result = NULL; - *ncoeff = 0; - *racc = 10*acc; - -/* Check inherited status */ - if( !astOK ) return result; - -/* Number of coefficients per poly. */ - ncof = order; - -/* Initialise the elements of the structure. */ - data.order = order; - data.nsamp = nsamp; - data.init_jac = 1; - data.xp1 = astMalloc( nsamp*order*sizeof( double ) ); - data.xp2 = NULL; - data.y[ 0 ] = table[ 1 ]; - data.y[ 1 ] = NULL; - -/* Work space to hold coefficients. */ - coeffs = astMalloc( ncof*sizeof( double ) ); - -/* Other work space. */ - work1 = astMalloc( nsamp*sizeof( double ) ); - work2 = astMalloc( ncof*nsamp*sizeof( double ) ); - work3 = astMalloc( ncof*sizeof( int ) ); - work4 = astMalloc( (5*ncof+nsamp)*sizeof( double ) ); - if( astOK ) { - -/* Find all the required powers of x1 and store them in the "xp1" - component of the data structure. The required initialisation is done - differently for different subclasses of PolyMap, so we need to wrap - it up in a virtual function. */ - astFitPoly1DInit( this, forward, table, &data, scales ); - -/* The initial guess at the coefficient values represents a unit - transformation in (normalised) tabulated (x,y) values. Using normalised - values means that we are, effectively, including a guess at the linear - scaling factor between input and output of the PolyMap (e.g. the PolyMap - may have inputs in mm and outputs in radians). */ - for( k = 0; k < ncof; k++ ) coeffs[ k ] = 0.0; - coeffs[ 1 ] = 1.0; - -/* Find the best coefficients */ - info = lmder1( MPFunc1D, &data, nsamp, ncof, coeffs, work1, work2, nsamp, - sqrt(DBL_EPSILON), work3, work4, (5*ncof+nsamp) ); - if( info == 0 ) astError( AST__MNPCK, "astPolyMap(PolyTran): Minpack error " - "detected (possible programming error).", status ); - -/* Return the achieved accuracy. The "work1" array holds the normalised Y - residuals at each tabulated point. */ - pr = work1; - tv = 0.0; - for( k = 0; k < nsamp; k++,pr++ ) tv += (*pr)*(*pr); - *racc = scales[ 1 ]*sqrt( tv/nsamp ); - -/* The best fitting polynomial coefficient found above relate to the - polynomial between the scaled positions stored in "table". These - scaled positions are related to PolyMap input/output axis values via - the scale factors supplied in "scales". Find the initial factor for the - current output. */ - f1 = scales[ 1 ]; - f2 = 1.0; - -/* Look at each coefficient. */ - pc = coeffs; - for( w1 = 0; w1 < order; w1++,pc++ ) { - -/* Get a pointer to the powers of X1 appropriate for the current coefficient, - at the first sample. */ - pxp1 = data.xp1 + w1; - -/* We find the contribution which this coefficient makes to the total - polynomial value. Find the maximum contribution made at any sample - points. */ - maxterm = 0.0; - for( k = 0; k < nsamp; k++ ) { - -/* Get the absolute value of the polynomial term that uses the current - coefficient. */ - term = fabs( ( *pc )*( *pxp1 ) ); - -/* Update the maximum term found at any sample. */ - if( term > maxterm ) maxterm = term; - -/* Increment the pointers to refer to the next sample. */ - pxp1 += order; - } - -/* If the maximum contribution made by this term is less than the - required accuracy, set the coefficient value to zero. */ - if( maxterm*f1 < acc ) { - *pc = 0.0; - -/* Scale the best fitting polynomial coefficient found above to take - account of the fact that the tabulated input and output positions in - "table" were are not actual PolyMap input and output axis values, but - are scaled by the factors stored in "scales". */ - } else { - *pc *= f1/f2; - } - - f2 *= scales[ 0 ]; - } - -/* Convert the array of coefficients into PolyMap form. */ - result = astMalloc( ncof*3*sizeof( double ) ); - - *ncoeff = 0; - pr = result; - pc = coeffs; - for( w1 = 0; w1 < order; w1++,pc++ ) { - if( *pc != 0.0 ) { - *(pr++) = *pc; - *(pr++) = 1; - *(pr++) = w1; - (*ncoeff)++; - } - } - -/* Truncate the returned array. */ - result = astRealloc( result, (*ncoeff)*3*sizeof( double ) ); - } - -/* Free resources. */ - coeffs = astFree( coeffs ); - data.xp1 = astFree( data.xp1 ); - work1 = astFree( work1 ); - work2 = astFree( work2 ); - work3 = astFree( work3 ); - work4 = astFree( work4 ); - -/* Return the coefficient array. */ - return result; - -} - -static void FitPoly1DInit( AstPolyMap *this, int forward, double **table, - AstMinPackData *data, double *scales, int *status ){ -/* -*+ -* Name: -* astFitPoly1DInit - -* Purpose: -* Perform initialisation needed for FitPoly1D - -* Type: -* Protected function. - -* Synopsis: -* #include "polymap.h" -* void astFitPoly1DInit( AstPolyMap *this, int forward, double **table, -* AstMinPackData *data, double *scales, -* int *status ) - -* Class Membership: -* PolyMap virtual function. - -* Description: -* This function performs initialisation needed for FitPoly1D. - -* Parameters: -* this -* Pointer to the PolyMap. -* forward -* Non-zero if the forward transformation of "this" is being -* replaced. Zero if the inverse transformation of "this" is being -* replaced. -* table -* Pointer to an array of 2 pointers. Each of these pointers points -* to an array of "nsamp" doubles, being the scaled and sampled values -* for x1 and y1 in that order. -* data -* Pointer to a structure holding information to pass the the -* service function invoked by the minimisation function. -* scales -* Array holding the scaling factors for the two columns of the table. -* Multiplying the table values by the scale factor produces PolyMap -* input or output axis values. -*- -*/ - -/* Local Variables; */ - double *px1; - double *pxp1; - double tv; - double x1; - int k; - int w1; - -/* Check the local error status. */ - if ( !astOK ) return; - -/* Get pointers to the supplied x1 values. */ - px1 = table[ 0 ]; - -/* Get pointers to the location for the next power of x1. */ - pxp1 = data->xp1; - -/* Loop round all samples. */ - for( k = 0; k < data->nsamp; k++ ) { - -/* Get the current x1 value. */ - x1 = *(px1++); - -/* Find all the required powers of x1 and store them in the "xp1" - component of the data structure. */ - tv = 1.0; - for( w1 = 0; w1 < data->order; w1++ ) { - *(pxp1++) = tv; - tv *= x1; - } - } -} - -static double *FitPoly2D( AstPolyMap *this, int forward, int nsamp, double acc, - int order, double **table, double scales[4], - int *ncoeff, double *racc, int *status ){ -/* -* Name: -* FitPoly2D - -* Purpose: -* Fit a (2-in,2-out) polynomial to a supplied set of data. - -* Type: -* Private function. - -* Synopsis: -* double *FitPoly2D( AstPolyMap *this, int forward, int nsamp, double acc, -* int order, double **table, double scales[4], -* int *ncoeff, double *racc, int *status ) - -* Description: -* This function fits a pair of least squares 2D polynomial surfaces -* to the positions in a supplied table, and returns the coefficients -* of a PolyMap to describe the fit. For the purposes of this function, -* the inputs to the fit is refered to as (x1,x2) and the output as -* (y1,y2). So the returned coefficients describe a PolyMap with forward -* transformations: -* -* y1 = P1( x1, x2 ) -* y2 = P2( x1, x2 ) -* -* P1 and P2 have the same maximum power on each input (specified by -* the "order" parameter). - -* Parameters: -* this -* Pointer to the PolyMap. -* forward -* Non-zero if the forward transformation of "this" is being -* replaced. Zero if the inverse transformation of "this" is being -* replaced. -* nsamp -* The number of (x1,x2,y1,y2) positions in the supplied table. -* acc -* The required accuracy, expressed as a geodesic distance within -* the polynomials output space (not the normalised tabular space). -* order -* The maximum power (plus one) of x1 or x2 within P1 and P2. So for -* instance, a value of 3 refers to a quadratic polynomial. Note, cross -* terms with total powers greater than or equal to "order" are not -* included in the fit. So the maximum number of terms in the fitted -* polynomial is order*(order+1)/2. -* table -* Pointer to an array of 4 pointers. Each of these pointers points -* to an array of "nsamp" doubles, being the normalised and sampled -* values for x1, x2, y1 or y2 in that order. -* scales -* Array holding the scaling factors used to produced the normalised -* values in the four columns of the table. Multiplying the normalised -* table values by the scale factor produces input or output axis -* values for the returned PolyMap -* ncoeff -* Pointer to an ant in which to return the number of coefficients -* described by the returned array. -* racc -* Pointer to a double in which to return the achieved accuracy -* (which may be greater than "acc"). -* status -* Pointer to the inherited status variable. - -* Returned Value: -* A pointer to an array of doubles defining the polynomial in the -* form required by the PolyMap contructor. The number of coefficients -* is returned via "ncoeff". If the polynomial could not be found with -* sufficient accuracy , then NULL is returned. The returned pointer -* should be freed using astFree when no longer needed. - -*/ - -/* Local Variables: */ - AstMinPackData data; - double *coeffs; - double *pc; - double *pr; - double *pxp1; - double *pxp2; - double *result; - double *work1; - double *work2; - double *work4; - double f1; - double f20; - double f2; - double f3; - double facc; - double maxterm; - double term; - double tv; - int *work3; - int info; - int iout; - int k; - int ncof; - int w12; - int w1; - int w2; - -/* Initialise returned value */ - result = NULL; - *ncoeff = 0; - *racc = 10*acc; - -/* Check inherited status */ - if( !astOK ) return result; - -/* Number of coefficients per poly. */ - ncof = order*( order + 1 )/2; - -/* Initialise the elements of the structure. */ - data.order = order; - data.nsamp = nsamp; - data.init_jac = 1; - data.xp1 = astMalloc( nsamp*order*sizeof( double ) ); - data.xp2 = astMalloc( nsamp*order*sizeof( double ) ); - data.y[ 0 ] = table[ 2 ]; - data.y[ 1 ] = table[ 3 ]; - -/* Work space to hold coefficients. */ - coeffs = astMalloc( 2*ncof*sizeof( double ) ); - -/* Other work space. */ - work1 = astMalloc( 2*nsamp*sizeof( double ) ); - work2 = astMalloc( 4*ncof*nsamp*sizeof( double ) ); - work3 = astMalloc( 2*ncof*sizeof( int ) ); - work4 = astMalloc( 2*(5*ncof+nsamp)*sizeof( double ) ); - if( astOK ) { - -/* Find all the required powers of x1 and x2 and store them in the "xp1" - and "xp2" components of the data structure. The required initialisation - is done differently for different subclasses of PolyMap, so we need to - wrap it up in a virtual function. */ - astFitPoly2DInit( this, forward, table, &data, scales ); - -/* The initial guess at the coefficient values represents a unit - transformation in (normalised) tabulated (x,y) values. Using normalised - values means that we are, effectively, including a guess at the linear - scaling factor between input and output of the PolyMap (e.g. the PolyMap - may have inputs in mm and outputs in radians). */ - for( k = 0; k < 2*ncof; k++ ) coeffs[ k ] = 0.0; - coeffs[ 1 ] = 1.0; - coeffs[ 5 ] = 1.0; - -/* Find the best coefficients */ - info = lmder1( MPFunc2D, &data, 2*nsamp, 2*ncof, coeffs, work1, work2, - 2*nsamp, sqrt(DBL_EPSILON), work3, work4, 2*(5*ncof+nsamp) ); - if( info == 0 ) astError( AST__MNPCK, "astPolyMap(PolyTran): Minpack error " - "detected (possible programming error).", status ); - -/* Return the achieved accuracy. */ - pr = work1; - tv = 0.0; - for( k = 0; k < 2*nsamp; k++,pr++ ) tv += (*pr)*(*pr); - facc = 1.0/(scales[2]*scales[2]) + 1.0/(scales[3]*scales[3]); - *racc = sqrt( tv/(2*nsamp*facc) ); - -/* Pointer to the first coefficient. */ - pc = coeffs; - -/* Look at coefficients for each output in turn. */ - for( iout = 0; iout < 2 && astOK; iout++ ) { - -/* The best fitting polynomial coefficient found above relate to the - polynomial between the scaled positions stored in "table". These - scaled positions are related to PolyMap input/output axis values via - the scale factors supplied in "scales". Find the initial factor for the - current output. */ - f1 = scales[ 2 + iout ]; - -/* Look at each coefficient for the current output. */ - f20 = 1.0; - for( w12 = 0; w12 < order; w12++ ) { - f3 = 1.0; - f2 = f20; - for( w2 = 0; w2 <= w12; w2++,pc++ ) { - w1 = w12 - w2; - -/* Get pointers to the powers of X1 and X2 appropriate for the current - coefficient, at the first sample. */ - pxp1 = data.xp1 + w1; - pxp2 = data.xp2 + w2; - -/* We find the contribution which this coefficient makes to the total - polynomial value. Find the maximum contribution made at any sample - points. */ - maxterm = 0.0; - for( k = 0; k < nsamp; k++ ) { - -/* Get the absolute value of the polynomial term that uses the current - coefficient. */ - term = fabs( ( *pc )*( *pxp1 )*( *pxp2 ) ); - -/* Update the maximum term found at any sample. */ - if( term > maxterm ) maxterm = term; - -/* Increment the pointers to refer to the next sample. */ - pxp1 += order; - pxp2 += order; - } - -/* If the maximum contribution made by this term is less than the - required accuracy, set the coefficient value to zero. */ - if( maxterm*f1 < acc ) { - *pc = 0.0; - -/* Scale the best fitting polynomial coefficient found above to take - account of the fact that the tabulated input and output positions in - "table" were are not actual PolyMap input and output axis values, but - are scaled by the factors stored in "scales". */ - } else { - *pc *= f1/( f2*f3 ); - } - - f2 /= scales[ 0 ]; - f3 *= scales[ 1 ]; - } - - f20 *= scales[ 0 ]; - } - } - -/* Convert the array of coefficients into PolyMap form. */ - result = astMalloc( 2*ncof*4*sizeof( double ) ); - - *ncoeff = 0; - pr = result; - pc = coeffs; - for( iout = 0; iout < 2 && astOK; iout++ ) { - for( w12 = 0; w12 < order; w12++ ) { - for( w2 = 0; w2 <= w12; w2++,pc++ ) { - w1 = w12 - w2; - if( *pc != 0.0 ) { - *(pr++) = *pc; - *(pr++) = iout + 1; - *(pr++) = w1; - *(pr++) = w2; - (*ncoeff)++; - } - } - } - } - -/* Truncate the returned array. */ - result = astRealloc( result, (*ncoeff)*4*sizeof( double ) ); - } - -/* Free resources. */ - coeffs = astFree( coeffs ); - data.xp1 = astFree( data.xp1 ); - data.xp2 = astFree( data.xp2 ); - work1 = astFree( work1 ); - work2 = astFree( work2 ); - work3 = astFree( work3 ); - work4 = astFree( work4 ); - -/* Return the coefficient array. */ - return result; - -} - -static void FitPoly2DInit( AstPolyMap *this, int forward, double **table, - AstMinPackData *data, double *scales, int *status ){ -/* -*+ -* Name: -* astFitPoly2DInit - -* Purpose: -* Perform initialisation needed for FitPoly2D - -* Type: -* Protected function. - -* Synopsis: -* #include "polymap.h" -* void astFitPoly2DInit( AstPolyMap *this, int forward, double **table, -* AstMinPackData *data, double *scales, -* int *status ) - -* Class Membership: -* PolyMap virtual function. - -* Description: -* This function performs initialisation needed for FitPoly2D. - -* Parameters: -* this -* Pointer to the PolyMap. -* forward -* Non-zero if the forward transformation of "this" is being -* replaced. Zero if the inverse transformation of "this" is being -* replaced. -* table -* Pointer to an array of 4 pointers. Each of these pointers points -* to an array of "nsamp" doubles, being the scaled and sampled values -* for x1, x2, y1 or y2 in that order. -* data -* Pointer to a structure holding information to pass the the -* service function invoked by the minimisation function. -* scales -* Array holding the scaling factors for the four columns of the table. -* Multiplying the table values by the scale factor produces PolyMap -* input or output axis values. -*- -*/ - -/* Local Variables; */ - double *px1; - double *px2; - double *pxp1; - double *pxp2; - double tv; - double x1; - double x2; - int k; - int w1; - int w2; - -/* Check the local error status. */ - if ( !astOK ) return; - -/* Get pointers to the supplied x1 and x2 values. */ - px1 = table[ 0 ]; - px2 = table[ 1 ]; - -/* Get pointers to the location for the next power of x1 and x2. */ - pxp1 = data->xp1; - pxp2 = data->xp2; - -/* Loop round all samples. */ - for( k = 0; k < data->nsamp; k++ ) { - -/* Get the current x1 and x2 values. */ - x1 = *(px1++); - x2 = *(px2++); - -/* Find all the required powers of x1 and store them in the "xp1" - component of the data structure. */ - tv = 1.0; - for( w1 = 0; w1 < data->order; w1++ ) { - *(pxp1++) = tv; - tv *= x1; - } - -/* Find all the required powers of x2 and store them in the "xp2" - comonent of the data structure. */ - tv = 1.0; - for( w2 = 0; w2 < data->order; w2++ ) { - *(pxp2++) = tv; - tv *= x2; - } - } -} - -static void FreeArrays( AstPolyMap *this, int forward, int *status ) { -/* -* Name: -* FreeArrays - -* Purpose: -* Free the dynamic arrays contained within a PolyMap - -* Type: -* Private function. - -* Synopsis: -* void FreeArrays( AstPolyMap *this, int forward, int *status ) - -* Description: -* This function frees all the dynamic arrays allocated as part of a -* PolyMap. - -* Parameters: -* this -* Pointer to the PolyMap. -* forward -* If non-zero, the arrays for the forward transformation are freed. -* Otherwise, the arrays for the inverse transformation are freed. -* 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: */ - int nin; /* Number of inputs */ - int nout; /* Number of outputs */ - int i; /* Loop count */ - int j; /* Loop count */ - -/* Get the number of inputs and outputs of the uninverted Mapping. */ - nin = ( (AstMapping *) this )->nin; - nout = ( (AstMapping *) this )->nout; - -/* Free the dynamic arrays for the forward transformation. */ - if( forward ) { - - if( this->coeff_f ) { - for( i = 0; i < nout; i++ ) { - this->coeff_f[ i ] = astFree( this->coeff_f[ i ] ); - } - this->coeff_f = astFree( this->coeff_f ); - } - - if( this->power_f ) { - for( i = 0; i < nout; i++ ) { - if( this->ncoeff_f && this->power_f[ i ] ) { - for( j = 0; j < this->ncoeff_f[ i ]; j++ ) { - this->power_f[ i ][ j ] = astFree( this->power_f[ i ][ j ] ); - } - } - this->power_f[ i ] = astFree( this->power_f[ i ] ); - } - this->power_f = astFree( this->power_f ); - } - - this->ncoeff_f = astFree( this->ncoeff_f ); - this->mxpow_f = astFree( this->mxpow_f ); - -/* Free the dynamic arrays for the inverse transformation. */ - } else { - - if( this->coeff_i ) { - for( i = 0; i < nin; i++ ) { - this->coeff_i[ i ] = astFree( this->coeff_i[ i ] ); - } - this->coeff_i = astFree( this->coeff_i ); - } - - if(this->power_i ) { - for( i = 0; i < nin; i++ ) { - if( this->ncoeff_i && this->power_i[ i ] ) { - for( j = 0; j < this->ncoeff_i[ i ]; j++ ) { - this->power_i[ i ][ j ] = astFree( this->power_i[ i ][ j ] ); - } - } - this->power_i[ i ] = astFree( this->power_i[ i ] ); - } - this->power_i = astFree( this->power_i ); - } - - this->ncoeff_i = astFree( this->ncoeff_i ); - this->mxpow_i = astFree( this->mxpow_i ); - } -} - -static const char *GetAttrib( AstObject *this_object, const char *attrib, int *status ) { -/* -* Name: -* GetAttrib - -* Purpose: -* Get the value of a specified attribute for a PolyMap. - -* Type: -* Private function. - -* Synopsis: -* #include "polymap.h" -* const char *GetAttrib( AstObject *this, const char *attrib, int *status ) - -* Class Membership: -* PolyMap member function (over-rides the protected astGetAttrib -* method inherited from the Mapping class). - -* Description: -* This function returns a pointer to the value of a specified -* attribute for a PolyMap, formatted as a character string. - -* Parameters: -* this -* Pointer to the PolyMap. -* attrib -* Pointer to a null-terminated string containing the name of -* the attribute whose value is required. This name should be in -* lower case, with all white space removed. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* - Pointer to a null-terminated string containing the attribute -* value. - -* Notes: -* - The returned string pointer may point at memory allocated -* within the PolyMap, or at static memory. The contents of the -* string may be over-written or the pointer may become invalid -* following a further invocation of the same function or any -* modification of the PolyMap. A copy of the string should -* therefore be made if necessary. -* - 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: */ - astDECLARE_GLOBALS /* Pointer to thread-specific global data */ - AstPolyMap *this; /* Pointer to the PolyMap structure */ - const char *result; /* Pointer value to return */ - double dval; /* Floating point attribute value */ - int ival; /* Integer attribute value */ - -/* Initialise. */ - result = NULL; - -/* Check the global error status. */ - if ( !astOK ) return result; - -/* Get a pointer to the thread specific global data structure. */ - astGET_GLOBALS(this_object); - -/* Obtain a pointer to the PolyMap structure. */ - this = (AstPolyMap *) this_object; - -/* Compare "attrib" with each recognised attribute name in turn, - obtaining the value of the required attribute. If necessary, write - the value into "getattrib_buff" as a null-terminated string in an appropriate - format. Set "result" to point at the result string. */ - -/* IterInverse. */ -/* ------------ */ - if ( !strcmp( attrib, "iterinverse" ) ) { - ival = astGetIterInverse( this ); - if ( astOK ) { - (void) sprintf( getattrib_buff, "%d", ival ); - result = getattrib_buff; - } - -/* NiterInverse. */ -/* ------------- */ - } else if ( !strcmp( attrib, "niterinverse" ) ) { - ival = astGetNiterInverse( this ); - if ( astOK ) { - (void) sprintf( getattrib_buff, "%d", ival ); - result = getattrib_buff; - } - -/* TolInverse. */ -/* ----------- */ - } else if ( !strcmp( attrib, "tolinverse" ) ) { - dval = astGetTolInverse( this ); - if ( astOK ) { - (void) sprintf( getattrib_buff, "%.*g", AST__DBL_DIG, dval ); - result = getattrib_buff; - } - -/* If the attribute name was not recognised, pass it on to the parent - method for further interpretation. */ - } else { - result = (*parent_getattrib)( this_object, attrib, status ); - } - -/* Return the result. */ - return result; - -} - -static AstPolyMap **GetJacobian( AstPolyMap *this, int *status ){ -/* -* Name: -* GetJacobian - -* Purpose: -* Get a description of a Jacobian matrix for the fwd transformation -* of a PolyMap. - -* Type: -* Private function. - -* Synopsis: -* AstPolyMap **GetJacobian( AstPolyMap *this, int *status ) - -* Description: -* This function returns a set of PolyMaps which define the Jacobian -* matrix of the forward transformation of the supplied PolyMap. -* -* The Jacobian matrix has "nout" rows and "nin" columns, where "nin" -* and "nout" are the number of inputs and outputs of the supplied PolyMap. -* Row "i", column "j" of the matrix holds the rate of change of the -* i'th PolyMap output with respect to the j'th PolyMap input. -* -* Since the values in the Jacobian matrix vary across the input space -* of the PolyMap, the matrix is returned in the form of a set of new -* PolyMaps which generate the elements of the Jacobian for any given -* position in the input space. The "nout" values in a single column of -* the Jacobian matrix are generated by the "nout" outputs from a single -* new PolyMap. The whole matrix is described by "nin" PolyMaps. -* -* The returned PolyMaps are cached in the supplied PolyMap object in -* order to speed up subsequent calls to this function. - -* Parameters: -* this -* The PolyMap for which the Jacbian is required. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* A pointer to an array of "nin" PolyMap pointers, where "nin" is the number -* of inputs for the sipplied PolyMap. The returned array should not be changed -* in any way, and the PolyMaps should not be freed (they will be freed when -* the supplied PolyMap is deleted). - -*/ - -/* Local Variables: */ - double *coeffs; - double *pc; - int icof; - int icol; - int iin; - int irow; - int ncof; - int ncof_row; - int ncof_total; - int nin; - int nout; - int power; - -/* Check inherited status */ - if( !astOK ) return NULL; - -/* Ensure there is a Jacobian stored in the PolyMap. */ - if( !this->jacobian ) { - -/* Get the number of inputs and outputs. */ - nin = astGetNin( this ); - nout = astGetNout( this ); - -/* Allocate memory to hold pointers to the PolyMaps used to describe the - Jacobian matrix. */ - this->jacobian = astCalloc( nin, sizeof(AstPolyMap *) ); - -/* Find the total number of coefficients used to describe the supplied - PolyMap (forward transformation) and allocate work space to hold the - coefficients for a single new PolyMap forward transformation. */ - ncof = 0; - for( irow = 0; irow < nout; irow++ ) { - ncof += this->ncoeff_f[ irow ]; - } - coeffs = astMalloc( ncof*( 2 + nin )*sizeof( double ) ); - -/* Check pointers can be used safely. */ - if( astOK ) { - -/* The Jacobian matrix has "nout" rows and "nin" columns. The "nout" values - in a single column of the Jacobian matrix corresponds to the "nout" outputs - from a single PolyMap. The whole matrix is described by "nin" PolyMaps. - Loop over each column of the Matrix, creating the corresponding PolyMap - for each. */ - for( icol = 0; icol < nin; icol++ ) { - -/* Initialise the total number of coefficients used to describe the - element of the PolyMap. */ - ncof_total = 0; - -/* Loop over each row of the Jacobian matrix (i.e. each PolyMap output). */ - pc = coeffs; - for( irow = 0; irow < nout; irow++ ) { - -/* Loop over each coefficient used in the polynomial that generates the - current PolyMap output. */ - ncof_row = this->ncoeff_f[ irow ]; - for( icof = 0; icof < ncof_row; icof++ ) { - -/* Get the power of input "icol" associated with the current coefficient. */ - power = (int)( this->power_f[ irow ][ icof ][ icol ] + 0.5 ); - -/* We can skip the coefficient if the power is zero. */ - if( power > 0 ) { - ncof_total++; - -/* Store the coefficient value, modified so that it describes a - polynomial that has been differentiated with respect to input "icol". */ - *(pc++) = this->coeff_f[ irow ][ icof ]*power; - -/* Store the output PolyMap to which the coeff relates. */ - *(pc++) = irow + 1; - -/* Store the powers of the inputs associated with the coeff. These are - the same as the original powers, except that the power of "icol" - (the input with respect to which the output has been differentiated) - is reduced by one. */ - for( iin = 0; iin < nin; iin++ ) { - if( iin != icol ) { - *(pc++) = this->power_f[ irow ][ icof ][ iin ]; - } else { - *(pc++) = this->power_f[ irow ][ icof ][ iin ] - 1; - } - } - } - } - } - -/* Create the PolyMap and store a pointer to it in the jacobian array in - the supplied PolyMap. */ - (this->jacobian)[ icol ] = astPolyMap( nin, nout, ncof_total, coeffs, - 0, NULL, " ", status ); - } - } - -/* Free resources */ - coeffs = astFree( coeffs ); - } - -/* Return the Jacobian. */ - return this->jacobian; -} - -static int GetObjSize( AstObject *this_object, int *status ) { -/* -* Name: -* GetObjSize - -* Purpose: -* Return the in-memory size of an Object. - -* Type: -* Private function. - -* Synopsis: -* #include "polymap.h" -* int GetObjSize( AstObject *this, int *status ) - -* Class Membership: -* PolyMap member function (over-rides the astGetObjSize protected -* method inherited from the parent class). - -* Description: -* This function returns the in-memory size of the supplied PolyMap, -* in bytes. - -* Parameters: -* this -* Pointer to the PolyMap. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* The Object size, in bytes. - -* 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: */ - AstPolyMap *this; - int ic; - int nc; - int result; - -/* Initialise. */ - result = 0; - -/* Check the global error status. */ - if ( !astOK ) return result; - -/* Obtain a pointers to the PolyMap structure. */ - this = (AstPolyMap *) this_object; - -/* Invoke the GetObjSize method inherited from the parent class, and then - add on any components of the class structure defined by this class - which are stored in dynamically allocated memory. */ - result = (*parent_getobjsize)( this_object, status ); - - if( this->jacobian ) { - nc = astGetNin( this ); - for( ic = 0; ic < nc; ic++ ) { - result += astGetObjSize( (this->jacobian)[ ic ] ); - } - result += sizeof( AstPolyMap * )*nc; - } - -/* If an error occurred, clear the result value. */ - if ( !astOK ) result = 0; - -/* Return the result, */ - return result; -} - -static int GetTranForward( AstMapping *this, int *status ) { -/* -* -* Name: -* GetTranForward - -* Purpose: -* Determine if a PolyMap defines a forward coordinate transformation. - -* Type: -* Private function. - -* Synopsis: -* #include "mapping.h" -* int GetTranForward( AstMapping *this, int *status ) - -* Class Membership: -* PolyMap member function (over-rides the astGetTranForward method -* inherited from the Mapping class). - -* Description: -* This function returns a value indicating if the PolyMap is able -* to perform a forward coordinate transformation. - -* Parameters: -* this -* Pointer to the PolyMap. -* 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: */ - AstPolyMap *map; /* Pointer to PolyMap to be queried */ - int result; /* The returned value */ - -/* Check the global error status. */ - if ( !astOK ) return 0; - -/* Obtain a pointer to the PolyMap. */ - map = (AstPolyMap *) this; - -/* First deal with cases where the PolyMap has not been inverted. */ - if( ! astGetInvert( this ) ) { - -/* The PolyMap has a defined forward transformation if one or more - coefficients values were supplied for the original forward - transformation. It is not possible to replace the original forward - transformation with an iterative algorithm. */ - result = map->ncoeff_f ? 1 : 0; - -/* Now deal with cases where the PolyMap has been inverted. */ - } else { - -/* The PolyMap has a defined forward transformation if one or more - coefficients values were supplied for the original inverse - transformation, or if the original inverse transformation is being - approximated using an iterative algorithm. */ - result = ( map->ncoeff_i || astGetIterInverse( map ) ) ? 1 : 0; - } - -/* Return the result. */ - return result; -} - -static int GetTranInverse( AstMapping *this, int *status ) { -/* -* -* Name: -* GetTranInverse - -* Purpose: -* Determine if a PolyMap defines an inverse coordinate transformation. - -* Type: -* Private function. - -* Synopsis: -* #include "matrixmap.h" -* int GetTranInverse( AstMapping *this, int *status ) - -* Class Membership: -* PolyMap member function (over-rides the astGetTranInverse method -* inherited from the Mapping class). - -* Description: -* This function returns a value indicating if the PolyMap is able -* to perform an inverse coordinate transformation. - -* Parameters: -* this -* Pointer to the PolyMap. -* 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: */ - AstPolyMap *map; /* Pointer to PolyMap to be queried */ - int result; /* The returned value */ - -/* Check the global error status. */ - if ( !astOK ) return 0; - -/* Obtain a pointer to the PolyMap. */ - map = (AstPolyMap *) this; - -/* First deal with cases where the PolyMap has not been inverted. */ - if( ! astGetInvert( this ) ) { - -/* The PolyMap has a defined inverse transformation if one or more - coefficients values were supplied for the original inverse - transformation, or if the original inverse transformation is being - approximated using an iterative algorithm. */ - result = ( map->ncoeff_i || astGetIterInverse( map ) ) ? 1 : 0; - -/* Now deal with cases where the PolyMap has been inverted. */ - } else { - -/* The PolyMap has a defined inverse transformation if one or more - coefficients values were supplied for the original forward - transformation. It is not possible to replace the original forward - transformation with an iterative algorithm. */ - result = map->ncoeff_f ? 1 : 0; - } - -/* Return the result. */ - return result; -} - -void astInitPolyMapVtab_( AstPolyMapVtab *vtab, const char *name, int *status ) { -/* -*+ -* Name: -* astInitPolyMapVtab - -* Purpose: -* Initialise a virtual function table for a PolyMap. - -* Type: -* Protected function. - -* Synopsis: -* #include "polymap.h" -* void astInitPolyMapVtab( AstPolyMapVtab *vtab, const char *name ) - -* Class Membership: -* PolyMap vtab initialiser. - -* Description: -* This function initialises the component of a virtual function -* table which is used by the PolyMap 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 astIsAPolyMap) 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->PolyPowers = PolyPowers; - vtab->FitPoly1DInit = FitPoly1DInit; - vtab->FitPoly2DInit = FitPoly2DInit; - vtab->PolyTran = PolyTran; - vtab->PolyCoeffs = PolyCoeffs; - - vtab->ClearIterInverse = ClearIterInverse; - vtab->GetIterInverse = GetIterInverse; - vtab->SetIterInverse = SetIterInverse; - vtab->TestIterInverse = TestIterInverse; - - vtab->ClearNiterInverse = ClearNiterInverse; - vtab->GetNiterInverse = GetNiterInverse; - vtab->SetNiterInverse = SetNiterInverse; - vtab->TestNiterInverse = TestNiterInverse; - - vtab->ClearTolInverse = ClearTolInverse; - vtab->GetTolInverse = GetTolInverse; - vtab->SetTolInverse = SetTolInverse; - vtab->TestTolInverse = TestTolInverse; - -/* 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_getobjsize = object->GetObjSize; - object->GetObjSize = GetObjSize; - -#if defined(THREAD_SAFE) - parent_managelock = object->ManageLock; - object->ManageLock = ManageLock; -#endif - - parent_clearattrib = object->ClearAttrib; - object->ClearAttrib = ClearAttrib; - parent_getattrib = object->GetAttrib; - object->GetAttrib = GetAttrib; - parent_setattrib = object->SetAttrib; - object->SetAttrib = SetAttrib; - parent_testattrib = object->TestAttrib; - object->TestAttrib = TestAttrib; - - parent_transform = mapping->Transform; - mapping->Transform = Transform; - mapping->GetTranForward = GetTranForward; - mapping->GetTranInverse = GetTranInverse; - -/* Store replacement pointers for methods which will be over-ridden by - new member functions implemented here. */ - object->Equal = Equal; - mapping->MapMerge = MapMerge; - -/* Declare the destructor and copy constructor. */ - astSetDelete( (AstObjectVtab *) vtab, Delete ); - astSetCopy( (AstObjectVtab *) vtab, Copy ); - -/* Declare the class dump function. */ - astSetDump( vtab, Dump, "PolyMap", "Polynomial 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 void IterInverse( AstPolyMap *this, AstPointSet *out, AstPointSet *result, - int *status ){ -/* -* Name: -* IterInverse - -* Purpose: -* Use an iterative method to evaluate the inverse transformation of a -* PolyMap at a set of output positions. - -* Type: -* Private function. - -* Synopsis: -* void IterInverse( AstPolyMap *this, AstPointSet *out, AstPointSet *result, -* int *status ) - -* Description: -* This function transforms a set of PolyMap output positions using -* the inverse transformation of the PolyMap, to generate the corresponding -* input positions. An iterative Newton-Raphson method is used which -* only required the forward transformation of the PolyMap to be deifned. - -* Parameters: -* this -* The PolyMap. -* out -* A PointSet holding the PolyMap output positions that are to be -* transformed using the inverse transformation. -* result -* A PointSet into which the generated PolyMap input positions are to be -* stored. -* status -* Pointer to the inherited status variable. - -*/ - -/* Local Variables: */ - AstPointSet *work; - AstPointSet **ps_jac; - AstPolyMap **jacob; - double *vec; - double *pb; - double **ptr_work; - double ***ptr_jac; - double *mat; - double **ptr_out; - double **ptr_in; - double *pa; - double det; - double maxerr; - double vlensq; - double xlensq; - double xx; - int *flags; - int *iw; - int fwd; - int icol; - int icoord; - int ipoint; - int irow; - int iter; - int maxiter; - int nconv; - int ncoord; - int npoint; - int sing; - -/* Check inherited status */ - if( !astOK ) return; - -/* Check the PolyMap has equal numbers of inputs and outputs. */ - ncoord = astGetNin( this ); - if( ncoord != astGetNout( this ) ) { - astError( AST__INTER, "astTransform(%s): Supplied %s has unequal numbers" - " of inputs and outputs and therefore an iterative inverse " - "cannot be used (internal AST Programming errpr).", status, - astGetClass(this), astGetClass(this) ); - } - -/* Get information about the Jacobian matrix for the forward polynomial - transformation. This matrix is a ncoord X ncoord matrix, in which - element (row=I,col=J) is the rate of change of output coord I with - respect to input coord J, of the supplied PolyMap's forward transformation. - The numerical values of the matrix vary depending on where it is - evaluated within the input space of the PolyMap. For this reason, the - "jacob" variable holds a vector of "ncoord" PolyMaps. The outputs of - each of these PolyMaps corresponds to a single column in the Jacobian - matrix. */ - jacob = GetJacobian( this, status ); - -/* Get the number of points to be transformed. */ - npoint = astGetNpoint( out ); - -/* Get another PointSet to hold intermediate results. */ - work = astPointSet( npoint, ncoord, " ", status ); - -/* See if the PolyMap has been inverted. */ - fwd = !astGetInvert( this ); - -/* Get pointers to the data arrays for all PointSets. Note, here "in" and - "out" refer to inputs and outputs of the PolyMap (i.e. the forward - transformation). These are respectively *outputs* and *inputs* of the - inverse transformation. */ - ptr_in = astGetPoints( result ); /* Returned input positions */ - ptr_out = astGetPoints( out ); /* Supplied output positions */ - ptr_work = astGetPoints( work ); /* Work space */ - -/* Allocate an array of PointSets to hold the elements of the Jacobian - matrix. */ - ptr_jac = astMalloc( sizeof( double ** )*ncoord ); - ps_jac = astCalloc( ncoord, sizeof( AstPointSet * ) ); - if( astOK ) { - for( icoord = 0; icoord < ncoord; icoord++ ) { - ps_jac[ icoord ] = astPointSet( npoint, ncoord, " ", status ); - ptr_jac[ icoord ] = astGetPoints( ps_jac[ icoord ] ); - } - } - -/* Allocate an array to hold flags indicating if each position has - converged. Initialise it to hold zero at every element. */ - flags = astCalloc( npoint, sizeof( int ) ); - -/* Allocate memory to hold the Jacobian matrix at a single point. */ - mat = astMalloc( sizeof( double )*ncoord*ncoord ); - -/* Allocate memory to hold the offset vector. */ - vec = astMalloc( sizeof( double )*ncoord ); - -/* Allocate memory to hold work space for palDmat. */ - iw = astMalloc( sizeof( int )*ncoord ); - -/* Check pointers can be used safely. */ - if( astOK ) { - -/* Store the initial guess at the required input positions. We assume initially - that the inverse transformation is a unit mapping, and so we just copy - the supplied outputs positions to the results PointSet holding the - corresponding input positions. */ - for( icoord = 0; icoord < ncoord; icoord++ ) { - memcpy( ptr_in[ icoord ], ptr_out[ icoord ], sizeof( double )*npoint ); - } - -/* Get the maximum number of iterations to perform. */ - maxiter = astGetNiterInverse( this ); - -/* Get the target relative error for the returned input axis values, and - square it. */ - maxerr = astGetTolInverse( this ); - maxerr *= maxerr; - -/* Initialise the number of positions which have reached the required - accuracy. */ - nconv = 0; - -/* Loop round doing iterations of a Newton-Raphson algorithm, until - all points have achieved the required relative error, or the - maximum number of iterations have been performed. */ - for( iter = 0; iter < maxiter && nconv < npoint && astOK; iter++ ) { - -/* Use the forward transformation of the supplied PolyMap to transform - the current guesses at the required input positions into the - corresponding output positions. Store the results in the "work" - PointSet. */ - (void) astTransform( this, result, fwd, work ); - -/* Modify the work PointSet so that it holds the offsets from the output - positions produced by the current input position guesses, and the - required output positions. */ - for( icoord = 0; icoord < ncoord; icoord++ ) { - pa = ptr_out[ icoord ]; - pb = ptr_work[ icoord ]; - for( ipoint = 0; ipoint< npoint; ipoint++,pa++,pb++ ) { - if( *pa != AST__BAD && *pb != AST__BAD ){ - *pb = *pa - *pb; - } else { - *pb = AST__BAD; - } - } - } - -/* Evaluate the elements of the Jacobian matrix at the current input - position guesses. */ - for( icoord = 0; icoord < ncoord; icoord++ ) { - (void) astTransform( jacob[ icoord ], result, 1, ps_jac[ icoord ] ); - } - -/* For each position, we now invert the matrix equation - - Dy = Jacobian.Dx - - to find a guess at the vector (dx) holding the offsets from the - current input positions guesses to their required values. Loop over all - points. */ - for( ipoint = 0; ipoint < npoint; ipoint++ ) { - -/* Do not change positions that have already converged. */ - if( !flags[ ipoint ] ) { - -/* Get the numerical values for the elements of the Jacobian matrix at - the current point. */ - pa = mat; - for( irow = 0; irow < ncoord; irow++ ) { - for( icol = 0; icol < ncoord; icol++ ) { - *(pa++) = ptr_jac[ icol ][ irow ][ ipoint ]; - } - -/* Store the offset from the current output position to the required - output position. */ - vec[ irow ] = ptr_work[ irow ][ ipoint ]; - } - -/* Find the corresponding offset from the current input position to the required - input position. */ - palDmat( ncoord, mat, vec, &det, &sing, iw ); - -/* If the matrix was singular, the input position cannot be evaluated so - store a bad value for it and indicate it has converged. */ - if( sing ) { - for( icoord = 0; icoord < ncoord; icoord++ ) { - ptr_in[ icoord ][ ipoint ] = AST__BAD; - } - flags[ ipoint ] = 1; - nconv++; - -/* Otherwise, update the input position guess. */ - } else { - vlensq = 0.0; - xlensq = 0.0; - pa = vec; - for( icoord = 0; icoord < ncoord; icoord++,pa++ ) { - xx = ptr_in[ icoord ][ ipoint ] + (*pa); - ptr_in[ icoord ][ ipoint ] = xx; - xlensq += xx*xx; - vlensq += (*pa)*(*pa); - } - -/* Check for convergence. */ - if( vlensq < maxerr*xlensq ) { - flags[ ipoint ] = 1; - nconv++; - } - } - } - } - } - } - -/* Free resources. */ - vec = astFree( vec ); - iw = astFree( iw ); - mat = astFree( mat ); - flags = astFree( flags ); - work = astAnnul( work ); - - if( ps_jac ) { - for( icoord = 0; icoord < ncoord; icoord++ ) { - ps_jac[ icoord ] = astAnnul( ps_jac[ icoord ] ); - } - ps_jac = astFree( ps_jac ); - } - - ptr_jac = astFree( ptr_jac ); - -} - -static void LMFunc1D( const double *p, double *hx, int m, int n, void *adata ){ -/* -* Name: -* LMFunc1D - -* Purpose: -* Evaluate a test 1D polynomal. - -* Type: -* Private function. - -* Synopsis: -* void LMFunc1D( const double *p, double *hx, int m, int n, void *adata ) - -* Description: -* This function finds the residuals implied by a supplied set of -* candidate polynomial coefficients. Each residual is a candidate -* polynomial evaluated at one of the sample points, minus the -* supplied target value for the polynomial at that test point. -* -* The minimisation process minimises the sum of the squared residuals. - -* Parameters: -* p -* An array of "m" coefficients for the candidate polynomial. The -* coefficients are ordered C0, C1, C2, etc. -* hx -* An array in which to return the "n" residuals. The residual at -* sample "k" is returned in element (k). -* m -* The length of the "p" array. This should be equal to order. -* n -* The length of the "hx" array. This should be equal to nsamp. -* adata -* Pointer to a structure holding the sample positions and values, -* and other information. - -*/ - -/* Local Variables: */ - AstMinPackData *data; - double *px1; - double *py; - const double *vp; - double *vr; - double res; - int k; - int w1; - -/* Get a pointer to the data structure. */ - data = (AstMinPackData *) adata; - -/* Initialise a pointer to the current returned residual value. */ - vr = hx; - -/* Initialise a pointer to the sampled Y values for the polynomial output. */ - py = data->y[ 0 ]; - -/* Initialise a pointer to the powers of the input X values at the curent - (i.e. first) sample. */ - px1 = data->xp1; - -/* Loop over the index of the sample to which this residual refers. */ - for( k = 0; k < data->nsamp; k++ ) { - -/* Initialise a pointer to the first coefficient (the constant term) for the - polynomial output coordinate. */ - vp = p; - -/* Initialise this residual to hold the sampled Y value. Increment the - pointer to the next sampled value for the current polynomial output. */ - res = -( *(py++) ); - -/* Loop over the coefficients. */ - for( w1 = 0; w1 < data->order; w1++ ) { - -/* Increment the residual by the value of the current term Cw1*(x1^w1). - Increment the pointer to the next coefficient (C). Also increment the - pointer to the next higher power of X1. */ - res += ( *(vp++) )*( *(px1++) ); - } - -/* Store the complete residual in the returned array, and increment the - pointer to the next residual. */ - *(vr++) = res; - } -} - -static void LMFunc2D( const double *p, double *hx, int m, int n, void *adata ){ -/* -* Name: -* LMFunc2D - -* Purpose: -* Evaluate a test 2D polynomal. - -* Type: -* Private function. - -* Synopsis: -* void LMFunc2D( const double *p, double *hx, int m, int n, void *adata ) - -* Description: -* This function finds the residuals implied by a supplied set of -* candidate polynomial coefficients. Each residual is a candidate -* polynomial (either P1 or P2) evaluated at one of the sample points -* (x1,x2), minus the supplied target value for the polynomial at that -* test point. -* -* The minimisation process minimises the sum of the squared residuals. - -* Parameters: -* p -* An array of "m" coefficients for the candidate polynomial. All the -* coefficients for polynomial P1 come first, followed by those for P2. -* Within each each polynomial, the coefficients are order C00, C10, -* C01, C20, C11, C02, C30, C21, C12, C03, etc. So the coefficient -* of (x1^j*x2^k) (=Cjk) for polynomial Pi is stored in element -* [k + (j + k)*(j + k + 1)/2 + i*order*(order+1)/2] of the "p" array. -* hx -* An array in which to return the "n" residuals. The residual at -* sample "k" for polynomial "i" is returned in element (k + nsamp*i). -* m -* The length of the "p" array. This should be equal to order*(order+1). -* n -* The length of the "hx" array. This should be equal to 2*nsamp. -* adata -* Pointer to a structure holding the sample positions and values, -* and other information. - -*/ - -/* Local Variables: */ - AstMinPackData *data; - const double *vp0; - const double *vp; - double *py; - double *px1; - double *px10; - double *px20; - double *vr; - double res; - double *px2; - int iout; - int k; - int w12; - int w2; - -/* Get a pointer to the data structure. */ - data = (AstMinPackData *) adata; - -/* Initialise a pointer to the current returned residual value. */ - vr = hx; - -/* Initialise a pointer to the first coefficient (the constant term) for the - current (i.e. first) polynomial output coordinate. */ - vp0 = p; - -/* Loop over each polynomial output coordinate. */ - for( iout = 0; iout < 2; iout++ ) { - -/* Initialise a pointer to the sampled Y values for the first polynomial - output. */ - py = data->y[ iout ]; - -/* Initialise pointers to the powers of the input X values at the curent - (i.e. first) sample. */ - px10 = data->xp1; - px20 = data->xp2; - -/* Loop over the index of the sample to which this residual refers. */ - for( k = 0; k < data->nsamp; k++ ) { - -/* Reset the pointer to the first coefficient (the constant term) - for the current polynomial output coordinate. */ - vp = vp0; - -/* Initialise this residual to hold the sampled Y value. Increment the - pointer to the next sampled value for the current polynomial output. */ - res = -( *(py++) ); - -/* The w12 value is the sum of the powers of X1 and X2. So w12=0 - corresponds to the constant term in the polynomial, and (e.g.) w12=6 - corresponds to all the terms for which the sum of the powerss (w1+w2) - is 6. Loop over all possible w12 values. */ - for( w12 = 0; w12 < data->order; w12++ ) { - -/* The next coeff refers to (x1^w12)*(x2^0). Get pointers to the values - holding x1^w12 and x2^0. */ - px1 = px10++; - px2 = px20; - -/* Loop over powers of x2. The corresponding power of x1 is "w12-w2", but - is not explicitly needed here. So x1 moves down from w12 to zero, as - w2 moves up from zero to w12. */ - for( w2 = 0; w2 <= w12; w2++ ) { - -/* Increment the residual by the value of the current term Cw1w2*(x1^w1)*(x2^w2). - Increment the pointer tio the next coefficient (C). Also decrement the - pointer to the next lower power of X1, and increment the pointer to the next - higher power of X2. */ - res += ( *(vp++) )*( *(px1--) )*( *(px2++) ); - } - } - -/* Move on to the x2 powers for the next sample. Don't need to do this - for x1 since px10 is incremented within the above loop. */ - px20 += data->order; - -/* Store the complete residual in the returned array, and increment the - pointer to the next residual. */ - *(vr++) = res; - } - -/* Get a pointer to the first coefficient (the constant term) for the - next polynomial output coordinate. */ - vp0 += data->order*( 1 + data->order )/2; - } -} - -static void LMJacob1D( const double *p, double *jac, int m, int n, void *adata ){ -/* -* Name: -* LMJacob1D - -* Purpose: -* Evaluate the Jacobian matrix of a test 1D polynomal. - -* Type: -* Private function. - -* Synopsis: -* void LMJacob1D( const double *p, double *jac, int m, int n, void *adata ) - -* Description: -* This function finds the Jacobian matrix that describes the rate of -* change of every residual with respect to every polynomial coefficient. -* Each residual is a candidate polynomial evaluated at one of the sample -* points minus the supplied target value for the polynomial at that test -* point. -* -* For a polynomial the Jacobian matrix is constant (i.e. does not -* depend on the values of the polynomial coefficients). So we only -* evaluate it on the first call. - -* Parameters: -* p -* An array of "m" coefficients for the candidate polynomial. -* jac -* An array in which to return the "m*n" elements of the Jacobian -* matrix. The rate of change of residual "r" with respect to -* coefficient "c" is returned in element "r + c*n". The residual -* at sample "k" of polynomial Pi has an "r" index of (k + nsamp*i). -* The coefficient of (x1^j) for polynomial Pi has a "c" index -* of j. -* m -* The length of the "p" array. This should be equal to order. -* n -* The number of residuals. This should be equal to nsamp. -* adata -* Pointer to a structure holding the sample positions and values, -* and other information. - -*/ - -/* Local Variables: */ - AstMinPackData *data; - int k; - int w1; - -/* Get a pointer to the data structure. */ - data = (AstMinPackData *) adata; - -/* The Jacobian of the residuals with respect to the polynomial - coefficients is constant (i.e. does not depend on the values of the - polynomial coefficients). So we only need to calculate it once. If - this is the first call, calculate the Jacobian and return it in "jac". - otherwise, just return immediately retaining the supplied "jac" values - (which will be the values returned by the previous call to this - function). */ - if( data->init_jac ) { - data->init_jac = 0; - -/* Loop over all residuals. */ - for( k = 0; k < n; k++ ) { - -/* Loop over all parameters (i.e. polynomial coefficients). */ - for( w1 = 0; w1 < m; w1++ ) { - -/* Store the Jacobian. */ - jac[ k + w1*n ] = (data->xp1)[ w1 + k*data->order ]; - } - } - } -} - -static void LMJacob2D( const double *p, double *jac, int m, int n, void *adata ){ -/* -* Name: -* LMJacob2D - -* Purpose: -* Evaluate the Jacobian matrix of a test 2D polynomal. - -* Type: -* Private function. - -* Synopsis: -* void LMJacob2D( const double *p, double *jac, int m, int n, void *adata ) - -* Description: -* This function finds the Jacobian matrix that describes the rate of -* change of every residual with respect to every polynomial coefficient. -* Each residual is a candidate polynomial (either P1 or P2) evaluated -* at one of the sample points (x1,x2), minus the supplied target value -* for the polynomial at that test point. -* -* For a polynomial the Jacobian matrix is constant (i.e. does not -* depend on the values of the polynomial coefficients). So we only -* evaluate it on the first call. - -* Parameters: -* p -* An array of "m" coefficients for the candidate polynomial. All the -* coefficients for polynomial P1 come first, followed by those for P2. -* Within each each polynomial, the coefficients are order C00, C10, -* C01, C20, C11, C02, C30, C21, C12, C03, etc. So the coefficient -* of (x1^j*x2^k) (=Cjk) for polynomial Pi is stored in element -* [k + (j + k)*(j + k + 1)/2 + i*order*(order+1)/2] of the "p" array. -* jac -* An array in which to return the "m*n" elements of the Jacobian -* matrix. The rate of change of residual "r" with respect to -* coefficient "c" is returned in element "r + c*n". The residual -* at sample "k" of polynomial Pi has an "r" index of (k + nsamp*i). -* The coefficient of (x1^j*x2^k) for polynomial Pi has a "c" index -* of [k + (j + k)*(j + k + 1)/2 + i*order*(order+1)/2]. -* m -* The length of the "p" array. This should be equal to order*(order+1). -* n -* The number of residuals. This should be equal to 2*nsamp. -* adata -* Pointer to a structure holdin gthe sample positions and values, -* and other information. - -*/ - -/* Local Variables: */ - AstMinPackData *data; - int iout; - int k; - int ncof; - int vp; - int vr; - int w1; - int w12; - int w2; - -/* Get a pointer to the data structure. */ - data = (AstMinPackData *) adata; - -/* The Jacobian of the residuals with respect to the polynomial - coefficients is constant (i.e. does not depend on the values of the - polynomial coefficients). So we only need to calculate it once. If - this is the first call, calculate the Jacobian and return it in "jac". - otherwise, just return immediately retaining the supplied "jac" values - (which will be the values returned by the previous call to this - function). */ - if( data->init_jac ) { - data->init_jac = 0; - -/* Store the number of coefficients in one polynomial. */ - ncof = data->order*( 1 + data->order )/2; - -/* Loop over all residuals. */ - for( vr = 0; vr < n; vr++ ) { - -/* Calculate the polynomial output index, and sample index, that creates - the current residual. */ - iout = vr/data->nsamp; - k = vr - iout*data->nsamp; - -/* Loop over all parameters (i.e. polynomial coefficients). */ - for( vp = 0; vp < m; vp++ ) { - -/* If this coefficient is not used in the creation of the current - polynomial output value, then the Jacobian value is zero. */ - if( vp/ncof != iout ) { - jac[ vr + vp*n ] = 0.0; - -/* Otherwise, get the powers of the two polynomial inputs, to which - the current coefficient relates. */ - } else { - w12 = ( -1.0 + sqrt( 1.0 + 8.0*(vp - iout*ncof) ) )/2.0; - w2 = vp - iout*ncof - w12*( w12 + 1 )/2; - w1 = w12 - w2; - -/* Store the Jacobian. */ - jac[ vr + vp*n ] = (data->xp1)[ w1 + k*data->order ]* - (data->xp2)[ w2 + k*data->order ]; - } - } - } - } -} - -#if defined(THREAD_SAFE) -static int ManageLock( AstObject *this_object, int mode, int extra, - AstObject **fail, int *status ) { -/* -* Name: -* ManageLock - -* Purpose: -* Manage the thread lock on an Object. - -* Type: -* Private function. - -* Synopsis: -* #include "object.h" -* AstObject *ManageLock( AstObject *this, int mode, int extra, -* AstObject **fail, int *status ) - -* Class Membership: -* PolyMap member function (over-rides the astManageLock protected -* method inherited from the parent class). - -* Description: -* This function manages the thread lock on the supplied Object. The -* lock can be locked, unlocked or checked by this function as -* deteremined by parameter "mode". See astLock for details of the way -* these locks are used. - -* Parameters: -* this -* Pointer to the Object. -* mode -* An integer flag indicating what the function should do: -* -* AST__LOCK: Lock the Object for exclusive use by the calling -* thread. The "extra" value indicates what should be done if the -* Object is already locked (wait or report an error - see astLock). -* -* AST__UNLOCK: Unlock the Object for use by other threads. -* -* AST__CHECKLOCK: Check that the object is locked for use by the -* calling thread (report an error if not). -* extra -* Extra mode-specific information. -* fail -* If a non-zero function value is returned, a pointer to the -* Object that caused the failure is returned at "*fail". This may -* be "this" or it may be an Object contained within "this". Note, -* the Object's reference count is not incremented, and so the -* returned pointer should not be annulled. A NULL pointer is -* returned if this function returns a value of zero. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* A local status value: -* 0 - Success -* 1 - Could not lock or unlock the object because it was already -* locked by another thread. -* 2 - Failed to lock a POSIX mutex -* 3 - Failed to unlock a POSIX mutex -* 4 - Bad "mode" value supplied. - -* Notes: -* - This function attempts to execute even if an error has already -* occurred. -*/ - -/* Local Variables: */ - AstPolyMap *this; - int ic; - int result; - int nc; - -/* Initialise */ - result = 0; - -/* Check the supplied pointer is not NULL. */ - if( !this_object ) return result; - -/* Obtain a pointer to the PolyMap structure. */ - this = (AstPolyMap *) this_object; - -/* Invoke the ManageLock method inherited from the parent class. */ - if( !result ) result = (*parent_managelock)( this_object, mode, extra, - fail, status ); - -/* Invoke the astManageLock method on any Objects contained within - the supplied Object. */ - if( this->jacobian ) { - nc = astGetNin( this ); - for( ic = 0; ic < nc && result; ic++ ) { - result = astManageLock( (this->jacobian)[ ic ], mode, - extra, fail ); - } - } - - return result; - -} -#endif - -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 PolyMap. - -* 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: -* PolyMap 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 PolyMap 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 PolyMap 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 PolyMap which is to be merged with -* its neighbours. This should be a cloned copy of the PolyMap -* pointer contained in the array element "(*map_list)[where]" -* (see below). This pointer will not be annulled, and the -* PolyMap it identifies will not be modified by this function. -* where -* Index in the "*map_list" array (below) at which the pointer -* to the nominated PolyMap 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: */ - AstPolyMap *pmap0; /* Nominated PolyMap */ - AstPolyMap *pmap1; /* Neighbouring PolyMap */ - int i; /* Index of neighbour */ - int nin; /* Number of input coordinates for nominated PolyMap */ - int nout; /* Number of output coordinates for nominated PolyMap */ - int ok; /* Are PolyMaps equivalent? */ - int oldinv0; /* Original Invert value in pmap0 */ - int oldinv1; /* Original Invert value in pmap1 */ - int result; /* Result value to return */ - -/* Initialise. */ - result = -1; - -/* Check the global error status. */ - if ( !astOK ) return result; - -/* Save a pointer to the nominated PolyMap. */ - pmap0 = (AstPolyMap *) ( *map_list )[ where ]; - -/* The only simplification which can currently be performed is to merge a PolyMap - with its own inverse. This can only be done in series. Obviously, - there are potentially other simplications which could be performed, but - time does not currently allow these to be coded. */ - if( series ) { - -/* Temporarily set the Invert flag of the nominated PolyMap to the - required value, first saving the original value so that it can be - re-instated later. */ - oldinv0 = astGetInvert( pmap0 ); - astSetInvert( pmap0, ( *invert_list )[ where ] ); - -/* Get the number of inputs and outputs to the nominated PolyMap. */ - nin = astGetNin( pmap0 ); - nout = astGetNout( pmap0 ); - -/* Check each neighbour. */ - for( i = where - 1; i <= where + 1; i += 2 ) { - -/* Continue with the next pass if the neighbour does not exist. */ - if( i < 0 || i >= *nmap ) continue; - -/* Continue with the next pass if this neighbour is not a PermMap. */ - if( ! astIsAPolyMap( ( *map_list )[ i ] ) ) continue; - -/* Get a pointer to it. */ - pmap1 = (AstPolyMap *) ( *map_list )[ i ]; - -/* Check it is used in the opposite direction to the nominated PolyMap. */ - if( ( *invert_list )[ i ] == ( *invert_list )[ where ] ) continue; - -/* We use the astEqual method to check that the two PolyMaps are equal. - But at the moment they may not be equal because they may have - different Invert flags. Therefore, temporarily set the invert flag - of the neighbour so that it is the same as the nominated PolyMap, - first saving the original value so that it can be re-instated later. - Note, we have already checked that the two PolyMaps are used in opposite - directions within the CmpMap. */ - oldinv1 = astGetInvert( pmap1 ); - astSetInvert( pmap1, ( *invert_list )[ where ] ); - -/* Use astEqual to check that the coefficients etc are equal in the two - PolyMaps. */ - ok = astEqual( pmap0, pmap1 ); - -/* Re-instate the original value of the Invert flag in the neighbour. */ - astSetInvert( pmap1, oldinv1 ); - -/* Pass on to the next neighbour if the current neighbour differs from - the nominated PolyMap. */ - if( !ok ) continue; - -/* If we get this far, then the nominated PolyMap and the current - neighbour cancel each other out, so replace each by a UnitMap. */ - pmap0 = astAnnul( pmap0 ); - pmap1 = astAnnul( pmap1 ); - if( i < where ) { - ( *map_list )[ where ] = (AstMapping *) astUnitMap( nout, "", status ); - ( *map_list )[ i ] = (AstMapping *) astUnitMap( nout, "", status ); - ( *invert_list )[ where ] = 0; - ( *invert_list )[ i ] = 0; - result = i; - } else { - ( *map_list )[ where ] = (AstMapping *) astUnitMap( nin, "", status ); - ( *map_list )[ i ] = (AstMapping *) astUnitMap( nin, "", status ); - ( *invert_list )[ where ] = 0; - ( *invert_list )[ i ] = 0; - result = where; - } - -/* Leave the loop. */ - break; - } - -/* If the nominated PolyMap was not replaced by a UnitMap, then - re-instate its original value for the Invert flag. */ - if( pmap0 ) astSetInvert( pmap0, oldinv0 ); - } - -/* Return the result. */ - return result; -} - -static int MPFunc1D( void *p, int m, int n, const double *x, double *fvec, - double *fjac, int ldfjac, int iflag ) { -/* -* Name: -* MPFunc1D - -* Purpose: -* Evaluate a test 1D polynomal or its Jacobian. - -* Type: -* Private function. - -* Synopsis: -* int MPFunc1D( void *p, int m, int n, const double *x, double *fvec, -* double *fjac, int ldfjac, int iflag ) - -* Description: -* This function finds the residuals or Jacobian implied by a supplied -* set of candidate polynomial coefficients. Each residual is a candidate -* polynomial evaluated at one of the sample points, minus the -* supplied target value for the polynomial at that test point. -* -* The minimisation process minimises the sum of the squared residuals. - -* Parameters: -* p -* A pointer to data needed to evaluate the required residuals or -* Jacobians. -* m -* The number of residuals. -* n -* The number of coefficients. -* x -* An array of "n" coefficients for the candidate polynomial. -* fvec -* An array in which to return the "m" residuals. -* fjac -* An array in which to return the "mxn" Jacobian values. -* ldflac -* Unused (should always be equal to "m"). -* iflag -* If 1 calculate the residuals, otherwise calculate the Jacobians. - -*/ - -/* If required, calculate the function values at X, and return this - vector in fvec. Do not alter fjac. */ - if( iflag == 1 ) { - LMFunc1D( x, fvec, n, m, p ); - -/* Otherwise, calculate the Jacobian values at X, and return this - matrix in fjac. Do not alter fvec. */ - } else { - LMJacob1D( x, fjac, n, m, p ); - } - - return 0; -} - -static int MPFunc2D( void *p, int m, int n, const double *x, double *fvec, - double *fjac, int ldfjac, int iflag ) { -/* -* Name: -* MPFunc1D - -* Purpose: -* Evaluate a test 2D polynomal or its Jacobian. - -* Type: -* Private function. - -* Synopsis: -* int MPFunc2D( void *p, int m, int n, const double *x, double *fvec, -* double *fjac, int ldfjac, int iflag ) - -* Description: -* This function finds the residuals or Jacobian implied by a supplied -* set of candidate polynomial coefficients. Each residual is a candidate -* polynomial evaluated at one of the sample points, minus the -* supplied target value for the polynomial at that test point. -* -* The minimisation process minimises the sum of the squared residuals. - -* Parameters: -* p -* A pointer to data needed to evaluate the required residuals or -* Jacobians. -* m -* The number of residuals. -* n -* The number of coefficients. -* x -* An array of "n" coefficients for the candidate polynomial. -* fvec -* An array in which to return the "m" residuals. -* fjac -* An array in which to return the "mxn" Jacobian values. -* ldflac -* Unused (should always be equal to "m"). -* iflag -* If 1 calculate the residuals, otherwise calculate the Jacobians. - -*/ - -/* If required, calculate the function values at X, and return this - vector in fvec. Do not alter fjac. */ - if( iflag == 1 ) { - LMFunc2D( x, fvec, n, m, p ); - -/* Otherwise, calculate the Jacobian values at X, and return this - matrix in fjac. Do not alter fvec. */ - } else { - LMJacob2D( x, fjac, n, m, p ); - - } - - return 0; -} - -static void PolyCoeffs( AstPolyMap *this, int forward, int nel, double *coeffs, - int *ncoeff, int *status ){ -/* -*++ -* Name: -c astPolyCoeffs -f AST_POLYCOEFFS - -* Purpose: -* Retrieve the coefficient values used by a PolyMap. - -* Type: -* Public function. - -* Synopsis: -c #include "polymap.h" -c void astPolyCoeffs( AstPolyMap *this, int forward, int nel, double *coeffs, -c int *ncoeff ) -f CALL AST_POLYCOEFFS( THIS, FORWARD, NEL, COEFFS, NCOEFF, STATUS ) - -* Class Membership: -* PolyMap method. - -* Description: -* This function returns the coefficient values used by either the -* forward or inverse transformation of a PolyMap, in the same form -* that they are supplied to the PolyMap constructor. -* -* Usually, you should call this method first with -c "nel" -f NEL -* set to zero to determine the number of coefficients used by the -* PolyMap. This allows you to allocate an array of the correct size to -* hold all coefficient data. You should then call this method a -* second time to get the coefficient data. - -* Parameters: -c this -f THIS = INTEGER (Given) -* Pointer to the original Mapping. -c forward -f FORWARD = LOGICAL (Given) -c If non-zero, -f If .TRUE., -* the coefficients of the forward PolyMap transformation are -* returned. Otherwise the inverse transformation coefficients are -* returned. -c nel -f NEL = INTEGER (Given) -* The length of the supplied -c "coeffs" -f COEFFS -* array. It should be at least "ncoeff*( nin + 2 )" if -c "foward" is non-zero, -f FORWARD is .TRUE., -* and "ncoeff*( nout + 2 )" otherwise, where "ncoeff" is the -* number of coefficients to be returned. If a value of zero -* is supplied, no coefficient values are returned, but the -* number of coefficients used by the transformation is still -* returned in -c "ncoeff". -f NCOEFF. -c coeffs -f COEFFS( NEL ) = DOUBLE PRECISION (Returned) -* An array in which to return the coefficients used by the -* requested transformation of the PolyMap. Ignored if -c "nel" is zero. -f NEL is zero. -* The coefficient data is returned in the form in which it is -* supplied to the PolyMap constructor. That is, each group of -* "2 + nin" or "2 + nout" adjacent elements describe a single -* coefficient of the forward or inverse transformation. See the -* PolyMap constructor documentation for further details. -* -* If the supplied array is too short to hold all the coefficients, -* trailing coefficients are excluded. If the supplied array is -* longer than needed to hold all the coefficients, trailing -* elements are filled with zeros. -c ncoeff -f NCOEFF = INTEGER (Returned) -* The number of coefficients used by the requested transformation. -* A value of zero is returned if the transformation does not -* have any defining polynomials. A value is returned for this argument -* even if -c "nel" is zero. -f NEL is zero. -f STATUS = INTEGER (Given and Returned) -f The global status. - -*-- -*/ - -/* Local Variables: */ - double **coeff; - int ***power; - int *nco; - int *pp; - int iax; - int icoeff; - int iel; - int ipoly; - int nax; - int npoly; - -/* Initialise */ - *ncoeff = 0; - -/* Check the inherited status. */ - if ( !astOK ) return; - -/* Fill any supplied array with zeros. */ - if( nel ) memset( coeffs, 0, nel*sizeof( *coeffs ) ); - -/* Get the values to use, taking account of whether the PolyMap has been - inverted or not. */ - if( forward != astGetInvert( this ) ){ - nco = this->ncoeff_f; - power = this->power_f; - coeff = this->coeff_f; - npoly = astGetNout( this ); - nax = astGetNin( this ); - } else { - nco = this->ncoeff_i; - power = this->power_i; - coeff = this->coeff_i; - npoly = astGetNin( this ); - nax = astGetNout( this ); - } - -/* Notheg to do if there are no coeffs. */ - if( nco && power && coeff ) { - -/* Initialise index of next value to store in returned array. */ - iel = 0; - -/* Loop round each 1D polynomial. */ - for( ipoly = 0; ipoly < npoly; ipoly++ ){ - -/* Loop round coefficients. */ - for( icoeff = 0; icoeff < nco[ ipoly ]; icoeff++ ) { - -/* Store the coefficient value in the next element of the returned array, - if the array is not already full. Increment the pointer to the next - element. */ - if( iel < nel ) coeffs[ iel++ ] = coeff[ipoly][icoeff]; - -/* Next store the integer index of the PolyMap input or output which uses - the coefficient within its defining polynomial (the first axis has - index 1). */ - if( iel < nel ) coeffs[ iel++ ] = ipoly + 1; - -/* The remaining elements of the group give the integer powers to use - with each input or output coordinate value. */ - pp = power[ipoly][icoeff]; - for( iax = 0; iax < nax; iax++,pp++ ) { - if( iel < nel ) coeffs[ iel++ ] = *pp; - } - } - -/* Increment the total number of coefficients used by the transformation. */ - *ncoeff += nco[ ipoly ]; - } - } -} - -static void PolyPowers( AstPolyMap *this, double **work, int ncoord, - const int *mxpow, double **ptr, int point, - int fwd, int *status ){ -/* -*+ -* Name: -* astPolyPowers - -* Purpose: -* Find the required powers of the input axis values. - -* Type: -* Protected function. - -* Synopsis: -* #include "polymap.h" -* void astPolyPowers( AstPolyMap *this, double **work, int ncoord, -* const int *mxpow, double **ptr, int point, -* int fwd ) - -* Class Membership: -* PolyMap virtual function. - -* Description: -* This function is used by astTransform to calculate the powers of -* the axis values for a single input position. In the case of -* sub-classes, the powers may not be simply powers of the supplied -* axis values but may be more complex quantities such as a Chebyshev -* polynomial of the required degree evaluated at the input axis values. - -* Parameters: -* this -* Pointer to the PolyMap. -* work -* An array of "ncoord" pointers, each pointing to an array of -* length "max(2,mxpow)". The required values are placed in this -* array on exit. -* ncoord -* The number of axes. -* mxpow -* Pointer to an array holding the maximum power required of each -* axis value. Should have "ncoord" elements. -* ptr -* An array of "ncoord" pointers, each pointing to an array holding -* the axis values. Each of these arrays of axis values must have -* at least "point+1" elements. -* point -* The zero based index of the point within "ptr" that holds the -* axis values to be exponentiated. -* fwd -* Do the supplied coefficients define the foward transformation of -* the PolyMap? -*- -*/ - -/* Local Variables; */ - double *pwork; - double x; - int coord; - int ip; - -/* Check the local error status. */ - if ( !astOK ) return; - -/* For the base PolyMap class, this method simply raises each input axis - value to the required power. Loop over all input axes. */ - for( coord = 0; coord < ncoord; coord++ ) { - -/* Get a pointer to the array in which the powers of the current axis - value are to be returned. */ - pwork = work[ coord ]; - -/* Anything to the power zero is 1.0. */ - pwork[ 0 ] = 1.0; - -/* Get the input axis value. If it is bad, store bad values for all - remaining powers. */ - x = ptr[ coord ][ point ]; - if( x == AST__BAD ) { - for( ip = 1; ip <= mxpow[ coord ]; ip++ ) pwork[ ip ] = AST__BAD; - -/* Otherwise, form and store the required powers of the input axis value. */ - } else { - for( ip = 1; ip <= mxpow[ coord ]; ip++ ) { - pwork[ ip ] = pwork[ ip - 1 ]*x; - } - } - } -} - -static AstPolyMap *PolyTran( AstPolyMap *this, int forward, double acc, - double maxacc, int maxorder, const double *lbnd, - const double *ubnd, int *status ){ -/* -*++ -* Name: -c astPolyTran -f AST_POLYTRAN - -* Purpose: -* Fit a PolyMap inverse or forward transformation. - -* Type: -* Public function. - -* Synopsis: -c #include "polymap.h" -c AstPolyMap *astPolyTran( AstPolyMap *this, int forward, double acc, -c double maxacc, int maxorder, const double *lbnd, -c const double *ubnd ) -f RESULT = AST_POLYTRAN( THIS, FORWARD, ACC, MAXACC, MAXORDER, LBND, -f UBND, STATUS ) - -* Class Membership: -* PolyMap method. - -* Description: -* This function creates a new PolyMap which is a copy of the supplied -* PolyMap, in which a specified transformation (forward or inverse) -* has been replaced by a new polynomial transformation. The -* coefficients of the new transformation are estimated by sampling -* the other transformation and performing a least squares polynomial -* fit in the opposite direction to the sampled positions and values. -* -* This method can only be used on (1-input,1-output) or (2-input,2-output) -* PolyMaps. -* -* The transformation to create is specified by the -c "forward" parameter. -f FORWARD parameter. -* In what follows "X" refers to the inputs of the PolyMap, and "Y" to -* the outputs of the PolyMap. The forward transformation transforms -* input values (X) into output values (Y), and the inverse transformation -* transforms output values (Y) into input values (X). Within a PolyMap, -* each transformation is represented by an independent set of -* polynomials, P_f or P_i: Y=P_f(X) for the forward transformation and -* X=P_i(Y) for the inverse transformation. -* -c The "forward" -f The FORWARD -* parameter specifies the transformation to be replaced. If it is -c non-zero, -f is .TRUE., -* a new forward transformation is created -* by first finding the input values (X) using the inverse transformation -* (which must be available) at a regular grid of points (Y) covering a -* rectangular region of the PolyMap's output space. The coefficients of -* the required forward polynomial, Y=P_f(X), are chosen in order to -* minimise the sum of the squared residuals between the sampled values -* of Y and P_f(X). -* -c If "forward" is zero (probably the most likely case), -f If FORWARD is .FALSE. (probably the most likely case), -* a new inverse transformation is created by -* first finding the output values (Y) using the forward transformation -* (which must be available) at a regular grid of points (X) covering a -* rectangular region of the PolyMap's input space. The coefficients of -* the required inverse polynomial, X=P_i(Y), are chosen in order to -* minimise the sum of the squared residuals between the sampled values -* of X and P_i(Y). -* -* This fitting process is performed repeatedly with increasing -* polynomial orders (starting with linear) until the target -* accuracy is achieved, or a specified maximum order is reached. If -* the target accuracy cannot be achieved even with this maximum-order -* polynomial, the best fitting maximum-order polynomial is returned so -* long as its accuracy is better than -c "maxacc". -f MAXACC. -* If it is not, a NULL pointer is returned but no error is reported. - -* Parameters: -c this -f THIS = INTEGER (Given) -* Pointer to the original Mapping. -c forward -f FORWARD = LOGICAL (Given) -c If non-zero, -f If .TRUE., -* the forward PolyMap transformation is replaced. Otherwise the -* inverse transformation is replaced. -c acc -f ACC = DOUBLE (Given) -* The target accuracy, expressed as a geodesic distance within -* the PolyMap's input space (if -c "forward" is zero) or output space (if "forward" is non-zero). -f FORWARD is .FALSE.) or output space (if FORWARD is .TRUE.). -c maxacc -f MAXACC = DOUBLE (Given) -* The maximum allowed accuracy for an acceptable polynomial, -* expressed as a geodesic distance within the PolyMap's input -* space (if -c "forward" is zero) or output space (if "forward" is non-zero). -f FORWARD is .FALSE.) or output space (if FORWARD is .TRUE.). -c maxorder -f MAXORDER = INTEGER (Given) -* The maximum allowed polynomial order. This is one more than the -* maximum power of either input axis. So for instance, a value of -* 3 refers to a quadratic polynomial. Note, cross terms with total -* powers greater than or equal to -c maxorder -f MAXORDER -* are not inlcuded in the fit. So the maximum number of terms in -* each of the fitted polynomials is -c maxorder*(maxorder+1)/2. -f MAXORDER*(MAXORDER+1)/2. -c lbnd -f LBND( * ) = DOUBLE PRECISION (Given) -c Pointer to an -f An -* array holding the lower bounds of a rectangular region within -* the PolyMap's input space (if -c "forward" is zero) or output space (if "forward" is non-zero). -f FORWARD is .FALSE.) or output space (if FORWARD is .TRUE.). -* The new polynomial will be evaluated over this rectangle. The -* length of this array should equal the value of the PolyMap's Nin -* or Nout attribute, depending on -c "forward". -f FORWARD. -c ubnd -f UBND( * ) = DOUBLE PRECISION (Given) -c Pointer to an -f An -* array holding the upper bounds of a rectangular region within -* the PolyMap's input space (if -c "forward" is zero) or output space (if "forward" is non-zero). -f FORWARD is .FALSE.) or output space (if FORWARD is .TRUE.). -* The new polynomial will be evaluated over this rectangle. The -* length of this array should equal the value of the PolyMap's Nin -* or Nout attribute, depending on -c "forward". -f FORWARD. -f STATUS = INTEGER (Given and Returned) -f The global status. - -* Returned Value: -c astPolyTran() -f AST_POLYTRAN = INTEGER -* A pointer to the new PolyMap. -c A NULL pointer -f AST__NULL -* will be returned if the fit fails to achieve the accuracy -* specified by -c "maxacc", -f MAXACC, -* but no error will be reported. - -* Applicability: -* PolyMap -* All PolyMaps have this method. -* ChebyMap -c The ChebyMap implementation of this method allows -c NULL pointers to be supplied for "lbnd" and/or "ubnd", -c in which case the corresponding bounds supplied when the ChebyMap -c was created are used. -* The returned PolyMap will be a ChebyMap, and the new transformation -* will be defined as a weighted sum of Chebyshev functions of the -* first kind. - -* Notes: -* - This function can only be used on 1D or 2D PolyMaps which have -* the same number of inputs and outputs. -* - 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. -*-- -*/ - -/* Local Variables: */ - AstPolyMap *result; - int ok; - -/* Initialise. */ - result = NULL; - -/* Check the inherited status. */ - if ( !astOK ) return result; - -/* Take a copy of the supplied PolyMap. */ - result = astCopy( this ); - -/* Replace the required transformation. */ - ok = ReplaceTransformation( result, forward, acc, maxacc, maxorder, lbnd, - ubnd, status ); - -/* If an error occurred, or the fit was not good enough, annul the returned - PolyMap. */ - if ( !ok || !astOK ) result = astAnnul( result ); - -/* Return the result. */ - return result; -} - -static int ReplaceTransformation( AstPolyMap *this, int forward, double acc, - double maxacc, int maxorder, const double *lbnd, - const double *ubnd, int *status ){ -/* -* Name: -* ReplaceTransformation - -* Purpose: -* Create a new inverse or forward transformation for a PolyMap. - -* Type: -* Private function. - -* Synopsis: -* int ReplaceTransformation( AstPolyMap *this, int forward, double acc, -* double maxacc, int maxorder, const double *lbnd, -* const double *ubnd, int *status ) - -* Description: -* This function creates a new forward or inverse transformation for -* the supplied PolyMap (replacing any existing transformation), by -* sampling the other transformation and performing a least squares -* polynomial fit to the sample positions and values. -* -* The transformation to create is specified by the "forward" parameter. -* In what follows "X" refers to the inputs of the PolyMap, and "Y" to -* the outputs of the PolyMap. The forward transformation transforms -* input values (X) into output values (Y), and the inverse transformation -* transforms output values (Y) into input values (X). Within a PolyMap, -* each transformation is represented by an independent set of -* polynomials: Y=P_f(X) for the forward transformation and X=P_i(Y) -* for the inverse transformation. -* -* If "forward" is zero then a new inverse transformation is created by -* first finding the output values (Y) using the forward transformation -* (which must be available) at a regular grid of points (X) covering a -* rectangular region of the PolyMap's input space. The coefficients of -* the required inverse polynomial, X=P_i(Y), are chosen in order to -* minimise the sum of the squared residuals between the sampled values -* of X and P_i(Y). -* -* If "forward" is non-zero then a new forward transformation is created -* by first finding the input values (X) using the inverse transformation -* (which must be available) at a regular grid of points (Y) covering a -* rectangular region of the PolyMap's output space. The coefficients of -* the required forward polynomial, Y=P_f(X), are chosen in order to -* minimise the sum of the squared residuals between the sampled values -* of Y and P_f(X). -* -* This fitting process is performed repeatedly with increasing -* polynomial orders (starting with linear) until the target -* accuracy is achieved, or a specified maximum order is reached. If -* the target accuracy cannot be achieved even with this maximum-order -* polynomial, the best fitting maximum-order polynomial is returned so -* long as its accuracy is better than "maxacc". - -* Parameters: -* this -* The PolyMap. -* forward -* If non-zero, then the forward PolyMap transformation is -* replaced. Otherwise the inverse transformation is replaced. -* acc -* The target accuracy, expressed as a geodesic distance within -* the PolyMap's input space (if "forward" is zero) or output -* space (if "forward" is non-zero). -* maxacc -* The maximum allowed accuracy for an acceptable polynomial, -* expressed as a geodesic distance within the PolyMap's input -* space (if "forward" is zero) or output space (if "forward" is -* non-zero). -* maxorder -* The maximum allowed polynomial order. This is one more than the -* maximum power of either input axis. So for instance, a value of -* 3 refers to a quadratic polynomial. Note, cross terms with total -* powers greater than or equal to maxorder are not inlcuded in the -* fit. So the maximum number of terms in each of the fitted polynomials -* is maxorder*(maxorder+1)/2. -* lbnd -* An array holding the lower bounds of a rectangular region within -* the PolyMap's input space (if "forward" is zero) or output -* space (if "forward" is non-zero). The new polynomial will be -* evaluated over this rectangle. -* ubnd -* An array holding the upper bounds of a rectangular region within -* the PolyMap's input space (if "forward" is zero) or output -* space (if "forward" is non-zero). The new polynomial will be -* evaluated over this rectangle. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* - Non-zero if a fit was performed succesfully, to at least the -* maximum allowed accuracy. Zero if the fit failed, or an error -* occurred. - -* Notes: -* - No error is reported if the fit fails to achieve the required -* maximum accuracy. -* - An error is reported if the transformation that is not being -* replaced is not defined. -* - An error is reported if the PolyMap does not have equal numbers -* of inputs and outputs. -* - An error is reported if the PolyMap has more than 2 inputs or outputs. - -*/ - -/* Local Variables: */ - double **table; - double *cofs; - double racc; - double scales[ 4 ]; - int ndim; - int ncof; - int nsamp; - int order; - int result; - -/* Check inherited status */ - if( !astOK ) return 0; - -/* Check the PolyMap can be used. */ - ndim = astGetNin( this ); - if( astGetNout( this ) != ndim && astOK ) { - astError( AST__BADNI, "astPolyTran(%s): Supplied %s has " - "different number of inputs (%d) and outputs (%d).", - status, astGetClass( this ), astGetClass( this ), ndim, - astGetNout( this ) ); - - } else if( ndim > 2 && astOK ) { - astError( AST__BADNI, "astPolyTran(%s): Supplied %s has " - "too many inputs and outputs (%d) - must be 1 or 2.", - status, astGetClass( this ), astGetClass( this ), ndim ); - } - - if( forward != astGetInvert( this ) ){ - if( ! this->ncoeff_i && astOK ) { - astError( AST__NODEF, "astPolyTran(%s): Supplied %s has " - "no inverse transformation.", status, astGetClass( this ), - astGetClass( this ) ); - } - } else { - if( ! this->ncoeff_f && astOK ) { - astError( AST__NODEF, "astPolyTran(%s): Supplied %s has " - "no forward transformation.", status, astGetClass( this ), - astGetClass( this ) ); - } - } - -/* Check the bounds can be used. */ - if( !lbnd || !ubnd ) { - astError( AST__NODEF, "astPolyTran(%s): No upper and/or lower bounds " - "supplied.", status, astGetClass( this ) ); - } else if( lbnd[ 0 ] >= ubnd[ 0 ] && astOK ) { - astError( AST__NODEF, "astPolyTran(%s): Supplied upper " - "bound for the first axis (%g) is less than or equal to the " - "supplied lower bound (%g).", status, astGetClass( this ), - lbnd[ 0 ], ubnd[ 0 ] ); - } else if( ndim == 2 && lbnd[ 1 ] >= ubnd[ 1 ] && astOK ) { - astError( AST__NODEF, "astPolyTran(%s): Supplied upper " - "bound for the second axis (%g) is less than or equal to the " - "supplied lower bound (%g).", status, astGetClass( this ), - lbnd[ 1 ], ubnd[ 1 ] ); - } - -/* Initialise pointer to work space. */ - table = NULL; - -/* Loop over increasing polynomial orders until the required accuracy is - achieved, up to the specified maximum order. The "order" value is one more - than the maximum power in the polynomial (so a quadratic has "order" 3). */ - if( maxorder < 2 ) maxorder = 2; - for( order = 2; order <= maxorder; order++ ) { - -/* First do 2D PolyMaps. */ - if( ndim == 2 ) { - -/* Sample the requested polynomial transformation at a grid of points. This - grid covers the user-supplied region, using 2*order points on each - axis. */ - table = SamplePoly2D( this, !forward, table, lbnd, ubnd, 2*order, - &nsamp, scales, status ); - -/* Fit the polynomial. Always fit a linear polynomial ("order" 2) to any - dummy second axis. If successfull, replace the PolyMap transformation - and break out of the order loop. */ - cofs = FitPoly2D( this, forward, nsamp, acc, order, table, scales, - &ncof, &racc, status ); - -/* Now do 1D PolyMaps. */ - } else { - table = SamplePoly1D( this, !forward, table, lbnd[ 0 ], ubnd[ 0 ], - 2*order, &nsamp, scales, status ); - cofs = FitPoly1D( this, forward, nsamp, acc, order, table, scales, - &ncof, &racc, status ); - } - -/* If the fit was succesful, replace the PolyMap transformation and break - out of the order loop. */ - if( cofs && ( racc < acc || ( racc < maxacc && order == maxorder ) ) ) { - StoreArrays( this, forward, ncof, cofs, status ); - break; - } else { - cofs = astFree( cofs ); - } - } - -/* If no fit was produced, return zero. */ - result = cofs ? 1 : 0; - -/* Free resources. */ - cofs = astFree( cofs ); - table = astFreeDouble( table ); - -/* Return the result. */ - return result; -} - -static double **SamplePoly1D( AstPolyMap *this, int forward, double **table, - double lbnd, double ubnd, int npoint, int *nsamp, - double scales[2], int *status ){ -/* -* Name: -* SamplePoly1D - -* Purpose: -* Create a table of input and output positions for a 1D PolyMap. - -* Type: -* Private function. - -* Synopsis: -* double **SamplePoly1D( AstPolyMap *this, int forward, double **table, -* double lbnd, double ubnd, int npoint, int *nsamp, -* double scales[2], int *status ) - -* Description: -* This function creates a table containing samples of the requested -* polynomial transformation at a grid of input points. This grid covers -* the user-supplied region, using "npoint" points. - -* Parameters: -* this -* The PolyMap. -* forward -* If non-zero, then the forward PolyMap transformation is sampled. -* Otherwise the inverse transformation is sampled. -* table -* Pointer to a previous table created by this function, which is -* to be re-used, or NULL. -* lbnd -* The lower bounds of the region within the PolyMap's 1D input space -* (if "forward" is non-zero) or output space (if "forward" is zero). -* The new polynomial will be evaluated over this region. -* ubnd -* The upper bounds of the region within the PolyMap's 1D input space -* (if "forward" is non-zero) or output space (if "forward" is zero). -* The new polynomial will be evaluated over this region. -* npoint -* The number of points to use. -* nsamp -* Address of an int in which to return the total number of samples -* in the returned table. -* scales -* Array in which to return the scaling factors for the two columns -* of the returned table. Multiplying the returned table values by -* the scale factor produces PolyMap input or output axis values. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* Pointer to an array of 2 pointers. Each of these pointers points -* to an array of "nsamp" doubles, being the sampled values for y1 -* and x1 in that order. Here x1 is the input value for the sampled -* transformation (these are spaced on the regular grid specified -* by lbnd, ubnd and npoint), and y1 is the output position produced -* by the sampled transformation. The returned values are scaled so -* that each column has an RMS value of 1.0. The scaling factors that -* convert scaled values into original values are returned in "scales". -* The returned pointer should be freed using astFreeDouble when no -* longer needed. - -*/ - -/* Local Variables: */ - AstPointSet *ps1; - AstPointSet *ps2; - double **result; - double *p0; - double *p1; - double *ptr1[ 1 ]; - double *ptr2[ 1 ]; - double delta0; - double rms; - double sum; - double val0; - int i; - int icol; - -/* Initialise returned value */ - result = table; - *nsamp = 0; - -/* Check inherited status */ - if( !astOK ) return result; - -/* Ensure we have a table of the correct size. */ - *nsamp = npoint; - if( !result ) result = astCalloc( 2, sizeof( double * ) ); - if( result ) { - for( i = 0; i < 2; i++ ) { - result[ i ] = astRealloc( result[ i ] , (*nsamp)*sizeof( double ) ); - } - } - -/* Work out the step sizes for the grid. */ - delta0 = ( ubnd - lbnd )/( npoint - 1 ); - -/* Create a PointSet to hold the grid of input positions. Use column 1 - of the table to hold the PointSet values. */ - ps1 = astPointSet( *nsamp, 1, " ", status ); - ptr1[ 0 ] = result[ 1 ]; - astSetPoints( ps1, ptr1 ); - -/* Create a PointSet to hold the grid of output positions. Use column 0 - of the table to hold the PointSet values. */ - ps2 = astPointSet( *nsamp, 1, " ", status ); - ptr2[ 0 ] = result[ 0 ]; - astSetPoints( ps2, ptr2 ); - if( astOK ) { - -/* Calculate the grid of input positions and store in the PointSet and - therefore also in the returned table. Rounding error may cause the - final point to be just over the upper bound, so we leave the final - point out of the loop and set it explicitly to the upper bound - afterwards. */ - val0 = lbnd; - p0 = ptr1[ 0 ]; - for( i = 0; i < npoint-1; i++ ) { - *(p0++) = val0; - val0 += delta0; - } - *p0 = ubnd; - -/* Transform the input grid to get the output grid. */ - (void) astTransform( this, ps1, forward, ps2 ); - -/* Scale each column in turn. */ - for( icol = 0; icol < 2; icol++ ) { - -/* Find the RMS of the values in the column. */ - sum = 0.0; - p0 = result[ icol ]; - p1 = p0 + (*nsamp); - for( ; p0 < p1; p0++ ) sum += ( *p0 )*( *p0 ); - rms = sqrt( sum/(*nsamp) ); - -/* Divide the table values by the RMS. */ - p0 = result[ icol ]; - p1 = p0 + (*nsamp); - for( ; p0 < p1; p0++ ) *p0 /= rms; - -/* Return the RMS as the scale factor. */ - scales[ icol ] = rms; - } - } - -/* Free resources */ - ps1 = astAnnul( ps1 ); - ps2 = astAnnul( ps2 ); - -/* If an error occurred, free the returned array. */ - if( !astOK ) result = astFreeDouble( result ); - -/* Return a pointer to the table. */ - return result; -} - -static double **SamplePoly2D( AstPolyMap *this, int forward, double **table, - const double *lbnd, const double *ubnd, int npoint, - int *nsamp, double scales[4], int *status ){ -/* -* Name: -* SamplePoly2D - -* Purpose: -* Create a table of input and output positions for a 2D PolyMap. - -* Type: -* Private function. - -* Synopsis: -* double **SamplePoly2D( AstPolyMap *this, int forward, double **table, -* const double *lbnd, const double *ubnd, int npoint, -* int *nsamp, double scales[4], int *status ) - -* Description: -* This function creates a table containing samples of the requested -* polynomial transformation at a grid of input points. This grid covers -* the user-supplied region, using "npoint" points on each axis. - -* Parameters: -* this -* The PolyMap. -* forward -* If non-zero, then the forward PolyMap transformation is sampled. -* Otherwise the inverse transformation is sampled. -* table -* Pointer to a previous table created by this function, which is -* to be re-used, or NULL. -* lbnd -* An array holding the lower bounds of a rectangular region within -* the PolyMap's input space (if "forward" is non-zero) or output -* space (if "forward" is zero). The new polynomial will be -* evaluated over this rectangle. -* ubnd -* An array holding the upper bounds of a rectangular region within -* the PolyMap's input space (if "forward" is non-zero) or output -* space (if "forward" is zero). The new polynomial will be -* evaluated over this rectangle. -* npoint -* The number of points along each edge of the grid. -* nsamp -* Address of an int in which to return the total number of samples -* in the returned table. -* scales -* Array in which to return the scaling factors for the four -* columns of the returned table. Multiplying the returned table -* values by the scale factor produces PolyMap input or output axis -* values. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* Pointer to an array of 4 pointers. Each of these pointers points -* to an array of "nsamp" doubles, being the sampled values for y1, -* y2, x1 or x2 in that order. Here (x1,x2) are the input values -* for the sampled transformation (these are spaced on the regular -* grid specified by lbnd, ubnd and npoint), and (y1,y2) are the -* output positions produced by the sampled transformation. The -* returned values are scaled so that each column has an RMS value -* of 1.0. The scaling factors that convert scaled values into -* original values are returned in "scales". The returned pointer -* should be freed using astFreeDouble when no longer needed. - -*/ - -/* Local Variables: */ - AstPointSet *ps1; - AstPointSet *ps2; - double **result; - double *p0; - double *p1; - double *ptr1[ 2 ]; - double *ptr2[ 2 ]; - double delta1; - double delta0; - double rms; - double sum; - double val0; - double val1; - int i; - int icol; - int j; - -/* Initialise returned value */ - result = table; - *nsamp = 0; - -/* Check inherited status */ - if( !astOK ) return result; - -/* Ensure we have a table of the correct size. */ - *nsamp = npoint*npoint; - if( !result ) result = astCalloc( 4, sizeof( double * ) ); - if( result ) { - for( i = 0; i < 4; i++ ) { - result[ i ] = astRealloc( result[ i ] , (*nsamp)*sizeof( double ) ); - } - } - -/* Work out the step sizes for the grid. */ - delta0 = ( ubnd[ 0 ] - lbnd[ 0 ] )/( npoint - 1 ); - delta1 = ( ubnd[ 1 ] - lbnd[ 1 ] )/( npoint - 1 ); - -/* Create a PointSet to hold the grid of input positions. Use columns 2 - and 3 of the table to hold the PointSet values. */ - ps1 = astPointSet( *nsamp, 2, " ", status ); - ptr1[ 0 ] = result[ 2 ]; - ptr1[ 1 ] = result[ 3 ]; - astSetPoints( ps1, ptr1 ); - -/* Create a PointSet to hold the grid of output positions. Use columns 0 - and 1 of the table to hold the PointSet values. */ - ps2 = astPointSet( *nsamp, 2, " ", status ); - ptr2[ 0 ] = result[ 0 ]; - ptr2[ 1 ] = result[ 1 ]; - astSetPoints( ps2, ptr2 ); - if( astOK ) { - -/* Calculate the grid of input positions and store in the PointSet and - therefore also in the returned table. Rounding error may cause the - final point on each axis to be just over the upper bound, so we leave - the final point out of the loop and set it explicitly to the upper bound - afterwards. */ - val0 = lbnd[ 0 ]; - p0 = ptr1[ 0 ]; - p1 = ptr1[ 1 ]; - for( i = 0; i < npoint - 1; i++ ) { - val1 = lbnd[ 1 ]; - for( j = 0; j < npoint-1; j++ ) { - *(p0++) = val0; - *(p1++) = val1; - val1 += delta1; - } - *(p0++) = val0; - *(p1++) = ubnd[ 1 ]; - val0 += delta0; - } - - - val1 = lbnd[ 1 ]; - for( j = 0; j < npoint-1; j++ ) { - *(p0++) = ubnd[ 0 ]; - *(p1++) = val1; - val1 += delta1; - } - *(p0++) = ubnd[ 0 ]; - *(p1++) = ubnd[ 1 ]; - - -/* Transform the input grid to get the output grid. */ - (void) astTransform( this, ps1, forward, ps2 ); - -/* Scale each pair of columns in turn. Use the same scale factor for - each axis in order to ensure an isotropic metric. */ - for( icol = 0; icol < 4; icol += 2 ) { - -/* Find the RMS of the values in the two columns. */ - sum = 0.0; - p0 = result[ icol ]; - p1 = p0 + (*nsamp); - for( ; p0 < p1; p0++ ) sum += ( *p0 )*( *p0 ); - - p0 = result[ icol + 1 ]; - p1 = p0 + (*nsamp); - for( ; p0 < p1; p0++ ) sum += ( *p0 )*( *p0 ); - - rms = sqrt( sum/(2*(*nsamp)) ); - -/* Divide the table values by the RMS. */ - p0 = result[ icol ]; - p1 = p0 + (*nsamp); - for( ; p0 < p1; p0++ ) *p0 /= rms; - - p0 = result[ icol + 1 ]; - p1 = p0 + (*nsamp); - for( ; p0 < p1; p0++ ) *p0 /= rms; - -/* Return the RMS as the scale factor. */ - scales[ icol ] = rms; - scales[ icol + 1 ] = rms; - } - } - -/* Free resources */ - ps1 = astAnnul( ps1 ); - ps2 = astAnnul( ps2 ); - -/* If an error occurred, free the returned array. */ - if( !astOK ) result = astFreeDouble( result ); - -/* Return a pointer to the table. */ - return result; -} - -static void SetAttrib( AstObject *this_object, const char *setting, int *status ) { -/* -* Name: -* SetAttrib - -* Purpose: -* Set an attribute value for a PolyMap. - -* Type: -* Private function. - -* Synopsis: -* #include "polymap.h" -* void SetAttrib( AstObject *this, const char *setting ) - -* Class Membership: -* PolyMap member function (over-rides the astSetAttrib protected -* method inherited from the Mapping class). - -* Description: -* This function assigns an attribute value for a PolyMap, the -* attribute and its value being specified by means of a string of -* the form: -* -* "attribute= value " -* -* Here, "attribute" specifies the attribute name and should be in -* lower case with no white space present. The value to the right -* of the "=" should be a suitable textual representation of the -* value to be assigned and this will be interpreted according to -* the attribute's data type. White space surrounding the value is -* only significant for string attributes. - -* Parameters: -* this -* Pointer to the PolyMap. -* setting -* Pointer to a null-terminated string specifying the new attribute -* value. -*/ - -/* Local Variables: */ - AstPolyMap *this; /* Pointer to the PolyMap structure */ - double dval; /* Floating point attribute value */ - int ival; /* Integer attribute value */ - int len; /* Length of setting string */ - int nc; /* Number of characters read by astSscanf */ - -/* Check the global error status. */ - if ( !astOK ) return; - -/* Obtain a pointer to the PolyMap structure. */ - this = (AstPolyMap *) this_object; - -/* Obtain the length of the setting string. */ - len = (int) strlen( setting ); - -/* Test for each recognised attribute in turn, using "astSscanf" to parse - the setting string and extract the attribute value (or an offset to - it in the case of string values). In each case, use the value set - in "nc" to check that the entire string was matched. Once a value - has been obtained, use the appropriate method to set it. */ - -/* IterInverse. */ -/* ------------ */ - if ( nc = 0, - ( 1 == astSscanf( setting, "iterinverse= %d %n", &ival, &nc ) ) - && ( nc >= len ) ) { - astSetIterInverse( this, ival ); - -/* NiterInverse. */ -/* ------------- */ - } else if ( nc = 0, - ( 1 == astSscanf( setting, "niterinverse= %d %n", &ival, &nc ) ) - && ( nc >= len ) ) { - astSetNiterInverse( this, ival ); - -/* TolInverse. */ -/* ----------- */ - } else if ( nc = 0, - ( 1 == astSscanf( setting, "tolinverse= %lg %n", &dval, &nc ) ) - && ( nc >= len ) ) { - astSetTolInverse( this, dval ); - -/* If the attribute is still not recognised, pass it on to the parent - method for further interpretation. */ - } else { - (*parent_setattrib)( this_object, setting, status ); - } -} - -static void StoreArrays( AstPolyMap *this, int forward, int ncoeff, - const double *coeff, int *status ){ -/* -* Name: -* StoreArrays - -* Purpose: -* Store the dynamic arrays for a single transformation within a PolyMap - -* Type: -* Private function. - -* Synopsis: -* #include "polymap.h" -* void StoreArrays( AstPolyMap *this, int forward, int ncoeff, -* const double *coeff, int *status ) - -* Class Membership: -* PolyMap initialiser. - -* Description: -* This function sets up the arrays within a PolyMap structure that -* describes either the forward or inverse transformation. - -* Parameters: -* this -* The PolyMap. -* forward -* If non-zero, replace the forward transformation. Otherwise, -* replace the inverse transformation. -* ncoeff -* The number of non-zero coefficients necessary to define the -* specified transformation of the PolyMap. If zero is supplied, the -* transformation will be undefined. -* coeff -* An array containing "ncof*( 2 + nin )" elements. Each group -* of "2 + nin" adjacent elements describe a single coefficient of -* the transformation. Within each such group, the first -* element is the coefficient value; the next element is the -* integer index of the PolyMap output which uses the coefficient -* within its defining polynomial (the first output has index 1); -* the remaining elements of the group give the integer powers to -* use with each input coordinate value (powers must not be -* negative). -* status -* Pointer to inherited status. -*/ - -/* Local Variables: */ - const double *group; /* Pointer to start of next coeff. description */ - int *pows; /* Pointer to powers for current coeff. */ - int gsize; /* Length of each coeff. description */ - int i; /* Loop count */ - int ico; /* Index of next coeff. for current input or output */ - int iin; /* Input index extracted from coeff. description */ - int iout; /* Output index extracted from coeff. description */ - int j; /* Loop count */ - int nin; /* Number of inputs */ - int nout; /* Number of outputs */ - int pow; /* Power extracted from coeff. description */ - -/* Check the global status. */ - if ( !astOK ) return; - -/* Get the number of inputs and outputs. */ - nin = astGetNin( this ); - nout = astGetNout( this ); - -/* First Free any existing arrays. */ - FreeArrays( this, forward, status ); - -/* Now initialise the forward transformation, if required. */ - if( forward && ncoeff > 0 ) { - -/* Create the arrays decribing the forward transformation. */ - this->ncoeff_f = astMalloc( sizeof( int )*(size_t) nout ); - this->mxpow_f = astMalloc( sizeof( int )*(size_t) nin ); - this->power_f = astMalloc( sizeof( int ** )*(size_t) nout ); - this->coeff_f = astMalloc( sizeof( double * )*(size_t) nout ); - if( astOK ) { - -/* Initialise the count of coefficients for each output coordinate to zero. */ - for( i = 0; i < nout; i++ ) this->ncoeff_f[ i ] = 0; - -/* Initialise max power for each input coordinate to zero. */ - for( j = 0; j < nin; j++ ) this->mxpow_f[ j ] = 0; - -/* Scan through the supplied forward coefficient array, counting the - number of coefficients which relate to each output. Also find the - highest power used for each input axis. Report errors if any unusable - values are found in the supplied array. */ - group = coeff; - gsize = 2 + nin; - for( i = 0; i < ncoeff && astOK; i++, group += gsize ) { - - iout = floor( group[ 1 ] + 0.5 ); - if( iout < 1 || iout > nout ) { - astError( AST__BADCI, "astInitPolyMap(%s): Forward " - "coefficient %d referred to an illegal output " - "coordinate %d.", status, astGetClass( this ), i + 1, - iout ); - astError( AST__BADCI, "This number should be in the " - "range 1 to %d.", status, nout ); - break; - } - - this->ncoeff_f[ iout - 1 ]++; - - for( j = 0; j < nin; j++ ) { - pow = floor( group[ 2 + j ] + 0.5 ); - if( pow < 0 ) { - astError( AST__BADPW, "astInitPolyMap(%s): Forward " - "coefficient %d has a negative power (%d) " - "for input coordinate %d.", status, - astGetClass( this ), i + 1, pow, j + 1 ); - astError( AST__BADPW, "All powers should be zero or " - "positive." , status); - break; - } - if( pow > this->mxpow_f[ j ] ) this->mxpow_f[ j ] = pow; - } - } - -/* Allocate the arrays to store the input powers associated with each - coefficient, and the coefficient values. Reset the coefficient count - for each axis to zero afterwards so that we can use the array as an index - to the next vacant slot withint he following loop. */ - for( i = 0; i < nout; i++ ) { - this->power_f[ i ] = astMalloc( sizeof( int * )* - (size_t) this->ncoeff_f[ i ] ); - this->coeff_f[ i ] = astMalloc( sizeof( double )* - (size_t) this->ncoeff_f[ i ] ); - this->ncoeff_f[ i ] = 0; - } - - if( astOK ) { - -/* Extract the coefficient values and powers form the supplied array and - store them in the arrays created above. */ - group = coeff; - for( i = 0; i < ncoeff && astOK; i++, group += gsize ) { - iout = floor( group[ 1 ] + 0.5 ) - 1; - ico = ( this->ncoeff_f[ iout ] )++; - this->coeff_f[ iout ][ ico ] = group[ 0 ]; - - pows = astMalloc( sizeof( int )*(size_t) nin ); - this->power_f[ iout ][ ico ] = pows; - if( astOK ) { - for( j = 0; j < nin; j++ ) { - pows[ j ] = floor( group[ 2 + j ] + 0.5 ); - } - } - } - } - } - } - -/* Now initialise the inverse transformation, if required. */ - if( !forward && ncoeff > 0 ) { - -/* Create the arrays decribing the inverse transformation. */ - this->ncoeff_i = astMalloc( sizeof( int )*(size_t) nin ); - this->mxpow_i = astMalloc( sizeof( int )*(size_t) nout ); - this->power_i = astMalloc( sizeof( int ** )*(size_t) nin ); - this->coeff_i = astMalloc( sizeof( double * )*(size_t) nin ); - if( astOK ) { - -/* Initialise the count of coefficients for each input coordinate to zero. */ - for( i = 0; i < nin; i++ ) this->ncoeff_i[ i ] = 0; - -/* Initialise max power for each output coordinate to zero. */ - for( j = 0; j < nout; j++ ) this->mxpow_i[ j ] = 0; - -/* Scan through the supplied inverse coefficient array, counting the - number of coefficients which relate to each input. Also find the - highest power used for each output axis. Report errors if any unusable - values are found in the supplied array. */ - group = coeff; - - gsize = 2 + nout; - for( i = 0; i < ncoeff && astOK; i++, group += gsize ) { - - iin = floor( group[ 1 ] + 0.5 ); - if( iin < 1 || iin > nin ) { - astError( AST__BADCI, "astInitPolyMap(%s): Inverse " - "coefficient %d referred to an illegal input " - "coordinate %d.", status, astGetClass( this ), - i + 1, iin ); - astError( AST__BADCI, "This number should be in the " - "range 1 to %d.", status, nin ); - break; - } - - this->ncoeff_i[ iin - 1 ]++; - - for( j = 0; j < nout; j++ ) { - pow = floor( group[ 2 + j ] + 0.5 ); - if( pow < 0 ) { - astError( AST__BADPW, "astInitPolyMap(%s): Inverse " - "coefficient %d has a negative power (%d) " - "for output coordinate %d.", status, - astGetClass( this ), i + 1, pow, j + 1 ); - astError( AST__BADPW, "All powers should be zero or " - "positive." , status); - break; - } - if( pow > this->mxpow_i[ j ] ) this->mxpow_i[ j ] = pow; - } - } - -/* Allocate the arrays to store the output powers associated with each - coefficient, and the coefficient values. Reset the coefficient count - for each axis to zero afterwards so that we can use the array as an index - to the next vacant slot within the following loop. */ - for( i = 0; i < nin; i++ ) { - this->power_i[ i ] = astMalloc( sizeof( int * )* - (size_t) this->ncoeff_i[ i ] ); - this->coeff_i[ i ] = astMalloc( sizeof( double )* - (size_t) this->ncoeff_i[ i ] ); - this->ncoeff_i[ i ] = 0; - } - - if( astOK ) { - -/* Extract the coefficient values and powers form the supplied array and - store them in the arrays created above. */ - group = coeff; - for( i = 0; i < ncoeff && astOK; i++, group += gsize ) { - iin = floor( group[ 1 ] + 0.5 ) - 1; - ico = ( this->ncoeff_i[ iin ] )++; - this->coeff_i[ iin ][ ico ] = group[ 0 ]; - - pows = astMalloc( sizeof( int )*(size_t) nout ); - this->power_i[ iin ][ ico ] = pows; - if( astOK ) { - for( j = 0; j < nout; j++ ) { - pows[ j ] = floor( group[ 2 + j ] + 0.5 ); - } - } - } - } - } - } -} - -static int TestAttrib( AstObject *this_object, const char *attrib, int *status ) { -/* -* Name: -* TestAttrib - -* Purpose: -* Test if a specified attribute value is set for a PolyMap. - -* Type: -* Private function. - -* Synopsis: -* #include "polymap.h" -* int TestAttrib( AstObject *this, const char *attrib, int *status ) - -* Class Membership: -* PolyMap member function (over-rides the astTestAttrib protected -* method inherited from the Mapping class). - -* Description: -* This function returns a boolean result (0 or 1) to indicate whether -* a value has been set for one of a PolyMap's attributes. - -* Parameters: -* this -* Pointer to the PolyMap. -* attrib -* Pointer to a null-terminated string specifying the attribute -* name. This should be in lower case with no surrounding white -* space. -* status -* Pointer to the inherited status variable. - -* Returned Value: -* One if a value has been set, otherwise zero. - -* 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: */ - AstPolyMap *this; /* Pointer to the PolyMap structure */ - int result; /* Result value to return */ - -/* Initialise. */ - result = 0; - -/* Check the global error status. */ - if ( !astOK ) return result; - -/* Obtain a pointer to the PolyMap structure. */ - this = (AstPolyMap *) this_object; - -/* Check the attribute name and test the appropriate attribute. */ - -/* IterInverse. */ -/* ------------ */ - if ( !strcmp( attrib, "iterinverse" ) ) { - result = astTestIterInverse( this ); - -/* NiterInverse. */ -/* ------------- */ - } else if ( !strcmp( attrib, "niterinverse" ) ) { - result = astTestNiterInverse( this ); - -/* TolInverse. */ -/* ----------- */ - } else if ( !strcmp( attrib, "tolinverse" ) ) { - result = astTestTolInverse( this ); - -/* If the attribute is still not recognised, pass it on to the parent - method for further interpretation. */ - } else { - result = (*parent_testattrib)( this_object, attrib, status ); - } - -/* Return the result, */ - return result; -} - -static AstPointSet *Transform( AstMapping *this, AstPointSet *in, - int forward, AstPointSet *out, int *status ) { -/* -* Name: -* Transform - -* Purpose: -* Apply a PolyMap to transform a set of points. - -* Type: -* Private function. - -* Synopsis: -* #include "polymap.h" -* AstPointSet *Transform( AstMapping *this, AstPointSet *in, -* int forward, AstPointSet *out, int *status ) - -* Class Membership: -* PolyMap member function (over-rides the astTransform protected -* method inherited from the Mapping class). - -* Description: -* This function takes a PolyMap and a set of points encapsulated in a -* PointSet and transforms the points. - -* Parameters: -* this -* Pointer to the PolyMap. -* 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 PolyMap being applied. -* - The number of coordinate values per point in the output PointSet will -* equal the number of rows in the PolyMap 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 */ - AstPolyMap *map; /* Pointer to PolyMap to be applied */ - double **coeff; /* Pointer to coefficient value arrays */ - double **ptr_in; /* Pointer to input coordinate data */ - double **ptr_out; /* Pointer to output coordinate data */ - double **work; /* Pointer to exponentiated axis values */ - double *outcof; /* Pointer to next coefficient value */ - double outval; /* Output axis value */ - double term; /* Term to be added to output value */ - double xp; /* Exponentiated input axis value */ - int ***power; /* Pointer to coefficient power arrays */ - int **outpow; /* Pointer to next set of axis powers */ - int *mxpow; /* Pointer to max used power for each input */ - int *ncoeff; /* Pointer to no. of coefficients */ - int in_coord; /* Index of output coordinate */ - int ico; /* Coefficient index */ - int nc; /* No. of coefficients in polynomial */ - 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 */ - int pow; /* Next axis power */ - -/* Check the global error status. */ - if ( !astOK ) return NULL; - -/* Obtain a pointer to the PolyMap. */ - map = (AstPolyMap *) 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 ); - -/* Determine whether to apply the original forward or inverse mapping, - according to the direction specified and whether the mapping has been - inverted. */ - if ( astGetInvert( map ) ) forward = !forward; - -/* We will now extend the parent astTransform method by performing the - calculations needed to generate the output coordinate values. */ - -/* If we are using the original inverse transformatiom, and the IterInverse - attribute is non-zero, use an iterative inverse algorithm rather than any - inverse transformation defined within the PolyMap. */ - if( !forward && astGetIterInverse(map) ) { - IterInverse( map, in, result, status ); - -/* Otherwise, 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. */ - } else { - ncoord_in = astGetNcoord( in ); - ncoord_out = astGetNcoord( result ); - npoint = astGetNpoint( in ); - ptr_in = astGetPoints( in ); - ptr_out = astGetPoints( result ); - -/* Get a pointer to the arrays holding the required coefficient - values and powers, according to the direction of mapping required. */ - if ( forward ) { - ncoeff = map->ncoeff_f; - coeff = map->coeff_f; - power = map->power_f; - mxpow = map->mxpow_f; - } else { - ncoeff = map->ncoeff_i; - coeff = map->coeff_i; - power = map->power_i; - mxpow = map->mxpow_i; - } - -/* Allocate memory to hold the required powers of the input axis values. */ - work = astMalloc( sizeof( double * )*(size_t) ncoord_in ); - for( in_coord = 0; in_coord < ncoord_in; in_coord++ ) { - work[ in_coord ] = astMalloc( sizeof( double )* - (size_t) ( astMAX( 2, mxpow[in_coord]+1 ) ) ); - } - -/* Perform coordinate arithmetic. */ -/* ------------------------------ */ - if ( astOK ) { - -/* Loop to apply the polynomial to each point in turn.*/ - for ( point = 0; point < npoint; point++ ) { - -/* Find the required powers of the input axis values and store them - in the work array. Note, using a virtual method here slows the PolyMap - Transform function down by about 5%, compared to doing the equivalent - calculations in-line. But we need some way to allow the ChebyMap class - to over-ride the calculation of the powers, so we must do something - like this. If the 5% slow-down is too much, it can be reduced down to - about 2% by replacing the invocation of the astPolyPowers_ interface - function with a direct call to the implementation function itself. - This involves replacing the astPolyPowers call below with this: - - (**astMEMBER(this,PolyMap,PolyPowers))( (AstPolyMap *) this, work, ncoord_in, - mxpow, ptr_in, point, forward, status ); - - The above could be wrapped up in an alternative implementation of the - astPolyPowers macro, so that it looks the same as the existing code. - In fact, this scheme could be more widely used to speed up invocation - of virtual functions within AST. The disadvantage is that the interface - functions for some virtual methods includes some extra processing, - over and above simply invoking the implementation function. - - Of course the other way to get rid of the 5% slow down, is to - revert to using in-line code below, and then replicate this entire - function in the ChebyMap class, making suitable changes to use - Chebyshev functions in place of simple powers. But that is bad - structuring... */ - astPolyPowers( this, work, ncoord_in, mxpow, ptr_in, point, - forward ); - -/* Loop round each output. */ - for( out_coord = 0; out_coord < ncoord_out; out_coord++ ) { - -/* Initialise the output value. */ - outval = 0.0; - -/* Get pointers to the coefficients and powers for this output. */ - outcof = coeff[ out_coord ]; - outpow = power[ out_coord ]; - -/* Loop round all polynomial coefficients.*/ - nc = ncoeff[ out_coord ]; - for ( ico = 0; ico < nc && outval != AST__BAD; - ico++, outcof++, outpow++ ) { - -/* Initialise the current term to be equal to the value of the coefficient. - If it is bad, store a bad output value. */ - term = *outcof; - if( term == AST__BAD ) { - outval = AST__BAD; - -/* Otherwise, loop round all inputs */ - } else { - for( in_coord = 0; in_coord < ncoord_in; in_coord++ ) { - -/* Get the power of the current input axis value used by the current - coefficient. If it is zero, pass on. */ - pow = (*outpow)[ in_coord ]; - if( pow > 0 ) { - -/* Get the axis value raised to the appropriate power. */ - xp = work[ in_coord ][ pow ]; - -/* If bad, set the output value bad and break. */ - if( xp == AST__BAD ) { - outval = AST__BAD; - break; - -/* Otherwise multiply the current term by the exponentiated axis value. */ - } else { - term *= xp; - } - } - } - } - -/* Increment the output value by the current term of the polynomial. */ - if( outval != AST__BAD ) outval += term; - - } - -/* Store the output value. */ - ptr_out[ out_coord ][ point ] = outval; - - } - } - } - -/* Free resources. */ - for( in_coord = 0; in_coord < ncoord_in; in_coord++ ) { - work[ in_coord ] = astFree( work[ in_coord ] ); - } - work = astFree( work ); - } - -/* Return a pointer to the output PointSet. */ - 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). */ - -/* -*att++ -* Name: -* IterInverse - -* Purpose: -* Provide an iterative inverse transformation? - -* Type: -* Public attribute. - -* Synopsis: -* Integer (boolean). - -* Description: -* This attribute indicates whether the inverse transformation of -* the PolyMap should be implemented via an iterative Newton-Raphson -* approximation that uses the forward transformation to transform -* candidate input positions until an output position is found which -* is close to the required output position. By default, an iterative -* inverse is provided if, and only if, no inverse polynomial was supplied -* when the PolyMap was constructed. -* -* The NiterInverse and TolInverse attributes provide parameters that -* control the behaviour of the inverse approximation method. -* -* The iterative inverse returns AST__BAD axis values at positions -* for which the inverse transformation is undefined. For instance, -* if the forward transformation is y = x*x, the iterative inverse -* will return x = AST__BAD at y = -1. If the inverse transformation -* is multiply defined, the position returned by the iterative inverse -* will be the position of the solution that is closest to the -* supplied position. For instance, using the above example, y = x*x, -* the iterative inverse will return x = +2 at y = 4, because x = +2 -* is the closest solution to 4 (the other solution is x = -2). - -* Applicability: -* PolyMap -* All PolyMaps have this attribute. -* ChebyMap -* The ChebyMap class does not currently provide an option for an -* iterative inverse, and so the IterInverse value is always zero. -* Setting or clearing the IterInverse attribute of a ChebyMap has -* no effect. - -* Notes: -* - The transformation replaced by the iterative algorithm is the -* transformation from the original PolyMap output space to the -* original PolyMap input space (i.e. the input and output spaces as -* defined by the arguments of the PolyMap constructor). This is still -* the case even if the PolyMap has subsequently been inverted. In -* other words if a PolyMap is created and then inverted, setting -* the IterInverse to a non-zero value will replace the forward -* transformation of the inverted PolyMap (i.e. the inverse -* transformation of the original PolyMap). It is not possible to -* replace the other transformation (i.e. from the original PolyMap -* input space to the original PolyMap output space) with an iterative -* algorithm. -* - If a PolyMap that has an iterative inverse transformation is -* subsequently inverted, the inverted PolyMap will have an iterative -* forward transformation. -* - An iterative inverse can only be used if the PolyMap has equal -* numbers of inputs and outputs, as given by the Nin and Nout -* attributes. An error will be reported if IterInverse is set non-zero -* for a PolyMap that does not meet this requirement. - -*att-- -*/ -astMAKE_CLEAR(PolyMap,IterInverse,iterinverse,-INT_MAX) -astMAKE_GET(PolyMap,IterInverse,int,0,( ( this->iterinverse == -INT_MAX ) ? - (this->ncoeff_i == 0) : this->iterinverse )) -astMAKE_SET(PolyMap,IterInverse,int,iterinverse, - (((astGetNin(this)==astGetNout(this))||!value)?((value?1:0)):(astError(AST__ATTIN,"astSetIterInverse(%s):" - "Cannot use an iterative inverse because the %s has unequal numbers of " - "inputs and outputs.", status, astGetClass(this),astGetClass(this)),this->iterinverse))) -astMAKE_TEST(PolyMap,IterInverse,( this->iterinverse != -INT_MAX )) - -/* NiterInverse. */ -/* --------- */ -/* -*att++ -* Name: -* NiterInverse - -* Purpose: -* Maximum number of iterations for the iterative inverse transformation. - -* Type: -* Public attribute. - -* Synopsis: -* Integer. - -* Description: -* This attribute controls the iterative inverse transformation -* used if the IterInverse attribute is non-zero. -* -* Its value gives the maximum number of iterations of the -* Newton-Raphson algorithm to be used for each transformed position. -* The default value is 4. See also attribute TolInverse. - -* Applicability: -* PolyMap -* All PolyMaps have this attribute. - -*att-- -*/ -astMAKE_CLEAR(PolyMap,NiterInverse,niterinverse,-INT_MAX) -astMAKE_GET(PolyMap,NiterInverse,int,0,( this->niterinverse == -INT_MAX ? 4 : this->niterinverse)) -astMAKE_SET(PolyMap,NiterInverse,int,niterinverse,value) -astMAKE_TEST(PolyMap,NiterInverse,( this->niterinverse != -INT_MAX )) - -/* TolInverse. */ -/* ----------- */ -/* -*att++ -* Name: -* TolInverse - -* Purpose: -* Target relative error for the iterative inverse transformation. - -* Type: -* Public attribute. - -* Synopsis: -* Floating point. - -* Description: -* This attribute controls the iterative inverse transformation -* used if the IterInverse attribute is non-zero. -* -* Its value gives the target relative error in the axis values of -* each transformed position. Further iterations will be performed -* until the target relative error is reached, or the maximum number -* of iterations given by attribute NiterInverse is reached. - -* The default value is 1.0E-6. - -* Applicability: -* PolyMap -* All PolyMaps have this attribute. -*att-- -*/ -astMAKE_CLEAR(PolyMap,TolInverse,tolinverse,AST__BAD) -astMAKE_GET(PolyMap,TolInverse,double,0.0,( this->tolinverse == AST__BAD ? 1.0E-6 : this->tolinverse)) -astMAKE_SET(PolyMap,TolInverse,double,tolinverse,value) -astMAKE_TEST(PolyMap,TolInverse,( this->tolinverse != AST__BAD )) - -/* Copy constructor. */ -/* ----------------- */ -static void Copy( const AstObject *objin, AstObject *objout, int *status ) { -/* -* Name: -* Copy - -* Purpose: -* Copy constructor for PolyMap objects. - -* Type: -* Private function. - -* Synopsis: -* void Copy( const AstObject *objin, AstObject *objout, int *status ) - -* Description: -* This function implements the copy constructor for PolyMap 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 -* coefficients associated with the input PolyMap. -*/ - - -/* Local Variables: */ - AstPolyMap *in; /* Pointer to input PolyMap */ - AstPolyMap *out; /* Pointer to output PolyMap */ - int nin; /* No. of input coordinates */ - int nout; /* No. of output coordinates */ - int i; /* Loop count */ - int j; /* Loop count */ - -/* Check the global error status. */ - if ( !astOK ) return; - -/* Obtain pointers to the input and output PolyMaps. */ - in = (AstPolyMap *) objin; - out = (AstPolyMap *) 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->ncoeff_f = NULL; - out->power_f = NULL; - out->coeff_f = NULL; - out->mxpow_f = NULL; - - out->ncoeff_i = NULL; - out->power_i = NULL; - out->coeff_i = NULL; - out->mxpow_i = NULL; - - out->jacobian = NULL; - -/* Get the number of inputs and outputs of the uninverted Mapping. */ - nin = ( (AstMapping *) in )->nin; - nout = ( (AstMapping *) in )->nout; - -/* Copy the number of coefficients associated with each output of the forward - transformation. */ - if( in->ncoeff_f ) { - out->ncoeff_f = (int *) astStore( NULL, (void *) in->ncoeff_f, - sizeof( int )*(size_t) nout ); - -/* Copy the maximum power of each input axis value used by the forward - transformation. */ - out->mxpow_f = (int *) astStore( NULL, (void *) in->mxpow_f, - sizeof( int )*(size_t) nin ); - -/* Copy the coefficient values used by the forward transformation. */ - if( in->coeff_f ) { - out->coeff_f = astMalloc( sizeof( double * )*(size_t) nout ); - if( astOK ) { - for( i = 0; i < nout; i++ ) { - out->coeff_f[ i ] = (double *) astStore( NULL, (void *) in->coeff_f[ i ], - sizeof( double )*(size_t) in->ncoeff_f[ i ] ); - } - } - } - -/* Copy the input axis powers associated with each coefficient of the forward - transformation. */ - if( in->power_f ) { - out->power_f = astMalloc( sizeof( int ** )*(size_t) nout ); - if( astOK ) { - for( i = 0; i < nout; i++ ) { - out->power_f[ i ] = astMalloc( sizeof( int * )*(size_t) in->ncoeff_f[ i ] ); - if( astOK ) { - for( j = 0; j < in->ncoeff_f[ i ]; j++ ) { - out->power_f[ i ][ j ] = (int *) astStore( NULL, (void *) in->power_f[ i ][ j ], - sizeof( int )*(size_t) nin ); - } - } - } - } - } - } - -/* Do the same for the inverse transformation. */ - if( in->ncoeff_i ) { - out->ncoeff_i = (int *) astStore( NULL, (void *) in->ncoeff_i, - sizeof( int )*(size_t) nin ); - - out->mxpow_i = (int *) astStore( NULL, (void *) in->mxpow_i, - sizeof( int )*(size_t) nout ); - - if( in->coeff_i ) { - out->coeff_i = astMalloc( sizeof( double * )*(size_t) nin ); - if( astOK ) { - for( i = 0; i < nin; i++ ) { - out->coeff_i[ i ] = (double *) astStore( NULL, (void *) in->coeff_i[ i ], - sizeof( double )*(size_t) in->ncoeff_i[ i ] ); - } - } - } - - if( in->power_i ) { - out->power_i = astMalloc( sizeof( int ** )*(size_t) nin ); - if( astOK ) { - for( i = 0; i < nin; i++ ) { - out->power_i[ i ] = astMalloc( sizeof( int * )*(size_t) in->ncoeff_i[ i ] ); - if( astOK ) { - for( j = 0; j < in->ncoeff_i[ i ]; j++ ) { - out->power_i[ i ][ j ] = (int *) astStore( NULL, (void *) in->power_i[ i ][ j ], - sizeof( int )*(size_t) nout ); - } - } - } - } - } - } - -/* If an error has occurred, free al the resources allocated above. */ - if( !astOK ) { - FreeArrays( out, 1, status ); - FreeArrays( out, 0, status ); - } - - return; - -} - -/* Destructor. */ -/* ----------- */ -static void Delete( AstObject *obj, int *status ) { -/* -* Name: -* Delete - -* Purpose: -* Destructor for PolyMap objects. - -* Type: -* Private function. - -* Synopsis: -* void Delete( AstObject *obj, int *status ) - -* Description: -* This function implements the destructor for PolyMap 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: */ - AstPolyMap *this; - int nc; - int ic; - int lstat; - int error; - -/* Obtain a pointer to the PolyMap structure. */ - this = (AstPolyMap *) obj; - -/* Free the arrays. */ - FreeArrays( this, 1, status ); - FreeArrays( this, 0, status ); - -/* Free the resources used to store the Jacobian of the forward - transformation. */ - if( this->jacobian ) { - -/* Get the number of PolyMap inputs. We need to clear any error status - first since astGetNin returns zero if an error has occurred. The - Jacobian will only be non-NULL if the number of inputs and outputs - are equal. */ - error = !astOK; - if( error ) { - lstat = astStatus; - astClearStatus; - } - nc = astGetNin( this ); - if( error ) astSetStatus( lstat ); - - for( ic = 0; ic < nc; ic++ ) { - (this->jacobian)[ ic ] = astAnnul( (this->jacobian)[ ic ] ); - } - this->jacobian = astFree( this->jacobian ); - } -} - -/* Dump function. */ -/* -------------- */ -static void Dump( AstObject *this_object, AstChannel *channel, int *status ) { -/* -* Name: -* Dump - -* Purpose: -* Dump function for PolyMap 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 PolyMap class to an output Channel. - -* Parameters: -* this -* Pointer to the PolyMap 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: */ - AstPolyMap *this; /* Pointer to the PolyMap structure */ - char buff[ KEY_LEN + 1 ]; /* Buffer for keyword string */ - char comm[ 100 ]; /* Buffer for comment string */ - double dval; /* Floating point attribute value */ - int i; /* Loop index */ - int iv; /* Vectorised keyword index */ - int ival; /* Integer value */ - int j; /* Loop index */ - int k; /* Loop index */ - int nin; /* No. of input coords */ - int nout; /* No. of output coords */ - int set; /* Attribute value set? */ - -/* Check the global error status. */ - if ( !astOK ) return; - -/* Obtain a pointer to the PolyMap structure. */ - this = (AstPolyMap *) this_object; - -/* Find the number of inputs and outputs of the uninverted Mapping. */ - nin = ( (AstMapping *) this )->nin; - nout = ( (AstMapping *) this )->nout; - -/* Write out values representing the instance variables for the - PolyMap class. */ - -/* First do the forward transformation arrays. Check they are used. */ - if( this->ncoeff_f ) { - -/* Store the maximum power of each input axis value used by the forward - transformation. */ - for( i = 0; i < nin; i++ ){ - (void) sprintf( buff, "MPF%d", i + 1 ); - (void) sprintf( comm, "Max. power of input %d in any forward polynomial", i + 1 ); - astWriteInt( channel, buff, 1, 1, (this->mxpow_f)[ i ], comm ); - } - -/* Store the number of coefficients associated with each output of the forward - transformation. */ - for( i = 0; i < nout; i++ ){ - (void) sprintf( buff, "NCF%d", i + 1 ); - (void) sprintf( comm, "No. of coeff.s for forward polynomial %d", i + 1 ); - astWriteInt( channel, buff, 1, 1, (this->ncoeff_f)[ i ], comm ); - } - -/* Store the coefficient values used by the forward transformation. */ - iv = 1; - for( i = 0; i < nout; i++ ){ - for( j = 0; j < this->ncoeff_f[ i ]; j++, iv++ ){ - if( (this->coeff_f)[ i ][ j ] != AST__BAD ) { - (void) sprintf( buff, "CF%d", iv ); - (void) sprintf( comm, "Coeff %d of forward polynomial %d", j + 1, i + 1 ); - astWriteDouble( channel, buff, 1, 1, (this->coeff_f)[ i ][ j ], comm ); - } - } - } - -/* Store the input axis powers associated with each coefficient of the forward - transformation. */ - iv = 1; - for( i = 0; i < nout; i++ ){ - for( j = 0; j < this->ncoeff_f[ i ]; j++ ){ - for( k = 0; k < nin; k++, iv++ ){ - if( (this->power_f)[ i ][ j ][ k ] > 0 ) { - (void) sprintf( buff, "PF%d", iv ); - (void) sprintf( comm, "Power of i/p %d for coeff %d of fwd poly %d", k + 1, j + 1, i + 1 ); - astWriteDouble( channel, buff, 1, 1, (this->power_f)[ i ][ j ][ k ], comm ); - } - } - } - } - } - -/* Now do the inverse transformation arrays. Check they are used. */ - if( this->ncoeff_i ) { - -/* Store the maximum power of each output axis value used by the inverse - transformation. */ - for( i = 0; i < nout; i++ ){ - (void) sprintf( buff, "MPI%d", i + 1 ); - (void) sprintf( comm, "Max. power of output %d in any inverse polynomial", i + 1 ); - astWriteInt( channel, buff, 1, 1, (this->mxpow_i)[ i ], comm ); - } - -/* Store the number of coefficients associated with each input of the inverse - transformation. */ - for( i = 0; i < nin; i++ ){ - (void) sprintf( buff, "NCI%d", i + 1 ); - (void) sprintf( comm, "No. of coeff.s for inverse polynomial %d", i + 1 ); - astWriteInt( channel, buff, 1, 1, (this->ncoeff_i)[ i ], comm ); - } - -/* Store the coefficient values used by the inverse transformation. */ - iv = 1; - for( i = 0; i < nin; i++ ){ - for( j = 0; j < this->ncoeff_i[ i ]; j++, iv++ ){ - if( (this->coeff_i)[ i ][ j ] != AST__BAD ) { - (void) sprintf( buff, "CI%d", iv ); - (void) sprintf( comm, "Coeff %d of inverse polynomial %d", j + 1, i + 1 ); - astWriteDouble( channel, buff, 1, 1, (this->coeff_i)[ i ][ j ], comm ); - } - } - } - -/* Store the output axis powers associated with each coefficient of the inverse - transformation. */ - iv = 1; - for( i = 0; i < nin; i++ ){ - for( j = 0; j < this->ncoeff_i[ i ]; j++ ){ - for( k = 0; k < nout; k++, iv++ ){ - if( (this->power_i)[ i ][ j ][ k ] > 0 ) { - (void) sprintf( buff, "PI%d", iv ); - (void) sprintf( comm, "Power of o/p %d for coeff %d of inv poly %d", k + 1, j + 1, i + 1 ); - astWriteDouble( channel, buff, 1, 1, (this->power_i)[ i ][ j ][ k ], comm ); - } - } - } - } - } - -/* Use an iterative inverse? */ - set = TestIterInverse( this, status ); - ival = set ? GetIterInverse( this, status ) : astGetIterInverse( this ); - astWriteInt( channel, "IterInv", set, 0, ival, ival ? "Use an iterative inverse" : "Do not use an iterative inverse" ); - -/* Max number of iterations for iterative inverse. */ - set = TestNiterInverse( this, status ); - ival = set ? GetNiterInverse( this, status ) : astGetNiterInverse( this ); - astWriteInt( channel, "NiterInv", set, 0, ival, "Max number of iterations for iterative inverse" ); - -/* Target relative error for iterative inverse. */ - set = TestTolInverse( this, status ); - dval = set ? GetTolInverse( this, status ) : astGetTolInverse( this ); - astWriteDouble( channel, "TolInv", set, 0, dval, "Target relative error for iterative inverse" ); - -/* Undefine macros local to this function. */ -#undef KEY_LEN -} - -/* Standard class functions. */ -/* ========================= */ -/* Implement the astIsAPolyMap and astCheckPolyMap functions using the macros - defined for this purpose in the "object.h" header file. */ -astMAKE_ISA(PolyMap,Mapping) -astMAKE_CHECK(PolyMap) - -AstPolyMap *astPolyMap_( int nin, int nout, int ncoeff_f, const double coeff_f[], - int ncoeff_i, const double coeff_i[], const char *options, int *status, ...){ -/* -*++ -* Name: -c astPolyMap -f AST_POLYMAP - -* Purpose: -* Create a PolyMap. - -* Type: -* Public function. - -* Synopsis: -c #include "polymap.h" -c AstPolyMap *astPolyMap( int nin, int nout, int ncoeff_f, const double coeff_f[], -c int ncoeff_i, const double coeff_i[], -c const char *options, ... ) -f RESULT = AST_POLYMAP( NIN, NOUT, NCOEFF_F, COEFF_F, NCOEFF_I, COEFF_I, -f OPTIONS, STATUS ) - -* Class Membership: -* PolyMap constructor. - -* Description: -* This function creates a new PolyMap and optionally initialises -* its attributes. -* -* A PolyMap is a form of Mapping which performs a general polynomial -* transformation. Each output coordinate is a polynomial function of -* all the input coordinates. The coefficients are specified separately -* for each output coordinate. The forward and inverse transformations -* are defined independantly by separate sets of coefficients. If no -* inverse transformation is supplied, the default behaviour is to use -* an iterative method to evaluate the inverse based only on the forward -* transformation (see attribute IterInverse). - -* Parameters: -c nin -f NIN = INTEGER (Given) -* The number of input coordinates. -c nout -f NOUT = INTEGER (Given) -* The number of output coordinates. -c ncoeff_f -f NCOEFF_F = INTEGER (Given) -* The number of non-zero coefficients necessary to define the -* forward transformation of the PolyMap. If zero is supplied, the -* forward transformation will be undefined. -c coeff_f -f COEFF_F( * ) = DOUBLE PRECISION (Given) -* An array containing -c "ncoeff_f*( 2 + nin )" elements. Each group of "2 + nin" -f "NCOEFF_F*( 2 + NIN )" elements. Each group of "2 + NIN" -* adjacent elements describe a single coefficient of the forward -* transformation. Within each such group, the first element is the -* coefficient value; the next element is the integer index of the -* PolyMap output which uses the coefficient within its defining -* polynomial (the first output has index 1); the remaining elements -* of the group give the integer powers to use with each input -* coordinate value (powers must not be negative, and floating -* point values are rounded to the nearest integer). -c If "ncoeff_f" is zero, a NULL pointer may be supplied for "coeff_f". -* -* For instance, if the PolyMap has 3 inputs and 2 outputs, each group -* consisting of 5 elements, A groups such as "(1.2, 2.0, 1.0, 3.0, 0.0)" -* describes a coefficient with value 1.2 which is used within the -* definition of output 2. The output value is incremented by the -* product of the coefficient value, the value of input coordinate -* 1 raised to the power 1, and the value of input coordinate 2 raised -* to the power 3. Input coordinate 3 is not used since its power is -* specified as zero. As another example, the group "(-1.0, 1.0, -* 0.0, 0.0, 0.0 )" describes adds a constant value -1.0 onto -* output 1 (it is a constant value since the power for every input -* axis is given as zero). -* -c Each final output coordinate value is the sum of the "ncoeff_f" terms -c described by the "ncoeff_f" groups within the supplied array. -f Each final output coordinate value is the sum of the "NCOEFF_F" terms -f described by the "NCOEFF_F" groups within the supplied array. -c ncoeff_i -f NCOEFF_I = INTEGER (Given) -* The number of non-zero coefficients necessary to define the -* inverse transformation of the PolyMap. If zero is supplied, the -* default behaviour is to use an iterative method to evaluate the -* inverse based only on the forward transformation (see attribute -* IterInverse). -c coeff_i -f COEFF_I( * ) = DOUBLE PRECISION (Given) -* An array containing -c "ncoeff_i*( 2 + nout )" elements. Each group of "2 + nout" -f "NCOEFF_I*( 2 + NOUT )" elements. Each group of "2 + NOUT" -* adjacent elements describe a single coefficient of the inverse -c transformation, using the same schame as "coeff_f", -f transformation, using the same schame as "COEFF_F", -* except that "inputs" and "outputs" are transposed. -c If "ncoeff_i" is zero, a NULL pointer may be supplied for "coeff_i". -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 PolyMap. 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 PolyMap. 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 astPolyMap() -f AST_POLYMAP = INTEGER -* A pointer to the new PolyMap. - -* Notes: -* - 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. -*-- -*/ - -/* Local Variables: */ - astDECLARE_GLOBALS /* Pointer to thread-specific global data */ - AstPolyMap *new; /* Pointer to new PolyMap */ - 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 PolyMap, allocating memory and initialising the - virtual function table as well if necessary. */ - new = astInitPolyMap( NULL, sizeof( AstPolyMap ), !class_init, - &class_vtab, "PolyMap", nin, nout, - ncoeff_f, coeff_f, ncoeff_i, coeff_i ); - -/* 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 PolyMap'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 PolyMap. */ - return new; -} - -AstPolyMap *astPolyMapId_( int nin, int nout, int ncoeff_f, const double coeff_f[], - int ncoeff_i, const double coeff_i[], const char *options, ... ){ -/* -* Name: -* astPolyMapId_ - -* Purpose: -* Create a PolyMap. - -* Type: -* Private function. - -* Synopsis: -* #include "polymap.h" -* AstPolyMap *astPolyMap( int nin, int nout, int ncoeff_f, const double coeff_f[], -* int ncoeff_i, const double coeff_i[], -* const char *options, ... ) - -* Class Membership: -* PolyMap constructor. - -* Description: -* This function implements the external (public) interface to the -* astPolyMap constructor function. It returns an ID value (instead -* of a true C pointer) to external users, and must be provided -* because astPolyMap_ 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 astPolyMap_ 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 astPolyMap_. - -* Returned Value: -* The ID value associated with the new PolyMap. -*/ - -/* Local Variables: */ - astDECLARE_GLOBALS /* Pointer to thread-specific global data */ - AstPolyMap *new; /* Pointer to new PolyMap */ - 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 PolyMap, allocating memory and initialising the - virtual function table as well if necessary. */ - new = astInitPolyMap( NULL, sizeof( AstPolyMap ), !class_init, - &class_vtab, "PolyMap", nin, nout, - ncoeff_f, coeff_f, ncoeff_i, coeff_i ); - -/* 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 PolyMap'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 PolyMap. */ - return astMakeId( new ); -} - -AstPolyMap *astInitPolyMap_( void *mem, size_t size, int init, - AstPolyMapVtab *vtab, const char *name, - int nin, int nout, int ncoeff_f, const double coeff_f[], - int ncoeff_i, const double coeff_i[], int *status ){ -/* -*+ -* Name: -* astInitPolyMap - -* Purpose: -* Initialise a PolyMap. - -* Type: -* Protected function. - -* Synopsis: -* #include "polymap.h" -* AstPolyMap *astInitPolyMap( void *mem, size_t size, int init, -* AstPolyMapVtab *vtab, const char *name, -* int nin, int nout, int ncoeff_f, -* const double coeff_f[], int ncoeff_i, -* const double coeff_i[] ) - -* Class Membership: -* PolyMap initialiser. - -* Description: -* This function is provided for use by class implementations to initialise -* a new PolyMap object. It allocates memory (if necessary) to accommodate -* the PolyMap plus any additional data associated with the derived class. -* It then initialises a PolyMap 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 PolyMap at the start of the memory passed via the -* "vtab" parameter. - -* Parameters: -* mem -* A pointer to the memory in which the PolyMap is to be initialised. -* This must be of sufficient size to accommodate the PolyMap data -* (sizeof(PolyMap)) 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 PolyMap (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 PolyMap -* structure, so a valid value must be supplied even if not required for -* allocating memory. -* init -* A logical flag indicating if the PolyMap'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 PolyMap. -* 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. -* ncoeff_f -* The number of non-zero coefficients necessary to define the -* forward transformation of the PolyMap. If zero is supplied, the -* forward transformation will be undefined. -* coeff_f -* An array containing "ncoeff_f*( 2 + nin )" elements. Each group -* of "2 + nin" adjacent elements describe a single coefficient of -* the forward transformation. Within each such group, the first -* element is the coefficient value; the next element is the -* integer index of the PolyMap output which uses the coefficient -* within its defining polynomial (the first output has index 1); -* the remaining elements of the group give the integer powers to -* use with each input coordinate value (powers must not be -* negative) -* -* For instance, if the PolyMap has 3 inputs and 2 outputs, each group -* consisting of 5 elements, A groups such as "(1.2, 2.0, 1.0, 3.0, 0.0)" -* describes a coefficient with value 1.2 which is used within the -* definition of output 2. The output value is incremented by the -* product of the coefficient value, the value of input coordinate -* 1 raised to the power 1, and the value of input coordinate 2 raised -* to the power 3. Input coordinate 3 is not used since its power is -* specified as zero. As another example, the group "(-1.0, 1.0, -* 0.0, 0.0, 0.0 )" describes adds a constant value -1.0 onto -* output 1 (it is a constant value since the power for every input -* axis is given as zero). -* -* Each final output coordinate value is the sum of the "ncoeff_f" terms -* described by the "ncoeff_f" groups within the supplied array. -* ncoeff_i -* The number of non-zero coefficients necessary to define the -* inverse transformation of the PolyMap. If zero is supplied, the -* inverse transformation will be undefined. -* coeff_i -* An array containing -* "ncoeff_i*( 2 + nout )" elements. Each group of "2 + nout" -* adjacent elements describe a single coefficient of the inverse -* transformation, using the same schame as "coeff_f", except that -* "inputs" and "outputs" are transposed. - -* Returned Value: -* A pointer to the new PolyMap. - -* 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: */ - AstPolyMap *new; /* Pointer to new PolyMap */ - -/* Check the global status. */ - if ( !astOK ) return NULL; - -/* If necessary, initialise the virtual function table. */ - if ( init ) astInitPolyMapVtab( vtab, name ); - -/* Initialise a Mapping structure (the parent class) as the first component - within the PolyMap structure, allocating memory if necessary. Specify that - the Mapping should be defined in both the forward and inverse directions. */ - new = (AstPolyMap *) astInitMapping( mem, size, 0, - (AstMappingVtab *) vtab, name, - nin, nout, 1, 1 ); - if ( astOK ) { - -/* Initialise the PolyMap data. */ -/* ---------------------------- */ - -/* First initialise the pointers in case of errors. */ - new->ncoeff_f = NULL; - new->power_f = NULL; - new->coeff_f = NULL; - new->mxpow_f = NULL; - - new->ncoeff_i = NULL; - new->power_i = NULL; - new->coeff_i = NULL; - new->mxpow_i = NULL; - -/* Store the forward transformation. */ - StoreArrays( new, 1, ncoeff_f, coeff_f, status ); - -/* Store the inverse transformation. */ - StoreArrays( new, 0, ncoeff_i, coeff_i, status ); - -/* Other class attributes. */ - new->iterinverse = -INT_MAX; - new->niterinverse = -INT_MAX; - new->tolinverse = AST__BAD; - new->jacobian = NULL; - -/* If an error occurred, clean up by deleting the new PolyMap. */ - if ( !astOK ) new = astDelete( new ); - } - -/* Return a pointer to the new PolyMap. */ - return new; -} - -AstPolyMap *astLoadPolyMap_( void *mem, size_t size, - AstPolyMapVtab *vtab, const char *name, - AstChannel *channel, int *status ) { -/* -*+ -* Name: -* astLoadPolyMap - -* Purpose: -* Load a PolyMap. - -* Type: -* Protected function. - -* Synopsis: -* #include "polymap.h" -* AstPolyMap *astLoadPolyMap( void *mem, size_t size, -* AstPolyMapVtab *vtab, const char *name, -* AstChannel *channel ) - -* Class Membership: -* PolyMap loader. - -* Description: -* This function is provided to load a new PolyMap 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 -* PolyMap 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 PolyMap at the start of the memory -* passed via the "vtab" parameter. - - -* Parameters: -* mem -* A pointer to the memory into which the PolyMap is to be -* loaded. This must be of sufficient size to accommodate the -* PolyMap data (sizeof(PolyMap)) 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 PolyMap (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 PolyMap 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(AstPolyMap) is used instead. -* vtab -* Pointer to the start of the virtual function table to be -* associated with the new PolyMap. If this is NULL, a pointer -* to the (static) virtual function table for the PolyMap 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 "PolyMap" is used instead. - -* Returned Value: -* A pointer to the new PolyMap. - -* 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: */ - AstPolyMap *new; /* Pointer to the new PolyMap */ - char buff[ KEY_LEN + 1 ]; /* Buffer for keyword string */ - int i; /* Loop index */ - int iv; /* Vectorised keyword index */ - int j; /* Loop index */ - int k; /* Loop index */ - int nin; /* No. of input coords */ - int nout; /* No. of output coords */ - int undef; /* Is the transformation undefined? */ - -/* 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 PolyMap. In this case the - PolyMap belongs to this class, so supply appropriate values to be - passed to the parent class loader (and its parent, etc.). */ - if ( !vtab ) { - size = sizeof( AstPolyMap ); - vtab = &class_vtab; - name = "PolyMap"; - -/* If required, initialise the virtual function table for this class. */ - if ( !class_init ) { - astInitPolyMapVtab( 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 PolyMap. */ - new = astLoadMapping( mem, size, (AstMappingVtab *) vtab, name, - channel ); - - if ( astOK ) { - -/* Get the number of inputs and outputs for the uninverted Mapping. */ - nin = ( (AstMapping *) new )->nin; - nout = ( (AstMapping *) new )->nout; - -/* Read input data. */ -/* ================ */ -/* Request the input Channel to read all the input data appropriate to - this class into the internal "values list". */ - astReadClassData( channel, "PolyMap" ); - -/* Allocate memory to hold the forward arrays. */ - new->ncoeff_f = astMalloc( sizeof( int )*(size_t) nout ); - new->mxpow_f = astMalloc( sizeof( int )*(size_t) nin ); - new->power_f = astMalloc( sizeof( int ** )*(size_t) nout ); - new->coeff_f = astMalloc( sizeof( double * )*(size_t) nout ); - if( astOK ) { - -/* Assume the forward transformation is defined. */ - undef = 0; - -/* Get the maximum power of each input axis value used by the forward - transformation. Set a flag "undef" if no values relating to the - forward transformation are found (this indicates that the forward - transformation is not defined). */ - for( i = 0; i < nin && !undef; i++ ){ - (void) sprintf( buff, "mpf%d", i + 1 ); - (new->mxpow_f)[ i ] = astReadInt( channel, buff, INT_MAX ); - if( (new->mxpow_f)[ i ] == INT_MAX ) undef = 1; - } - -/* Get the number of coefficients associated with each output of the forward - transformation. */ - for( i = 0; i < nout && !undef; i++ ){ - (void) sprintf( buff, "ncf%d", i + 1 ); - (new->ncoeff_f)[ i ] = astReadInt( channel, buff, INT_MAX ); - if( (new->ncoeff_f)[ i ] == INT_MAX ) undef = 1; - } - -/* Get the coefficient values used by the forward transformation. This - uses new style vectorised key names if available. Otherwise it uses - old style indexed names (which were superceded by vectorised names - because they are shorter and so work better with FitsChans). */ - iv = 0; - for( i = 0; i < nout && !undef; i++ ){ - - (new->coeff_f)[ i ] = astMalloc( sizeof( double )* - (size_t) new->ncoeff_f[ i ] ); - if( astOK ) { - for( j = 0; j < new->ncoeff_f[ i ]; j++ ){ - (void) sprintf( buff, "cf%d", ++iv ); - (new->coeff_f)[ i ][ j ] = astReadDouble( channel, buff, AST__BAD ); - if( (new->coeff_f)[ i ][ j ] == AST__BAD ) { - (void) sprintf( buff, "cf%d_%d", i + 1, j + 1 ); - (new->coeff_f)[ i ][ j ] = astReadDouble( channel, buff, AST__BAD ); - } - } - } - } - -/* Get the input axis powers associated with each coefficient of the forward - transformation. */ - iv = 0; - for( i = 0; i < nout && !undef; i++ ){ - (new->power_f)[ i ] = astMalloc( sizeof( int * )* - (size_t) new->ncoeff_f[ i ] ); - if( astOK ) { - for( j = 0; j < new->ncoeff_f[ i ]; j++ ){ - (new->power_f)[ i ][ j ] = astMalloc( sizeof( int )* (size_t) nin ); - if( astOK ) { - for( k = 0; k < nin; k++ ){ - (void) sprintf( buff, "pf%d", ++iv ); - (new->power_f)[ i ][ j ][ k ] = astReadInt( channel, buff, 0 ); - if( (new->power_f)[ i ][ j ][ k ] == 0 ) { - (void) sprintf( buff, "pf%d_%d_%d", i + 1, j + 1, k + 1 ); - (new->power_f)[ i ][ j ][ k ] = astReadInt( channel, buff, 0 ); - } - } - } - } - } - } - -/* Free the arrays if the forward transformation is undefined. */ - if( undef ) { - new->ncoeff_f = astFree( new->ncoeff_f ); - new->mxpow_f = astFree( new->mxpow_f ); - new->power_f = astFree( new->power_f ); - new->coeff_f = astFree( new->coeff_f ); - } - } - -/* Allocate memory to hold the inverse arrays. */ - new->ncoeff_i = astMalloc( sizeof( int )*(size_t) nin ); - new->mxpow_i = astMalloc( sizeof( int )*(size_t) nout ); - new->power_i = astMalloc( sizeof( int ** )*(size_t) nin ); - new->coeff_i = astMalloc( sizeof( double * )*(size_t) nin ); - if( astOK ) { - -/* Assume the inverse transformation is defined. */ - undef = 0; - -/* Get the maximum power of each output axis value used by the inverse - transformation. Set a flag "undef" if no values relating to the - inverse transformation are found (this indicates that the inverse - transformation is not defined). */ - for( i = 0; i < nout && !undef; i++ ){ - (void) sprintf( buff, "mpi%d", i + 1 ); - (new->mxpow_i)[ i ] = astReadInt( channel, buff, INT_MAX ); - if( (new->mxpow_i)[ i ] == INT_MAX ) undef = 1; - } - -/* Get the number of coefficients associated with each input of the inverse - transformation. */ - for( i = 0; i < nin && !undef; i++ ){ - (void) sprintf( buff, "nci%d", i + 1 ); - (new->ncoeff_i)[ i ] = astReadInt( channel, buff, INT_MAX ); - if( (new->ncoeff_i)[ i ] == INT_MAX ) undef = 1; - } - -/* Get the coefficient values used by the inverse transformation. */ - iv = 0; - for( i = 0; i < nin && !undef; i++ ){ - - (new->coeff_i)[ i ] = astMalloc( sizeof( double )* - (size_t) new->ncoeff_i[ i ] ); - if( astOK ) { - for( j = 0; j < new->ncoeff_i[ i ]; j++ ){ - (void) sprintf( buff, "ci%d", ++iv ); - (new->coeff_i)[ i ][ j ] = astReadDouble( channel, buff, AST__BAD ); - if( (new->coeff_i)[ i ][ j ] == AST__BAD ) { - (void) sprintf( buff, "ci%d_%d", i + 1, j + 1 ); - (new->coeff_i)[ i ][ j ] = astReadDouble( channel, buff, AST__BAD ); - } - } - } - } - -/* Get the output axis powers associated with each coefficient of the inverse - transformation. */ - iv = 0; - for( i = 0; i < nin && !undef; i++ ){ - (new->power_i)[ i ] = astMalloc( sizeof( int * )* - (size_t) new->ncoeff_i[ i ] ); - if( astOK ) { - for( j = 0; j < new->ncoeff_i[ i ]; j++ ){ - (new->power_i)[ i ][ j ] = astMalloc( sizeof( int )* (size_t) nout ); - if( astOK ) { - for( k = 0; k < nout; k++ ){ - (void) sprintf( buff, "pi%d", ++iv ); - (new->power_i)[ i ][ j ][ k ] = astReadInt( channel, buff, 0 ); - if( (new->power_i)[ i ][ j ][ k ] == 0 ) { - (void) sprintf( buff, "pi%d_%d_%d", i + 1, j + 1, k + 1 ); - (new->power_i)[ i ][ j ][ k ] = astReadInt( channel, buff, 0 ); - } - } - } - } - } - } - -/* Free the arrays if the inverse transformation is undefined. */ - if( undef ) { - new->ncoeff_i = astFree( new->ncoeff_i ); - new->mxpow_i = astFree( new->mxpow_i ); - new->power_i = astFree( new->power_i ); - new->coeff_i = astFree( new->coeff_i ); - } - } - -/* Whether to use an iterative inverse transformation. */ - new->iterinverse = astReadInt( channel, "iterinv", -INT_MAX ); - if ( TestIterInverse( new, status ) ) SetIterInverse( new, new->iterinverse, status ); - -/* Max number of iterations for iterative inverse transformation. */ - new->niterinverse = astReadInt( channel, "niterinv", -INT_MAX ); - if ( TestNiterInverse( new, status ) ) SetNiterInverse( new, new->niterinverse, status ); - -/* Target relative error for iterative inverse transformation. */ - new->tolinverse = astReadDouble( channel, "tolinv", AST__BAD ); - if ( TestTolInverse( new, status ) ) SetTolInverse( new, new->tolinverse, status ); - -/* The Jacobian of the PolyMap's forward transformation has not yet been - found. */ - new->jacobian = NULL; - -/* If an error occurred, clean up by deleting the new PolyMap. */ - if ( !astOK ) new = astDelete( new ); - } - -/* Return the new PolyMap 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. */ - -void astPolyPowers_( AstPolyMap *this, double **work, int ncoord, - const int *mxpow, double **ptr, int point, int fwd, - int *status ){ - if ( !astOK ) return; - (**astMEMBER(this,PolyMap,PolyPowers))( this, work, ncoord, mxpow, ptr, - point, fwd, status ); -} - -AstPolyMap *astPolyTran_( AstPolyMap *this, int forward, double acc, - double maxacc, int maxorder, const double *lbnd, - const double *ubnd, int *status ){ - if ( !astOK ) return NULL; - return (**astMEMBER(this,PolyMap,PolyTran))( this, forward, acc, - maxacc, maxorder, lbnd, - ubnd, status ); -} - -void astPolyCoeffs_( AstPolyMap *this, int forward, int nel, double *array, - int *ncoeff, int *status ){ - if ( !astOK ) return; - (**astMEMBER(this,PolyMap,PolyCoeffs))( this, forward, nel, - array, ncoeff, status ); -} - -void astFitPoly1DInit_( AstPolyMap *this, int forward, double **table, - AstMinPackData *data, double *scales, - int *status ){ - if ( !astOK ) return; - (**astMEMBER(this,PolyMap,FitPoly1DInit))( this, forward, table, data, scales, - status ); -} - - -void astFitPoly2DInit_( AstPolyMap *this, int forward, double **table, - AstMinPackData *data, double *scales, - int *status ){ - if ( !astOK ) return; - (**astMEMBER(this,PolyMap,FitPoly2DInit))( this, forward, table, data, scales, - status ); -} - - - - |