diff options
author | William Joye <wjoye@cfa.harvard.edu> | 2016-10-26 21:13:00 (GMT) |
---|---|---|
committer | William Joye <wjoye@cfa.harvard.edu> | 2016-10-26 21:13:00 (GMT) |
commit | da2e3d212171bbe64c1af39114fd067308656990 (patch) | |
tree | 9601f7ed15fa1394762124630c12a792bc073ec2 /funtools/wcs/lin.c | |
parent | 76b109ad6d97d19ab835596dc70149ef379f3733 (diff) | |
download | blt-da2e3d212171bbe64c1af39114fd067308656990.zip blt-da2e3d212171bbe64c1af39114fd067308656990.tar.gz blt-da2e3d212171bbe64c1af39114fd067308656990.tar.bz2 |
rm funtools for update
Diffstat (limited to 'funtools/wcs/lin.c')
-rw-r--r-- | funtools/wcs/lin.c | 448 |
1 files changed, 0 insertions, 448 deletions
diff --git a/funtools/wcs/lin.c b/funtools/wcs/lin.c deleted file mode 100644 index 999a9d0..0000000 --- a/funtools/wcs/lin.c +++ /dev/null @@ -1,448 +0,0 @@ -/*============================================================================= -* -* WCSLIB - an implementation of the FITS WCS proposal. -* Copyright (C) 1995-2002, Mark Calabretta -* -* This library 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 2 of the License, or (at your option) any later version. -* -* This library 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 Public -* License along with this library; if not, write to the Free Software -* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -* -* Correspondence concerning WCSLIB may be directed to: -* Internet email: mcalabre@atnf.csiro.au -* Postal address: Dr. Mark Calabretta, -* Australia Telescope National Facility, -* P.O. Box 76, -* Epping, NSW, 2121, -* AUSTRALIA -* -*============================================================================= -* -* C routines which implement the FITS World Coordinate System (WCS) -* convention. -* -* Summary of routines -* ------------------- -* These utility routines apply the linear transformation defined by the WCS -* FITS header cards. There are separate routines for the image-to-pixel, -* linfwd(), and pixel-to-image, linrev(), transformations. -* -* An initialization routine, linset(), computes intermediate values from -* the transformation parameters but need not be called explicitly - see the -* explanation of lin.flag below. -* -* An auxiliary matrix inversion routine, matinv(), is included. It uses -* LU-triangular factorization with scaled partial pivoting. -* -* -* Initialization routine; linset() -* -------------------------------- -* Initializes members of a linprm data structure which hold intermediate -* values. Note that this routine need not be called directly; it will be -* invoked by linfwd() and linrev() if the "flag" structure member is -* anything other than a predefined magic value. -* -* Given and/or returned: -* lin linprm* Linear transformation parameters (see below). -* -* Function return value: -* int Error status -* 0: Success. -* 1: Memory allocation error. -* 2: PC matrix is singular. -* -* Forward transformation; linfwd() -* -------------------------------- -* Compute pixel coordinates from image coordinates. Note that where -* celestial coordinate systems are concerned the image coordinates -* correspond to (x,y) in the plane of projection, not celestial (lng,lat). -* -* Given: -* imgcrd const double[] -* Image (world) coordinate. -* -* Given and returned: -* lin linprm* Linear transformation parameters (see below). -* -* Returned: -* pixcrd d[] Pixel coordinate. -* -* Function return value: -* int Error status -* 0: Success. -* 1: The transformation is not invertible. -* -* Reverse transformation; linrev() -* -------------------------------- -* Compute image coordinates from pixel coordinates. Note that where -* celestial coordinate systems are concerned the image coordinates -* correspond to (x,y) in the plane of projection, not celestial (lng,lat). -* -* Given: -* pixcrd const double[] -* Pixel coordinate. -* -* Given and/or returned: -* lin linprm* Linear transformation parameters (see below). -* -* Returned: -* imgcrd d[] Image (world) coordinate. -* -* Function return value: -* int Error status -* 0: Success. -* 1: Error. -* -* Linear transformation parameters -* -------------------------------- -* The linprm struct consists of the following: -* -* int flag -* This flag must be set to zero whenever any of the following members -* are set or modified. This signals the initialization routine, -* linset(), to recompute intermediaries. -* int naxis -* Number of image axes. -* double *crpix -* Pointer to the first element of an array of double containing the -* coordinate reference pixel, CRPIXn. -* double *pc -* Pointer to the first element of the PC (pixel coordinate) -* transformation matrix. The expected order is -* -* lin.pc = {PC1_1, PC1_2, PC2_1, PC2_2}; -* -* This may be conveniently constructed from a two-dimensional array -* via -* -* double m[2][2] = {{PC1_1, PC1_2}, -* {PC2_1, PC2_2}}; -* -* which is equivalent to, -* -* double m[2][2]; -* m[0][0] = PC1_1; -* m[0][1] = PC1_2; -* m[1][0] = PC2_1; -* m[1][1] = PC2_2; -* -* for which the storage order is -* -* PC1_1, PC1_2, PC2_1, PC2_2 -* -* so it would be legitimate to set lin.pc = *m. -* double *cdelt -* Pointer to the first element of an array of double containing the -* coordinate increments, CDELTn. -* -* The remaining members of the linprm struct are maintained by the -* initialization routine and should not be modified. -* -* double *piximg -* Pointer to the first element of the matrix containing the product -* of the CDELTn diagonal matrix and the PC matrix. -* double *imgpix -* Pointer to the first element of the inverse of the piximg matrix. -* -* linset allocates storage for the above arrays using malloc(). Note, -* however, that these routines do not free this storage so if a linprm -* variable has itself been malloc'd then these structure members must be -* explicitly freed before the linprm variable is free'd otherwise a memory -* leak will result. -* -* Author: Mark Calabretta, Australia Telescope National Facility -* $Id: lin.c,v 2.8 2002/01/30 06:04:03 mcalabre Exp $ -*===========================================================================*/ - -#include <stdlib.h> -#include <math.h> -#include "wcslib.h" - -/* Map error number to error message for each function. */ -const char *linset_errmsg[] = { - 0, - "Memory allocation error", - "PC matrix is singular"}; - -const char *linfwd_errmsg[] = { - 0, - "Memory allocation error", - "PC matrix is singular"}; - -const char *linrev_errmsg[] = { - 0, - "Memory allocation error", - "PC matrix is singular"}; - -int linset(lin) - -struct linprm *lin; - -{ - int i, ij, j, mem, n; - - n = lin->naxis; - - /* Allocate memory for internal arrays. */ - mem = n * n * sizeof(double); - lin->piximg = (double*)malloc(mem); - if (lin->piximg == (double*)0) return 1; - - lin->imgpix = (double*)malloc(mem); - if (lin->imgpix == (double*)0) { - free(lin->piximg); - return 1; - } - - /* Compute the pixel-to-image transformation matrix. */ - for (i = 0, ij = 0; i < n; i++) { - for (j = 0; j < n; j++, ij++) { - lin->piximg[ij] = lin->cdelt[i] * lin->pc[ij]; - } - } - - /* Compute the image-to-pixel transformation matrix. */ - if (matinv(n, lin->piximg, lin->imgpix)) return 2; - - lin->flag = LINSET; - - return 0; -} - -/*--------------------------------------------------------------------------*/ - -int linfwd(imgcrd, lin, pixcrd) - -const double imgcrd[]; -struct linprm *lin; -double pixcrd[]; - -{ - int i, ij, j, n; - - n = lin->naxis; - - if (lin->flag != LINSET) { - if (linset(lin)) return 1; - } - - for (i = 0, ij = 0; i < n; i++) { - pixcrd[i] = 0.0; - for (j = 0; j < n; j++, ij++) { - pixcrd[i] += lin->imgpix[ij] * imgcrd[j]; - } - } - - for (j = 0; j < n; j++) { - pixcrd[j] += lin->crpix[j]; - } - - return 0; -} - -/*--------------------------------------------------------------------------*/ - -int linrev(pixcrd, lin, imgcrd) - -const double pixcrd[]; -struct linprm *lin; -double imgcrd[]; - -{ - int i, ij, j, n; - double temp; - - n = lin->naxis; - - if (lin->flag != LINSET) { - if (linset(lin)) return 1; - } - - for (i = 0; i < n; i++) { - imgcrd[i] = 0.0; - } - - for (j = 0; j < n; j++) { - temp = pixcrd[j] - lin->crpix[j]; - for (i = 0, ij = j; i < n; i++, ij+=n) { - imgcrd[i] += lin->piximg[ij] * temp; - } - } - - return 0; -} - -/*--------------------------------------------------------------------------*/ - -int matinv(n, mat, inv) - -const int n; -const double mat[]; -double inv[]; - -{ - register int i, ij, ik, j, k, kj, pj; - int itemp, mem, *mxl, *lxm, pivot; - double colmax, *lu, *rowmax, dtemp; - - - /* Allocate memory for internal arrays. */ - mem = n * sizeof(int); - if ((mxl = (int*)malloc(mem)) == (int*)0) return 1; - if ((lxm = (int*)malloc(mem)) == (int*)0) { - free(mxl); - return 1; - } - - mem = n * sizeof(double); - if ((rowmax = (double*)malloc(mem)) == (double*)0) { - free(mxl); - free(lxm); - return 1; - } - - mem *= n; - if ((lu = (double*)malloc(mem)) == (double*)0) { - free(mxl); - free(lxm); - free(rowmax); - return 1; - } - - - /* Initialize arrays. */ - for (i = 0, ij = 0; i < n; i++) { - /* Vector which records row interchanges. */ - mxl[i] = i; - - rowmax[i] = 0.0; - - for (j = 0; j < n; j++, ij++) { - dtemp = fabs(mat[ij]); - if (dtemp > rowmax[i]) rowmax[i] = dtemp; - - lu[ij] = mat[ij]; - } - - /* A row of zeroes indicates a singular matrix. */ - if (rowmax[i] == 0.0) { - free(mxl); - free(lxm); - free(rowmax); - free(lu); - return 2; - } - } - - - /* Form the LU triangular factorization using scaled partial pivoting. */ - for (k = 0; k < n; k++) { - /* Decide whether to pivot. */ - colmax = fabs(lu[k*n+k]) / rowmax[k]; - pivot = k; - - for (i = k+1; i < n; i++) { - ik = i*n + k; - dtemp = fabs(lu[ik]) / rowmax[i]; - if (dtemp > colmax) { - colmax = dtemp; - pivot = i; - } - } - - if (pivot > k) { - /* We must pivot, interchange the rows of the design matrix. */ - for (j = 0, pj = pivot*n, kj = k*n; j < n; j++, pj++, kj++) { - dtemp = lu[pj]; - lu[pj] = lu[kj]; - lu[kj] = dtemp; - } - - /* Amend the vector of row maxima. */ - dtemp = rowmax[pivot]; - rowmax[pivot] = rowmax[k]; - rowmax[k] = dtemp; - - /* Record the interchange for later use. */ - itemp = mxl[pivot]; - mxl[pivot] = mxl[k]; - mxl[k] = itemp; - } - - /* Gaussian elimination. */ - for (i = k+1; i < n; i++) { - ik = i*n + k; - - /* Nothing to do if lu[ik] is zero. */ - if (lu[ik] != 0.0) { - /* Save the scaling factor. */ - lu[ik] /= lu[k*n+k]; - - /* Subtract rows. */ - for (j = k+1; j < n; j++) { - lu[i*n+j] -= lu[ik]*lu[k*n+j]; - } - } - } - } - - - /* mxl[i] records which row of mat corresponds to row i of lu. */ - /* lxm[i] records which row of lu corresponds to row i of mat. */ - for (i = 0; i < n; i++) { - lxm[mxl[i]] = i; - } - - - /* Determine the inverse matrix. */ - for (i = 0, ij = 0; i < n; i++) { - for (j = 0; j < n; j++, ij++) { - inv[ij] = 0.0; - } - } - - for (k = 0; k < n; k++) { - inv[lxm[k]*n+k] = 1.0; - - /* Forward substitution. */ - for (i = lxm[k]+1; i < n; i++) { - for (j = lxm[k]; j < i; j++) { - inv[i*n+k] -= lu[i*n+j]*inv[j*n+k]; - } - } - - /* Backward substitution. */ - for (i = n-1; i >= 0; i--) { - for (j = i+1; j < n; j++) { - inv[i*n+k] -= lu[i*n+j]*inv[j*n+k]; - } - inv[i*n+k] /= lu[i*n+i]; - } - } - - free(mxl); - free(lxm); - free(rowmax); - free(lu); - - return 0; -} -/* Dec 20 1999 Doug Mink - Include wcslib.h, which includes lin.h - * - * Feb 15 2001 Doug Mink - Add comments for WCSLIB 2.6; no code changes - * Sep 19 2001 Doug Mink - Add above change to WCSLIB 2.7 code - * Nov 20 2001 Doug Mink - Always include stdlib.h - * - * Jan 15 2002 Bill Joye - Add ifdef so this compiles on MacOS/X - * - * Nov 18 2003 Doug Mink - Include stdlib.h instead of malloc.h - */ |