IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ticket #349: psFits.c

File psFits.c, 50.0 KB (added by robert.desonia@…, 21 years ago)

replacement src/fileUtils/psFits.c

Line 
1/** @file psFits.c
2 *
3 * @brief Contains Fits I/O routines
4 *
5 * @ingroup FileIO
6 *
7 * @author Robert DeSonia, MHPCC
8 *
9 * @version $Revision: 1.25 $ $Name: $
10 * @date $Date: 2005/04/05 19:33:30 $
11 *
12 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
13 */
14
15#include <unistd.h>
16
17#include "psFits.h"
18#include "string.h"
19#include "psError.h"
20#include "psFileUtilsErrors.h"
21#include "psImageExtraction.h"
22#include "psMemory.h"
23#include "psString.h"
24#include "psLogMsg.h"
25#include "psTrace.h"
26
27#define MAX_STRING_LENGTH 256 // maximum length string for FITS routines
28
29static char* standardFitsKeys[] = {
30 NULL
31 };
32
33static psElemType convertFitsToPsType(int datatype)
34{
35 switch (datatype) {
36 case TBYTE:
37 return PS_TYPE_U8;
38 case TSBYTE:
39 return PS_TYPE_S8;
40 case TSHORT:
41 return PS_TYPE_S16;
42 case TUSHORT:
43 return PS_TYPE_U16;
44 case TLONG:
45 if (sizeof(long) == 8) {
46 return PS_TYPE_S64;
47 }
48 // no break
49 case TINT:
50 return PS_TYPE_S32;
51 case TULONG:
52 if (sizeof(unsigned long) == 8) {
53 return PS_TYPE_U64;
54 }
55 // no break
56 case TUINT:
57 return PS_TYPE_U32;
58 case TLONGLONG:
59 return PS_TYPE_S64;
60 case TFLOAT:
61 return PS_TYPE_F32;
62 case TDOUBLE:
63 return PS_TYPE_F64;
64 case TCOMPLEX:
65 return PS_TYPE_C32;
66 case TDBLCOMPLEX:
67 return PS_TYPE_C64;
68 case TLOGICAL:
69 return PS_TYPE_BOOL;
70 default:
71 psError(PS_ERR_IO, true,
72 "Unknown FITS datatype, %d.",
73 datatype);
74 return PS_TYPE_PTR;
75 }
76}
77
78static bool convertPsTypeToFits(psElemType type, int* bitPix, double* bZero, int* dataType)
79{
80
81 int bitpix;
82 int datatype;
83 double bzero = 0.0;
84
85 switch (type) {
86
87 case PS_TYPE_U8:
88 bitpix = BYTE_IMG;
89 datatype = TBYTE;
90 break;
91
92 case PS_TYPE_BOOL:
93 case PS_TYPE_S8:
94 bitpix = BYTE_IMG;
95 bzero = INT8_MIN;
96 datatype = TSBYTE;
97 break;
98
99 case PS_TYPE_U16:
100 bitpix = SHORT_IMG;
101 bzero = -1.0 * INT16_MIN;
102 datatype = TUSHORT;
103 break;
104
105 case PS_TYPE_S16:
106 bitpix = SHORT_IMG;
107 datatype = TSHORT;
108 break;
109
110 case PS_TYPE_U32:
111 bitpix = LONG_IMG;
112 bzero = -1.0 * INT32_MIN;
113 datatype = TUINT;
114 break;
115
116 case PS_TYPE_S32:
117 bitpix = LONG_IMG;
118 datatype = TINT;
119 break;
120
121 case PS_TYPE_F32:
122 bitpix = FLOAT_IMG;
123 datatype = TFLOAT;
124 break;
125
126 case PS_TYPE_F64:
127 bitpix = DOUBLE_IMG;
128 datatype = TDOUBLE;
129 break;
130
131 default: {
132 char* typeStr;
133 PS_TYPE_NAME(typeStr,type);
134 psError(PS_ERR_BAD_PARAMETER_TYPE, true,
135 PS_ERRORTEXT_psFits_TYPE_UNSUPPORTED,
136 typeStr);
137 return false;
138 }
139 }
140
141 // pass the requested parameters (NULL parameters are not set, of course).
142 if (bitPix != NULL) {
143 *bitPix = bitpix;
144 }
145
146 if (dataType != NULL) {
147 *dataType = datatype;
148 }
149
150 if (bZero != NULL) {
151 *bZero = bzero;
152 }
153
154 return true;
155}
156
157static bool convertMetadataTypeToBinaryTForm(psMetadataType type, char** fitsType)
158{
159 switch (type) {
160 case PS_META_BOOL:
161 *fitsType = psStringCopy("1L");
162 break;
163 case PS_META_S32:
164 *fitsType = psStringCopy("1J");
165 break;
166 case PS_META_F32:
167 *fitsType = psStringCopy("1E");
168 break;
169 case PS_META_F64:
170 *fitsType = psStringCopy("1D");
171 break;
172 // XXX: Handle other types, e.g., Vectors, etc.
173 default:
174 return false;
175 }
176
177 return true;
178}
179
180static bool isHDUEmpty(const psFits* fits)
181{
182 /* check for keys - no keys means this is really an empty HDU */
183 int keysexist = -1;
184 int morekeys;
185 int status = 0;
186
187 fits_get_hdrspace(fits->p_fd, &keysexist, &morekeys, &status);
188
189 // if no keys exist and not primary HDU, this really is an empty HDU
190 if (keysexist == 0) {
191 return true;
192 }
193
194 return false;
195
196}
197
198static void fitsFree(psFits* fits)
199{
200 int status = 0;
201
202 if (fits != NULL) {
203 (void)fits_close_file(fits->p_fd, &status);
204 psFree((void*)fits->filename);
205 }
206}
207
208psFits* psFitsAlloc(const char* name)
209{
210 int status = 0;
211 fitsfile *fptr = NULL; /* Pointer to the FITS file */
212
213 if (name == NULL) {
214 psError(PS_ERR_BAD_PARAMETER_NULL, true,
215 PS_ERRORTEXT_psFits_FILENAME_NULL);
216 return NULL;
217 }
218
219 /* Open/Create the FITS file */
220 if (access(name, F_OK) == 0) { // file exists
221 (void)fits_open_file(&fptr, name, READWRITE, &status);
222 if (fptr == NULL) { // if failed, try openning as just read-only
223 status = 0;
224 (void)fits_open_file(&fptr, name, READONLY, &status);
225 }
226 if (fptr == NULL || status != 0) {
227 char fitsErr[MAX_STRING_LENGTH];
228 fits_get_errstatus(status, fitsErr);
229 psError(PS_ERR_BAD_PARAMETER_VALUE, true,
230 PS_ERRORTEXT_psFits_FILENAME_INVALID,
231 name, fitsErr);
232 return NULL;
233 }
234 } else { // file does not exist, so create.
235 (void)fits_create_file(&fptr, name, &status);
236 if (fptr == NULL || status != 0) {
237 char fitsErr[MAX_STRING_LENGTH];
238 fits_get_errstatus(status, fitsErr);
239 psError(PS_ERR_IO, true,
240 PS_ERRORTEXT_psFits_FILENAME_CREATE_FAILED,
241 name, fitsErr);
242 return NULL;
243 }
244 }
245
246 psFits* fits = psAlloc(sizeof(psFits));
247 fits->filename = psAlloc(strlen(name)+1);
248 fits->p_fd = fptr;
249 strcpy((char*)fits->filename,name);
250 p_psMemSetDeallocator(fits,(psFreeFcn)fitsFree);
251
252 return fits;
253}
254
255bool psFitsMoveExtName(const psFits* fits,
256 const char* extname)
257{
258 int status = 0;
259
260 if (fits == NULL) {
261 psError(PS_ERR_BAD_PARAMETER_NULL, true,
262 PS_ERRORTEXT_psFits_NULL);
263 return false;
264 }
265
266 if (extname == NULL) {
267 psError(PS_ERR_BAD_PARAMETER_NULL, true,
268 PS_ERRORTEXT_psFits_EXTNAME_NULL);
269 return false;
270 }
271
272
273 if (fits_movnam_hdu(fits->p_fd, ANY_HDU, (char*)extname, 0, &status) != 0) {
274 char fitsErr[MAX_STRING_LENGTH];
275 fits_get_errstatus(status, fitsErr);
276 psError(PS_ERR_LOCATION_INVALID, true,
277 PS_ERRORTEXT_psFits_EXTNAME_INVALID,
278 extname, fits->filename, fitsErr);
279 return false;
280 }
281
282 return true;
283}
284
285bool psFitsMoveExtNum(const psFits* fits,
286 int extnum,
287 bool relative)
288{
289 if (fits == NULL) {
290 psError(PS_ERR_BAD_PARAMETER_NULL, true,
291 PS_ERRORTEXT_psFits_NULL);
292 return false;
293 }
294
295 int status = 0;
296 int hdutype = 0;
297
298 if (relative) {
299 fits_movrel_hdu(fits->p_fd, extnum, &hdutype, &status);
300 if (status != 0) {
301 char fitsErr[MAX_STRING_LENGTH];
302 fits_get_errstatus(status, fitsErr);
303 psError(PS_ERR_LOCATION_INVALID, true,
304 PS_ERRORTEXT_psFits_EXTNUM_REL_MOVE_FAILED,
305 extnum, fits->filename, fitsErr);
306 return false;
307 }
308 } else {
309 fits_movabs_hdu(fits->p_fd, extnum+1, &hdutype, &status);
310 if (status != 0) {
311 char fitsErr[MAX_STRING_LENGTH];
312 fits_get_errstatus(status, fitsErr);
313 psError(PS_ERR_LOCATION_INVALID, true,
314 PS_ERRORTEXT_psFits_EXTNUM_ABS_MOVE_FAILED,
315 extnum, fits->filename, fitsErr);
316 return false;
317 }
318 }
319
320 return true;
321}
322
323int psFitsGetExtNum(const psFits* fits)
324{
325 int hdunum;
326
327 if (fits == NULL) {
328 psError(PS_ERR_BAD_PARAMETER_NULL, true,
329 PS_ERRORTEXT_psFits_NULL);
330 return PS_FITS_TYPE_NONE;
331 }
332
333
334 return fits_get_hdu_num(fits->p_fd,&hdunum) - 1;
335}
336
337char* psFitsGetExtName(const psFits* fits)
338{
339 if (fits == NULL) {
340 psError(PS_ERR_BAD_PARAMETER_NULL, true,
341 PS_ERRORTEXT_psFits_NULL);
342 return NULL;
343 }
344
345 int status = 0;
346 char name[MAX_STRING_LENGTH];
347
348 if (fits_read_key_str(fits->p_fd, "EXTNAME", name, NULL, &status) != 0) {
349 status = 0;
350 if (fits_read_key_str(fits->p_fd, "HDUNAME", name, NULL, &status) != 0) {
351 int num = psFitsGetExtNum(fits);
352 snprintf(name, MAX_STRING_LENGTH, "EXT-%3d",num);
353 }
354 }
355 return psStringCopy(name);
356}
357
358bool psFitsSetExtName(const psFits* fits, const char* name)
359{
360 if (fits == NULL) {
361 psError(PS_ERR_BAD_PARAMETER_NULL, true,
362 PS_ERRORTEXT_psFits_NULL);
363 return false;
364 }
365
366 if (name == NULL) {
367 psError(PS_ERR_BAD_PARAMETER_NULL, true,
368 PS_ERRORTEXT_psFits_EXTNAME_NULL);
369 return false;
370 }
371
372 int status = 0;
373
374 if (fits_update_key_str(fits->p_fd, "EXTNAME", (char*)name, NULL, &status) != 0) {
375 char fitsErr[MAX_STRING_LENGTH];
376 (void)fits_get_errstatus(status, fitsErr);
377 psError(PS_ERR_IO, true,
378 PS_ERRORTEXT_psFits_WRITE_FAILED,
379 fits->filename, fitsErr);
380 return false;
381 }
382
383 return true;
384}
385
386int psFitsGetSize(const psFits* fits)
387{
388 if (fits == NULL) {
389 psError(PS_ERR_BAD_PARAMETER_NULL, true,
390 PS_ERRORTEXT_psFits_NULL);
391 return 0;
392 }
393
394 int num = 0;
395 int status = 0;
396
397 if (fits_get_num_hdus(fits->p_fd, &num, &status) != 0) {
398 char fitsErr[MAX_STRING_LENGTH];
399 fits_get_errstatus(status, fitsErr);
400 psError(PS_ERR_LOCATION_INVALID, true,
401 PS_ERRORTEXT_psFits_GETNUMHDUS_FAILED,
402 fits->filename, fitsErr);
403 return 0;
404 }
405
406 return num;
407}
408
409psFitsType psFitsGetExtType(const psFits* fits)
410{
411 if (fits == NULL) {
412 psError(PS_ERR_BAD_PARAMETER_NULL, true,
413 PS_ERRORTEXT_psFits_NULL);
414 return PS_FITS_TYPE_NONE;
415 }
416
417 int status = 0;
418 int hdutype = PS_FITS_TYPE_NONE;
419
420 if (fits_get_hdu_type(fits->p_fd, &hdutype, &status) != 0) {
421 char fitsErr[MAX_STRING_LENGTH];
422 fits_get_errstatus(status, fitsErr);
423 psError(PS_ERR_LOCATION_INVALID, true,
424 PS_ERRORTEXT_psFits_GETHDUTYPE_FAILED,
425 fits->filename, fitsErr);
426 return PS_FITS_TYPE_NONE;
427 }
428
429 if (hdutype == PS_FITS_TYPE_IMAGE &&
430 psFitsGetExtNum(fits) > 0 &&
431 isHDUEmpty(fits)) {
432 return PS_FITS_TYPE_ANY;
433 }
434
435 return hdutype;
436}
437
438psMetadata* psFitsReadHeader(psMetadata* out,
439 const psFits* fits)
440{
441 if (fits == NULL) {
442 psError(PS_ERR_BAD_PARAMETER_NULL, true,
443 PS_ERRORTEXT_psFits_NULL);
444 return NULL;
445 }
446
447 if (out == NULL) {
448 out = psMetadataAlloc();
449 if (out == NULL) {
450 psError(PS_ERR_UNKNOWN, false,
451 "Failed to allocate a new psMetadata container.");
452 return NULL;
453 }
454 }
455
456 // Get number of key names
457 int numKeys = 0;
458 int keyNum = 0;
459 int status = 0;
460 fits_get_hdrpos(fits->p_fd, &numKeys, &keyNum, &status);
461
462 // Get each key name. Keywords start at one.
463 char keyType;
464 char keyName[MAX_STRING_LENGTH];
465 char keyValue[MAX_STRING_LENGTH];
466 char keyComment[MAX_STRING_LENGTH];
467 psBool tempBool;
468 psBool success;
469 psBool stdKey;
470 for (int i = 1; i <= numKeys; i++) {
471
472 fits_read_keyn(fits->p_fd, i, keyName, keyValue, keyComment, &status);
473
474 stdKey = false;
475
476 int stdKeyIdx = 0;
477 while (standardFitsKeys[stdKeyIdx] != NULL && ! stdKey) {
478 if (strcmp(keyName,standardFitsKeys[stdKeyIdx++]) == 0) {
479 stdKey = true;
480 }
481 }
482
483 if (keyValue[0] != 0) { // blank values are not handled by fits_get_keytype
484 fits_get_keytype(keyValue, &keyType, &status);
485 } else {
486 keyType = 'C';
487 }
488 if (status != 0) {
489 break;
490 }
491
492 if (! stdKey) {
493 switch (keyType) {
494 case 'X': // bit
495 case 'I': // short int.
496 case 'J': // int.
497 case 'B': // byte
498 success = psMetadataAdd(out,
499 PS_LIST_TAIL,
500 keyName,
501 PS_META_S32 | PS_META_DUPLICATE_OK,
502 keyComment,
503 atoi(keyValue));
504 break;
505 case 'U': // unsigned int. may not fit in a psS32
506 case 'K': // long int. can't all fit in a psS32
507 case 'F':
508 success = psMetadataAdd(out,
509 PS_LIST_TAIL,
510 keyName,
511 PS_META_F64 | PS_META_DUPLICATE_OK,
512 keyComment,
513 atof(keyValue));
514 break;
515 case 'C':
516 // remove the single-quotes at front/end
517 if (keyValue[0] == '\'' && keyValue[strlen(keyValue)-1] == '\'') {
518 keyValue[strlen(keyValue)-1] = '\0';
519 success = psMetadataAdd(out,
520 PS_LIST_TAIL,
521 keyName,
522 PS_META_STR | PS_META_DUPLICATE_OK,
523 keyComment,
524 keyValue+1);
525 } else {
526 success = psMetadataAdd(out,
527 PS_LIST_TAIL,
528 keyName,
529 PS_META_STR | PS_META_DUPLICATE_OK,
530 keyComment,
531 keyValue);
532 }
533 break;
534 case 'L':
535 tempBool = (keyValue[0] == 'T') ? 1 : 0;
536 success = psMetadataAdd(out,
537 PS_LIST_TAIL,
538 keyName,
539 PS_META_BOOL | PS_META_DUPLICATE_OK,
540 keyComment,
541 tempBool);
542 break;
543 default:
544 psError(PS_ERR_IO, true,
545 PS_ERRORTEXT_psFits_METATYPE_INVALID,
546 keyType);
547 return out;
548 }
549
550 if (!success) {
551 psError(PS_ERR_UNKNOWN, false,
552 PS_ERRORTEXT_psFits_METADATA_ADD_FAILED,
553 keyName);
554 return out;
555 }
556 }
557
558 }
559
560 if ( status != 0) {
561 char fitsErr[MAX_STRING_LENGTH];
562 (void)fits_get_errstatus(status, fitsErr);
563 psError(PS_ERR_IO, true,
564 PS_ERRORTEXT_psFits_METADATA_ADD_FAILED,
565 fitsErr);
566 return false;
567 }
568
569 return out;
570}
571
572psHash* psFitsReadHeaderSet(psHash* out,
573 const psFits* fits)
574{
575 if (fits == NULL) {
576 psError(PS_ERR_BAD_PARAMETER_NULL, true,
577 PS_ERRORTEXT_psFits_NULL);
578 psFree(out);
579 return NULL;
580 }
581
582 if (out == NULL) {
583 out = psHashAlloc(10);
584 // XXX: what is the appropriate number of buckets? Got 10 from psMetadataAlloc.
585 if (out == NULL) {
586 psError(PS_ERR_UNKNOWN, false,
587 "Failed to allocate a new psHash container.");
588 return NULL;
589 }
590 }
591
592 int size = psFitsGetSize(fits);
593
594 int origPosition = psFitsGetExtNum(fits);
595
596 for (int lcv=0; lcv < size; lcv++) {
597 psFitsMoveExtNum(fits, lcv, false);
598
599 char* name = NULL;
600 if (lcv == 0) {
601 name = psStringCopy("PHU");
602 } else {
603 name = psFitsGetExtName(fits);
604 }
605
606 psMetadata* header = psFitsReadHeader(NULL, fits);
607 if (name != NULL && header != NULL) {
608 psHashAdd(out, name, header);
609 } else { // XXX: is this a warning or error?
610 psLogMsg(__func__, PS_LOG_WARN,
611 "Failed to read HDU#%d header data.",
612 lcv);
613 }
614
615 psFree(name);
616 psFree(header);
617 }
618
619 // reposition to the original position
620 psFitsMoveExtNum(fits, origPosition, false);
621
622 return out;
623}
624
625psImage* psFitsReadImage(psImage* output, // a psImage to recycle.
626 const psFits* fits, // the psFits object
627 psRegion region, // the region in the FITS image to read
628 int z) // the z-plane in the FITS image cube to read
629{
630 psS32 status = 0; /* CFITSIO file vars */
631 psS32 nAxis = 0;
632 psS32 anynull = 0;
633 psS32 bitPix = 0; /* Pixel type */
634 long nAxes[3];
635 long firstPixel[3]; /* lower-left corner of image subset */
636 long lastPixel[3]; /* upper-right corner of image subset */
637 long increment[3]; /* increment for image subset */
638 char fitsErr[80] = ""; /* CFITSIO error message string */
639 psS32 fitsDatatype = 0;
640 psS32 datatype = 0;
641
642 if (fits == NULL) {
643 psError(PS_ERR_BAD_PARAMETER_NULL, true,
644 PS_ERRORTEXT_psFits_NULL);
645 psFree(output);
646 return NULL;
647 }
648
649 // check to see if we even are positioned on an image HDU
650 int hdutype;
651 if ( fits_get_hdu_type(fits->p_fd,&hdutype, &status) != 0) {
652 char fitsErr[MAX_STRING_LENGTH];
653 (void)fits_get_errstatus(status, fitsErr);
654 psError(PS_ERR_IO, true,
655 PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED,
656 fitsErr);
657 return NULL;
658 }
659 if (hdutype != IMAGE_HDU) {
660 psError(PS_ERR_IO, true,
661 PS_ERRORTEXT_psFits_NOT_IMAGE_TYPE);
662 return NULL;
663 }
664
665 /* Get the data type 'bitPix' from the FITS image */
666 if (fits_get_img_equivtype(fits->p_fd, &bitPix, &status) != 0) {
667 fits_get_errstatus(status, fitsErr);
668 psError(PS_ERR_IO, true,
669 PS_ERRORTEXT_psFits_DATATYPE_UNKNOWN,
670 fitsErr);
671 psFree(output);
672 return NULL;
673 }
674
675 /* Get the dimensions 'nAxis' from the FITS image */
676 if (fits_get_img_dim(fits->p_fd, &nAxis, &status) != 0) {
677 (void)fits_get_errstatus(status, fitsErr);
678 psError(PS_ERR_IO, true,
679 PS_ERRORTEXT_psFits_IMAGE_DIM_UNKNOWN,
680 fitsErr);
681 psFree(output);
682 return NULL;
683 }
684
685 /* Validate the number of axis */
686 if ((nAxis < 2) || (nAxis > 3)) {
687 psError(PS_ERR_IO, true,
688 PS_ERRORTEXT_psFits_IMAGE_DIMENSION_UNSUPPORTED,
689 nAxis);
690 psFree(output);
691 return NULL;
692 }
693
694 /* Get the Image size from the FITS file */
695 if (fits_get_img_size(fits->p_fd, nAxis, nAxes, &status) != 0) {
696 (void)fits_get_errstatus(status, fitsErr);
697 psError(PS_ERR_IO, true,
698 PS_ERRORTEXT_psFits_IMAGE_SIZE_UNKNOWN,
699 fitsErr);
700 psFree(output);
701 return NULL;
702 }
703
704 firstPixel[0] = region.x0 + 1;
705 firstPixel[1] = region.y0 + 1;
706 firstPixel[2] = z + 1;
707
708 if (region.x1 > 0) {
709 lastPixel[0] = region.x1;
710 } else {
711 lastPixel[0] = nAxes[0] + region.x1; // n.b., region.x1 < 0
712 }
713 if (region.y1 > 0) {
714 lastPixel[1] = region.y1;
715 } else {
716 lastPixel[1] = nAxes[1] + region.y1; // n.b., region.y1 < 0
717 }
718 lastPixel[2] = z + 1;
719
720 increment[0] = 1;
721 increment[1] = 1;
722 increment[2] = 1;
723
724 switch (bitPix) {
725 case BYTE_IMG:
726 datatype = PS_TYPE_U8;
727 fitsDatatype = TBYTE;
728 break;
729 case SBYTE_IMG:
730 datatype = PS_TYPE_S8;
731 fitsDatatype = TSBYTE;
732 break;
733 case USHORT_IMG:
734 datatype = PS_TYPE_U16;
735 fitsDatatype = TUSHORT;
736 break;
737 case SHORT_IMG:
738 datatype = PS_TYPE_S16;
739 fitsDatatype = TSHORT;
740 break;
741 case ULONG_IMG:
742 datatype = PS_TYPE_U32;
743 fitsDatatype = TUINT;
744 break;
745 case LONG_IMG:
746 datatype = PS_TYPE_S32;
747 fitsDatatype = TINT;
748 break;
749 case LONGLONG_IMG:
750 datatype = PS_TYPE_S64;
751 fitsDatatype = TLONGLONG;
752 break;
753 case FLOAT_IMG:
754 datatype = PS_TYPE_F32;
755 fitsDatatype = TFLOAT;
756 break;
757 case DOUBLE_IMG:
758 datatype = PS_TYPE_F64;
759 fitsDatatype = TDOUBLE;
760 break;
761 default:
762 psError(PS_ERR_IO, true,
763 PS_ERRORTEXT_psFits_FITS_TYPE_UNSUPPORTED,
764 bitPix);
765 psFree(output);
766 return NULL;
767 }
768
769 output = psImageRecycle(output,
770 lastPixel[0]-firstPixel[0]+1,
771 lastPixel[1]-firstPixel[1]+1,
772 datatype);
773
774 if (output == NULL) {
775 psError(PS_ERR_UNKNOWN, false,
776 "Failed to allocate a properly sized image.");
777 return false;
778 }
779
780 // n.b., this assumes contiguous image buffer
781 if (fits_read_subset(fits->p_fd, fitsDatatype, firstPixel, lastPixel, increment,
782 NULL, output->data.V[0], &anynull, &status) != 0) {
783 psFree(output);
784 (void)fits_get_errstatus(status, fitsErr);
785 psError(PS_ERR_IO, true,
786 PS_ERRORTEXT_psFits_READ_FAILED,
787 fitsErr);
788 return NULL;
789 }
790
791 return output;
792
793}
794
795bool psFitsWriteImage(const psFits* fits,
796 const psMetadata* header,
797 const psImage* input,
798 int numZPlanes,
799 char* extname)
800{
801
802 if (fits == NULL) {
803 psError(PS_ERR_BAD_PARAMETER_NULL, true,
804 PS_ERRORTEXT_psFits_NULL);
805 return false;
806 }
807
808 if (input == NULL) {
809 psError(PS_ERR_BAD_PARAMETER_NULL, true,
810 PS_ERRORTEXT_psFits_IMAGE_NULL);
811 return false;
812 }
813 int numCols = input->numCols;
814 int numRows = input->numRows;
815
816 int status = 0;
817
818 // determine the FITS-equivalent parameters
819 int bitPix;
820 double bZero;
821 int dataType;
822 if (! convertPsTypeToFits(input->type.type, &bitPix, &bZero, &dataType) ) {
823 return false;
824 }
825
826 int naxis = 3;
827 long naxes[3];
828
829 naxes[0] = numCols;
830 naxes[1] = numRows;
831 naxes[2] = numZPlanes;
832
833 if (numZPlanes < 2) {
834 naxis = 2;
835 }
836
837 fits_create_img(fits->p_fd, bitPix, naxis, naxes, &status);
838
839 if (extname != NULL) {
840 fits_update_key_str(fits->p_fd, "EXTNAME", (char*)extname, NULL, &status);
841 }
842
843 if (bZero != 0) { // set the bscale/bzero
844 fits_write_key_dbl(fits->p_fd, "BZERO", bZero, 12, "Pixel Value Offset", &status);
845 fits_write_key_dbl(fits->p_fd, "BSCALE", 1.0, 12, "Pixel Value Scale", &status);
846 fits_set_bscale(fits->p_fd, 1.0, bZero, &status);
847 }
848
849 // write the header, if any.
850 if (header != NULL) {
851 psFitsWriteHeader(header, fits);
852 }
853
854 if (input->parent == NULL) { // if no parent, assume that the image data is contiguous
855 fits_write_img(fits->p_fd,
856 dataType, // datatype
857 1, // writing to the first z-plane
858 numCols*numRows, // number of elements to write, i.e., the whole image
859 input->data.V[0], // the data
860 &status);
861 } else { // image data may not be contiguous; write one row at a time
862 int firstPixel = 1;
863 for (int row = 0; row < numRows; row++) {
864 fits_write_img(fits->p_fd,
865 dataType, // datatype
866 firstPixel,
867 numCols, // number of elements to write, i.e., one row's worth
868 input->data.V[row],// the raw row data
869 &status);
870 firstPixel += numCols; // move to next row
871 }
872 }
873
874 if ( status != 0) {
875 char fitsErr[MAX_STRING_LENGTH];
876 (void)fits_get_errstatus(status, fitsErr);
877 psError(PS_ERR_IO, true,
878 PS_ERRORTEXT_psFits_WRITE_FAILED,
879 fits->filename, fitsErr);
880 return false;
881 }
882
883 return true;
884
885}
886
887bool psFitsUpdateImage(const psFits* fits,
888 const psImage* input,
889 psRegion region,
890 int z)
891{
892 int status = 0;
893
894 if (fits == NULL) {
895 psError(PS_ERR_BAD_PARAMETER_NULL, true,
896 PS_ERRORTEXT_psFits_NULL);
897 return false;
898 }
899
900 if (input == NULL) {
901 psError(PS_ERR_BAD_PARAMETER_NULL, true,
902 PS_ERRORTEXT_psFits_IMAGE_NULL);
903 return false;
904 }
905
906 // check to see if we are positioned on an image HDU
907 int hdutype;
908 if ( fits_get_hdu_type(fits->p_fd,&hdutype, &status) != 0) {
909 char fitsErr[MAX_STRING_LENGTH];
910 (void)fits_get_errstatus(status, fitsErr);
911 psError(PS_ERR_IO, true,
912 PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED,
913 fitsErr);
914 return NULL;
915 }
916 if (hdutype != IMAGE_HDU) {
917 psError(PS_ERR_IO, true,
918 PS_ERRORTEXT_psFits_NOT_IMAGE_TYPE);
919 return NULL;
920 }
921
922 int numCols = input->numCols;
923 int numRows = input->numRows;
924
925 // determine the FITS-equivalent parameters
926 int bitPix;
927 double bZero;
928 int dataType;
929 if (! convertPsTypeToFits(input->type.type, &bitPix, &bZero, &dataType) ) {
930 return false;
931 }
932
933 //check to see if the HDU has the same datatype
934 int fileBitpix;
935 int naxis;
936 long nAxes[3];
937 nAxes[2] = 1;
938 fits_get_img_param(fits->p_fd, 3, &fileBitpix, &naxis, nAxes, &status);
939
940 //check to see if the HDU has the same datatype
941 if (bitPix != fileBitpix) {
942 char* fitsTypeStr;
943 char* imageTypeStr;
944 PS_TYPE_NAME(fitsTypeStr,fileBitpix);
945 PS_TYPE_NAME(imageTypeStr,input->type.type);
946 psError(PS_ERR_IO, true,
947 PS_ERRORTEXT_psFits_IMAGE_UPDATE_TYPE_MISMATCH,
948 fitsTypeStr, imageTypeStr);
949 return false;
950 }
951
952 //check if the HDU has the z-plane requested
953 if (z >= nAxes[2]) {
954 psError(PS_ERR_BAD_PARAMETER_SIZE, true,
955 PS_ERRORTEXT_psFits_FITS_Z_SMALL,
956 nAxes[2],z);
957 return false;
958 }
959
960 // determine the region in the FITS file domain
961 long firstPixel[3];
962 long lastPixel[3];
963
964 firstPixel[0] = region.x0 + 1;
965 firstPixel[1] = region.y0 + 1;
966 firstPixel[2] = z + 1;
967
968 if (region.x1 > 0) {
969 lastPixel[0] = region.x1;
970 } else {
971 lastPixel[0] = nAxes[0] + region.x1; // n.b., region.x1 < 0
972 }
973 if (region.y1 > 0) {
974 lastPixel[1] = region.y1;
975 } else {
976 lastPixel[1] = nAxes[1] + region.y1; // n.b., region.y1 < 0
977 }
978 lastPixel[2] = z + 1;
979
980 if (firstPixel[0] < 1 || firstPixel[0] > nAxes[0] ||
981 firstPixel[1] < 1 || firstPixel[1] > nAxes[1] ||
982 lastPixel[0] < 1 || lastPixel[0] > nAxes[0] ||
983 lastPixel[1] < 1 || lastPixel[1] > nAxes[1]) {
984 psError(PS_ERR_BAD_PARAMETER_VALUE, true,
985 "Specified region [%d:%d,%d:%d], is not valid given the %dx%d FITS image.",
986 region.y0,region.y1-1,region.x0,region.x1-1);
987 return false;
988
989 }
990
991 int dx = lastPixel[0] - firstPixel[0];
992 int dy = lastPixel[1] - firstPixel[1];
993 if (dx > numCols ||
994 dy > numRows) {
995 psError(PS_ERR_BAD_PARAMETER_VALUE, true,
996 "The region [%d:%d,%d:%d], is not valid given the input %dx%d image.",
997 firstPixel[1]-1,lastPixel[1]-1,
998 firstPixel[0]-1,lastPixel[0]-1,
999 numCols, numRows);
1000 return false;
1001 }
1002
1003 psImage* subset;
1004 if (dx != numCols || dy != numRows) {
1005 // the input image needs to be subsetted
1006 subset = psImageSubset((psImage*)input,0,0,dx+1,dy+1);
1007 } else {
1008 subset = psMemIncrRefCounter((psImage*)input);
1009 }
1010
1011 fits_write_subset(fits->p_fd, dataType, firstPixel, lastPixel, subset->data.V[0], &status);
1012
1013 if ( status != 0) {
1014 char fitsErr[MAX_STRING_LENGTH];
1015 (void)fits_get_errstatus(status, fitsErr);
1016 psError(PS_ERR_IO, true,
1017 PS_ERRORTEXT_psFits_WRITE_FAILED,
1018 fits->filename, fitsErr);
1019 return false;
1020 }
1021
1022 return true;
1023}
1024
1025bool psFitsWriteHeader(const psMetadata* header,
1026 const psFits* fits)
1027{
1028
1029 if (fits == NULL) {
1030 psError(PS_ERR_BAD_PARAMETER_NULL, true,
1031 PS_ERRORTEXT_psFits_NULL);
1032 return false;
1033 }
1034
1035 if (header == NULL) {
1036 psError(PS_ERR_BAD_PARAMETER_NULL, true,
1037 PS_ERRORTEXT_psFits_METADATA_NULL);
1038 return false;
1039 }
1040
1041 int status = 0;
1042
1043 //transverse the metadata list and add each key.
1044
1045 psListIterator* iter = psListIteratorAlloc(header->list,PS_LIST_HEAD,true);
1046 psMetadataItem* item;
1047 while ( (item=psListGetAndIncrement(iter)) != NULL ) {
1048 switch (item->type) {
1049 case PS_META_BOOL: {
1050 int value = item->data.B;
1051 fits_update_key(fits->p_fd,
1052 TLOGICAL,
1053 item->name,
1054 &value,
1055 item->comment,
1056 &status);
1057 break;
1058 }
1059 case PS_META_S32:
1060 fits_update_key(fits->p_fd,
1061 TINT,
1062 item->name,
1063 &item->data.S32,
1064 item->comment,
1065 &status);
1066 break;
1067 case PS_META_F32:
1068 fits_update_key(fits->p_fd,
1069 TFLOAT,
1070 item->name,
1071 &item->data.F32,
1072 item->comment,
1073 &status);
1074 break;
1075 case PS_META_F64:
1076 fits_update_key(fits->p_fd,
1077 TDOUBLE,
1078 item->name,
1079 &item->data.F64,
1080 item->comment,
1081 &status);
1082 break;
1083 case PS_META_STR:
1084 fits_update_key(fits->p_fd,
1085 TSTRING,
1086 item->name,
1087 item->data.V,
1088 item->comment,
1089 &status);
1090 break;
1091 default: // all other META types are ignored
1092 break;
1093 }
1094
1095 if ( status != 0) {
1096 char fitsErr[MAX_STRING_LENGTH];
1097 (void)fits_get_errstatus(status, fitsErr);
1098 psError(PS_ERR_IO, true,
1099 PS_ERRORTEXT_psFits_WRITE_FAILED,
1100 fits->filename, fitsErr);
1101 return false;
1102 }
1103 }
1104
1105 return true;
1106}
1107
1108psMetadata* psFitsReadTableRow(const psFits* fits,
1109 int row)
1110{
1111 long numRows;
1112 int numCols;
1113 int status = 0;
1114
1115 if (fits == NULL) {
1116 psError(PS_ERR_BAD_PARAMETER_NULL, true,
1117 PS_ERRORTEXT_psFits_NULL);
1118 return NULL;
1119 }
1120
1121 // check to see if we even are positioned on a table HDU
1122 int hdutype;
1123 fits_get_hdu_type(fits->p_fd,&hdutype, &status);
1124 if ( status != 0) {
1125 char fitsErr[MAX_STRING_LENGTH];
1126 (void)fits_get_errstatus(status, fitsErr);
1127 psError(PS_ERR_IO, true,
1128 PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED,
1129 fitsErr);
1130 return NULL;
1131 }
1132 if (hdutype != ASCII_TBL && hdutype != BINARY_TBL) {
1133 psError(PS_ERR_IO, true,
1134 PS_ERRORTEXT_psFits_NOT_TABLE_TYPE);
1135 return NULL;
1136 }
1137
1138 // get the size of the FITS table
1139 fits_get_num_rows(fits->p_fd, &numRows, &status);
1140 fits_get_num_cols(fits->p_fd, &numCols, &status);
1141 if ( status != 0) {
1142 char fitsErr[MAX_STRING_LENGTH];
1143 (void)fits_get_errstatus(status, fitsErr);
1144 psError(PS_ERR_IO, true,
1145 PS_ERRORTEXT_psFits_GET_TABLE_SIZE_FAILED,
1146 fitsErr);
1147 return NULL;
1148 }
1149
1150 // the row parameter in the proper range?
1151 if (row < 0 || row >= numRows) {
1152 psError(PS_ERR_BAD_PARAMETER_VALUE, true,
1153 PS_ERRORTEXT_psFits_ROW_INVALID,
1154 row,numRows);
1155 return NULL;
1156 }
1157
1158 psMetadata* data = psMetadataAlloc();
1159
1160 int typecode;
1161 long repeat;
1162 long width;
1163 char name[60];
1164 for (int col = 1; col <= numCols; col++) {
1165 // get the column name
1166 if (hdutype == BINARY_TBL) {
1167 fits_get_bcolparms(fits->p_fd, col, name,
1168 NULL, NULL, NULL, NULL, NULL, NULL, NULL, &status);
1169 } else {
1170 fits_get_acolparms(fits->p_fd, col, name,
1171 NULL, NULL, NULL, NULL, NULL, NULL, NULL, &status);
1172 }
1173
1174 // get the column type
1175 fits_get_coltype(fits->p_fd, col, &typecode, &repeat, &width, &status);
1176
1177 if (status == 0) {
1178
1179 #define READ_TABLE_ROW_CASE(FITSTYPE, NATIVETYPE, TYPE) \
1180 case FITSTYPE: { \
1181 NATIVETYPE value = 0; \
1182 int anynul = 0; \
1183 fits_read_col(fits->p_fd, FITSTYPE, col,row+1, \
1184 1, 1, NULL, &value, &anynul, &status); \
1185 psMetadataAdd(data,PS_LIST_TAIL, name, \
1186 PS_META_##TYPE, \
1187 "", (ps##TYPE)value); \
1188 break; \
1189 }
1190
1191 switch (typecode) {
1192 case TBYTE:
1193 case TSHORT:
1194 case TLONGLONG:
1195 READ_TABLE_ROW_CASE(TLONG, long, S32)
1196 READ_TABLE_ROW_CASE(TFLOAT, float, F32)
1197 READ_TABLE_ROW_CASE(TDOUBLE, double, F64)
1198 READ_TABLE_ROW_CASE(TLOGICAL, bool, BOOL);
1199 case TSTRING: {
1200 char* value;
1201 int anynul = 0;
1202 fits_read_col(fits->p_fd, TSTRING, col,row+1,
1203 1, 1, NULL, &value, &anynul, &status);
1204 if (anynul == 0) {
1205 psMetadataAdd(data,PS_LIST_TAIL, name,
1206 PS_META_STR,
1207 "", value);
1208 }
1209 break;
1210 }
1211 default:
1212 psTrace("psFits.psFitsReadTableRow", 2,
1213 "Column %d or row %d was of a non primitive type, %d",
1214 col, row, typecode);
1215 }
1216 }
1217
1218 if ( status != 0) {
1219 char fitsErr[MAX_STRING_LENGTH];
1220 (void)fits_get_errstatus(status, fitsErr);
1221 psError(PS_ERR_IO, true,
1222 PS_ERRORTEXT_psFits_GET_TABLE_ELEMENT,
1223 col,row,fitsErr);
1224 psFree(data);
1225 return NULL;
1226 }
1227
1228 }
1229
1230 return data;
1231}
1232
1233psArray* psFitsReadTableColumn(const psFits* fits,
1234 const char* colname)
1235{
1236 int colnum = 0;
1237 int status = 0;
1238
1239 if (fits == NULL) {
1240 psError(PS_ERR_BAD_PARAMETER_NULL, true,
1241 PS_ERRORTEXT_psFits_NULL);
1242 return NULL;
1243 }
1244
1245 // check to see if we even are positioned on a table HDU
1246 int hdutype;
1247 if ( fits_get_hdu_type(fits->p_fd,&hdutype, &status) != 0) {
1248 char fitsErr[MAX_STRING_LENGTH];
1249 (void)fits_get_errstatus(status, fitsErr);
1250 psError(PS_ERR_IO, true,
1251 PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED,
1252 fitsErr);
1253 return NULL;
1254 }
1255 if (hdutype != ASCII_TBL && hdutype != BINARY_TBL) {
1256 psError(PS_ERR_IO, true,
1257 PS_ERRORTEXT_psFits_NOT_TABLE_TYPE);
1258 return NULL;
1259 }
1260
1261 // find the column by name
1262 if ( fits_get_colnum(fits->p_fd, CASESEN, (char*)colname, &colnum, &status) != 0) {
1263 char fitsErr[MAX_STRING_LENGTH];
1264 (void)fits_get_errstatus(status, fitsErr);
1265 psError(PS_ERR_IO, true,
1266 PS_ERRORTEXT_psFits_FIND_COLUMN,
1267 colname, fitsErr);
1268 return NULL;
1269 }
1270
1271 // get the number of rows
1272 long numRows = 0;
1273 fits_get_num_rows(fits->p_fd, &numRows, &status);
1274
1275 // get the column length.
1276 int width;
1277 if ( fits_get_col_display_width(fits->p_fd, colnum, &width, &status) != 0) {
1278 char fitsErr[MAX_STRING_LENGTH];
1279 (void)fits_get_errstatus(status, fitsErr);
1280 psError(PS_ERR_IO, true,
1281 PS_ERRORTEXT_psFits_GET_COLTYPE,
1282 fitsErr);
1283 return NULL;
1284 }
1285
1286 // allocate the buffers
1287 psArray* result = psArrayAlloc(numRows);
1288 for (int row = 0; row < numRows; row++) {
1289 result->data[row] = psAlloc((width+1)*sizeof(char));
1290 }
1291 result->n = numRows;
1292
1293 fits_read_col_str(fits->p_fd,
1294 colnum,
1295 1, // firstrow
1296 1, // firestelem
1297 numRows,
1298 "", // nulstr
1299 (char**)result->data,
1300 NULL,
1301 &status);
1302
1303 if ( status != 0) {
1304 char fitsErr[MAX_STRING_LENGTH];
1305 (void)fits_get_errstatus(status, fitsErr);
1306 psError(PS_ERR_IO, true,
1307 PS_ERRORTEXT_psFits_TABLE_READ_COL,
1308 fitsErr);
1309 return NULL;
1310 }
1311
1312 return result;
1313}
1314
1315psVector* psFitsReadTableColumnNum(const psFits* fits,
1316 const char* colname)
1317{
1318 int status = 0;
1319 int colnum = 0;
1320
1321 if (fits == NULL) {
1322 psError(PS_ERR_BAD_PARAMETER_NULL, true,
1323 PS_ERRORTEXT_psFits_NULL);
1324 return NULL;
1325 }
1326
1327 // check to see if we even are positioned on a table HDU
1328 int hdutype;
1329 if ( fits_get_hdu_type(fits->p_fd,&hdutype, &status) != 0) {
1330 char fitsErr[MAX_STRING_LENGTH];
1331 (void)fits_get_errstatus(status, fitsErr);
1332 psError(PS_ERR_IO, true,
1333 PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED,
1334 fitsErr);
1335 return NULL;
1336 }
1337 if (hdutype != ASCII_TBL && hdutype != BINARY_TBL) {
1338 psError(PS_ERR_IO, true,
1339 PS_ERRORTEXT_psFits_NOT_TABLE_TYPE);
1340 return NULL;
1341 }
1342
1343 // find the column by name
1344 if ( fits_get_colnum(fits->p_fd, CASESEN, (char*)colname, &colnum, &status) != 0) {
1345 char fitsErr[MAX_STRING_LENGTH];
1346 (void)fits_get_errstatus(status, fitsErr);
1347 psError(PS_ERR_IO, true,
1348 PS_ERRORTEXT_psFits_FIND_COLUMN,
1349 colname, fitsErr);
1350 return NULL;
1351 }
1352
1353 // get the number of rows
1354 long numRows = 0;
1355 fits_get_num_rows(fits->p_fd,
1356 &numRows,
1357 &status);
1358
1359 // get the column datatype.
1360 int typecode;
1361 long repeat;
1362 long width;
1363 if ( fits_get_eqcoltype(fits->p_fd, colnum, &typecode, &repeat, &width, &status) != 0) {
1364 char fitsErr[MAX_STRING_LENGTH];
1365 (void)fits_get_errstatus(status, fitsErr);
1366 psError(PS_ERR_IO, true,
1367 PS_ERRORTEXT_psFits_GET_COLTYPE,
1368 fitsErr);
1369 return NULL;
1370 }
1371
1372 psVector* result = psVectorAlloc(numRows, convertFitsToPsType(typecode));
1373
1374 fits_read_col(fits->p_fd,
1375 typecode,
1376 colnum,
1377 1 /* firstrow */,
1378 1 /* firstelem */,
1379 numRows,
1380 NULL,
1381 (psPtr)(result->data.U8),
1382 NULL,
1383 &status);
1384
1385 if ( status != 0) {
1386 char fitsErr[MAX_STRING_LENGTH];
1387 (void)fits_get_errstatus(status, fitsErr);
1388 psError(PS_ERR_IO, true,
1389 PS_ERRORTEXT_psFits_TABLE_READ_COL,
1390 fitsErr);
1391 return NULL;
1392 }
1393
1394 return result;
1395}
1396
1397
1398psArray* psFitsReadTable(const psFits* fits)
1399{
1400 int status = 0;
1401
1402 if (fits == NULL) {
1403 psError(PS_ERR_BAD_PARAMETER_NULL, true,
1404 PS_ERRORTEXT_psFits_NULL);
1405 return NULL;
1406 }
1407
1408 // check to see if we even are positioned on a table HDU
1409 int hdutype;
1410 if ( fits_get_hdu_type(fits->p_fd,&hdutype, &status) != 0) {
1411 char fitsErr[MAX_STRING_LENGTH];
1412 (void)fits_get_errstatus(status, fitsErr);
1413 psError(PS_ERR_IO, true,
1414 PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED,
1415 fitsErr);
1416 return NULL;
1417 }
1418 if (hdutype != ASCII_TBL && hdutype != BINARY_TBL) {
1419 psError(PS_ERR_IO, true,
1420 PS_ERRORTEXT_psFits_NOT_TABLE_TYPE);
1421 return NULL;
1422 }
1423
1424 // get the size of the FITS table
1425 long numRows = 0;
1426 fits_get_num_rows(fits->p_fd, &numRows, &status);
1427 if ( status != 0) {
1428 char fitsErr[MAX_STRING_LENGTH];
1429 (void)fits_get_errstatus(status, fitsErr);
1430 psError(PS_ERR_IO, true,
1431 PS_ERRORTEXT_psFits_GET_TABLE_SIZE_FAILED,
1432 fitsErr);
1433 return NULL;
1434 }
1435
1436 psArray* table = psArrayAlloc(numRows);
1437
1438 for (int row = 0; row < numRows; row++) {
1439 table->data[row] = psFitsReadTableRow(fits,row);
1440 }
1441
1442 return table;
1443}
1444
1445bool psFitsWriteTable(const psFits* fits,
1446 psMetadata* header,
1447 psArray* table,
1448 char* extname)
1449{
1450 int status = 0;
1451 psMetadataItem* item;
1452
1453 if (fits == NULL) {
1454 psError(PS_ERR_BAD_PARAMETER_NULL, true,
1455 PS_ERRORTEXT_psFits_NULL);
1456 return false;
1457 }
1458
1459 if (table == NULL) {
1460 psError(PS_ERR_BAD_PARAMETER_NULL, true,
1461 PS_ERRORTEXT_psFits_IMAGE_NULL);
1462 return false;
1463 }
1464
1465 int rows = table->n;
1466 if (rows < 1) {
1467 // no table data, what can I do?
1468 psError(PS_ERR_BAD_PARAMETER_SIZE, true,
1469 PS_ERRORTEXT_psFits_TABLE_EMPTY);
1470 return false;
1471 }
1472
1473 // find all the columns needed
1474 psArray* columns = psArrayAlloc(((psMetadata*)table->data[0])->list->size);
1475 columns->n=0;
1476
1477 // find the unique items in the array of metadata 'rows'
1478 for (int row=0; row < rows; row++) {
1479 psMetadata* rowMeta = table->data[row];
1480 if (rowMeta != NULL) {
1481 psListIterator* iter = psListIteratorAlloc(rowMeta->list,
1482 PS_LIST_HEAD,true);
1483 while ( (item=psListGetAndIncrement(iter)) != NULL) {
1484 if (PS_META_IS_PRIMITIVE(item->type)) {
1485 bool found = false;
1486 for (int n=0; n < columns->n && ! found; n++) {
1487 if (strcmp(item->name,
1488 ((psMetadataItem*)(columns->data[n]))->name) == 0) {
1489 found = true;
1490 }
1491 }
1492 if (! found) {
1493 psArrayAdd(columns, columns->nalloc, item);
1494 }
1495 }
1496 }
1497 psFree(iter);
1498 }
1499 }
1500
1501 if (columns->n == 0) { // no table columns found
1502 // XXX: Error?
1503 return false;
1504 }
1505
1506 //create list of column names and types.
1507 psArray* columnNames = psArrayAlloc(columns->n);
1508 psArray* columnTypes = psArrayAlloc(columns->n);
1509 for (int n=0; n < columns->n; n++) {
1510 char* fitsType;
1511 columnNames->data[n] = psMemIncrRefCounter(((psMetadataItem*)columns->data[n])->name);
1512 if ( ! convertMetadataTypeToBinaryTForm(((psMetadataItem*)columns->data[n])->type,
1513 &fitsType)) {
1514 // XXX: error message
1515 return false;
1516 }
1517 columnTypes->data[n] = fitsType;
1518 }
1519
1520 fits_create_tbl(fits->p_fd,
1521 BINARY_TBL,
1522 table->n, // number of rows in table
1523 columns->n, // number of columns in table
1524 (char**)columnNames->data, // names of the columns
1525 (char**)columnTypes->data, // format of the columns
1526 NULL, // physical unit of columns
1527 extname, // extension name
1528 &status);
1529
1530 psFree(columnNames);
1531 psFree(columnTypes);
1532
1533 // fill in the table elements with data
1534 for (int n = 0; n < columns->n; n++) {
1535 int row;
1536 item = columns->data[n];
1537 if (PS_META_IS_PRIMITIVE(item->type)) {
1538 psVector* col = NULL;
1539 switch (item->type) {
1540 case PS_META_S32:
1541 col = psVectorAlloc(table->n, PS_TYPE_S32);
1542 for (row = 0; row < table->n; row++) {
1543 col->data.S32[row] = psMetadataLookupS32(NULL,
1544 table->data[row],
1545 item->name);
1546 }
1547 fits_write_col_int(fits->p_fd,
1548 n+1, // column number
1549 1, // firstrow
1550 1, // firstelem
1551 table->n, // nelements
1552 col->data.S32,
1553 &status);
1554 break;
1555 case PS_META_F32:
1556 col = psVectorAlloc(table->n, PS_TYPE_F32);
1557 for (row = 0; row < table->n; row++) {
1558 col->data.F32[row] = psMetadataLookupF32(NULL,
1559 table->data[row],
1560 item->name);
1561 }
1562 fits_write_col_flt(fits->p_fd,
1563 n+1, // column number
1564 1, // firstrow
1565 1, // firstelem
1566 table->n, // nelements
1567 col->data.F32,
1568 &status);
1569 break;
1570 case PS_META_F64:
1571 col = psVectorAlloc(table->n, PS_TYPE_F64);
1572 for (row = 0; row < table->n; row++) {
1573 col->data.F64[row] = psMetadataLookupF64(NULL,
1574 table->data[row],
1575 item->name);
1576 }
1577 fits_write_col_dbl(fits->p_fd,
1578 n+1, // column number
1579 1, // firstrow
1580 1, // firstelem
1581 table->n, // nelements
1582 col->data.F64,
1583 &status);
1584 break;
1585 case PS_META_BOOL:
1586 col = psVectorAlloc(table->n, PS_TYPE_BOOL);
1587 for (row = 0; row < table->n; row++) {
1588 col->data.S8[row] = psMetadataLookupBool(NULL,
1589 table->data[row],
1590 item->name);
1591 }
1592 fits_write_col_log(fits->p_fd,
1593 n+1, // column number
1594 1, // firstrow
1595 1, // firstelem
1596 table->n, // nelements
1597 col->data.S8,
1598 &status);
1599 break;
1600 default:
1601 // XXX: error message?
1602 break;
1603 }
1604 psFree(col);
1605 } else if (item->type == PS_META_STR) {
1606 psArray* col = psArrayAlloc(table->n);
1607 for (row = 0; row < table->n; row++) {
1608 col->data[row] = item->data.V;
1609 }
1610 fits_write_col_str(fits->p_fd,
1611 n, // column number
1612 1, // firstrow
1613 1, // firstelem
1614 table->n, // nelements
1615 (char**)col->data,
1616 &status);
1617 psFree(col);
1618 }
1619 }
1620
1621 psFree(columns);
1622
1623 return true;
1624}
1625
1626bool psFitsUpdateTable(const psFits* fits,
1627 psMetadata* data,
1628 int row)
1629{
1630 int status = 0;
1631
1632 if (fits == NULL) {
1633 psError(PS_ERR_BAD_PARAMETER_NULL, true,
1634 PS_ERRORTEXT_psFits_NULL);
1635 return false;
1636 }
1637
1638 if (data == NULL) {
1639 psError(PS_ERR_BAD_PARAMETER_NULL, true,
1640 PS_ERRORTEXT_psFits_IMAGE_NULL);
1641 return false;
1642 }
1643
1644 // check to see if we even are positioned on a table HDU
1645 int hdutype;
1646 if ( fits_get_hdu_type(fits->p_fd,&hdutype, &status) != 0) {
1647 char fitsErr[MAX_STRING_LENGTH];
1648 (void)fits_get_errstatus(status, fitsErr);
1649 psError(PS_ERR_IO, true,
1650 PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED,
1651 fitsErr);
1652 return false;
1653 }
1654 if (hdutype != ASCII_TBL && hdutype != BINARY_TBL) {
1655 psError(PS_ERR_IO, true,
1656 PS_ERRORTEXT_psFits_NOT_TABLE_TYPE);
1657 return false;
1658 }
1659
1660 psMetadataIterator* iter = psMetadataIteratorAlloc(data,PS_LIST_HEAD,NULL);
1661
1662 psMetadataItem* item;
1663
1664 while ( (item=psMetadataGetAndIncrement(iter)) != NULL) {
1665 if (PS_META_IS_PRIMITIVE(item->type)) {
1666 // operating on primitive data type, i.e., not a complex object
1667 int colnum = 0;
1668
1669 if ( fits_get_colnum(fits->p_fd, CASESEN, item->name, &colnum, &status) == 0) {
1670 // cooresponding column found in table
1671 int dataType = 0;
1672 convertPsTypeToFits(item->type, NULL, NULL, &dataType);
1673
1674 if (fits_write_col(fits->p_fd, dataType, colnum, row+1, 1, 1, &item->data,&status) != 0) {
1675 char fitsErr[MAX_STRING_LENGTH];
1676 (void)fits_get_errstatus(status, fitsErr);
1677 psError(PS_ERR_IO, true,
1678 PS_ERRORTEXT_psFits_WRITE_FAILED,
1679 fits->filename, fitsErr);
1680 psFree(iter);
1681 return false;
1682 }
1683 } else {
1684 // the column was not found.
1685 psWarning("No column with the name '%s' exists in the table.",
1686 item->name);
1687 }
1688 }
1689 }
1690
1691 psFree(iter);
1692
1693 return true;
1694}