/* $Id: cpl_wcs.c,v 1.17 2008/02/26 10:36:43 scastro Exp $ * * This file is part of the ESO Common Pipeline Library * Copyright (C) 2001-2004 European Southern Observatory * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 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 General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ /* * $Author: scastro $ * $Date: 2008/02/26 10:36:43 $ * $Revision: 1.17 $ * $Name: $ */ #ifdef HAVE_CONFIG_H #include #endif /*---------------------------------------------------------------------------- Includes ----------------------------------------------------------------------------*/ #include #include #include #include #include #include #include #include #include #include "cpl_wcs.h" #ifdef CPL_WCS_INSTALLED /* If WCS is installed */ #include #endif /* End If WCS is installed */ /* Remove this when cfitsio is upgraded */ #include /*---------------------------------------------------------------------------*/ /** * @defgroup cpl_wcs World Coordinate System * * This module provides functions to manipulate FITS World Coordinate Systems * * A @em cpl_wcs is an object containing a pointer to the WCSLIB structure * and the physical dimensions of the image from which the WCS was read. * The functionality provided includes general transformations between physical * and world coordinates as well as a few convenience routines for * x,y <=> RA,Dec transformations. * * @par Synopsis: * @code * #include "cpl_wcs.h" * @endcode */ /*---------------------------------------------------------------------------*/ /**@{*/ /*---------------------------------------------------------------------------- Type definition ----------------------------------------------------------------------------*/ #ifdef CPL_WCS_INSTALLED /* If WCS is installed */ struct _cpl_wcs_ { struct wcsprm *wcsptr; /* WCSLIB structure */ int naxis; /* Number of dimensions of the image */ int *dims; /* Dimensions of image */ }; /* Static routine declarations */ static cpl_wcs *cpl_wcs_init(void); /* Fix when cfitsio is upgraded */ /* static char *cpl_wcs_plist2fitsstr(const cpl_propertylist *self); */ static char *cpl_wcs_plist2fitsstr(const cpl_propertylist *self, int *nkeys); static cpl_propertylist *cpl_wcs_fitsstr2plist(char *fitsstr); static void cpl_wcs_platesol_4(cpl_matrix *xy, cpl_matrix *std, cpl_array *bad, cpl_array **plateconsts); static void cpl_wcs_platesol_6(cpl_matrix *xy, cpl_matrix *std, cpl_array *bad, cpl_array **plateconsts); /* Remove when cfitsio is upgraded */ static int cpl_wcs_ffhdr2str(fitsfile *fptr, int exclude_comm, char **exclist, int nexc, char **header, int *nkeys, int *status); /* Useful conversion factors */ /* Using CPL_MATH_2PI instead */ /* #define M_TWOPI 6.28218530717958647692 */ /* FITS END card */ static const char *endcard = "END "; /* WCSLIB error messages */ #define WCSLIB_ERRCODE_MAX 9 static const char *wcslib_errmsgs[WCSLIB_ERRCODE_MAX+1] = { "", "WCSLIB undefined input structure pointer", "WCSLIB unable to allocate required memory", "WCSLIB linear transformation matrix is singular", "WCSLIB invalid coordinate axis types", "WCSLIB invalid parameter value", "WCSLIB invalid coordinate transformation parameters", "WCSLIB Ill-conditioned coordinate transformation parameters", "WCSLIB One or more input coordinates invalid", "WCSLIB One or more input coordinates invalid"}; #endif /* End If WCS is installed */ /*-------------------------------------------------------------------------*/ /*-------------------------------------------------------------------------*/ /** * @brief * Create a wcs structure by parsing a propertylist. * * @param plist The input propertylist * * @return * The newly created and filled cpl_wcs object or NULL if it could not be * created. In the latter case an appropriate error code is set. * * @error * * * * * * * * * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The parameter plist is a NULL pointer. *
CPL_ERROR_TYPE_MISMATCH * NAXIS information in image propertylist is not an integer *
CPL_ERROR_DATA_NOT_FOUND * Error in getting NAXIS information for image propertylists *
CPL_ERROR_NO_WCS * The WCS sub library is not available. *
* @enderror * * The function allocates memory for a WCS structure. A pointer to the WCSLIB * header information is created by parsing the FITS WCS keywords from the * header of a file. A few ancillary items are also filled in. * * The returned property must be destroyed using the wcs destructor * @b cpl_wcs_delete(). * * @see cpl_wcs_delete() */ cpl_wcs *cpl_wcs_new_from_propertylist(const cpl_propertylist *plist) { #ifdef CPL_WCS_INSTALLED /* If WCS is installed */ const char *_id = "cpl_wcs_new"; char *shdr,nax[9]; cpl_wcs *wcs; int retval,nrej,nwcs,np,i; struct wcsprm *wwcs = NULL; /* Check to see if the propertylist exists */ if (plist == NULL) { cpl_error_set(_id,CPL_ERROR_NULL_INPUT); return(NULL); } /* Get the size of the propertylist */ /* np = cpl_propertylist_get_size(plist); */ /* Get a cpl_wcs structure and initialise it */ if ((wcs = cpl_wcs_init()) == NULL) return(NULL); /* Convert the propertylist into a string */ /* shdr = cpl_wcs_plist2fitsstr(plist); */ shdr = cpl_wcs_plist2fitsstr(plist,&np); if (! shdr) { cpl_error_set_where(_id); cpl_wcs_delete(wcs); return(NULL); } /* Parse the header to get the wcslib structure */ retval = wcspih(shdr,np,0,0,&nrej,&nwcs,&wwcs); if (retval != 0) { cpl_msg_error(_id,"Cannot parse WCS header"); cpl_free(shdr); wcsvfree(&nwcs,&wwcs); cpl_wcs_delete(wcs); return(NULL); } free(shdr); /* Now extract the one you want and ditch the rest */ wcs->wcsptr = (struct wcsprm *)cpl_calloc(1,sizeof(struct wcsprm)); (wcs->wcsptr)->flag = -1; wcscopy(1,wwcs,wcs->wcsptr); wcsset(wcs->wcsptr); wcsvfree(&nwcs,&wwcs); /* Fill the rest of the stuff */ if (cpl_propertylist_has(plist,"NAXIS")) { wcs->naxis = cpl_propertylist_get_int(plist,"NAXIS"); if (cpl_error_get_code() != CPL_ERROR_NONE) { cpl_wcs_delete(wcs); return(NULL); } wcs->dims = cpl_calloc(wcs->naxis,sizeof(int)); for (i = 1; i <= wcs->naxis; i++) { (void)snprintf(nax,9,"NAXIS%d",i); wcs->dims[i-1] = cpl_propertylist_get_int(plist,nax); if (cpl_error_get_code() != CPL_ERROR_NONE) { cpl_wcs_delete(wcs); return(NULL); } } } /* Get out of here */ return(wcs); #else cpl_ensure(0, CPL_ERROR_NO_WCS, NULL); #endif /* End If WCS is installed */ } /** * @brief * Destroy a WCS structure * * @param wcs The WCS structure to destroy * * @return * Nothing. * * @error * * * * * *
CPL_ERROR_NO_WCS * The WCS sub library is not available. *
* @enderror * The function destroys the WCS structure @em wcs and its whole * contents. */ void cpl_wcs_delete(cpl_wcs *wcs) { #ifdef CPL_WCS_INSTALLED /* If WCS is installed */ /* Free the workspace */ if (wcs == NULL) return; if (wcs->wcsptr != NULL) { (void)wcsfree(wcs->wcsptr); cpl_free(wcs->wcsptr); } if (wcs->dims != NULL) cpl_free(wcs->dims); cpl_free(wcs); #else cpl_error_set(cpl_func, CPL_ERROR_NO_WCS); #endif /* End If WCS is installed */ } /** * @brief * Convert between physical and world coordiantes. * * @param wcs The input cpl_wcs structure * @param from The input coordinate matrix * @param to The output coordinate matrix * @param status The output status array * @param transform The transformation mode * * @return * An appropriate error code. * * @error * * * * * * * * * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The parameter wcs is a NULL pointer, the parameter * from is a NULL pointer or wcs is missing * some of its information. *
CPL_ERROR_UNSPECIFIED * No rows or columns in the input matrix, or an unspecified error * has occurred in the WCSLIB routine *
CPL_ERROR_UNSUPPORTED_MODE * The input conversion mode is not supported *
CPL_ERROR_NO_WCS * The WCS sub library is not available. *
* @enderror * * This function converts between several types of coordinates. These include: * -- physical coordinates: The physical location on a detector (i.e. * pixel coordinates) * -- world coordinates: The real astronomical coordinate system for the * observations. This may be spectral, celestial, * time, etc. * -- standard coordinates: These are an intermediate relative coordinate * representation, defined as a distance from * a reference point in the natural units of the * world coordinate system. Any defined projection * geometry will have already been included in * the definition of standard coordinates. * * The supported conversion modes are: * -- CPL_WCS_PHYS2WORLD: Converts input physical to world coordinates * -- CPL_WCS_WORLD2PHYS: Converts input world to physical coordinates * -- CPL_WCS_WORLD2STD: Converts input world to standard coordinates * -- CPL_WCS_PHYS2STD: Converts input physical to standard coordinates * The output matrix and status arrays will be allocated here, and thus will * need to be freed by the calling routine. The status array is used to flag * input coordinates where there has been some sort of failure in the * transformation. * */ cpl_error_code cpl_wcs_convert(const cpl_wcs *wcs, const cpl_matrix *from, cpl_matrix **to, cpl_array **status, cpl_wcs_trans_mode transform) { #ifdef CPL_WCS_INSTALLED /* If WCS is installed */ int nrows,ncols,*sdata,retval; cpl_matrix *x1; cpl_array *x2,*x3; cpl_error_code code; double *tdata,*x1data,*x2data,*x3data; const double *fdata; const char *_id = "cpl_wcs_convert"; /* Initialise output */ *to = NULL; *status = NULL; /* Basic checks on the input pointers */ if (wcs == NULL || wcs->wcsptr == NULL || from == NULL) { cpl_error_set(_id,CPL_ERROR_NULL_INPUT); return(CPL_ERROR_NULL_INPUT); } /* Check the number of rows/columns for the input matrix */ nrows = cpl_matrix_get_nrow(from); ncols = cpl_matrix_get_ncol(from); if (nrows == 0 || ncols == 0) { cpl_error_set(_id,CPL_ERROR_UNSPECIFIED); return(CPL_ERROR_UNSPECIFIED); } /* Get the output memory and some memory for wcslib to use */ *to = cpl_matrix_new(nrows,ncols); *status = cpl_array_new(nrows,CPL_TYPE_INT); x1 = cpl_matrix_new(nrows,ncols); x2 = cpl_array_new(nrows,CPL_TYPE_DOUBLE); x3 = cpl_array_new(nrows,CPL_TYPE_DOUBLE); /* Now get the pointers for the data arrays */ fdata = cpl_matrix_get_data_const(from); tdata = cpl_matrix_get_data(*to); sdata = cpl_array_get_data_int(*status); x1data = cpl_matrix_get_data(x1); x2data = cpl_array_get_data_double(x2); x3data = cpl_array_get_data_double(x3); /* Right, now switch for the transform type. First physical to world coordinates */ switch (transform) { case CPL_WCS_PHYS2WORLD: retval = wcsp2s(wcs->wcsptr,nrows,wcs->naxis,fdata,x1data,x2data, x3data,tdata,sdata); break; case CPL_WCS_WORLD2PHYS: retval = wcss2p(wcs->wcsptr,nrows,wcs->naxis,fdata,x2data,x3data, x1data,tdata,sdata); break; case CPL_WCS_WORLD2STD: retval = wcss2p(wcs->wcsptr,nrows,wcs->naxis,fdata,x2data,x3data, tdata,x1data,sdata); break; case CPL_WCS_PHYS2STD: retval = wcsp2s(wcs->wcsptr,nrows,wcs->naxis,fdata,tdata,x2data, x3data,x1data,sdata); break; default: cpl_error_set(_id,CPL_ERROR_UNSUPPORTED_MODE); cpl_matrix_delete(*to); *to = NULL; cpl_array_delete(*status); *status = NULL; cpl_matrix_delete(x1); return(CPL_ERROR_UNSUPPORTED_MODE); } /* Ditch the intermediate coordinate results */ cpl_matrix_delete(x1); cpl_array_delete(x2); cpl_array_delete(x3); /* Now interpret the return value from wcslib */ switch (retval) { case 0: code = CPL_ERROR_NONE; break; case 1: code = CPL_ERROR_NULL_INPUT; cpl_msg_error(_id,wcslib_errmsgs[1]); cpl_error_set(_id,code); break; case 2: case 3: case 4: case 5: case 6: case 7: case 8: case 9: code = CPL_ERROR_UNSPECIFIED; cpl_error_set(_id,code); cpl_msg_error(_id,wcslib_errmsgs[retval]); break; default: code = CPL_ERROR_UNSPECIFIED; cpl_error_set(_id,code); cpl_msg_error(_id,"WCSLIB found an unspecified error: %d", retval); break; } return(code); #else cpl_ensure_code(0, CPL_ERROR_NO_WCS); #endif /* End If WCS is installed */ } /** * @brief * Do a 2d plate solution given physical and celestial coordinates * * @param ilist The input property list containing the first pass WCS * @param cel The celestial coordinate matrix * @param xy The physical coordinate matrix * @param niter The number of fitting iterations * @param thresh The threshold for the fitting rejection cycle * @param fitmode The fitting mode (see below) * @param outmode The output mode (see below) * @param olist The output propertylist containing the new WCS * * @return * An appropriate error code. * * @error * * * * * * * * * * * * * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The parameter cel is a NULL pointer, the parameter * xy is a NULL pointer or ilist is a * NULL pointer. *
CPL_ERROR_UNSPECIFIED * Unable to parse the input propertylist into a proper FITS WCS or * there are too few points in the input matrices for a fit. *
CPL_ERROR_INCOMPATIBLE_INPUT * The matrices cel and xy have different sizes. *
CPL_ERROR_UNSUPPORTED_MODE * Either fitmode or outmode are specified incorrectly *
CPL_ERROR_NO_WCS * The WCS sub library is not available. *
* @enderror * * This function allows for the following type of fits: * -- CPL_WCS_PLATESOL_4: Fit for zero point, 1 scale and 1 rotation. * -- CPL_WCS_PLATESOL_6: Fit for zero point, 2 scales, 1 rotation, 1 shear. * * This fucntion allows the zeropoint to be defined by shifting either the * physical or the celestial coordinates of the reference point: * -- CPL_WCS_MV_CRVAL: Keeps the physical point fixed and shifts the celestial * -- CPL_WCS_MV_CRPIX: Keeps the celestial point fixed and shifts the physical * * The output propertylist contains WCS relevant information only. * */ cpl_error_code cpl_wcs_platesol(cpl_propertylist *ilist, cpl_matrix *cel, cpl_matrix *xy, int niter, float thresh, cpl_wcs_platesol_fitmode fitmode, cpl_wcs_platesol_outmode outmode, cpl_propertylist **olist) { #ifdef CPL_WCS_INSTALLED /* If WCS is installed */ int npts,iter,n,i,nrejnew,*isbad,ngood,ntotrej,retval,stat,ival; double xifit,etafit,*work,mederr,mederr_xi,mederr_eta,*xydata,*stddata; double *pc,crpix1,crpix2,crval1,crval2,phi,theta,dval; const char *_id = "cpl_wcs_platesol"; char *outstr,*o; struct wcsprm *wp; cpl_wcs *wcs; cpl_matrix *std; cpl_array *status,*bad,*plateconsts; cpl_vector *vwork; cpl_property *p; cpl_propertylist *test; /* Initialise the output pointer */ *olist = NULL; /* Basic checks on the input pointers */ if (cel == NULL || xy == NULL || ilist == NULL) { cpl_error_set(_id,CPL_ERROR_NULL_INPUT); return(CPL_ERROR_NULL_INPUT); } /* Open the cpl_wcs structure */ wcs = cpl_wcs_new_from_propertylist(ilist); if (wcs == NULL) { cpl_msg_error(_id,"Unable to parse header"); cpl_error_set(_id,CPL_ERROR_UNSPECIFIED); return(CPL_ERROR_UNSPECIFIED); } /* Get the number of celestial points and compare this with the size of the matrix with the xy coordinates. Also look at the total number of points available */ npts = cpl_matrix_get_nrow(cel); if (npts != cpl_matrix_get_nrow(xy)) { cpl_wcs_delete(wcs); cpl_error_set(_id,CPL_ERROR_INCOMPATIBLE_INPUT); return(CPL_ERROR_INCOMPATIBLE_INPUT); } if (npts < 2) { cpl_wcs_delete(wcs); cpl_msg_error(_id,"Insufficient points for a fit"); cpl_error_set(_id,CPL_ERROR_UNSPECIFIED); return(CPL_ERROR_UNSPECIFIED); } /* Convert the celestial coordinates to standard coordinates */ cpl_wcs_convert(wcs,cel,&std,&status,CPL_WCS_WORLD2STD); cpl_array_delete(status); /* Get an array to flag bad pairs */ bad = cpl_array_new(npts,CPL_TYPE_INT); isbad = cpl_array_get_data_int(bad); for (i = 0; i < npts; i++) isbad[i] = 0; /* Get some workspace for rejection algorithm */ if ((work = cpl_malloc(2*npts*sizeof(double))) == NULL) { cpl_wcs_delete(wcs); cpl_array_delete(bad); return(cpl_error_get_code()); } /* A couple of convenience variables */ xydata = cpl_matrix_get_data(xy); stddata = cpl_matrix_get_data(std); /* Iterative loop */ iter = 1; plateconsts = NULL; while (iter <= niter) { /* Check the number of good points left in this iteration */ ngood = 0; for (i = 0; i < npts; i++) ngood += (isbad[i] == 0); if (ngood < 2) break; if (plateconsts != NULL) cpl_array_delete(plateconsts); /* Do a plate solution */ switch (fitmode) { case CPL_WCS_PLATESOL_6: cpl_wcs_platesol_6(xy,std,bad,&plateconsts); break; case CPL_WCS_PLATESOL_4: cpl_wcs_platesol_4(xy,std,bad,&plateconsts); break; default: cpl_error_set(_id,CPL_ERROR_UNSUPPORTED_MODE); cpl_matrix_delete(std); cpl_wcs_delete(wcs); cpl_array_delete(bad); return(CPL_ERROR_UNSUPPORTED_MODE); } /* Get the fit residuals */ pc = cpl_array_get_data_double(plateconsts); n = 0; for (i = 0; i < npts; i++) { if (isbad[i] == 1) continue; xifit = xydata[2*i]*pc[0] + xydata[2*i+1]*pc[1] + pc[2]; etafit = xydata[2*i]*pc[3] + xydata[2*i+1]*pc[4] + pc[5]; work[n++] = fabs(xifit - stddata[2*i]); work[n++] = fabs(etafit - stddata[2*i+1]); } /* Get the median of the array */ vwork = cpl_vector_wrap(n,work); mederr = 1.48*cpl_vector_get_median(vwork); cpl_vector_unwrap(vwork); /* Get out of here if it's the last iteration */ if (iter == niter) break; /* Now reject the bad ones... */ ntotrej = 0; nrejnew = 0; for (i = 0; i < npts; i++) { if (isbad[i] == 1) { ntotrej++; continue; } if (work[i*2] > thresh*mederr || work[i*2+1] > thresh*mederr) { isbad[i] = 1; nrejnew++; ntotrej++; } } /* Were any more stars rejected in this iteration */ if (nrejnew == 0) break; iter++; } /* Now work out the median error in each axis */ n = 0; for (i = 0; i < npts; i++) { if (isbad[i] == 1) continue; xifit = xydata[2*i]*pc[0] + xydata[2*i+1]*pc[1] + pc[2]; work[n++] = fabs(xifit - stddata[2*i]); } vwork = cpl_vector_wrap(n,work); mederr_xi = 1.48*cpl_vector_get_median(vwork); cpl_vector_unwrap(vwork); n = 0; for (i = 0; i < npts; i++) { if (isbad[i] == 1) continue; etafit = xydata[2*i]*pc[3] + xydata[2*i+1]*pc[4] + pc[5]; work[n++] = fabs(etafit - stddata[2*i+1]); } vwork = cpl_vector_wrap(n,work); mederr_eta = 1.48*cpl_vector_get_median(vwork); /* Do some intermediate tidying */ cpl_vector_unwrap(vwork); cpl_free(work); cpl_array_delete(bad); cpl_matrix_delete(std); /* Define the reference point result */ wp = wcs->wcsptr; switch (outmode) { case CPL_WCS_MV_CRPIX: crpix1 = (pc[4]*pc[2] - pc[1]*pc[5])/(pc[3]*pc[1] - pc[4]*pc[0]); crpix2 = (pc[0]*pc[5] - pc[3]*pc[2])/(pc[3]*pc[1] - pc[4]*pc[0]); crval1 = (wp->crval)[0]; crval2 = (wp->crval)[1]; break; case CPL_WCS_MV_CRVAL: crpix1 = (wp->crpix)[0]; crpix2 = (wp->crpix)[1]; xifit = crpix1*pc[0] + crpix2*pc[1] + pc[2]; etafit = crpix1*pc[3] + crpix2*pc[4] + pc[5]; retval = celx2s(&(wp->cel),1,1,2,2,&xifit,&etafit,&phi,&theta, &crval1,&crval2,&stat); break; } /* Now update the wcs structure */ (wp->crval)[0] = crval1; (wp->crval)[1] = crval2; (wp->crpix)[0] = crpix1; (wp->crpix)[1] = crpix2; (wp->pc)[0] = pc[0]; (wp->pc)[1] = pc[1]; (wp->pc)[2] = pc[3]; (wp->pc)[3] = pc[4]; for (i = 0; i < 4; i++) (wp->cd)[i] = (wp->pc)[i]; (wp->cdelt)[0] = 1.0; (wp->cdelt)[1] = 1.0; (wp->csyer)[0] = mederr_xi; (wp->csyer)[1] = mederr_eta; /* Make a FITS string and convert it to a propertylist */ (void)wcshdo(0,wp,&i,&o); outstr = cpl_calloc((i+1)*80+1,1); strncpy(outstr,o,i*80); strncat(outstr,endcard,80); *olist = cpl_wcs_fitsstr2plist(outstr); free(o); cpl_free(outstr); /* Ok, massage the result a bit to make it look a bit nicer */ p = cpl_propertylist_get_property(*olist,"PC1_1"); cpl_property_set_name(p,"CD1_1"); p = cpl_propertylist_get_property(*olist,"PC1_2"); cpl_property_set_name(p,"CD1_2"); p = cpl_propertylist_get_property(*olist,"PC2_1"); cpl_property_set_name(p,"CD2_1"); p = cpl_propertylist_get_property(*olist,"PC2_2"); cpl_property_set_name(p,"CD2_2"); cpl_propertylist_erase(*olist,"CDELT1"); cpl_propertylist_erase(*olist,"CDELT2"); cpl_propertylist_erase(*olist,"RESTFRQ"); cpl_propertylist_erase(*olist,"RESTWAV"); cpl_propertylist_erase(*olist,"END"); /* Fix boo-boo in wcslib */ test = cpl_propertylist_new(); cpl_propertylist_copy_property_regexp(test,(const cpl_propertylist *)*olist, "PV",0); n = cpl_propertylist_get_size(test); for (i = 0; i < n; i++) { p = cpl_propertylist_get(test,i); if (cpl_property_get_type(p) == CPL_TYPE_INT) { ival = cpl_property_get_int(p); dval = (double)ival; cpl_propertylist_erase(*olist,cpl_property_get_name(p)); cpl_propertylist_append_double(*olist,cpl_property_get_name(p),dval); } } cpl_propertylist_delete(test); /* Tidy and exit */ cpl_array_delete(plateconsts); cpl_wcs_delete(wcs); return CPL_ERROR_NONE; #else cpl_ensure_code(0, CPL_ERROR_NO_WCS); #endif /* End If WCS is installed */ } /**@}*/ #ifdef CPL_WCS_INSTALLED /* If WCS is installed */ /*-------------------------------------------------------------------------*/ /** * @brief * Create an empty wcs structure. * * @return * The output wcs structure * * @error * None * * This is a static routine that creates an empty cpl_wcs structure. */ static cpl_wcs *cpl_wcs_init(void) { cpl_wcs *wcs = NULL; /* Get the main structure workspace */ wcs = cpl_malloc(sizeof(cpl_wcs)); if (wcs == NULL) return(NULL); /* Initialise the output structure */ wcs->wcsptr = NULL; wcs->naxis = 0; wcs->dims = NULL; return(wcs); } /** * @brief * Convert a propertylist to a FITS string * * @param self The input propertylist * * @return * The output character string with the properties formatted as in a FITS * header. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The parameter self is a NULL pointer. *
* @enderror * * This converts a propertylist into a single string with all properties formatted * as FITS cards. This is needed for wcspih. The output string must be freed by * the calling routine. */ /* static char *cpl_wcs_plist2fitsstr(const cpl_propertylist *self) { */ static char *cpl_wcs_plist2fitsstr(const cpl_propertylist *self, int *nkeys) { const char *_id = "cpl_wcs_plist2fitsstr"; char *header; int n,status; fitsfile *fptr; /* Sanity testing of input propertylist */ if (self == NULL) { cpl_error_set(_id,CPL_ERROR_NULL_INPUT); return(NULL); } /* Find out how big the propertylist is */ n = cpl_propertylist_get_size(self); /* Open a memory file with CFITSIO */ status = 0; (void)fits_create_file(&fptr,"mem://",&status); /* (void)fits_create_img(fptr,BYTE_IMG,0,NULL,&status); */ /* Add the properties into the FITS file */ cpl_propertylist_to_fitsfile(fptr,self,NULL); /* Parse the header */ /* (void)fits_hdr2str(fptr,1,NULL,0,&header,&ival,&status); */ (void)cpl_wcs_ffhdr2str(fptr,1,NULL,0,&header,nkeys,&status); (void)fits_close_file(fptr,&status); /* Get out of here */ return(header); } /** * @brief * Convert a FITS string to a propertylist * * @param fitsstr The input FITS header string * * @return * The output propertylist. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The parameter fitsstr is a NULL pointer. *
* @enderror * * This converts a single string formatted with FITS cards into a propertylist. * This is needed for wcspih. The output propertylist must be freed by the * calling routine */ static cpl_propertylist *cpl_wcs_fitsstr2plist(char *fitsstr) { int status; fitsfile *fptr; char *f; const char *_id="cpl_wcs_fitsstr2plist"; cpl_propertylist *p; /* Check input */ if (fitsstr == NULL) { cpl_error_set(_id,CPL_ERROR_NULL_INPUT); return(NULL); } /* Create a memory file using CFITSIO. We can then add individual cards to the header and allow CFITSIO to decide on data types etc. */ status = 0; (void)fits_create_file(&fptr,"mem://",&status); /* (void)fits_create_img(fptr,BYTE_IMG,0,NULL,&status); */ /* Loop for all the cards in the FITS string and add them to the memory FITS header */ f = fitsstr; while (strncmp(f,"END ",8)) { (void)fits_insert_card(fptr,f,&status); f += (FLEN_CARD - 1); } (void)fits_insert_card(fptr,f,&status); /* Now create the propertylist */ p = cpl_propertylist_from_fitsfile(fptr); /* Close the memory file and get out of here */ (void)fits_close_file(fptr,&status); return(p); } /** * @brief * Do a 6 plate constant fit * * @param xy The matrix of physical coordinates * @param std The matrix of celestial coordinates * @param bad An array to flag points to be ignored from the fit * @param plateconsts The output array of plate constants * * @return * Nothing * * This routine fits the constants a,b,c,d,e,f to the equations: * xi = ax + by + c * eta = dx + ey + f * The values of these coefficients are passed back in the plateconsts * array. * */ static void cpl_wcs_platesol_6(cpl_matrix *xy, cpl_matrix *std, cpl_array *bad, cpl_array **plateconsts) { double sx1sq,sy1sq,sx1y1,sx1x2,sy1x2,*xydata,*stddata,*pc; double sy1y2,sx1y2,xposmean,yposmean,ximean,etamean,xx1,yy1,xx2,yy2; int i,ngood,*isbad,nstds; /* Get some convenience variables */ xydata = cpl_matrix_get_data(xy); stddata = cpl_matrix_get_data(std); isbad = cpl_array_get_data_int(bad); nstds = cpl_array_get_size(bad); /* Initialise all the counters and summations */ sx1sq = 0.0; sy1sq = 0.0; sx1y1 = 0.0; sx1x2 = 0.0; sy1x2 = 0.0; sy1y2 = 0.0; sx1y2 = 0.0; xposmean = 0.0; yposmean = 0.0; ximean = 0.0; etamean = 0.0; /* Find means in each coordinate system */ ngood = 0; for (i = 0; i < nstds; i++) { if (isbad[i] != 0) continue; xposmean += xydata[2*i]; yposmean += xydata[2*i+1]; ximean += stddata[2*i]; etamean += stddata[2*i+1]; ngood++; } xposmean /= (double)ngood; yposmean /= (double)ngood; ximean /= (double)ngood; etamean /= (double)ngood; /* Now accumulate the sums */ for (i = 0; i < nstds; i++) { if (! isbad[i]) { xx1 = xydata[2*i] - xposmean; yy1 = xydata[2*i+1] - yposmean; xx2 = stddata[2*i] - ximean; yy2 = stddata[2*i+1] - etamean; sx1sq += xx1*xx1; sy1sq += yy1*yy1; sx1y1 += xx1*yy1; sx1x2 += xx1*xx2; sy1x2 += yy1*xx2; sy1y2 += yy1*yy2; sx1y2 += xx1*yy2; } } /* Get an output array for the results */ *plateconsts = cpl_array_new(6,CPL_TYPE_DOUBLE); pc = cpl_array_get_data_double(*plateconsts); /* Do solution for X */ pc[0] = (sx1y1*sy1x2 - sx1x2*sy1sq)/(sx1y1*sx1y1 - sx1sq*sy1sq); pc[1] = (sx1x2*sx1y1 - sx1sq*sy1x2)/(sx1y1*sx1y1 - sx1sq*sy1sq); pc[2] = -xposmean*pc[0] - yposmean*pc[1] + ximean; /* Now the solution for Y */ pc[3] = (sy1y2*sx1y1 - sy1sq*sx1y2)/(sx1y1*sx1y1 - sy1sq*sx1sq); pc[4] = (sx1y1*sx1y2 - sy1y2*sx1sq)/(sx1y1*sx1y1 - sy1sq*sx1sq); pc[5] = -xposmean*pc[3] - yposmean*pc[4] + etamean; } /** * @brief * Do a 4 plate constant fit * * @param xy The matrix of physical coordinates * @param std The matrix of celestial coordinates * @param bad An array to flag points to be ignored from the fit * @param plateconsts The output array of plate constants * * @return * Nothing * * This routine fits the constants a,b,c,d,e,f to the equations: * xi = ax + by + c * eta = dx + ey + f * but where the scale and rotation implied by the coefficients a,b,d,e are * constrained to be the same for each axis. The 6 coefficients are passed * back in the plateconsts array. * */ static void cpl_wcs_platesol_4(cpl_matrix *xy, cpl_matrix *std, cpl_array *bad, cpl_array **plateconsts) { double sx1sq,sy1sq,sx1x2,sy1x2,sy1y2,sx1y2,xposmean,yposmean,*pc; double ximean,etamean,xx1,yy1,xx2,yy2,det,num,denom,theta,mag; double stheta,ctheta,*xydata,*stddata; int i,ngood,*isbad,nstds; /* Get some convenience variables */ xydata = cpl_matrix_get_data(xy); stddata = cpl_matrix_get_data(std); isbad = cpl_array_get_data_int(bad); nstds = cpl_array_get_size(bad); /* Initialise all the counters and summations */ sx1sq = 0.0; sy1sq = 0.0; sx1x2 = 0.0; sy1x2 = 0.0; sy1y2 = 0.0; sx1y2 = 0.0; xposmean = 0.0; yposmean = 0.0; ximean = 0.0; etamean = 0.0; /* Find means in each coordinate system */ ngood = 0; for (i = 0; i < nstds; i++) { if (isbad[i] != 0) continue; xposmean += xydata[2*i]; yposmean += xydata[2*i+1]; ximean += stddata[2*i]; etamean += stddata[2*i+1]; ngood++; } xposmean /= (double)ngood; yposmean /= (double)ngood; ximean /= (double)ngood; etamean /= (double)ngood; /* Now accumulate the sums */ for (i = 0; i < nstds; i++) { if (! isbad[i]) { xx1 = xydata[2*i] - xposmean; yy1 = xydata[2*i+1] - yposmean; xx2 = stddata[2*i] - ximean; yy2 = stddata[2*i+1] - etamean; sx1sq += xx1*xx1; sy1sq += yy1*yy1; sx1x2 += xx1*xx2; sy1x2 += yy1*xx2; sy1y2 += yy1*yy2; sx1y2 += xx1*yy2; } } /* Compute the rotation angle */ det = sx1x2*sy1y2 - sy1x2*sx1y2; if (det < 0.0) { num = sy1x2 + sx1y2; denom = -sx1x2 + sy1y2; } else { num = sy1x2 - sx1y2; denom = sx1x2 + sy1y2; } if (num == 0.0 && denom == 0.0) { theta = 0.0; } else { theta = atan2(num,denom); if (theta < 0.0) theta += CPL_MATH_2PI; } /* Compute magnification factor */ ctheta = cos(theta); stheta = sin(theta); num = denom*ctheta + num*stheta; denom = sx1sq + sy1sq; if (denom <= 0.0) { mag = 1.0; } else { mag = num/denom; } /* Get an output array for the results */ *plateconsts = cpl_array_new(6,CPL_TYPE_DOUBLE); pc = cpl_array_get_data_double(*plateconsts); /* Compute coeffs */ if (det < 0.0) { pc[0] = -mag*ctheta; pc[3] = mag*stheta; } else { pc[0] = mag*ctheta; pc[3] = -mag*stheta; } pc[1] = mag*stheta; pc[4] = mag*ctheta; pc[2] = -xposmean*pc[0] - yposmean*pc[1] + ximean; pc[5] = -xposmean*pc[3] - yposmean*pc[4] + etamean; } /*-------------------------------------------------------------------------*/ /* * NOTE: * Work around for cfitsio-2.5.10 memory errors in ffhdr2str(). This * problem has been fixed in cfitsio > 3.0.3 */ /* read header keywords into a long string of chars. This routine allocates memory for the string, so the calling routine must eventually free the memory when it is not needed any more. If exclude_comm is TRUE, then all the COMMENT, HISTORY, and keywords will be excluded from the output string of keywords. Any other list of keywords to be excluded may be specified with the exclist parameter. */ static int cpl_wcs_ffhdr2str( fitsfile *fptr, /* I - FITS file pointer */ int exclude_comm, /* I - if TRUE, exclude commentary keywords */ char **exclist, /* I - list of excluded keyword names */ int nexc, /* I - number of names in exclist */ char **header, /* O - returned header string */ int *nkeys, /* O - returned number of 80-char keywords */ int *status) /* IO - error status */ { int casesn, match, exact, totkeys; long ii, jj; char keybuf[162], keyname[FLEN_KEYWORD], *headptr; /* NOTE: Use the actual CFITSIO function if the version number is more recent than 3.0.3 */ float version; fits_get_version(&version); if (version >= 3.030) return ffhdr2str(fptr, exclude_comm, exclist, nexc, header, nkeys, status); *nkeys = 0; if (*status > 0) return(*status); /* get number of keywords in the header (doesn't include END) */ if (ffghsp(fptr, &totkeys, NULL, status) > 0) return(*status); /* allocate memory for all the keywords (multiple of 2880 bytes) */ *header = (char *) calloc ( (totkeys + 1) * 80 + 1, 1); if (!(*header)) { *status = MEMORY_ALLOCATION; ffpmsg("failed to allocate memory to hold all the header keywords"); return(*status); } headptr = *header; casesn = FALSE; /* read every keyword */ for (ii = 1; ii <= totkeys; ii++) { ffgrec(fptr, ii, keybuf, status); /* pad record with blanks so that it is at least 80 chars long */ strcat(keybuf, " "); keyname[0] = '\0'; strncat(keyname, keybuf, 8); /* copy the keyword name */ if (exclude_comm) { if (!FSTRCMP("COMMENT ", keyname) || !FSTRCMP("HISTORY ", keyname) || !FSTRCMP(" ", keyname) ) continue; /* skip this commentary keyword */ } /* does keyword match any names in the exclusion list? */ for (jj = 0; jj < nexc; jj++ ) { ffcmps(exclist[jj], keyname, casesn, &match, &exact); if (match) break; } if (jj == nexc) { /* not in exclusion list, add this keyword to the string */ strcpy(headptr, keybuf); headptr += 80; (*nkeys)++; } } /* add the END keyword */ strcpy(headptr, "END "); headptr += 80; (*nkeys)++; *headptr = '\0'; /* terminate the header string */ /* * This is the actual fix for this function */ *header = realloc(*header, (*nkeys *80) + 1); /* minimize the allocated memory */ return(*status); } #endif /* End If WCS is installed */