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 nullVal = LONG_MIN;
645 else if( typecode==TLONGLONG )
646 nullVal = LONGLONG_MIN;
647
648 if( nullVal ) {
649 ffpkyj( outfptr, nullKwd, nullVal, "Null value", status );
650 fits_set_btblnull( outfptr, colNo, nullVal, status );
651 newNullKwd = 1;
652 }
653 } else if( gParse.hdutype==ASCII_TBL ) {
654 ffpkys( outfptr, nullKwd, "NULL", "Null value string", status );
655 fits_set_atblnull( outfptr, colNo, "NULL", status );
656 newNullKwd = 1;
657 }
658 }
659
660 }
661
662 } else if( *status ) {
663 ffcprs();
664 FFUNLOCK;
665 return( *status );
666 } else {
667
668 /********************************************************/
669 /* Check if a TDIM keyword should be written/updated. */
670 /********************************************************/
671
672 ffkeyn("TDIM", colNo, tdimKwd, status);
673 ffgcrd( outfptr, tdimKwd, card, status );
674 if( *status==0 ) {
675 /* TDIM exists, so update it with result's dimension */
676 ffptdm( outfptr, colNo, naxis, naxes, status );
677 } else if( *status==KEY_NO_EXIST ) {
678 /* TDIM does not exist, so clear error stack and */
679 /* write a TDIM only if result is multi-dimensional */
680 *status = 0;
681 ffcmsg();
682 if( naxis>1 )
683 ffptdm( outfptr, colNo, naxis, naxes, status );
684 }
685 if( *status ) {
686 /* Either some other error happened in ffgcrd */
687 /* or one happened in ffptdm */
688 ffcprs();
689 FFUNLOCK;
690 return( *status );
691 }
692
693 }
694
695 if( colNo>0 ) {
696
697 /* Output column exists (now)... put results into it */
698
699 int anyNull = 0;
700 int nPerLp, i;
701 long totaln;
702
703 ffgkyj(infptr, "NAXIS2", &totaln, 0, status);
704
705 /*************************************/
706 /* Create new iterator Output Column */
707 /*************************************/
708
709 col_cnt = gParse.nCols;
710 if( allocateCol( col_cnt, status ) ) {
711 ffcprs();
712 FFUNLOCK;
713 return( *status );
714 }
715
716 fits_iter_set_by_num( gParse.colData+col_cnt, outfptr,
717 colNo, 0, OutputCol );
718 gParse.nCols++;
719
720 for( i=0; i<nRngs; i++ ) {
721 Info.dataPtr = NULL;
722 Info.maxRows = end[i]-start[i]+1;
723
724 /*
725 If there is only 1 range, and it includes all the rows,
726 and there are 10 or more rows, then set nPerLp = 0 so
727 that the iterator function will dynamically choose the
728 most efficient number of rows to process in each loop.
729 Otherwise, set nPerLp to the number of rows in this range.
730 */
731
732 if( (Info.maxRows >= 10) && (nRngs == 1) &&
733 (start[0] == 1) && (end[0] == totaln))
734 nPerLp = 0;
735 else
736 nPerLp = Info.maxRows;
737
738 if( ffiter( gParse.nCols, gParse.colData, start[i]-1,
739 nPerLp, parse_data, (void*)&Info, status ) == -1 )
740 *status = 0;
741 else if( *status ) {
742 ffcprs();
743 FFUNLOCK;
744 return( *status );
745 }
746 if( Info.anyNull ) anyNull = 1;
747 }
748
749 if( newNullKwd && !anyNull ) {
750 ffdkey( outfptr, nullKwd, status );
751 }
752
753 } else {
754
755 /* Put constant result into keyword */
756
757 result = gParse.Nodes + gParse.resultNode;
758 switch( Info.datatype ) {
759 case TDOUBLE:
760 ffukyd( outfptr, parName, result->value.data.dbl, 15,
761 parInfo, status );
762 break;
763 case TLONG:
764 ffukyj( outfptr, parName, result->value.data.lng, parInfo, status );
765 break;
766 case TLOGICAL:
767 ffukyl( outfptr, parName, result->value.data.log, parInfo, status );
768 break;
769 case TBIT:
770 case TSTRING:
771 if (fits_strcasecmp(parName,"HISTORY") == 0) {
772 ffphis( outfptr, result->value.data.str, status);
773 } else if (fits_strcasecmp(parName,"COMMENT") == 0) {
774 ffpcom( outfptr, result->value.data.str, status);
775 } else {
776 ffukys( outfptr, parName, result->value.data.str, parInfo, status );
777 }
778 break;
779 }
780 }
781
782 ffcprs();
783 FFUNLOCK;
784 return( *status );
785 }
786
787 /*--------------------------------------------------------------------------*/
fftexp(fitsfile * fptr,char * expr,int maxdim,int * datatype,long * nelem,int * naxis,long * naxes,int * status)788 int fftexp( fitsfile *fptr, /* I - Input FITS file */
789 char *expr, /* I - Arithmetic expression */
790 int maxdim, /* I - Max Dimension of naxes */
791 int *datatype, /* O - Data type of result */
792 long *nelem, /* O - Vector length of result */
793 int *naxis, /* O - # of dimensions of result */
794 long *naxes, /* O - Size of each dimension */
795 int *status ) /* O - Error status */
796 /* */
797 /* Evaluate the given expression and return information on the result. */
798 /*--------------------------------------------------------------------------*/
799 {
800 FFLOCK;
801 ffiprs( fptr, 0, expr, maxdim, datatype, nelem, naxis, naxes, status );
802 ffcprs();
803 FFUNLOCK;
804 return( *status );
805 }
806
807 /*--------------------------------------------------------------------------*/
ffiprs(fitsfile * fptr,int compressed,char * expr,int maxdim,int * datatype,long * nelem,int * naxis,long * naxes,int * status)808 int ffiprs( fitsfile *fptr, /* I - Input FITS file */
809 int compressed, /* I - Is FITS file hkunexpanded? */
810 char *expr, /* I - Arithmetic expression */
811 int maxdim, /* I - Max Dimension of naxes */
812 int *datatype, /* O - Data type of result */
813 long *nelem, /* O - Vector length of result */
814 int *naxis, /* O - # of dimensions of result */
815 long *naxes, /* O - Size of each dimension */
816 int *status ) /* O - Error status */
817 /* */
818 /* Initialize the parser and determine what type of result the expression */
819 /* produces. */
820 /*--------------------------------------------------------------------------*/
821 {
822 Node *result;
823 int i,lexpr, tstatus = 0;
824 int xaxis, bitpix;
825 long xaxes[9];
826 static iteratorCol dmyCol;
827
828 if( *status ) return( *status );
829
830 /* make sure all internal structures for this HDU are current */
831 if ( ffrdef(fptr, status) ) return(*status);
832
833 /* Initialize the Parser structure */
834
835 gParse.def_fptr = fptr;
836 gParse.compressed = compressed;
837 gParse.nCols = 0;
838 gParse.colData = NULL;
839 gParse.varData = NULL;
840 gParse.getData = find_column;
841 gParse.loadData = load_column;
842 gParse.Nodes = NULL;
843 gParse.nNodesAlloc= 0;
844 gParse.nNodes = 0;
845 gParse.hdutype = 0;
846 gParse.status = 0;
847
848 fits_get_hdu_type(fptr, &gParse.hdutype, status );
849
850 if (gParse.hdutype == IMAGE_HDU) {
851
852 fits_get_img_param(fptr, 9, &bitpix, &xaxis, xaxes, status);
853 if (*status) {
854 ffpmsg("ffiprs: unable to get image dimensions");
855 return( *status );
856 }
857 gParse.totalRows = xaxis > 0 ? 1 : 0;
858 for (i = 0; i < xaxis; ++i)
859 gParse.totalRows *= xaxes[i];
860 if (DEBUG_PIXFILTER)
861 printf("naxis=%d, gParse.totalRows=%ld\n", xaxis, gParse.totalRows);
862 }
863 else if( ffgkyj(fptr, "NAXIS2", &gParse.totalRows, 0, &tstatus) )
864 {
865 /* this might be a 1D or null image with no NAXIS2 keyword */
866 gParse.totalRows = 0;
867 }
868
869
870 /* Copy expression into parser... read from file if necessary */
871
872
873 if( expr[0]=='@' ) {
874 if( ffimport_file( expr+1, &gParse.expr, status ) ) return( *status );
875 lexpr = strlen(gParse.expr);
876 } else {
877 lexpr = strlen(expr);
878 gParse.expr = (char*)malloc( (2+lexpr)*sizeof(char));
879 strcpy(gParse.expr,expr);
880 }
881 strcat(gParse.expr + lexpr,"\n");
882 gParse.index = 0;
883 gParse.is_eobuf = 0;
884
885 /* Parse the expression, building the Nodes and determing */
886 /* which columns are needed and what data type is returned */
887
888 ffrestart(NULL);
889 if( ffparse() ) {
890 return( *status = PARSE_SYNTAX_ERR );
891 }
892 /* Check results */
893
894 *status = gParse.status;
895 if( *status ) return(*status);
896
897 if( !gParse.nNodes ) {
898 ffpmsg("Blank expression");
899 return( *status = PARSE_SYNTAX_ERR );
900 }
901 if( !gParse.nCols ) {
902 dmyCol.fptr = fptr; /* This allows iterator to know value of */
903 gParse.colData = &dmyCol; /* fptr when no columns are referenced */
904 }
905
906 result = gParse.Nodes + gParse.resultNode;
907
908 *naxis = result->value.naxis;
909 *nelem = result->value.nelem;
910 for( i=0; i<*naxis && i<maxdim; i++ )
911 naxes[i] = result->value.naxes[i];
912
913 switch( result->type ) {
914 case BOOLEAN:
915 *datatype = TLOGICAL;
916 break;
917 case LONG:
918 *datatype = TLONG;
919 break;
920 case DOUBLE:
921 *datatype = TDOUBLE;
922 break;
923 case BITSTR:
924 *datatype = TBIT;
925 break;
926 case STRING:
927 *datatype = TSTRING;
928 break;
929 default:
930 *datatype = 0;
931 ffpmsg("Bad return data type");
932 *status = gParse.status = PARSE_BAD_TYPE;
933 break;
934 }
935 gParse.datatype = *datatype;
936 FREE(gParse.expr);
937
938 if( result->operation==CONST_OP ) *nelem = - *nelem;
939 return(*status);
940 }
941
942 /*--------------------------------------------------------------------------*/
ffcprs(void)943 void ffcprs( void ) /* No parameters */
944 /* */
945 /* Clear the parser, making it ready to accept a new expression. */
946 /*--------------------------------------------------------------------------*/
947 {
948 int col, node, i;
949
950 if( gParse.nCols > 0 ) {
951 FREE( gParse.colData );
952 for( col=0; col<gParse.nCols; col++ ) {
953 if( gParse.varData[col].undef == NULL ) continue;
954 if( gParse.varData[col].type == BITSTR )
955 FREE( ((char**)gParse.varData[col].data)[0] );
956 free( gParse.varData[col].undef );
957 }
958 FREE( gParse.varData );
959 gParse.nCols = 0;
960 }
961
962 if( gParse.nNodes > 0 ) {
963 node = gParse.nNodes;
964 while( node-- ) {
965 if( gParse.Nodes[node].operation==gtifilt_fct ) {
966 i = gParse.Nodes[node].SubNodes[0];
967 if (gParse.Nodes[ i ].value.data.ptr)
968 FREE( gParse.Nodes[ i ].value.data.ptr );
969 }
970 else if( gParse.Nodes[node].operation==regfilt_fct ) {
971 i = gParse.Nodes[node].SubNodes[0];
972 fits_free_region( (SAORegion *)gParse.Nodes[ i ].value.data.ptr );
973 }
974 }
975 gParse.nNodes = 0;
976 }
977 if( gParse.Nodes ) free( gParse.Nodes );
978 gParse.Nodes = NULL;
979
980 gParse.hdutype = ANY_HDU;
981 gParse.pixFilter = 0;
982 }
983
984 /*---------------------------------------------------------------------------*/
parse_data(long totalrows,long offset,long firstrow,long nrows,int nCols,iteratorCol * colData,void * userPtr)985 int parse_data( long totalrows, /* I - Total rows to be processed */
986 long offset, /* I - Number of rows skipped at start*/
987 long firstrow, /* I - First row of this iteration */
988 long nrows, /* I - Number of rows in this iter */
989 int nCols, /* I - Number of columns in use */
990 iteratorCol *colData, /* IO- Column information/data */
991 void *userPtr ) /* I - Data handling instructions */
992 /* */
993 /* Iterator work function which calls the parser and copies the results */
994 /* into either an OutputCol or a data pointer supplied in the userPtr */
995 /* structure. */
996 /*---------------------------------------------------------------------------*/
997 {
998 int status, constant=0, anyNullThisTime=0;
999 long jj, kk, idx, remain, ntodo;
1000 Node *result;
1001 iteratorCol * outcol;
1002
1003 /* declare variables static to preserve their values between calls */
1004 static void *Data, *Null;
1005 static int datasize;
1006 static long lastRow, repeat, resDataSize;
1007 static LONGLONG jnull;
1008 static parseInfo *userInfo;
1009 static long zeros[4] = {0,0,0,0};
1010
1011 if (DEBUG_PIXFILTER)
1012 printf("parse_data(total=%ld, offset=%ld, first=%ld, rows=%ld, cols=%d)\n",
1013 totalrows, offset, firstrow, nrows, nCols);
1014 /*--------------------------------------------------------*/
1015 /* Initialization procedures: execute on the first call */
1016 /*--------------------------------------------------------*/
1017 outcol = colData + (nCols - 1);
1018 if (firstrow == offset+1)
1019 {
1020 userInfo = (parseInfo*)userPtr;
1021 userInfo->anyNull = 0;
1022
1023 if( userInfo->maxRows>0 )
1024 userInfo->maxRows = minvalue(totalrows,userInfo->maxRows);
1025 else if( userInfo->maxRows<0 )
1026 userInfo->maxRows = totalrows;
1027 else
1028 userInfo->maxRows = nrows;
1029
1030 lastRow = firstrow + userInfo->maxRows - 1;
1031
1032 if( userInfo->dataPtr==NULL ) {
1033
1034 if( outcol->iotype == InputCol ) {
1035 ffpmsg("Output column for parser results not found!");
1036 return( PARSE_NO_OUTPUT );
1037 }
1038 /* Data gets set later */
1039 Null = outcol->array;
1040 userInfo->datatype = outcol->datatype;
1041
1042 /* Check for a TNULL/BLANK keyword for output column/image */
1043
1044 status = 0;
1045 jnull = 0;
1046 if (gParse.hdutype == IMAGE_HDU) {
1047 if (gParse.pixFilter->blank)
1048 jnull = (LONGLONG) gParse.pixFilter->blank;
1049 }
1050 else {
1051 ffgknjj( outcol->fptr, "TNULL", outcol->colnum,
1052 1, &jnull, (int*)&jj, &status );
1053
1054 if( status==BAD_INTKEY ) {
1055 /* Probably ASCII table with text TNULL keyword */
1056 switch( userInfo->datatype ) {
1057 case TSHORT: jnull = (LONGLONG) SHRT_MIN; break;
1058 case TINT: jnull = (LONGLONG) INT_MIN; break;
1059 case TLONG: jnull = (LONGLONG) LONG_MIN; break;
1060 }
1061 }
1062 }
1063 repeat = outcol->repeat;
1064 /*
1065 if (DEBUG_PIXFILTER)
1066 printf("parse_data: using null value %ld\n", jnull);
1067 */
1068 } else {
1069
1070 Data = userInfo->dataPtr;
1071 Null = (userInfo->nullPtr ? userInfo->nullPtr : zeros);
1072 repeat = gParse.Nodes[gParse.resultNode].value.nelem;
1073
1074 }
1075
1076 /* Determine the size of each element of the returned result */
1077
1078 switch( userInfo->datatype ) {
1079 case TBIT: /* Fall through to TBYTE */
1080 case TLOGICAL: /* Fall through to TBYTE */
1081 case TBYTE: datasize = sizeof(char); break;
1082 case TSHORT: datasize = sizeof(short); break;
1083 case TINT: datasize = sizeof(int); break;
1084 case TLONG: datasize = sizeof(long); break;
1085 case TLONGLONG: datasize = sizeof(LONGLONG); break;
1086 case TFLOAT: datasize = sizeof(float); break;
1087 case TDOUBLE: datasize = sizeof(double); break;
1088 case TSTRING: datasize = sizeof(char*); break;
1089 }
1090
1091 /* Determine the size of each element of the calculated result */
1092 /* (only matters for numeric/logical data) */
1093
1094 switch( gParse.Nodes[gParse.resultNode].type ) {
1095 case BOOLEAN: resDataSize = sizeof(char); break;
1096 case LONG: resDataSize = sizeof(long); break;
1097 case DOUBLE: resDataSize = sizeof(double); break;
1098 }
1099 }
1100
1101 /*-------------------------------------------*/
1102 /* Main loop: process all the rows of data */
1103 /*-------------------------------------------*/
1104
1105 /* If writing to output column, set first element to appropriate */
1106 /* null value. If no NULLs encounter, zero out before returning. */
1107 /*
1108 if (DEBUG_PIXFILTER)
1109 printf("parse_data: using null value %ld\n", jnull);
1110 */
1111
1112 if( userInfo->dataPtr == NULL ) {
1113 /* First, reset Data pointer to start of output array */
1114 Data = (char*) outcol->array + datasize;
1115
1116 switch( userInfo->datatype ) {
1117 case TLOGICAL: *(char *)Null = 'U'; break;
1118 case TBYTE: *(char *)Null = (char )jnull; break;
1119 case TSHORT: *(short *)Null = (short)jnull; break;
1120 case TINT: *(int *)Null = (int )jnull; break;
1121 case TLONG: *(long *)Null = (long )jnull; break;
1122 case TLONGLONG: *(LONGLONG *)Null = (LONGLONG )jnull; break;
1123 case TFLOAT: *(float *)Null = FLOATNULLVALUE; break;
1124 case TDOUBLE: *(double*)Null = DOUBLENULLVALUE; break;
1125 case TSTRING: (*(char **)Null)[0] = '\1';
1126 (*(char **)Null)[1] = '\0'; break;
1127 }
1128 }
1129
1130 /* Alter nrows in case calling routine didn't want to do all rows */
1131
1132 nrows = minvalue(nrows,lastRow-firstrow+1);
1133
1134 Setup_DataArrays( nCols, colData, firstrow, nrows );
1135
1136 /* Parser allocates arrays for each column and calculation it performs. */
1137 /* Limit number of rows processed during each pass to reduce memory */
1138 /* requirements... In most cases, iterator will limit rows to less */
1139 /* than 2500 rows per iteration, so this is really only relevant for */
1140 /* hk-compressed files which must be decompressed in memory and sent */
1141 /* whole to parse_data in a single iteration. */
1142
1143 remain = nrows;
1144 while( remain ) {
1145 ntodo = minvalue(remain,2500);
1146 Evaluate_Parser ( firstrow, ntodo );
1147 if( gParse.status ) break;
1148
1149 firstrow += ntodo;
1150 remain -= ntodo;
1151
1152 /* Copy results into data array */
1153
1154 result = gParse.Nodes + gParse.resultNode;
1155 if( result->operation==CONST_OP ) constant = 1;
1156
1157 switch( result->type ) {
1158
1159 case BOOLEAN:
1160 case LONG:
1161 case DOUBLE:
1162 if( constant ) {
1163 char undef=0;
1164 for( kk=0; kk<ntodo; kk++ )
1165 for( jj=0; jj<repeat; jj++ )
1166 ffcvtn( gParse.datatype,
1167 &(result->value.data),
1168 &undef, result->value.nelem /* 1 */,
1169 userInfo->datatype, Null,
1170 (char*)Data + (kk*repeat+jj)*datasize,
1171 &anyNullThisTime, &gParse.status );
1172 } else {
1173 if ( repeat == result->value.nelem ) {
1174 ffcvtn( gParse.datatype,
1175 result->value.data.ptr,
1176 result->value.undef,
1177 result->value.nelem*ntodo,
1178 userInfo->datatype, Null, Data,
1179 &anyNullThisTime, &gParse.status );
1180 } else if( result->value.nelem == 1 ) {
1181 for( kk=0; kk<ntodo; kk++ )
1182 for( jj=0; jj<repeat; jj++ ) {
1183 ffcvtn( gParse.datatype,
1184 (char*)result->value.data.ptr + kk*resDataSize,
1185 (char*)result->value.undef + kk,
1186 1, userInfo->datatype, Null,
1187 (char*)Data + (kk*repeat+jj)*datasize,
1188 &anyNullThisTime, &gParse.status );
1189 }
1190 } else {
1191 int nCopy;
1192 nCopy = minvalue( repeat, result->value.nelem );
1193 for( kk=0; kk<ntodo; kk++ ) {
1194 ffcvtn( gParse.datatype,
1195 (char*)result->value.data.ptr
1196 + kk*result->value.nelem*resDataSize,
1197 (char*)result->value.undef
1198 + kk*result->value.nelem,
1199 nCopy, userInfo->datatype, Null,
1200 (char*)Data + (kk*repeat)*datasize,
1201 &anyNullThisTime, &gParse.status );
1202 if( nCopy < repeat ) {
1203 memset( (char*)Data + (kk*repeat+nCopy)*datasize,
1204 0, (repeat-nCopy)*datasize);
1205 }
1206 }
1207
1208 }
1209 if( result->operation>0 ) {
1210 FREE( result->value.data.ptr );
1211 }
1212 }
1213 if( gParse.status==OVERFLOW_ERR ) {
1214 gParse.status = NUM_OVERFLOW;
1215 ffpmsg("Numerical overflow while converting expression to necessary datatype");
1216 }
1217 break;
1218
1219 case BITSTR:
1220 switch( userInfo->datatype ) {
1221 case TBYTE:
1222 idx = -1;
1223 for( kk=0; kk<ntodo; kk++ ) {
1224 for( jj=0; jj<result->value.nelem; jj++ ) {
1225 if( jj%8 == 0 )
1226 ((char*)Data)[++idx] = 0;
1227 if( constant ) {
1228 if( result->value.data.str[jj]=='1' )
1229 ((char*)Data)[idx] |= 128>>(jj%8);
1230 } else {
1231 if( result->value.data.strptr[kk][jj]=='1' )
1232 ((char*)Data)[idx] |= 128>>(jj%8);
1233 }
1234 }
1235 }
1236 break;
1237 case TBIT:
1238 case TLOGICAL:
1239 if( constant ) {
1240 for( kk=0; kk<ntodo; kk++ )
1241 for( jj=0; jj<result->value.nelem; jj++ ) {
1242 ((char*)Data)[ jj+kk*result->value.nelem ] =
1243 ( result->value.data.str[jj]=='1' );
1244 }
1245 } else {
1246 for( kk=0; kk<ntodo; kk++ )
1247 for( jj=0; jj<result->value.nelem; jj++ ) {
1248 ((char*)Data)[ jj+kk*result->value.nelem ] =
1249 ( result->value.data.strptr[kk][jj]=='1' );
1250 }
1251 }
1252 break;
1253 case TSTRING:
1254 if( constant ) {
1255 for( jj=0; jj<ntodo; jj++ ) {
1256 strcpy( ((char**)Data)[jj], result->value.data.str );
1257 }
1258 } else {
1259 for( jj=0; jj<ntodo; jj++ ) {
1260 strcpy( ((char**)Data)[jj], result->value.data.strptr[jj] );
1261 }
1262 }
1263 break;
1264 default:
1265 ffpmsg("Cannot convert bit expression to desired type.");
1266 gParse.status = PARSE_BAD_TYPE;
1267 break;
1268 }
1269 if( result->operation>0 ) {
1270 FREE( result->value.data.strptr[0] );
1271 FREE( result->value.data.strptr );
1272 }
1273 break;
1274
1275 case STRING:
1276 if( userInfo->datatype==TSTRING ) {
1277 if( constant ) {
1278 for( jj=0; jj<ntodo; jj++ )
1279 strcpy( ((char**)Data)[jj], result->value.data.str );
1280 } else {
1281 for( jj=0; jj<ntodo; jj++ )
1282 if( result->value.undef[jj] ) {
1283 anyNullThisTime = 1;
1284 strcpy( ((char**)Data)[jj],
1285 *(char **)Null );
1286 } else {
1287 strcpy( ((char**)Data)[jj],
1288 result->value.data.strptr[jj] );
1289 }
1290 }
1291 } else {
1292 ffpmsg("Cannot convert string expression to desired type.");
1293 gParse.status = PARSE_BAD_TYPE;
1294 }
1295 if( result->operation>0 ) {
1296 FREE( result->value.data.strptr[0] );
1297 FREE( result->value.data.strptr );
1298 }
1299 break;
1300 }
1301
1302 if( gParse.status ) break;
1303
1304 /* Increment Data to point to where the next block should go */
1305
1306 if( result->type==BITSTR && userInfo->datatype==TBYTE )
1307 Data = (char*)Data
1308 + datasize * ( (result->value.nelem+7)/8 ) * ntodo;
1309 else if( result->type==STRING )
1310 Data = (char*)Data + datasize * ntodo;
1311 else
1312 Data = (char*)Data + datasize * ntodo * repeat;
1313 }
1314
1315 /* If no NULLs encountered during this pass, set Null value to */
1316 /* zero to make the writing of the output column data faster */
1317
1318 if( anyNullThisTime )
1319 userInfo->anyNull = 1;
1320 else if( userInfo->dataPtr == NULL ) {
1321 if( userInfo->datatype == TSTRING )
1322 memcpy( *(char **)Null, zeros, 2 );
1323 else
1324 memcpy( Null, zeros, datasize );
1325 }
1326
1327 /*-------------------------------------------------------*/
1328 /* Clean up procedures: after processing all the rows */
1329 /*-------------------------------------------------------*/
1330
1331 /* if the calling routine specified that only a limited number */
1332 /* of rows in the table should be processed, return a value of -1 */
1333 /* once all the rows have been done, if no other error occurred. */
1334
1335 if (gParse.hdutype != IMAGE_HDU && firstrow - 1 == lastRow) {
1336 if (!gParse.status && userInfo->maxRows<totalrows) {
1337 return (-1);
1338 }
1339 }
1340
1341 return(gParse.status); /* return successful status */
1342 }
1343
Setup_DataArrays(int nCols,iteratorCol * cols,long fRow,long nRows)1344 static void Setup_DataArrays( int nCols, iteratorCol *cols,
1345 long fRow, long nRows )
1346 /***********************************************************************/
1347 /* Setup the varData array in gParse to contain the fits column data. */
1348 /* Then, allocate and initialize the necessary UNDEF arrays for each */
1349 /* column used by the parser. */
1350 /***********************************************************************/
1351 {
1352 int i;
1353 long nelem, len, row, idx;
1354 char **bitStrs;
1355 char **sptr;
1356 char *barray;
1357 long *iarray;
1358 double *rarray;
1359 char msg[80];
1360
1361 gParse.firstDataRow = fRow;
1362 gParse.nDataRows = nRows;
1363
1364 /* Resize and fill in UNDEF arrays for each column */
1365
1366 for( i=0; i<nCols; i++ ) {
1367
1368 iteratorCol *icol = cols + i;
1369 DataInfo *varData = gParse.varData + i;
1370
1371 if( icol->iotype == OutputCol ) continue;
1372
1373 nelem = varData->nelem;
1374 len = nelem * nRows;
1375
1376 switch ( varData->type ) {
1377
1378 case BITSTR:
1379 /* No need for UNDEF array, but must make string DATA array */
1380 len = (nelem+1)*nRows; /* Count '\0' */
1381 bitStrs = (char**)varData->data;
1382 if( bitStrs ) FREE( bitStrs[0] );
1383 free( bitStrs );
1384 bitStrs = (char**)malloc( nRows*sizeof(char*) );
1385 if( bitStrs==NULL ) {
1386 varData->data = varData->undef = NULL;
1387 gParse.status = MEMORY_ALLOCATION;
1388 break;
1389 }
1390 bitStrs[0] = (char*)malloc( len*sizeof(char) );
1391 if( bitStrs[0]==NULL ) {
1392 free( bitStrs );
1393 varData->data = varData->undef = NULL;
1394 gParse.status = MEMORY_ALLOCATION;
1395 break;
1396 }
1397
1398 for( row=0; row<nRows; row++ ) {
1399 bitStrs[row] = bitStrs[0] + row*(nelem+1);
1400 idx = (row)*( (nelem+7)/8 ) + 1;
1401 for(len=0; len<nelem; len++) {
1402 if( ((char*)icol->array)[idx] & (1<<(7-len%8)) )
1403 bitStrs[row][len] = '1';
1404 else
1405 bitStrs[row][len] = '0';
1406 if( len%8==7 ) idx++;
1407 }
1408 bitStrs[row][len] = '\0';
1409 }
1410 varData->undef = (char*)bitStrs;
1411 varData->data = (char*)bitStrs;
1412 break;
1413
1414 case STRING:
1415 sptr = (char**)icol->array;
1416 if (varData->undef)
1417 free( varData->undef );
1418 varData->undef = (char*)malloc( nRows*sizeof(char) );
1419 if( varData->undef==NULL ) {
1420 gParse.status = MEMORY_ALLOCATION;
1421 break;
1422 }
1423 row = nRows;
1424 while( row-- )
1425 varData->undef[row] =
1426 ( **sptr != '\0' && FSTRCMP( sptr[0], sptr[row+1] )==0 );
1427 varData->data = sptr + 1;
1428 break;
1429
1430 case BOOLEAN:
1431 barray = (char*)icol->array;
1432 if (varData->undef)
1433 free( varData->undef );
1434 varData->undef = (char*)malloc( len*sizeof(char) );
1435 if( varData->undef==NULL ) {
1436 gParse.status = MEMORY_ALLOCATION;
1437 break;
1438 }
1439 while( len-- ) {
1440 varData->undef[len] =
1441 ( barray[0]!=0 && barray[0]==barray[len+1] );
1442 }
1443 varData->data = barray + 1;
1444 break;
1445
1446 case LONG:
1447 iarray = (long*)icol->array;
1448 if (varData->undef)
1449 free( varData->undef );
1450 varData->undef = (char*)malloc( len*sizeof(char) );
1451 if( varData->undef==NULL ) {
1452 gParse.status = MEMORY_ALLOCATION;
1453 break;
1454 }
1455 while( len-- ) {
1456 varData->undef[len] =
1457 ( iarray[0]!=0L && iarray[0]==iarray[len+1] );
1458 }
1459 varData->data = iarray + 1;
1460 break;
1461
1462 case DOUBLE:
1463 rarray = (double*)icol->array;
1464 if (varData->undef)
1465 free( varData->undef );
1466 varData->undef = (char*)malloc( len*sizeof(char) );
1467 if( varData->undef==NULL ) {
1468 gParse.status = MEMORY_ALLOCATION;
1469 break;
1470 }
1471 while( len-- ) {
1472 varData->undef[len] =
1473 ( rarray[0]!=0.0 && rarray[0]==rarray[len+1]);
1474 }
1475 varData->data = rarray + 1;
1476 break;
1477
1478 default:
1479 snprintf(msg, 80, "SetupDataArrays, unhandled type %d\n",
1480 varData->type);
1481 ffpmsg(msg);
1482 }
1483
1484 if( gParse.status ) { /* Deallocate NULL arrays of previous columns */
1485 while( i-- ) {
1486 varData = gParse.varData + i;
1487 if( varData->type==BITSTR )
1488 FREE( ((char**)varData->data)[0] );
1489 FREE( varData->undef );
1490 varData->undef = NULL;
1491 }
1492 return;
1493 }
1494 }
1495 }
1496
1497 /*--------------------------------------------------------------------------*/
ffcvtn(int inputType,void * input,char * undef,long ntodo,int outputType,void * nulval,void * output,int * anynull,int * status)1498 int ffcvtn( int inputType, /* I - Data type of input array */
1499 void *input, /* I - Input array of type inputType */
1500 char *undef, /* I - Array of flags indicating UNDEF elems */
1501 long ntodo, /* I - Number of elements to process */
1502 int outputType, /* I - Data type of output array */
1503 void *nulval, /* I - Ptr to value to use for UNDEF elements */
1504 void *output, /* O - Output array of type outputType */
1505 int *anynull, /* O - Any nulls flagged? */
1506 int *status ) /* O - Error status */
1507 /* */
1508 /* Convert an array of any input data type to an array of any output */
1509 /* data type, using an array of UNDEF flags to assign nulvals to */
1510 /*--------------------------------------------------------------------------*/
1511 {
1512 long i;
1513
1514 switch( outputType ) {
1515
1516 case TLOGICAL:
1517 switch( inputType ) {
1518 case TLOGICAL:
1519 case TBYTE:
1520 for( i=0; i<ntodo; i++ )
1521 if( ((unsigned char*)input)[i] )
1522 ((unsigned char*)output)[i] = 1;
1523 else
1524 ((unsigned char*)output)[i] = 0;
1525 break;
1526 case TSHORT:
1527 for( i=0; i<ntodo; i++ )
1528 if( ((short*)input)[i] )
1529 ((unsigned char*)output)[i] = 1;
1530 else
1531 ((unsigned char*)output)[i] = 0;
1532 break;
1533 case TLONG:
1534 for( i=0; i<ntodo; i++ )
1535 if( ((long*)input)[i] )
1536 ((unsigned char*)output)[i] = 1;
1537 else
1538 ((unsigned char*)output)[i] = 0;
1539 break;
1540 case TFLOAT:
1541 for( i=0; i<ntodo; i++ )
1542 if( ((float*)input)[i] )
1543 ((unsigned char*)output)[i] = 1;
1544 else
1545 ((unsigned char*)output)[i] = 0;
1546 break;
1547 case TDOUBLE:
1548 for( i=0; i<ntodo; i++ )
1549 if( ((double*)input)[i] )
1550 ((unsigned char*)output)[i] = 1;
1551 else
1552 ((unsigned char*)output)[i] = 0;
1553 break;
1554 default:
1555 *status = BAD_DATATYPE;
1556 break;
1557 }
1558 for(i=0;i<ntodo;i++) {
1559 if( undef[i] ) {
1560 ((unsigned char*)output)[i] = *(unsigned char*)nulval;
1561 *anynull = 1;
1562 }
1563 }
1564 break;
1565
1566 case TBYTE:
1567 switch( inputType ) {
1568 case TLOGICAL:
1569 case TBYTE:
1570 for( i=0; i<ntodo; i++ )
1571 ((unsigned char*)output)[i] = ((unsigned char*)input)[i];
1572 break;
1573 case TSHORT:
1574 fffi2i1((short*)input,ntodo,1.,0.,0,0,0,NULL,NULL,(unsigned char*)output,status);
1575 break;
1576 case TLONG:
1577 for (i = 0; i < ntodo; i++) {
1578 if( undef[i] ) {
1579 ((unsigned char*)output)[i] = *(unsigned char*)nulval;
1580 *anynull = 1;
1581 } else {
1582 if( ((long*)input)[i] < 0 ) {
1583 *status = OVERFLOW_ERR;
1584 ((unsigned char*)output)[i] = 0;
1585 } else if( ((long*)input)[i] > UCHAR_MAX ) {
1586 *status = OVERFLOW_ERR;
1587 ((unsigned char*)output)[i] = UCHAR_MAX;
1588 } else
1589 ((unsigned char*)output)[i] =
1590 (unsigned char) ((long*)input)[i];
1591 }
1592 }
1593 return( *status );
1594 case TFLOAT:
1595 fffr4i1((float*)input,ntodo,1.,0.,0,0,NULL,NULL,
1596 (unsigned char*)output,status);
1597 break;
1598 case TDOUBLE:
1599 fffr8i1((double*)input,ntodo,1.,0.,0,0,NULL,NULL,
1600 (unsigned char*)output,status);
1601 break;
1602 default:
1603 *status = BAD_DATATYPE;
1604 break;
1605 }
1606 for(i=0;i<ntodo;i++) {
1607 if( undef[i] ) {
1608 ((unsigned char*)output)[i] = *(unsigned char*)nulval;
1609 *anynull = 1;
1610 }
1611 }
1612 break;
1613
1614 case TSHORT:
1615 switch( inputType ) {
1616 case TLOGICAL:
1617 case TBYTE:
1618 for( i=0; i<ntodo; i++ )
1619 ((short*)output)[i] = ((unsigned char*)input)[i];
1620 break;
1621 case TSHORT:
1622 for( i=0; i<ntodo; i++ )
1623 ((short*)output)[i] = ((short*)input)[i];
1624 break;
1625 case TLONG:
1626 for (i = 0; i < ntodo; i++) {
1627 if( undef[i] ) {
1628 ((short*)output)[i] = *(short*)nulval;
1629 *anynull = 1;
1630 } else {
1631 if( ((long*)input)[i] < SHRT_MIN ) {
1632 *status = OVERFLOW_ERR;
1633 ((short*)output)[i] = SHRT_MIN;
1634 } else if ( ((long*)input)[i] > SHRT_MAX ) {
1635 *status = OVERFLOW_ERR;
1636 ((short*)output)[i] = SHRT_MAX;
1637 } else
1638 ((short*)output)[i] = (short) ((long*)input)[i];
1639 }
1640 }
1641 return( *status );
1642 case TFLOAT:
1643 fffr4i2((float*)input,ntodo,1.,0.,0,0,NULL,NULL,
1644 (short*)output,status);
1645 break;
1646 case TDOUBLE:
1647 fffr8i2((double*)input,ntodo,1.,0.,0,0,NULL,NULL,
1648 (short*)output,status);
1649 break;
1650 default:
1651 *status = BAD_DATATYPE;
1652 break;
1653 }
1654 for(i=0;i<ntodo;i++) {
1655 if( undef[i] ) {
1656 ((short*)output)[i] = *(short*)nulval;
1657 *anynull = 1;
1658 }
1659 }
1660 break;
1661
1662 case TINT:
1663 switch( inputType ) {
1664 case TLOGICAL:
1665 case TBYTE:
1666 for( i=0; i<ntodo; i++ )
1667 ((int*)output)[i] = ((unsigned char*)input)[i];
1668 break;
1669 case TSHORT:
1670 for( i=0; i<ntodo; i++ )
1671 ((int*)output)[i] = ((short*)input)[i];
1672 break;
1673 case TLONG:
1674 for( i=0; i<ntodo; i++ )
1675 ((int*)output)[i] = ((long*)input)[i];
1676 break;
1677 case TFLOAT:
1678 fffr4int((float*)input,ntodo,1.,0.,0,0,NULL,NULL,
1679 (int*)output,status);
1680 break;
1681 case TDOUBLE:
1682 fffr8int((double*)input,ntodo,1.,0.,0,0,NULL,NULL,
1683 (int*)output,status);
1684 break;
1685 default:
1686 *status = BAD_DATATYPE;
1687 break;
1688 }
1689 for(i=0;i<ntodo;i++) {
1690 if( undef[i] ) {
1691 ((int*)output)[i] = *(int*)nulval;
1692 *anynull = 1;
1693 }
1694 }
1695 break;
1696
1697 case TLONG:
1698 switch( inputType ) {
1699 case TLOGICAL:
1700 case TBYTE:
1701 for( i=0; i<ntodo; i++ )
1702 ((long*)output)[i] = ((unsigned char*)input)[i];
1703 break;
1704 case TSHORT:
1705 for( i=0; i<ntodo; i++ )
1706 ((long*)output)[i] = ((short*)input)[i];
1707 break;
1708 case TLONG:
1709 for( i=0; i<ntodo; i++ )
1710 ((long*)output)[i] = ((long*)input)[i];
1711 break;
1712 case TFLOAT:
1713 fffr4i4((float*)input,ntodo,1.,0.,0,0,NULL,NULL,
1714 (long*)output,status);
1715 break;
1716 case TDOUBLE:
1717 fffr8i4((double*)input,ntodo,1.,0.,0,0,NULL,NULL,
1718 (long*)output,status);
1719 break;
1720 default:
1721 *status = BAD_DATATYPE;
1722 break;
1723 }
1724 for(i=0;i<ntodo;i++) {
1725 if( undef[i] ) {
1726 ((long*)output)[i] = *(long*)nulval;
1727 *anynull = 1;
1728 }
1729 }
1730 break;
1731
1732 case TLONGLONG:
1733 switch( inputType ) {
1734 case TLOGICAL:
1735 case TBYTE:
1736 for( i=0; i<ntodo; i++ )
1737 ((LONGLONG*)output)[i] = ((unsigned char*)input)[i];
1738 break;
1739 case TSHORT:
1740 for( i=0; i<ntodo; i++ )
1741 ((LONGLONG*)output)[i] = ((short*)input)[i];
1742 break;
1743 case TLONG:
1744 for( i=0; i<ntodo; i++ )
1745 ((LONGLONG*)output)[i] = ((long*)input)[i];
1746 break;
1747 case TFLOAT:
1748 fffr4i8((float*)input,ntodo,1.,0.,0,0,NULL,NULL,
1749 (LONGLONG*)output,status);
1750 break;
1751 case TDOUBLE:
1752 fffr8i8((double*)input,ntodo,1.,0.,0,0,NULL,NULL,
1753 (LONGLONG*)output,status);
1754
1755 break;
1756 default:
1757 *status = BAD_DATATYPE;
1758 break;
1759 }
1760 for(i=0;i<ntodo;i++) {
1761 if( undef[i] ) {
1762 ((LONGLONG*)output)[i] = *(LONGLONG*)nulval;
1763 *anynull = 1;
1764 }
1765 }
1766 break;
1767
1768 case TFLOAT:
1769 switch( inputType ) {
1770 case TLOGICAL:
1771 case TBYTE:
1772 for( i=0; i<ntodo; i++ )
1773 ((float*)output)[i] = ((unsigned char*)input)[i];
1774 break;
1775 case TSHORT:
1776 for( i=0; i<ntodo; i++ )
1777 ((float*)output)[i] = ((short*)input)[i];
1778 break;
1779 case TLONG:
1780 for( i=0; i<ntodo; i++ )
1781 ((float*)output)[i] = (float) ((long*)input)[i];
1782 break;
1783 case TFLOAT:
1784 for( i=0; i<ntodo; i++ )
1785 ((float*)output)[i] = ((float*)input)[i];
1786 break;
1787 case TDOUBLE:
1788 fffr8r4((double*)input,ntodo,1.,0.,0,0,NULL,NULL,
1789 (float*)output,status);
1790 break;
1791 default:
1792 *status = BAD_DATATYPE;
1793 break;
1794 }
1795 for(i=0;i<ntodo;i++) {
1796 if( undef[i] ) {
1797 ((float*)output)[i] = *(float*)nulval;
1798 *anynull = 1;
1799 }
1800 }
1801 break;
1802
1803 case TDOUBLE:
1804 switch( inputType ) {
1805 case TLOGICAL:
1806 case TBYTE:
1807 for( i=0; i<ntodo; i++ )
1808 ((double*)output)[i] = ((unsigned char*)input)[i];
1809 break;
1810 case TSHORT:
1811 for( i=0; i<ntodo; i++ )
1812 ((double*)output)[i] = ((short*)input)[i];
1813 break;
1814 case TLONG:
1815 for( i=0; i<ntodo; i++ )
1816 ((double*)output)[i] = ((long*)input)[i];
1817 break;
1818 case TFLOAT:
1819 for( i=0; i<ntodo; i++ )
1820 ((double*)output)[i] = ((float*)input)[i];
1821 break;
1822 case TDOUBLE:
1823 for( i=0; i<ntodo; i++ )
1824 ((double*)output)[i] = ((double*)input)[i];
1825 break;
1826 default:
1827 *status = BAD_DATATYPE;
1828 break;
1829 }
1830 for(i=0;i<ntodo;i++) {
1831 if( undef[i] ) {
1832 ((double*)output)[i] = *(double*)nulval;
1833 *anynull = 1;
1834 }
1835 }
1836 break;
1837
1838 default:
1839 *status = BAD_DATATYPE;
1840 break;
1841 }
1842
1843 return ( *status );
1844 }
1845
1846 /*---------------------------------------------------------------------------*/
fffrwc(fitsfile * fptr,char * expr,char * timeCol,char * parCol,char * valCol,long ntimes,double * times,char * time_status,int * status)1847 int fffrwc( fitsfile *fptr, /* I - Input FITS file */
1848 char *expr, /* I - Boolean expression */
1849 char *timeCol, /* I - Name of time column */
1850 char *parCol, /* I - Name of parameter column */
1851 char *valCol, /* I - Name of value column */
1852 long ntimes, /* I - Number of distinct times in file */
1853 double *times, /* O - Array of times in file */
1854 char *time_status, /* O - Array of boolean results */
1855 int *status ) /* O - Error status */
1856 /* */
1857 /* Evaluate a boolean expression for each time in a compressed file, */
1858 /* returning an array of flags indicating which times evaluated to TRUE/FALSE*/
1859 /*---------------------------------------------------------------------------*/
1860 {
1861 parseInfo Info;
1862 long alen, width;
1863 int parNo, typecode;
1864 int naxis, constant, nCol=0;
1865 long nelem, naxes[MAXDIMS], elem;
1866 char result;
1867
1868 if( *status ) return( *status );
1869
1870 fits_get_colnum( fptr, CASEINSEN, timeCol, &gParse.timeCol, status );
1871 fits_get_colnum( fptr, CASEINSEN, parCol, &gParse.parCol , status );
1872 fits_get_colnum( fptr, CASEINSEN, valCol, &gParse.valCol, status );
1873 if( *status ) return( *status );
1874
1875 if( ffiprs( fptr, 1, expr, MAXDIMS, &Info.datatype, &nelem,
1876 &naxis, naxes, status ) ) {
1877 ffcprs();
1878 return( *status );
1879 }
1880 if( nelem<0 ) {
1881 constant = 1;
1882 nelem = -nelem;
1883 nCol = gParse.nCols;
1884 gParse.nCols = 0; /* Ignore all column references */
1885 } else
1886 constant = 0;
1887
1888 if( Info.datatype!=TLOGICAL || nelem!=1 ) {
1889 ffcprs();
1890 ffpmsg("Expression does not evaluate to a logical scalar.");
1891 return( *status = PARSE_BAD_TYPE );
1892 }
1893
1894 /*******************************************/
1895 /* Allocate data arrays for each parameter */
1896 /*******************************************/
1897
1898 parNo = gParse.nCols;
1899 while( parNo-- ) {
1900 switch( gParse.colData[parNo].datatype ) {
1901 case TLONG:
1902 if( (gParse.colData[parNo].array =
1903 (long *)malloc( (ntimes+1)*sizeof(long) )) )
1904 ((long*)gParse.colData[parNo].array)[0] = 1234554321;
1905 else
1906 *status = MEMORY_ALLOCATION;
1907 break;
1908 case TDOUBLE:
1909 if( (gParse.colData[parNo].array =
1910 (double *)malloc( (ntimes+1)*sizeof(double) )) )
1911 ((double*)gParse.colData[parNo].array)[0] = DOUBLENULLVALUE;
1912 else
1913 *status = MEMORY_ALLOCATION;
1914 break;
1915 case TSTRING:
1916 if( !fits_get_coltype( fptr, gParse.valCol, &typecode,
1917 &alen, &width, status ) ) {
1918 alen++;
1919 if( (gParse.colData[parNo].array =
1920 (char **)malloc( (ntimes+1)*sizeof(char*) )) ) {
1921 if( (((char **)gParse.colData[parNo].array)[0] =
1922 (char *)malloc( (ntimes+1)*sizeof(char)*alen )) ) {
1923 for( elem=1; elem<=ntimes; elem++ )
1924 ((char **)gParse.colData[parNo].array)[elem] =
1925 ((char **)gParse.colData[parNo].array)[elem-1]+alen;
1926 ((char **)gParse.colData[parNo].array)[0][0] = '\0';
1927 } else {
1928 free( gParse.colData[parNo].array );
1929 *status = MEMORY_ALLOCATION;
1930 }
1931 } else {
1932 *status = MEMORY_ALLOCATION;
1933 }
1934 }
1935 break;
1936 }
1937 if( *status ) {
1938 while( parNo-- ) {
1939 if( gParse.colData[parNo].datatype==TSTRING )
1940 FREE( ((char **)gParse.colData[parNo].array)[0] );
1941 FREE( gParse.colData[parNo].array );
1942 }
1943 return( *status );
1944 }
1945 }
1946
1947 /**********************************************************************/
1948 /* Read data from columns needed for the expression and then parse it */
1949 /**********************************************************************/
1950
1951 if( !uncompress_hkdata( fptr, ntimes, times, status ) ) {
1952 if( constant ) {
1953 result = gParse.Nodes[gParse.resultNode].value.data.log;
1954 elem = ntimes;
1955 while( elem-- ) time_status[elem] = result;
1956 } else {
1957 Info.dataPtr = time_status;
1958 Info.nullPtr = NULL;
1959 Info.maxRows = ntimes;
1960 *status = parse_data( ntimes, 0, 1, ntimes, gParse.nCols,
1961 gParse.colData, (void*)&Info );
1962 }
1963 }
1964
1965 /************/
1966 /* Clean up */
1967 /************/
1968
1969 parNo = gParse.nCols;
1970 while ( parNo-- ) {
1971 if( gParse.colData[parNo].datatype==TSTRING )
1972 FREE( ((char **)gParse.colData[parNo].array)[0] );
1973 FREE( gParse.colData[parNo].array );
1974 }
1975
1976 if( constant ) gParse.nCols = nCol;
1977
1978 ffcprs();
1979 return(*status);
1980 }
1981
1982 /*---------------------------------------------------------------------------*/
uncompress_hkdata(fitsfile * fptr,long ntimes,double * times,int * status)1983 int uncompress_hkdata( fitsfile *fptr,
1984 long ntimes,
1985 double *times,
1986 int *status )
1987 /* */
1988 /* description */
1989 /*---------------------------------------------------------------------------*/
1990 {
1991 char parName[256], *sPtr[1], found[1000];
1992 int parNo, anynul;
1993 long naxis2, row, currelem;
1994 double currtime, newtime;
1995
1996 sPtr[0] = parName;
1997 currelem = 0;
1998 currtime = -1e38;
1999
2000 parNo=gParse.nCols;
2001 while( parNo-- ) found[parNo] = 0;
2002
2003 if( ffgkyj( fptr, "NAXIS2", &naxis2, NULL, status ) ) return( *status );
2004
2005 for( row=1; row<=naxis2; row++ ) {
2006 if( ffgcvd( fptr, gParse.timeCol, row, 1L, 1L, 0.0,
2007 &newtime, &anynul, status ) ) return( *status );
2008 if( newtime != currtime ) {
2009 /* New time encountered... propogate parameters to next row */
2010 if( currelem==ntimes ) {
2011 ffpmsg("Found more unique time stamps than caller indicated");
2012 return( *status = PARSE_BAD_COL );
2013 }
2014 times[currelem++] = currtime = newtime;
2015 parNo = gParse.nCols;
2016 while( parNo-- ) {
2017 switch( gParse.colData[parNo].datatype ) {
2018 case TLONG:
2019 ((long*)gParse.colData[parNo].array)[currelem] =
2020 ((long*)gParse.colData[parNo].array)[currelem-1];
2021 break;
2022 case TDOUBLE:
2023 ((double*)gParse.colData[parNo].array)[currelem] =
2024 ((double*)gParse.colData[parNo].array)[currelem-1];
2025 break;
2026 case TSTRING:
2027 strcpy( ((char **)gParse.colData[parNo].array)[currelem],
2028 ((char **)gParse.colData[parNo].array)[currelem-1] );
2029 break;
2030 }
2031 }
2032 }
2033
2034 if( ffgcvs( fptr, gParse.parCol, row, 1L, 1L, "",
2035 sPtr, &anynul, status ) ) return( *status );
2036 parNo = gParse.nCols;
2037 while( parNo-- )
2038 if( !fits_strcasecmp( parName, gParse.varData[parNo].name ) ) break;
2039
2040 if( parNo>=0 ) {
2041 found[parNo] = 1; /* Flag this parameter as found */
2042 switch( gParse.colData[parNo].datatype ) {
2043 case TLONG:
2044 ffgcvj( fptr, gParse.valCol, row, 1L, 1L,
2045 ((long*)gParse.colData[parNo].array)[0],
2046 ((long*)gParse.colData[parNo].array)+currelem,
2047 &anynul, status );
2048 break;
2049 case TDOUBLE:
2050 ffgcvd( fptr, gParse.valCol, row, 1L, 1L,
2051 ((double*)gParse.colData[parNo].array)[0],
2052 ((double*)gParse.colData[parNo].array)+currelem,
2053 &anynul, status );
2054 break;
2055 case TSTRING:
2056 ffgcvs( fptr, gParse.valCol, row, 1L, 1L,
2057 ((char**)gParse.colData[parNo].array)[0],
2058 ((char**)gParse.colData[parNo].array)+currelem,
2059 &anynul, status );
2060 break;
2061 }
2062 if( *status ) return( *status );
2063 }
2064 }
2065
2066 if( currelem<ntimes ) {
2067 ffpmsg("Found fewer unique time stamps than caller indicated");
2068 return( *status = PARSE_BAD_COL );
2069 }
2070
2071 /* Check for any parameters which were not located in the table */
2072 parNo = gParse.nCols;
2073 while( parNo-- )
2074 if( !found[parNo] ) {
2075 snprintf( parName, 256, "Parameter not found: %-30s",
2076 gParse.varData[parNo].name );
2077 ffpmsg( parName );
2078 *status = PARSE_SYNTAX_ERR;
2079 }
2080 return( *status );
2081 }
2082
2083 /*---------------------------------------------------------------------------*/
ffffrw(fitsfile * fptr,char * expr,long * rownum,int * status)2084 int ffffrw( fitsfile *fptr, /* I - Input FITS file */
2085 char *expr, /* I - Boolean expression */
2086 long *rownum, /* O - First row of table to eval to T */
2087 int *status ) /* O - Error status */
2088 /* */
2089 /* Evaluate a boolean expression, returning the row number of the first */
2090 /* row which evaluates to TRUE */
2091 /*---------------------------------------------------------------------------*/
2092 {
2093 int naxis, constant, dtype;
2094 long nelem, naxes[MAXDIMS];
2095 char result;
2096
2097 if( *status ) return( *status );
2098
2099 FFLOCK;
2100 if( ffiprs( fptr, 0, expr, MAXDIMS, &dtype, &nelem, &naxis,
2101 naxes, status ) ) {
2102 ffcprs();
2103 FFUNLOCK;
2104 return( *status );
2105 }
2106 if( nelem<0 ) {
2107 constant = 1;
2108 nelem = -nelem;
2109 } else
2110 constant = 0;
2111
2112 if( dtype!=TLOGICAL || nelem!=1 ) {
2113 ffcprs();
2114 ffpmsg("Expression does not evaluate to a logical scalar.");
2115 FFUNLOCK;
2116 return( *status = PARSE_BAD_TYPE );
2117 }
2118
2119 *rownum = 0;
2120 if( constant ) { /* No need to call parser... have result from ffiprs */
2121 result = gParse.Nodes[gParse.resultNode].value.data.log;
2122 if( result ) {
2123 /* Make sure there is at least 1 row in table */
2124 ffgnrw( fptr, &nelem, status );
2125 if( nelem )
2126 *rownum = 1;
2127 }
2128 } else {
2129 if( ffiter( gParse.nCols, gParse.colData, 0, 0,
2130 ffffrw_work, (void*)rownum, status ) == -1 )
2131 *status = 0; /* -1 indicates exitted without error before end... OK */
2132 }
2133
2134 ffcprs();
2135 FFUNLOCK;
2136 return(*status);
2137 }
2138
2139 /*---------------------------------------------------------------------------*/
ffffrw_work(long totalrows,long offset,long firstrow,long nrows,int nCols,iteratorCol * colData,void * userPtr)2140 int ffffrw_work(long totalrows, /* I - Total rows to be processed */
2141 long offset, /* I - Number of rows skipped at start*/
2142 long firstrow, /* I - First row of this iteration */
2143 long nrows, /* I - Number of rows in this iter */
2144 int nCols, /* I - Number of columns in use */
2145 iteratorCol *colData, /* IO- Column information/data */
2146 void *userPtr ) /* I - Data handling instructions */
2147 /* */
2148 /* Iterator work function which calls the parser and searches for the */
2149 /* first row which evaluates to TRUE. */
2150 /*---------------------------------------------------------------------------*/
2151 {
2152 long idx;
2153 Node *result;
2154
2155 Evaluate_Parser( firstrow, nrows );
2156
2157 if( !gParse.status ) {
2158
2159 result = gParse.Nodes + gParse.resultNode;
2160 if( result->operation==CONST_OP ) {
2161
2162 if( result->value.data.log ) {
2163 *(long*)userPtr = firstrow;
2164 return( -1 );
2165 }
2166
2167 } else {
2168
2169 for( idx=0; idx<nrows; idx++ )
2170 if( result->value.data.logptr[idx] && !result->value.undef[idx] ) {
2171 *(long*)userPtr = firstrow + idx;
2172 return( -1 );
2173 }
2174 }
2175 }
2176
2177 return( gParse.status );
2178 }
2179
2180
set_image_col_types(fitsfile * fptr,const char * name,int bitpix,DataInfo * varInfo,iteratorCol * colIter)2181 static int set_image_col_types (fitsfile * fptr, const char * name, int bitpix,
2182 DataInfo * varInfo, iteratorCol *colIter) {
2183
2184 int istatus;
2185 double tscale, tzero;
2186 char temp[80];
2187
2188 switch (bitpix) {
2189 case BYTE_IMG:
2190 case SHORT_IMG:
2191 case LONG_IMG:
2192 istatus = 0;
2193 if (fits_read_key(fptr, TDOUBLE, "BZERO", &tzero, NULL, &istatus))
2194 tzero = 0.0;
2195
2196 istatus = 0;
2197 if (fits_read_key(fptr, TDOUBLE, "BSCALE", &tscale, NULL, &istatus))
2198 tscale = 1.0;
2199
2200 if (tscale == 1.0 && (tzero == 0.0 || tzero == 32768.0 )) {
2201 varInfo->type = LONG;
2202 colIter->datatype = TLONG;
2203 }
2204 else {
2205 varInfo->type = DOUBLE;
2206 colIter->datatype = TDOUBLE;
2207 if (DEBUG_PIXFILTER)
2208 printf("use DOUBLE for %s with BSCALE=%g/BZERO=%g\n",
2209 name, tscale, tzero);
2210 }
2211 break;
2212
2213 case LONGLONG_IMG:
2214 case FLOAT_IMG:
2215 case DOUBLE_IMG:
2216 varInfo->type = DOUBLE;
2217 colIter->datatype = TDOUBLE;
2218 break;
2219 default:
2220 snprintf(temp, 80,"set_image_col_types: unrecognized image bitpix [%d]\n",
2221 bitpix);
2222 ffpmsg(temp);
2223 return gParse.status = PARSE_BAD_TYPE;
2224 }
2225 return 0;
2226 }
2227
2228
2229 /*************************************************************************
2230
2231 Functions used by the evaluator to access FITS data
2232 (find_column, find_keywd, allocateCol, load_column)
2233
2234 *************************************************************************/
2235
find_column(char * colName,void * itslval)2236 static int find_column( char *colName, void *itslval )
2237 {
2238 FFSTYPE *thelval = (FFSTYPE*)itslval;
2239 int col_cnt, status;
2240 int colnum, typecode, type;
2241 long repeat, width;
2242 fitsfile *fptr;
2243 char temp[80];
2244 double tzero,tscale;
2245 int istatus;
2246 DataInfo *varInfo;
2247 iteratorCol *colIter;
2248
2249 if (DEBUG_PIXFILTER)
2250 printf("find_column(%s)\n", colName);
2251
2252 if( *colName == '#' )
2253 return( find_keywd( colName + 1, itslval ) );
2254
2255 fptr = gParse.def_fptr;
2256
2257 status = 0;
2258 col_cnt = gParse.nCols;
2259
2260 if (gParse.hdutype == IMAGE_HDU) {
2261 int i;
2262 if (!gParse.pixFilter) {
2263 gParse.status = COL_NOT_FOUND;
2264 ffpmsg("find_column: IMAGE_HDU but no PixelFilter");
2265 return pERROR;
2266 }
2267
2268 colnum = -1;
2269 for (i = 0; i < gParse.pixFilter->count; ++i) {
2270 if (!fits_strcasecmp(colName, gParse.pixFilter->tag[i]))
2271 colnum = i;
2272 }
2273 if (colnum < 0) {
2274 snprintf(temp, 80, "find_column: PixelFilter tag %s not found", colName);
2275 ffpmsg(temp);
2276 gParse.status = COL_NOT_FOUND;
2277 return pERROR;
2278 }
2279
2280 if( allocateCol( col_cnt, &gParse.status ) ) return pERROR;
2281
2282 varInfo = gParse.varData + col_cnt;
2283 colIter = gParse.colData + col_cnt;
2284
2285 fptr = gParse.pixFilter->ifptr[colnum];
2286 fits_get_img_param(fptr,
2287 MAXDIMS,
2288 &typecode, /* actually bitpix */
2289 &varInfo->naxis,
2290 &varInfo->naxes[0],
2291 &status);
2292 varInfo->nelem = 1;
2293 type = COLUMN;
2294 if (set_image_col_types(fptr, colName, typecode, varInfo, colIter))
2295 return pERROR;
2296 colIter->fptr = fptr;
2297 colIter->iotype = InputCol;
2298 }
2299 else { /* HDU holds a table */
2300 if( gParse.compressed )
2301 colnum = gParse.valCol;
2302 else
2303 if( fits_get_colnum( fptr, CASEINSEN, colName, &colnum, &status ) ) {
2304 if( status == COL_NOT_FOUND ) {
2305 type = find_keywd( colName, itslval );
2306 if( type != pERROR ) ffcmsg();
2307 return( type );
2308 }
2309 gParse.status = status;
2310 return pERROR;
2311 }
2312
2313 if( fits_get_coltype( fptr, colnum, &typecode,
2314 &repeat, &width, &status ) ) {
2315 gParse.status = status;
2316 return pERROR;
2317 }
2318
2319 if( allocateCol( col_cnt, &gParse.status ) ) return pERROR;
2320
2321 varInfo = gParse.varData + col_cnt;
2322 colIter = gParse.colData + col_cnt;
2323
2324 fits_iter_set_by_num( colIter, fptr, colnum, 0, InputCol );
2325 }
2326
2327 /* Make sure we don't overflow variable name array */
2328 strncpy(varInfo->name,colName,MAXVARNAME);
2329 varInfo->name[MAXVARNAME] = '\0';
2330
2331 if (gParse.hdutype != IMAGE_HDU) {
2332 switch( typecode ) {
2333 case TBIT:
2334 varInfo->type = BITSTR;
2335 colIter->datatype = TBYTE;
2336 type = BITCOL;
2337 break;
2338 case TBYTE:
2339 case TSHORT:
2340 case TLONG:
2341 /* The datatype of column with TZERO and TSCALE keywords might be
2342 float or double.
2343 */
2344 snprintf(temp,80,"TZERO%d",colnum);
2345 istatus = 0;
2346 if(fits_read_key(fptr,TDOUBLE,temp,&tzero,NULL,&istatus)) {
2347 tzero = 0.0;
2348 }
2349 snprintf(temp,80,"TSCAL%d",colnum);
2350 istatus = 0;
2351 if(fits_read_key(fptr,TDOUBLE,temp,&tscale,NULL,&istatus)) {
2352 tscale = 1.0;
2353 }
2354 if (tscale == 1.0 && (tzero == 0.0 || tzero == 32768.0 )) {
2355 varInfo->type = LONG;
2356 colIter->datatype = TLONG;
2357 /* Reading an unsigned long column as a long can cause overflow errors.
2358 Treat the column as a double instead.
2359 } else if (tscale == 1.0 && tzero == 2147483648.0 ) {
2360 varInfo->type = LONG;
2361 colIter->datatype = TULONG;
2362 */
2363
2364 }
2365 else {
2366 varInfo->type = DOUBLE;
2367 colIter->datatype = TDOUBLE;
2368 }
2369 type = COLUMN;
2370 break;
2371 /*
2372 For now, treat 8-byte integer columns as type double.
2373 This can lose precision, so the better long term solution
2374 will be to add support for TLONGLONG as a separate datatype.
2375 */
2376 case TLONGLONG:
2377 case TFLOAT:
2378 case TDOUBLE:
2379 varInfo->type = DOUBLE;
2380 colIter->datatype = TDOUBLE;
2381 type = COLUMN;
2382 break;
2383 case TLOGICAL:
2384 varInfo->type = BOOLEAN;
2385 colIter->datatype = TLOGICAL;
2386 type = BCOLUMN;
2387 break;
2388 case TSTRING:
2389 varInfo->type = STRING;
2390 colIter->datatype = TSTRING;
2391 type = SCOLUMN;
2392 if ( width >= MAX_STRLEN ) {
2393 snprintf(temp, 80, "column %d is wider than maximum %d characters",
2394 colnum, MAX_STRLEN-1);
2395 ffpmsg(temp);
2396 gParse.status = PARSE_LRG_VECTOR;
2397 return pERROR;
2398 }
2399 if( gParse.hdutype == ASCII_TBL ) repeat = width;
2400 break;
2401 default:
2402 if (typecode < 0) {
2403 snprintf(temp, 80,"variable-length array columns are not supported. typecode = %d", typecode);
2404 ffpmsg(temp);
2405 }
2406 gParse.status = PARSE_BAD_TYPE;
2407 return pERROR;
2408 }
2409 varInfo->nelem = repeat;
2410 if( repeat>1 && typecode!=TSTRING ) {
2411 if( fits_read_tdim( fptr, colnum, MAXDIMS,
2412 &varInfo->naxis,
2413 &varInfo->naxes[0], &status )
2414 ) {
2415 gParse.status = status;
2416 return pERROR;
2417 }
2418 } else {
2419 varInfo->naxis = 1;
2420 varInfo->naxes[0] = 1;
2421 }
2422 }
2423 gParse.nCols++;
2424 thelval->lng = col_cnt;
2425
2426 return( type );
2427 }
2428
find_keywd(char * keyname,void * itslval)2429 static int find_keywd(char *keyname, void *itslval )
2430 {
2431 FFSTYPE *thelval = (FFSTYPE*)itslval;
2432 int status, type;
2433 char keyvalue[FLEN_VALUE], dtype;
2434 fitsfile *fptr;
2435 double rval;
2436 int bval;
2437 long ival;
2438
2439 status = 0;
2440 fptr = gParse.def_fptr;
2441 if( fits_read_keyword( fptr, keyname, keyvalue, NULL, &status ) ) {
2442 if( status == KEY_NO_EXIST ) {
2443 /* Do this since ffgkey doesn't put an error message on stack */
2444 snprintf(keyvalue,FLEN_VALUE, "ffgkey could not find keyword: %s",keyname);
2445 ffpmsg(keyvalue);
2446 }
2447 gParse.status = status;
2448 return( pERROR );
2449 }
2450
2451 if( fits_get_keytype( keyvalue, &dtype, &status ) ) {
2452 gParse.status = status;
2453 return( pERROR );
2454 }
2455
2456 switch( dtype ) {
2457 case 'C':
2458 fits_read_key_str( fptr, keyname, keyvalue, NULL, &status );
2459 type = STRING;
2460 strcpy( thelval->str , keyvalue );
2461 break;
2462 case 'L':
2463 fits_read_key_log( fptr, keyname, &bval, NULL, &status );
2464 type = BOOLEAN;
2465 thelval->log = bval;
2466 break;
2467 case 'I':
2468 fits_read_key_lng( fptr, keyname, &ival, NULL, &status );
2469 type = LONG;
2470 thelval->lng = ival;
2471 break;
2472 case 'F':
2473 fits_read_key_dbl( fptr, keyname, &rval, NULL, &status );
2474 type = DOUBLE;
2475 thelval->dbl = rval;
2476 break;
2477 default:
2478 type = pERROR;
2479 break;
2480 }
2481
2482 if( status ) {
2483 gParse.status=status;
2484 return pERROR;
2485 }
2486
2487 return( type );
2488 }
2489
allocateCol(int nCol,int * status)2490 static int allocateCol( int nCol, int *status )
2491 {
2492 if( (nCol%25)==0 ) {
2493 if( nCol ) {
2494 gParse.colData = (iteratorCol*) realloc( gParse.colData,
2495 (nCol+25)*sizeof(iteratorCol) );
2496 gParse.varData = (DataInfo *) realloc( gParse.varData,
2497 (nCol+25)*sizeof(DataInfo) );
2498 } else {
2499 gParse.colData = (iteratorCol*) malloc( 25*sizeof(iteratorCol) );
2500 gParse.varData = (DataInfo *) malloc( 25*sizeof(DataInfo) );
2501 }
2502 if( gParse.colData == NULL
2503 || gParse.varData == NULL ) {
2504 if( gParse.colData ) free(gParse.colData);
2505 if( gParse.varData ) free(gParse.varData);
2506 gParse.colData = NULL;
2507 gParse.varData = NULL;
2508 return( *status = MEMORY_ALLOCATION );
2509 }
2510 }
2511 gParse.varData[nCol].data = NULL;
2512 gParse.varData[nCol].undef = NULL;
2513 return 0;
2514 }
2515
load_column(int varNum,long fRow,long nRows,void * data,char * undef)2516 static int load_column( int varNum, long fRow, long nRows,
2517 void *data, char *undef )
2518 {
2519 iteratorCol *var = gParse.colData+varNum;
2520 long nelem,nbytes,row,len,idx;
2521 char **bitStrs, msg[80];
2522 unsigned char *bytes;
2523 int status = 0, anynul;
2524
2525 if (gParse.hdutype == IMAGE_HDU) {
2526 /* This test would need to be on a per varNum basis to support
2527 * cross HDU operations */
2528 fits_read_imgnull(var->fptr, var->datatype, fRow, nRows,
2529 data, undef, &anynul, &status);
2530 if (DEBUG_PIXFILTER)
2531 printf("load_column: IMAGE_HDU fRow=%ld, nRows=%ld => %d\n",
2532 fRow, nRows, status);
2533 } else {
2534
2535 nelem = nRows * var->repeat;
2536
2537 switch( var->datatype ) {
2538 case TBYTE:
2539 nbytes = ((var->repeat+7)/8) * nRows;
2540 bytes = (unsigned char *)malloc( nbytes * sizeof(char) );
2541
2542 ffgcvb(var->fptr, var->colnum, fRow, 1L, nbytes,
2543 0, bytes, &anynul, &status);
2544
2545 nelem = var->repeat;
2546 bitStrs = (char **)data;
2547 for( row=0; row<nRows; row++ ) {
2548 idx = (row)*( (nelem+7)/8 ) + 1;
2549 for(len=0; len<nelem; len++) {
2550 if( bytes[idx] & (1<<(7-len%8)) )
2551 bitStrs[row][len] = '1';
2552 else
2553 bitStrs[row][len] = '0';
2554 if( len%8==7 ) idx++;
2555 }
2556 bitStrs[row][len] = '\0';
2557 }
2558
2559 FREE( (char *)bytes );
2560 break;
2561 case TSTRING:
2562 ffgcfs(var->fptr, var->colnum, fRow, 1L, nRows,
2563 (char **)data, undef, &anynul, &status);
2564 break;
2565 case TLOGICAL:
2566 ffgcfl(var->fptr, var->colnum, fRow, 1L, nelem,
2567 (char *)data, undef, &anynul, &status);
2568 break;
2569 case TLONG:
2570 ffgcfj(var->fptr, var->colnum, fRow, 1L, nelem,
2571 (long *)data, undef, &anynul, &status);
2572 break;
2573 case TDOUBLE:
2574 ffgcfd(var->fptr, var->colnum, fRow, 1L, nelem,
2575 (double *)data, undef, &anynul, &status);
2576 break;
2577 default:
2578 snprintf(msg,80,"load_column: unexpected datatype %d", var->datatype);
2579 ffpmsg(msg);
2580 }
2581 }
2582 if( status ) {
2583 gParse.status = status;
2584 return pERROR;
2585 }
2586
2587 return 0;
2588 }
2589
2590
2591 /*--------------------------------------------------------------------------*/
fits_pixel_filter(PixelFilter * filter,int * status)2592 int fits_pixel_filter (PixelFilter * filter, int * status)
2593 /* Evaluate an expression using the data in the input FITS file(s) */
2594 /*--------------------------------------------------------------------------*/
2595 {
2596 parseInfo Info = { 0 };
2597 int naxis, bitpix;
2598 long nelem, naxes[MAXDIMS];
2599 int col_cnt;
2600 Node *result;
2601 int datatype;
2602 fitsfile * infptr;
2603 fitsfile * outfptr;
2604 char * DEFAULT_TAGS[] = { "X" };
2605 char msg[256];
2606 int writeBlankKwd = 0; /* write BLANK if any output nulls? */
2607
2608 DEBUG_PIXFILTER = getenv("DEBUG_PIXFILTER") ? 1 : 0;
2609
2610 if (*status)
2611 return (*status);
2612
2613 FFLOCK;
2614 if (!filter->tag || !filter->tag[0] || !filter->tag[0][0]) {
2615 filter->tag = DEFAULT_TAGS;
2616 if (DEBUG_PIXFILTER)
2617 printf("using default tag '%s'\n", filter->tag[0]);
2618 }
2619
2620 infptr = filter->ifptr[0];
2621 outfptr = filter->ofptr;
2622 gParse.pixFilter = filter;
2623
2624 if (ffiprs(infptr, 0, filter->expression, MAXDIMS,
2625 &Info.datatype, &nelem, &naxis, naxes, status)) {
2626 goto CLEANUP;
2627 }
2628
2629 if (nelem < 0) {
2630 nelem = -nelem;
2631 }
2632
2633 {
2634 /* validate result type */
2635 const char * type = 0;
2636 switch (Info.datatype) {
2637 case TLOGICAL: type = "LOGICAL"; break;
2638 case TLONG: type = "LONG"; break;
2639 case TDOUBLE: type = "DOUBLE"; break;
2640 case TSTRING: type = "STRING";
2641 *status = pERROR;
2642 ffpmsg("pixel_filter: cannot have string image");
2643 case TBIT: type = "BIT";
2644 if (DEBUG_PIXFILTER)
2645 printf("hmm, image from bits?\n");
2646 break;
2647 default: type = "UNKNOWN?!";
2648 *status = pERROR;
2649 ffpmsg("pixel_filter: unexpected result datatype");
2650 }
2651 if (DEBUG_PIXFILTER)
2652 printf("result type is %s [%d]\n", type, Info.datatype);
2653 if (*status)
2654 goto CLEANUP;
2655 }
2656
2657 if (fits_get_img_param(infptr, MAXDIMS,
2658 &bitpix, &naxis, &naxes[0], status)) {
2659 ffpmsg("pixel_filter: unable to read input image parameters");
2660 goto CLEANUP;
2661 }
2662
2663 if (DEBUG_PIXFILTER)
2664 printf("input bitpix %d\n", bitpix);
2665
2666 if (Info.datatype == TDOUBLE) {
2667 /* for floating point expressions, set the default output image to
2668 bitpix = -32 (float) unless the default is already a double */
2669 if (bitpix != DOUBLE_IMG)
2670 bitpix = FLOAT_IMG;
2671 }
2672
2673 /* override output image bitpix if specified by caller */
2674 if (filter->bitpix)
2675 bitpix = filter->bitpix;
2676 if (DEBUG_PIXFILTER)
2677 printf("output bitpix %d\n", bitpix);
2678
2679 if (fits_create_img(outfptr, bitpix, naxis, naxes, status)) {
2680 ffpmsg("pixel_filter: unable to create output image");
2681 goto CLEANUP;
2682 }
2683
2684 /* transfer keycards */
2685 {
2686 int i, ncards, more;
2687 if (fits_get_hdrspace(infptr, &ncards, &more, status)) {
2688 ffpmsg("pixel_filter: unable to determine number of keycards");
2689 goto CLEANUP;
2690 }
2691
2692 for (i = 1; i <= ncards; ++i) {
2693
2694 int keyclass;
2695 char card[FLEN_CARD];
2696
2697 if (fits_read_record(infptr, i, card, status)) {
2698 snprintf(msg, 256,"pixel_filter: unable to read keycard %d", i);
2699 ffpmsg(msg);
2700 goto CLEANUP;
2701 }
2702
2703 keyclass = fits_get_keyclass(card);
2704 if (keyclass == TYP_STRUC_KEY) {
2705 /* output structure defined by fits_create_img */
2706 }
2707 else if (keyclass == TYP_COMM_KEY && i < 12) {
2708 /* assume this is one of the FITS standard comments */
2709 }
2710 else if (keyclass == TYP_NULL_KEY && bitpix < 0) {
2711 /* do not transfer BLANK to real output image */
2712 }
2713 else if (keyclass == TYP_SCAL_KEY && bitpix < 0) {
2714 /* do not transfer BZERO, BSCALE to real output image */
2715 }
2716 else if (fits_write_record(outfptr, card, status)) {
2717 snprintf(msg,256, "pixel_filter: unable to write keycard '%s' [%d]\n",
2718 card, *status);
2719 ffpmsg(msg);
2720 goto CLEANUP;
2721 }
2722 }
2723 }
2724
2725 switch (bitpix) {
2726 case BYTE_IMG: datatype = TLONG; Info.datatype = TBYTE; break;
2727 case SHORT_IMG: datatype = TLONG; Info.datatype = TSHORT; break;
2728 case LONG_IMG: datatype = TLONG; Info.datatype = TLONG; break;
2729 case FLOAT_IMG: datatype = TDOUBLE; Info.datatype = TFLOAT; break;
2730 case DOUBLE_IMG: datatype = TDOUBLE; Info.datatype = TDOUBLE; break;
2731
2732 default:
2733 snprintf(msg, 256,"pixel_filter: unexpected output bitpix %d\n", bitpix);
2734 ffpmsg(msg);
2735 *status = pERROR;
2736 goto CLEANUP;
2737 }
2738
2739 if (bitpix > 0) { /* arrange for NULLs in output */
2740 long nullVal = filter->blank;
2741 if (!filter->blank) {
2742 int tstatus = 0;
2743 if (fits_read_key_lng(infptr, "BLANK", &nullVal, 0, &tstatus)) {
2744
2745 writeBlankKwd = 1;
2746
2747 if (bitpix == BYTE_IMG)
2748 nullVal = UCHAR_MAX;
2749 else if (bitpix == SHORT_IMG)
2750 nullVal = SHRT_MIN;
2751 else if (bitpix == LONG_IMG)
2752 nullVal = LONG_MIN;
2753 else
2754 printf("unhandled positive output BITPIX %d\n", bitpix);
2755 }
2756
2757 filter->blank = nullVal;
2758 }
2759
2760 fits_set_imgnull(outfptr, filter->blank, status);
2761 if (DEBUG_PIXFILTER)
2762 printf("using blank %ld\n", nullVal);
2763
2764 }
2765
2766 if (!filter->keyword[0]) {
2767 iteratorCol * colIter;
2768 DataInfo * varInfo;
2769
2770 /*************************************/
2771 /* Create new iterator Output Column */
2772 /*************************************/
2773 col_cnt = gParse.nCols;
2774 if (allocateCol(col_cnt, status))
2775 goto CLEANUP;
2776 gParse.nCols++;
2777
2778 colIter = &gParse.colData[col_cnt];
2779 colIter->fptr = filter->ofptr;
2780 colIter->iotype = OutputCol;
2781 varInfo = &gParse.varData[col_cnt];
2782 set_image_col_types(colIter->fptr, "CREATED", bitpix, varInfo, colIter);
2783
2784 Info.maxRows = -1;
2785
2786 if (ffiter(gParse.nCols, gParse.colData, 0,
2787 0, parse_data, &Info, status) == -1)
2788 *status = 0;
2789 else if (*status)
2790 goto CLEANUP;
2791
2792 if (Info.anyNull) {
2793 if (writeBlankKwd) {
2794 fits_update_key_lng(outfptr, "BLANK", filter->blank, "NULL pixel value", status);
2795 if (*status)
2796 ffpmsg("pixel_filter: unable to write BLANK keyword");
2797 if (DEBUG_PIXFILTER) {
2798 printf("output has NULLs\n");
2799 printf("wrote blank [%d]\n", *status);
2800 }
2801 }
2802 }
2803 else if (bitpix > 0) /* never used a null */
2804 if (fits_set_imgnull(outfptr, -1234554321, status))
2805 ffpmsg("pixel_filter: unable to reset imgnull");
2806 }
2807 else {
2808
2809 /* Put constant result into keyword */
2810 char * parName = filter->keyword;
2811 char * parInfo = filter->comment;
2812
2813 result = gParse.Nodes + gParse.resultNode;
2814 switch (Info.datatype) {
2815 case TDOUBLE:
2816 ffukyd(outfptr, parName, result->value.data.dbl, 15, parInfo, status);
2817 break;
2818 case TLONG:
2819 ffukyj(outfptr, parName, result->value.data.lng, parInfo, status);
2820 break;
2821 case TLOGICAL:
2822 ffukyl(outfptr, parName, result->value.data.log, parInfo, status);
2823 break;
2824 case TBIT:
2825 case TSTRING:
2826 ffukys(outfptr, parName, result->value.data.str, parInfo, status);
2827 break;
2828 default:
2829 snprintf(msg, 256,"pixel_filter: unexpected constant result type [%d]\n",
2830 Info.datatype);
2831 ffpmsg(msg);
2832 }
2833 }
2834
2835 CLEANUP:
2836 ffcprs();
2837 FFUNLOCK;
2838 return (*status);
2839 }
2840