1 /************************************************************************/
2 /* */
3 /* CFITSIO Lexical Parser */
4 /* */
5 /* This file is one of 3 files containing code which parses an */
6 /* arithmetic expression and evaluates it in the context of an input */
7 /* FITS file table extension. The CFITSIO lexical parser is divided */
8 /* into the following 3 parts/files: the CFITSIO "front-end", */
9 /* eval_f.c, contains the interface between the user/CFITSIO and the */
10 /* real core of the parser; the FLEX interpreter, eval_l.c, takes the */
11 /* input string and parses it into tokens and identifies the FITS */
12 /* information required to evaluate the expression (ie, keywords and */
13 /* columns); and, the BISON grammar and evaluation routines, eval_y.c, */
14 /* receives the FLEX output and determines and performs the actual */
15 /* operations. The files eval_l.c and eval_y.c are produced from */
16 /* running flex and bison on the files eval.l and eval.y, respectively. */
17 /* (flex and bison are available from any GNU archive: see www.gnu.org) */
18 /* */
19 /* The grammar rules, rather than evaluating the expression in situ, */
20 /* builds a tree, or Nodal, structure mapping out the order of */
21 /* operations and expression dependencies. This "compilation" process */
22 /* allows for much faster processing of multiple rows. This technique */
23 /* was developed by Uwe Lammers of the XMM Science Analysis System, */
24 /* although the CFITSIO implementation is entirely code original. */
25 /* */
26 /* */
27 /* Modification History: */
28 /* */
29 /* Kent Blackburn c1992 Original parser code developed for the */
30 /* FTOOLS software package, in particular, */
31 /* the fselect task. */
32 /* Kent Blackburn c1995 BIT column support added */
33 /* Peter D Wilson Feb 1998 Vector column support added */
34 /* Peter D Wilson May 1998 Ported to CFITSIO library. User */
35 /* interface routines written, in essence */
36 /* making fselect, fcalc, and maketime */
37 /* capabilities available to all tools */
38 /* via single function calls. */
39 /* Peter D Wilson Jun 1998 Major rewrite of parser core, so as to */
40 /* create a run-time evaluation tree, */
41 /* inspired by the work of Uwe Lammers, */
42 /* resulting in a speed increase of */
43 /* 10-100 times. */
44 /* Peter D Wilson Jul 1998 gtifilter(a,b,c,d) function added */
45 /* Peter D Wilson Aug 1998 regfilter(a,b,c,d) function added */
46 /* Peter D Wilson Jul 1999 Make parser fitsfile-independent, */
47 /* allowing a purely vector-based usage */
48 /* Peter D Wilson Aug 1999 Add row-offset capability */
49 /* Peter D Wilson Sep 1999 Add row-range capability to ffcalc_rng */
50 /* */
51 /************************************************************************/
52
53 #include <limits.h>
54 #include <ctype.h>
55 #include "eval_defs.h"
56 #include "region.h"
57
58 typedef struct {
59 int datatype; /* Data type to cast parse results into for user */
60 void *dataPtr; /* Pointer to array of results, NULL if to use iterCol */
61 void *nullPtr; /* Pointer to nulval, use zero if NULL */
62 long maxRows; /* Max No. of rows to process, -1=all, 0=1 iteration */
63 int anyNull; /* Flag indicating at least 1 undef value encountered */
64 } parseInfo;
65
66 /* Internal routines needed to allow the evaluator to operate on FITS data */
67
68 static void Setup_DataArrays( int nCols, iteratorCol *cols,
69 long fRow, long nRows );
70 static int find_column( char *colName, void *itslval );
71 static int find_keywd ( char *key, void *itslval );
72 static int allocateCol( int nCol, int *status );
73 static int load_column( int varNum, long fRow, long nRows,
74 void *data, char *undef );
75
76 static int DEBUG_PIXFILTER;
77
78 #define FREE(x) { if (x) free(x); else printf("invalid free(" #x ") at %s:%d\n", __FILE__, __LINE__); }
79
80 /*---------------------------------------------------------------------------*/
fffrow(fitsfile * fptr,char * expr,long firstrow,long nrows,long * n_good_rows,char * row_status,int * status)81 int fffrow( fitsfile *fptr, /* I - Input FITS file */
82 char *expr, /* I - Boolean expression */
83 long firstrow, /* I - First row of table to eval */
84 long nrows, /* I - Number of rows to evaluate */
85 long *n_good_rows, /* O - Number of rows eval to True */
86 char *row_status, /* O - Array of boolean results */
87 int *status ) /* O - Error status */
88 /* */
89 /* Evaluate a boolean expression using the indicated rows, returning an */
90 /* array of flags indicating which rows evaluated to TRUE/FALSE */
91 /*---------------------------------------------------------------------------*/
92 {
93 parseInfo Info;
94 int naxis, constant;
95 long nelem, naxes[MAXDIMS], elem;
96 char result;
97
98 if( *status ) return( *status );
99
100 FFLOCK;
101 if( ffiprs( fptr, 0, expr, MAXDIMS, &Info.datatype, &nelem, &naxis,
102 naxes, status ) ) {
103 ffcprs();
104 FFUNLOCK;
105 return( *status );
106 }
107 if( nelem<0 ) {
108 constant = 1;
109 nelem = -nelem;
110 } else
111 constant = 0;
112
113 if( Info.datatype!=TLOGICAL || nelem!=1 ) {
114 ffcprs();
115 ffpmsg("Expression does not evaluate to a logical scalar.");
116 FFUNLOCK;
117 return( *status = PARSE_BAD_TYPE );
118 }
119
120 if( constant ) { /* No need to call parser... have result from ffiprs */
121 result = gParse.Nodes[gParse.resultNode].value.data.log;
122 *n_good_rows = nrows;
123 for( elem=0; elem<nrows; elem++ )
124 row_status[elem] = result;
125 } else {
126 firstrow = (firstrow>1 ? firstrow : 1);
127 Info.dataPtr = row_status;
128 Info.nullPtr = NULL;
129 Info.maxRows = nrows;
130
131 if( ffiter( gParse.nCols, gParse.colData, firstrow-1, 0,
132 parse_data, (void*)&Info, status ) == -1 )
133 *status = 0; /* -1 indicates exitted without error before end... OK */
134
135 if( *status ) {
136
137 /***********************/
138 /* Error... Do nothing */
139 /***********************/
140
141 } else {
142
143 /***********************************/
144 /* Count number of good rows found */
145 /***********************************/
146
147 *n_good_rows = 0L;
148 for( elem=0; elem<Info.maxRows; elem++ ) {
149 if( row_status[elem]==1 ) ++*n_good_rows;
150 }
151 }
152 }
153
154 ffcprs();
155 FFUNLOCK;
156 return(*status);
157 }
158
159 /*--------------------------------------------------------------------------*/
ffsrow(fitsfile * infptr,fitsfile * outfptr,char * expr,int * status)160 int ffsrow( fitsfile *infptr, /* I - Input FITS file */
161 fitsfile *outfptr, /* I - Output FITS file */
162 char *expr, /* I - Boolean expression */
163 int *status ) /* O - Error status */
164 /* */
165 /* Evaluate an expression on all rows of a table. If the input and output */
166 /* files are not the same, copy the TRUE rows to the output file. If the */
167 /* files are the same, delete the FALSE rows (preserve the TRUE rows). */
168 /* Can copy rows between extensions of the same file, *BUT* if output */
169 /* extension is before the input extension, the second extension *MUST* be */
170 /* opened using ffreopen, so that CFITSIO can handle changing file lengths. */
171 /*--------------------------------------------------------------------------*/
172 {
173 parseInfo Info;
174 int naxis, constant;
175 long nelem, rdlen, naxes[MAXDIMS], maxrows, nbuff, nGood, inloc, outloc;
176 LONGLONG ntodo, inbyteloc, outbyteloc, hsize;
177 long freespace;
178 unsigned char *buffer, result;
179 struct {
180 LONGLONG rowLength, numRows, heapSize;
181 LONGLONG dataStart, heapStart;
182 } inExt, outExt;
183
184 if( *status ) return( *status );
185
186 FFLOCK;
187 if( ffiprs( infptr, 0, expr, MAXDIMS, &Info.datatype, &nelem, &naxis,
188 naxes, status ) ) {
189 ffcprs();
190 FFUNLOCK;
191 return( *status );
192 }
193
194 if( nelem<0 ) {
195 constant = 1;
196 nelem = -nelem;
197 } else
198 constant = 0;
199
200 /**********************************************************************/
201 /* Make sure expression evaluates to the right type... logical scalar */
202 /**********************************************************************/
203
204 if( Info.datatype!=TLOGICAL || nelem!=1 ) {
205 ffcprs();
206 ffpmsg("Expression does not evaluate to a logical scalar.");
207 FFUNLOCK;
208 return( *status = PARSE_BAD_TYPE );
209 }
210
211 /***********************************************************/
212 /* Extract various table information from each extension */
213 /***********************************************************/
214
215 if( infptr->HDUposition != (infptr->Fptr)->curhdu )
216 ffmahd( infptr, (infptr->HDUposition) + 1, NULL, status );
217 if( *status ) {
218 ffcprs();
219 FFUNLOCK;
220 return( *status );
221 }
222 inExt.rowLength = (long) (infptr->Fptr)->rowlength;
223 inExt.numRows = (infptr->Fptr)->numrows;
224 inExt.heapSize = (infptr->Fptr)->heapsize;
225 if( inExt.numRows == 0 ) { /* Nothing to copy */
226 ffcprs();
227 FFUNLOCK;
228 return( *status );
229 }
230
231 if( outfptr->HDUposition != (outfptr->Fptr)->curhdu )
232 ffmahd( outfptr, (outfptr->HDUposition) + 1, NULL, status );
233 if( (outfptr->Fptr)->datastart < 0 )
234 ffrdef( outfptr, status );
235 if( *status ) {
236 ffcprs();
237 FFUNLOCK;
238 return( *status );
239 }
240 outExt.rowLength = (long) (outfptr->Fptr)->rowlength;
241 outExt.numRows = (outfptr->Fptr)->numrows;
242 if( !outExt.numRows )
243 (outfptr->Fptr)->heapsize = 0L;
244 outExt.heapSize = (outfptr->Fptr)->heapsize;
245
246 if( inExt.rowLength != outExt.rowLength ) {
247 ffpmsg("Output table has different row length from input");
248 ffcprs();
249 FFUNLOCK;
250 return( *status = PARSE_BAD_OUTPUT );
251 }
252
253 /***********************************/
254 /* Fill out Info data for parser */
255 /***********************************/
256
257 Info.dataPtr = (char *)malloc( (size_t) ((inExt.numRows + 1) * sizeof(char)) );
258 Info.nullPtr = NULL;
259 Info.maxRows = (long) inExt.numRows;
260 if( !Info.dataPtr ) {
261 ffpmsg("Unable to allocate memory for row selection");
262 ffcprs();
263 FFUNLOCK;
264 return( *status = MEMORY_ALLOCATION );
265 }
266
267 /* make sure array is zero terminated */
268 ((char*)Info.dataPtr)[inExt.numRows] = 0;
269
270 if( constant ) { /* Set all rows to the same value from constant result */
271
272 result = gParse.Nodes[gParse.resultNode].value.data.log;
273 for( ntodo = 0; ntodo<inExt.numRows; ntodo++ )
274 ((char*)Info.dataPtr)[ntodo] = result;
275 nGood = (long) (result ? inExt.numRows : 0);
276
277 } else {
278
279 ffiter( gParse.nCols, gParse.colData, 0L, 0L,
280 parse_data, (void*)&Info, status );
281
282 nGood = 0;
283 for( ntodo = 0; ntodo<inExt.numRows; ntodo++ )
284 if( ((char*)Info.dataPtr)[ntodo] ) nGood++;
285 }
286
287 if( *status ) {
288 /* Error... Do nothing */
289 } else {
290 rdlen = (long) inExt.rowLength;
291 buffer = (unsigned char *)malloc(maxvalue(500000,rdlen) * sizeof(char) );
292 if( buffer==NULL ) {
293 ffcprs();
294 FFUNLOCK;
295 return( *status=MEMORY_ALLOCATION );
296 }
297 maxrows = maxvalue( (500000L/rdlen), 1);
298 nbuff = 0;
299 inloc = 1;
300 if( infptr==outfptr ) { /* Skip initial good rows if input==output file */
301 while( ((char*)Info.dataPtr)[inloc-1] ) inloc++;
302 outloc = inloc;
303 } else {
304 outloc = (long) (outExt.numRows + 1);
305 if (outloc > 1)
306 ffirow( outfptr, outExt.numRows, nGood, status );
307 }
308
309 do {
310 if( ((char*)Info.dataPtr)[inloc-1] ) {
311 ffgtbb( infptr, inloc, 1L, rdlen, buffer+rdlen*nbuff, status );
312 nbuff++;
313 if( nbuff==maxrows ) {
314 ffptbb( outfptr, outloc, 1L, rdlen*nbuff, buffer, status );
315 outloc += nbuff;
316 nbuff = 0;
317 }
318 }
319 inloc++;
320 } while( !*status && inloc<=inExt.numRows );
321
322 if( nbuff ) {
323 ffptbb( outfptr, outloc, 1L, rdlen*nbuff, buffer, status );
324 outloc += nbuff;
325 }
326
327 if( infptr==outfptr ) {
328
329 if( outloc<=inExt.numRows )
330 ffdrow( infptr, outloc, inExt.numRows-outloc+1, status );
331
332 } else if( inExt.heapSize && nGood ) {
333
334 /* Copy heap, if it exists and at least one row copied */
335
336 /********************************************************/
337 /* Get location information from the output extension */
338 /********************************************************/
339
340 if( outfptr->HDUposition != (outfptr->Fptr)->curhdu )
341 ffmahd( outfptr, (outfptr->HDUposition) + 1, NULL, status );
342 outExt.dataStart = (outfptr->Fptr)->datastart;
343 outExt.heapStart = (outfptr->Fptr)->heapstart;
344
345 /*************************************************/
346 /* Insert more space into outfptr if necessary */
347 /*************************************************/
348
349 hsize = outExt.heapStart + outExt.heapSize;
350 freespace = (long) (( ( (hsize + 2879) / 2880) * 2880) - hsize);
351 ntodo = inExt.heapSize;
352
353 if ( (freespace - ntodo) < 0) { /* not enough existing space? */
354 ntodo = (ntodo - freespace + 2879) / 2880; /* number of blocks */
355 ffiblk(outfptr, (long) ntodo, 1, status); /* insert the blocks */
356 }
357 ffukyj( outfptr, "PCOUNT", inExt.heapSize+outExt.heapSize,
358 NULL, status );
359
360 /*******************************************************/
361 /* Get location information from the input extension */
362 /*******************************************************/
363
364 if( infptr->HDUposition != (infptr->Fptr)->curhdu )
365 ffmahd( infptr, (infptr->HDUposition) + 1, NULL, status );
366 inExt.dataStart = (infptr->Fptr)->datastart;
367 inExt.heapStart = (infptr->Fptr)->heapstart;
368
369 /**********************************/
370 /* Finally copy heap to outfptr */
371 /**********************************/
372
373 ntodo = inExt.heapSize;
374 inbyteloc = inExt.heapStart + inExt.dataStart;
375 outbyteloc = outExt.heapStart + outExt.dataStart + outExt.heapSize;
376
377 while ( ntodo && !*status ) {
378 rdlen = (long) minvalue(ntodo,500000);
379 ffmbyt( infptr, inbyteloc, REPORT_EOF, status );
380 ffgbyt( infptr, rdlen, buffer, status );
381 ffmbyt( outfptr, outbyteloc, IGNORE_EOF, status );
382 ffpbyt( outfptr, rdlen, buffer, status );
383 inbyteloc += rdlen;
384 outbyteloc += rdlen;
385 ntodo -= rdlen;
386 }
387
388 /***********************************************************/
389 /* But must update DES if data is being appended to a */
390 /* pre-existing heap space. Edit each new entry in file */
391 /***********************************************************/
392
393 if( outExt.heapSize ) {
394 LONGLONG repeat, offset, j;
395 int i;
396 for( i=1; i<=(outfptr->Fptr)->tfield; i++ ) {
397 if( (outfptr->Fptr)->tableptr[i-1].tdatatype<0 ) {
398 for( j=outExt.numRows+1; j<=outExt.numRows+nGood; j++ ) {
399 ffgdesll( outfptr, i, j, &repeat, &offset, status );
400 offset += outExt.heapSize;
401 ffpdes( outfptr, i, j, repeat, offset, status );
402 }
403 }
404 }
405 }
406
407 } /* End of HEAP copy */
408
409 FREE(buffer);
410 }
411
412 FREE(Info.dataPtr);
413 ffcprs();
414
415 ffcmph(outfptr, status); /* compress heap, deleting any orphaned data */
416 FFUNLOCK;
417 return(*status);
418 }
419
420 /*---------------------------------------------------------------------------*/
ffcrow(fitsfile * fptr,int datatype,char * expr,long firstrow,long nelements,void * nulval,void * array,int * anynul,int * status)421 int ffcrow( fitsfile *fptr, /* I - Input FITS file */
422 int datatype, /* I - Datatype to return results as */
423 char *expr, /* I - Arithmetic expression */
424 long firstrow, /* I - First row to evaluate */
425 long nelements, /* I - Number of elements to return */
426 void *nulval, /* I - Ptr to value to use as UNDEF */
427 void *array, /* O - Array of results */
428 int *anynul, /* O - Were any UNDEFs encountered? */
429 int *status ) /* O - Error status */
430 /* */
431 /* Calculate an expression for the indicated rows of a table, returning */
432 /* the results, cast as datatype (TSHORT, TDOUBLE, etc), in array. If */
433 /* nulval==NULL, UNDEFs will be zeroed out. For vector results, the number */
434 /* of elements returned may be less than nelements if nelements is not an */
435 /* even multiple of the result dimension. Call fftexp to obtain the */
436 /* dimensions of the results. */
437 /*---------------------------------------------------------------------------*/
438 {
439 parseInfo Info;
440 int naxis;
441 long nelem1, naxes[MAXDIMS];
442
443 if( *status ) return( *status );
444
445 FFLOCK;
446 if( ffiprs( fptr, 0, expr, MAXDIMS, &Info.datatype, &nelem1, &naxis,
447 naxes, status ) ) {
448 ffcprs();
449 FFUNLOCK;
450 return( *status );
451 }
452 if( nelem1<0 ) nelem1 = - nelem1;
453
454 if( nelements<nelem1 ) {
455 ffcprs();
456 ffpmsg("Array not large enough to hold at least one row of data.");
457 FFUNLOCK;
458 return( *status = PARSE_LRG_VECTOR );
459 }
460
461 firstrow = (firstrow>1 ? firstrow : 1);
462
463 if( datatype ) Info.datatype = datatype;
464
465 Info.dataPtr = array;
466 Info.nullPtr = nulval;
467 Info.maxRows = nelements / nelem1;
468
469 if( ffiter( gParse.nCols, gParse.colData, firstrow-1, 0,
470 parse_data, (void*)&Info, status ) == -1 )
471 *status=0; /* -1 indicates exitted without error before end... OK */
472
473 *anynul = Info.anyNull;
474 ffcprs();
475 FFUNLOCK;
476 return( *status );
477 }
478
479 /*--------------------------------------------------------------------------*/
ffcalc(fitsfile * infptr,char * expr,fitsfile * outfptr,char * parName,char * parInfo,int * status)480 int ffcalc( fitsfile *infptr, /* I - Input FITS file */
481 char *expr, /* I - Arithmetic expression */
482 fitsfile *outfptr, /* I - Output fits file */
483 char *parName, /* I - Name of output parameter */
484 char *parInfo, /* I - Extra information on parameter */
485 int *status ) /* O - Error status */
486 /* */
487 /* Evaluate an expression for all rows of a table. Call ffcalc_rng with */
488 /* a row range of 1-MAX. */
489 {
490 long start=1, end=LONG_MAX;
491
492 return ffcalc_rng( infptr, expr, outfptr, parName, parInfo,
493 1, &start, &end, status );
494 }
495
496 /*--------------------------------------------------------------------------*/
ffcalc_rng(fitsfile * infptr,char * expr,fitsfile * outfptr,char * parName,char * parInfo,int nRngs,long * start,long * end,int * status)497 int ffcalc_rng( fitsfile *infptr, /* I - Input FITS file */
498 char *expr, /* I - Arithmetic expression */
499 fitsfile *outfptr, /* I - Output fits file */
500 char *parName, /* I - Name of output parameter */
501 char *parInfo, /* I - Extra information on parameter */
502 int nRngs, /* I - Row range info */
503 long *start, /* I - Row range info */
504 long *end, /* I - Row range info */
505 int *status ) /* O - Error status */
506 /* */
507 /* Evaluate an expression using the data in the input FITS file and place */
508 /* the results into either a column or keyword in the output fits file, */
509 /* depending on the value of parName (keywords normally prefixed with '#') */
510 /* and whether the expression evaluates to a constant or a table column. */
511 /* The logic is as follows: */
512 /* (1) If a column exists with name, parName, put results there. */
513 /* (2) If parName starts with '#', as in #NAXIS, put result there, */
514 /* with parInfo used as the comment. If expression does not evaluate */
515 /* to a constant, flag an error. */
516 /* (3) If a keyword exists with name, parName, and expression is a */
517 /* constant, put result there, using parInfo as the new comment. */
518 /* (4) Else, create a new column with name parName and TFORM parInfo. */
519 /* If parInfo is NULL, use a default data type for the column. */
520 /*--------------------------------------------------------------------------*/
521 {
522 parseInfo Info;
523 int naxis, constant, typecode, newNullKwd=0;
524 long nelem, naxes[MAXDIMS], repeat, width;
525 int col_cnt, colNo;
526 Node *result;
527 char card[81], tform[16], nullKwd[9], tdimKwd[9];
528
529 if( *status ) return( *status );
530
531 FFLOCK;
532 if( ffiprs( infptr, 0, expr, MAXDIMS, &Info.datatype, &nelem, &naxis,
533 naxes, status ) ) {
534
535 ffcprs();
536 FFUNLOCK;
537 return( *status );
538 }
539 if( nelem<0 ) {
540 constant = 1;
541 nelem = -nelem;
542 } else
543 constant = 0;
544
545 /* Case (1): If column exists put it there */
546
547 colNo = 0;
548 if( ffgcno( outfptr, CASEINSEN, parName, &colNo, status )==COL_NOT_FOUND ) {
549
550 /* Output column doesn't exist. Test for keyword. */
551
552 /* Case (2): Does parName indicate result should be put into keyword */
553
554 *status = 0;
555 if( parName[0]=='#' ) {
556 if( ! constant ) {
557 ffcprs();
558 ffpmsg( "Cannot put tabular result into keyword (ffcalc)" );
559 FFUNLOCK;
560 return( *status = PARSE_BAD_TYPE );
561 }
562 parName++; /* Advance past '#' */
563 if ( (fits_strcasecmp(parName,"HISTORY") == 0 || fits_strcasecmp(parName,"COMMENT") == 0) &&
564 Info.datatype != TSTRING ) {
565 ffcprs();
566 ffpmsg( "HISTORY and COMMENT values must be strings (ffcalc)" );
567 FFUNLOCK;
568 return( *status = PARSE_BAD_TYPE );
569 }
570
571 } else if( constant ) {
572
573 /* Case (3): Does a keyword named parName already exist */
574
575 if( ffgcrd( outfptr, parName, card, status )==KEY_NO_EXIST ) {
576 colNo = -1;
577 } else if( *status ) {
578 ffcprs();
579 FFUNLOCK;
580 return( *status );
581 }
582
583 } else
584 colNo = -1;
585
586 if( colNo<0 ) {
587
588 /* Case (4): Create new column */
589
590 *status = 0;
591 ffgncl( outfptr, &colNo, status );
592 colNo++;
593 if( parInfo==NULL || *parInfo=='\0' ) {
594 /* Figure out best default column type */
595 if( gParse.hdutype==BINARY_TBL ) {
596 snprintf(tform,15,"%ld",nelem);
597 switch( Info.datatype ) {
598 case TLOGICAL: strcat(tform,"L"); break;
599 case TLONG: strcat(tform,"J"); break;
600 case TDOUBLE: strcat(tform,"D"); break;
601 case TSTRING: strcat(tform,"A"); break;
602 case TBIT: strcat(tform,"X"); break;
603 case TLONGLONG: strcat(tform,"K"); break;
604 }
605 } else {
606 switch( Info.datatype ) {
607 case TLOGICAL:
608 ffcprs();
609 ffpmsg("Cannot create LOGICAL column in ASCII table");
610 FFUNLOCK;
611 return( *status = NOT_BTABLE );
612 case TLONG: strcpy(tform,"I11"); break;
613 case TDOUBLE: strcpy(tform,"D23.15"); break;
614 case TSTRING:
615 case TBIT: snprintf(tform,16,"A%ld",nelem); break;
616 }
617 }
618 parInfo = tform;
619 } else if( !(isdigit((int) *parInfo)) && gParse.hdutype==BINARY_TBL ) {
620 if( Info.datatype==TBIT && *parInfo=='B' )
621 nelem = (nelem+7)/8;
622 snprintf(tform,16,"%ld%s",nelem,parInfo);
623 parInfo = tform;
624 }
625 fficol( outfptr, colNo, parName, parInfo, status );
626 if( naxis>1 )
627 ffptdm( outfptr, colNo, naxis, naxes, status );
628
629 /* Setup TNULLn keyword in case NULLs are encountered */
630
631 ffkeyn("TNULL", colNo, nullKwd, status);
632 if( ffgcrd( outfptr, nullKwd, card, status )==KEY_NO_EXIST ) {
633 *status = 0;
634 if( gParse.hdutype==BINARY_TBL ) {
635 LONGLONG nullVal=0;
636 fits_binary_tform( parInfo, &typecode, &repeat, &width, status );
637 if( typecode==TBYTE )
638 nullVal = UCHAR_MAX;
639 else if( typecode==TSHORT )
640 nullVal = SHRT_MIN;
641 else if( typecode==TINT )
642 nullVal = INT_MIN;
643 else if( typecode==TLONG ) {
644 if (sizeof(long) == 8 && sizeof(int) == 4)
645 nullVal = INT_MIN;
646 else
647 nullVal = LONG_MIN;
648 }
649 else if( typecode==TLONGLONG )
650 nullVal = LONGLONG_MIN;
651
652 if( nullVal ) {
653 ffpkyj( outfptr, nullKwd, nullVal, "Null value", status );
654 fits_set_btblnull( outfptr, colNo, nullVal, status );
655 newNullKwd = 1;
656 }
657 } else if( gParse.hdutype==ASCII_TBL ) {
658 ffpkys( outfptr, nullKwd, "NULL", "Null value string", status );
659 fits_set_atblnull( outfptr, colNo, "NULL", status );
660 newNullKwd = 1;
661 }
662 }
663
664 }
665
666 } else if( *status ) {
667 ffcprs();
668 FFUNLOCK;
669 return( *status );
670 } else {
671
672 /********************************************************/
673 /* Check if a TDIM keyword should be written/updated. */
674 /********************************************************/
675
676 ffkeyn("TDIM", colNo, tdimKwd, status);
677 ffgcrd( outfptr, tdimKwd, card, status );
678 if( *status==0 ) {
679 /* TDIM exists, so update it with result's dimension */
680 ffptdm( outfptr, colNo, naxis, naxes, status );
681 } else if( *status==KEY_NO_EXIST ) {
682 /* TDIM does not exist, so clear error stack and */
683 /* write a TDIM only if result is multi-dimensional */
684 *status = 0;
685 ffcmsg();
686 if( naxis>1 )
687 ffptdm( outfptr, colNo, naxis, naxes, status );
688 }
689 if( *status ) {
690 /* Either some other error happened in ffgcrd */
691 /* or one happened in ffptdm */
692 ffcprs();
693 FFUNLOCK;
694 return( *status );
695 }
696
697 }
698
699 if( colNo>0 ) {
700
701 /* Output column exists (now)... put results into it */
702
703 int anyNull = 0;
704 int nPerLp, i;
705 long totaln;
706
707 ffgkyj(infptr, "NAXIS2", &totaln, 0, status);
708
709 /*************************************/
710 /* Create new iterator Output Column */
711 /*************************************/
712
713 col_cnt = gParse.nCols;
714 if( allocateCol( col_cnt, status ) ) {
715 ffcprs();
716 FFUNLOCK;
717 return( *status );
718 }
719
720 fits_iter_set_by_num( gParse.colData+col_cnt, outfptr,
721 colNo, 0, OutputCol );
722 gParse.nCols++;
723
724 for( i=0; i<nRngs; i++ ) {
725 Info.dataPtr = NULL;
726 Info.maxRows = end[i]-start[i]+1;
727
728 /*
729 If there is only 1 range, and it includes all the rows,
730 and there are 10 or more rows, then set nPerLp = 0 so
731 that the iterator function will dynamically choose the
732 most efficient number of rows to process in each loop.
733 Otherwise, set nPerLp to the number of rows in this range.
734 */
735
736 if( (Info.maxRows >= 10) && (nRngs == 1) &&
737 (start[0] == 1) && (end[0] == totaln))
738 nPerLp = 0;
739 else
740 nPerLp = Info.maxRows;
741
742 if( ffiter( gParse.nCols, gParse.colData, start[i]-1,
743 nPerLp, parse_data, (void*)&Info, status ) == -1 )
744 *status = 0;
745 else if( *status ) {
746 ffcprs();
747 FFUNLOCK;
748 return( *status );
749 }
750 if( Info.anyNull ) anyNull = 1;
751 }
752
753 if( newNullKwd && !anyNull ) {
754 ffdkey( outfptr, nullKwd, status );
755 }
756
757 } else {
758
759 /* Put constant result into keyword */
760
761 result = gParse.Nodes + gParse.resultNode;
762 switch( Info.datatype ) {
763 case TDOUBLE:
764 ffukyd( outfptr, parName, result->value.data.dbl, 15,
765 parInfo, status );
766 break;
767 case TLONG:
768 ffukyj( outfptr, parName, result->value.data.lng, parInfo, status );
769 break;
770 case TLOGICAL:
771 ffukyl( outfptr, parName, result->value.data.log, parInfo, status );
772 break;
773 case TBIT:
774 case TSTRING:
775 if (fits_strcasecmp(parName,"HISTORY") == 0) {
776 ffphis( outfptr, result->value.data.str, status);
777 } else if (fits_strcasecmp(parName,"COMMENT") == 0) {
778 ffpcom( outfptr, result->value.data.str, status);
779 } else {
780 ffukys( outfptr, parName, result->value.data.str, parInfo, status );
781 }
782 break;
783 }
784 }
785
786 ffcprs();
787 FFUNLOCK;
788 return( *status );
789 }
790
791 /*--------------------------------------------------------------------------*/
fftexp(fitsfile * fptr,char * expr,int maxdim,int * datatype,long * nelem,int * naxis,long * naxes,int * status)792 int fftexp( fitsfile *fptr, /* I - Input FITS file */
793 char *expr, /* I - Arithmetic expression */
794 int maxdim, /* I - Max Dimension of naxes */
795 int *datatype, /* O - Data type of result */
796 long *nelem, /* O - Vector length of result */
797 int *naxis, /* O - # of dimensions of result */
798 long *naxes, /* O - Size of each dimension */
799 int *status ) /* O - Error status */
800 /* */
801 /* Evaluate the given expression and return information on the result. */
802 /*--------------------------------------------------------------------------*/
803 {
804 FFLOCK;
805 ffiprs( fptr, 0, expr, maxdim, datatype, nelem, naxis, naxes, status );
806 ffcprs();
807 FFUNLOCK;
808 return( *status );
809 }
810
811 /*--------------------------------------------------------------------------*/
ffiprs(fitsfile * fptr,int compressed,char * expr,int maxdim,int * datatype,long * nelem,int * naxis,long * naxes,int * status)812 int ffiprs( fitsfile *fptr, /* I - Input FITS file */
813 int compressed, /* I - Is FITS file hkunexpanded? */
814 char *expr, /* I - Arithmetic expression */
815 int maxdim, /* I - Max Dimension of naxes */
816 int *datatype, /* O - Data type of result */
817 long *nelem, /* O - Vector length of result */
818 int *naxis, /* O - # of dimensions of result */
819 long *naxes, /* O - Size of each dimension */
820 int *status ) /* O - Error status */
821 /* */
822 /* Initialize the parser and determine what type of result the expression */
823 /* produces. */
824 /*--------------------------------------------------------------------------*/
825 {
826 Node *result;
827 int i,lexpr, tstatus = 0;
828 int xaxis, bitpix;
829 long xaxes[9];
830 static iteratorCol dmyCol;
831
832 if( *status ) return( *status );
833
834 /* make sure all internal structures for this HDU are current */
835 if ( ffrdef(fptr, status) ) return(*status);
836
837 /* Initialize the Parser structure */
838
839 gParse.def_fptr = fptr;
840 gParse.compressed = compressed;
841 gParse.nCols = 0;
842 gParse.colData = NULL;
843 gParse.varData = NULL;
844 gParse.getData = find_column;
845 gParse.loadData = load_column;
846 gParse.Nodes = NULL;
847 gParse.nNodesAlloc= 0;
848 gParse.nNodes = 0;
849 gParse.hdutype = 0;
850 gParse.status = 0;
851
852 fits_get_hdu_type(fptr, &gParse.hdutype, status );
853
854 if (gParse.hdutype == IMAGE_HDU) {
855
856 fits_get_img_param(fptr, 9, &bitpix, &xaxis, xaxes, status);
857 if (*status) {
858 ffpmsg("ffiprs: unable to get image dimensions");
859 return( *status );
860 }
861 gParse.totalRows = xaxis > 0 ? 1 : 0;
862 for (i = 0; i < xaxis; ++i)
863 gParse.totalRows *= xaxes[i];
864 if (DEBUG_PIXFILTER)
865 printf("naxis=%d, gParse.totalRows=%ld\n", xaxis, gParse.totalRows);
866 }
867 else if( ffgkyj(fptr, "NAXIS2", &gParse.totalRows, 0, &tstatus) )
868 {
869 /* this might be a 1D or null image with no NAXIS2 keyword */
870 gParse.totalRows = 0;
871 }
872
873
874 /* Copy expression into parser... read from file if necessary */
875
876
877 if( expr[0]=='@' ) {
878 if( ffimport_file( expr+1, &gParse.expr, status ) ) return( *status );
879 lexpr = strlen(gParse.expr);
880 } else {
881 lexpr = strlen(expr);
882 gParse.expr = (char*)malloc( (2+lexpr)*sizeof(char));
883 strcpy(gParse.expr,expr);
884 }
885 strcat(gParse.expr + lexpr,"\n");
886 gParse.index = 0;
887 gParse.is_eobuf = 0;
888
889 /* Parse the expression, building the Nodes and determing */
890 /* which columns are needed and what data type is returned */
891
892 ffrestart(NULL);
893 if( ffparse() ) {
894 return( *status = PARSE_SYNTAX_ERR );
895 }
896 /* Check results */
897
898 *status = gParse.status;
899 if( *status ) return(*status);
900
901 if( !gParse.nNodes ) {
902 ffpmsg("Blank expression");
903 return( *status = PARSE_SYNTAX_ERR );
904 }
905 if( !gParse.nCols ) {
906 dmyCol.fptr = fptr; /* This allows iterator to know value of */
907 gParse.colData = &dmyCol; /* fptr when no columns are referenced */
908 }
909
910 result = gParse.Nodes + gParse.resultNode;
911
912 *naxis = result->value.naxis;
913 *nelem = result->value.nelem;
914 for( i=0; i<*naxis && i<maxdim; i++ )
915 naxes[i] = result->value.naxes[i];
916
917 switch( result->type ) {
918 case BOOLEAN:
919 *datatype = TLOGICAL;
920 break;
921 case LONG:
922 *datatype = TLONG;
923 break;
924 case DOUBLE:
925 *datatype = TDOUBLE;
926 break;
927 case BITSTR:
928 *datatype = TBIT;
929 break;
930 case STRING:
931 *datatype = TSTRING;
932 break;
933 default:
934 *datatype = 0;
935 ffpmsg("Bad return data type");
936 *status = gParse.status = PARSE_BAD_TYPE;
937 break;
938 }
939 gParse.datatype = *datatype;
940 FREE(gParse.expr);
941
942 if( result->operation==CONST_OP ) *nelem = - *nelem;
943 return(*status);
944 }
945
946 /*--------------------------------------------------------------------------*/
ffcprs(void)947 void ffcprs( void ) /* No parameters */
948 /* */
949 /* Clear the parser, making it ready to accept a new expression. */
950 /*--------------------------------------------------------------------------*/
951 {
952 int col, node, i;
953
954 if( gParse.nCols > 0 ) {
955 FREE( gParse.colData );
956 for( col=0; col<gParse.nCols; col++ ) {
957 if( gParse.varData[col].undef == NULL ) continue;
958 if( gParse.varData[col].type == BITSTR )
959 FREE( ((char**)gParse.varData[col].data)[0] );
960 free( gParse.varData[col].undef );
961 }
962 FREE( gParse.varData );
963 gParse.nCols = 0;
964 }
965
966 if( gParse.nNodes > 0 ) {
967 node = gParse.nNodes;
968 while( node-- ) {
969 if( gParse.Nodes[node].operation==gtifilt_fct ) {
970 i = gParse.Nodes[node].SubNodes[0];
971 if (gParse.Nodes[ i ].value.data.ptr)
972 FREE( gParse.Nodes[ i ].value.data.ptr );
973 }
974 else if( gParse.Nodes[node].operation==regfilt_fct ) {
975 i = gParse.Nodes[node].SubNodes[0];
976 fits_free_region( (SAORegion *)gParse.Nodes[ i ].value.data.ptr );
977 }
978 }
979 gParse.nNodes = 0;
980 }
981 if( gParse.Nodes ) free( gParse.Nodes );
982 gParse.Nodes = NULL;
983
984 gParse.hdutype = ANY_HDU;
985 gParse.pixFilter = 0;
986 }
987
988 /*---------------------------------------------------------------------------*/
parse_data(long totalrows,long offset,long firstrow,long nrows,int nCols,iteratorCol * colData,void * userPtr)989 int parse_data( long totalrows, /* I - Total rows to be processed */
990 long offset, /* I - Number of rows skipped at start*/
991 long firstrow, /* I - First row of this iteration */
992 long nrows, /* I - Number of rows in this iter */
993 int nCols, /* I - Number of columns in use */
994 iteratorCol *colData, /* IO- Column information/data */
995 void *userPtr ) /* I - Data handling instructions */
996 /* */
997 /* Iterator work function which calls the parser and copies the results */
998 /* into either an OutputCol or a data pointer supplied in the userPtr */
999 /* structure. */
1000 /*---------------------------------------------------------------------------*/
1001 {
1002 int status, constant=0, anyNullThisTime=0;
1003 long jj, kk, idx, remain, ntodo;
1004 Node *result;
1005 iteratorCol * outcol;
1006
1007 /* declare variables static to preserve their values between calls */
1008 static void *Data, *Null;
1009 static int datasize;
1010 static long lastRow, repeat, resDataSize;
1011 static LONGLONG jnull;
1012 static parseInfo *userInfo;
1013 static long zeros[4] = {0,0,0,0};
1014
1015 if (DEBUG_PIXFILTER)
1016 printf("parse_data(total=%ld, offset=%ld, first=%ld, rows=%ld, cols=%d)\n",
1017 totalrows, offset, firstrow, nrows, nCols);
1018 /*--------------------------------------------------------*/
1019 /* Initialization procedures: execute on the first call */
1020 /*--------------------------------------------------------*/
1021 outcol = colData + (nCols - 1);
1022 if (firstrow == offset+1)
1023 {
1024 userInfo = (parseInfo*)userPtr;
1025 userInfo->anyNull = 0;
1026
1027 if( userInfo->maxRows>0 )
1028 userInfo->maxRows = minvalue(totalrows,userInfo->maxRows);
1029 else if( userInfo->maxRows<0 )
1030 userInfo->maxRows = totalrows;
1031 else
1032 userInfo->maxRows = nrows;
1033
1034 lastRow = firstrow + userInfo->maxRows - 1;
1035
1036 if( userInfo->dataPtr==NULL ) {
1037
1038 if( outcol->iotype == InputCol ) {
1039 ffpmsg("Output column for parser results not found!");
1040 return( PARSE_NO_OUTPUT );
1041 }
1042 /* Data gets set later */
1043 Null = outcol->array;
1044 userInfo->datatype = outcol->datatype;
1045
1046 /* Check for a TNULL/BLANK keyword for output column/image */
1047
1048 status = 0;
1049 jnull = 0;
1050 if (gParse.hdutype == IMAGE_HDU) {
1051 if (gParse.pixFilter->blank)
1052 jnull = (LONGLONG) gParse.pixFilter->blank;
1053 }
1054 else {
1055 ffgknjj( outcol->fptr, "TNULL", outcol->colnum,
1056 1, &jnull, (int*)&jj, &status );
1057
1058 if( status==BAD_INTKEY ) {
1059 /* Probably ASCII table with text TNULL keyword */
1060 switch( userInfo->datatype ) {
1061 case TSHORT: jnull = (LONGLONG) SHRT_MIN; break;
1062 case TINT: jnull = (LONGLONG) INT_MIN; break;
1063 case TLONG: jnull = (LONGLONG) LONG_MIN; break;
1064 }
1065 }
1066 }
1067 repeat = outcol->repeat;
1068 /*
1069 if (DEBUG_PIXFILTER)
1070 printf("parse_data: using null value %ld\n", jnull);
1071 */
1072 } else {
1073
1074 Data = userInfo->dataPtr;
1075 Null = (userInfo->nullPtr ? userInfo->nullPtr : zeros);
1076 repeat = gParse.Nodes[gParse.resultNode].value.nelem;
1077
1078 }
1079
1080 /* Determine the size of each element of the returned result */
1081
1082 switch( userInfo->datatype ) {
1083 case TBIT: /* Fall through to TBYTE */
1084 case TLOGICAL: /* Fall through to TBYTE */
1085 case TBYTE: datasize = sizeof(char); break;
1086 case TSHORT: datasize = sizeof(short); break;
1087 case TINT: datasize = sizeof(int); break;
1088 case TLONG: datasize = sizeof(long); break;
1089 case TLONGLONG: datasize = sizeof(LONGLONG); break;
1090 case TFLOAT: datasize = sizeof(float); break;
1091 case TDOUBLE: datasize = sizeof(double); break;
1092 case TSTRING: datasize = sizeof(char*); break;
1093 }
1094
1095 /* Determine the size of each element of the calculated result */
1096 /* (only matters for numeric/logical data) */
1097
1098 switch( gParse.Nodes[gParse.resultNode].type ) {
1099 case BOOLEAN: resDataSize = sizeof(char); break;
1100 case LONG: resDataSize = sizeof(long); break;
1101 case DOUBLE: resDataSize = sizeof(double); break;
1102 }
1103 }
1104
1105 /*-------------------------------------------*/
1106 /* Main loop: process all the rows of data */
1107 /*-------------------------------------------*/
1108
1109 /* If writing to output column, set first element to appropriate */
1110 /* null value. If no NULLs encounter, zero out before returning. */
1111 /*
1112 if (DEBUG_PIXFILTER)
1113 printf("parse_data: using null value %ld\n", jnull);
1114 */
1115
1116 if( userInfo->dataPtr == NULL ) {
1117 /* First, reset Data pointer to start of output array */
1118 Data = (char*) outcol->array + datasize;
1119
1120 switch( userInfo->datatype ) {
1121 case TLOGICAL: *(char *)Null = 'U'; break;
1122 case TBYTE: *(char *)Null = (char )jnull; break;
1123 case TSHORT: *(short *)Null = (short)jnull; break;
1124 case TINT: *(int *)Null = (int )jnull; break;
1125 case TLONG: *(long *)Null = (long )jnull; break;
1126 case TLONGLONG: *(LONGLONG *)Null = (LONGLONG )jnull; break;
1127 case TFLOAT: *(float *)Null = FLOATNULLVALUE; break;
1128 case TDOUBLE: *(double*)Null = DOUBLENULLVALUE; break;
1129 case TSTRING: (*(char **)Null)[0] = '\1';
1130 (*(char **)Null)[1] = '\0'; break;
1131 }
1132 }
1133
1134 /* Alter nrows in case calling routine didn't want to do all rows */
1135
1136 nrows = minvalue(nrows,lastRow-firstrow+1);
1137
1138 Setup_DataArrays( nCols, colData, firstrow, nrows );
1139
1140 /* Parser allocates arrays for each column and calculation it performs. */
1141 /* Limit number of rows processed during each pass to reduce memory */
1142 /* requirements... In most cases, iterator will limit rows to less */
1143 /* than 2500 rows per iteration, so this is really only relevant for */
1144 /* hk-compressed files which must be decompressed in memory and sent */
1145 /* whole to parse_data in a single iteration. */
1146
1147 remain = nrows;
1148 while( remain ) {
1149 ntodo = minvalue(remain,2500);
1150 Evaluate_Parser ( firstrow, ntodo );
1151 if( gParse.status ) break;
1152
1153 firstrow += ntodo;
1154 remain -= ntodo;
1155
1156 /* Copy results into data array */
1157
1158 result = gParse.Nodes + gParse.resultNode;
1159 if( result->operation==CONST_OP ) constant = 1;
1160
1161 switch( result->type ) {
1162
1163 case BOOLEAN:
1164 case LONG:
1165 case DOUBLE:
1166 if( constant ) {
1167 char undef=0;
1168 for( kk=0; kk<ntodo; kk++ )
1169 for( jj=0; jj<repeat; jj++ )
1170 ffcvtn( gParse.datatype,
1171 &(result->value.data),
1172 &undef, result->value.nelem /* 1 */,
1173 userInfo->datatype, Null,
1174 (char*)Data + (kk*repeat+jj)*datasize,
1175 &anyNullThisTime, &gParse.status );
1176 } else {
1177 if ( repeat == result->value.nelem ) {
1178 ffcvtn( gParse.datatype,
1179 result->value.data.ptr,
1180 result->value.undef,
1181 result->value.nelem*ntodo,
1182 userInfo->datatype, Null, Data,
1183 &anyNullThisTime, &gParse.status );
1184 } else if( result->value.nelem == 1 ) {
1185 for( kk=0; kk<ntodo; kk++ )
1186 for( jj=0; jj<repeat; jj++ ) {
1187 ffcvtn( gParse.datatype,
1188 (char*)result->value.data.ptr + kk*resDataSize,
1189 (char*)result->value.undef + kk,
1190 1, userInfo->datatype, Null,
1191 (char*)Data + (kk*repeat+jj)*datasize,
1192 &anyNullThisTime, &gParse.status );
1193 }
1194 } else {
1195 int nCopy;
1196 nCopy = minvalue( repeat, result->value.nelem );
1197 for( kk=0; kk<ntodo; kk++ ) {
1198 ffcvtn( gParse.datatype,
1199 (char*)result->value.data.ptr
1200 + kk*result->value.nelem*resDataSize,
1201 (char*)result->value.undef
1202 + kk*result->value.nelem,
1203 nCopy, userInfo->datatype, Null,
1204 (char*)Data + (kk*repeat)*datasize,
1205 &anyNullThisTime, &gParse.status );
1206 if( nCopy < repeat ) {
1207 memset( (char*)Data + (kk*repeat+nCopy)*datasize,
1208 0, (repeat-nCopy)*datasize);
1209 }
1210 }
1211
1212 }
1213 if( result->operation>0 ) {
1214 FREE( result->value.data.ptr );
1215 }
1216 }
1217 if( gParse.status==OVERFLOW_ERR ) {
1218 gParse.status = NUM_OVERFLOW;
1219 ffpmsg("Numerical overflow while converting expression to necessary datatype");
1220 }
1221 break;
1222
1223 case BITSTR:
1224 switch( userInfo->datatype ) {
1225 case TBYTE:
1226 idx = -1;
1227 for( kk=0; kk<ntodo; kk++ ) {
1228 for( jj=0; jj<result->value.nelem; jj++ ) {
1229 if( jj%8 == 0 )
1230 ((char*)Data)[++idx] = 0;
1231 if( constant ) {
1232 if( result->value.data.str[jj]=='1' )
1233 ((char*)Data)[idx] |= 128>>(jj%8);
1234 } else {
1235 if( result->value.data.strptr[kk][jj]=='1' )
1236 ((char*)Data)[idx] |= 128>>(jj%8);
1237 }
1238 }
1239 }
1240 break;
1241 case TBIT:
1242 case TLOGICAL:
1243 if( constant ) {
1244 for( kk=0; kk<ntodo; kk++ )
1245 for( jj=0; jj<result->value.nelem; jj++ ) {
1246 ((char*)Data)[ jj+kk*result->value.nelem ] =
1247 ( result->value.data.str[jj]=='1' );
1248 }
1249 } else {
1250 for( kk=0; kk<ntodo; kk++ )
1251 for( jj=0; jj<result->value.nelem; jj++ ) {
1252 ((char*)Data)[ jj+kk*result->value.nelem ] =
1253 ( result->value.data.strptr[kk][jj]=='1' );
1254 }
1255 }
1256 break;
1257 case TSTRING:
1258 if( constant ) {
1259 for( jj=0; jj<ntodo; jj++ ) {
1260 strcpy( ((char**)Data)[jj], result->value.data.str );
1261 }
1262 } else {
1263 for( jj=0; jj<ntodo; jj++ ) {
1264 strcpy( ((char**)Data)[jj], result->value.data.strptr[jj] );
1265 }
1266 }
1267 break;
1268 default:
1269 ffpmsg("Cannot convert bit expression to desired type.");
1270 gParse.status = PARSE_BAD_TYPE;
1271 break;
1272 }
1273 if( result->operation>0 ) {
1274 FREE( result->value.data.strptr[0] );
1275 FREE( result->value.data.strptr );
1276 }
1277 break;
1278
1279 case STRING:
1280 if( userInfo->datatype==TSTRING ) {
1281 if( constant ) {
1282 for( jj=0; jj<ntodo; jj++ )
1283 strcpy( ((char**)Data)[jj], result->value.data.str );
1284 } else {
1285 for( jj=0; jj<ntodo; jj++ )
1286 if( result->value.undef[jj] ) {
1287 anyNullThisTime = 1;
1288 strcpy( ((char**)Data)[jj],
1289 *(char **)Null );
1290 } else {
1291 strcpy( ((char**)Data)[jj],
1292 result->value.data.strptr[jj] );
1293 }
1294 }
1295 } else {
1296 ffpmsg("Cannot convert string expression to desired type.");
1297 gParse.status = PARSE_BAD_TYPE;
1298 }
1299 if( result->operation>0 ) {
1300 FREE( result->value.data.strptr[0] );
1301 FREE( result->value.data.strptr );
1302 }
1303 break;
1304 }
1305
1306 if( gParse.status ) break;
1307
1308 /* Increment Data to point to where the next block should go */
1309
1310 if( result->type==BITSTR && userInfo->datatype==TBYTE )
1311 Data = (char*)Data
1312 + datasize * ( (result->value.nelem+7)/8 ) * ntodo;
1313 else if( result->type==STRING )
1314 Data = (char*)Data + datasize * ntodo;
1315 else
1316 Data = (char*)Data + datasize * ntodo * repeat;
1317 }
1318
1319 /* If no NULLs encountered during this pass, set Null value to */
1320 /* zero to make the writing of the output column data faster */
1321
1322 if( anyNullThisTime )
1323 userInfo->anyNull = 1;
1324 else if( userInfo->dataPtr == NULL ) {
1325 if( userInfo->datatype == TSTRING )
1326 memcpy( *(char **)Null, zeros, 2 );
1327 else
1328 memcpy( Null, zeros, datasize );
1329 }
1330
1331 /*-------------------------------------------------------*/
1332 /* Clean up procedures: after processing all the rows */
1333 /*-------------------------------------------------------*/
1334
1335 /* if the calling routine specified that only a limited number */
1336 /* of rows in the table should be processed, return a value of -1 */
1337 /* once all the rows have been done, if no other error occurred. */
1338
1339 if (gParse.hdutype != IMAGE_HDU && firstrow - 1 == lastRow) {
1340 if (!gParse.status && userInfo->maxRows<totalrows) {
1341 return (-1);
1342 }
1343 }
1344
1345 return(gParse.status); /* return successful status */
1346 }
1347
Setup_DataArrays(int nCols,iteratorCol * cols,long fRow,long nRows)1348 static void Setup_DataArrays( int nCols, iteratorCol *cols,
1349 long fRow, long nRows )
1350 /***********************************************************************/
1351 /* Setup the varData array in gParse to contain the fits column data. */
1352 /* Then, allocate and initialize the necessary UNDEF arrays for each */
1353 /* column used by the parser. */
1354 /***********************************************************************/
1355 {
1356 int i;
1357 long nelem, len, row, idx;
1358 char **bitStrs;
1359 char **sptr;
1360 char *barray;
1361 long *iarray;
1362 double *rarray;
1363 char msg[80];
1364
1365 gParse.firstDataRow = fRow;
1366 gParse.nDataRows = nRows;
1367
1368 /* Resize and fill in UNDEF arrays for each column */
1369
1370 for( i=0; i<nCols; i++ ) {
1371
1372 iteratorCol *icol = cols + i;
1373 DataInfo *varData = gParse.varData + i;
1374
1375 if( icol->iotype == OutputCol ) continue;
1376
1377 nelem = varData->nelem;
1378 len = nelem * nRows;
1379
1380 switch ( varData->type ) {
1381
1382 case BITSTR:
1383 /* No need for UNDEF array, but must make string DATA array */
1384 len = (nelem+1)*nRows; /* Count '\0' */
1385 bitStrs = (char**)varData->data;
1386 if( bitStrs ) FREE( bitStrs[0] );
1387 free( bitStrs );
1388 bitStrs = (char**)malloc( nRows*sizeof(char*) );
1389 if( bitStrs==NULL ) {
1390 varData->data = varData->undef = NULL;
1391 gParse.status = MEMORY_ALLOCATION;
1392 break;
1393 }
1394 bitStrs[0] = (char*)malloc( len*sizeof(char) );
1395 if( bitStrs[0]==NULL ) {
1396 free( bitStrs );
1397 varData->data = varData->undef = NULL;
1398 gParse.status = MEMORY_ALLOCATION;
1399 break;
1400 }
1401
1402 for( row=0; row<nRows; row++ ) {
1403 bitStrs[row] = bitStrs[0] + row*(nelem+1);
1404 idx = (row)*( (nelem+7)/8 ) + 1;
1405 for(len=0; len<nelem; len++) {
1406 if( ((char*)icol->array)[idx] & (1<<(7-len%8)) )
1407 bitStrs[row][len] = '1';
1408 else
1409 bitStrs[row][len] = '0';
1410 if( len%8==7 ) idx++;
1411 }
1412 bitStrs[row][len] = '\0';
1413 }
1414 varData->undef = (char*)bitStrs;
1415 varData->data = (char*)bitStrs;
1416 break;
1417
1418 case STRING:
1419 sptr = (char**)icol->array;
1420 if (varData->undef)
1421 free( varData->undef );
1422 varData->undef = (char*)malloc( nRows*sizeof(char) );
1423 if( varData->undef==NULL ) {
1424 gParse.status = MEMORY_ALLOCATION;
1425 break;
1426 }
1427 row = nRows;
1428 while( row-- )
1429 varData->undef[row] =
1430 ( **sptr != '\0' && FSTRCMP( sptr[0], sptr[row+1] )==0 );
1431 varData->data = sptr + 1;
1432 break;
1433
1434 case BOOLEAN:
1435 barray = (char*)icol->array;
1436 if (varData->undef)
1437 free( varData->undef );
1438 varData->undef = (char*)malloc( len*sizeof(char) );
1439 if( varData->undef==NULL ) {
1440 gParse.status = MEMORY_ALLOCATION;
1441 break;
1442 }
1443 while( len-- ) {
1444 varData->undef[len] =
1445 ( barray[0]!=0 && barray[0]==barray[len+1] );
1446 }
1447 varData->data = barray + 1;
1448 break;
1449
1450 case LONG:
1451 iarray = (long*)icol->array;
1452 if (varData->undef)
1453 free( varData->undef );
1454 varData->undef = (char*)malloc( len*sizeof(char) );
1455 if( varData->undef==NULL ) {
1456 gParse.status = MEMORY_ALLOCATION;
1457 break;
1458 }
1459 while( len-- ) {
1460 varData->undef[len] =
1461 ( iarray[0]!=0L && iarray[0]==iarray[len+1] );
1462 }
1463 varData->data = iarray + 1;
1464 break;
1465
1466 case DOUBLE:
1467 rarray = (double*)icol->array;
1468 if (varData->undef)
1469 free( varData->undef );
1470 varData->undef = (char*)malloc( len*sizeof(char) );
1471 if( varData->undef==NULL ) {
1472 gParse.status = MEMORY_ALLOCATION;
1473 break;
1474 }
1475 while( len-- ) {
1476 varData->undef[len] =
1477 ( rarray[0]!=0.0 && rarray[0]==rarray[len+1]);
1478 }
1479 varData->data = rarray + 1;
1480 break;
1481
1482 default:
1483 snprintf(msg, 80, "SetupDataArrays, unhandled type %d\n",
1484 varData->type);
1485 ffpmsg(msg);
1486 }
1487
1488 if( gParse.status ) { /* Deallocate NULL arrays of previous columns */
1489 while( i-- ) {
1490 varData = gParse.varData + i;
1491 if( varData->type==BITSTR )
1492 FREE( ((char**)varData->data)[0] );
1493 FREE( varData->undef );
1494 varData->undef = NULL;
1495 }
1496 return;
1497 }
1498 }
1499 }
1500
1501 /*--------------------------------------------------------------------------*/
ffcvtn(int inputType,void * input,char * undef,long ntodo,int outputType,void * nulval,void * output,int * anynull,int * status)1502 int ffcvtn( int inputType, /* I - Data type of input array */
1503 void *input, /* I - Input array of type inputType */
1504 char *undef, /* I - Array of flags indicating UNDEF elems */
1505 long ntodo, /* I - Number of elements to process */
1506 int outputType, /* I - Data type of output array */
1507 void *nulval, /* I - Ptr to value to use for UNDEF elements */
1508 void *output, /* O - Output array of type outputType */
1509 int *anynull, /* O - Any nulls flagged? */
1510 int *status ) /* O - Error status */
1511 /* */
1512 /* Convert an array of any input data type to an array of any output */
1513 /* data type, using an array of UNDEF flags to assign nulvals to */
1514 /*--------------------------------------------------------------------------*/
1515 {
1516 long i;
1517
1518 switch( outputType ) {
1519
1520 case TLOGICAL:
1521 switch( inputType ) {
1522 case TLOGICAL:
1523 case TBYTE:
1524 for( i=0; i<ntodo; i++ )
1525 if( ((unsigned char*)input)[i] )
1526 ((unsigned char*)output)[i] = 1;
1527 else
1528 ((unsigned char*)output)[i] = 0;
1529 break;
1530 case TSHORT:
1531 for( i=0; i<ntodo; i++ )
1532 if( ((short*)input)[i] )
1533 ((unsigned char*)output)[i] = 1;
1534 else
1535 ((unsigned char*)output)[i] = 0;
1536 break;
1537 case TLONG:
1538 for( i=0; i<ntodo; i++ )
1539 if( ((long*)input)[i] )
1540 ((unsigned char*)output)[i] = 1;
1541 else
1542 ((unsigned char*)output)[i] = 0;
1543 break;
1544 case TFLOAT:
1545 for( i=0; i<ntodo; i++ )
1546 if( ((float*)input)[i] )
1547 ((unsigned char*)output)[i] = 1;
1548 else
1549 ((unsigned char*)output)[i] = 0;
1550 break;
1551 case TDOUBLE:
1552 for( i=0; i<ntodo; i++ )
1553 if( ((double*)input)[i] )
1554 ((unsigned char*)output)[i] = 1;
1555 else
1556 ((unsigned char*)output)[i] = 0;
1557 break;
1558 default:
1559 *status = BAD_DATATYPE;
1560 break;
1561 }
1562 for(i=0;i<ntodo;i++) {
1563 if( undef[i] ) {
1564 ((unsigned char*)output)[i] = *(unsigned char*)nulval;
1565 *anynull = 1;
1566 }
1567 }
1568 break;
1569
1570 case TBYTE:
1571 switch( inputType ) {
1572 case TLOGICAL:
1573 case TBYTE:
1574 for( i=0; i<ntodo; i++ )
1575 ((unsigned char*)output)[i] = ((unsigned char*)input)[i];
1576 break;
1577 case TSHORT:
1578 fffi2i1((short*)input,ntodo,1.,0.,0,0,0,NULL,NULL,(unsigned char*)output,status);
1579 break;
1580 case TLONG:
1581 for (i = 0; i < ntodo; i++) {
1582 if( undef[i] ) {
1583 ((unsigned char*)output)[i] = *(unsigned char*)nulval;
1584 *anynull = 1;
1585 } else {
1586 if( ((long*)input)[i] < 0 ) {
1587 *status = OVERFLOW_ERR;
1588 ((unsigned char*)output)[i] = 0;
1589 } else if( ((long*)input)[i] > UCHAR_MAX ) {
1590 *status = OVERFLOW_ERR;
1591 ((unsigned char*)output)[i] = UCHAR_MAX;
1592 } else
1593 ((unsigned char*)output)[i] =
1594 (unsigned char) ((long*)input)[i];
1595 }
1596 }
1597 return( *status );
1598 case TFLOAT:
1599 fffr4i1((float*)input,ntodo,1.,0.,0,0,NULL,NULL,
1600 (unsigned char*)output,status);
1601 break;
1602 case TDOUBLE:
1603 fffr8i1((double*)input,ntodo,1.,0.,0,0,NULL,NULL,
1604 (unsigned char*)output,status);
1605 break;
1606 default:
1607 *status = BAD_DATATYPE;
1608 break;
1609 }
1610 for(i=0;i<ntodo;i++) {
1611 if( undef[i] ) {
1612 ((unsigned char*)output)[i] = *(unsigned char*)nulval;
1613 *anynull = 1;
1614 }
1615 }
1616 break;
1617
1618 case TSHORT:
1619 switch( inputType ) {
1620 case TLOGICAL:
1621 case TBYTE:
1622 for( i=0; i<ntodo; i++ )
1623 ((short*)output)[i] = ((unsigned char*)input)[i];
1624 break;
1625 case TSHORT:
1626 for( i=0; i<ntodo; i++ )
1627 ((short*)output)[i] = ((short*)input)[i];
1628 break;
1629 case TLONG:
1630 for (i = 0; i < ntodo; i++) {
1631 if( undef[i] ) {
1632 ((short*)output)[i] = *(short*)nulval;
1633 *anynull = 1;
1634 } else {
1635 if( ((long*)input)[i] < SHRT_MIN ) {
1636 *status = OVERFLOW_ERR;
1637 ((short*)output)[i] = SHRT_MIN;
1638 } else if ( ((long*)input)[i] > SHRT_MAX ) {
1639 *status = OVERFLOW_ERR;
1640 ((short*)output)[i] = SHRT_MAX;
1641 } else
1642 ((short*)output)[i] = (short) ((long*)input)[i];
1643 }
1644 }
1645 return( *status );
1646 case TFLOAT:
1647 fffr4i2((float*)input,ntodo,1.,0.,0,0,NULL,NULL,
1648 (short*)output,status);
1649 break;
1650 case TDOUBLE:
1651 fffr8i2((double*)input,ntodo,1.,0.,0,0,NULL,NULL,
1652 (short*)output,status);
1653 break;
1654 default:
1655 *status = BAD_DATATYPE;
1656 break;
1657 }
1658 for(i=0;i<ntodo;i++) {
1659 if( undef[i] ) {
1660 ((short*)output)[i] = *(short*)nulval;
1661 *anynull = 1;
1662 }
1663 }
1664 break;
1665
1666 case TINT:
1667 switch( inputType ) {
1668 case TLOGICAL:
1669 case TBYTE:
1670 for( i=0; i<ntodo; i++ )
1671 ((int*)output)[i] = ((unsigned char*)input)[i];
1672 break;
1673 case TSHORT:
1674 for( i=0; i<ntodo; i++ )
1675 ((int*)output)[i] = ((short*)input)[i];
1676 break;
1677 case TLONG:
1678 for( i=0; i<ntodo; i++ )
1679 ((int*)output)[i] = ((long*)input)[i];
1680 break;
1681 case TFLOAT:
1682 fffr4int((float*)input,ntodo,1.,0.,0,0,NULL,NULL,
1683 (int*)output,status);
1684 break;
1685 case TDOUBLE:
1686 fffr8int((double*)input,ntodo,1.,0.,0,0,NULL,NULL,
1687 (int*)output,status);
1688 break;
1689 default:
1690 *status = BAD_DATATYPE;
1691 break;
1692 }
1693 for(i=0;i<ntodo;i++) {
1694 if( undef[i] ) {
1695 ((int*)output)[i] = *(int*)nulval;
1696 *anynull = 1;
1697 }
1698 }
1699 break;
1700
1701 case TLONG:
1702 switch( inputType ) {
1703 case TLOGICAL:
1704 case TBYTE:
1705 for( i=0; i<ntodo; i++ )
1706 ((long*)output)[i] = ((unsigned char*)input)[i];
1707 break;
1708 case TSHORT:
1709 for( i=0; i<ntodo; i++ )
1710 ((long*)output)[i] = ((short*)input)[i];
1711 break;
1712 case TLONG:
1713 for( i=0; i<ntodo; i++ )
1714 ((long*)output)[i] = ((long*)input)[i];
1715 break;
1716 case TFLOAT:
1717 fffr4i4((float*)input,ntodo,1.,0.,0,0,NULL,NULL,
1718 (long*)output,status);
1719 break;
1720 case TDOUBLE:
1721 fffr8i4((double*)input,ntodo,1.,0.,0,0,NULL,NULL,
1722 (long*)output,status);
1723 break;
1724 default:
1725 *status = BAD_DATATYPE;
1726 break;
1727 }
1728 for(i=0;i<ntodo;i++) {
1729 if( undef[i] ) {
1730 ((long*)output)[i] = *(long*)nulval;
1731 *anynull = 1;
1732 }
1733 }
1734 break;
1735
1736 case TLONGLONG:
1737 switch( inputType ) {
1738 case TLOGICAL:
1739 case TBYTE:
1740 for( i=0; i<ntodo; i++ )
1741 ((LONGLONG*)output)[i] = ((unsigned char*)input)[i];
1742 break;
1743 case TSHORT:
1744 for( i=0; i<ntodo; i++ )
1745 ((LONGLONG*)output)[i] = ((short*)input)[i];
1746 break;
1747 case TLONG:
1748 for( i=0; i<ntodo; i++ )
1749 ((LONGLONG*)output)[i] = ((long*)input)[i];
1750 break;
1751 case TFLOAT:
1752 fffr4i8((float*)input,ntodo,1.,0.,0,0,NULL,NULL,
1753 (LONGLONG*)output,status);
1754 break;
1755 case TDOUBLE:
1756 fffr8i8((double*)input,ntodo,1.,0.,0,0,NULL,NULL,
1757 (LONGLONG*)output,status);
1758
1759 break;
1760 default:
1761 *status = BAD_DATATYPE;
1762 break;
1763 }
1764 for(i=0;i<ntodo;i++) {
1765 if( undef[i] ) {
1766 ((LONGLONG*)output)[i] = *(LONGLONG*)nulval;
1767 *anynull = 1;
1768 }
1769 }
1770 break;
1771
1772 case TFLOAT:
1773 switch( inputType ) {
1774 case TLOGICAL:
1775 case TBYTE:
1776 for( i=0; i<ntodo; i++ )
1777 ((float*)output)[i] = ((unsigned char*)input)[i];
1778 break;
1779 case TSHORT:
1780 for( i=0; i<ntodo; i++ )
1781 ((float*)output)[i] = ((short*)input)[i];
1782 break;
1783 case TLONG:
1784 for( i=0; i<ntodo; i++ )
1785 ((float*)output)[i] = (float) ((long*)input)[i];
1786 break;
1787 case TFLOAT:
1788 for( i=0; i<ntodo; i++ )
1789 ((float*)output)[i] = ((float*)input)[i];
1790 break;
1791 case TDOUBLE:
1792 fffr8r4((double*)input,ntodo,1.,0.,0,0,NULL,NULL,
1793 (float*)output,status);
1794 break;
1795 default:
1796 *status = BAD_DATATYPE;
1797 break;
1798 }
1799 for(i=0;i<ntodo;i++) {
1800 if( undef[i] ) {
1801 ((float*)output)[i] = *(float*)nulval;
1802 *anynull = 1;
1803 }
1804 }
1805 break;
1806
1807 case TDOUBLE:
1808 switch( inputType ) {
1809 case TLOGICAL:
1810 case TBYTE:
1811 for( i=0; i<ntodo; i++ )
1812 ((double*)output)[i] = ((unsigned char*)input)[i];
1813 break;
1814 case TSHORT:
1815 for( i=0; i<ntodo; i++ )
1816 ((double*)output)[i] = ((short*)input)[i];
1817 break;
1818 case TLONG:
1819 for( i=0; i<ntodo; i++ )
1820 ((double*)output)[i] = ((long*)input)[i];
1821 break;
1822 case TFLOAT:
1823 for( i=0; i<ntodo; i++ )
1824 ((double*)output)[i] = ((float*)input)[i];
1825 break;
1826 case TDOUBLE:
1827 for( i=0; i<ntodo; i++ )
1828 ((double*)output)[i] = ((double*)input)[i];
1829 break;
1830 default:
1831 *status = BAD_DATATYPE;
1832 break;
1833 }
1834 for(i=0;i<ntodo;i++) {
1835 if( undef[i] ) {
1836 ((double*)output)[i] = *(double*)nulval;
1837 *anynull = 1;
1838 }
1839 }
1840 break;
1841
1842 default:
1843 *status = BAD_DATATYPE;
1844 break;
1845 }
1846
1847 return ( *status );
1848 }
1849
1850 /*---------------------------------------------------------------------------*/
fffrwc(fitsfile * fptr,char * expr,char * timeCol,char * parCol,char * valCol,long ntimes,double * times,char * time_status,int * status)1851 int fffrwc( fitsfile *fptr, /* I - Input FITS file */
1852 char *expr, /* I - Boolean expression */
1853 char *timeCol, /* I - Name of time column */
1854 char *parCol, /* I - Name of parameter column */
1855 char *valCol, /* I - Name of value column */
1856 long ntimes, /* I - Number of distinct times in file */
1857 double *times, /* O - Array of times in file */
1858 char *time_status, /* O - Array of boolean results */
1859 int *status ) /* O - Error status */
1860 /* */
1861 /* Evaluate a boolean expression for each time in a compressed file, */
1862 /* returning an array of flags indicating which times evaluated to TRUE/FALSE*/
1863 /*---------------------------------------------------------------------------*/
1864 {
1865 parseInfo Info;
1866 long alen, width;
1867 int parNo, typecode;
1868 int naxis, constant, nCol=0;
1869 long nelem, naxes[MAXDIMS], elem;
1870 char result;
1871
1872 if( *status ) return( *status );
1873
1874 fits_get_colnum( fptr, CASEINSEN, timeCol, &gParse.timeCol, status );
1875 fits_get_colnum( fptr, CASEINSEN, parCol, &gParse.parCol , status );
1876 fits_get_colnum( fptr, CASEINSEN, valCol, &gParse.valCol, status );
1877 if( *status ) return( *status );
1878
1879 if( ffiprs( fptr, 1, expr, MAXDIMS, &Info.datatype, &nelem,
1880 &naxis, naxes, status ) ) {
1881 ffcprs();
1882 return( *status );
1883 }
1884 if( nelem<0 ) {
1885 constant = 1;
1886 nelem = -nelem;
1887 nCol = gParse.nCols;
1888 gParse.nCols = 0; /* Ignore all column references */
1889 } else
1890 constant = 0;
1891
1892 if( Info.datatype!=TLOGICAL || nelem!=1 ) {
1893 ffcprs();
1894 ffpmsg("Expression does not evaluate to a logical scalar.");
1895 return( *status = PARSE_BAD_TYPE );
1896 }
1897
1898 /*******************************************/
1899 /* Allocate data arrays for each parameter */
1900 /*******************************************/
1901
1902 parNo = gParse.nCols;
1903 while( parNo-- ) {
1904 switch( gParse.colData[parNo].datatype ) {
1905 case TLONG:
1906 if( (gParse.colData[parNo].array =
1907 (long *)malloc( (ntimes+1)*sizeof(long) )) )
1908 ((long*)gParse.colData[parNo].array)[0] = 1234554321;
1909 else
1910 *status = MEMORY_ALLOCATION;
1911 break;
1912 case TDOUBLE:
1913 if( (gParse.colData[parNo].array =
1914 (double *)malloc( (ntimes+1)*sizeof(double) )) )
1915 ((double*)gParse.colData[parNo].array)[0] = DOUBLENULLVALUE;
1916 else
1917 *status = MEMORY_ALLOCATION;
1918 break;
1919 case TSTRING:
1920 if( !fits_get_coltype( fptr, gParse.valCol, &typecode,
1921 &alen, &width, status ) ) {
1922 alen++;
1923 if( (gParse.colData[parNo].array =
1924 (char **)malloc( (ntimes+1)*sizeof(char*) )) ) {
1925 if( (((char **)gParse.colData[parNo].array)[0] =
1926 (char *)malloc( (ntimes+1)*sizeof(char)*alen )) ) {
1927 for( elem=1; elem<=ntimes; elem++ )
1928 ((char **)gParse.colData[parNo].array)[elem] =
1929 ((char **)gParse.colData[parNo].array)[elem-1]+alen;
1930 ((char **)gParse.colData[parNo].array)[0][0] = '\0';
1931 } else {
1932 free( gParse.colData[parNo].array );
1933 *status = MEMORY_ALLOCATION;
1934 }
1935 } else {
1936 *status = MEMORY_ALLOCATION;
1937 }
1938 }
1939 break;
1940 }
1941 if( *status ) {
1942 while( parNo-- ) {
1943 if( gParse.colData[parNo].datatype==TSTRING )
1944 FREE( ((char **)gParse.colData[parNo].array)[0] );
1945 FREE( gParse.colData[parNo].array );
1946 }
1947 return( *status );
1948 }
1949 }
1950
1951 /**********************************************************************/
1952 /* Read data from columns needed for the expression and then parse it */
1953 /**********************************************************************/
1954
1955 if( !uncompress_hkdata( fptr, ntimes, times, status ) ) {
1956 if( constant ) {
1957 result = gParse.Nodes[gParse.resultNode].value.data.log;
1958 elem = ntimes;
1959 while( elem-- ) time_status[elem] = result;
1960 } else {
1961 Info.dataPtr = time_status;
1962 Info.nullPtr = NULL;
1963 Info.maxRows = ntimes;
1964 *status = parse_data( ntimes, 0, 1, ntimes, gParse.nCols,
1965 gParse.colData, (void*)&Info );
1966 }
1967 }
1968
1969 /************/
1970 /* Clean up */
1971 /************/
1972
1973 parNo = gParse.nCols;
1974 while ( parNo-- ) {
1975 if( gParse.colData[parNo].datatype==TSTRING )
1976 FREE( ((char **)gParse.colData[parNo].array)[0] );
1977 FREE( gParse.colData[parNo].array );
1978 }
1979
1980 if( constant ) gParse.nCols = nCol;
1981
1982 ffcprs();
1983 return(*status);
1984 }
1985
1986 /*---------------------------------------------------------------------------*/
uncompress_hkdata(fitsfile * fptr,long ntimes,double * times,int * status)1987 int uncompress_hkdata( fitsfile *fptr,
1988 long ntimes,
1989 double *times,
1990 int *status )
1991 /* */
1992 /* description */
1993 /*---------------------------------------------------------------------------*/
1994 {
1995 char parName[256], *sPtr[1], found[1000];
1996 int parNo, anynul;
1997 long naxis2, row, currelem;
1998 double currtime, newtime;
1999
2000 sPtr[0] = parName;
2001 currelem = 0;
2002 currtime = -1e38;
2003
2004 parNo=gParse.nCols;
2005 while( parNo-- ) found[parNo] = 0;
2006
2007 if( ffgkyj( fptr, "NAXIS2", &naxis2, NULL, status ) ) return( *status );
2008
2009 for( row=1; row<=naxis2; row++ ) {
2010 if( ffgcvd( fptr, gParse.timeCol, row, 1L, 1L, 0.0,
2011 &newtime, &anynul, status ) ) return( *status );
2012 if( newtime != currtime ) {
2013 /* New time encountered... propogate parameters to next row */
2014 if( currelem==ntimes ) {
2015 ffpmsg("Found more unique time stamps than caller indicated");
2016 return( *status = PARSE_BAD_COL );
2017 }
2018 times[currelem++] = currtime = newtime;
2019 parNo = gParse.nCols;
2020 while( parNo-- ) {
2021 switch( gParse.colData[parNo].datatype ) {
2022 case TLONG:
2023 ((long*)gParse.colData[parNo].array)[currelem] =
2024 ((long*)gParse.colData[parNo].array)[currelem-1];
2025 break;
2026 case TDOUBLE:
2027 ((double*)gParse.colData[parNo].array)[currelem] =
2028 ((double*)gParse.colData[parNo].array)[currelem-1];
2029 break;
2030 case TSTRING:
2031 strcpy( ((char **)gParse.colData[parNo].array)[currelem],
2032 ((char **)gParse.colData[parNo].array)[currelem-1] );
2033 break;
2034 }
2035 }
2036 }
2037
2038 if( ffgcvs( fptr, gParse.parCol, row, 1L, 1L, "",
2039 sPtr, &anynul, status ) ) return( *status );
2040 parNo = gParse.nCols;
2041 while( parNo-- )
2042 if( !fits_strcasecmp( parName, gParse.varData[parNo].name ) ) break;
2043
2044 if( parNo>=0 ) {
2045 found[parNo] = 1; /* Flag this parameter as found */
2046 switch( gParse.colData[parNo].datatype ) {
2047 case TLONG:
2048 ffgcvj( fptr, gParse.valCol, row, 1L, 1L,
2049 ((long*)gParse.colData[parNo].array)[0],
2050 ((long*)gParse.colData[parNo].array)+currelem,
2051 &anynul, status );
2052 break;
2053 case TDOUBLE:
2054 ffgcvd( fptr, gParse.valCol, row, 1L, 1L,
2055 ((double*)gParse.colData[parNo].array)[0],
2056 ((double*)gParse.colData[parNo].array)+currelem,
2057 &anynul, status );
2058 break;
2059 case TSTRING:
2060 ffgcvs( fptr, gParse.valCol, row, 1L, 1L,
2061 ((char**)gParse.colData[parNo].array)[0],
2062 ((char**)gParse.colData[parNo].array)+currelem,
2063 &anynul, status );
2064 break;
2065 }
2066 if( *status ) return( *status );
2067 }
2068 }
2069
2070 if( currelem<ntimes ) {
2071 ffpmsg("Found fewer unique time stamps than caller indicated");
2072 return( *status = PARSE_BAD_COL );
2073 }
2074
2075 /* Check for any parameters which were not located in the table */
2076 parNo = gParse.nCols;
2077 while( parNo-- )
2078 if( !found[parNo] ) {
2079 snprintf( parName, 256, "Parameter not found: %-30s",
2080 gParse.varData[parNo].name );
2081 ffpmsg( parName );
2082 *status = PARSE_SYNTAX_ERR;
2083 }
2084 return( *status );
2085 }
2086
2087 /*---------------------------------------------------------------------------*/
ffffrw(fitsfile * fptr,char * expr,long * rownum,int * status)2088 int ffffrw( fitsfile *fptr, /* I - Input FITS file */
2089 char *expr, /* I - Boolean expression */
2090 long *rownum, /* O - First row of table to eval to T */
2091 int *status ) /* O - Error status */
2092 /* */
2093 /* Evaluate a boolean expression, returning the row number of the first */
2094 /* row which evaluates to TRUE */
2095 /*---------------------------------------------------------------------------*/
2096 {
2097 int naxis, constant, dtype;
2098 long nelem, naxes[MAXDIMS];
2099 char result;
2100
2101 if( *status ) return( *status );
2102
2103 FFLOCK;
2104 if( ffiprs( fptr, 0, expr, MAXDIMS, &dtype, &nelem, &naxis,
2105 naxes, status ) ) {
2106 ffcprs();
2107 FFUNLOCK;
2108 return( *status );
2109 }
2110 if( nelem<0 ) {
2111 constant = 1;
2112 nelem = -nelem;
2113 } else
2114 constant = 0;
2115
2116 if( dtype!=TLOGICAL || nelem!=1 ) {
2117 ffcprs();
2118 ffpmsg("Expression does not evaluate to a logical scalar.");
2119 FFUNLOCK;
2120 return( *status = PARSE_BAD_TYPE );
2121 }
2122
2123 *rownum = 0;
2124 if( constant ) { /* No need to call parser... have result from ffiprs */
2125 result = gParse.Nodes[gParse.resultNode].value.data.log;
2126 if( result ) {
2127 /* Make sure there is at least 1 row in table */
2128 ffgnrw( fptr, &nelem, status );
2129 if( nelem )
2130 *rownum = 1;
2131 }
2132 } else {
2133 if( ffiter( gParse.nCols, gParse.colData, 0, 0,
2134 ffffrw_work, (void*)rownum, status ) == -1 )
2135 *status = 0; /* -1 indicates exitted without error before end... OK */
2136 }
2137
2138 ffcprs();
2139 FFUNLOCK;
2140 return(*status);
2141 }
2142
2143 /*---------------------------------------------------------------------------*/
ffffrw_work(long totalrows,long offset,long firstrow,long nrows,int nCols,iteratorCol * colData,void * userPtr)2144 int ffffrw_work(long totalrows, /* I - Total rows to be processed */
2145 long offset, /* I - Number of rows skipped at start*/
2146 long firstrow, /* I - First row of this iteration */
2147 long nrows, /* I - Number of rows in this iter */
2148 int nCols, /* I - Number of columns in use */
2149 iteratorCol *colData, /* IO- Column information/data */
2150 void *userPtr ) /* I - Data handling instructions */
2151 /* */
2152 /* Iterator work function which calls the parser and searches for the */
2153 /* first row which evaluates to TRUE. */
2154 /*---------------------------------------------------------------------------*/
2155 {
2156 long idx;
2157 Node *result;
2158
2159 Evaluate_Parser( firstrow, nrows );
2160
2161 if( !gParse.status ) {
2162
2163 result = gParse.Nodes + gParse.resultNode;
2164 if( result->operation==CONST_OP ) {
2165
2166 if( result->value.data.log ) {
2167 *(long*)userPtr = firstrow;
2168 return( -1 );
2169 }
2170
2171 } else {
2172
2173 for( idx=0; idx<nrows; idx++ )
2174 if( result->value.data.logptr[idx] && !result->value.undef[idx] ) {
2175 *(long*)userPtr = firstrow + idx;
2176 return( -1 );
2177 }
2178 }
2179 }
2180
2181 return( gParse.status );
2182 }
2183
2184
set_image_col_types(fitsfile * fptr,const char * name,int bitpix,DataInfo * varInfo,iteratorCol * colIter)2185 static int set_image_col_types (fitsfile * fptr, const char * name, int bitpix,
2186 DataInfo * varInfo, iteratorCol *colIter) {
2187
2188 int istatus;
2189 double tscale, tzero;
2190 char temp[80];
2191
2192 switch (bitpix) {
2193 case BYTE_IMG:
2194 case SHORT_IMG:
2195 case LONG_IMG:
2196 istatus = 0;
2197 if (fits_read_key(fptr, TDOUBLE, "BZERO", &tzero, NULL, &istatus))
2198 tzero = 0.0;
2199
2200 istatus = 0;
2201 if (fits_read_key(fptr, TDOUBLE, "BSCALE", &tscale, NULL, &istatus))
2202 tscale = 1.0;
2203
2204 if (tscale == 1.0 && (tzero == 0.0 || tzero == 32768.0 )) {
2205 varInfo->type = LONG;
2206 colIter->datatype = TLONG;
2207 }
2208 else {
2209 varInfo->type = DOUBLE;
2210 colIter->datatype = TDOUBLE;
2211 if (DEBUG_PIXFILTER)
2212 printf("use DOUBLE for %s with BSCALE=%g/BZERO=%g\n",
2213 name, tscale, tzero);
2214 }
2215 break;
2216
2217 case LONGLONG_IMG:
2218 case FLOAT_IMG:
2219 case DOUBLE_IMG:
2220 varInfo->type = DOUBLE;
2221 colIter->datatype = TDOUBLE;
2222 break;
2223 default:
2224 snprintf(temp, 80,"set_image_col_types: unrecognized image bitpix [%d]\n",
2225 bitpix);
2226 ffpmsg(temp);
2227 return gParse.status = PARSE_BAD_TYPE;
2228 }
2229 return 0;
2230 }
2231
2232
2233 /*************************************************************************
2234
2235 Functions used by the evaluator to access FITS data
2236 (find_column, find_keywd, allocateCol, load_column)
2237
2238 *************************************************************************/
2239
find_column(char * colName,void * itslval)2240 static int find_column( char *colName, void *itslval )
2241 {
2242 FFSTYPE *thelval = (FFSTYPE*)itslval;
2243 int col_cnt, status;
2244 int colnum, typecode, type;
2245 long repeat, width;
2246 fitsfile *fptr;
2247 char temp[80];
2248 double tzero,tscale;
2249 int istatus;
2250 DataInfo *varInfo;
2251 iteratorCol *colIter;
2252
2253 if (DEBUG_PIXFILTER)
2254 printf("find_column(%s)\n", colName);
2255
2256 if( *colName == '#' )
2257 return( find_keywd( colName + 1, itslval ) );
2258
2259 fptr = gParse.def_fptr;
2260
2261 status = 0;
2262 col_cnt = gParse.nCols;
2263
2264 if (gParse.hdutype == IMAGE_HDU) {
2265 int i;
2266 if (!gParse.pixFilter) {
2267 gParse.status = COL_NOT_FOUND;
2268 ffpmsg("find_column: IMAGE_HDU but no PixelFilter");
2269 return pERROR;
2270 }
2271
2272 colnum = -1;
2273 for (i = 0; i < gParse.pixFilter->count; ++i) {
2274 if (!fits_strcasecmp(colName, gParse.pixFilter->tag[i]))
2275 colnum = i;
2276 }
2277 if (colnum < 0) {
2278 snprintf(temp, 80, "find_column: PixelFilter tag %s not found", colName);
2279 ffpmsg(temp);
2280 gParse.status = COL_NOT_FOUND;
2281 return pERROR;
2282 }
2283
2284 if( allocateCol( col_cnt, &gParse.status ) ) return pERROR;
2285
2286 varInfo = gParse.varData + col_cnt;
2287 colIter = gParse.colData + col_cnt;
2288
2289 fptr = gParse.pixFilter->ifptr[colnum];
2290 fits_get_img_param(fptr,
2291 MAXDIMS,
2292 &typecode, /* actually bitpix */
2293 &varInfo->naxis,
2294 &varInfo->naxes[0],
2295 &status);
2296 varInfo->nelem = 1;
2297 type = COLUMN;
2298 if (set_image_col_types(fptr, colName, typecode, varInfo, colIter))
2299 return pERROR;
2300 colIter->fptr = fptr;
2301 colIter->iotype = InputCol;
2302 }
2303 else { /* HDU holds a table */
2304 if( gParse.compressed )
2305 colnum = gParse.valCol;
2306 else
2307 if( fits_get_colnum( fptr, CASEINSEN, colName, &colnum, &status ) ) {
2308 if( status == COL_NOT_FOUND ) {
2309 type = find_keywd( colName, itslval );
2310 if( type != pERROR ) ffcmsg();
2311 return( type );
2312 }
2313 gParse.status = status;
2314 return pERROR;
2315 }
2316
2317 if( fits_get_coltype( fptr, colnum, &typecode,
2318 &repeat, &width, &status ) ) {
2319 gParse.status = status;
2320 return pERROR;
2321 }
2322
2323 if( allocateCol( col_cnt, &gParse.status ) ) return pERROR;
2324
2325 varInfo = gParse.varData + col_cnt;
2326 colIter = gParse.colData + col_cnt;
2327
2328 fits_iter_set_by_num( colIter, fptr, colnum, 0, InputCol );
2329 }
2330
2331 /* Make sure we don't overflow variable name array */
2332 strncpy(varInfo->name,colName,MAXVARNAME);
2333 varInfo->name[MAXVARNAME] = '\0';
2334
2335 if (gParse.hdutype != IMAGE_HDU) {
2336 switch( typecode ) {
2337 case TBIT:
2338 varInfo->type = BITSTR;
2339 colIter->datatype = TBYTE;
2340 type = BITCOL;
2341 break;
2342 case TBYTE:
2343 case TSHORT:
2344 case TLONG:
2345 /* The datatype of column with TZERO and TSCALE keywords might be
2346 float or double.
2347 */
2348 snprintf(temp,80,"TZERO%d",colnum);
2349 istatus = 0;
2350 if(fits_read_key(fptr,TDOUBLE,temp,&tzero,NULL,&istatus)) {
2351 tzero = 0.0;
2352 }
2353 snprintf(temp,80,"TSCAL%d",colnum);
2354 istatus = 0;
2355 if(fits_read_key(fptr,TDOUBLE,temp,&tscale,NULL,&istatus)) {
2356 tscale = 1.0;
2357 }
2358 if (tscale == 1.0 && (tzero == 0.0 || tzero == 32768.0 )) {
2359 varInfo->type = LONG;
2360 colIter->datatype = TLONG;
2361 /* Reading an unsigned long column as a long can cause overflow errors.
2362 Treat the column as a double instead.
2363 } else if (tscale == 1.0 && tzero == 2147483648.0 ) {
2364 varInfo->type = LONG;
2365 colIter->datatype = TULONG;
2366 */
2367
2368 }
2369 else {
2370 varInfo->type = DOUBLE;
2371 colIter->datatype = TDOUBLE;
2372 }
2373 type = COLUMN;
2374 break;
2375 /*
2376 For now, treat 8-byte integer columns as type double.
2377 This can lose precision, so the better long term solution
2378 will be to add support for TLONGLONG as a separate datatype.
2379 */
2380 case TLONGLONG:
2381 case TFLOAT:
2382 case TDOUBLE:
2383 varInfo->type = DOUBLE;
2384 colIter->datatype = TDOUBLE;
2385 type = COLUMN;
2386 break;
2387 case TLOGICAL:
2388 varInfo->type = BOOLEAN;
2389 colIter->datatype = TLOGICAL;
2390 type = BCOLUMN;
2391 break;
2392 case TSTRING:
2393 varInfo->type = STRING;
2394 colIter->datatype = TSTRING;
2395 type = SCOLUMN;
2396 if ( width >= MAX_STRLEN ) {
2397 snprintf(temp, 80, "column %d is wider than maximum %d characters",
2398 colnum, MAX_STRLEN-1);
2399 ffpmsg(temp);
2400 gParse.status = PARSE_LRG_VECTOR;
2401 return pERROR;
2402 }
2403 if( gParse.hdutype == ASCII_TBL ) repeat = width;
2404 break;
2405 default:
2406 if (typecode < 0) {
2407 snprintf(temp, 80,"variable-length array columns are not supported. typecode = %d", typecode);
2408 ffpmsg(temp);
2409 }
2410 gParse.status = PARSE_BAD_TYPE;
2411 return pERROR;
2412 }
2413 varInfo->nelem = repeat;
2414 if( repeat>1 && typecode!=TSTRING ) {
2415 if( fits_read_tdim( fptr, colnum, MAXDIMS,
2416 &varInfo->naxis,
2417 &varInfo->naxes[0], &status )
2418 ) {
2419 gParse.status = status;
2420 return pERROR;
2421 }
2422 } else {
2423 varInfo->naxis = 1;
2424 varInfo->naxes[0] = 1;
2425 }
2426 }
2427 gParse.nCols++;
2428 thelval->lng = col_cnt;
2429
2430 return( type );
2431 }
2432
find_keywd(char * keyname,void * itslval)2433 static int find_keywd(char *keyname, void *itslval )
2434 {
2435 FFSTYPE *thelval = (FFSTYPE*)itslval;
2436 int status, type;
2437 char keyvalue[FLEN_VALUE], dtype;
2438 fitsfile *fptr;
2439 double rval;
2440 int bval;
2441 long ival;
2442
2443 status = 0;
2444 fptr = gParse.def_fptr;
2445 if( fits_read_keyword( fptr, keyname, keyvalue, NULL, &status ) ) {
2446 if( status == KEY_NO_EXIST ) {
2447 /* Do this since ffgkey doesn't put an error message on stack */
2448 snprintf(keyvalue,FLEN_VALUE, "ffgkey could not find keyword: %s",keyname);
2449 ffpmsg(keyvalue);
2450 }
2451 gParse.status = status;
2452 return( pERROR );
2453 }
2454
2455 if( fits_get_keytype( keyvalue, &dtype, &status ) ) {
2456 gParse.status = status;
2457 return( pERROR );
2458 }
2459
2460 switch( dtype ) {
2461 case 'C':
2462 fits_read_key_str( fptr, keyname, keyvalue, NULL, &status );
2463 type = STRING;
2464 strcpy( thelval->str , keyvalue );
2465 break;
2466 case 'L':
2467 fits_read_key_log( fptr, keyname, &bval, NULL, &status );
2468 type = BOOLEAN;
2469 thelval->log = bval;
2470 break;
2471 case 'I':
2472 fits_read_key_lng( fptr, keyname, &ival, NULL, &status );
2473 type = LONG;
2474 thelval->lng = ival;
2475 break;
2476 case 'F':
2477 fits_read_key_dbl( fptr, keyname, &rval, NULL, &status );
2478 type = DOUBLE;
2479 thelval->dbl = rval;
2480 break;
2481 default:
2482 type = pERROR;
2483 break;
2484 }
2485
2486 if( status ) {
2487 gParse.status=status;
2488 return pERROR;
2489 }
2490
2491 return( type );
2492 }
2493
allocateCol(int nCol,int * status)2494 static int allocateCol( int nCol, int *status )
2495 {
2496 if( (nCol%25)==0 ) {
2497 if( nCol ) {
2498 gParse.colData = (iteratorCol*) realloc( gParse.colData,
2499 (nCol+25)*sizeof(iteratorCol) );
2500 gParse.varData = (DataInfo *) realloc( gParse.varData,
2501 (nCol+25)*sizeof(DataInfo) );
2502 } else {
2503 gParse.colData = (iteratorCol*) malloc( 25*sizeof(iteratorCol) );
2504 gParse.varData = (DataInfo *) malloc( 25*sizeof(DataInfo) );
2505 }
2506 if( gParse.colData == NULL
2507 || gParse.varData == NULL ) {
2508 if( gParse.colData ) free(gParse.colData);
2509 if( gParse.varData ) free(gParse.varData);
2510 gParse.colData = NULL;
2511 gParse.varData = NULL;
2512 return( *status = MEMORY_ALLOCATION );
2513 }
2514 }
2515 gParse.varData[nCol].data = NULL;
2516 gParse.varData[nCol].undef = NULL;
2517 return 0;
2518 }
2519
load_column(int varNum,long fRow,long nRows,void * data,char * undef)2520 static int load_column( int varNum, long fRow, long nRows,
2521 void *data, char *undef )
2522 {
2523 iteratorCol *var = gParse.colData+varNum;
2524 long nelem,nbytes,row,len,idx;
2525 char **bitStrs, msg[80];
2526 unsigned char *bytes;
2527 int status = 0, anynul;
2528
2529 if (gParse.hdutype == IMAGE_HDU) {
2530 /* This test would need to be on a per varNum basis to support
2531 * cross HDU operations */
2532 fits_read_imgnull(var->fptr, var->datatype, fRow, nRows,
2533 data, undef, &anynul, &status);
2534 if (DEBUG_PIXFILTER)
2535 printf("load_column: IMAGE_HDU fRow=%ld, nRows=%ld => %d\n",
2536 fRow, nRows, status);
2537 } else {
2538
2539 nelem = nRows * var->repeat;
2540
2541 switch( var->datatype ) {
2542 case TBYTE:
2543 nbytes = ((var->repeat+7)/8) * nRows;
2544 bytes = (unsigned char *)malloc( nbytes * sizeof(char) );
2545
2546 ffgcvb(var->fptr, var->colnum, fRow, 1L, nbytes,
2547 0, bytes, &anynul, &status);
2548
2549 nelem = var->repeat;
2550 bitStrs = (char **)data;
2551 for( row=0; row<nRows; row++ ) {
2552 idx = (row)*( (nelem+7)/8 ) + 1;
2553 for(len=0; len<nelem; len++) {
2554 if( bytes[idx] & (1<<(7-len%8)) )
2555 bitStrs[row][len] = '1';
2556 else
2557 bitStrs[row][len] = '0';
2558 if( len%8==7 ) idx++;
2559 }
2560 bitStrs[row][len] = '\0';
2561 }
2562
2563 FREE( (char *)bytes );
2564 break;
2565 case TSTRING:
2566 ffgcfs(var->fptr, var->colnum, fRow, 1L, nRows,
2567 (char **)data, undef, &anynul, &status);
2568 break;
2569 case TLOGICAL:
2570 ffgcfl(var->fptr, var->colnum, fRow, 1L, nelem,
2571 (char *)data, undef, &anynul, &status);
2572 break;
2573 case TLONG:
2574 ffgcfj(var->fptr, var->colnum, fRow, 1L, nelem,
2575 (long *)data, undef, &anynul, &status);
2576 break;
2577 case TDOUBLE:
2578 ffgcfd(var->fptr, var->colnum, fRow, 1L, nelem,
2579 (double *)data, undef, &anynul, &status);
2580 break;
2581 default:
2582 snprintf(msg,80,"load_column: unexpected datatype %d", var->datatype);
2583 ffpmsg(msg);
2584 }
2585 }
2586 if( status ) {
2587 gParse.status = status;
2588 return pERROR;
2589 }
2590
2591 return 0;
2592 }
2593
2594
2595 /*--------------------------------------------------------------------------*/
fits_pixel_filter(PixelFilter * filter,int * status)2596 int fits_pixel_filter (PixelFilter * filter, int * status)
2597 /* Evaluate an expression using the data in the input FITS file(s) */
2598 /*--------------------------------------------------------------------------*/
2599 {
2600 parseInfo Info = { 0 };
2601 int naxis, bitpix;
2602 long nelem, naxes[MAXDIMS];
2603 int col_cnt;
2604 Node *result;
2605 int datatype;
2606 fitsfile * infptr;
2607 fitsfile * outfptr;
2608 char * DEFAULT_TAGS[] = { "X" };
2609 char msg[256];
2610 int writeBlankKwd = 0; /* write BLANK if any output nulls? */
2611
2612 DEBUG_PIXFILTER = getenv("DEBUG_PIXFILTER") ? 1 : 0;
2613
2614 if (*status)
2615 return (*status);
2616
2617 FFLOCK;
2618 if (!filter->tag || !filter->tag[0] || !filter->tag[0][0]) {
2619 filter->tag = DEFAULT_TAGS;
2620 if (DEBUG_PIXFILTER)
2621 printf("using default tag '%s'\n", filter->tag[0]);
2622 }
2623
2624 infptr = filter->ifptr[0];
2625 outfptr = filter->ofptr;
2626 gParse.pixFilter = filter;
2627
2628 if (ffiprs(infptr, 0, filter->expression, MAXDIMS,
2629 &Info.datatype, &nelem, &naxis, naxes, status)) {
2630 goto CLEANUP;
2631 }
2632
2633 if (nelem < 0) {
2634 nelem = -nelem;
2635 }
2636
2637 {
2638 /* validate result type */
2639 const char * type = 0;
2640 switch (Info.datatype) {
2641 case TLOGICAL: type = "LOGICAL"; break;
2642 case TLONG: type = "LONG"; break;
2643 case TDOUBLE: type = "DOUBLE"; break;
2644 case TSTRING: type = "STRING";
2645 *status = pERROR;
2646 ffpmsg("pixel_filter: cannot have string image");
2647 case TBIT: type = "BIT";
2648 if (DEBUG_PIXFILTER)
2649 printf("hmm, image from bits?\n");
2650 break;
2651 default: type = "UNKNOWN?!";
2652 *status = pERROR;
2653 ffpmsg("pixel_filter: unexpected result datatype");
2654 }
2655 if (DEBUG_PIXFILTER)
2656 printf("result type is %s [%d]\n", type, Info.datatype);
2657 if (*status)
2658 goto CLEANUP;
2659 }
2660
2661 if (fits_get_img_param(infptr, MAXDIMS,
2662 &bitpix, &naxis, &naxes[0], status)) {
2663 ffpmsg("pixel_filter: unable to read input image parameters");
2664 goto CLEANUP;
2665 }
2666
2667 if (DEBUG_PIXFILTER)
2668 printf("input bitpix %d\n", bitpix);
2669
2670 if (Info.datatype == TDOUBLE) {
2671 /* for floating point expressions, set the default output image to
2672 bitpix = -32 (float) unless the default is already a double */
2673 if (bitpix != DOUBLE_IMG)
2674 bitpix = FLOAT_IMG;
2675 }
2676
2677 /* override output image bitpix if specified by caller */
2678 if (filter->bitpix)
2679 bitpix = filter->bitpix;
2680 if (DEBUG_PIXFILTER)
2681 printf("output bitpix %d\n", bitpix);
2682
2683 if (fits_create_img(outfptr, bitpix, naxis, naxes, status)) {
2684 ffpmsg("pixel_filter: unable to create output image");
2685 goto CLEANUP;
2686 }
2687
2688 /* transfer keycards */
2689 {
2690 int i, ncards, more;
2691 if (fits_get_hdrspace(infptr, &ncards, &more, status)) {
2692 ffpmsg("pixel_filter: unable to determine number of keycards");
2693 goto CLEANUP;
2694 }
2695
2696 for (i = 1; i <= ncards; ++i) {
2697
2698 int keyclass;
2699 char card[FLEN_CARD];
2700
2701 if (fits_read_record(infptr, i, card, status)) {
2702 snprintf(msg, 256,"pixel_filter: unable to read keycard %d", i);
2703 ffpmsg(msg);
2704 goto CLEANUP;
2705 }
2706
2707 keyclass = fits_get_keyclass(card);
2708 if (keyclass == TYP_STRUC_KEY) {
2709 /* output structure defined by fits_create_img */
2710 }
2711 else if (keyclass == TYP_COMM_KEY && i < 12) {
2712 /* assume this is one of the FITS standard comments */
2713 }
2714 else if (keyclass == TYP_NULL_KEY && bitpix < 0) {
2715 /* do not transfer BLANK to real output image */
2716 }
2717 else if (keyclass == TYP_SCAL_KEY && bitpix < 0) {
2718 /* do not transfer BZERO, BSCALE to real output image */
2719 }
2720 else if (fits_write_record(outfptr, card, status)) {
2721 snprintf(msg,256, "pixel_filter: unable to write keycard '%s' [%d]\n",
2722 card, *status);
2723 ffpmsg(msg);
2724 goto CLEANUP;
2725 }
2726 }
2727 }
2728
2729 switch (bitpix) {
2730 case BYTE_IMG: datatype = TLONG; Info.datatype = TBYTE; break;
2731 case SHORT_IMG: datatype = TLONG; Info.datatype = TSHORT; break;
2732 case LONG_IMG: datatype = TLONG; Info.datatype = TLONG; break;
2733 case FLOAT_IMG: datatype = TDOUBLE; Info.datatype = TFLOAT; break;
2734 case DOUBLE_IMG: datatype = TDOUBLE; Info.datatype = TDOUBLE; break;
2735
2736 default:
2737 snprintf(msg, 256,"pixel_filter: unexpected output bitpix %d\n", bitpix);
2738 ffpmsg(msg);
2739 *status = pERROR;
2740 goto CLEANUP;
2741 }
2742
2743 if (bitpix > 0) { /* arrange for NULLs in output */
2744 long nullVal = filter->blank;
2745 if (!filter->blank) {
2746 int tstatus = 0;
2747 if (fits_read_key_lng(infptr, "BLANK", &nullVal, 0, &tstatus)) {
2748
2749 writeBlankKwd = 1;
2750
2751 if (bitpix == BYTE_IMG)
2752 nullVal = UCHAR_MAX;
2753 else if (bitpix == SHORT_IMG)
2754 nullVal = SHRT_MIN;
2755 else if (bitpix == LONG_IMG) {
2756 if (sizeof(long) == 8 && sizeof(int) == 4)
2757 nullVal = INT_MIN;
2758 else
2759 nullVal = LONG_MIN;
2760 }
2761 else
2762 printf("unhandled positive output BITPIX %d\n", bitpix);
2763 }
2764
2765 filter->blank = nullVal;
2766 }
2767
2768 fits_set_imgnull(outfptr, filter->blank, status);
2769 if (DEBUG_PIXFILTER)
2770 printf("using blank %ld\n", nullVal);
2771
2772 }
2773
2774 if (!filter->keyword[0]) {
2775 iteratorCol * colIter;
2776 DataInfo * varInfo;
2777
2778 /*************************************/
2779 /* Create new iterator Output Column */
2780 /*************************************/
2781 col_cnt = gParse.nCols;
2782 if (allocateCol(col_cnt, status))
2783 goto CLEANUP;
2784 gParse.nCols++;
2785
2786 colIter = &gParse.colData[col_cnt];
2787 colIter->fptr = filter->ofptr;
2788 colIter->iotype = OutputCol;
2789 varInfo = &gParse.varData[col_cnt];
2790 set_image_col_types(colIter->fptr, "CREATED", bitpix, varInfo, colIter);
2791
2792 Info.maxRows = -1;
2793
2794 if (ffiter(gParse.nCols, gParse.colData, 0,
2795 0, parse_data, &Info, status) == -1)
2796 *status = 0;
2797 else if (*status)
2798 goto CLEANUP;
2799
2800 if (Info.anyNull) {
2801 if (writeBlankKwd) {
2802 fits_update_key_lng(outfptr, "BLANK", filter->blank, "NULL pixel value", status);
2803 if (*status)
2804 ffpmsg("pixel_filter: unable to write BLANK keyword");
2805 if (DEBUG_PIXFILTER) {
2806 printf("output has NULLs\n");
2807 printf("wrote blank [%d]\n", *status);
2808 }
2809 }
2810 }
2811 else if (bitpix > 0) /* never used a null */
2812 if (fits_set_imgnull(outfptr, -1234554321, status))
2813 ffpmsg("pixel_filter: unable to reset imgnull");
2814 }
2815 else {
2816
2817 /* Put constant result into keyword */
2818 char * parName = filter->keyword;
2819 char * parInfo = filter->comment;
2820
2821 result = gParse.Nodes + gParse.resultNode;
2822 switch (Info.datatype) {
2823 case TDOUBLE:
2824 ffukyd(outfptr, parName, result->value.data.dbl, 15, parInfo, status);
2825 break;
2826 case TLONG:
2827 ffukyj(outfptr, parName, result->value.data.lng, parInfo, status);
2828 break;
2829 case TLOGICAL:
2830 ffukyl(outfptr, parName, result->value.data.log, parInfo, status);
2831 break;
2832 case TBIT:
2833 case TSTRING:
2834 ffukys(outfptr, parName, result->value.data.str, parInfo, status);
2835 break;
2836 default:
2837 snprintf(msg, 256,"pixel_filter: unexpected constant result type [%d]\n",
2838 Info.datatype);
2839 ffpmsg(msg);
2840 }
2841 }
2842
2843 CLEANUP:
2844 ffcprs();
2845 FFUNLOCK;
2846 return (*status);
2847 }
2848