/*** 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()
*/