?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
moved kdeaccessibility kdeaddons kdeadmin kdeartwork kdebindings kdeedu kdegames...
[kdeedu.git] / kstars / kstars / indi / fitsrw.c
blobda6084f5e5c6f92d067a97163a402e2b5c310aa0
1 /******************************************************************************/
2 /* Peter Kirchgessner */
3 /* e-mail: pkirchg@aol.com */
4 /* WWW : http://members.aol.com/pkirchg */
5 /******************************************************************************/
6 /* #BEG-HDR */
7 /* */
8 /* Package : FITS reading/writing library */
9 /* Modul-Name : fitsrw.c */
10 /* Description : Support of reading/writing FITS-files */
11 /* Function(s) : fits_new_filestruct - (local) initialize file structure*/
12 /* fits_new_hdulist - (local) initialize hdulist struct*/
13 /* fits_delete_filestruct - (local) delete file structure */
14 /* fits_delete_recordlist - (local) delete record list */
15 /* fits_delete_hdulist - (local) delete hdu list */
16 /* fits_nan_32 - (local) check IEEE NaN values */
17 /* fits_nan_64 - (local) check IEEE NaN values */
18 /* fits_get_error - get error message */
19 /* fits_set_error - (local) set error message */
20 /* fits_drop_error - (local) remove an error message */
21 /* fits_open - open a FITS file */
22 /* fits_close - close a FITS file */
23 /* fits_add_hdu - add a HDU to a FITS file */
24 /* fits_add_card - add a card to the HDU */
25 /* fits_print_header - print a single FITS header */
26 /* fits_read_header - (local) read in FITS header */
27 /* fits_write_header - write a FITS header */
28 /* fits_decode_header - (local) decode a header */
29 /* fits_eval_pixrange - (local) evaluate range of pixels */
30 /* fits_decode_card - decode a card */
31 /* fits_search_card - search a card in a record list */
32 /* fits_image_info - get information about image */
33 /* fits_seek_image - position to an image */
34 /* fits_read_pixel - read pixel values from file */
35 /* fits_to_pgmraw - convert FITS-file to PGM-file */
36 /* pgmraw_to_fits - convert PGM-file to FITS-file */
37 /* */
38 /* Author : P. Kirchgessner */
39 /* Date of Gen. : 12-Apr-97 */
40 /* Last modified : 20-Dec-97 */
41 /* Version : 0.11 */
42 /* Compiler Opt. : */
43 /* Changes : */
44 /* #MOD-0001, nn, 20-Dec-97, Initialize some variables */
45 /* */
46 /* #END-HDR */
47 /******************************************************************************/
48 /* References: */
49 /* - NOST, Definition of the Flexible Image Transport System (FITS), */
50 /* September 29, 1995, Standard, NOST 100-1.1 */
51 /* (ftp://nssdc.gsfc.nasa.gov/pub/fits/fits_standard_ps.Z) */
52 /* - The FITS IMAGE Extension. A Proposal. J.D. Ponz, R.W. Thompson, */
53 /* J.R. Munoz, Feb. 7, 1992 */
54 /* (ftp://www.cv.nrao.edu/fits/documents/standards/image.ps.gz) */
55 /* */
56 /******************************************************************************/
57
58 #define VERSIO 0.11
59 /* Identifikation: "@(#) <product> <ver> <dd-mmm-yy>" */
60 static char ident[] = "@(#) libfits.c 0.11 20-Dec-97 (%I%)";
61
62 /******************************************************************************/
63 /* FITS reading/writing library */
64 /* Copyright (C) 1997 Peter Kirchgessner */
65 /* (email: pkirchg@aol.com, WWW: http://members.aol.com/pkirchg) */
66 /* The library was developed for a FITS-plug-in to GIMP, the GNU Image */
67 /* Manipulation Program. But it is completely independant to that. If someone */
68 /* finds it useful for other purposes, try to keep it independant from your */
69 /* application. */
70 /******************************************************************************/
71 /* */
72 /* This program is free software; you can redistribute it and/or modify */
73 /* it under the terms of the GNU General Public License as published by */
74 /* the Free Software Foundation; either version 2 of the License, or */
75 /* (at your option) any later version. */
76 /* */
77 /* This program is distributed in the hope that it will be useful, */
78 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
79 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
80 /* GNU General Public License for more details. */
81 /* */
82 /* You should have received a copy of the GNU General Public License */
83 /* along with this program; if not, write to the Free Software */
84 /* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
85 /******************************************************************************/
86
87 #include <stdio.h>
88 #include <stdlib.h>
89 #include <string.h>
90
91 #include "fitsrw.h"
92
93
94 /* Declaration of local funtions */
95
96 static FITS_FILE *fits_new_filestruct (void);
97 static FITS_HDU_LIST *fits_new_hdulist (void);
98 static void fits_delete_filestruct (FITS_FILE *ff);
99 static void fits_delete_recordlist (FITS_RECORD_LIST *rl);
100 static void fits_delete_hdulist (FITS_HDU_LIST *hl);
101 int fits_nan_32 (unsigned char *value);
102 int fits_nan_64 (unsigned char *value);
103 static void fits_set_error (const char *errmsg);
104 static void fits_drop_error (void);
105 static FITS_RECORD_LIST *fits_read_header (FILE *fp, int *nrec);
106 static FITS_HDU_LIST *fits_decode_header (FITS_RECORD_LIST *hdr,
107 long hdr_offset, long dat_offset);
108 static int fits_eval_pixrange (FILE *fp, FITS_HDU_LIST *hdu);
109
110
111 /* Error handling like a FIFO */
112 #define FITS_MAX_ERROR 16
113 #define FITS_ERROR_LENGTH 256
114 static int fits_n_error = 0;
115 static char fits_error[FITS_MAX_ERROR][FITS_ERROR_LENGTH];
116
117 /* What byte ordering for IEEE-format we are running on ? */
118 int fits_ieee32_intel = 0;
119 int fits_ieee32_motorola = 0;
120 int fits_ieee64_intel = 0;
121 int fits_ieee64_motorola = 0;
122
123 /* Macros */
124 #define FITS_RETURN(msg, val) { fits_set_error (msg); return (val); }
125 #define FITS_VRETURN(msg) { fits_set_error (msg); return; }
126
127 /* Get pixel values from memory. p must be an (unsigned char *) */
128 #define FITS_GETBITPIX16(p,val) val = ((p[0] << 8) | (p[1]))
129 #define FITS_GETBITPIX32(p,val) val = \
130 ((p[0] << 24) | (p[1] << 16) | (p[2] << 8) | p[3])
131
132 /* Get floating point values from memory. p must be an (unsigned char *). */
133 /* The floating point values must directly correspond */
134 /* to machine representation. Otherwise it does not work. */
135 #define FITS_GETBITPIXM32(p,val) \
136 { if (fits_ieee32_intel) {unsigned char uc[4]; \
137 uc[0] = p[3]; uc[1] = p[2]; uc[2] = p[1]; uc[3] = p[0]; \
138 val = *(FITS_BITPIXM32 *)uc; } \
139 else if (fits_ieee32_motorola) { val = *(FITS_BITPIXM32 *)p; } \
140 else if (fits_ieee64_motorola) {FITS_BITPIXM64 m64; \
141 unsigned char *uc= (unsigned char *)&m64; \
142 uc[0]=p[0]; uc[1]=p[1]; uc[2]=p[2]; uc[3]=p[3]; uc[4]=uc[5]=uc[6]=uc[7]=0; \
143 val = (FITS_BITPIXM32)m64; } \
144 else if (fits_ieee64_intel) {FITS_BITPIXM64 i64; \
145 unsigned char *uc= (unsigned char *)&i64; \
146 uc[0]=uc[1]=uc[2]=uc[3]=0; uc[7]=p[0]; uc[6]=p[1]; uc[5]=p[2]; uc[4]=p[3]; \
147 val = (FITS_BITPIXM32)i64;}\
148 }
149
150 #define FITS_GETBITPIXM64(p,val) \
151 { if (fits_ieee64_intel) {unsigned char uc[8]; \
152 uc[0] = p[7]; uc[1] = p[6]; uc[2] = p[5]; uc[3] = p[4]; \
153 uc[4] = p[3]; uc[5] = p[2]; uc[6] = p[1]; uc[7] = p[0]; \
154 val = *(FITS_BITPIXM64 *)uc; } else val = *(FITS_BITPIXM64 *)p; }
155
156 #define FITS_WRITE_BOOLCARD(fp,key,value) \
157 {char card[81]; \
158 snprintf (card, sizeof( card ), "%-8.8s= %20s%50s", key, value ? "T" : "F", " "); \
159 fwrite (card, 1, 80, fp); }
160
161 #define FITS_WRITE_LONGCARD(fp,key,value) \
162 {char card[81]; \
163 snprintf (card, sizeof( card ), "%-8.8s= %20ld%50s", key, (long)value, " "); \
164 fwrite (card, 1, 80, fp); }
165
166 #define FITS_WRITE_DOUBLECARD(fp,key,value) \
167 {char card[81], dbl[21], *istr; \
168 snprintf (dbl, sizeof( dbl ), "%20f", (double)value); istr = strstr (dbl, "e"); \
169 if (istr) *istr = 'E'; \
170 snprintf (card, sizeof( card ), "%-8.8s= %20.20s%50s", key, dbl, " "); \
171 fwrite (card, 1, 80, fp); }
172
173 #define FITS_WRITE_STRINGCARD(fp,key,value) \
174 {char card[81]; int k;\
175 snprintf (card, sizeof( card ), "%-8.8s= \'%s", key, value); \
176 for (k = strlen (card); k < 81; k++) card[k] = ' '; \
177 k = strlen (key); if (k < 8) card[19] = '\''; else card[11+k] = '\''; \
178 fwrite (card, 1, 80, fp); }
179
180 #define FITS_WRITE_CARD(fp,value) \
181 {char card[81]; \
182 snprintf (card, sizeof( card ), "%-80.80s", value); \
183 fwrite (card, 1, 80, fp); }
184
185
186 /*****************************************************************************/
187 /* #BEG-PAR */
188 /* */
189 /* Function : fits_new_filestruct - (local) initialize file structure */
190 /* */
191 /* Parameters: */
192 /* -none- */
193 /* */
194 /* Returns a pointer to an initialized fits file structure. */
195 /* On failure, a NULL-pointer is returned. */
196 /* */
197 /* #END-PAR */
198 /*****************************************************************************/
199
200 static FITS_FILE *fits_new_filestruct (void)
201
202 {FITS_FILE *ff;
203
204 ff = (FITS_FILE *)malloc (sizeof (FITS_FILE));
205 if (ff == NULL) return (NULL);
206
207 memset ((char *)ff, 0, sizeof (*ff));
208 return (ff);
209 }
210
211
212 /*****************************************************************************/
213 /* #BEG-PAR */
214 /* */
215 /* Function : fits_new_hdulist - (local) initialize hdulist structure */
216 /* */
217 /* Parameters: */
218 /* -none- */
219 /* */
220 /* Returns a pointer to an initialized hdulist structure. */
221 /* On failure, a NULL-pointer is returned. */
222 /* */
223 /* #END-PAR */
224 /*****************************************************************************/
225
226 static FITS_HDU_LIST *fits_new_hdulist (void)
227
228 {FITS_HDU_LIST *hdl;
229
230 hdl = (FITS_HDU_LIST *)malloc (sizeof (FITS_HDU_LIST));
231 if (hdl == NULL) return (NULL);
232
233 memset ((char *)hdl, 0, sizeof (*hdl));
234 hdl->pixmin = hdl->pixmax = hdl->datamin = hdl->datamax = 0.0;
235 hdl->bzero = hdl->bscale = 0.0;
236
237 return (hdl);
238 }
239
240
241 /*****************************************************************************/
242 /* #BEG-PAR */
243 /* */
244 /* Function : fits_delete_filestruct - (local) delete file structure */
245 /* */
246 /* Parameters: */
247 /* FITS_FILE *ff [I] : pointer to fits file structure */
248 /* ( mode : I=input, O=output, I/O=input/output ) */
249 /* */
250 /* Frees all memory allocated by the file structure. */
251 /* */
252 /* #END-PAR */
253 /*****************************************************************************/
254
255 static void fits_delete_filestruct (FITS_FILE *ff)
256
257 {
258 if (ff == NULL) return;
259
260 fits_delete_hdulist (ff->hdu_list);
261 ff->hdu_list = NULL;
262
263 ff->fp = NULL;
264 free ((char *)ff);
265 }
266
267
268 /*****************************************************************************/
269 /* #BEG-PAR */
270 /* */
271 /* Function : fits_delete_recordlist - (local) delete record list */
272 /* */
273 /* Parameters: */
274 /* FITS_RECORD_LIST *rl [I] : record list to delete */
275 /* ( mode : I=input, O=output, I/O=input/output ) */
276 /* */
277 /* #END-PAR */
278 /*****************************************************************************/
279
280 static void fits_delete_recordlist (FITS_RECORD_LIST *rl)
281
282 {FITS_RECORD_LIST *next;
283
284 while (rl != NULL)
285 {
286 next = rl->next_record;
287 rl->next_record = NULL;
288 free ((char *)rl);
289 rl = next;
290 }
291 }
292
293
294 /*****************************************************************************/
295 /* #BEG-PAR */
296 /* */
297 /* Function : fits_delete_hdulist - (local) delete hdu list */
298 /* */
299 /* Parameters: */
300 /* FITS_HDU_LIST *hl [I] : hdu list to delete */
301 /* ( mode : I=input, O=output, I/O=input/output ) */
302 /* */
303 /* #END-PAR */
304 /*****************************************************************************/
305
306 static void fits_delete_hdulist (FITS_HDU_LIST *hl)
307
308 {FITS_HDU_LIST *next;
309
310 while (hl != NULL)
311 {
312 fits_delete_recordlist (hl->header_record_list);
313 next = hl->next_hdu;
314 hl->next_hdu = NULL;
315 free ((char *)hl);
316 hl = next;
317 }
318 }
319
320
321 /*****************************************************************************/
322 /* #BEG-PAR */
323 /* */
324 /* Function : fits_nan_32 - (local) check for IEEE NaN values (32 bit) */
325 /* */
326 /* Parameters: */
327 /* unsigned char *v [I] : value to check */
328 /* ( mode : I=input, O=output, I/O=input/output ) */
329 /* */
330 /* The function returns 1 if the value is a NaN. Otherwise 0 is returned. */
331 /* The byte sequence at v must start with the sign/eponent byte. */
332 /* */
333 /* #END-PAR */
334 /*****************************************************************************/
335
336 int fits_nan_32 (unsigned char *v)
337
338 {register unsigned long k;
339
340 k = (v[0] << 24) | (v[1] << 16) | (v[2] << 8) | v[3];
341 k &= 0x7fffffff; /* Dont care about the sign bit */
342
343 /* See NOST Definition of the Flexible Image Transport System (FITS), */
344 /* Appendix F, IEEE special formats. */
345 return ( ((k >= 0x7f7fffff) && (k <= 0x7fffffff))
346 || ((k >= 0x00000001) && (k <= 0x00800000)));
347 }
348
349
350 /*****************************************************************************/
351 /* #BEG-PAR */
352 /* */
353 /* Function : fits_nan_64 - (local) check for IEEE NaN values (64 bit) */
354 /* */
355 /* Parameters: */
356 /* unsigned char *v [I] : value to check */
357 /* ( mode : I=input, O=output, I/O=input/output ) */
358 /* */
359 /* The function returns 1 if the value is a NaN. Otherwise 0 is returned. */
360 /* The byte sequence at v must start with the sign/eponent byte. */
361 /* (currently we ignore the low order 4 bytes of the mantissa. Therefore */
362 /* this function is the same as for 32 bits). */
363 /* */
364 /* #END-PAR */
365 /*****************************************************************************/
366
367 int fits_nan_64 (unsigned char *v)
368
369 {register unsigned long k;
370
371 k = (v[0] << 24) | (v[1] << 16) | (v[2] << 8) | v[3];
372 k &= 0x7fffffff; /* Dont care about the sign bit */
373
374 /* See NOST Definition of the Flexible Image Transport System (FITS), */
375 /* Appendix F, IEEE special formats. */
376 return ( ((k >= 0x7f7fffff) && (k <= 0x7fffffff))
377 || ((k >= 0x00000001) && (k <= 0x00800000)));
378 }
379
380
381 /*****************************************************************************/
382 /* #BEG-PAR */
383 /* */
384 /* Function : fits_get_error - get an error message */
385 /* */
386 /* Parameters: */
387 /* -none- */
388 /* */
389 /* If an error message has been set, a pointer to the message is returned. */
390 /* Otherwise a NULL pointer is returned. */
391 /* An inquired error message is removed from the error FIFO. */
392 /* */
393 /* #END-PAR */
394 /*****************************************************************************/
395
396 char *fits_get_error (void)
397
398 {static char errmsg[FITS_ERROR_LENGTH];
399 int k;
400
401 if (fits_n_error <= 0) return (NULL);
402 strcpy (errmsg, fits_error[0]);
403
404 for (k = 1; k < fits_n_error; k++)
405 strcpy (fits_error[k-1], fits_error[k]);
406
407 fits_n_error--;
408
409 return (errmsg);
410 }
411
412
413 /*****************************************************************************/
414 /* #BEG-PAR */
415 /* */
416 /* Function : fits_set_error - (local) set an error message */
417 /* */
418 /* Parameters: */
419 /* char *errmsg [I] : Error message to set */
420 /* ( mode : I=input, O=output, I/O=input/output ) */
421 /* */
422 /* Places the error message in the FIFO. If the FIFO is full, */
423 /* the message is discarded. */
424 /* */
425 /* #END-PAR */
426 /*****************************************************************************/
427
428 static void fits_set_error (const char *errmsg)
429
430 {
431 if (fits_n_error < FITS_MAX_ERROR)
432 {
433 strncpy (fits_error[fits_n_error], errmsg, FITS_ERROR_LENGTH);
434 fits_error[fits_n_error++][FITS_ERROR_LENGTH-1] = '\0';
435 }
436 }
437
438
439 /*****************************************************************************/
440 /* #BEG-PAR */
441 /* */
442 /* Function : fits_drop_error - (local) remove an error message */
443 /* */
444 /* Parameters: */
445 /* -none- */
446 /* */
447 /* Removes the last error message from the error message FIFO */
448 /* */
449 /* #END-PAR */
450 /*****************************************************************************/
451
452 static void fits_drop_error (void)
453
454 {
455 if (fits_n_error > 0) fits_n_error--;
456 }
457
458
459 /*****************************************************************************/
460 /* #BEG-PAR */
461 /* */
462 /* Function : fits_open - open a FITS file */
463 /* */
464 /* Parameters: */
465 /* char *filename [I] : name of file to open */
466 /* char *openmode [I] : mode to open the file ("r", "w") */
467 /* ( mode : I=input, O=output, I/O=input/output ) */
468 /* */
469 /* On success, a FITS_FILE-pointer is returned. On failure, a NULL- */
470 /* pointer is returned. */
471 /* The functions scans through the file loading each header and analyzing */
472 /* them. */
473 /* */
474 /* #END-PAR */
475 /*****************************************************************************/
476
477
478 FITS_FILE *fits_open (const char* filename, const char *openmode)
479
480 {int reading, writing, n_rec, n_hdr;
481 long fpos_header, fpos_data;
482 FILE *fp;
483 FITS_FILE *ff;
484 FITS_RECORD_LIST *hdrlist;
485 FITS_HDU_LIST *hdulist, *last_hdulist;
486
487 if ((filename == NULL) || (*filename == '\0') || (openmode == NULL))
488 FITS_RETURN ("fits_open: Invalid parameters", NULL);
489
490 /* initialize */
491 hdulist = NULL;
492 last_hdulist = NULL;
493
494 /* Check the IEEE-format we are running on */
495 {float one32 = 1.0;
496 double one64 = 1.0;
497 unsigned char *op32 = (unsigned char *)&one32;
498 unsigned char *op64 = (unsigned char *)&one64;
499
500 if (sizeof (float) == 4)
501 {
502 fits_ieee32_intel = (op32[3] == 0x3f);
503 fits_ieee32_motorola = (op32[0] == 0x3f);
504 }
505 if (sizeof (double) == 8)
506 {
507 fits_ieee64_intel = (op64[7] == 0x3f);
508 fits_ieee64_motorola = (op64[0] == 0x3f);
509 }
510 }
511
512 reading = (strcmp (openmode, "r") == 0);
513 writing = (strcmp (openmode, "w") == 0);
514 if ((!reading) && (!writing))
515 FITS_RETURN ("fits_open: Invalid openmode", NULL);
516
517 fp = fopen (filename, reading ? "rb" : "wb");
518 if (fp == NULL) FITS_RETURN ("fits_open: fopen() failed", NULL);
519
520 ff = fits_new_filestruct ();
521 if (ff == NULL)
522 {
523 fclose (fp);
524 FITS_RETURN ("fits_open: No more memory", NULL);
525 }
526
527 ff->fp = fp;
528 ff->openmode = *openmode;
529
530 if (writing) return (ff);
531
532 for (n_hdr = 0; ; n_hdr++) /* Read through all HDUs */
533 {
534 fpos_header = ftell (fp); /* Save file position of header */
535 hdrlist = fits_read_header (fp, &n_rec);
536
537 if (hdrlist == NULL)
538 {
539 if (n_hdr > 0) /* At least one header must be present. */
540 fits_drop_error (); /* If we got a header already, drop the error */
541 break;
542 }
543 fpos_data = ftell (fp); /* Save file position of data */
544
545 /* Decode the header */
546 hdulist = fits_decode_header (hdrlist, fpos_header, fpos_data);
547 if (hdulist == NULL)
548 {
549 fits_delete_recordlist (hdrlist);
550 break;
551 }
552 ff->n_hdu++;
553 ff->n_pic += hdulist->numpic;
554
555 if (hdulist->used.blank_value) ff->blank_used = 1;
556 if (hdulist->used.nan_value) ff->nan_used = 1;
557
558 if (n_hdr == 0)
559 ff->hdu_list = hdulist;
560 else
561 last_hdulist->next_hdu = hdulist;
562 last_hdulist = hdulist;
563 /* Evaluate the range of pixel data */
564 fits_eval_pixrange (fp, hdulist);
565
566 /* Reposition to start of next header */
567 if (fseek (fp, hdulist->data_offset+hdulist->data_size, SEEK_SET) < 0)
568 break;
569 }
570
571 return (ff);
572 }
573
574
575 /*****************************************************************************/
576 /* #BEG-PAR */
577 /* */
578 /* Function : fits_close - close a FITS file */
579 /* */
580 /* Parameters: */
581 /* FITS_FILE *ff [I] : FITS file pointer */
582 /* ( mode : I=input, O=output, I/O=input/output ) */
583 /* */
584 /* #END-PAR */
585 /*****************************************************************************/
586
587 void fits_close (FITS_FILE *ff)
588
589 {
590 if (ff == NULL) FITS_VRETURN ("fits_close: Invalid parameter");
591
592 fclose (ff->fp);
593
594 fits_delete_filestruct (ff);
595 }
596
597
598 /*****************************************************************************/
599 /* #BEG-PAR */
600 /* */
601 /* Function : fits_add_hdu - add a HDU to the file */
602 /* */
603 /* Parameters: */
604 /* FITS_FILE *ff [I] : FITS file pointer */
605 /* ( mode : I=input, O=output, I/O=input/output ) */
606 /* */
607 /* Adds a new HDU to the list kept in ff. A pointer to the new HDU is */
608 /* returned. On failure, a NULL-pointer is returned. */
609 /* */
610 /* #END-PAR */
611 /*****************************************************************************/
612
613 FITS_HDU_LIST *fits_add_hdu (FITS_FILE *ff)
614
615 {FITS_HDU_LIST *newhdu, *hdu;
616
617 if (ff->openmode != 'w')
618 FITS_RETURN ("fits_add_hdu: file not open for writing", NULL);
619
620 newhdu = fits_new_hdulist ();
621 if (newhdu == NULL) return (NULL);
622
623 if (ff->hdu_list == NULL)
624 {
625 ff->hdu_list = newhdu;
626 }
627 else
628 {
629 hdu = ff->hdu_list;
630 while (hdu->next_hdu != NULL)
631 hdu = hdu->next_hdu;
632 hdu->next_hdu = newhdu;
633 }
634
635 return (newhdu);
636 }
637
638
639 /*****************************************************************************/
640 /* #BEG-PAR */
641 /* */
642 /* Function : fits_add_card - add a card to the HDU */
643 /* */
644 /* Parameters: */
645 /* FITS_HDU_LIST *hdulist [I] : HDU listr */
646 /* char *card [I] : card to add */
647 /* ( mode : I=input, O=output, I/O=input/output ) */
648 /* */
649 /* The card must follow the standards of FITS. The card must not use a */
650 /* keyword that is written using *hdulist itself. On success 0 is returned. */
651 /* On failure -1 is returned. */
652 /* */
653 /* #END-PAR */
654 /*****************************************************************************/
655
656 int fits_add_card (FITS_HDU_LIST *hdulist, char *card)
657
658 {int k;
659
660 if (hdulist->naddcards >= FITS_NADD_CARDS) return (-1);
661
662 k = strlen (card);
663 if (k < FITS_CARD_SIZE)
664 {
665 memset (&(hdulist->addcards[hdulist->naddcards][k]), ' ', FITS_CARD_SIZE-k);
666 memcpy (hdulist->addcards[(hdulist->naddcards)++], card, k);
667 }
668 else
669 {
670 memcpy (hdulist->addcards[(hdulist->naddcards)++], card, FITS_CARD_SIZE);
671 }
672 return (0);
673 }
674
675
676 /*****************************************************************************/
677 /* #BEG-PAR */
678 /* */
679 /* Function : fits_print_header - print the internal representation */
680 /* of a single header */
681 /* Parameters: */
682 /* FITS_HDU_LIST *hdr [I] : pointer to the header */
683 /* ( mode : I=input, O=output, I/O=input/output ) */
684 /* */
685 /* #END-PAR */
686 /*****************************************************************************/
687
688 void fits_print_header (FITS_HDU_LIST *hdr)
689
690 {int k;
691
692 if (hdr->used.simple)
693 printf ("Content of SIMPLE-header:\n");
694 else
695 printf ("Content of XTENSION-header %s:\n", hdr->xtension);
696 printf ("header_offset : %ld\n", hdr->header_offset);
697 printf ("data_offset : %ld\n", hdr->data_offset);
698 printf ("data_size : %ld\n", hdr->data_size);
699 printf ("used data_size: %ld\n", hdr->udata_size);
700 printf ("bytes p.pixel : %d\n", hdr->bpp);
701 printf ("pixmin : %f\n", hdr->pixmin);
702 printf ("pixmax : %f\n", hdr->pixmax);
703
704 printf ("naxis : %d\n", hdr->naxis);
705 for (k = 1; k <= hdr->naxis; k++)
706 printf ("naxis%-3d : %d\n", k, hdr->naxisn[k-1]);
707
708 printf ("bitpix : %d\n", hdr->bitpix);
709
710 if (hdr->used.blank)
711 printf ("blank : %ld\n", hdr->blank);
712 else
713 printf ("blank : not used\n");
714
715 if (hdr->used.datamin)
716 printf ("datamin : %f\n", hdr->datamin);
717 else
718 printf ("datamin : not used\n");
719 if (hdr->used.datamax)
720 printf ("datamax : %f\n", hdr->datamax);
721 else
722 printf ("datamax : not used\n");
723
724 if (hdr->used.gcount)
725 printf ("gcount : %ld\n", hdr->gcount);
726 else
727 printf ("gcount : not used\n");
728 if (hdr->used.pcount)
729 printf ("pcount : %ld\n", hdr->pcount);
730 else
731 printf ("pcount : not used\n");
732
733 if (hdr->used.bscale)
734 printf ("bscale : %f\n", hdr->bscale);
735 else
736 printf ("bscale : not used\n");
737 if (hdr->used.bzero)
738 printf ("bzero : %f\n", hdr->bzero);
739 else
740 printf ("bzero : not used\n");
741 }
742
743
744 /*****************************************************************************/
745 /* #BEG-PAR */
746 /* */
747 /* Function : fits_read_header - (local) read FITS header */
748 /* */
749 /* Parameters: */
750 /* FILE *fp [I] : file pointer */
751 /* int *nrec [O] : number of records read */
752 /* ( mode : I=input, O=output, I/O=input/output ) */
753 /* */
754 /* Reads in all header records up to the record that keeps the END-card. */
755 /* A pointer to the record list is returned on success. */
756 /* On failure, a NULL-pointer is returned. */
757 /* */
758 /* #END-PAR */
759 /*****************************************************************************/
760
761 static FITS_RECORD_LIST *fits_read_header (FILE *fp, int *nrec)
762
763 {unsigned char record[FITS_RECORD_SIZE];
764 FITS_RECORD_LIST *start_list = NULL, *cu_record = NULL, *new_record;
765 FITS_DATA *fdat;
766 int k, simple, xtension;
767
768 *nrec = 0;
769
770 k = fread (record, 1, FITS_RECORD_SIZE, fp);
771 if (k != FITS_RECORD_SIZE)
772 FITS_RETURN ("fits_read_header: Error in read of first record", NULL);
773
774 simple = (strncmp (record, "SIMPLE ", 8) == 0);
775 xtension = (strncmp (record, "XTENSION", 8) == 0);
776 if ((!simple) && (!xtension))
777 FITS_RETURN ("fits_read_header: Missing keyword SIMPLE or XTENSION", NULL);
778
779 if (simple)
780 {
781 fdat = fits_decode_card (record, typ_fbool);
782 if (fdat && !fdat->fbool)
783 fits_set_error ("fits_read_header (warning): keyword SIMPLE does not have\
784 value T");
785 }
786
787 for (;;) /* Process all header records */
788 {
789 new_record = (FITS_RECORD_LIST *)malloc (sizeof (FITS_RECORD_LIST));
790 if (new_record == NULL)
791 {
792 fits_delete_recordlist (start_list);
793 FITS_RETURN ("fits_read_header: Not enough memory", NULL);
794 }
795 memcpy (new_record->data, record, FITS_RECORD_SIZE);
796 new_record->next_record = NULL;
797 (*nrec)++;
798
799 if (start_list == NULL) /* Add new record to the list */
800 start_list = new_record;
801 else
802 cu_record->next_record = new_record;
803
804 cu_record = new_record;
805 /* Was this the last record ? */
806 if (fits_search_card (cu_record, "END") != NULL) break;
807
808 k = fread (record, 1, FITS_RECORD_SIZE, fp);
809 if (k != FITS_RECORD_SIZE)
810 FITS_RETURN ("fits_read_header: Error in read of record", NULL);
811 }
812 return (start_list);
813 }
814
815
816 /*****************************************************************************/
817 /* #BEG-PAR */
818 /* */
819 /* Function : fits_write_header - write a FITS header */
820 /* */
821 /* Parameters: */
822 /* FITS_FILE *ff [I] : FITS-file pointer */
823 /* FITS_HDU_LIST [I] : pointer to header */
824 /* ( mode : I=input, O=output, I/O=input/output ) */
825 /* */
826 /* Writes a header to the file. On success, 0 is returned. On failure, */
827 /* -1 is returned. */
828 /* */
829 /* #END-PAR */
830 /*****************************************************************************/
831
832 int fits_write_header (FITS_FILE *ff, FITS_HDU_LIST *hdulist)
833
834 {int numcards;
835 int r;
836
837 if (ff->openmode != 'w')
838 FITS_RETURN ("fits_write_header: file not open for writing", -1);
839
840 numcards = 0;
841
842 if (hdulist->used.simple)
843 {
844 FITS_WRITE_BOOLCARD (ff->fp, "SIMPLE", 1);
845 numcards++;
846 }
847 else if (hdulist->used.xtension)
848 {
849 FITS_WRITE_STRINGCARD (ff->fp, "XTENSION", hdulist->xtension);
850 numcards++;
851 }
852
853 FITS_WRITE_LONGCARD (ff->fp, "BITPIX", hdulist->bitpix);
854 numcards++;
855
856 FITS_WRITE_LONGCARD (ff->fp, "NAXIS", hdulist->naxis);
857 numcards++;
858
859 for (r = 0; r < hdulist->naxis; r++)
860 {char naxisn[10];
861 snprintf (naxisn, sizeof( naxisn ), "NAXIS%d", r+1);
862 FITS_WRITE_LONGCARD (ff->fp, naxisn, hdulist->naxisn[r]);
863 numcards++;
864 }
865
866 if (hdulist->used.extend)
867 {
868 FITS_WRITE_BOOLCARD (ff->fp, "EXTEND", hdulist->extend);
869 numcards++;
870 }
871
872 if (hdulist->used.groups)
873 {
874 FITS_WRITE_BOOLCARD (ff->fp, "GROUPS", hdulist->groups);
875 numcards++;
876 }
877
878 if (hdulist->used.pcount)
879 {
880 FITS_WRITE_LONGCARD (ff->fp, "PCOUNT", hdulist->pcount);
881 numcards++;
882 }
883 if (hdulist->used.gcount)
884 {
885 FITS_WRITE_LONGCARD (ff->fp, "GCOUNT", hdulist->gcount);
886 numcards++;
887 }
888
889 if (hdulist->used.bzero)
890 {
891 FITS_WRITE_DOUBLECARD (ff->fp, "BZERO", hdulist->bzero);
892 numcards++;
893 }
894 if (hdulist->used.bscale)
895 {
896 FITS_WRITE_DOUBLECARD (ff->fp, "BSCALE", hdulist->bscale);
897 numcards++;
898 }
899
900 if (hdulist->used.datamin)
901 {
902 FITS_WRITE_DOUBLECARD (ff->fp, "DATAMIN", hdulist->datamin);
903 numcards++;
904 }
905 if (hdulist->used.datamax)
906 {
907 FITS_WRITE_DOUBLECARD (ff->fp, "DATAMAX", hdulist->datamax);
908 numcards++;
909 }
910
911 if (hdulist->used.blank)
912 {
913 FITS_WRITE_LONGCARD (ff->fp, "BLANK", hdulist->blank);
914 numcards++;
915 }
916
917 /* Write additional cards */
918 if (hdulist->naddcards > 0)
919 {
920 fwrite (hdulist->addcards, FITS_CARD_SIZE, hdulist->naddcards, ff->fp);
921 numcards += hdulist->naddcards;
922 }
923
924 FITS_WRITE_CARD (ff->fp, "END");
925 numcards++;
926
927 r = (numcards*FITS_CARD_SIZE) % FITS_RECORD_SIZE;
928 if (r) /* Must the record be filled up ? */
929 {
930 while (r++ < FITS_RECORD_SIZE)
931 putc (' ', ff->fp);
932 }
933
934 return (ferror (ff->fp) ? -1 : 0);
935 }
936
937
938 /*****************************************************************************/
939 /* #BEG-PAR */
940 /* */
941 /* Function : fits_decode_header - (local) decode a header */
942 /* */
943 /* Parameters: */
944 /* FITS_RECORD_LIST *hdr [I] : the header record list */
945 /* long hdr_offset [I] : fileposition of header */
946 /* long dat_offset [I] : fileposition of data (end of header) */
947 /* ( mode : I=input, O=output, I/O=input/output ) */
948 /* */
949 /* The function decodes the mostly used data within the header and generates */
950 /* a FITS_HDU_LIST-entry. On failure, a NULL-pointer is returned. */
951 /* */
952 /* #END-PAR */
953 /*****************************************************************************/
954
955 static FITS_HDU_LIST *fits_decode_header (FITS_RECORD_LIST *hdr,
956 long hdr_offset, long dat_offset)
957
958 {FITS_HDU_LIST *hdulist;
959 FITS_DATA *fdat;
960 char errmsg[80], key[9];
961 int k, bpp, random_groups;
962 long mul_axis, data_size, bitpix_supported;
963
964 #define FITS_DECODE_CARD(mhdr,mkey,mfdat,mtyp) \
965 {strcpy (key, mkey); \
966 mfdat = fits_decode_card (fits_search_card (mhdr, mkey), mtyp); \
967 if (mfdat == NULL) goto err_missing; }
968
969 #define FITS_TRY_CARD(mhdr,mhdu,mkey,mvar,mtyp,unionvar) \
970 {FITS_DATA *mfdat = fits_decode_card (fits_search_card (mhdr,mkey), mtyp); \
971 mhdu->used.mvar = (mfdat != NULL); \
972 if (mhdu->used.mvar) mhdu->mvar = mfdat->unionvar; }
973
974 hdulist = fits_new_hdulist ();
975 if (hdulist == NULL)
976 FITS_RETURN ("fits_decode_header: Not enough memory", NULL);
977
978 /* Initialize the header data */
979 hdulist->header_offset = hdr_offset;
980 hdulist->data_offset = dat_offset;
981
982 hdulist->used.simple = (strncmp (hdr->data, "SIMPLE ", 8) == 0);
983 hdulist->used.xtension = (strncmp (hdr->data, "XTENSION", 8) == 0);
984 if (hdulist->used.xtension)
985 {
986 fdat = fits_decode_card (fits_search_card (hdr, "XTENSION"), typ_fstring);
987 strcpy (hdulist->xtension, fdat->fstring);
988 }
989
990 FITS_DECODE_CARD (hdr, "NAXIS", fdat, typ_flong);
991 hdulist->naxis = fdat->flong;
992
993 FITS_DECODE_CARD (hdr, "BITPIX", fdat, typ_flong);
994 bpp = hdulist->bitpix = (int)fdat->flong;
995 if ( (bpp != 8) && (bpp != 16) && (bpp != 32)
996 && (bpp != -32) && (bpp != -64))
997 {
998 strcpy (errmsg, "fits_decode_header: Invalid BITPIX-value");
999 goto err_return;
1000 }
1001 if (bpp < 0) bpp = -bpp;
1002 bpp /= 8;
1003 hdulist->bpp = bpp;
1004
1005 FITS_TRY_CARD (hdr, hdulist, "GCOUNT", gcount, typ_flong, flong);
1006 FITS_TRY_CARD (hdr, hdulist, "PCOUNT", pcount, typ_flong, flong);
1007
1008 FITS_TRY_CARD (hdr, hdulist, "GROUPS", groups, typ_fbool, fbool);
1009 random_groups = hdulist->used.groups && hdulist->groups;
1010
1011 FITS_TRY_CARD (hdr, hdulist, "EXTEND", extend, typ_fbool, fbool);
1012
1013 if (hdulist->used.xtension) /* Extension requires GCOUNT and PCOUNT */
1014 {
1015 if ((!hdulist->used.gcount) || (!hdulist->used.pcount))
1016 {
1017 strcpy (errmsg, "fits_decode_header: Missing GCOUNT/PCOUNT for XTENSION");
1018 goto err_return;
1019 }
1020 }
1021
1022 mul_axis = 1;
1023
1024 /* Find all NAXISx-cards */
1025 for (k = 1; k <= FITS_MAX_AXIS; k++)
1026 {char naxisn[9];
1027
1028 snprintf (naxisn, sizeof( naxisn ), "NAXIS%-3d", k);
1029 fdat = fits_decode_card (fits_search_card (hdr, naxisn), typ_flong);
1030 if (fdat == NULL)
1031 {
1032 k--; /* Save the last NAXISk read */
1033 break;
1034 }
1035 hdulist->naxisn[k-1] = (int)fdat->flong;
1036 if (hdulist->naxisn[k-1] < 0)
1037 {
1038 strcpy (errmsg, "fits_decode_header: Negative value in NAXISn");
1039 goto err_return;
1040 }
1041 if ((k == 1) && (random_groups))
1042 {
1043 if (hdulist->naxisn[0] != 0)
1044 {
1045 strcpy (errmsg, "fits_decode_header: Random groups with NAXIS1 != 0");
1046 goto err_return;
1047 }
1048 }
1049 else
1050 mul_axis *= hdulist->naxisn[k-1];
1051 }
1052
1053 if ((hdulist->naxis > 0) && (k < hdulist->naxis))
1054 {
1055 strcpy (errmsg, "fits_decode_card: Not enough NAXISn-cards");
1056 goto err_return;
1057 }
1058
1059 /* If we have only one dimension, just set the second to size one. */
1060 /* So we dont have to check for naxis < 2 in some places. */
1061 if (hdulist->naxis < 2)
1062 hdulist->naxisn[1] = 1;
1063 if (hdulist->naxis < 1)
1064 {
1065 mul_axis = 0;
1066 hdulist->naxisn[0] = 1;
1067 }
1068
1069 if (hdulist->used.xtension)
1070 data_size = bpp*hdulist->gcount*(hdulist->pcount + mul_axis);
1071 else
1072 data_size = bpp*mul_axis;
1073 hdulist->udata_size = data_size; /* Used data size without padding */
1074
1075 /* Datasize must be a multiple of the FITS logical record size */
1076 data_size = (data_size + FITS_RECORD_SIZE - 1) / FITS_RECORD_SIZE;
1077 data_size *= FITS_RECORD_SIZE;
1078 hdulist->data_size = data_size;
1079
1080
1081 FITS_TRY_CARD (hdr, hdulist, "BLANK", blank, typ_flong, flong);
1082
1083 FITS_TRY_CARD (hdr, hdulist, "DATAMIN", datamin, typ_fdouble, fdouble);
1084 FITS_TRY_CARD (hdr, hdulist, "DATAMAX", datamax, typ_fdouble, fdouble);
1085
1086 FITS_TRY_CARD (hdr, hdulist, "BZERO", bzero, typ_fdouble, fdouble);
1087 FITS_TRY_CARD (hdr, hdulist, "BSCALE", bscale, typ_fdouble, fdouble);
1088
1089 /* Evaluate number of interpretable images for this HDU */
1090 hdulist->numpic = 0;
1091
1092 /* We must support this format */
1093 bitpix_supported = (hdulist->bitpix > 0)
1094 || ( (hdulist->bitpix == -64)
1095 && (fits_ieee64_intel || fits_ieee64_motorola))
1096 || ( (hdulist->bitpix == -32)
1097 && ( fits_ieee32_intel || fits_ieee32_motorola
1098 || fits_ieee64_intel || fits_ieee64_motorola));
1099
1100 if (bitpix_supported)
1101 {
1102 if (hdulist->used.simple)
1103 {
1104 if (hdulist->naxis > 0)
1105 {
1106 hdulist->numpic = 1;
1107 for (k = 3; k <= hdulist->naxis; k++)
1108 hdulist->numpic *= hdulist->naxisn[k-1];
1109 }
1110 }
1111 else if ( hdulist->used.xtension
1112 && (strncmp (hdulist->xtension, "IMAGE", 5) == 0))
1113 {
1114 if (hdulist->naxis > 0)
1115 {
1116 hdulist->numpic = 1;
1117 for (k = 3; k <= hdulist->naxis; k++)
1118 hdulist->numpic *= hdulist->naxisn[k-1];
1119 }
1120 }
1121 }
1122 else
1123 {char msg[160];
1124 snprintf (msg, sizeof( msg ), "fits_decode_header: IEEE floating point format required for\
1125 BITPIX=%d\nis not supported on this machine", hdulist->bitpix);
1126 fits_set_error (msg);
1127 }
1128
1129 hdulist->header_record_list = hdr; /* Add header records to the list */
1130 return (hdulist);
1131
1132 err_missing:
1133 snprintf (errmsg, sizeof(errmsg), "fits_decode_header: missing/invalid %.50s card", key);
1134
1135 err_return:
1136 fits_delete_hdulist (hdulist);
1137 fits_set_error (errmsg);
1138 return (NULL);
1139
1140 #undef FITS_DECODE_CARD
1141 }
1142
1143
1144 /*****************************************************************************/
1145 /* #BEG-PAR */
1146 /* */
1147 /* Function : fits_eval_pixrange - (local) evaluate range of pixel data */
1148 /* */
1149 /* Parameters: */
1150 /* FILE *fp [I] : file pointer */
1151 /* FITS_HDU_LIST *hdu [I] : pointer to header */
1152 /* ( mode : I=input, O=output, I/O=input/output ) */
1153 /* */
1154 /* The Function sets the values hdu->pixmin and hdu->pixmax. On success 0 */
1155 /* is returned. On failure, -1 is returned. */
1156 /* */
1157 /* #END-PAR */
1158 /*****************************************************************************/
1159
1160 static int fits_eval_pixrange (FILE *fp, FITS_HDU_LIST *hdu)
1161
1162 {register unsigned int maxelem;
1163 #define FITSNPIX 4096
1164 unsigned char pixdat[FITSNPIX];
1165 unsigned int nelem, bpp;
1166 int blank_found = 0, nan_found = 0;
1167
1168 if (fseek (fp, hdu->data_offset, SEEK_SET) < 0)
1169 FITS_RETURN ("fits_eval_pixrange: cant position file", -1);
1170
1171 bpp = hdu->bpp; /* Number of bytes per pixel */
1172 nelem = hdu->udata_size / bpp; /* Number of data elements */
1173
1174 switch (hdu->bitpix)
1175 {
1176 case 8: {
1177 register FITS_BITPIX8 pixval;
1178 register unsigned char *ptr;
1179 FITS_BITPIX8 minval = 255, maxval = 0;
1180 FITS_BITPIX8 blankval;
1181
1182 while (nelem > 0)
1183 {
1184 maxelem = sizeof (pixdat)/bpp;
1185 if (nelem < maxelem) maxelem = nelem;
1186 nelem -= maxelem;
1187 if (fread ((char *)pixdat, bpp, maxelem, fp) != maxelem)
1188 FITS_RETURN ("fits_eval_pixrange: error on read bitpix 8 data", -1);
1189
1190 ptr = pixdat;
1191 if (hdu->used.blank)
1192 {
1193 blankval = (FITS_BITPIX8)hdu->blank;
1194 while (maxelem-- > 0)
1195 {
1196 pixval = (FITS_BITPIX8)*(ptr++);
1197 if (pixval != blankval)
1198 {
1199 if (pixval < minval) minval = pixval;
1200 else if (pixval > maxval) maxval = pixval;
1201 }
1202 else blank_found = 1;
1203 }
1204 }
1205 else
1206 {
1207 while (maxelem-- > 0)
1208 {
1209 pixval = (FITS_BITPIX8)*(ptr++);
1210 if (pixval < minval) minval = pixval;
1211 else if (pixval > maxval) maxval = pixval;
1212 }
1213 }
1214 }
1215 hdu->pixmin = minval;
1216 hdu->pixmax = maxval;
1217 break; }
1218
1219 case 16: {
1220 register FITS_BITPIX16 pixval;
1221 register unsigned char *ptr;
1222 FITS_BITPIX16 minval = 0x7fff, maxval = ~0x7fff;
1223
1224 while (nelem > 0)
1225 {
1226 maxelem = sizeof (pixdat)/bpp;
1227 if (nelem < maxelem) maxelem = nelem;
1228 nelem -= maxelem;
1229 if (fread ((char *)pixdat, bpp, maxelem, fp) != maxelem)
1230 FITS_RETURN ("fits_eval_pixrange: error on read bitpix 16 data", -1);
1231
1232 ptr = pixdat;
1233 if (hdu->used.blank)
1234 {FITS_BITPIX16 blankval = (FITS_BITPIX16)hdu->blank;
1235
1236 while (maxelem-- > 0)
1237 {
1238 FITS_GETBITPIX16 (ptr, pixval);
1239 ptr += 2;
1240 if (pixval != blankval)
1241 {
1242 if (pixval < minval) minval = pixval;
1243 else if (pixval > maxval) maxval = pixval;
1244 }
1245 else blank_found = 1;
1246 }
1247 }
1248 else
1249 {
1250 while (maxelem-- > 0)
1251 {
1252 FITS_GETBITPIX16 (ptr, pixval);
1253 ptr += 2;
1254 if (pixval < minval) minval = pixval;
1255 else if (pixval > maxval) maxval = pixval;
1256 }
1257 }
1258 }
1259 hdu->pixmin = minval;
1260 hdu->pixmax = maxval;
1261 break; }
1262
1263
1264 case 32: {
1265 register FITS_BITPIX32 pixval;
1266 register unsigned char *ptr;
1267 FITS_BITPIX32 minval = 0x7fffffff, maxval = ~0x7fffffff;
1268
1269 while (nelem > 0)
1270 {
1271 maxelem = sizeof (pixdat)/bpp;
1272 if (nelem < maxelem) maxelem = nelem;
1273 nelem -= maxelem;
1274 if (fread ((char *)pixdat, bpp, maxelem, fp) != maxelem)
1275 FITS_RETURN ("fits_eval_pixrange: error on read bitpix 32 data", -1);
1276
1277 ptr = pixdat;
1278 if (hdu->used.blank)
1279 {FITS_BITPIX32 blankval = (FITS_BITPIX32)hdu->blank;
1280
1281 while (maxelem-- > 0)
1282 {
1283 FITS_GETBITPIX32 (ptr, pixval);
1284 ptr += 4;
1285 if (pixval != blankval)
1286 {
1287 if (pixval < minval) minval = pixval;
1288 else if (pixval > maxval) maxval = pixval;
1289 }
1290 else blank_found = 1;
1291 }
1292 }
1293 else
1294 {
1295 while (maxelem-- > 0)
1296 {
1297 FITS_GETBITPIX32 (ptr, pixval);
1298 ptr += 4;
1299 if (pixval < minval) minval = pixval;
1300 else if (pixval > maxval) maxval = pixval;
1301 }
1302 }
1303 }
1304 hdu->pixmin = minval;
1305 hdu->pixmax = maxval;
1306 break; }
1307
1308 case -32: {
1309 register FITS_BITPIXM32 pixval;
1310 register unsigned char *ptr;
1311 FITS_BITPIXM32 minval, maxval;
1312 int first = 1;
1313
1314 /* initialize */
1315
1316 pixval = 0;
1317 minval = 0;
1318 maxval = 0;
1319
1320 while (nelem > 0)
1321 {
1322 maxelem = sizeof (pixdat)/bpp;
1323 if (nelem < maxelem) maxelem = nelem;
1324 nelem -= maxelem;
1325 if (fread ((char *)pixdat, bpp, maxelem, fp) != maxelem)
1326 FITS_RETURN ("fits_eval_pixrange: error on read bitpix -32 data", -1);
1327
1328 ptr = pixdat;
1329 while (maxelem-- > 0)
1330 {
1331 if (!fits_nan_32 (ptr))
1332 {
1333 FITS_GETBITPIXM32 (ptr, pixval);
1334 ptr += 4;
1335 if (first)
1336 {
1337 first = 0;
1338 minval = maxval = pixval;
1339 }
1340 else if (pixval < minval) { minval = pixval; }
1341 else if (pixval > maxval) { maxval = pixval; }
1342 }
1343 else nan_found = 1;
1344 }
1345 }
1346 hdu->pixmin = minval;
1347 hdu->pixmax = maxval;
1348 break; }
1349
1350 case -64: {
1351 register FITS_BITPIXM64 pixval;
1352 register unsigned char *ptr;
1353 FITS_BITPIXM64 minval, maxval;
1354 int first = 1;
1355
1356 /* initialize */
1357
1358 minval = 0;
1359 maxval = 0;
1360
1361 while (nelem > 0)
1362 {
1363 maxelem = sizeof (pixdat)/bpp;
1364 if (nelem < maxelem) maxelem = nelem;
1365 nelem -= maxelem;
1366 if (fread ((char *)pixdat, bpp, maxelem, fp) != maxelem)
1367 FITS_RETURN ("fits_eval_pixrange: error on read bitpix -64 data", -1);
1368
1369 ptr = pixdat;
1370 while (maxelem-- > 0)
1371 {
1372 if (!fits_nan_64 (ptr))
1373 {
1374 FITS_GETBITPIXM64 (ptr, pixval);
1375 ptr += 8;
1376 if (first)
1377 {
1378 first = 0;
1379 minval = maxval = pixval;
1380 }
1381 else if (pixval < minval) { minval = pixval; }
1382 else if (pixval > maxval) { maxval = pixval; }
1383 }
1384 else nan_found = 1;
1385 }
1386 }
1387 hdu->pixmin = minval;
1388 hdu->pixmax = maxval;
1389 break; }
1390 }
1391 if (nan_found) hdu->used.nan_value = 1;
1392 if (blank_found) hdu->used.blank_value = 1;
1393
1394 return (0);
1395 }
1396
1397
1398 /*****************************************************************************/
1399 /* #BEG-PAR */
1400 /* */
1401 /* Function : fits_decode_card - decode a card */
1402 /* */
1403 /* Parameters: */
1404 /* const char *card [I] : pointer to card image */
1405 /* FITS_DATA_TYPES data_type [I] : datatype to decode */
1406 /* ( mode : I=input, O=output, I/O=input/output ) */
1407 /* */
1408 /* Decodes a card and returns a pointer to the union, keeping the data. */
1409 /* If card is NULL or on failure, a NULL-pointer is returned. */
1410 /* If the card does not have the value indicator, an error is generated, */
1411 /* but its tried to decode the card. The data is only valid up to the next */
1412 /* call of the function. */
1413 /* */
1414 /* #END-PAR */
1415 /*****************************************************************************/
1416
1417 FITS_DATA *fits_decode_card (const char *card, FITS_DATA_TYPES data_type)
1418
1419 {static FITS_DATA data;
1420 long l_long;
1421 double l_double;
1422 char l_card[FITS_CARD_SIZE+1], msg[256];
1423 char *cp = ident, *dst, *end;
1424
1425 if (card == NULL) return (NULL);
1426
1427 memcpy (l_card, card, FITS_CARD_SIZE);
1428 l_card[FITS_CARD_SIZE] = '\0';
1429
1430 if (strncmp (card+8, "= ", 2) != 0)
1431 {
1432 snprintf (msg, sizeof( msg ), "fits_decode_card (warning): Missing value indicator\
1433 '= ' for %8.8s", l_card);
1434 fits_set_error (msg);
1435 }
1436
1437 switch (data_type)
1438 {
1439 case typ_bitpix8:
1440 data.bitpix8 = (FITS_BITPIX8)(l_card[10]);
1441 break;
1442
1443 case typ_bitpix16:
1444 if (sscanf (l_card+10, "%ld", &l_long) != 1)
1445 FITS_RETURN ("fits_decode_card: error decoding typ_bitpix16", NULL);
1446 data.bitpix16 = (FITS_BITPIX16)l_long;
1447 break;
1448
1449 case typ_bitpix32:
1450 if (sscanf (l_card+10, "%ld", &l_long) != 1)
1451 FITS_RETURN ("fits_decode_card: error decoding typ_bitpix32", NULL);
1452 data.bitpix32 = (FITS_BITPIX32)l_long;
1453 break;
1454
1455 case typ_bitpixm32:
1456 if (sscanf (l_card+10, "%lf", &l_double) != 1)
1457 FITS_RETURN ("fits_decode_card: error decoding typ_bitpixm32", NULL);
1458 data.bitpixm32 = (FITS_BITPIXM32)l_double;
1459 break;
1460
1461 case typ_bitpixm64:
1462 if (sscanf (l_card+10, "%lf", &l_double) != 1)
1463 FITS_RETURN ("fits_decode_card: error decoding typ_bitpixm64", NULL);
1464 data.bitpixm64 = (FITS_BITPIXM64)l_double;
1465 break;
1466
1467 case typ_fbool:
1468 cp = l_card+10;
1469 while (*cp == ' ') cp++;
1470 if (*cp == 'T') data.fbool = 1;
1471 else if (*cp == 'F') data.fbool = 0;
1472 else FITS_RETURN ("fits_decode_card: error decoding typ_fbool", NULL);
1473 break;
1474
1475 case typ_flong:
1476 if (sscanf (l_card+10, "%ld", &l_long) != 1)
1477 FITS_RETURN ("fits_decode_card: error decoding typ_flong", NULL);
1478 data.flong = (FITS_BITPIX32)l_long;
1479 break;
1480
1481 case typ_fdouble:
1482 if (sscanf (l_card+10, "%lf", &l_double) != 1)
1483 FITS_RETURN ("fits_decode_card: error decoding typ_fdouble", NULL);
1484 data.fdouble = (FITS_BITPIXM32)l_double;
1485 break;
1486
1487 case typ_fstring:
1488 cp = l_card+10;
1489 if (*cp != '\'')
1490 FITS_RETURN ("fits_decode_card: missing \' decoding typ_fstring", NULL);
1491
1492 dst = data.fstring;
1493 cp++;
1494 end = l_card+FITS_CARD_SIZE-1;
1495 for (;;) /* Search for trailing quote */
1496 {
1497 if (*cp != '\'') /* All characters but quote are used. */
1498 {
1499 *(dst++) = *cp;
1500 }
1501 else /* Maybe there is a quote in the string */
1502 {
1503 if (cp >= end) break; /* End of card ? finished */
1504 if (*(cp+1) != '\'') break;
1505 *(dst++) = *(cp++);
1506 }
1507 if (cp >= end) break;
1508 cp++;
1509 }
1510 *dst = '\0';
1511 break;
1512 }
1513 return (&data);
1514 }
1515
1516
1517 /*****************************************************************************/
1518 /* #BEG-PAR */
1519 /* */
1520 /* Function : fits_search_card - search a card in the record list */
1521 /* */
1522 /* Parameters: */
1523 /* FITS_RECORD_LIST *rl [I] : record list to search */
1524 /* char *keyword [I] : keyword identifying the card */
1525 /* ( mode : I=input, O=output, I/O=input/output ) */
1526 /* */
1527 /* A card is searched in the reord list. Only the first eight characters of */
1528 /* keyword are significant. If keyword is less than 8 characters, its filled */
1529 /* with blanks. */
1530 /* If the card is found, a pointer to the card is returned. */
1531 /* The pointer does not point to a null-terminated string. Only the next */
1532 /* 80 bytes are allowed to be read. */
1533 /* On failure a NULL-pointer is returned. */
1534 /* */
1535 /* #END-PAR */
1536 /*****************************************************************************/
1537
1538 char *fits_search_card (FITS_RECORD_LIST *rl, const char *keyword)
1539
1540 {int key_len, k;
1541 char *card;
1542 char key[9];
1543
1544 key_len = strlen (keyword);
1545 if (key_len > 8) key_len = 8;
1546 if (key_len == 0)
1547 FITS_RETURN ("fits_search_card: Invalid parameter", NULL);
1548
1549 strcpy (key, " ");
1550 memcpy (key, keyword, key_len);
1551
1552 while (rl != NULL)
1553 {
1554 card = (char *)rl->data;
1555 for (k = 0; k < FITS_RECORD_SIZE / FITS_CARD_SIZE; k++)
1556 {
1557 if (strncmp (card, key, 8) == 0) return (card);
1558 card += FITS_CARD_SIZE;
1559 }
1560 rl = rl->next_record;
1561 }
1562 return (NULL);
1563 }
1564
1565
1566 /*****************************************************************************/
1567 /* #BEG-PAR */
1568 /* */
1569 /* Function : fits_image_info - get information about an image */
1570 /* */
1571 /* Parameters: */
1572 /* FITS_FILE *ff [I] : FITS file structure */
1573 /* int picind [I] : Index of picture in file (1,2,...) */
1574 /* int *hdupicind [O] : Index of picture in HDU (1,2,...) */
1575 /* ( mode : I=input, O=output, I/O=input/output ) */
1576 /* */
1577 /* The function returns on success a pointer to a FITS_HDU_LIST. hdupicind */
1578 /* then gives the index of the image within the HDU. */
1579 /* On failure, NULL is returned. */
1580 /* */
1581 /* #END-PAR */
1582 /*****************************************************************************/
1583
1584 FITS_HDU_LIST *fits_image_info (FITS_FILE *ff, int picind, int *hdupicind)
1585
1586 {FITS_HDU_LIST *hdulist;
1587 int firstpic, lastpic;
1588
1589 if (ff == NULL)
1590 FITS_RETURN ("fits_image_info: ff is NULL", NULL);
1591
1592 if (ff->openmode != 'r')
1593 FITS_RETURN ("fits_image_info: file not open for reading", NULL);
1594
1595 if ((picind < 1) || (picind > ff->n_pic))
1596 FITS_RETURN ("fits_image_info: picind out of range", NULL);
1597
1598 firstpic = 1;
1599 for (hdulist = ff->hdu_list; hdulist != NULL; hdulist = hdulist->next_hdu)
1600 {
1601 if (hdulist->numpic <= 0) continue;
1602 lastpic = firstpic+hdulist->numpic-1;
1603 if (picind <= lastpic) /* Found image in current HDU ? */
1604 break;
1605
1606 firstpic = lastpic+1;
1607 }
1608 *hdupicind = picind - firstpic + 1;
1609 return (hdulist);
1610 }
1611
1612
1613 /*****************************************************************************/
1614 /* #BEG-PAR */
1615 /* */
1616 /* Function : fits_seek_image - position to a specific image */
1617 /* */
1618 /* Parameters: */
1619 /* FITS_FILE *ff [I] : FITS file structure */
1620 /* int picind [I] : Index of picture to seek (1,2,...) */
1621 /* ( mode : I=input, O=output, I/O=input/output ) */
1622 /* */
1623 /* The function positions the file pointer to a specified image. */
1624 /* The function returns on success a pointer to a FITS_HDU_LIST. This pointer*/
1625 /* must also be used when reading data from the image. */
1626 /* On failure, NULL is returned. */
1627 /* */
1628 /* #END-PAR */
1629 /*****************************************************************************/
1630
1631 FITS_HDU_LIST *fits_seek_image (FITS_FILE *ff, int picind)
1632
1633 {FITS_HDU_LIST *hdulist;
1634 int hdupicind;
1635 long offset, pic_size;
1636
1637 hdulist = fits_image_info (ff, picind, &hdupicind);
1638 if (hdulist == NULL) return (NULL);
1639
1640 pic_size = hdulist->bpp * hdulist->naxisn[0] * hdulist->naxisn[1];
1641 offset = hdulist->data_offset + (hdupicind-1)*pic_size;
1642 if (fseek (ff->fp, offset, SEEK_SET) < 0)
1643 FITS_RETURN ("fits_seek_image: Unable to position to image", NULL);
1644
1645 return (hdulist);
1646 }
1647
1648
1649 /*****************************************************************************/
1650 /* #BEG-PAR */
1651 /* */
1652 /* Function : fits_read_pixel - read pixel values from a file */
1653 /* */
1654 /* Parameters: */
1655 /* FITS_FILE *ff [I] : FITS file structure */
1656 /* FITS_HDU_LIST *hdulist [I] : pointer to hdulist that describes image */
1657 /* int npix [I] : number of pixel values to read */
1658 /* FITS_PIX_TRANSFORM *trans [I]: pixel transformation */
1659 /* void *buf [O] : buffer where to place transformed pixels */
1660 /* ( mode : I=input, O=output, I/O=input/output ) */
1661 /* */
1662 /* The function reads npix pixel values from the file, transforms them */
1663 /* checking for blank/NaN pixels and stores the transformed values in buf. */
1664 /* The number of transformed pixels is returned. If the returned value is */
1665 /* less than npix (or even -1), an error has occured. */
1666 /* hdulist must be a pointer returned by fits_seek_image(). Before starting */
1667 /* to read an image, fits_seek_image() must be called. Even for successive */
1668 /* images. */
1669 /* */
1670 /* #END-PAR */
1671 /*****************************************************************************/
1672
1673 int fits_read_pixel (FITS_FILE *ff, FITS_HDU_LIST *hdulist, int npix,
1674 FITS_PIX_TRANSFORM *trans, void *buf)
1675
1676 {double offs, scale;
1677 double datadiff, pixdiff;
1678 unsigned char pixbuffer[4096], *pix, *cdata;
1679 unsigned char creplace;
1680 int transcount = 0;
1681 long tdata, tmin, tmax;
1682 int maxelem;
1683 FITS_BITPIX8 bp8, bp8blank;
1684 FITS_BITPIX16 bp16, bp16blank;
1685 FITS_BITPIX32 bp32, bp32blank;
1686 FITS_BITPIXM32 bpm32;
1687 FITS_BITPIXM64 bpm64;
1688
1689 /* initialize */
1690
1691 bpm32 = 0;
1692
1693 if (ff->openmode != 'r') return (-1); /* Not open for reading */
1694 if (trans->dsttyp != 'c') return (-1); /* Currently we only return chars */
1695 if (npix <= 0) return (npix);
1696
1697 datadiff = trans->datamax - trans->datamin;
1698 pixdiff = trans->pixmax - trans->pixmin;
1699
1700 offs = trans->datamin - trans->pixmin*datadiff/pixdiff;
1701 scale = datadiff / pixdiff;
1702
1703 tmin = (long)trans->datamin;
1704 tmax = (long)trans->datamax;
1705 if (tmin < 0) tmin = 0; else if (tmin > 255) tmin = 255;
1706 if (tmax < 0) tmax = 0; else if (tmax > 255) tmax = 255;
1707
1708 cdata = (unsigned char *)buf;
1709 creplace = (unsigned char)trans->replacement;
1710
1711 switch (hdulist->bitpix)
1712 {
1713 case 8:
1714 while (npix > 0) /* For all pixels to read */
1715 {
1716 maxelem = sizeof (pixbuffer) / hdulist->bpp;
1717 if (maxelem > npix) maxelem = npix;
1718 if (fread ((char *)pixbuffer, hdulist->bpp, maxelem, ff->fp) != (unsigned) maxelem)
1719 return (-1);
1720 npix -= maxelem;
1721
1722 pix = pixbuffer;
1723 if (hdulist->used.blank)
1724 {
1725 bp8blank = (FITS_BITPIX8)hdulist->blank;
1726 while (maxelem--)
1727 {
1728 bp8 = (FITS_BITPIX8)*(pix++);
1729 if (bp8 == bp8blank) /* Is it a blank pixel ? */
1730 *(cdata++) = creplace;
1731 else /* Do transform */
1732 {
1733 tdata = (long)(bp8 * scale + offs);
1734 if (tdata < tmin) tdata = tmin;
1735 else if (tdata > tmax) tdata = tmax;
1736 *(cdata++) = (unsigned char)tdata;
1737 }
1738 transcount++;
1739 }
1740 }
1741 else
1742 {
1743 while (maxelem--)
1744 {
1745 bp8 = (FITS_BITPIX8)*(pix++);
1746 tdata = (long)(bp8 * scale + offs);
1747 if (tdata < tmin) tdata = tmin;
1748 else if (tdata > tmax) tdata = tmax;
1749 *(cdata++) = (unsigned char)tdata;
1750 transcount++;
1751 }
1752 }
1753 }
1754 break;
1755
1756 case 16:
1757 while (npix > 0) /* For all pixels to read */
1758 {
1759 maxelem = sizeof (pixbuffer) / hdulist->bpp;
1760 if (maxelem > npix) maxelem = npix;
1761 if (fread ((char *)pixbuffer, hdulist->bpp, maxelem, ff->fp) != (unsigned) maxelem)
1762 return (-1);
1763 npix -= maxelem;
1764
1765 pix = pixbuffer;
1766 if (hdulist->used.blank)
1767 {
1768 bp16blank = (FITS_BITPIX16)hdulist->blank;
1769 while (maxelem--)
1770 {
1771 FITS_GETBITPIX16 (pix, bp16);
1772 if (bp16 == bp16blank)
1773 *(cdata++) = creplace;
1774 else
1775 {
1776 tdata = (long)(bp16 * scale + offs);
1777 if (tdata < tmin) tdata = tmin;
1778 else if (tdata > tmax) tdata = tmax;
1779 *(cdata++) = (unsigned char)tdata;
1780 }
1781 transcount++;
1782 pix += 2;
1783 }
1784 }
1785 else
1786 {
1787 while (maxelem--)
1788 {
1789 FITS_GETBITPIX16 (pix, bp16);
1790 tdata = (long)(bp16 * scale + offs);
1791 if (tdata < tmin) tdata = tmin;
1792 else if (tdata > tmax) tdata = tmax;
1793 *(cdata++) = (unsigned char)tdata;
1794 transcount++;
1795 pix += 2;
1796 }
1797 }
1798 }
1799 break;
1800
1801 case 32:
1802 while (npix > 0) /* For all pixels to read */
1803 {
1804 maxelem = sizeof (pixbuffer) / hdulist->bpp;
1805 if (maxelem > npix) maxelem = npix;
1806 if (fread ((char *)pixbuffer, hdulist->bpp, maxelem, ff->fp) != (unsigned) maxelem)
1807 return (-1);
1808 npix -= maxelem;
1809
1810 pix = pixbuffer;
1811 if (hdulist->used.blank)
1812 {
1813 bp32blank = (FITS_BITPIX32)hdulist->blank;
1814 while (maxelem--)
1815 {
1816 FITS_GETBITPIX32 (pix, bp32);
1817 if (bp32 == bp32blank)
1818 *(cdata++) = creplace;
1819 else
1820 {
1821 tdata = (long)(bp32 * scale + offs);
1822 if (tdata < tmin) tdata = tmin;
1823 else if (tdata > tmax) tdata = tmax;
1824 *(cdata++) = (unsigned char)tdata;
1825 }
1826 transcount++;
1827 pix += 4;
1828 }
1829 }
1830 else
1831 {
1832 while (maxelem--)
1833 {
1834 FITS_GETBITPIX32 (pix, bp32);
1835 tdata = (long)(bp32 * scale + offs);
1836 if (tdata < tmin) tdata = tmin;
1837 else if (tdata > tmax) tdata = tmax;
1838 *(cdata++) = (unsigned char)tdata;
1839 transcount++;
1840 pix += 4;
1841 }
1842 }
1843 }
1844 break;
1845
1846 case -32:
1847 while (npix > 0) /* For all pixels to read */
1848 {
1849 maxelem = sizeof (pixbuffer) / hdulist->bpp;
1850 if (maxelem > npix) maxelem = npix;
1851 if (fread ((char *)pixbuffer, hdulist->bpp, maxelem, ff->fp) != (unsigned) maxelem)
1852 return (-1);
1853 npix -= maxelem;
1854
1855 pix = pixbuffer;
1856 while (maxelem--)
1857 {
1858 if (fits_nan_32 (pix)) /* An IEEE special value ? */
1859 *(cdata++) = creplace;
1860 else /* Do transform */
1861 {
1862 FITS_GETBITPIXM32 (pix, bpm32);
1863 tdata = (long)(bpm32 * scale + offs);
1864 if (tdata < tmin) tdata = tmin;
1865 else if (tdata > tmax) tdata = tmax;
1866 *(cdata++) = (unsigned char)tdata;
1867 }
1868 transcount++;
1869 pix += 4;
1870 }
1871 }
1872 break;
1873
1874 case -64:
1875 while (npix > 0) /* For all pixels to read */
1876 {
1877 maxelem = sizeof (pixbuffer) / hdulist->bpp;
1878 if (maxelem > npix) maxelem = npix;
1879 if (fread ((char *)pixbuffer, hdulist->bpp, maxelem, ff->fp) != (unsigned) maxelem)
1880 return (-1);
1881 npix -= maxelem;
1882
1883 pix = pixbuffer;
1884 while (maxelem--)
1885 {
1886 if (fits_nan_64 (pix))
1887 *(cdata++) = creplace;
1888 else
1889 {
1890 FITS_GETBITPIXM64 (pix, bpm64);
1891 tdata = (long)(bpm64 * scale + offs);
1892 if (tdata < tmin) tdata = tmin;
1893 else if (tdata > tmax) tdata = tmax;
1894 *(cdata++) = (unsigned char)tdata;
1895 }
1896 transcount++;
1897 pix += 8;
1898 }
1899 }
1900 break;
1901 }
1902 return (transcount);
1903 }
1904
1905
1906
1907 #ifndef FITS_NO_DEMO
1908 /*****************************************************************************/
1909 /* #BEG-PAR */
1910 /* */
1911 /* Function : fits_to_pgmraw - convert FITS-file to raw PGM-file */
1912 /* */
1913 /* Parameters: */
1914 /* char *fitsfile [I] : name of fitsfile */
1915 /* char *pgmfile [I] : name of pgmfile */
1916 /* ( mode : I=input, O=output, I/O=input/output ) */
1917 /* */
1918 /* The function converts a FITS-file to a raw PGM-file. The PGM-file will */
1919 /* be upside down, because the orientation for storing the image is */
1920 /* different. On success, 0 is returned. On failure, -1 is returned. */
1921 /* */
1922 /* #END-PAR */
1923 /*****************************************************************************/
1924
1925 int fits_to_pgmraw (char *fitsfile, char *pgmfile)
1926
1927 {FITS_FILE *fitsin = NULL;
1928 FILE *pgmout = NULL;
1929 FITS_HDU_LIST *hdu;
1930 FITS_PIX_TRANSFORM trans;
1931 int retval = -1, nbytes, maxbytes;
1932 char buffer[1024];
1933
1934 fitsin = fits_open (fitsfile, "r"); /* Open FITS-file for reading */
1935 if (fitsin == NULL) goto err_return;
1936
1937 if (fitsin->n_pic < 1) goto err_return; /* Any picture in it ? */
1938
1939 hdu = fits_seek_image (fitsin, 1); /* Position to the first image */
1940 if (hdu == NULL) goto err_return;
1941 if (hdu->naxis < 2) goto err_return; /* Enough dimensions ? */
1942
1943 pgmout = fopen (pgmfile, "wb");
1944 if (pgmout == NULL) goto err_return;
1945
1946 /* Write PGM header with width/height */
1947 fprintf (pgmout, "P5\n%d %d\n255\n", hdu->naxisn[0], hdu->naxisn[1]);
1948
1949 /* Set up transformation for FITS pixel values to 0...255 */
1950 /* It maps trans.pixmin to trans.datamin and trans.pixmax to trans.datamax. */
1951 /* Values out of range [datamin, datamax] are clamped */
1952 trans.pixmin = hdu->pixmin;
1953 trans.pixmax = hdu->pixmax;
1954 trans.datamin = 0.0;
1955 trans.datamax = 255.0;
1956 trans.replacement = 0.0; /* Blank/NaN replacement value */
1957 trans.dsttyp = 'c'; /* Output type is character */
1958
1959 nbytes = hdu->naxisn[0]*hdu->naxisn[1];
1960 while (nbytes > 0)
1961 {
1962 maxbytes = sizeof (buffer);
1963 if (maxbytes > nbytes) maxbytes = nbytes;
1964
1965 /* Read pixels and transform them */
1966 if (fits_read_pixel (fitsin, hdu, maxbytes, &trans, buffer) != maxbytes)
1967 goto err_return;
1968
1969 if (fwrite (buffer, 1, maxbytes, pgmout) != maxbytes)
1970 goto err_return;
1971
1972 nbytes -= maxbytes;
1973 }
1974 retval = 0;
1975
1976 err_return:
1977
1978 if (fitsin) fits_close (fitsin);
1979 if (pgmout) fclose (pgmout);
1980
1981 return (retval);
1982 }
1983
1984
1985 /*****************************************************************************/
1986 /* #BEG-PAR */
1987 /* */
1988 /* Function : pgmraw to fits - convert raw PGM-file to FITS-file */
1989 /* */
1990 /* Parameters: */
1991 /* char *pgmfile [I] : name of pgmfile */
1992 /* char *fitsfile [I] : name of fitsfile */
1993 /* ( mode : I=input, O=output, I/O=input/output ) */
1994 /* */
1995 /* The function converts a raw PGM-file to a FITS-file. The FITS-file will */
1996 /* be upside down, because the orientation for storing the image is */
1997 /* different. On success, 0 is returned. On failure, -1 is returned. */
1998 /* */
1999 /* #END-PAR */
2000 /*****************************************************************************/
2001
2002 int pgmraw_to_fits (char *pgmfile, char *fitsfile)
2003
2004 {FITS_FILE *fitsout = NULL;
2005 FILE *pgmin = NULL;
2006 FITS_HDU_LIST *hdu;
2007 char buffer[1024];
2008 int width, height, numbytes, maxbytes;
2009 int retval = -1;
2010
2011 fitsout = fits_open (fitsfile, "w");
2012 if (fitsout == NULL) goto err_return;
2013
2014 pgmin = fopen (pgmfile, "r");
2015 if (pgmin == NULL) goto err_return;
2016
2017 /* Read signature of PGM file */
2018 if (fgets (buffer, sizeof (buffer), pgmin) == NULL) goto err_return;
2019 if ((buffer[0] != 'P') || (buffer[1] != '5')) goto err_return;
2020
2021 /* Skip comments upto width/height */
2022 do
2023 {
2024 if (fgets (buffer, sizeof (buffer), pgmin) == NULL) goto err_return;
2025 } while (buffer[0] == '#');
2026
2027 if (sscanf (buffer, "%d%d", &width, &height) != 2) goto err_return;
2028 if ((width < 1) || (height < 1)) goto err_return;
2029
2030 /* Skip comments upto maxval */
2031 do
2032 {
2033 if (fgets (buffer, sizeof (buffer), pgmin) == NULL) goto err_return;
2034 } while (buffer[0] == '#');
2035 /* Ignore maxval */
2036
2037 hdu = fits_add_hdu (fitsout); /* Create a HDU for the FITS file */
2038 if (hdu == NULL) goto err_return;
2039
2040 hdu->used.simple = 1; /* Set proper values */
2041 hdu->bitpix = 8;
2042 hdu->naxis = 2;
2043 hdu->naxisn[0] = width;
2044 hdu->naxisn[1] = height;
2045 hdu->used.datamin = 1;
2046 hdu->datamin = 0.0;
2047 hdu->used.datamax = 1;
2048 hdu->datamax = 255.0;
2049 hdu->used.bzero = 1;
2050 hdu->bzero = 0.0;
2051 hdu->used.bscale = 1;
2052 hdu->bscale = 1.0;
2053
2054 fits_add_card (hdu, "");
2055 fits_add_card (hdu, "HISTORY THIS FITS FILE WAS GENERATED BY FITSRW\
2056 USING PGMRAW_TO_FITS");
2057
2058 /* Write the header. Blocking is done automatically */
2059 if (fits_write_header (fitsout, hdu) < 0) goto err_return;
2060
2061 /* The primary array plus blocking must be written manually */
2062 numbytes = width * height;
2063
2064 while (numbytes > 0)
2065 {
2066 maxbytes = sizeof (buffer);
2067 if (maxbytes > numbytes) maxbytes = numbytes;
2068
2069 if (fread (buffer, 1, maxbytes, pgmin) != maxbytes) goto err_return;
2070 if (fwrite (buffer, 1, maxbytes, fitsout->fp) != maxbytes) goto err_return;
2071
2072 numbytes -= maxbytes;
2073 }
2074
2075 /* Do blocking */
2076 numbytes = (width * height) % FITS_RECORD_SIZE;
2077 if (numbytes)
2078 {
2079 while (numbytes++ < FITS_RECORD_SIZE)
2080 if (putc (0, fitsout->fp) == EOF) goto err_return;
2081 }
2082 retval = 0;
2083
2084 err_return:
2085
2086 if (fitsout) fits_close (fitsout);
2087 if (pgmin) fclose (pgmin);
2088
2089 return (retval);
2090 }
2091
2092 #endif
Experimental kdeedu stuff