1 /*
2 * xvfits.c - load/save routines for 'fits' format pictures
3 *
4 * LoadFITS(fname, pinfo, quick) - loads a FITS file
5 * WriteFITS(fp,pic,ptype,w,h,r,g,b,numcols,style,raw,cmt,comment)
6 */
7
8 /*
9 * Copyright 1992, 1993, 1994 by David Robinson.
10 *
11 * Permission to use, copy, modify, and distribute this software and its
12 * documentation for any purpose and without fee is hereby granted, provided
13 * that the above copyright notice appear in all copies. This software is
14 * provided "as is" without express or implied warranty.
15 */
16
17 #define NEEDSDIR /* for S_IRUSR|S_IWUSR */
18 #include "xv.h"
19
20 #define NCARDS (36)
21 #define BLOCKSIZE (2880)
22
23 /* data types */
24 typedef enum datatype { T_INT, T_LOG, T_STR, T_NOVAL } DATTYPE;
25
26 typedef struct {
27 FILE *fp; /* file pointer */
28 int bitpix; /* number of bits per pixel */
29 int size; /* Size of each pixel, in bytes */
30 int naxis; /* number of axes */
31 long int axes[3]; /* size of each axis */
32 long int ndata; /* number of elements in data */
33 long int cpos; /* current position in data file */
34 char *comment; /* malloc'ed comment string, or NULL if none */
35 } FITS;
36
37
38 /* block workspace of BLOCKSIZE bytes */
39 static char *fits_block=NULL;
40
41
42 static int splitfits PARM((byte *, char *, int, int, int, char *));
43 static const char *ftopen3d PARM((FITS *, char *, int *, int *, int *, int *));
44 static void ftclose PARM((FITS *));
45 static int ftgbyte PARM((FITS *, byte *, int));
46 static const char *rdheader PARM((FITS *));
47 static const char *wrheader PARM((FILE *, int, int, char *));
48 static const char *rdcard PARM((char *, const char *, DATTYPE, long int *));
49 static void wrcard PARM((char *, const char *, DATTYPE, int, char *));
50 static int ftgdata PARM((FITS *, void *, int));
51 static void ftfixdata PARM((FITS *, void *, int));
52 static void flip PARM((byte *, int, int));
53
54
55
56 /*******************************************/
LoadFITS(fname,pinfo,quick)57 int LoadFITS(fname, pinfo, quick)
58 char *fname;
59 PICINFO *pinfo;
60 int quick;
61 /*******************************************/
62 {
63 /* returns '1' on success */
64
65 FITS fs;
66 int i, nx, ny, nz, bitpix, nrd, ioerror, npixels, bufsize;
67 byte *image;
68 const char *error;
69 char basename[64];
70
71 if (fits_block == NULL) {
72 fits_block = (char *) malloc((size_t) BLOCKSIZE);
73 if (!fits_block) FatalError("Insufficient memory for FITS block buffer");
74 }
75
76 error = ftopen3d(&fs, fname, &nx, &ny, &nz, &bitpix);
77 if (error) {
78 SetISTR(ISTR_WARNING, "%s", error);
79 return 0;
80 }
81
82 if (quick) nz = 1; /* only load first plane */
83 npixels = nx * ny;
84 bufsize = nz * npixels;
85 if (nx <= 0 || ny <= 0 || npixels/nx != ny || bufsize/nz != npixels) {
86 SetISTR(ISTR_WARNING, "FITS image dimensions out of range (%dx%dx%d)",
87 nx, ny, nz);
88 return 0;
89 }
90
91 image = (byte *) malloc((size_t) bufsize);
92 if (!image) FatalError("Insufficient memory for image");
93
94 /*
95 * Read in image. For three dimensional images, read it in in one go
96 * to ensure that we get that same scaling for all planes.
97 */
98
99 nrd = ftgbyte(&fs, image, bufsize);
100 ioerror = ferror(fs.fp);
101 ftclose(&fs);
102
103 if (nrd == 0) { /* didn't read any data at all */
104 if (ioerror)
105 SetISTR(ISTR_WARNING, "%s", "I/O error reading FITS file");
106 else
107 SetISTR(ISTR_WARNING, "%s", "Unexpected EOF reading FITS file");
108
109 free(image);
110 return 0;
111 }
112
113 else if (nrd < bufsize) { /* read partial image */
114 if (ioerror)
115 SetISTR(ISTR_WARNING, "%s", "Truncated FITS file due to I/O error");
116 else
117 SetISTR(ISTR_WARNING, "%s", "Truncated FITS file");
118
119 { byte *foo;
120 for (foo=image+nrd; foo<image+bufsize; foo++) *foo=0x80; /* pad with grey */
121 }
122 }
123
124 if (nz > 1) {
125 /* how many planes do we actually have? */
126 nz = (nrd-1)/(npixels) + 1;
127
128 /* returns how many sub-files created */
129 nz = splitfits(image, fs.comment, nx, ny, nz, basename);
130 image = (byte *)realloc(image, (size_t) npixels); /* toss all but first */
131 }
132
133 /* There seems to be a convention that fits files be displayed using
134 * a cartesian coordinate system. Thus the first pixel is in the lower left
135 * corner. Fix this by reflecting in the line y=ny/2.
136 */
137 flip(image, nx, ny);
138
139 /* Success! */
140 pinfo->pic = image;
141 pinfo->type = PIC8;
142 pinfo->w = pinfo->normw = nx;
143 pinfo->h = pinfo->normh = ny;
144
145 for (i=0; i < 256; i++) pinfo->r[i] = pinfo->g[i] = pinfo->b[i] = i;
146 pinfo->frmType = F_FITS;
147 pinfo->colType = F_GREYSCALE;
148
149 sprintf(pinfo->fullInfo, "FITS, bitpix: %d", bitpix);
150 sprintf(pinfo->shrtInfo, "%dx%d FITS.", nx, ny);
151 pinfo->comment = fs.comment;
152
153 if (nz > 1) {
154 pinfo->numpages = nz;
155 strcpy(pinfo->pagebname, basename);
156 }
157
158 return 1;
159 }
160
161
162
163 /*******************************************/
WriteFITS(fp,pic,ptype,w,h,rmap,gmap,bmap,numcols,colorstyle,comment)164 int WriteFITS(fp,pic,ptype,w,h,rmap,gmap,bmap,numcols,colorstyle,comment)
165 FILE *fp;
166 byte *pic;
167 int ptype, w,h;
168 byte *rmap, *gmap, *bmap;
169 int numcols, colorstyle;
170 char *comment;
171 {
172 int i, j, npixels, nend;
173 byte *ptr;
174 const char *error;
175 byte rgb[256];
176
177 if (!fits_block) {
178 fits_block = (char *) malloc((size_t) BLOCKSIZE);
179 if (!fits_block) FatalError("Insufficient memory for FITS block buffer");
180 }
181
182 error = wrheader(fp, w, h, comment);
183 if (error) {
184 SetISTR(ISTR_WARNING, "%s", error);
185 return -1;
186 }
187
188 if (ptype == PIC8) {
189 /* If ptype is PIC8, we need to calculate the greyscale colourmap */
190 for (i=0; i < numcols; i++) rgb[i] = MONO(rmap[i], gmap[i], bmap[i]);
191
192 for (i=h-1; i >= 0; i--) { /* flip line ordering when writing out */
193 ptr = &pic[i*w];
194 for (j=0; j < w; j++, ptr++) putc(rgb[*ptr], fp);
195 }
196 }
197
198 else { /* PIC24 */
199 for (i=h-1; i >= 0; i--) { /* flip line ordering when writing out */
200 ptr = &pic[i*w*3];
201 for (j=0; j < w; j++, ptr += 3) putc(MONO(ptr[0], ptr[1], ptr[2]), fp);
202 }
203 }
204
205 npixels = w*h;
206
207 /* nend is the number of padding characters at the end of the last block */
208 nend = ((npixels+BLOCKSIZE-1)/BLOCKSIZE)*BLOCKSIZE - npixels;
209 if (nend) for (i=0; i<nend; i++) putc('\0', fp);
210
211 return 0;
212 }
213
214
215
216 /************************************/
splitfits(image,comment,nx,ny,nz,basename)217 static int splitfits(image, comment, nx, ny, nz, basename)
218 byte *image;
219 char *comment;
220 int nx, ny, nz;
221 char *basename;
222 {
223 /*
224 * Given a 3-dimensional FITS image, this splits it up into nz 2-d files.
225 * It returns the number of files actually stored.
226 * If only one file could be written, then no split files are created.
227 * It returns the basename of the split files in bname.
228 * If there was a problem writing files, then a error message will be set.
229 */
230
231 int i, npixels=nx * ny, nwrt, tmpfd;
232 FILE *fp;
233 const char *error;
234 char filename[70];
235
236 #ifndef VMS
237 sprintf(basename, "%s/xvpgXXXXXX", tmpdir);
238 #else
239 sprintf(basename, "Sys$Disk:[]xvpgXXXXXX");
240 #endif
241
242 #ifdef USE_MKSTEMP
243 close(mkstemp(basename));
244 #else
245 mktemp(basename);
246 #endif
247 if (basename[0] == '\0') {
248 SetISTR(ISTR_WARNING, "%s", "Unable to build temporary filename");
249 return 1;
250 }
251
252 strcat(basename, ".");
253 error = NULL;
254
255 for (i=0; i < nz && !error; i++) {
256 sprintf(filename, "%s%d", basename, i+1);
257 tmpfd = open(filename,O_WRONLY|O_CREAT|O_EXCL,S_IRWUSR);
258 if (tmpfd < 0) {
259 error = "Unable to open temporary file";
260 break;
261 }
262 fp = fdopen(tmpfd, "w");
263 if (!fp) {
264 error = "Unable to open temporary file";
265 break;
266 }
267
268 if (wrheader(fp, nx, ny, comment)) {
269 error = "I/O error writing temporary file";
270 fflush(fp);
271 fclose(fp);
272 unlink(filename);
273 close(tmpfd);
274 break;
275 }
276
277 nwrt = fwrite(image+i*npixels, sizeof(byte), (size_t) npixels, fp);
278 fflush(fp);
279 fclose(fp);
280 close(tmpfd);
281
282 if (nwrt == 0) { /* failed to write any data */
283 error = "I/O error writing temporary file";
284 unlink(filename);
285 break;
286 }
287 else if (nwrt < npixels)
288 error = "I/O error writing temporary file";
289 }
290
291
292 /* at this point, i is the number of files created */
293 if (i == 0) return 1;
294 if (i == 1) {
295 sprintf(filename, "%s1", basename);
296 unlink(filename);
297 }
298
299 if (error) SetISTR(ISTR_WARNING, "%s", error);
300 return i;
301 }
302
303
304 /************************************/
wrheader(fp,nx,ny,comment)305 static const char *wrheader(fp, nx, ny, comment)
306 FILE *fp;
307 int nx, ny;
308 char *comment;
309 {
310 /* Writes a minimalist FITS file header */
311
312 char *block = fits_block, *bp;
313 int i, j, lenhist;
314 char history[80];
315
316 for (i=0, bp=block; i<BLOCKSIZE; i++, bp++) *bp = ' ';
317
318 sprintf(history, "Written by XV %s", VERSTR);
319 lenhist = strlen(history);
320
321 i = 0;
322 wrcard(&block[80*i++], "SIMPLE", T_LOG, 1, NULL); /* write SIMPLE card */
323 wrcard(&block[80*i++], "BITPIX", T_INT, 8, NULL); /* write BITPIX card */
324 wrcard(&block[80*i++], "NAXIS", T_INT, 2, NULL); /* write NAXIS card */
325 wrcard(&block[80*i++], "NAXIS1", T_INT, nx, NULL); /* write NAXIS1 card */
326 wrcard(&block[80*i++], "NAXIS2", T_INT, ny, NULL); /* write NAXIS2 card */
327
328 /* Write HISTORY keyword */
329 wrcard(&block[80*i++], "HISTORY", T_STR, lenhist, history);
330
331 if (comment && *comment != '\0') {
332 while (*comment == '\n') comment++; /* Skip any blank lines */
333 while (*comment != '\0') {
334 for (j=0; j<72; j++)
335 if (comment[j] == '\0' || comment[j] == '\n') break;
336
337 /*
338 * Check to see if it is an xv history record; if so, then avoid
339 * duplicating it.
340 */
341 if (j != lenhist || xvbcmp(comment, history, (size_t) j) != 0)
342 wrcard(&block[80*i++], "COMMENT", T_STR, j, comment);
343
344 if (i == NCARDS) { /* Filled up a block */
345 i = fwrite(block, sizeof(char), (size_t) BLOCKSIZE, fp);
346 if (i != BLOCKSIZE) return("Error writing FITS file");
347 for (i=0, bp=block; i<BLOCKSIZE; i++, bp++) *bp = ' ';
348 i = 0;
349 }
350
351 comment += j;
352 while (*comment == '\n') comment++; /* Skip any blank lines */
353 }
354 }
355
356 wrcard(&block[80*i++], "END", T_NOVAL, 0, NULL); /* write END keyword */
357 i = fwrite(block, sizeof(char), (size_t) BLOCKSIZE, fp);
358
359 if (i != BLOCKSIZE) return "Error writing FITS file";
360 return NULL;
361 }
362
363
364
365 /************************************/
ftopen3d(fs,file,nx,ny,nz,bitpix)366 static const char *ftopen3d(fs, file, nx, ny, nz, bitpix)
367 FITS *fs;
368 char *file;
369 int *nx, *ny, *nz, *bitpix;
370 {
371 /* open a 2 or 3-dimensional fits file.
372 * Stores the dimensions of the file in nx, ny and nz, and updates the FITS
373 * structure passed in fs.
374 * If successful, returns NULL otherwise returns an error message.
375 * Will return an error message if the primary data unit is not a
376 * 2 or 3-dimensional array.
377 */
378
379 FILE *fp;
380 int naxis, i;
381 const char *error;
382
383 fp = xv_fopen(file, "r");
384 if (!fp) return "Unable to open FITS file";
385
386 fs->fp = fp;
387 fs->bitpix = 0;
388 fs->naxis = 0;
389 fs->cpos = 0;
390
391 /* read header */
392 error = rdheader(fs);
393 if (error) {
394 ftclose(fs);
395 return error;
396 }
397
398 naxis = fs->naxis;
399
400 /* get number of data */
401 fs->ndata = 1;
402 for (i=0; i<naxis; i++)
403 fs->ndata = fs->ndata * fs->axes[i];
404
405 *nx = fs->axes[0];
406 *ny = fs->axes[1];
407 if (naxis == 2) *nz = 1;
408 else *nz = fs->axes[2];
409
410 *bitpix = fs->bitpix;
411
412 return NULL;
413 }
414
415
416 /************************************/
ftclose(fs)417 static void ftclose(fs)
418 FITS *fs;
419 {
420 if (fs == NULL) return;
421 if (fs->fp != NULL) fclose(fs->fp);
422 }
423
424
425 /************************************/
rdheader(fs)426 static const char *rdheader(fs)
427 FITS *fs;
428 {
429 /* reads the fits header, and updates the FITS structure fs.
430 * Returns NULL on success, or an error message otherwise.
431 */
432
433 int i, j, res, commlen, commsize;
434 char name[9];
435 char *block=fits_block, *p;
436 const char *error;
437 long int val; /* the value */
438
439 fs->comment = NULL;
440 commlen = 0;
441 commsize = 256;
442
443 res = fread(block, sizeof(char), (size_t) BLOCKSIZE, fs->fp);
444 if (res != BLOCKSIZE) return "Error reading FITS file";
445 i = 0;
446
447 /* read SIMPLE key */
448 error = rdcard(block, "SIMPLE", T_LOG, &val);
449 if (error) return error;
450 if (val == 0) return "Not a SIMPLE FITS file";
451 i++;
452
453 /* read BITPIX key */
454 error = rdcard(&block[80], "BITPIX", T_INT, &val);
455 if (error) return error;
456
457 if (val != 8 && val != 16 && val != 32 && val != 64 && val != -32 &&
458 val != -64)
459 return "Bad BITPIX value in FITS file";
460
461 j = fs->bitpix = val;
462 if (j<0) j = -j;
463 fs->size = j/8;
464 i++;
465
466 /* read NAXIS key */
467 error = rdcard(&block[2*80], "NAXIS", T_INT, &val);
468 if (error) return error;
469 if (val < 0 || val > 999) return "Bad NAXIS value in FITS file";
470 if (val == 0) return "FITS file does not contain an image";
471 if (val < 2) return "FITS file has fewer than two dimensions";
472 fs->naxis = val;
473 i++;
474
475 /* read NAXISnnn keys.
476 * We allow NAXIS to be > 3 iff the dimensions of the extra axes are 1
477 */
478 for (j=0; j < fs->naxis; j++) {
479 if (i == NCARDS) {
480 res = fread(block, sizeof(char), (size_t) BLOCKSIZE, fs->fp);
481 if (res != BLOCKSIZE) return "Error reading FITS file";
482 i = 0;
483 }
484
485 sprintf(name, "NAXIS%d", j+1);
486 error = rdcard(&block[i*80], name, T_INT, &val);
487 if (error) return error;
488 if (val < 0) return "Bad NAXISn value in FITS file";
489 if (val == 0) return "FITS file does not contain an image";
490
491 if (j < 3) fs->axes[j] = val;
492 else if (val != 1) return "FITS file has more than three dimensions";
493 i++;
494 }
495 if (fs->naxis > 3) fs->naxis = 3;
496
497 /* do remainder, looking for comment cards */
498 /* Section 5.2.2.4 of the NOST standard groups COMMENT and HISTORY keywords
499 * under 'Commentary Keywords'. Many applications write HISTORY records in
500 * preference to COMMENT records, so we recognise both types here.
501 */
502 for (;;) {
503 if (i == NCARDS) {
504 res = fread(block, sizeof(char), (size_t) BLOCKSIZE, fs->fp);
505 if (res != BLOCKSIZE) return "Unexpected eof in FITS file";
506 i = 0;
507 }
508
509 p = &block[i*80];
510 if (strncmp(p, "END ", (size_t) 8) == 0) break;
511 if (strncmp(p, "HISTORY ", (size_t) 8) == 0 ||
512 strncmp(p, "COMMENT ", (size_t) 8) == 0) {
513 p += 8; /* skip keyword */
514 for (j=71; j >= 0; j--) if (p[j] != ' ') break;
515 j++; /* make j length of comment */
516 if (j > 0) { /* skip blank comment cards */
517 if (fs->comment == NULL) {
518 fs->comment = (char *) malloc((size_t) commsize); /* initially 256 */
519 if (fs->comment == NULL)
520 FatalError("Insufficient memory for comment buffer");
521 }
522
523 if (commlen + j + 2 > commsize) { /* if too small */
524 char *new;
525 commsize += commsize; /* double size of array */
526 new = (char *) malloc((size_t) commsize);
527
528 if (new == NULL)
529 FatalError("Insufficient memory for comment buffer");
530
531 if (commlen) xvbcopy(fs->comment, new, (size_t) commlen);
532 free(fs->comment);
533 fs->comment = new;
534 }
535
536 xvbcopy(p, &fs->comment[commlen], (size_t) j); /* add string */
537 commlen += j;
538 fs->comment[commlen++] = '\n'; /* with trailing cr */
539 fs->comment[commlen] = '\0';
540 }
541 }
542 i++;
543 }
544
545 return NULL;
546 }
547
548
549 /************************************/
wrcard(card,name,dtype,kvalue,svalue)550 static void wrcard(card, name, dtype, kvalue, svalue)
551 char *card;
552 const char *name;
553 DATTYPE dtype; /* type of value */
554 int kvalue;
555 char *svalue;
556 {
557 /* write a header record into the 80 byte buffer card.
558 * The keyword name is passed in name. The value type is in dtype; this
559 * can have the following values:
560 * dtype = T_NOVAL
561 * no keyword value is written
562 * dtype = T_LOG
563 * a logical value, either 'T' or 'F' in column 30 is written
564 * dtype = T_INT
565 * an integer is written, right justified in columns 11-30
566 * dtype = T_STR
567 * a string 'svalue' of length kvalue is written, in columns 9-80.
568 */
569
570 int l;
571 char *sp;
572
573 for (l=0, sp=card; l<80; l++,sp++) *sp=' ';
574
575 l = strlen(name);
576 if (l) xvbcopy(name, card, (size_t) l); /* copy name */
577
578 if (dtype == T_NOVAL) return;
579
580 if (dtype == T_STR) {
581 l = kvalue;
582 if (l <= 0) return;
583 if (l > 72) l = 72;
584 xvbcopy(svalue, &card[8], (size_t) l);
585 return;
586 }
587
588 card[8] = '=';
589
590 if (dtype == T_LOG)
591 card[29] = kvalue ? 'T' : 'F';
592 else { /* T_INT */
593 sprintf(&card[10], "%20d", kvalue);
594 card[30] = ' ';
595 }
596 }
597
598
599 /************************************/
rdcard(card,name,dtype,kvalue)600 static const char *rdcard(card, name, dtype, kvalue)
601 char *card;
602 const char *name;
603 DATTYPE dtype; /* type of value */
604 long int *kvalue;
605 {
606 /* Read a header record, from the 80 byte buffer card.
607 * the keyword name must match 'name'; and parse its value according to
608 * dtype. This can have the following values:
609 * dtype = T_LOG
610 * value is logical, either 'T' or 'F' in column 30.
611 * dtype = T_INT
612 * value is an integer, right justified in columns 11-30.
613 *
614 * The value is stored in kvalue.
615 * It returns NULL on success, or an error message otherwise.
616 */
617
618 int i, ptr;
619 char namestr[9];
620 static char error[45];
621
622 xvbcopy(card, namestr, (size_t) 8);
623
624 for (i=7; i>=0 && namestr[i] == ' '; i--);
625 namestr[i+1] = '\0';
626
627 if (strcmp(namestr, name) != 0) {
628 sprintf(error, "Keyword %s not found in FITS file", name);
629 return error;
630 }
631
632
633 /* get start of value */
634 ptr = 10;
635 while (ptr < 80 && card[ptr] == ' ') ptr++;
636 if (ptr == 80) return "FITS file has missing keyword value"; /* no value */
637
638 if (dtype == T_LOG) {
639 if (ptr != 29 || (card[29] != 'T' && card[29] != 'F'))
640 return "Keyword has bad logical value in FITS file";
641 *kvalue = (card[29] == 'T');
642 }
643
644 else { /* an integer */
645 int j;
646 long int ival;
647 char num[21];
648
649 if (ptr > 29) return "Keyword has bad integer value in FITS file";
650 xvbcopy(&card[ptr], num, (size_t) (30-ptr));
651 num[30-ptr] = '\0';
652 j = sscanf(num, "%ld", &ival);
653 if (j != 1) return "Keyword has bad integer value in FITS file";
654 *kvalue = ival;
655 }
656
657 return NULL;
658 }
659
660
661 /************************************/
ftgdata(fs,buffer,nelem)662 static int ftgdata(fs, buffer, nelem)
663 FITS *fs;
664 void *buffer;
665 int nelem;
666 {
667 /* reads nelem values into the buffer.
668 * returns NULL for success or an error message.
669 * Copes with the fact that the last 2880 byte record of the FITS file
670 * may be truncated, and should be padded out with zeros.
671 * bitpix type of data
672 * 8 byte
673 * 16 short int (NOT 2-byte integer)
674 * 32 int (NOT 4-byte integer)
675 * -32 float
676 * -64 double
677 *
678 * Returns the number of elements actually read.
679 */
680
681 int res;
682
683 if (nelem == 0) return 0;
684
685 res = fread(buffer, (size_t) fs->size, (size_t) nelem, fs->fp);
686 /* if failed to read all the data because at end of file */
687 if (res != nelem && feof(fs->fp)) {
688 /* nblock is the number of elements in a record.
689 size is always a factor of BLOCKSIZE */
690
691 int loffs, nblock=BLOCKSIZE/fs->size;
692
693 /*
694 * the last record might be short; check this.
695 * loffs is the offset of the start of the last record from the current
696 * position.
697 */
698
699 loffs = ((fs->ndata + nblock - 1) / nblock - 1) * nblock - fs->cpos;
700
701 /* if we read to the end of the penultimate record */
702 if (res >= loffs) {
703 /* pad with zeros */
704 xvbzero((char *)buffer+res*fs->size, (size_t) ((nelem-res)*fs->size));
705 res = nelem;
706 }
707 }
708
709 fs->cpos += res;
710 if (res) ftfixdata(fs, buffer, res);
711 return res;
712 }
713
714
715 /************************************/
ftfixdata(fs,buffer,nelem)716 static void ftfixdata(fs, buffer, nelem)
717 FITS *fs;
718 void *buffer;
719 int nelem;
720 {
721 /* convert the raw data, as stored in the FITS file, to the format
722 * appropiate for the data representation of the host computer.
723 * Assumes that
724 * short int = 2 or more byte integer
725 * int = 4 or more byte integer
726 * float = 4 byte floating point, not necessarily IEEE.
727 * double = 8 byte floating point.
728 * This can exit for lack of memory on the first call with float or
729 * double data.
730 */
731
732 int i, n=nelem;
733 byte *ptr=buffer;
734
735 /*
736 * conversions. Although the data may be signed, reverse using unsigned
737 * variables.
738 * Because the native int types may be larger than the types in the file,
739 * we start from the end and work backwards to avoid overwriting data
740 * prematurely.
741 */
742 /* convert from big-endian two-byte signed integer to native form */
743
744 if (fs->bitpix == 16) {
745 unsigned short int *iptr=(unsigned short int *)ptr;
746 iptr += n-1; /* last short int */
747 ptr += (n-1)*2; /* last pair of bytes */
748 for (i=0; i < n; i++, ptr-=2, iptr--)
749 *iptr = (((int)*ptr) << 8) | (int)(ptr[1]);
750 }
751
752 /* convert from big-endian four-byte signed integer to native form */
753 else if (fs->bitpix == 32) {
754 unsigned int *iptr = (unsigned int *)ptr;
755 iptr += n-1; /* last integer */
756 ptr += (n-1)*4; /* last 4 bytes */
757 for (i=0; i < n; i++, ptr-=4, iptr--)
758 *iptr = ((unsigned int)ptr[0] << 24) |
759 ((unsigned int)ptr[1] << 16) |
760 ((unsigned int)ptr[2] << 8) |
761 ((unsigned int)ptr[3]);
762 }
763
764 /* convert from IEE 754 single precision to native form */
765 else if (fs->bitpix == -32) {
766 int j, k, expo;
767 static float *exps=NULL;
768
769 if (exps == NULL) {
770 exps = (float *)malloc(256 * sizeof(float));
771 if (exps == NULL) FatalError("Insufficient memory for exps store");
772 exps[150] = 1.;
773 for (i=151; i < 256; i++) exps[i] = 2.*exps[i-1];
774 for (i=149; i >= 0; i--) exps[i] = 0.5*exps[i+1];
775 }
776
777 for (i=0; i < n; i++, ptr+=4) {
778 k = (int)*ptr;
779 j = ((int)ptr[1] << 16) | ((int)ptr[2] << 8) | (int)ptr[3];
780 expo = ((k & 127) << 1) | (j >> 23);
781 if ((expo | j) == 0) *(float *)ptr = 0.;
782 else *(float *)ptr = exps[expo]*(float)(j | 0x800000);
783 if (k & 128) *(float *)ptr = - *(float *)ptr;
784 }
785
786 }
787
788 /* convert from IEE 754 double precision to native form */
789 else if (fs->bitpix == -64) {
790 int expo, k, l;
791 unsigned int j;
792 static double *exps=NULL;
793
794 if (exps == NULL) {
795 exps = (double *)malloc(2048 * sizeof(double));
796 if (exps == NULL) FatalError("Insufficient memory for exps store");
797 exps[1075] = 1.;
798 for (i=1076; i < 2048; i++) exps[i] = 2.*exps[i-1];
799 for (i=1074; i >= 0; i--) exps[i] = 0.5*exps[i+1];
800 }
801
802 for (i=0; i < n; i++, ptr+=8) {
803 k = (int)*ptr;
804 j = ((unsigned int)ptr[1] << 24) | ((unsigned int)ptr[2] << 16) |
805 ((unsigned int)ptr[3] << 8) | (unsigned int)ptr[4];
806 l = ((int)ptr[5] << 16) | ((int)ptr[6] << 8) | (int)ptr[7];
807 expo = ((k & 127) << 4) | (j >> 28);
808 if ((expo | j | l) == 0) *(double *)ptr = 0.;
809 else *(double *)ptr = exps[expo] * (16777216. *
810 (double)((j&0x0FFFFFFF)|0x10000000) + (double)l);
811 if (k & 128) *(double *)ptr = - *(double *)ptr;
812 }
813 }
814 }
815
816 #define maxmin(x, max, min) {\
817 maxmin_t=(x); if(maxmin_t > max) max=maxmin_t; \
818 if (maxmin_t<min) min=maxmin_t;}
819
820
821 /************************************/
ftgbyte(fs,cbuff,nelem)822 static int ftgbyte(fs, cbuff, nelem)
823 FITS *fs;
824 byte *cbuff;
825 int nelem;
826 {
827 /* Reads a byte image from the FITS file fs. The image contains nelem pixels.
828 * If bitpix = 8, then the image is loaded as stored in the file.
829 * Otherwise, it is rescaled so that the minimum value is stored as 0, and
830 * the maximum is stored as 255.
831 * Returns the number of pixels read.
832 */
833
834 void *voidbuff;
835 int i, n, nrd, bufsize, overflow=0;
836
837 /* if the data is byte, then read it directly */
838 if (fs->bitpix == 8)
839 return ftgdata(fs, cbuff, nelem);
840
841 /* allocate a buffer to store the image */
842 if (fs->bitpix == 16) {
843 bufsize = nelem * sizeof(short int);
844 if (bufsize/nelem != (int)sizeof(short int))
845 overflow = 1;
846 } else if (fs->bitpix == 32) {
847 bufsize = nelem * sizeof(int);
848 if (bufsize/nelem != (int)sizeof(short int))
849 overflow = 1;
850 } else {
851 bufsize = nelem * fs->size; /* float, double */
852 if (bufsize/nelem != fs->size)
853 overflow = 1;
854 }
855
856 if (overflow) {
857 SetISTR(ISTR_WARNING, "FITS image dimensions out of range");
858 return 0;
859 }
860
861 voidbuff = (void *)malloc((size_t) bufsize);
862 if (voidbuff == NULL) {
863 char emess[60];
864 sprintf(emess, "Insufficient memory for raw image of %d bytes",
865 nelem*fs->size);
866 FatalError(emess);
867 }
868
869 nrd = ftgdata(fs, voidbuff, nelem);
870 if (nrd == 0) return 0;
871 n = nrd;
872
873 /* convert short int to byte */
874 if (fs->bitpix == 16) {
875 short int *buffer=voidbuff;
876 int max, min, maxmin_t;
877 float scale;
878
879 min = max = buffer[0];
880 for (i=1; i < n; i++, buffer++) maxmin(*buffer, max, min);
881 scale = (max == min) ? 0. : 255./(float)(max-min);
882
883 /* rescale and convert */
884 for (i=0, buffer=voidbuff; i < n; i++)
885 cbuff[i] = (byte)(scale*(float)((int)buffer[i]-min));
886
887 /* convert long int to byte */
888 }
889
890 else if (fs->bitpix == 32) {
891 int *buffer=voidbuff;
892 int max, min, maxmin_t;
893 float scale, fmin;
894
895 min = max = buffer[0];
896 for (i=1; i < n; i++, buffer++) maxmin(*buffer, max, min);
897 scale = (max == min) ? 1. : 255./((double)max-(double)min);
898 fmin = (float)min;
899
900 /* rescale and convert */
901 if (scale < 255./2.1e9) /* is max-min too big for an int ? */
902 for (i=0, buffer=voidbuff; i < n; i++)
903 cbuff[i] = (byte)(scale*((float)buffer[i]-fmin));
904 else /* use integer subtraction */
905 for (i=0, buffer=voidbuff; i < n; i++)
906 cbuff[i] = (byte)(scale*(float)(buffer[i]-min));
907
908
909 }
910
911 /* convert float to byte */
912 else if (fs->bitpix == -32) {
913 float *buffer=voidbuff;
914 float max, min, maxmin_t, scale;
915
916 min = max = buffer[0];
917 for (i=1; i < n; i++, buffer++) maxmin(*buffer, max, min);
918 scale = (max == min) ? 0. : 255./(max-min);
919
920 /* rescale and convert */
921 for (i=0, buffer=voidbuff; i < n; i++)
922 cbuff[i] = (byte)(scale*(buffer[i]-min));
923
924 }
925
926 /* convert double to byte */
927 else if (fs->bitpix == -64) {
928 double *buffer=voidbuff;
929 double max, min, maxmin_t, scale;
930
931 min = max = buffer[0];
932 for (i=1; i < n; i++, buffer++) maxmin(*buffer, max, min);
933 scale = (max == min) ? 0. : 255./(max-min);
934
935 /* rescale and convert */
936 for (i=0, buffer=voidbuff; i < n; i++)
937 cbuff[i] = (byte)(scale*(buffer[i]-min));
938 }
939
940 free(voidbuff);
941 return n;
942 }
943
944 #undef maxmin
945
946
947 /************************************/
flip(buffer,nx,ny)948 static void flip(buffer, nx, ny)
949 byte *buffer;
950 int nx;
951 int ny;
952 {
953 /* reverse order of lines in image */
954
955 int i;
956 int j, v;
957 byte *buff1, *buff2;
958
959 for (i=0; i < ny/2; i++) {
960 buff1 = &buffer[i*nx];
961 buff2 = &buffer[(ny-1-i)*nx];
962 for (j=0; j < nx; j++) {
963 v = *buff1;
964 *(buff1++) = *buff2;
965 *(buff2++) = v;
966 }
967 }
968 }
969
970