1 /*
2 readimage.c
3
4 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 *
6 * Part of: SExtractor
7 *
8 * Author: E.BERTIN (IAP)
9 *
10 * Contents: functions for input of image data.
11 *
12 * Last modify: 13/07/2006
13 *
14 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15 */
16
17 #ifdef HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20
21 #include <math.h>
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <string.h>
25
26 #include "wcs/wcs.h"
27 #include "define.h"
28 #include "globals.h"
29 #include "prefs.h"
30 #include "check.h"
31 #include "field.h"
32 #include "fits/fitscat.h"
33 #include "interpolate.h"
34 #include "back.h"
35 #include "astrom.h"
36 #include "weight.h"
37 #include "wcs/tnx.h"
38
39 /******************************* loadstrip ***********************************/
40 /*
41 Load a new strip of pixel data into the buffer.
42 */
loadstrip(picstruct * field,picstruct * wfield)43 void *loadstrip(picstruct *field, picstruct *wfield)
44
45 {
46 checkstruct *check;
47 int y, w, flags, interpflag;
48 PIXTYPE *data, *wdata, *rmsdata;
49
50 w = field->width;
51 flags = field->flags;
52 interpflag = (wfield && wfield->interp_flag);
53 wdata = NULL; /* To avoid gcc -Wall warnings */
54
55 if (!field->y)
56 {
57 /*- First strip */
58 int nbpix;
59
60 nbpix = w*field->stripheight;
61
62 if (flags ^ FLAG_FIELD)
63 {
64 /*---- Allocate space for the frame-buffer */
65 if (!(field->strip=(PIXTYPE *)malloc(field->stripheight*field->width
66 *sizeof(PIXTYPE))))
67 error(EXIT_FAILURE,"Not enough memory for the image buffer of ",
68 field->rfilename);
69
70 data = field->strip;
71 /*---- We assume weight data have been read just before */
72 if (interpflag)
73 wdata = wfield->strip;
74 if (flags & BACKRMS_FIELD)
75 for (y=0, rmsdata=data; y<field->stripheight; y++, rmsdata += w)
76 backrmsline(field, y, rmsdata);
77 else if (flags & INTERP_FIELD)
78 copydata(field, 0, nbpix);
79 else
80 readdata(field, data, nbpix);
81 if (flags & (WEIGHT_FIELD|RMS_FIELD|BACKRMS_FIELD|VAR_FIELD))
82 weight_to_var(field, data, nbpix);
83 if ((flags & MEASURE_FIELD) && (check=prefs.check[CHECK_IDENTICAL]))
84 writecheck(check, data, nbpix);
85 for (y=0; y<field->stripheight; y++, data += w)
86 {
87 /*------ This is the only place where one can pick-up safely the current bkg */
88 if (flags & (MEASURE_FIELD|DETECT_FIELD))
89 subbackline(field, y, data);
90 /*------ Go to interpolation process */
91 if (interpflag)
92 {
93 interpolate(field,wfield, data, wdata);
94 wdata += w;
95 }
96 /*------ Check-image stuff */
97 if (prefs.check_flag)
98 {
99 if (flags & MEASURE_FIELD)
100 {
101 if ((check = prefs.check[CHECK_BACKGROUND]))
102 writecheck(check, field->backline, w);
103 if ((check = prefs.check[CHECK_SUBTRACTED]))
104 writecheck(check, data, w);
105 if ((check = prefs.check[CHECK_APERTURES]))
106 writecheck(check, data, w);
107 if ((check = prefs.check[CHECK_SUBPSFPROTOS]))
108 writecheck(check, data, w);
109 if ((check = prefs.check[CHECK_SUBPCPROTOS]))
110 writecheck(check, data, w);
111 }
112 if ((flags&DETECT_FIELD) && (check=prefs.check[CHECK_BACKRMS]))
113 {
114 backrmsline(field, y, (PIXTYPE *)check->pix);
115 writecheck(check, check->pix, w);
116 }
117 }
118 }
119 }
120 else
121 {
122 if (!(field->fstrip=(FLAGTYPE *)malloc(field->stripheight*field->width
123 *sizeof(FLAGTYPE))))
124 error(EXIT_FAILURE,"Not enough memory for the flag buffer of ",
125 field->rfilename);
126 readidata(field, field->fstrip, nbpix);
127 }
128
129 field->ymax = field->stripheight;
130 if (field->ymax < field->height)
131 field->stripysclim = field->stripheight - field->stripmargin;
132 }
133 else
134 {
135 /*- other strips */
136 if (flags ^ FLAG_FIELD)
137 {
138 data = field->strip + field->stripylim*w;
139 /*---- We assume weight data have been read just before */
140 if (interpflag)
141 wdata = wfield->strip + field->stripylim*w;
142
143 /*---- copy to Check-image the "oldest" line before it is replaced */
144 if ((flags & MEASURE_FIELD) && (check=prefs.check[CHECK_SUBOBJECTS]))
145 writecheck(check, data, w);
146
147 if (flags & BACKRMS_FIELD)
148 backrmsline(field, field->ymax, data);
149 else if (flags & INTERP_FIELD)
150 copydata(field, field->stripylim*w, w);
151 else
152 readdata(field, data, w);
153 if (flags & (WEIGHT_FIELD|RMS_FIELD|BACKRMS_FIELD|VAR_FIELD))
154 weight_to_var(field, data, w);
155
156 if ((flags & MEASURE_FIELD) && (check=prefs.check[CHECK_IDENTICAL]))
157 writecheck(check, data, w);
158 /*---- Interpolate and subtract the background at current line */
159 if (flags & (MEASURE_FIELD|DETECT_FIELD))
160 subbackline(field, field->ymax, data);
161 if (interpflag)
162 interpolate(field,wfield, data, wdata);
163 /*---- Check-image stuff */
164 if (prefs.check_flag)
165 {
166 if (flags & MEASURE_FIELD)
167 {
168 if ((check = prefs.check[CHECK_BACKGROUND]))
169 writecheck(check, field->backline, w);
170 if ((check = prefs.check[CHECK_SUBTRACTED]))
171 writecheck(check, data, w);
172 if ((check = prefs.check[CHECK_APERTURES]))
173 writecheck(check, data, w);
174 if ((check = prefs.check[CHECK_SUBPSFPROTOS]))
175 writecheck(check, data, w);
176 if ((check = prefs.check[CHECK_SUBPCPROTOS]))
177 writecheck(check, data, w);
178 }
179 if ((flags&DETECT_FIELD) && (check=prefs.check[CHECK_BACKRMS]))
180 {
181 backrmsline(field, field->ymax, (PIXTYPE *)check->pix);
182 writecheck(check, check->pix, w);
183 }
184 }
185 }
186 else
187 readidata(field, field->fstrip + field->stripylim*w, w);
188
189 field->stripylim = (++field->ymin)%field->stripheight;
190 if ((++field->ymax)<field->height)
191 field->stripysclim = (++field->stripysclim)%field->stripheight;
192 }
193
194 return (flags ^ FLAG_FIELD)?
195 (void *)(field->strip + field->stripy*w)
196 : (void *)(field->fstrip + field->stripy*w);
197 }
198
199
200 /******************************** copydata **********************************/
201 /*
202 Copy image data from one field to the other.
203 */
copydata(picstruct * field,int offset,int size)204 void copydata(picstruct *field, int offset, int size)
205 {
206 memcpy(field->strip+offset, field->reffield->strip+offset,
207 size*sizeof(PIXTYPE));
208 return;
209 }
210
211
212 /******************************** readdata **********************************/
213 /*
214 read and convert input data stream in PIXTYPE (float) format.
215 */
readdata(picstruct * field,PIXTYPE * ptr,int size)216 void readdata(picstruct *field, PIXTYPE *ptr, int size)
217 {
218 static char bufdata0[DATA_BUFSIZE];
219 char *bufdata;
220 short val16;
221 int i, bowl, spoonful, npix, curval, dval;
222 PIXTYPE bs,bz;
223
224 bs = (PIXTYPE)field->bscale;
225 bz = (PIXTYPE)field->bzero;
226 switch(field->compress_type)
227 {
228 /*-- Uncompressed image */
229 case ICOMPRESS_NONE:
230 bowl = DATA_BUFSIZE/field->bytepix;
231 spoonful = size<bowl?size:bowl;
232 for(; size>0; size -= spoonful)
233 {
234 if (spoonful>size)
235 spoonful = size;
236 bufdata = bufdata0;
237 QFREAD(bufdata, spoonful*field->bytepix, field->file,field->filename);
238 switch(field->bitpix)
239 {
240 case BP_BYTE:
241 if (field->bitsgn)
242 for (i=spoonful; i--;)
243 *(ptr++) = *(bufdata++)*bs + bz;
244 else
245 for (i=spoonful; i--;)
246 *(ptr++) = *((unsigned char *)bufdata++)*bs + bz;
247 break;
248
249 case BP_SHORT:
250 if (bswapflag)
251 swapbytes(bufdata, 2, spoonful);
252 if (field->bitsgn)
253 for (i=spoonful; i--; bufdata += sizeof(short))
254 *(ptr++) = *((short *)bufdata)*bs + bz;
255 else
256 for (i=spoonful; i--; bufdata += sizeof(USHORT))
257 *(ptr++) = *((USHORT *)bufdata)*bs + bz;
258 break;
259
260 case BP_LONG:
261 if (bswapflag)
262 swapbytes(bufdata, 4, spoonful);
263 if (field->bitsgn)
264 for (i=spoonful; i--; bufdata += sizeof(LONG))
265 *(ptr++) = *((LONG *)bufdata)*bs + bz;
266 else
267 for (i=spoonful; i--; bufdata += sizeof(ULONG))
268 *(ptr++) = *((ULONG *)bufdata)*bs + bz;
269 break;
270
271 case BP_FLOAT:
272 if (bswapflag)
273 swapbytes(bufdata, 4, spoonful);
274 for (i=spoonful; i--; bufdata += sizeof(float))
275 *(ptr++) = *((float *)bufdata)*bs + bz;
276 break;
277
278 case BP_DOUBLE:
279 if (bswapflag)
280 swapbytes(bufdata, 8, spoonful);
281 for (i=spoonful; i--; bufdata += sizeof(double))
282 *(ptr++) = *((double *)bufdata)*bs + bz;
283 break;
284
285 default:
286 error(EXIT_FAILURE,"*FATAL ERROR*: unknown BITPIX type in ",
287 "readdata()");
288 break;
289 }
290 }
291 break;
292
293 /*-- Compressed image */
294 case ICOMPRESS_BASEBYTE:
295 bufdata = field->compress_bufptr;
296 curval = field->compress_curval;
297 npix = field->compress_npix;
298 while (size--)
299 {
300 if (!(npix--))
301 {
302 if (curval != field->compress_checkval)
303 error(EXIT_FAILURE, "*Error*: invalid BASEBYTE checksum in ",
304 field->filename);
305 bufdata = field->compress_buf;
306 QFREAD(bufdata, FBSIZE, field->file, field->filename);
307 curval = 0;
308 if (bswapflag)
309 swapbytes(bufdata, 4, 1);
310 field->compress_checkval = *((int *)bufdata);
311 bufdata += 4;
312 if (bswapflag)
313 swapbytes(bufdata, 2, 1);
314 npix = (int)(*((short *)bufdata))-1;
315 bufdata+=2;
316 }
317 if ((dval=(int)*(bufdata++))==-128)
318 {
319 if (bswapflag)
320 swapbytes(bufdata, 2, 1);
321 memcpy(&val16, bufdata, 2);
322 dval = (int)val16;
323 if (dval==-32768)
324 {
325 bufdata += 2;
326 if (bswapflag)
327 swapbytes(bufdata, 4, 1);
328 memcpy(&dval,bufdata,4);
329 bufdata += 4;
330 }
331 else
332 bufdata += 2;
333 }
334 *(ptr++) = dval*bs + bz;
335 curval += dval;
336 }
337 field->compress_curval = curval;
338 field->compress_bufptr = bufdata;
339 field->compress_npix = npix;
340 break;
341
342 case ICOMPRESS_PREVPIX:
343 bufdata = field->compress_bufptr;
344 curval = field->compress_curval;
345 npix = field->compress_npix;
346 while (size--)
347 {
348 if (!(npix--))
349 {
350 if (curval != field->compress_checkval)
351 error(EXIT_FAILURE, "*Error*: invalid PREV_PIX checksum in ",
352 field->filename);
353 bufdata = field->compress_buf;
354 QFREAD(bufdata, FBSIZE, field->file, field->filename);
355 if (bswapflag)
356 swapbytes(bufdata, 2, 3);
357 curval = (int)*(short *)bufdata;
358 npix = (int)*(short *)(bufdata+=2)-1;
359 field->compress_checkval = (int)(*(short *)(bufdata+=2));
360 bufdata+=4;
361 }
362 if ((dval=(int)*(bufdata++))==-128)
363 {
364 if (bswapflag)
365 swapbytes(bufdata, 2, 1);
366 memcpy(&val16, bufdata, 2);
367 curval = (int)val16;
368 bufdata += 2;
369 }
370 else
371 curval += dval;
372 *(ptr++) = curval*bs + bz;
373 }
374 field->compress_curval = curval;
375 field->compress_bufptr = bufdata;
376 field->compress_npix = npix;
377 break;
378
379 default:
380 error(EXIT_FAILURE,"*Internal Error*: unknown compression mode in ",
381 "readdata()");
382 }
383
384 return;
385 }
386
387
388 /******************************** readidata *********************************/
389 /*
390 read and convert input data stream in FLAGTYPE (unsigned int) format.
391 */
readidata(picstruct * field,FLAGTYPE * ptr,int size)392 void readidata(picstruct *field, FLAGTYPE *ptr, int size)
393 {
394 static char bufdata0[DATA_BUFSIZE];
395 char *bufdata;
396 short val16;
397 int i, bowl, spoonful, npix, curval, dval;
398
399 switch(field->compress_type)
400 {
401 /*-- Uncompressed image */
402 case ICOMPRESS_NONE:
403 bowl = DATA_BUFSIZE/field->bytepix;
404 spoonful = size<bowl?size:bowl;
405 for(; size>0; size -= spoonful)
406 {
407 if (spoonful>size)
408 spoonful = size;
409 bufdata = bufdata0;
410 QFREAD(bufdata, spoonful*field->bytepix, field->file, field->filename);
411 switch(field->bitpix)
412 {
413 case BP_BYTE:
414 for (i=spoonful; i--;)
415 *(ptr++) = (FLAGTYPE)*((unsigned char *)bufdata++);
416 break;
417
418 case BP_SHORT:
419 if (bswapflag)
420 swapbytes(bufdata, 2, spoonful);
421 for (i=spoonful; i--; bufdata += sizeof(USHORT))
422 *(ptr++) = (FLAGTYPE)*((USHORT *)bufdata);
423 break;
424
425 case BP_LONG:
426 if (bswapflag)
427 swapbytes(bufdata, 4, spoonful);
428 for (i=spoonful; i--; bufdata += sizeof(ULONG))
429 *(ptr++) = (FLAGTYPE)*((ULONG *)bufdata);
430 break;
431
432 case BP_FLOAT:
433 case BP_DOUBLE:
434 error(EXIT_FAILURE,"*Error*: I was expecting integers in ",
435 field->filename);
436 break;
437 default:
438 error(EXIT_FAILURE,"*FATAL ERROR*: unknown BITPIX type in ",
439 "readdata()");
440 break;
441 }
442 }
443 break;
444
445 /*-- Compressed image */
446 case ICOMPRESS_BASEBYTE:
447 bufdata = field->compress_bufptr;
448 curval = field->compress_curval;
449 npix = field->compress_npix;
450 while (size--)
451 {
452 if (!(npix--))
453 {
454 if (curval != field->compress_checkval)
455 error(EXIT_FAILURE, "*Error*: invalid BASEBYTE checksum in ",
456 field->filename);
457 bufdata = field->compress_buf;
458 QFREAD(bufdata, FBSIZE, field->file, field->filename);
459 curval = 0;
460 if (bswapflag)
461 swapbytes(bufdata, 4, 1);
462 field->compress_checkval = *((int *)bufdata);
463 bufdata += 4;
464 if (bswapflag)
465 swapbytes(bufdata, 2, 1);
466 npix = (int)(*((short *)bufdata))-1;
467 bufdata+=2;
468 }
469 if ((dval=(int)*(bufdata++))==-128)
470 {
471 if (bswapflag)
472 swapbytes(bufdata, 2, 1);
473 memcpy(&val16, bufdata, 2);
474 dval = (int)val16;
475 if (dval==-32768)
476 {
477 bufdata += 2;
478 if (bswapflag)
479 swapbytes(bufdata, 4, 1);
480 memcpy(&dval,bufdata,4);
481 bufdata += 4;
482 }
483 else
484 bufdata += 2;
485 }
486 *(ptr++) = (FLAGTYPE)dval;
487 curval += dval;
488 }
489 field->compress_curval = curval;
490 field->compress_bufptr = bufdata;
491 field->compress_npix = npix;
492 break;
493
494 case ICOMPRESS_PREVPIX:
495 bufdata = field->compress_bufptr;
496 curval = field->compress_curval;
497 npix = field->compress_npix;
498 while (size--)
499 {
500 if (!(npix--))
501 {
502 if (curval != field->compress_checkval)
503 error(EXIT_FAILURE, "*Error*: invalid PREV_PIX checksum in ",
504 field->filename);
505 bufdata = field->compress_buf;
506 QFREAD(bufdata, FBSIZE, field->file, field->filename);
507 if (bswapflag)
508 swapbytes(bufdata, 2, 3);
509 curval = (int)*(short *)bufdata;
510 npix = (int)*(short *)(bufdata+=2)-1;
511 field->compress_checkval = (int)(*(short *)(bufdata+=2));
512 bufdata+=4;
513 }
514 if ((dval=(int)*(bufdata++))==-128)
515 {
516 if (bswapflag)
517 swapbytes(bufdata, 2, 1);
518 memcpy(&val16, bufdata, 2);
519 curval = (int)val16;
520 bufdata += 2;
521 }
522 else
523 curval += dval;
524 *(ptr++) = (FLAGTYPE)curval;
525 }
526 field->compress_curval = curval;
527 field->compress_bufptr = bufdata;
528 field->compress_npix = npix;
529 break;
530
531 default:
532 error(EXIT_FAILURE,"*Internal Error*: unknown compression mode in ",
533 "readdata()");
534 }
535
536 return;
537 }
538
539
540 /******************************* readimagehead *******************************/
541 /*
542 extract some data from the FITS-file header
543 */
readimagehead(picstruct * field)544 void readimagehead(picstruct *field)
545 {
546 int j,l, n;
547 char wstr1[TNX_MAXCHARS], wstr2[TNX_MAXCHARS],
548 st[80], str[80],
549 *buf, *point;
550
551 /* Open the file */
552 if (!(field->file = fopen(field->filename, "rb")))
553 error(EXIT_FAILURE,"*Error*: cannot open ", field->filename);
554 /* Go directly to the right extension */
555 if (field->mefpos)
556 {
557 QFSEEK(field->file, field->mefpos, SEEK_SET, field->filename);
558 }
559 buf = readfitshead(field->file, field->filename, &n);
560 if(FITSTOI("NAXIS ", 0) < 2)
561 error(EXIT_FAILURE, field->filename, " does NOT contain 2D-data!");
562
563 /*---------------------------- Basic keywords ------------------------------*/
564 field->bitpix = FITSTOI("BITPIX ", 0);
565 if (field->bitpix != BP_BYTE
566 && field->bitpix != BP_SHORT
567 && field->bitpix != BP_LONG
568 && field->bitpix != BP_FLOAT
569 && field->bitpix != BP_DOUBLE)
570 error(EXIT_FAILURE, "Sorry, I don't know that kind of data.", "");
571
572 field->bytepix = (field->bitpix>0?field->bitpix:-field->bitpix)>>3;
573 field->width = FITSTOI("NAXIS1 ", 0);
574 field->height = FITSTOI("NAXIS2 ", 0);
575 field->npix = (KINGSIZE_T)field->width*field->height;
576 field->bscale = FITSTOF("BSCALE ", 1.0);
577
578 field->bzero = FITSTOF("BZERO ", 0.0);
579 field->bitsgn = FITSTOI("BITSGN ", 1);
580 if (field->bitsgn && prefs.fitsunsigned_flag)
581 field->bitsgn = 0;
582
583 FITSTOS("OBJECT ", field->ident, "Unnamed");
584
585 /*----------------------------- Compression --------------------------------*/
586 if (fitsread(buf, "IMAGECOD", st, H_STRING, T_STRING)==RETURN_OK)
587 {
588 if (!strcmp(st, "NONE"))
589 field->compress_type = ICOMPRESS_NONE;
590 else if (!strcmp(st, "BASEBYTE"))
591 field->compress_type = ICOMPRESS_BASEBYTE;
592 else if (!strcmp(st, "PREV_PIX"))
593 field->compress_type = ICOMPRESS_PREVPIX;
594 else
595 warning("Compression skipped: unknown IMAGECOD parameter:", st);
596 }
597
598 /*----------------------------- Astrometry ---------------------------------*/
599 /* Presently, astrometry is done only on the measurement and detect images */
600 if (field->flags&(MEASURE_FIELD|DETECT_FIELD))
601 {
602 astromstruct *as;
603 double drota, s;
604 int naxis;
605
606 QCALLOC(as, astromstruct, 1);
607 field->astrom = as;
608
609 naxis = as->naxis = 2;
610 for (l=0; l<naxis; l++)
611 {
612 sprintf(str, "CTYPE%-3d", l+1);
613 FITSTOS(str, str, "");
614 strncpy(as->ctype[l], str, 8);
615 sprintf(str, "CUNIT%-3d", l+1);
616 FITSTOS(str, str, "deg");
617 strncpy(as->cunit[l], str, 32);
618 sprintf(str, "CRVAL%-3d", l+1);
619 as->crval[l] = FITSTOF(str, 0.0);
620 sprintf(str, "CRPIX%-3d", l+1);
621 as->crpix[l] = FITSTOF(str, 1.0);
622 sprintf(str, "CDELT%-3d", l+1);
623 as->cdelt[l] = FITSTOF(str, 1.0);
624 if (fabs(as->cdelt[l]) < 1/BIG)
625 error(EXIT_FAILURE, "*Error*: CDELT parameters out of range in ",
626 field->filename);
627 }
628 if (fitsnfind(buf, "CD1_1", n))
629 {
630 /*---- If CD keywords exist, use them for the linear mapping terms... */
631 for (l=0; l<naxis; l++)
632 for (j=0; j<naxis; j++)
633 {
634 sprintf(str, "CD%d_%d", l+1, j+1);
635 as->pc[l*naxis+j] = FITSTOF(str, l==j?1.0:0.0)/as->cdelt[l];
636 }
637 }
638 else if (fitsnfind(buf, "PC001001", n))
639 /*---- ...If PC keywords exist, use them for the linear mapping terms... */
640 for (l=0; l<naxis; l++)
641 for (j=0; j<naxis; j++)
642 {
643 sprintf(str, "PC%03d%03d", l+1, j+1);
644 as->pc[l*naxis+j] = FITSTOF(str, l==j?1.0:0.0);
645 }
646 else
647 {
648 /*---- ...otherwise take the obsolete CROTA2 parameter */
649 s = as->cdelt[1]/as->cdelt[0];
650 drota = FITSTOF("CROTA2 ", 0.0);
651 as->pc[3] = as->pc[0] = cos(drota*DEG);
652 as->pc[1] = -(as->pc[2] = sin(drota*DEG));
653 as->pc[1] *= s;
654 as->pc[2] /= s;
655 }
656
657 QMALLOC(as->wcs, struct wcsprm, 1);
658 /*-- Test if the WCS is recognized and a celestial pair is found */
659 if (prefs.world_flag
660 && !wcsset(as->naxis,(const char(*)[9])as->ctype, as->wcs)
661 && as->wcs->flag<999)
662 {
663 char *pstr;
664 double date;
665 int biss, dpar[3];
666
667 as->wcs_flag = 1;
668 /*---- Coordinate reference frame */
669 /*---- Search for an observation date expressed in Julian days */
670 date = FITSTOF("MJD-OBS ", -1.0);
671 /*---- Precession date (defined from Ephemerides du Bureau des Longitudes) */
672 /*---- in Julian years from 2000.0 */
673 if (date>0.0)
674 as->equinox = 2000.0 - (MJD2000 - date)/365.25;
675 else
676 {
677 /*------ Search for an observation date expressed in "civil" format */
678 FITSTOS("DATE-OBS", str, "");
679 if (*str)
680 {
681 /*-------- Decode DATE-OBS format: DD/MM/YY or YYYY-MM-DD */
682 for (l=0; l<3 && (pstr = strtok(l?NULL:str,"/- ")); l++)
683 dpar[l] = atoi(pstr);
684 if (l<3 || !dpar[0] || !dpar[1] || !dpar[2])
685 {
686 /*---------- If DATE-OBS value corrupted or incomplete, assume 2000-1-1 */
687 warning("Invalid DATE-OBS value in header: ", str);
688 dpar[0] = 2000; dpar[1] = 1; dpar[2] = 1;
689 }
690 else if (strchr(str, '/') && dpar[0]<32 && dpar[2]<100)
691 {
692 j = dpar[0];
693 dpar[0] = dpar[2]+1900;
694 dpar[2] = j;
695 }
696
697 biss = (dpar[0]%4)?0:1;
698 /*-------- Convert date to MJD */
699 date = -678956 + (365*dpar[0]+dpar[0]/4) - biss
700 + ((dpar[1]>2?((int)((dpar[1]+1)*30.6)-63+biss)
701 :((dpar[1]-1)*(63+biss))/2) + dpar[2]);
702 as->equinox = 2000.0 - (MJD2000 - date)/365.25;
703 }
704 else
705 /*-------- Well if really no date is found */
706 as->equinox = 2000.0;
707 }
708 if (field->flags&MEASURE_FIELD)
709 prefs.epoch = as->equinox;
710 FITSTOS("RADECSYS", str, as->equinox<1984.0?
711 "FK4": (as->equinox<1999.9999? "FK5" : "ICRS"));
712 if (!strcmp(str, "ICRS"))
713 {
714 as->radecsys = RDSYS_ICRS;
715 as->equinox = FITSTOF("EQUINOX ", 2000.0);
716 }
717 else if (!strcmp(str, "FK5"))
718 {
719 as->radecsys = RDSYS_FK5;
720 as->equinox = FITSTOF("EQUINOX ", FITSTOF("EPOCH ", 2000.0));
721 if (field->flags&MEASURE_FIELD)
722 sprintf(prefs.coosys, "eq_FK5");
723 }
724 else if (!strcmp(str, "FK4"))
725 {
726 if (as->equinox == 2000.0)
727 as->equinox = FITSTOF("EQUINOX ", FITSTOF("EPOCH ", 1950.0));
728 as->radecsys = RDSYS_FK4;
729 warning("FK4 precession formulae not yet implemented:\n",
730 " Astrometry may be slightly inaccurate");
731 }
732 else if (!strcmp(str, "FK4-NO-E"))
733 {
734 if (as->equinox == 2000.0)
735 as->equinox = FITSTOF("EQUINOX ", FITSTOF("EPOCH ", 1950.0));
736 as->radecsys = RDSYS_FK4_NO_E;
737 warning("FK4 precession formulae not yet implemented:\n",
738 " Astrometry may be slightly inaccurate");
739 }
740 else if (!strcmp(str, "GAPPT"))
741 {
742 as->radecsys = RDSYS_GAPPT;
743 warning("GAPPT reference frame not yet implemented:\n",
744 " Astrometry may be slightly inaccurate");
745 }
746 else
747 {
748 warning("Using ICRS instead of unknown astrometric reference frame: ",
749 str);
750 as->radecsys = RDSYS_ICRS;
751 as->equinox = FITSTOF("EQUINOX ", 2000.0);
752 }
753
754 /*---- Projection parameters */
755 if (!strcmp(as->wcs->pcode, "TNX"))
756 {
757 /*---- IRAF's TNX projection: decode these #$!?@#!! WAT parameters */
758 if (fitsfind(buf, "WAT?????") != RETURN_ERROR)
759 {
760 /*-------- First we need to concatenate strings */
761 pstr = wstr1;
762 sprintf(str, "WAT1_001");
763 for (j=2; fitsread(buf,str,pstr,H_STRINGS,T_STRING)==RETURN_OK; j++)
764 {
765 sprintf(str, "WAT1_%03d", j);
766 pstr += strlen(pstr);
767 }
768 pstr = wstr2;
769 sprintf(str, "WAT2_001");
770 for (j=2; fitsread(buf,str,pstr,H_STRINGS,T_STRING)==RETURN_OK; j++)
771 {
772 sprintf(str, "WAT2_%03d", j);
773 pstr += strlen(pstr);
774 }
775 /*-------- LONGPOLE defaulted to 180 deg if not found */
776 if ((pstr = strstr(wstr1, "longpole"))
777 || (pstr = strstr(wstr2, "longpole")))
778 pstr = strpbrk(pstr, "1234567890-+.");
779 as->longpole = pstr? atof(pstr) : 999.0;
780 as->latpole = 999.0;
781 /*-------- RO defaulted to 180/PI if not found */
782 if ((pstr = strstr(wstr1, "ro"))
783 || (pstr = strstr(wstr2, "ro")))
784 pstr = strpbrk(pstr, "1234567890-+.");
785 as->r0 = pstr? atof(pstr) : 0.0;
786 /*-------- Read the remaining TNX parameters */
787 if ((pstr = strstr(wstr1, "lngcor"))
788 || (pstr = strstr(wstr2, "lngcor")))
789 as->tnx_lngcor = read_tnxaxis(pstr);
790 if (!as->tnx_lngcor)
791 error(EXIT_FAILURE, "*Error*: incorrect TNX parameters in ",
792 field->filename);
793 if ((pstr = strstr(wstr1, "latcor"))
794 || (pstr = strstr(wstr2, "latcor")))
795 as->tnx_latcor = read_tnxaxis(pstr);
796 if (!as->tnx_latcor)
797 error(EXIT_FAILURE, "*Error*: incorrect TNX parameters in ",
798 field->filename);
799 }
800 }
801 else
802 {
803 as->longpole = FITSTOF("LONGPOLE", 999.0);
804 as->latpole = FITSTOF("LATPOLE ", 999.0);
805 if (fitsnfind(buf, "PROJP1 ", n))
806 for (l=0; l<10; l++)
807 {
808 sprintf(str, "PROJP%-3d", l);
809 as->projp[l] = FITSTOF(str, 0.0);
810 }
811 }
812 }
813 else
814 {
815 /*---- No need to keep memory allocated for a useless WCS structure */
816 free(as->wcs);
817 as->wcs_flag = 0;
818 }
819 }
820
821 /*---------------------------------------------------------------------------*/
822
823 field->fitshead = buf;
824 field->fitsheadsize = n*FBSIZE;
825
826 return;
827 }
828
829
830 /******************************* readfitshead ********************************/
831 /*
832 read data from the FITS-file header
833 */
readfitshead(FILE * file,char * filename,int * nblock)834 char *readfitshead(FILE *file, char *filename, int *nblock)
835
836 {
837 int n;
838 char *buf;
839
840 if (!(buf=(char *)malloc((size_t)FBSIZE)))
841 error(EXIT_FAILURE, "*Error*: Not enough memory in ", "readfitshead()");
842
843 /* Find the number of FITS blocks of the header while reading it */
844 QFREAD(buf, FBSIZE, file, filename);
845
846 if (strncmp(buf, "SIMPLE ", 8))
847 {
848 /* Ugly but necessary patch to handle this stupid DeNIS compressed format! */
849 if (strncmp(buf, "XTENSION", 8))
850 error(EXIT_FAILURE, filename, " is NOT a FITS file!");
851 /*
852 else
853 {
854 memset(buf, ' ', 80);
855 strncpy(buf,
856 "SIMPLE = T / Decompressed by SExtractor", 59);
857 }
858 */
859 }
860
861 for (n=1; !fitsnfind(buf,"END ", n); n++)
862 {
863 if (!(buf=(char *)realloc(buf, (size_t)(FBSIZE*(n+1)))))
864 error(EXIT_FAILURE, "*Error*: Not enough memory in ", "readfitshead()");
865 QFREAD(buf+n*FBSIZE, FBSIZE, file, filename);
866 }
867
868 *nblock = n;
869 return buf;
870 }
871
872
873