1 /*** File wcslib/imio.c
2  *** October 30, 2012
3  *** By Jessica Mink, jmink@cfa.harvard.edu
4  *** Harvard-Smithsonian Center for Astrophysics
5  *** Copyright (C) 1996-2012
6  *** Smithsonian Astrophysical Observatory, Cambridge, MA, USA
7 
8     This library is free software; you can redistribute it and/or
9     modify it under the terms of the GNU Lesser General Public
10     License as published by the Free Software Foundation; either
11     version 2 of the License, or (at your option) any later version.
12 
13     This library is distributed in the hope that it will be useful,
14     but WITHOUT ANY WARRANTY; without even the implied warranty of
15     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
16     Lesser General Public License for more details.
17 
18     You should have received a copy of the GNU Lesser General Public
19     License along with this library; if not, write to the Free Software
20     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
21 
22     Correspondence concerning WCSTools should be addressed as follows:
23            Internet email: jmink@cfa.harvard.edu
24            Postal address: Jessica Mink
25                            Smithsonian Astrophysical Observatory
26                            60 Garden St.
27                            Cambridge, MA 02138 USA
28 
29  * Module:      imio.c (image pixel manipulation)
30  * Purpose:     Read and write pixels from arbitrary data type 2D arrays
31  * Subroutine:	getpix (image, bitpix, w, h, bz, bs, x, y)
32  *		Read pixel from 2D image of any numeric type (0,0 lower left)
33  * Subroutine:	getpix1 (image, bitpix, w, h, bz, bs, x, y)
34  *		Read pixel from 2D image of any numeric type (1,1 lower left)
35  * Subroutine:	putpix (image, bitpix, w, h, bz, bs, x, y, dpix)
36  *		Write pixel into 2D image of any numeric type (0,0 lower left)
37  * Subroutine:	putpix1 (image, bitpix, w, h, bz, bs, x, y, dpix)
38  *		Write pixel into 2D image of any numeric type (1,1 lower left)
39  * Subroutine:	addpix (image, bitpix, w, h, bz, bs, x, y, dpix)
40  *		Copy pixel into 2D image of any numeric type (0,0 lower left)
41  * Subroutine:	addpix1 (image, bitpix, w, h, bz, bs, x, y, dpix)
42  *		Add pixel into 2D image of any numeric type (1,1 lower left)
43  * Subroutine:	maxvec (image, bitpix, bz, bs, pix1, npix)
44  *		Get maximum of vector from 2D image of any numeric type
45  * Subroutine:	minvec (image, bitpix, bz, bs, pix1, npix)
46  *		Get minimum of vector from 2D image of any numeric type
47  * Subroutine:	getvec (image, bitpix, bz, bs, pix1, npix, dvec)
48  *		Get vector from 2D image of any numeric type
49  * Subroutine:	putvec (image, bitpix, bz, bs, pix1, npix, dvec)
50  *		Copy pixel vector into a vector of any numeric type
51  * Subroutine:	addvec (image, bitpix, bz, bs, pix1, npix, dpix)
52  *		Add constant to pixel values in a vector
53  * Subroutine:	multvec (image, bitpix, bz, bs, pix1, npix, dpix)
54  *		Multiply pixel values in a vector by a constant
55  * Subroutine:	fillvec (image, bitpix, bz, bs, pix1, npix, dpix)
56  *		Copy pixel value in a vector of any numeric type
57  * Subroutine:	fillvec1 (image, bitpix, bz, bs, pix1, npix, dpix)
58  *		Copy pixel value int a vector of any numeric type
59  * Subroutine:	movepix (image1, bitpix, w1, x1, y1, image2, w2, x2, y2)
60  *		Copy pixel from one image location to another
61  * Subroutine:	imswap (bitpix,string,nbytes)
62  *		Swap bytes in string in place, with FITS bits/pixel code
63  * Subroutine:	imswap2 (string,nbytes)
64  *		Swap bytes in string in place
65  * Subroutine	imswap4 (string,nbytes)
66  *		Reverse bytes of Integer*4 or Real*4 vector in place
67  * Subroutine	imswap8 (string,nbytes)
68  *		Reverse bytes of Real*8 vector in place
69  * Subroutine	imswapped ()
70  *		Return 1 if PC/DEC byte order, else 0
71  */
72 
73 #include <stdlib.h>
74 #include <stdio.h>
75 #include "fitsfile.h"
76 
77 static int scale = 1;	/* If 0, skip scaling step */
78 void
setscale(scale0)79 setscale (scale0)
80 int scale0;
81 {scale = scale0; return;}
82 
83 /* GETPIX1 -- Get pixel from 2D FITS image of any numeric type */
84 
85 double
getpix1(image,bitpix,w,h,bzero,bscale,x,y)86 getpix1 (image, bitpix, w, h, bzero, bscale, x, y)
87 
88 char	*image;		/* Image array as 1-D vector */
89 int	bitpix;		/* FITS bits per pixel */
90 			/*  16 = short, -16 = unsigned short, 32 = int */
91 			/* -32 = float, -64 = double */
92 int	w;		/* Image width in pixels */
93 int	h;		/* Image height in pixels */
94 double  bzero;		/* Zero point for pixel scaling */
95 double  bscale;		/* Scale factor for pixel scaling */
96 int	x;		/* One-based horizontal pixel number */
97 int	y;		/* One-based vertical pixel number */
98 
99 {
100     return (getpix (image, bitpix, w, h, bzero, bscale, x-1, y-1));
101 }
102 
103 
104 /* GETPIX -- Get pixel from 2D image of any numeric type */
105 
106 double
getpix(image,bitpix,w,h,bzero,bscale,x,y)107 getpix (image, bitpix, w, h, bzero, bscale, x, y)
108 
109 char	*image;		/* Image array as 1-D vector */
110 int	bitpix;		/* FITS bits per pixel */
111 			/*  16 = short, -16 = unsigned short, 32 = int */
112 			/* -32 = float, -64 = double */
113 int	w;		/* Image width in pixels */
114 int	h;		/* Image height in pixels */
115 double  bzero;		/* Zero point for pixel scaling */
116 double  bscale;		/* Scale factor for pixel scaling */
117 int	x;		/* Zero-based horizontal pixel number */
118 int	y;		/* Zero-based vertical pixel number */
119 
120 {
121     short *im2;
122     int *im4;
123     unsigned char *im1;
124     unsigned short *imu;
125     float *imr;
126     double *imd;
127     double dpix;
128 
129 /* Return 0 if coordinates are not inside image */
130     if (x < 0 || x >= w)
131 	return (0.0);
132     if (y < 0 || y >= h)
133 	return (0.0);
134 
135 /* Extract pixel from appropriate type of array */
136     switch (bitpix) {
137 
138 	case 8:
139 	  im1 = (unsigned char *)image;
140 	  dpix = (double) im1[(y*w) + x];
141 	  break;
142 
143 	case 16:
144 	  im2 = (short *)image;
145 	  dpix = (double) im2[(y*w) + x];
146 	  break;
147 
148 	case 32:
149 	  im4 = (int *)image;
150 	  dpix = (double) im4[(y*w) + x];
151 	  break;
152 
153 	case -16:
154 	  imu = (unsigned short *)image;
155 	  dpix = (double) imu[(y*w) + x];
156 	  break;
157 
158 	case -32:
159 	  imr = (float *)image;
160 	  dpix = (double) imr[(y*w) + x];
161 	  break;
162 
163 	case -64:
164 	  imd = (double *)image;
165 	  dpix = imd[(y*w) + x];
166 	  break;
167 
168 	default:
169 	  dpix = 0.0;
170 	}
171     if (scale)
172 	return (bzero + (bscale * dpix));
173     else
174 	return (dpix);
175 }
176 
177 
178 /* PUTPIX1 -- Copy pixel into 2D FITS image of any numeric type */
179 
180 void
putpix1(image,bitpix,w,h,bzero,bscale,x,y,dpix)181 putpix1 (image, bitpix, w, h, bzero, bscale, x, y, dpix)
182 
183 char	*image;
184 int	bitpix;		/* Number of bits per pixel */
185 			/*  16 = short, -16 = unsigned short, 32 = int */
186 			/* -32 = float, -64 = double */
187 int	w;		/* Image width in pixels */
188 int	h;		/* Image height in pixels */
189 double  bzero;		/* Zero point for pixel scaling */
190 double  bscale;		/* Scale factor for pixel scaling */
191 int	x;		/* One-based horizontal pixel number */
192 int	y;		/* One-based vertical pixel number */
193 double	dpix;
194 
195 {
196     putpix (image, bitpix, w, h, bzero, bscale, x-1, y-1, dpix);
197     return;
198 }
199 
200 
201 /* PUTPIX -- Copy pixel into 2D image of any numeric type */
202 
203 void
putpix(image,bitpix,w,h,bzero,bscale,x,y,dpix)204 putpix (image, bitpix, w, h, bzero, bscale, x, y, dpix)
205 
206 char	*image;
207 int	bitpix;		/* Number of bits per pixel */
208 			/*  16 = short, -16 = unsigned short, 32 = int */
209 			/* -32 = float, -64 = double */
210 int	w;		/* Image width in pixels */
211 int	h;		/* Image height in pixels */
212 double  bzero;		/* Zero point for pixel scaling */
213 double  bscale;		/* Scale factor for pixel scaling */
214 int	x;
215 int	y;
216 double	dpix;
217 
218 {
219     double *imd;
220     float *imr;
221     int *im4;
222     short *im2;
223     unsigned short *imu;
224     unsigned char *im1;
225 
226 /* Return if coordinates are not inside image */
227     if (x < 0 || x >= w)
228 	return;
229     if (y < 0 || y >= h)
230 	return;
231 
232     if (scale)
233 	dpix = (dpix - bzero) / bscale;
234 
235     switch (bitpix) {
236 
237 	case 8:
238 	    im1 = (unsigned char *)image;
239 	    if (dpix < 0)
240 		im1[(y*w) + x] = (unsigned char) (dpix - 0.5);
241 	    else
242 		im1[(y*w) + x] = (unsigned char) (dpix + 0.5);
243 	    break;
244 
245 	case 16:
246 	    im2 = (short *)image;
247 	    if (dpix < 0)
248 		im2[(y*w) + x] = (short) (dpix - 0.5);
249 	    else
250 		im2[(y*w) + x] = (short) (dpix + 0.5);
251 	    break;
252 
253 	case 32:
254 	    im4 = (int *)image;
255 	    if (dpix < 0)
256 		im4[(y*w) + x] = (int) (dpix - 0.5);
257 	    else
258 		im4[(y*w) + x] = (int) (dpix + 0.5);
259 	    break;
260 
261 	case -16:
262 	    imu = (unsigned short *)image;
263 	    if (dpix < 0)
264 		imu[(y*w) + x] = (unsigned short) 0;
265 	    else
266 		imu[(y*w) + x] = (unsigned short) (dpix + 0.5);
267 	    break;
268 
269 	case -32:
270 	    imr = (float *)image;
271 	    imr[(y*w) + x] = (float) dpix;
272 	    break;
273 
274 	case -64:
275 	    imd = (double *)image;
276 	    imd[(y*w) + x] = dpix;
277 	    break;
278 
279 	}
280     return;
281 }
282 
283 
284 /* ADDPIX1 -- Add pixel value into 2D FITS image of any numeric type */
285 
286 void
addpix1(image,bitpix,w,h,bzero,bscale,x,y,dpix)287 addpix1 (image, bitpix, w, h, bzero, bscale, x, y, dpix)
288 
289 char	*image;
290 int	bitpix;		/* Number of bits per pixel */
291 			/*  16 = short, -16 = unsigned short, 32 = int */
292 			/* -32 = float, -64 = double */
293 int	w;		/* Image width in pixels */
294 int	h;		/* Image height in pixels */
295 double  bzero;		/* Zero point for pixel scaling */
296 double  bscale;		/* Scale factor for pixel scaling */
297 int	x;		/* One-based horizontal pixel number */
298 int	y;		/* One-based vertical pixel number */
299 double	dpix;		/* Value to add to pixel */
300 
301 {
302     addpix (image, bitpix, w, h, bzero, bscale, x-1, y-1, dpix);
303     return;
304 }
305 
306 
307 /* ADDPIX -- Add constant to pixel values in 2D image of any numeric type */
308 
309 void
addpix(image,bitpix,w,h,bzero,bscale,x,y,dpix)310 addpix (image, bitpix, w, h, bzero, bscale, x, y, dpix)
311 
312 char	*image;
313 int	bitpix;		/* Number of bits per pixel */
314 			/*  16 = short, -16 = unsigned short, 32 = int */
315 			/* -32 = float, -64 = double */
316 int	w;		/* Image width in pixels */
317 int	h;		/* Image height in pixels */
318 double  bzero;		/* Zero point for pixel scaling */
319 double  bscale;		/* Scale factor for pixel scaling */
320 int	x;		/* Zero-based horizontal pixel number */
321 int	y;		/* Zero-based vertical pixel number */
322 double	dpix;		/* Value to add to pixel */
323 
324 {
325     double *imd;
326     float *imr;
327     int *im4;
328     short *im2;
329     unsigned short *imu;
330     unsigned char *im1;
331     int ipix;
332 
333 /* Return if coordinates are not inside image */
334     if (x < 0 || x >= w)
335 	return;
336     if (y < 0 || y >= h)
337 	return;
338 
339     if (scale)
340 	dpix = (dpix - bzero) / bscale;
341     ipix = (y * w) + x;
342 
343     switch (bitpix) {
344 
345 	case 8:
346 	    im1 = (unsigned char *)image;
347 	    if (dpix < 0)
348 		image[ipix] = im1[ipix] + (unsigned char) (dpix - 0.5);
349 	    else
350 		image[ipix] = im1[ipix] + (unsigned char) (dpix + 0.5);
351 	    break;
352 
353 	case 16:
354 	    im2 = (short *)image;
355 	    if (dpix < 0)
356 		im2[ipix] = im2[ipix] + (short) (dpix - 0.5);
357 	    else
358 		im2[ipix] = im2[ipix] + (short) (dpix + 0.5);
359 	    break;
360 
361 	case 32:
362 	    im4 = (int *)image;
363 	    if (dpix < 0)
364 		im4[ipix] = im4[ipix] + (int) (dpix - 0.5);
365 	    else
366 		im4[ipix] = im4[ipix] + (int) (dpix + 0.5);
367 	    break;
368 
369 	case -16:
370 	    imu = (unsigned short *)image;
371 	    if (dpix > 0)
372 		imu[ipix] = imu[ipix] + (unsigned short) (dpix + 0.5);
373 	    break;
374 
375 	case -32:
376 	    imr = (float *)image;
377 	    imr[ipix] = imr[ipix] + (float) dpix;
378 	    break;
379 
380 	case -64:
381 	    imd = (double *)image;
382 	    imd[ipix] = imd[ipix] + dpix;
383 	    break;
384 
385 	}
386     return;
387 }
388 
389 
390 /* MOVEPIX -- Copy pixel between images */
391 
392 void
movepix(image1,bitpix1,w1,x1,y1,image2,bitpix2,w2,x2,y2)393 movepix (image1, bitpix1, w1, x1, y1, image2, bitpix2, w2, x2, y2)
394 
395 char	*image1;	/* Pointer to first pixel in input image */
396 int	bitpix1;	/* Bits per input pixel (FITS codes) */
397 			/*  16 = short, -16 = unsigned short, 32 = int */
398 			/* -32 = float, -64 = double */
399 int	w1;		/* Number of horizontal pixels in input image */
400 int	x1, y1;		/* Row and column for input pixel */
401 
402 char	*image2;	/* Pointer to first pixel in output image */
403 int	bitpix2;	/* Bits per output pixel (FITS codes) */
404 			/*  16 = short, -16 = unsigned short, 32 = int */
405 			/* -32 = float, -64 = double */
406 int	w2;		/* Number of horizontal pixels in output image */
407 int	x2, y2;		/* Row and column for output pixel */
408 
409 {
410     double dpix, *imd1, *imd2;
411     float rpix, *imr1, *imr2;
412     int *imi1, *imi2;
413     short *ims1, *ims2;
414     unsigned short *imu1, *imu2;
415     unsigned char *imc1, *imc2;
416 
417     if (x1 < 0 || x2 < 0 || x1 >= w1 || x2 >= w2)
418 	return;
419     if (y1 < 0 || y2 < 0)
420 	return;
421 
422     switch (bitpix1) {
423 
424 	case 8:
425 	    imc1 = (unsigned char *)image1;
426 	    switch (bitpix2) {
427 		case 8:
428 		    imc2 = (unsigned char *)image2;
429 		    imc2[(y2*w2) + x2] = imc1[(y1*w1) + x1];
430 		    break;
431 		case 16:
432 		    ims2 = (short *)image2;
433 		    ims2[(y2*w2) + x2] = (short) imc1[(y1*w1) + x1];
434 		    break;
435 		case 32:
436 		    imi2 = (int *)image2;
437 		    imi2[(y2*w2) + x2] = (int) imc1[(y1*w1) + x1];
438 		    break;
439 		case -16:
440 		    imu2 = (unsigned short *)image2;
441 		    imu2[(y2*w2) + x2] = (unsigned short) imc1[(y1*w1) + x1];
442 		    break;
443 		case -32:
444 		    imr2 = (float *)image2;
445 		    imr2[(y2*w2) + x2] = (float) imc1[(y1*w1) + x1];
446 		    break;
447 		case -64:
448 		    imd2 = (double *)image2;
449 		    imd2[(y2*w2) + x2] = (double) imc1[(y1*w1) + x1];
450 		    break;
451 		}
452 	    break;
453 
454 	case 16:
455 	    ims1 = (short *)image1;
456 	    switch (bitpix2) {
457 		case 8:
458 		    imc2 = (unsigned char *)image1;
459 		    imc2[(y2*w2) + x2] = (unsigned char) ims1[(y1*w1) + x1];
460 		    break;
461 		case 16:
462 		    ims2 = (short *)image2;
463 		    ims2[(y2*w2) + x2] = ims1[(y1*w1) + x1];
464 		    break;
465 		case 32:
466 		    imi2 = (int *)image2;
467 		    imi2[(y2*w2) + x2] = (int) ims1[(y1*w1) + x1];
468 		    break;
469 		case -16:
470 		    imu2 = (unsigned short *)image2;
471 		    imu2[(y2*w2) + x2] = (unsigned short) ims1[(y1*w1) + x1];
472 		    break;
473 		case -32:
474 		    imr2 = (float *)image2;
475 		    imr2[(y2*w2) + x2] = (float) ims1[(y1*w1) + x1];
476 		    break;
477 		case -64:
478 		    imd2 = (double *)image2;
479 		    imd2[(y2*w2) + x2] = (double) ims1[(y1*w1) + x1];
480 		    break;
481 		}
482 	    break;
483 
484 	case 32:
485 	    imi1 = (int *)image1;
486 	    switch (bitpix2) {
487 		case 8:
488 		    imc2 = (unsigned char *)image2;
489 		    imc2[(y2*w2) + x2] = (unsigned char) imi1[(y1*w1) + x1];
490 		    break;
491 		case 16:
492 		    ims2 = (short *)image2;
493 		    ims2[(y2*w2) + x2] = (short) imi1[(y1*w1) + x1];
494 		    break;
495 		case 32:
496 		    imi2 = (int *)image2;
497 		    imi2[(y2*w2) + x2] = imi1[(y1*w1) + x1];
498 		    break;
499 		case -16:
500 		    imu2 = (unsigned short *)image2;
501 		    imu2[(y2*w2) + x2] = (unsigned short) imi1[(y1*w1) + x1];
502 		    break;
503 		case -32:
504 		    imr2 = (float *)image2;
505 		    imr2[(y2*w2) + x2] = (float) imi1[(y1*w1) + x1];
506 		    break;
507 		case -64:
508 		    imd2 = (double *)image2;
509 		    imd2[(y2*w2) + x2] = (double) imi1[(y1*w1) + x1];
510 		    break;
511 		}
512 	    break;
513 
514 	case -16:
515 	    imu1 = (unsigned short *)image1;
516 	    switch (bitpix2) {
517 		case 8:
518 		    imc2 = (unsigned char *)image2;
519 		    imc2[(y2*w2) + x2] = (unsigned char) imu1[(y1*w1) + x1];
520 		    break;
521 		case 16:
522 		    ims2 = (short *)image2;
523 		    ims2[(y2*w2) + x2] = (short) imu1[(y1*w1) + x1];
524 		    break;
525 		case 32:
526 		    imi2 = (int *)image2;
527 		    imi2[(y2*w2) + x2] = (int) imu1[(y1*w1) + x1];
528 		    break;
529 		case -16:
530 		    imu2 = (unsigned short *)image2;
531 		    imu2[(y2*w2) + x2] = imu1[(y1*w1) + x1];
532 		    break;
533 		case -32:
534 		    imr2 = (float *)image2;
535 		    imr2[(y2*w2) + x2] = (float) imu1[(y1*w1) + x1];
536 		    break;
537 		case -64:
538 		    imd2 = (double *)image2;
539 		    imd2[(y2*w2) + x2] = (double) imu1[(y1*w1) + x1];
540 		    break;
541 		}
542 	    break;
543 
544 	case -32:
545 	    imr1 = (float *)image1;
546 	    rpix = imr1[(y1*w1) + x1];
547 	    switch (bitpix2) {
548 		case 8:
549 		    imc2 = (unsigned char *)image2;
550 		    if (rpix < 0.0)
551 			imc2[(y2*w2) + x2] = (unsigned char) 0;
552 		    else
553 			imc2[(y2*w2) + x2] = (unsigned char) (rpix + 0.5);
554 		    break;
555 		case 16:
556 		    ims2 = (short *)image2;
557 		    if (rpix < 0.0)
558 			ims2[(y2*w2) + x2] = (short) (rpix - 0.5);
559 		    else
560 			ims2[(y2*w2) + x2] = (short) (rpix + 0.5);
561 		    break;
562 		case 32:
563 		    imi2 = (int *)image2;
564 		    if (rpix < 0.0)
565 			imi2[(y2*w2) + x2] = (int) (rpix - 0.5);
566 		    else
567 			imi2[(y2*w2) + x2] = (int) (rpix + 0.5);
568 		    break;
569 		case -16:
570 		    imu2 = (unsigned short *)image2;
571 		    if (rpix < 0.0)
572 			imu2[(y2*w2) + x2] = (unsigned short) 0;
573 		    else
574 			imu2[(y2*w2) + x2] = (unsigned short) (rpix + 0.5);
575 		    break;
576 		case -32:
577 		    imr2 = (float *)image2;
578 		    imr2[(y2*w2) + x2] = rpix;
579 		    break;
580 		case -64:
581 		    imd2 = (double *)image2;
582 		    imd2[(y2*w2) + x2] = (double) rpix;
583 		    break;
584 		}
585 	    break;
586 
587 	case -64:
588 	    imd1 = (double *)image1;
589 	    dpix = imd1[(y1*w1) + x1];
590 	    switch (bitpix2) {
591 		case 8:
592 		    imc2 = (unsigned char *)image2;
593 		    if (dpix < 0.0)
594 			imc2[(y2*w2) + x2] = (unsigned char) 0;
595 		    else
596 			imc2[(y2*w2) + x2] = (unsigned char) (dpix + 0.5);
597 		    break;
598 		case 16:
599 		    ims2 = (short *)image2;
600 		    if (dpix < 0.0)
601 			ims2[(y2*w2) + x2] = (short) (dpix - 0.5);
602 		    else
603 			ims2[(y2*w2) + x2] = (short) (dpix + 0.5);
604 		    break;
605 		case 32:
606 		    imi2 = (int *)image2;
607 		    if (dpix < 0.0)
608 			imi2[(y2*w2) + x2] = (int) (dpix - 0.5);
609 		    else
610 			imi2[(y2*w2) + x2] = (int) (dpix + 0.5);
611 		    break;
612 		case -16:
613 		    imu2 = (unsigned short *)image2;
614 		    if (dpix < 0.0)
615 			imu2[(y2*w2) + x2] = (unsigned short) 0;
616 		    else
617 			imu2[(y2*w2) + x2] = (unsigned short) (dpix + 0.5);
618 		    break;
619 		case -32:
620 		    imr2 = (float *)image2;
621 		    imr2[(y2*w2) + x2] = (float) dpix;
622 		    break;
623 		case -64:
624 		    imd2 = (double *)image2;
625 		    imd2[(y2*w2) + x2] = dpix;
626 		    break;
627 		}
628 	    break;
629 	}
630     return;
631 }
632 
633 
634 /* MAXVEC -- Get maximum value in vector from 2D image of any numeric type */
635 
636 double
maxvec(image,bitpix,bzero,bscale,pix1,npix)637 maxvec (image, bitpix, bzero, bscale, pix1, npix)
638 
639 char	*image;		/* Image array from which to read vector */
640 int	bitpix;		/* Number of bits per pixel in image */
641 			/*  16 = short, -16 = unsigned short, 32 = int */
642 			/* -32 = float, -64 = double */
643 double  bzero;		/* Zero point for pixel scaling */
644 double  bscale;		/* Scale factor for pixel scaling */
645 int	pix1;		/* Offset of first pixel to check */
646 int	npix;		/* Number of pixels to check */
647 
648 {
649     short *im2, imax2, ip2;
650     int *im4, imax4, ip4;
651     unsigned short *imu, imaxu, ipu;
652     float *imr, imaxr, ipr;
653     double *imd;
654     double dmax = 0.0;
655     double ipd;
656     int ipix, pix2;
657     unsigned char *imc, imaxc, ipc;
658 
659     pix2 = pix1 + npix;
660 
661     switch (bitpix) {
662 
663 	case 8:
664 	    imc = (unsigned char *)(image);
665 	    imaxc = *(imc + pix1);
666 	    for (ipix = pix1; ipix < pix2; ipix++) {
667 		ipc = *(imc + ipix);
668 		if (ipc > imaxc)
669 		    imaxc = ipc;
670 		}
671 	    dmax = (double) imaxc;
672 	    break;
673 
674 	case 16:
675 	    im2 = (short *)image;
676 	    imax2 = *(im2 + pix1);
677 	    for (ipix = pix1; ipix < pix2; ipix++) {
678 		ip2 = *(im2 + ipix);
679 		if (ip2 > imax2)
680 		    imax2 = ip2;
681 		}
682 	    dmax = (double) imax2;
683 	    break;
684 
685 	case 32:
686 	    im4 = (int *)image;
687 	    imax4 = *(im4 + pix1);
688 	    for (ipix = pix1; ipix < pix2; ipix++) {
689 		ip4 = *(im4 + ipix);
690 		if (ip4 > imax4)
691 		    imax4 = ip4;
692 		}
693 	    dmax = (double) imax4;
694 	    break;
695 
696 	case -16:
697 	    imu = (unsigned short *)image;
698 	    imaxu = *(imu + pix1);
699 	    for (ipix = pix1; ipix < pix2; ipix++) {
700 		ipu = *(imu + ipix);
701 		if (ipu > imaxu)
702 		    imaxu = ipu;
703 		}
704 	    dmax = (double) imaxu;
705 	    break;
706 
707 	case -32:
708 	    imr = (float *)image;
709 	    imaxr = *(imr + pix1);
710 	    for (ipix = pix1; ipix < pix2; ipix++) {
711 		ipr = *(imr + ipix);
712 		if (ipr > imaxr)
713 		    imax2 = ipr;
714 		}
715 	    dmax = (double) imaxr;
716 	    break;
717 
718 	case -64:
719 	    imd = (double *)image;
720 	    dmax = *(imd + pix1);
721 	    for (ipix = pix1; ipix < pix2; ipix++) {
722 		ipd = *(imd + ipix);
723 		if (ipd > dmax)
724 		    dmax = ipd;
725 		}
726 	    break;
727 
728 	}
729 
730     /* Scale data if either BZERO or BSCALE keyword has been set */
731     if (scale && (bzero != 0.0 || bscale != 1.0))
732 	dmax = (dmax * bscale) + bzero;
733 
734     return (dmax);
735 }
736 
737 
738 /* MINVEC -- Get minimum value in vector from 2D image of any numeric type */
739 
740 double
minvec(image,bitpix,bzero,bscale,pix1,npix)741 minvec (image, bitpix, bzero, bscale, pix1, npix)
742 
743 char	*image;		/* Image array from which to read vector */
744 int	bitpix;		/* Number of bits per pixel in image */
745 			/*  16 = short, -16 = unsigned short, 32 = int */
746 			/* -32 = float, -64 = double */
747 double  bzero;		/* Zero point for pixel scaling */
748 double  bscale;		/* Scale factor for pixel scaling */
749 int	pix1;		/* Offset of first pixel to check */
750 int	npix;		/* Number of pixels to check */
751 
752 {
753     short *im2, imin2, ip2;
754     int *im4, imin4, ip4;
755     unsigned short *imu, iminu, ipu;
756     float *imr, iminr, ipr;
757     double *imd, ipd;
758     double dmin = 0.0;
759     int ipix, pix2;
760     unsigned char *imc, cmin, cp;
761 
762     pix2 = pix1 + npix;
763 
764     switch (bitpix) {
765 
766 	case 8:
767 	    imc = (unsigned char *)image;
768 	    cmin = *(imc + pix1);
769 	    for (ipix = pix1; ipix < pix2; ipix++) {
770 		cp = *(imc + ipix);
771 		if (cp < cmin)
772 		    cmin = cp;
773 		}
774 	    dmin = (double) cmin;
775 	    break;
776 
777 	case 16:
778 	    im2 = (short *)image + pix1;
779 	    imin2 = *im2;
780 	    for (ipix = pix1; ipix < pix2; ipix++) {
781 		ip2 = *(im2 + ipix);
782 		if (ip2 < imin2)
783 		    imin2 = ip2;
784 		}
785 	    dmin = (double) imin2;
786 	    break;
787 
788 	case 32:
789 	    im4 = (int *)image;
790 	    imin4 = *(im4 + pix1);
791 	    for (ipix = pix1; ipix < pix2; ipix++) {
792 		ip4 = *(im4 + ipix);
793 		if (ip4 < imin4)
794 		    imin4 = ip4;
795 		}
796 	    dmin = (double) imin4;
797 	    break;
798 
799 	case -16:
800 	    imu = (unsigned short *)image;
801 	    iminu = *(imu + pix1);
802 	    for (ipix = pix1; ipix < pix2; ipix++) {
803 		ipu = *(imu + ipix);
804 		if (ipu < iminu)
805 		    iminu = ipu;
806 		}
807 	    dmin = (double) iminu;
808 	    break;
809 
810 	case -32:
811 	    imr = (float *)image;
812 	    iminr = *(imr + pix1);
813 	    for (ipix = pix1; ipix < pix2; ipix++) {
814 		ipr = *(imr + ipix);
815 		if (ipr < iminr)
816 		    iminr = ipr;
817 		}
818 	    dmin = (double) iminr;
819 	    break;
820 
821 	case -64:
822 	    imd = (double *)image;
823 	    dmin = *(imd + pix1);
824 	    for (ipix = pix1; ipix < pix2; ipix++) {
825 		ipd = *(imd + ipix);
826 		if (ipd < dmin)
827 		    dmin = ipd;
828 		}
829 	    break;
830 
831 	}
832 
833     /* Scale data if either BZERO or BSCALE keyword has been set */
834     if (scale && (bzero != 0.0 || bscale != 1.0))
835 	dmin = (dmin * bscale) + bzero;
836 
837     return (dmin);
838 }
839 
840 
841 /* ADDVEC -- Add constant to pixel values in 2D image of any numeric type */
842 
843 void
addvec(image,bitpix,bzero,bscale,pix1,npix,dpix)844 addvec (image, bitpix, bzero, bscale, pix1, npix, dpix)
845 
846 char	*image;		/* Image array from which to extract vector */
847 int	bitpix;		/* Number of bits per pixel in image */
848 			/*  16 = short, -16 = unsigned short, 32 = int */
849 			/* -32 = float, -64 = double */
850 double  bzero;		/* Zero point for pixel scaling */
851 double  bscale;		/* Scale factor for pixel scaling */
852 int	pix1;		/* Offset of first pixel to extract */
853 int	npix;		/* Number of pixels to extract */
854 double	dpix;		/* Value to add to pixels */
855 
856 {
857     unsigned char *imc, ccon;
858     short *im2, jcon;
859     int *im4, icon;
860     unsigned short *imu, ucon;
861     float *imr, rcon;
862     double *imd;
863     int ipix, pix2;
864 
865     pix2 = pix1 + npix;
866 
867     if (scale)
868 	dpix = (dpix - bzero) / bscale;
869 
870     switch (bitpix) {
871 
872 	case 8:
873 	    imc = (unsigned char *) (image + pix1);
874 	    if (dpix < 0)
875 		ccon = (unsigned char) (dpix - 0.5);
876 	    else
877 		ccon = (unsigned char) (dpix + 0.5);
878 	    for (ipix = pix1; ipix < pix2; ipix++)
879 		*imc++ += ccon;
880 	    break;
881 
882 	case 16:
883 	    im2 = (short *) (image + pix1);
884 	    if (dpix < 0)
885 		jcon = (short) (dpix - 0.5);
886 	    else
887 		jcon = (short) (dpix + 0.5);
888 	    for (ipix = pix1; ipix < pix2; ipix++)
889 		*im2++ += jcon;
890 	    break;
891 
892 	case 32:
893 	    im4 = (int *) (image + pix1);
894 	    if (dpix < 0)
895 		icon = (int) (dpix - 0.5);
896 	    else
897 		icon = (int) (dpix + 0.5);
898 	    for (ipix = pix1; ipix < pix2; ipix++)
899 		*im4++ += icon;
900 	    break;
901 
902 	case -16:
903 	    imu = (unsigned short *) (image + pix1);
904 	    if (dpix > 0) {
905 		ucon = (unsigned short) (dpix + 0.5);
906 		imu = (unsigned short *) (image + pix1);
907 		for (ipix = pix1; ipix < pix2; ipix++)
908 		    *imu++ += ucon;
909 		}
910 	    else {
911 		icon = (int) (dpix - 0.5);
912 		imu = (unsigned short *) (image + pix1);
913 		for (ipix = pix1; ipix < pix2; ipix++) {
914 		    unsigned short tmp = (icon + (int) *imu);
915 		    *imu++ += tmp;
916 		    }
917 		}
918 	    break;
919 
920 	case -32:
921 	    rcon = (float) dpix;
922 	    imr = (float *) (image + pix1);
923 	    for (ipix = pix1; ipix < pix2; ipix++)
924 		*imr++ += rcon;
925 	    break;
926 
927 	case -64:
928 	    imd = (double *) (image + pix1);
929 	    for (ipix = pix1; ipix < pix2; ipix++)
930 		*imd++ += dpix;
931 	    break;
932 	}
933     return;
934 }
935 
936 
937 /* MULTVEC -- Multiply pixel values in place in 2D image of any numeric type */
938 
939 void
multvec(image,bitpix,bzero,bscale,pix1,npix,dpix)940 multvec (image, bitpix, bzero, bscale, pix1, npix, dpix)
941 
942 char	*image;		/* Image array from which to extract vector */
943 int	bitpix;		/* Number of bits per pixel in image */
944 			/*  16 = short, -16 = unsigned short, 32 = int */
945 			/* -32 = float, -64 = double */
946 double  bzero;		/* Zero point for pixel scaling */
947 double  bscale;		/* Scale factor for pixel scaling */
948 int	pix1;		/* Offset of first pixel to extract */
949 int	npix;		/* Number of pixels to extract */
950 double	dpix;		/* Value by which to multiply pixels */
951 
952 {
953     char *imc, ccon;
954     short *im2, jcon;
955     int *im4, icon, isint;
956     unsigned short *imu, ucon;
957     float *imr, rcon;
958     double *imd, dcon, dval;
959     int ipix, pix2;
960 
961     pix2 = pix1 + npix;
962 
963     if (scale)
964 	dpix = (dpix - bzero) / bscale;
965     ipix = (int) dpix;
966     dcon = (double) ipix;
967     if (dcon == dpix)
968 	isint = 1;
969     else
970 	isint = 0;
971 
972     switch (bitpix) {
973 
974 	case 8:
975 	    imc = image + pix1;
976 	    if (isint) {
977 		if (dpix < 0)
978 		    ccon = (char) (dpix - 0.5);
979 		else
980 		    ccon = (char) (dpix + 0.5);
981 		for (ipix = pix1; ipix < pix2; ipix++)
982 		    *imc++ *= ccon;
983 		}
984 	    else {
985 		for (ipix = pix1; ipix < pix2; ipix++) {
986 		    dval = ((double) *imc) * dpix;
987 		    if (dval < 256.0)
988 			*imc++ = (char) dval;
989 		    else
990 			*imc++ = (char) 255;
991 		    }
992 		}
993 	    break;
994 
995 	case 16:
996 	    im2 = (short *) (image + pix1);
997 	    if (isint) {
998 		im2 = (short *)image;
999 		if (dpix < 0)
1000 		    jcon = (short) (dpix - 0.5);
1001 		else
1002 		    jcon = (short) (dpix + 0.5);
1003 		for (ipix = pix1; ipix < pix2; ipix++)
1004 		    *im2++ *= jcon;
1005 		}
1006 	    else {
1007 		for (ipix = pix1; ipix < pix2; ipix++) {
1008 		    dval = ((double) *im2) * dpix;
1009 		    if (dval < 32768.0)
1010 			*im2++ = (short) dval;
1011 		    else
1012 			*im2++ = (short) 32767;
1013 		    }
1014 		}
1015 	    break;
1016 
1017 	case 32:
1018 	    im4 = (int *) (image + pix1);
1019 	    if (isint) {
1020 		if (dpix < 0)
1021 		    icon = (int) (dpix - 0.5);
1022 		else
1023 		    icon = (int) (dpix + 0.5);
1024 		for (ipix = pix1; ipix < pix2; ipix++)
1025 		    *im4++ *= icon;
1026 		}
1027 	    else {
1028 		for (ipix = pix1; ipix < pix2; ipix++) {
1029 		    dval = ((double) *im4) * dpix;
1030 		    if (dval < 32768.0)
1031 			*im4++ = (int) dval;
1032 		    else
1033 			*im4++ = (int) 32767;
1034 		    }
1035 		}
1036 	    break;
1037 
1038 	case -16:
1039 	    imu = (unsigned short *) (image + pix1);
1040 	    if (dpix > 0) {
1041 		ucon = (unsigned short) (dpix + 0.5);
1042 		imu = (unsigned short *) (image + pix1);
1043 		for (ipix = pix1; ipix < pix2; ipix++)
1044 		    *imu++ *= ucon;
1045 		}
1046 	    break;
1047 
1048 	case -32:
1049 	    rcon = (float) dpix;
1050 	    imr = (float *) (image + pix1);
1051 	    for (ipix = pix1; ipix < pix2; ipix++)
1052 		*imr++ *= rcon;
1053 	    break;
1054 
1055 	case -64:
1056 	    imd = (double *) (image + pix1);
1057 	    for (ipix = pix1; ipix < pix2; ipix++)
1058 		*imd++ *= dpix;
1059 	    break;
1060 
1061 	}
1062     return;
1063 }
1064 
1065 
1066 /* GETVEC -- Get vector from 2D image of any numeric type */
1067 
1068 void
getvec(image,bitpix,bzero,bscale,pix1,npix,dvec0)1069 getvec (image, bitpix, bzero, bscale, pix1, npix, dvec0)
1070 
1071 char	*image;		/* Image array from which to extract vector */
1072 int	bitpix;		/* Number of bits per pixel in image */
1073 			/*  16 = short, -16 = unsigned short, 32 = int */
1074 			/* -32 = float, -64 = double */
1075 double  bzero;		/* Zero point for pixel scaling */
1076 double  bscale;		/* Scale factor for pixel scaling */
1077 int	pix1;		/* Offset of first pixel to extract */
1078 int	npix;		/* Number of pixels to extract */
1079 double	*dvec0;		/* Vector of pixels (returned) */
1080 
1081 {
1082     short *im2;
1083     int *im4;
1084     unsigned short *imu;
1085     float *imr;
1086     double *imd;
1087     double *dvec;
1088     int ipix, pix2;
1089 
1090     pix2 = pix1 + npix;
1091     dvec = dvec0;
1092 
1093     switch (bitpix) {
1094 
1095 	case 8:
1096 	    for (ipix = pix1; ipix < pix2; ipix++)
1097 		*dvec++ = (double) *(image + ipix);
1098 	    break;
1099 
1100 	case 16:
1101 	    im2 = (short *)image;
1102 	    for (ipix = pix1; ipix < pix2; ipix++)
1103 		*dvec++ = (double) *(im2 + ipix);
1104 	    break;
1105 
1106 	case 32:
1107 	    im4 = (int *)image;
1108 	    for (ipix = pix1; ipix < pix2; ipix++)
1109 		*dvec++ = (double) *(im4 + ipix);
1110 	    break;
1111 
1112 	case -16:
1113 	    imu = (unsigned short *)image;
1114 	    for (ipix = pix1; ipix < pix2; ipix++)
1115 		*dvec++ = (double) *(imu + ipix);
1116 	    break;
1117 
1118 	case -32:
1119 	    imr = (float *)image;
1120 	    for (ipix = pix1; ipix < pix2; ipix++)
1121 		*dvec++ = (double) *(imr + ipix);
1122 	    break;
1123 
1124 	case -64:
1125 	    imd = (double *)image;
1126 	    for (ipix = pix1; ipix < pix2; ipix++)
1127 		*dvec++ = (double) *(imd + ipix);
1128 	    break;
1129 
1130 	}
1131 
1132     /* Scale data if either BZERO or BSCALE keyword has been set */
1133     if (scale && (bzero != 0.0 || bscale != 1.0)) {
1134 	dvec = dvec0;
1135 	for (ipix = pix1; ipix < pix2; ipix++) {
1136 	    *dvec = (*dvec * bscale) + bzero;
1137 	    dvec++;
1138 	    }
1139 	}
1140 
1141     return;
1142 }
1143 
1144 
1145 /* PUTVEC -- Copy pixel vector into 2D image of any numeric type */
1146 
1147 void
putvec(image,bitpix,bzero,bscale,pix1,npix,dvec)1148 putvec (image, bitpix, bzero, bscale, pix1, npix, dvec)
1149 
1150 char	*image;		/* Image into which to copy vector */
1151 int	bitpix;		/* Number of bits per pixel im image */
1152 			/*  16 = short, -16 = unsigned short, 32 = int */
1153 			/* -32 = float, -64 = double */
1154 double  bzero;		/* Zero point for pixel scaling */
1155 double  bscale;		/* Scale factor for pixel scaling */
1156 int	pix1;		/* Offset of first pixel of vector in image */
1157 int	npix;		/* Number of pixels to copy */
1158 double	*dvec;		/* Vector of pixels to copy */
1159 
1160 {
1161     short *im2;
1162     int *im4;
1163     unsigned short *imu;
1164     float *imr;
1165     double *imd;
1166     int ipix, pix2;
1167     double *dp = dvec;
1168 
1169     pix2 = pix1 + npix;
1170 
1171     /* Scale data if either BZERO or BSCALE keyword has been set */
1172     if (scale && (bzero != 0.0 || bscale != 1.0)) {
1173 	for (ipix = pix1; ipix < pix2; ipix++) {
1174 	    *dp = (*dp - bzero) / bscale;
1175 	    dp++;
1176 	    }
1177 	dp = dvec;
1178 	}
1179 
1180     switch (bitpix) {
1181 
1182 	case 8:
1183 	    for (ipix = pix1; ipix < pix2; ipix++)
1184 		*(image+ipix) = (char) *dp++;
1185 	    break;
1186 
1187 	case 16:
1188 	    im2 = (short *)image;
1189 	    for (ipix = pix1; ipix < pix2; ipix++) {
1190 		if (*dp < 0.0)
1191 		    *(im2+ipix) = (short) (*dp++ - 0.5);
1192 		else
1193 		    *(im2+ipix) = (short) (*dp++ + 0.5);
1194 		}
1195 	    break;
1196 
1197 	case 32:
1198 	    im4 = (int *)image;
1199 	    for (ipix = pix1; ipix < pix2; ipix++) {
1200 		if (*dp < 0.0)
1201 		    *(im4+ipix) = (int) (*dp++ - 0.5);
1202 		else
1203 		    *(im4+ipix) = (int) (*dp++ + 0.5);
1204 		}
1205 	    break;
1206 
1207 	case -16:
1208 	    imu = (unsigned short *)image;
1209 	    for (ipix = pix1; ipix < pix2; ipix++) {
1210 		if (*dp < 0.0)
1211 		    *(imu+ipix) = (unsigned short) 0;
1212 		else
1213 		    *(imu+ipix) = (unsigned short) (*dp++ + 0.5);
1214 		}
1215 	    break;
1216 
1217 	case -32:
1218 	    imr = (float *)image;
1219 	    for (ipix = pix1; ipix < pix2; ipix++)
1220 		*(imr+ipix) = (float) *dp++;
1221 	    break;
1222 
1223 	case -64:
1224 	    imd = (double *)image;
1225 	    for (ipix = pix1; ipix < pix2; ipix++)
1226 		*(imd+ipix) = (double) *dp++;
1227 	    break;
1228 	}
1229     return;
1230 }
1231 
1232 
1233 /* FILLVEC1 -- Copy single value into a vector of any numeric type */
1234 
1235 void
fillvec1(image,bitpix,bzero,bscale,pix1,npix,dpix)1236 fillvec1 (image, bitpix, bzero, bscale, pix1, npix, dpix)
1237 
1238 char	*image;		/* Vector to fill */
1239 int	bitpix;		/* Number of bits per pixel im image */
1240 			/*  16 = short, -16 = unsigned short, 32 = int */
1241 			/* -32 = float, -64 = double */
1242 double  bzero;		/* Zero point for pixel scaling */
1243 double  bscale;		/* Scale factor for pixel scaling */
1244 int	pix1;		/* First pixel to fill */
1245 int	npix;		/* Number of pixels to fill */
1246 double	dpix;		/* Value with which to fill pixels */
1247 {
1248     fillvec (image, bitpix, bzero, bscale, pix1-1, npix, dpix);
1249     return;
1250 }
1251 
1252 
1253 /* FILLVEC -- Copy single value into a vector of any numeric type */
1254 
1255 void
fillvec(image,bitpix,bzero,bscale,pix1,npix,dpix)1256 fillvec (image, bitpix, bzero, bscale, pix1, npix, dpix)
1257 
1258 char	*image;		/* Vector to fill */
1259 int	bitpix;		/* Number of bits per pixel im image */
1260 			/*  16 = short, -16 = unsigned short, 32 = int */
1261 			/* -32 = float, -64 = double */
1262 double  bzero;		/* Zero point for pixel scaling */
1263 double  bscale;		/* Scale factor for pixel scaling */
1264 int	pix1;		/* First pixel to fill */
1265 int	npix;		/* Number of pixels to fill */
1266 double	dpix;		/* Value with which to fill pixels */
1267 {
1268     char ipc;
1269     short *im2, ip2;
1270     int *im4, ip4;
1271     unsigned short *imu, ipu;
1272     float *imr, ipr;
1273     double *imd;
1274     int ipix, pix2;
1275     double dp;
1276 
1277     pix2 = pix1 + npix;
1278 
1279     /* Scale data if either BZERO or BSCALE keyword has been set */
1280     dp = dpix;
1281     if (scale && (bzero != 0.0 || bscale != 1.0))
1282 	dp = (dp - bzero) / bscale;
1283 
1284     switch (bitpix) {
1285 
1286 	case 8:
1287 	    if (dp < 0.0)
1288 		ipc = (char) (dp - 0.5);
1289 	    else
1290 		ipc = (char) (dp + 0.5);
1291 	    for (ipix = pix1; ipix < pix2; ipix++)
1292 		image[ipix] = ipc;
1293 	    break;
1294 
1295 	case 16:
1296 	    im2 = (short *)image;
1297 	    if (dp < 0.0)
1298 		ip2 = (short) (dp - 0.5);
1299 	    else
1300 		ip2 = (short) (dp + 0.5);
1301 	    for (ipix = pix1; ipix < pix2; ipix++)
1302 		im2[ipix] = ip2;
1303 	    break;
1304 
1305 	case 32:
1306 	    im4 = (int *)image;
1307 	    if (dp < 0.0)
1308 		ip4 = (int) (dp - 0.5);
1309 	    else
1310 		ip4 = (int) (dp + 0.5);
1311 	    for (ipix = pix1; ipix < pix2; ipix++)
1312 		im4[ipix] = ip4;
1313 	    break;
1314 
1315 	case -16:
1316 	    imu = (unsigned short *)image;
1317 	    if (dp < 0.0)
1318 		ipu = (unsigned short) (dp - 0.5);
1319 	    else
1320 		ipu = (unsigned short) (dp + 0.5);
1321 	    for (ipix = pix1; ipix < pix2; ipix++)
1322 		imu[ipix] = ipu;
1323 	    break;
1324 
1325 	case -32:
1326 	    imr = (float *)image;
1327 	    ipr = (float) dp;
1328 	    for (ipix = pix1; ipix < pix2; ipix++)
1329 		imr[ipix] = ipr;
1330 	    break;
1331 
1332 	case -64:
1333 	    imd = (double *)image;
1334 	    for (ipix = pix1; ipix < pix2; ipix++)
1335 		imd[ipix] = dp;
1336 	    break;
1337 	}
1338     return;
1339 }
1340 
1341 
1342 /* IMSWAP -- Reverse bytes of any type of vector in place */
1343 
1344 void
imswap(bitpix,string,nbytes)1345 imswap (bitpix, string, nbytes)
1346 
1347 int	bitpix;		/* Number of bits per pixel */
1348 			/*  16 = short, -16 = unsigned short, 32 = int */
1349 			/* -32 = float, -64 = double */
1350 char	*string;	/* Address of starting point of bytes to swap */
1351 int	nbytes;		/* Number of bytes to swap */
1352 
1353 {
1354     switch (bitpix) {
1355 
1356 	case 8:
1357 	    break;
1358 
1359 	case 16:
1360 	    if (nbytes < 2) return;
1361 	    imswap2 (string,nbytes);
1362 	    break;
1363 
1364 	case 32:
1365 	    if (nbytes < 4) return;
1366 	    imswap4 (string,nbytes);
1367 	    break;
1368 
1369 	case -16:
1370 	    if (nbytes < 2) return;
1371 	    imswap2 (string,nbytes);
1372 	    break;
1373 
1374 	case -32:
1375 	    if (nbytes < 4) return;
1376 	    imswap4 (string,nbytes);
1377 	    break;
1378 
1379 	case -64:
1380 	    if (nbytes < 8) return;
1381 	    imswap8 (string,nbytes);
1382 	    break;
1383 
1384 	}
1385     return;
1386 }
1387 
1388 
1389 /* IMSWAP2 -- Swap bytes in string in place */
1390 
1391 void
imswap2(string,nbytes)1392 imswap2 (string,nbytes)
1393 
1394 
1395 char *string;	/* Address of starting point of bytes to swap */
1396 int nbytes;	/* Number of bytes to swap */
1397 
1398 {
1399     char *sbyte, temp, *slast;
1400 
1401     slast = string + nbytes;
1402     sbyte = string;
1403     while (sbyte < slast) {
1404 	temp = sbyte[0];
1405 	sbyte[0] = sbyte[1];
1406 	sbyte[1] = temp;
1407 	sbyte= sbyte + 2;
1408 	}
1409     return;
1410 }
1411 
1412 
1413 /* IMSWAP4 -- Reverse bytes of Integer*4 or Real*4 vector in place */
1414 
1415 void
imswap4(string,nbytes)1416 imswap4 (string,nbytes)
1417 
1418 char *string;	/* Address of Integer*4 or Real*4 vector */
1419 int nbytes;	/* Number of bytes to reverse */
1420 
1421 {
1422     char *sbyte, *slast;
1423     char temp0, temp1, temp2, temp3;
1424 
1425     slast = string + nbytes;
1426     sbyte = string;
1427     while (sbyte < slast) {
1428 	temp3 = sbyte[0];
1429 	temp2 = sbyte[1];
1430 	temp1 = sbyte[2];
1431 	temp0 = sbyte[3];
1432 	sbyte[0] = temp0;
1433 	sbyte[1] = temp1;
1434 	sbyte[2] = temp2;
1435 	sbyte[3] = temp3;
1436 	sbyte = sbyte + 4;
1437 	}
1438 
1439     return;
1440 }
1441 
1442 
1443 /* IMSWAP8 -- Reverse bytes of Real*8 vector in place */
1444 
1445 void
imswap8(string,nbytes)1446 imswap8 (string,nbytes)
1447 
1448 char *string;	/* Address of Real*8 vector */
1449 int nbytes;	/* Number of bytes to reverse */
1450 
1451 {
1452     char *sbyte, *slast;
1453     char temp[8];
1454 
1455     slast = string + nbytes;
1456     sbyte = string;
1457     while (sbyte < slast) {
1458 	temp[7] = sbyte[0];
1459 	temp[6] = sbyte[1];
1460 	temp[5] = sbyte[2];
1461 	temp[4] = sbyte[3];
1462 	temp[3] = sbyte[4];
1463 	temp[2] = sbyte[5];
1464 	temp[1] = sbyte[6];
1465 	temp[0] = sbyte[7];
1466 	sbyte[0] = temp[0];
1467 	sbyte[1] = temp[1];
1468 	sbyte[2] = temp[2];
1469 	sbyte[3] = temp[3];
1470 	sbyte[4] = temp[4];
1471 	sbyte[5] = temp[5];
1472 	sbyte[6] = temp[6];
1473 	sbyte[7] = temp[7];
1474 	sbyte = sbyte + 8;
1475 	}
1476     return;
1477 }
1478 
1479 /* IMSWAPPED -- Returns 0 if big-endian (Sun,Mac),
1480 		1 if little-endian(PC,Alpha) */
1481 
1482 int
imswapped()1483 imswapped ()
1484 
1485 {
1486     char *ctest;
1487     int itest;
1488 
1489     itest = 1;
1490     ctest = (char *)&itest;
1491     if (*ctest)
1492 	return (1);
1493     else
1494 	return (0);
1495 }
1496 
1497 /* Apr 17 1996	New file
1498  * May 22 1996	Add H so that PUTPIX and GETPIX can check coordinates
1499  * Jun 11 1996	Simplify NEWIMAGE subroutine
1500  * Jun 12 1996	Add byte-swapping subroutines
1501  *
1502  * Jul 24 1997	Add 8-bit option to subroutines
1503  *
1504  * May 27 1998	Include imio.h instead of fitshead.h
1505  * Jun 17 1998	Fix bug, changing all unsigned int's to unsigned short's
1506  *
1507  * Apr 29 1999	Add scaling to getpix, putpix, getvec, and putvec
1508  * Apr 29 1999	Fix bug in getvec in dealing with 1-byte data
1509  * Sep 14 1999	Change dp incrementing so it works on Alpha compiler
1510  * Sep 27 1999	Add interface for 1-based (FITS) image access
1511  * Sep 27 1999	Add addpix() and addpix1()
1512  * Dec 14 1999	In putpix(), addpix(), putvec(), round when output is integer
1513  *
1514  * Sep 20 2000	In getvec(), scale only if necessary
1515  *
1516  * Nov 27 2001	In movepix(), add char to char move
1517  *
1518  * Jan 23 2002	Add global scale switch to turn off scaling
1519  * Jun  4 2002	In getvec() and putvec(), change dpix to dvec
1520  * Jun  4 2002	Add addvec() to add to a vector
1521  * Jul 19 2002	Fix getvec() bug rescaling scaled numbers
1522  *
1523  * May 20 2003	Declare scale0 in setscale()
1524  *
1525  * Jan 28 2004	Add image limit check to movepix()
1526  * Feb 27 2004	Add fillvec() and fillvec1() to set vector to a constant
1527  *
1528  * Jun 27 2005	Fix major bug in fillvec(); pass value dpix in fillvec1(), too
1529  * Aug 18 2005	Add maxvec(), addvec(), and multvec()
1530  *
1531  * Mar  1 2006	Fix bug of occasional double application of bscale in getvec()
1532  * Apr  3 2006	Fix bad cast in unisigned int section of addvec()
1533  * May  3 2006	Code fixes in addpix and multpix suggested by Robert Lupton
1534  * Jun  8 2006	Drop erroneous second im2 assignment without offset in addvec()
1535  * Jun 20 2006	Fix typos masquerading as unitialized variables
1536  *
1537  * Jan  8 2007	Include fitsfile.h instead of imio.h
1538  * Jun 11 2007	Add minvec() and speed up maxvec()
1539  *
1540  * Apr 12 2012	Fix 8-bit variables to be unsigned char
1541  * Oct 19 2012	Fix errors with character images in minvec() and maxvec()
1542  * Oct 31 2012	Fix errors with short images in minvec() and maxvec()
1543  * Oct 31 2012	Drop unused variable il2 from minvec()
1544  */
1545