diff options
Diffstat (limited to 'wcs/imio.c')
| -rw-r--r-- | wcs/imio.c | 1543 |
1 files changed, 1543 insertions, 0 deletions
diff --git a/wcs/imio.c b/wcs/imio.c new file mode 100644 index 0000000..9bbf5a4 --- /dev/null +++ b/wcs/imio.c @@ -0,0 +1,1543 @@ +/*** File wcslib/imio.c + *** October 30, 2012 + *** By Jessica Mink, jmink@cfa.harvard.edu + *** Harvard-Smithsonian Center for Astrophysics + *** Copyright (C) 1996-2012 + *** Smithsonian Astrophysical Observatory, Cambridge, MA, USA + + 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 WCSTools should be addressed as follows: + Internet email: jmink@cfa.harvard.edu + Postal address: Jessica Mink + Smithsonian Astrophysical Observatory + 60 Garden St. + Cambridge, MA 02138 USA + + * Module: imio.c (image pixel manipulation) + * Purpose: Read and write pixels from arbitrary data type 2D arrays + * Subroutine: getpix (image, bitpix, w, h, bz, bs, x, y) + * Read pixel from 2D image of any numeric type (0,0 lower left) + * Subroutine: getpix1 (image, bitpix, w, h, bz, bs, x, y) + * Read pixel from 2D image of any numeric type (1,1 lower left) + * Subroutine: putpix (image, bitpix, w, h, bz, bs, x, y, dpix) + * Write pixel into 2D image of any numeric type (0,0 lower left) + * Subroutine: putpix1 (image, bitpix, w, h, bz, bs, x, y, dpix) + * Write pixel into 2D image of any numeric type (1,1 lower left) + * Subroutine: addpix (image, bitpix, w, h, bz, bs, x, y, dpix) + * Copy pixel into 2D image of any numeric type (0,0 lower left) + * Subroutine: addpix1 (image, bitpix, w, h, bz, bs, x, y, dpix) + * Add pixel into 2D image of any numeric type (1,1 lower left) + * Subroutine: maxvec (image, bitpix, bz, bs, pix1, npix) + * Get maximum of vector from 2D image of any numeric type + * Subroutine: minvec (image, bitpix, bz, bs, pix1, npix) + * Get minimum of vector from 2D image of any numeric type + * Subroutine: getvec (image, bitpix, bz, bs, pix1, npix, dvec) + * Get vector from 2D image of any numeric type + * Subroutine: putvec (image, bitpix, bz, bs, pix1, npix, dvec) + * Copy pixel vector into a vector of any numeric type + * Subroutine: addvec (image, bitpix, bz, bs, pix1, npix, dpix) + * Add constant to pixel values in a vector + * Subroutine: multvec (image, bitpix, bz, bs, pix1, npix, dpix) + * Multiply pixel values in a vector by a constant + * Subroutine: fillvec (image, bitpix, bz, bs, pix1, npix, dpix) + * Copy pixel value in a vector of any numeric type + * Subroutine: fillvec1 (image, bitpix, bz, bs, pix1, npix, dpix) + * Copy pixel value int a vector of any numeric type + * Subroutine: movepix (image1, bitpix, w1, x1, y1, image2, w2, x2, y2) + * Copy pixel from one image location to another + * Subroutine: imswap (bitpix,string,nbytes) + * Swap bytes in string in place, with FITS bits/pixel code + * Subroutine: imswap2 (string,nbytes) + * Swap bytes in string in place + * Subroutine imswap4 (string,nbytes) + * Reverse bytes of Integer*4 or Real*4 vector in place + * Subroutine imswap8 (string,nbytes) + * Reverse bytes of Real*8 vector in place + * Subroutine imswapped () + * Return 1 if PC/DEC byte order, else 0 + */ + +#include <stdlib.h> +#include <stdio.h> +#include "fitsfile.h" + +static int scale = 1; /* If 0, skip scaling step */ +void +setscale (scale0) +int scale0; +{scale = scale0; return;} + +/* GETPIX1 -- Get pixel from 2D FITS image of any numeric type */ + +double +getpix1 (image, bitpix, w, h, bzero, bscale, x, y) + +char *image; /* Image array as 1-D vector */ +int bitpix; /* FITS bits per pixel */ + /* 16 = short, -16 = unsigned short, 32 = int */ + /* -32 = float, -64 = double */ +int w; /* Image width in pixels */ +int h; /* Image height in pixels */ +double bzero; /* Zero point for pixel scaling */ +double bscale; /* Scale factor for pixel scaling */ +int x; /* One-based horizontal pixel number */ +int y; /* One-based vertical pixel number */ + +{ + return (getpix (image, bitpix, w, h, bzero, bscale, x-1, y-1)); +} + + +/* GETPIX -- Get pixel from 2D image of any numeric type */ + +double +getpix (image, bitpix, w, h, bzero, bscale, x, y) + +char *image; /* Image array as 1-D vector */ +int bitpix; /* FITS bits per pixel */ + /* 16 = short, -16 = unsigned short, 32 = int */ + /* -32 = float, -64 = double */ +int w; /* Image width in pixels */ +int h; /* Image height in pixels */ +double bzero; /* Zero point for pixel scaling */ +double bscale; /* Scale factor for pixel scaling */ +int x; /* Zero-based horizontal pixel number */ +int y; /* Zero-based vertical pixel number */ + +{ + short *im2; + int *im4; + unsigned char *im1; + unsigned short *imu; + float *imr; + double *imd; + double dpix; + +/* Return 0 if coordinates are not inside image */ + if (x < 0 || x >= w) + return (0.0); + if (y < 0 || y >= h) + return (0.0); + +/* Extract pixel from appropriate type of array */ + switch (bitpix) { + + case 8: + im1 = (unsigned char *)image; + dpix = (double) im1[(y*w) + x]; + break; + + case 16: + im2 = (short *)image; + dpix = (double) im2[(y*w) + x]; + break; + + case 32: + im4 = (int *)image; + dpix = (double) im4[(y*w) + x]; + break; + + case -16: + imu = (unsigned short *)image; + dpix = (double) imu[(y*w) + x]; + break; + + case -32: + imr = (float *)image; + dpix = (double) imr[(y*w) + x]; + break; + + case -64: + imd = (double *)image; + dpix = imd[(y*w) + x]; + break; + + default: + dpix = 0.0; + } + if (scale) + return (bzero + (bscale * dpix)); + else + return (dpix); +} + + +/* PUTPIX1 -- Copy pixel into 2D FITS image of any numeric type */ + +void +putpix1 (image, bitpix, w, h, bzero, bscale, x, y, dpix) + +char *image; +int bitpix; /* Number of bits per pixel */ + /* 16 = short, -16 = unsigned short, 32 = int */ + /* -32 = float, -64 = double */ +int w; /* Image width in pixels */ +int h; /* Image height in pixels */ +double bzero; /* Zero point for pixel scaling */ +double bscale; /* Scale factor for pixel scaling */ +int x; /* One-based horizontal pixel number */ +int y; /* One-based vertical pixel number */ +double dpix; + +{ + putpix (image, bitpix, w, h, bzero, bscale, x-1, y-1, dpix); + return; +} + + +/* PUTPIX -- Copy pixel into 2D image of any numeric type */ + +void +putpix (image, bitpix, w, h, bzero, bscale, x, y, dpix) + +char *image; +int bitpix; /* Number of bits per pixel */ + /* 16 = short, -16 = unsigned short, 32 = int */ + /* -32 = float, -64 = double */ +int w; /* Image width in pixels */ +int h; /* Image height in pixels */ +double bzero; /* Zero point for pixel scaling */ +double bscale; /* Scale factor for pixel scaling */ +int x; +int y; +double dpix; + +{ + double *imd; + float *imr; + int *im4; + short *im2; + unsigned short *imu; + unsigned char *im1; + +/* Return if coordinates are not inside image */ + if (x < 0 || x >= w) + return; + if (y < 0 || y >= h) + return; + + if (scale) + dpix = (dpix - bzero) / bscale; + + switch (bitpix) { + + case 8: + im1 = (unsigned char *)image; + if (dpix < 0) + im1[(y*w) + x] = (unsigned char) (dpix - 0.5); + else + im1[(y*w) + x] = (unsigned char) (dpix + 0.5); + break; + + case 16: + im2 = (short *)image; + if (dpix < 0) + im2[(y*w) + x] = (short) (dpix - 0.5); + else + im2[(y*w) + x] = (short) (dpix + 0.5); + break; + + case 32: + im4 = (int *)image; + if (dpix < 0) + im4[(y*w) + x] = (int) (dpix - 0.5); + else + im4[(y*w) + x] = (int) (dpix + 0.5); + break; + + case -16: + imu = (unsigned short *)image; + if (dpix < 0) + imu[(y*w) + x] = (unsigned short) 0; + else + imu[(y*w) + x] = (unsigned short) (dpix + 0.5); + break; + + case -32: + imr = (float *)image; + imr[(y*w) + x] = (float) dpix; + break; + + case -64: + imd = (double *)image; + imd[(y*w) + x] = dpix; + break; + + } + return; +} + + +/* ADDPIX1 -- Add pixel value into 2D FITS image of any numeric type */ + +void +addpix1 (image, bitpix, w, h, bzero, bscale, x, y, dpix) + +char *image; +int bitpix; /* Number of bits per pixel */ + /* 16 = short, -16 = unsigned short, 32 = int */ + /* -32 = float, -64 = double */ +int w; /* Image width in pixels */ +int h; /* Image height in pixels */ +double bzero; /* Zero point for pixel scaling */ +double bscale; /* Scale factor for pixel scaling */ +int x; /* One-based horizontal pixel number */ +int y; /* One-based vertical pixel number */ +double dpix; /* Value to add to pixel */ + +{ + addpix (image, bitpix, w, h, bzero, bscale, x-1, y-1, dpix); + return; +} + + +/* ADDPIX -- Add constant to pixel values in 2D image of any numeric type */ + +void +addpix (image, bitpix, w, h, bzero, bscale, x, y, dpix) + +char *image; +int bitpix; /* Number of bits per pixel */ + /* 16 = short, -16 = unsigned short, 32 = int */ + /* -32 = float, -64 = double */ +int w; /* Image width in pixels */ +int h; /* Image height in pixels */ +double bzero; /* Zero point for pixel scaling */ +double bscale; /* Scale factor for pixel scaling */ +int x; /* Zero-based horizontal pixel number */ +int y; /* Zero-based vertical pixel number */ +double dpix; /* Value to add to pixel */ + +{ + double *imd; + float *imr; + int *im4; + short *im2; + unsigned short *imu; + unsigned char *im1; + int ipix; + +/* Return if coordinates are not inside image */ + if (x < 0 || x >= w) + return; + if (y < 0 || y >= h) + return; + + if (scale) + dpix = (dpix - bzero) / bscale; + ipix = (y * w) + x; + + switch (bitpix) { + + case 8: + im1 = (unsigned char *)image; + if (dpix < 0) + image[ipix] = im1[ipix] + (unsigned char) (dpix - 0.5); + else + image[ipix] = im1[ipix] + (unsigned char) (dpix + 0.5); + break; + + case 16: + im2 = (short *)image; + if (dpix < 0) + im2[ipix] = im2[ipix] + (short) (dpix - 0.5); + else + im2[ipix] = im2[ipix] + (short) (dpix + 0.5); + break; + + case 32: + im4 = (int *)image; + if (dpix < 0) + im4[ipix] = im4[ipix] + (int) (dpix - 0.5); + else + im4[ipix] = im4[ipix] + (int) (dpix + 0.5); + break; + + case -16: + imu = (unsigned short *)image; + if (dpix > 0) + imu[ipix] = imu[ipix] + (unsigned short) (dpix + 0.5); + break; + + case -32: + imr = (float *)image; + imr[ipix] = imr[ipix] + (float) dpix; + break; + + case -64: + imd = (double *)image; + imd[ipix] = imd[ipix] + dpix; + break; + + } + return; +} + + +/* MOVEPIX -- Copy pixel between images */ + +void +movepix (image1, bitpix1, w1, x1, y1, image2, bitpix2, w2, x2, y2) + +char *image1; /* Pointer to first pixel in input image */ +int bitpix1; /* Bits per input pixel (FITS codes) */ + /* 16 = short, -16 = unsigned short, 32 = int */ + /* -32 = float, -64 = double */ +int w1; /* Number of horizontal pixels in input image */ +int x1, y1; /* Row and column for input pixel */ + +char *image2; /* Pointer to first pixel in output image */ +int bitpix2; /* Bits per output pixel (FITS codes) */ + /* 16 = short, -16 = unsigned short, 32 = int */ + /* -32 = float, -64 = double */ +int w2; /* Number of horizontal pixels in output image */ +int x2, y2; /* Row and column for output pixel */ + +{ + double dpix, *imd1, *imd2; + float rpix, *imr1, *imr2; + int *imi1, *imi2; + short *ims1, *ims2; + unsigned short *imu1, *imu2; + unsigned char *imc1, *imc2; + + if (x1 < 0 || x2 < 0 || x1 >= w1 || x2 >= w2) + return; + if (y1 < 0 || y2 < 0) + return; + + switch (bitpix1) { + + case 8: + imc1 = (unsigned char *)image1; + switch (bitpix2) { + case 8: + imc2 = (unsigned char *)image2; + imc2[(y2*w2) + x2] = imc1[(y1*w1) + x1]; + break; + case 16: + ims2 = (short *)image2; + ims2[(y2*w2) + x2] = (short) imc1[(y1*w1) + x1]; + break; + case 32: + imi2 = (int *)image2; + imi2[(y2*w2) + x2] = (int) imc1[(y1*w1) + x1]; + break; + case -16: + imu2 = (unsigned short *)image2; + imu2[(y2*w2) + x2] = (unsigned short) imc1[(y1*w1) + x1]; + break; + case -32: + imr2 = (float *)image2; + imr2[(y2*w2) + x2] = (float) imc1[(y1*w1) + x1]; + break; + case -64: + imd2 = (double *)image2; + imd2[(y2*w2) + x2] = (double) imc1[(y1*w1) + x1]; + break; + } + break; + + case 16: + ims1 = (short *)image1; + switch (bitpix2) { + case 8: + imc2 = (unsigned char *)image1; + imc2[(y2*w2) + x2] = (unsigned char) ims1[(y1*w1) + x1]; + break; + case 16: + ims2 = (short *)image2; + ims2[(y2*w2) + x2] = ims1[(y1*w1) + x1]; + break; + case 32: + imi2 = (int *)image2; + imi2[(y2*w2) + x2] = (int) ims1[(y1*w1) + x1]; + break; + case -16: + imu2 = (unsigned short *)image2; + imu2[(y2*w2) + x2] = (unsigned short) ims1[(y1*w1) + x1]; + break; + case -32: + imr2 = (float *)image2; + imr2[(y2*w2) + x2] = (float) ims1[(y1*w1) + x1]; + break; + case -64: + imd2 = (double *)image2; + imd2[(y2*w2) + x2] = (double) ims1[(y1*w1) + x1]; + break; + } + break; + + case 32: + imi1 = (int *)image1; + switch (bitpix2) { + case 8: + imc2 = (unsigned char *)image2; + imc2[(y2*w2) + x2] = (unsigned char) imi1[(y1*w1) + x1]; + break; + case 16: + ims2 = (short *)image2; + ims2[(y2*w2) + x2] = (short) imi1[(y1*w1) + x1]; + break; + case 32: + imi2 = (int *)image2; + imi2[(y2*w2) + x2] = imi1[(y1*w1) + x1]; + break; + case -16: + imu2 = (unsigned short *)image2; + imu2[(y2*w2) + x2] = (unsigned short) imi1[(y1*w1) + x1]; + break; + case -32: + imr2 = (float *)image2; + imr2[(y2*w2) + x2] = (float) imi1[(y1*w1) + x1]; + break; + case -64: + imd2 = (double *)image2; + imd2[(y2*w2) + x2] = (double) imi1[(y1*w1) + x1]; + break; + } + break; + + case -16: + imu1 = (unsigned short *)image1; + switch (bitpix2) { + case 8: + imc2 = (unsigned char *)image2; + imc2[(y2*w2) + x2] = (unsigned char) imu1[(y1*w1) + x1]; + break; + case 16: + ims2 = (short *)image2; + ims2[(y2*w2) + x2] = (short) imu1[(y1*w1) + x1]; + break; + case 32: + imi2 = (int *)image2; + imi2[(y2*w2) + x2] = (int) imu1[(y1*w1) + x1]; + break; + case -16: + imu2 = (unsigned short *)image2; + imu2[(y2*w2) + x2] = imu1[(y1*w1) + x1]; + break; + case -32: + imr2 = (float *)image2; + imr2[(y2*w2) + x2] = (float) imu1[(y1*w1) + x1]; + break; + case -64: + imd2 = (double *)image2; + imd2[(y2*w2) + x2] = (double) imu1[(y1*w1) + x1]; + break; + } + break; + + case -32: + imr1 = (float *)image1; + rpix = imr1[(y1*w1) + x1]; + switch (bitpix2) { + case 8: + imc2 = (unsigned char *)image2; + if (rpix < 0.0) + imc2[(y2*w2) + x2] = (unsigned char) 0; + else + imc2[(y2*w2) + x2] = (unsigned char) (rpix + 0.5); + break; + case 16: + ims2 = (short *)image2; + if (rpix < 0.0) + ims2[(y2*w2) + x2] = (short) (rpix - 0.5); + else + ims2[(y2*w2) + x2] = (short) (rpix + 0.5); + break; + case 32: + imi2 = (int *)image2; + if (rpix < 0.0) + imi2[(y2*w2) + x2] = (int) (rpix - 0.5); + else + imi2[(y2*w2) + x2] = (int) (rpix + 0.5); + break; + case -16: + imu2 = (unsigned short *)image2; + if (rpix < 0.0) + imu2[(y2*w2) + x2] = (unsigned short) 0; + else + imu2[(y2*w2) + x2] = (unsigned short) (rpix + 0.5); + break; + case -32: + imr2 = (float *)image2; + imr2[(y2*w2) + x2] = rpix; + break; + case -64: + imd2 = (double *)image2; + imd2[(y2*w2) + x2] = (double) rpix; + break; + } + break; + + case -64: + imd1 = (double *)image1; + dpix = imd1[(y1*w1) + x1]; + switch (bitpix2) { + case 8: + imc2 = (unsigned char *)image2; + if (dpix < 0.0) + imc2[(y2*w2) + x2] = (unsigned char) 0; + else + imc2[(y2*w2) + x2] = (unsigned char) (dpix + 0.5); + break; + case 16: + ims2 = (short *)image2; + if (dpix < 0.0) + ims2[(y2*w2) + x2] = (short) (dpix - 0.5); + else + ims2[(y2*w2) + x2] = (short) (dpix + 0.5); + break; + case 32: + imi2 = (int *)image2; + if (dpix < 0.0) + imi2[(y2*w2) + x2] = (int) (dpix - 0.5); + else + imi2[(y2*w2) + x2] = (int) (dpix + 0.5); + break; + case -16: + imu2 = (unsigned short *)image2; + if (dpix < 0.0) + imu2[(y2*w2) + x2] = (unsigned short) 0; + else + imu2[(y2*w2) + x2] = (unsigned short) (dpix + 0.5); + break; + case -32: + imr2 = (float *)image2; + imr2[(y2*w2) + x2] = (float) dpix; + break; + case -64: + imd2 = (double *)image2; + imd2[(y2*w2) + x2] = dpix; + break; + } + break; + } + return; +} + + +/* MAXVEC -- Get maximum value in vector from 2D image of any numeric type */ + +double +maxvec (image, bitpix, bzero, bscale, pix1, npix) + +char *image; /* Image array from which to read vector */ +int bitpix; /* Number of bits per pixel in image */ + /* 16 = short, -16 = unsigned short, 32 = int */ + /* -32 = float, -64 = double */ +double bzero; /* Zero point for pixel scaling */ +double bscale; /* Scale factor for pixel scaling */ +int pix1; /* Offset of first pixel to check */ +int npix; /* Number of pixels to check */ + +{ + short *im2, imax2, ip2; + int *im4, imax4, ip4; + unsigned short *imu, imaxu, ipu; + float *imr, imaxr, ipr; + double *imd; + double dmax = 0.0; + double ipd; + int ipix, pix2; + unsigned char *imc, imaxc, ipc; + + pix2 = pix1 + npix; + + switch (bitpix) { + + case 8: + imc = (unsigned char *)(image); + imaxc = *(imc + pix1); + for (ipix = pix1; ipix < pix2; ipix++) { + ipc = *(imc + ipix); + if (ipc > imaxc) + imaxc = ipc; + } + dmax = (double) imaxc; + break; + + case 16: + im2 = (short *)image; + imax2 = *(im2 + pix1); + for (ipix = pix1; ipix < pix2; ipix++) { + ip2 = *(im2 + ipix); + if (ip2 > imax2) + imax2 = ip2; + } + dmax = (double) imax2; + break; + + case 32: + im4 = (int *)image; + imax4 = *(im4 + pix1); + for (ipix = pix1; ipix < pix2; ipix++) { + ip4 = *(im4 + ipix); + if (ip4 > imax4) + imax4 = ip4; + } + dmax = (double) imax4; + break; + + case -16: + imu = (unsigned short *)image; + imaxu = *(imu + pix1); + for (ipix = pix1; ipix < pix2; ipix++) { + ipu = *(imu + ipix); + if (ipu > imaxu) + imaxu = ipu; + } + dmax = (double) imaxu; + break; + + case -32: + imr = (float *)image; + imaxr = *(imr + pix1); + for (ipix = pix1; ipix < pix2; ipix++) { + ipr = *(imr + ipix); + if (ipr > imaxr) + imax2 = ipr; + } + dmax = (double) imaxr; + break; + + case -64: + imd = (double *)image; + dmax = *(imd + pix1); + for (ipix = pix1; ipix < pix2; ipix++) { + ipd = *(imd + ipix); + if (ipd > dmax) + dmax = ipd; + } + break; + + } + + /* Scale data if either BZERO or BSCALE keyword has been set */ + if (scale && (bzero != 0.0 || bscale != 1.0)) + dmax = (dmax * bscale) + bzero; + + return (dmax); +} + + +/* MINVEC -- Get minimum value in vector from 2D image of any numeric type */ + +double +minvec (image, bitpix, bzero, bscale, pix1, npix) + +char *image; /* Image array from which to read vector */ +int bitpix; /* Number of bits per pixel in image */ + /* 16 = short, -16 = unsigned short, 32 = int */ + /* -32 = float, -64 = double */ +double bzero; /* Zero point for pixel scaling */ +double bscale; /* Scale factor for pixel scaling */ +int pix1; /* Offset of first pixel to check */ +int npix; /* Number of pixels to check */ + +{ + short *im2, imin2, ip2; + int *im4, imin4, ip4; + unsigned short *imu, iminu, ipu; + float *imr, iminr, ipr; + double *imd, ipd; + double dmin = 0.0; + int ipix, pix2; + unsigned char *imc, cmin, cp; + + pix2 = pix1 + npix; + + switch (bitpix) { + + case 8: + imc = (unsigned char *)image; + cmin = *(imc + pix1); + for (ipix = pix1; ipix < pix2; ipix++) { + cp = *(imc + ipix); + if (cp < cmin) + cmin = cp; + } + dmin = (double) cmin; + break; + + case 16: + im2 = (short *)image + pix1; + imin2 = *im2; + for (ipix = pix1; ipix < pix2; ipix++) { + ip2 = *(im2 + ipix); + if (ip2 < imin2) + imin2 = ip2; + } + dmin = (double) imin2; + break; + + case 32: + im4 = (int *)image; + imin4 = *(im4 + pix1); + for (ipix = pix1; ipix < pix2; ipix++) { + ip4 = *(im4 + ipix); + if (ip4 < imin4) + imin4 = ip4; + } + dmin = (double) imin4; + break; + + case -16: + imu = (unsigned short *)image; + iminu = *(imu + pix1); + for (ipix = pix1; ipix < pix2; ipix++) { + ipu = *(imu + ipix); + if (ipu < iminu) + iminu = ipu; + } + dmin = (double) iminu; + break; + + case -32: + imr = (float *)image; + iminr = *(imr + pix1); + for (ipix = pix1; ipix < pix2; ipix++) { + ipr = *(imr + ipix); + if (ipr < iminr) + iminr = ipr; + } + dmin = (double) iminr; + break; + + case -64: + imd = (double *)image; + dmin = *(imd + pix1); + for (ipix = pix1; ipix < pix2; ipix++) { + ipd = *(imd + ipix); + if (ipd < dmin) + dmin = ipd; + } + break; + + } + + /* Scale data if either BZERO or BSCALE keyword has been set */ + if (scale && (bzero != 0.0 || bscale != 1.0)) + dmin = (dmin * bscale) + bzero; + + return (dmin); +} + + +/* ADDVEC -- Add constant to pixel values in 2D image of any numeric type */ + +void +addvec (image, bitpix, bzero, bscale, pix1, npix, dpix) + +char *image; /* Image array from which to extract vector */ +int bitpix; /* Number of bits per pixel in image */ + /* 16 = short, -16 = unsigned short, 32 = int */ + /* -32 = float, -64 = double */ +double bzero; /* Zero point for pixel scaling */ +double bscale; /* Scale factor for pixel scaling */ +int pix1; /* Offset of first pixel to extract */ +int npix; /* Number of pixels to extract */ +double dpix; /* Value to add to pixels */ + +{ + unsigned char *imc, ccon; + short *im2, jcon; + int *im4, icon; + unsigned short *imu, ucon; + float *imr, rcon; + double *imd; + int ipix, pix2; + + pix2 = pix1 + npix; + + if (scale) + dpix = (dpix - bzero) / bscale; + + switch (bitpix) { + + case 8: + imc = (unsigned char *) (image + pix1); + if (dpix < 0) + ccon = (unsigned char) (dpix - 0.5); + else + ccon = (unsigned char) (dpix + 0.5); + for (ipix = pix1; ipix < pix2; ipix++) + *imc++ += ccon; + break; + + case 16: + im2 = (short *) (image + pix1); + if (dpix < 0) + jcon = (short) (dpix - 0.5); + else + jcon = (short) (dpix + 0.5); + for (ipix = pix1; ipix < pix2; ipix++) + *im2++ += jcon; + break; + + case 32: + im4 = (int *) (image + pix1); + if (dpix < 0) + icon = (int) (dpix - 0.5); + else + icon = (int) (dpix + 0.5); + for (ipix = pix1; ipix < pix2; ipix++) + *im4++ += icon; + break; + + case -16: + imu = (unsigned short *) (image + pix1); + if (dpix > 0) { + ucon = (unsigned short) (dpix + 0.5); + imu = (unsigned short *) (image + pix1); + for (ipix = pix1; ipix < pix2; ipix++) + *imu++ += ucon; + } + else { + icon = (int) (dpix - 0.5); + imu = (unsigned short *) (image + pix1); + for (ipix = pix1; ipix < pix2; ipix++) { + unsigned short tmp = (icon + (int) *imu); + *imu++ += tmp; + } + } + break; + + case -32: + rcon = (float) dpix; + imr = (float *) (image + pix1); + for (ipix = pix1; ipix < pix2; ipix++) + *imr++ += rcon; + break; + + case -64: + imd = (double *) (image + pix1); + for (ipix = pix1; ipix < pix2; ipix++) + *imd++ += dpix; + break; + } + return; +} + + +/* MULTVEC -- Multiply pixel values in place in 2D image of any numeric type */ + +void +multvec (image, bitpix, bzero, bscale, pix1, npix, dpix) + +char *image; /* Image array from which to extract vector */ +int bitpix; /* Number of bits per pixel in image */ + /* 16 = short, -16 = unsigned short, 32 = int */ + /* -32 = float, -64 = double */ +double bzero; /* Zero point for pixel scaling */ +double bscale; /* Scale factor for pixel scaling */ +int pix1; /* Offset of first pixel to extract */ +int npix; /* Number of pixels to extract */ +double dpix; /* Value by which to multiply pixels */ + +{ + char *imc, ccon; + short *im2, jcon; + int *im4, icon, isint; + unsigned short *imu, ucon; + float *imr, rcon; + double *imd, dcon, dval; + int ipix, pix2; + + pix2 = pix1 + npix; + + if (scale) + dpix = (dpix - bzero) / bscale; + ipix = (int) dpix; + dcon = (double) ipix; + if (dcon == dpix) + isint = 1; + else + isint = 0; + + switch (bitpix) { + + case 8: + imc = image + pix1; + if (isint) { + if (dpix < 0) + ccon = (char) (dpix - 0.5); + else + ccon = (char) (dpix + 0.5); + for (ipix = pix1; ipix < pix2; ipix++) + *imc++ *= ccon; + } + else { + for (ipix = pix1; ipix < pix2; ipix++) { + dval = ((double) *imc) * dpix; + if (dval < 256.0) + *imc++ = (char) dval; + else + *imc++ = (char) 255; + } + } + break; + + case 16: + im2 = (short *) (image + pix1); + if (isint) { + im2 = (short *)image; + if (dpix < 0) + jcon = (short) (dpix - 0.5); + else + jcon = (short) (dpix + 0.5); + for (ipix = pix1; ipix < pix2; ipix++) + *im2++ *= jcon; + } + else { + for (ipix = pix1; ipix < pix2; ipix++) { + dval = ((double) *im2) * dpix; + if (dval < 32768.0) + *im2++ = (short) dval; + else + *im2++ = (short) 32767; + } + } + break; + + case 32: + im4 = (int *) (image + pix1); + if (isint) { + if (dpix < 0) + icon = (int) (dpix - 0.5); + else + icon = (int) (dpix + 0.5); + for (ipix = pix1; ipix < pix2; ipix++) + *im4++ *= icon; + } + else { + for (ipix = pix1; ipix < pix2; ipix++) { + dval = ((double) *im4) * dpix; + if (dval < 32768.0) + *im4++ = (int) dval; + else + *im4++ = (int) 32767; + } + } + break; + + case -16: + imu = (unsigned short *) (image + pix1); + if (dpix > 0) { + ucon = (unsigned short) (dpix + 0.5); + imu = (unsigned short *) (image + pix1); + for (ipix = pix1; ipix < pix2; ipix++) + *imu++ *= ucon; + } + break; + + case -32: + rcon = (float) dpix; + imr = (float *) (image + pix1); + for (ipix = pix1; ipix < pix2; ipix++) + *imr++ *= rcon; + break; + + case -64: + imd = (double *) (image + pix1); + for (ipix = pix1; ipix < pix2; ipix++) + *imd++ *= dpix; + break; + + } + return; +} + + +/* GETVEC -- Get vector from 2D image of any numeric type */ + +void +getvec (image, bitpix, bzero, bscale, pix1, npix, dvec0) + +char *image; /* Image array from which to extract vector */ +int bitpix; /* Number of bits per pixel in image */ + /* 16 = short, -16 = unsigned short, 32 = int */ + /* -32 = float, -64 = double */ +double bzero; /* Zero point for pixel scaling */ +double bscale; /* Scale factor for pixel scaling */ +int pix1; /* Offset of first pixel to extract */ +int npix; /* Number of pixels to extract */ +double *dvec0; /* Vector of pixels (returned) */ + +{ + short *im2; + int *im4; + unsigned short *imu; + float *imr; + double *imd; + double *dvec; + int ipix, pix2; + + pix2 = pix1 + npix; + dvec = dvec0; + + switch (bitpix) { + + case 8: + for (ipix = pix1; ipix < pix2; ipix++) + *dvec++ = (double) *(image + ipix); + break; + + case 16: + im2 = (short *)image; + for (ipix = pix1; ipix < pix2; ipix++) + *dvec++ = (double) *(im2 + ipix); + break; + + case 32: + im4 = (int *)image; + for (ipix = pix1; ipix < pix2; ipix++) + *dvec++ = (double) *(im4 + ipix); + break; + + case -16: + imu = (unsigned short *)image; + for (ipix = pix1; ipix < pix2; ipix++) + *dvec++ = (double) *(imu + ipix); + break; + + case -32: + imr = (float *)image; + for (ipix = pix1; ipix < pix2; ipix++) + *dvec++ = (double) *(imr + ipix); + break; + + case -64: + imd = (double *)image; + for (ipix = pix1; ipix < pix2; ipix++) + *dvec++ = (double) *(imd + ipix); + break; + + } + + /* Scale data if either BZERO or BSCALE keyword has been set */ + if (scale && (bzero != 0.0 || bscale != 1.0)) { + dvec = dvec0; + for (ipix = pix1; ipix < pix2; ipix++) { + *dvec = (*dvec * bscale) + bzero; + dvec++; + } + } + + return; +} + + +/* PUTVEC -- Copy pixel vector into 2D image of any numeric type */ + +void +putvec (image, bitpix, bzero, bscale, pix1, npix, dvec) + +char *image; /* Image into which to copy vector */ +int bitpix; /* Number of bits per pixel im image */ + /* 16 = short, -16 = unsigned short, 32 = int */ + /* -32 = float, -64 = double */ +double bzero; /* Zero point for pixel scaling */ +double bscale; /* Scale factor for pixel scaling */ +int pix1; /* Offset of first pixel of vector in image */ +int npix; /* Number of pixels to copy */ +double *dvec; /* Vector of pixels to copy */ + +{ + short *im2; + int *im4; + unsigned short *imu; + float *imr; + double *imd; + int ipix, pix2; + double *dp = dvec; + + pix2 = pix1 + npix; + + /* Scale data if either BZERO or BSCALE keyword has been set */ + if (scale && (bzero != 0.0 || bscale != 1.0)) { + for (ipix = pix1; ipix < pix2; ipix++) { + *dp = (*dp - bzero) / bscale; + dp++; + } + dp = dvec; + } + + switch (bitpix) { + + case 8: + for (ipix = pix1; ipix < pix2; ipix++) + *(image+ipix) = (char) *dp++; + break; + + case 16: + im2 = (short *)image; + for (ipix = pix1; ipix < pix2; ipix++) { + if (*dp < 0.0) + *(im2+ipix) = (short) (*dp++ - 0.5); + else + *(im2+ipix) = (short) (*dp++ + 0.5); + } + break; + + case 32: + im4 = (int *)image; + for (ipix = pix1; ipix < pix2; ipix++) { + if (*dp < 0.0) + *(im4+ipix) = (int) (*dp++ - 0.5); + else + *(im4+ipix) = (int) (*dp++ + 0.5); + } + break; + + case -16: + imu = (unsigned short *)image; + for (ipix = pix1; ipix < pix2; ipix++) { + if (*dp < 0.0) + *(imu+ipix) = (unsigned short) 0; + else + *(imu+ipix) = (unsigned short) (*dp++ + 0.5); + } + break; + + case -32: + imr = (float *)image; + for (ipix = pix1; ipix < pix2; ipix++) + *(imr+ipix) = (float) *dp++; + break; + + case -64: + imd = (double *)image; + for (ipix = pix1; ipix < pix2; ipix++) + *(imd+ipix) = (double) *dp++; + break; + } + return; +} + + +/* FILLVEC1 -- Copy single value into a vector of any numeric type */ + +void +fillvec1 (image, bitpix, bzero, bscale, pix1, npix, dpix) + +char *image; /* Vector to fill */ +int bitpix; /* Number of bits per pixel im image */ + /* 16 = short, -16 = unsigned short, 32 = int */ + /* -32 = float, -64 = double */ +double bzero; /* Zero point for pixel scaling */ +double bscale; /* Scale factor for pixel scaling */ +int pix1; /* First pixel to fill */ +int npix; /* Number of pixels to fill */ +double dpix; /* Value with which to fill pixels */ +{ + fillvec (image, bitpix, bzero, bscale, pix1-1, npix, dpix); + return; +} + + +/* FILLVEC -- Copy single value into a vector of any numeric type */ + +void +fillvec (image, bitpix, bzero, bscale, pix1, npix, dpix) + +char *image; /* Vector to fill */ +int bitpix; /* Number of bits per pixel im image */ + /* 16 = short, -16 = unsigned short, 32 = int */ + /* -32 = float, -64 = double */ +double bzero; /* Zero point for pixel scaling */ +double bscale; /* Scale factor for pixel scaling */ +int pix1; /* First pixel to fill */ +int npix; /* Number of pixels to fill */ +double dpix; /* Value with which to fill pixels */ +{ + char ipc; + short *im2, ip2; + int *im4, ip4; + unsigned short *imu, ipu; + float *imr, ipr; + double *imd; + int ipix, pix2; + double dp; + + pix2 = pix1 + npix; + + /* Scale data if either BZERO or BSCALE keyword has been set */ + dp = dpix; + if (scale && (bzero != 0.0 || bscale != 1.0)) + dp = (dp - bzero) / bscale; + + switch (bitpix) { + + case 8: + if (dp < 0.0) + ipc = (char) (dp - 0.5); + else + ipc = (char) (dp + 0.5); + for (ipix = pix1; ipix < pix2; ipix++) + image[ipix] = ipc; + break; + + case 16: + im2 = (short *)image; + if (dp < 0.0) + ip2 = (short) (dp - 0.5); + else + ip2 = (short) (dp + 0.5); + for (ipix = pix1; ipix < pix2; ipix++) + im2[ipix] = ip2; + break; + + case 32: + im4 = (int *)image; + if (dp < 0.0) + ip4 = (int) (dp - 0.5); + else + ip4 = (int) (dp + 0.5); + for (ipix = pix1; ipix < pix2; ipix++) + im4[ipix] = ip4; + break; + + case -16: + imu = (unsigned short *)image; + if (dp < 0.0) + ipu = (unsigned short) (dp - 0.5); + else + ipu = (unsigned short) (dp + 0.5); + for (ipix = pix1; ipix < pix2; ipix++) + imu[ipix] = ipu; + break; + + case -32: + imr = (float *)image; + ipr = (float) dp; + for (ipix = pix1; ipix < pix2; ipix++) + imr[ipix] = ipr; + break; + + case -64: + imd = (double *)image; + for (ipix = pix1; ipix < pix2; ipix++) + imd[ipix] = dp; + break; + } + return; +} + + +/* IMSWAP -- Reverse bytes of any type of vector in place */ + +void +imswap (bitpix, string, nbytes) + +int bitpix; /* Number of bits per pixel */ + /* 16 = short, -16 = unsigned short, 32 = int */ + /* -32 = float, -64 = double */ +char *string; /* Address of starting point of bytes to swap */ +int nbytes; /* Number of bytes to swap */ + +{ + switch (bitpix) { + + case 8: + break; + + case 16: + if (nbytes < 2) return; + imswap2 (string,nbytes); + break; + + case 32: + if (nbytes < 4) return; + imswap4 (string,nbytes); + break; + + case -16: + if (nbytes < 2) return; + imswap2 (string,nbytes); + break; + + case -32: + if (nbytes < 4) return; + imswap4 (string,nbytes); + break; + + case -64: + if (nbytes < 8) return; + imswap8 (string,nbytes); + break; + + } + return; +} + + +/* IMSWAP2 -- Swap bytes in string in place */ + +void +imswap2 (string,nbytes) + + +char *string; /* Address of starting point of bytes to swap */ +int nbytes; /* Number of bytes to swap */ + +{ + char *sbyte, temp, *slast; + + slast = string + nbytes; + sbyte = string; + while (sbyte < slast) { + temp = sbyte[0]; + sbyte[0] = sbyte[1]; + sbyte[1] = temp; + sbyte= sbyte + 2; + } + return; +} + + +/* IMSWAP4 -- Reverse bytes of Integer*4 or Real*4 vector in place */ + +void +imswap4 (string,nbytes) + +char *string; /* Address of Integer*4 or Real*4 vector */ +int nbytes; /* Number of bytes to reverse */ + +{ + char *sbyte, *slast; + char temp0, temp1, temp2, temp3; + + slast = string + nbytes; + sbyte = string; + while (sbyte < slast) { + temp3 = sbyte[0]; + temp2 = sbyte[1]; + temp1 = sbyte[2]; + temp0 = sbyte[3]; + sbyte[0] = temp0; + sbyte[1] = temp1; + sbyte[2] = temp2; + sbyte[3] = temp3; + sbyte = sbyte + 4; + } + + return; +} + + +/* IMSWAP8 -- Reverse bytes of Real*8 vector in place */ + +void +imswap8 (string,nbytes) + +char *string; /* Address of Real*8 vector */ +int nbytes; /* Number of bytes to reverse */ + +{ + char *sbyte, *slast; + char temp[8]; + + slast = string + nbytes; + sbyte = string; + while (sbyte < slast) { + temp[7] = sbyte[0]; + temp[6] = sbyte[1]; + temp[5] = sbyte[2]; + temp[4] = sbyte[3]; + temp[3] = sbyte[4]; + temp[2] = sbyte[5]; + temp[1] = sbyte[6]; + temp[0] = sbyte[7]; + sbyte[0] = temp[0]; + sbyte[1] = temp[1]; + sbyte[2] = temp[2]; + sbyte[3] = temp[3]; + sbyte[4] = temp[4]; + sbyte[5] = temp[5]; + sbyte[6] = temp[6]; + sbyte[7] = temp[7]; + sbyte = sbyte + 8; + } + return; +} + +/* IMSWAPPED -- Returns 0 if big-endian (Sun,Mac), + 1 if little-endian(PC,Alpha) */ + +int +imswapped () + +{ + char *ctest; + int itest; + + itest = 1; + ctest = (char *)&itest; + if (*ctest) + return (1); + else + return (0); +} + +/* Apr 17 1996 New file + * May 22 1996 Add H so that PUTPIX and GETPIX can check coordinates + * Jun 11 1996 Simplify NEWIMAGE subroutine + * Jun 12 1996 Add byte-swapping subroutines + * + * Jul 24 1997 Add 8-bit option to subroutines + * + * May 27 1998 Include imio.h instead of fitshead.h + * Jun 17 1998 Fix bug, changing all unsigned int's to unsigned short's + * + * Apr 29 1999 Add scaling to getpix, putpix, getvec, and putvec + * Apr 29 1999 Fix bug in getvec in dealing with 1-byte data + * Sep 14 1999 Change dp incrementing so it works on Alpha compiler + * Sep 27 1999 Add interface for 1-based (FITS) image access + * Sep 27 1999 Add addpix() and addpix1() + * Dec 14 1999 In putpix(), addpix(), putvec(), round when output is integer + * + * Sep 20 2000 In getvec(), scale only if necessary + * + * Nov 27 2001 In movepix(), add char to char move + * + * Jan 23 2002 Add global scale switch to turn off scaling + * Jun 4 2002 In getvec() and putvec(), change dpix to dvec + * Jun 4 2002 Add addvec() to add to a vector + * Jul 19 2002 Fix getvec() bug rescaling scaled numbers + * + * May 20 2003 Declare scale0 in setscale() + * + * Jan 28 2004 Add image limit check to movepix() + * Feb 27 2004 Add fillvec() and fillvec1() to set vector to a constant + * + * Jun 27 2005 Fix major bug in fillvec(); pass value dpix in fillvec1(), too + * Aug 18 2005 Add maxvec(), addvec(), and multvec() + * + * Mar 1 2006 Fix bug of occasional double application of bscale in getvec() + * Apr 3 2006 Fix bad cast in unisigned int section of addvec() + * May 3 2006 Code fixes in addpix and multpix suggested by Robert Lupton + * Jun 8 2006 Drop erroneous second im2 assignment without offset in addvec() + * Jun 20 2006 Fix typos masquerading as unitialized variables + * + * Jan 8 2007 Include fitsfile.h instead of imio.h + * Jun 11 2007 Add minvec() and speed up maxvec() + * + * Apr 12 2012 Fix 8-bit variables to be unsigned char + * Oct 19 2012 Fix errors with character images in minvec() and maxvec() + * Oct 30 2012 Fix errors with short images in minvec() and maxvec() + */ |