summaryrefslogtreecommitdiffstats
path: root/funtools/funim.c
diff options
context:
space:
mode:
Diffstat (limited to 'funtools/funim.c')
-rw-r--r--funtools/funim.c1680
1 files changed, 1680 insertions, 0 deletions
diff --git a/funtools/funim.c b/funtools/funim.c
new file mode 100644
index 0000000..f15e4f3
--- /dev/null
+++ b/funtools/funim.c
@@ -0,0 +1,1680 @@
+/*
+ * Copyright (c) 1999-2003 Smithsonian Astrophysical Observatory
+ */
+
+#include <funtoolsP.h>
+
+/*
+ *
+ * private routines
+ *
+ */
+
+#define MAXEV 10240
+
+/*
+ *
+ * _FunTableBin -- bin a blocked section of data from a table
+ *
+*/
+#ifdef ANSI_FUNC
+static void *
+_FunTableBin(Fun fun, void *buf, char *plist)
+#else
+static void *_FunTableBin(fun, buf, plist)
+ Fun fun;
+ void *buf;
+ char *plist;
+#endif
+{
+ int skip; /* temp for skipping data to end of ext */
+ int offset[2]; /* offset into event record for x, y */
+ int size[2]; /* size of x, y value in bytes */
+ int val[2]; /* x, y bin values of this event */
+ int convert; /* whether we have to convert to native */
+ int index; /* index into image array for this pixel */
+ int esize; /* max number of events to read at once */
+ int dofilt; /* true if we can filter */
+ int dovcol; /* true if we use a value column */
+ int voffset; /* offset into event record for vcol */
+ /* not used, so avoid gcc compiler warning */
+ /* int vsize; */ /* size of vcol in bytes */
+ int vtype; /* data type of event vcol */
+ int imdim; /* total dims of image buffer */
+ int itsnan; /* if value is a nan */
+ int transparent; /* allow everything through the filter? */
+ int *rbuf=NULL; /* valid event flags */
+ size_t i, j; /* loop counters */
+ size_t get; /* number of events to read */
+ size_t got; /* number of events read */
+ unsigned char *ebuf=NULL; /* event buffer */
+ unsigned char *eptr; /* current pointer into ebuf */
+ double block2; /* block * block */
+ double minval[2]; /* x, y min axis values */
+ double binsiz[2]; /* x, y bin values */
+ char *obuf; /* output image buffer */
+ char *pbuf; /* temp plist buf */
+ char tbuf[SZ_LINE]; /* ever-present temp buffer */
+ short sval; /* short val */
+ unsigned short usval; /* unsigned short val */
+ int ival; /* int val */
+ longlong lval; /* 64-bit int val */
+ unsigned int uival; /* unsigned int vala */
+ float fval; /* float val */
+ double dval; /* double val */
+ double vval; /* value of the event */
+
+ /* just in case ... */
+ if( !_FunValid(fun) ){
+ gerror(stderr, "invalid funtools handle\n");
+ return NULL;
+ }
+
+ /* if we delayed the open before, we have to open now */
+ if( !fun->header && strchr(fun->mode, 'r') )
+ _FunFITSOpen(fun, fun->fname, "r");
+
+ /* make sure something good can happen */
+ if( !fun->header || !fun->header->table ){
+ gerror(stderr, "not a binary table\n");
+ return NULL;
+ }
+
+ /* make sure we have binning parameters */
+ if( (fun->bin[0] < 0) || ((fun->dims > 1) && (fun->bin[1] < 0)) ){
+ gerror(stderr, "invalid or missing binning column(s) for table\n");
+ return NULL;
+ }
+
+ /* make sure we have valid output dimensions */
+ if( (fun->odim1 <= 0) || ((fun->dims > 1) && (fun->odim2 <= 0)) ){
+ gerror(stderr, "invalid or zero image dimension(s) for table (%s;%d,%d)\n",
+ fun->bincols ? fun->bincols : "invalid bincols?",
+ fun->odim1, fun->odim2);
+ return NULL;
+ }
+
+ /* we have to convert to native if the data is not the
+ same as the big-endian-ness of the machine */
+ convert = (fun->endian != is_bigendian());
+
+ /* set x and y info -- whatever columns they actually are */
+ offset[0] = fun->header->table->col[fun->bin[0]].offset;
+ size[0] = fun->header->table->col[fun->bin[0]].size;
+ minval[0] = fun->min1;
+ binsiz[0] = fun->binsiz1;
+ if( fun->dims > 1 ){
+ offset[1] = fun->header->table->col[fun->bin[1]].offset;
+ size[1] = fun->header->table->col[fun->bin[1]].size;
+ minval[1] = fun->min2;
+ binsiz[1] = fun->binsiz2;
+ }
+ /* optional value column for binning */
+ if( fun->vbin >= 0 ){
+ dovcol = 1;
+ voffset = fun->header->table->col[fun->vbin].offset;
+ /* not used, so avoid gcc compiler warning */
+ /* vsize = fun->header->table->col[fun->vbin].size; */
+ vtype = fun->header->table->col[fun->vbin].type;
+ }
+ else{
+ dovcol = 0;
+ }
+
+ /* init some convenience variables */
+ imdim = fun->odim1*MAX(fun->odim2,1);
+ block2 = (double)fun->block*fun->block;
+ transparent = 0;
+
+ /* see if default data type is to be over-ridden */
+ pbuf = xstrdup(plist);
+ if( _FunKeyword(pbuf, "bitpix", NULL, tbuf, SZ_LINE) ){
+ if( (i = atoi(tbuf)) ){
+ fun->obitpix = i;
+ fun->odtype = fun->obitpix / FT_WORDLEN;
+ }
+ }
+ else if( _FunKeyword(pbuf, "mask", "FUN_MASK", tbuf, SZ_LINE) ){
+ if( !strcasecmp(tbuf, "transparent") )
+ transparent = 1;
+ }
+ if( pbuf ) xfree(pbuf);
+
+ /* allocate output buffer, if necessary */
+ if( buf )
+ obuf = buf;
+ else{
+ obuf = (char *)xcalloc(imdim, ABS(fun->odtype));
+ if( !obuf ){
+ gerror(stderr, "can't allocate image buffer (%dx%dx%d)\n",
+ fun->odim1, MAX(fun->odim2,1), ABS(fun->odtype));
+ return NULL;
+ }
+ }
+
+ /* create the plist string */
+ if( convert )
+ strcpy(tbuf, "convert=true");
+ else
+ strcpy(tbuf, "convert=false");
+ /* add the binning key */
+ if( fun->bincols ){
+ strcat(tbuf, ",");
+ strcat(tbuf, fun->bincols);
+ }
+ /* add the index flag */
+ if( fun->idx == 1 ){
+ strcat(tbuf, ",idx=true");
+ }
+ else if( fun->idx == -1 ){
+ strcat(tbuf, ",idx=false");
+ }
+
+ /* open the filter, if its not already been done */
+ if( !fun->filt ){
+ if( fun->filter && *fun->filter ){
+ fun->filt = FilterOpen(fun->header, fun->filter, tbuf);
+ /* check for error */
+ if( fun->filt == NULL ){
+ fun->filt = NOFILTER;
+ gerror(stderr, "can't open filter: %s\n", fun->filter);
+ return NULL;
+ }
+ }
+ else{
+ fun->filt = NOFILTER;
+ fun->dofilt = 0;
+ }
+ }
+
+ /* allocate space for a pile of event records */
+ for(esize=MIN(MAXEV,(fun->left>0?fun->left:MAXEV)); esize>0; esize /=2){
+ if( ((ebuf = (unsigned char *)malloc(fun->rawsize*esize)) != NULL) &&
+ ((rbuf = (int *)malloc(sizeof(int)*esize)) != NULL) )
+ break;
+ }
+
+ /* while there are still events to read ... */
+ while( fun->left != 0 ){
+ /* figure out how many to read this time */
+ if( fun->left > 0 )
+ get = MIN(esize, fun->left);
+ /* read to EOF */
+ else
+ get = esize;
+ /* read in a pile of events */
+ _FunRead(fun, (void *)ebuf, fun->rawsize, get, &got);
+ /* if we got what we asked for ... */
+ if( got == get ){
+ /* if we are reading a set number of rows ... */
+ if( fun->left > 0 ){
+ /* subtract the number we just got */
+ fun->left -= got;
+ /* just in case! */
+ if( fun->left < 0 )
+ fun->left = 0;
+ }
+ }
+ /* didn't get what we asked for ... */
+ else{
+ /* if we are reading a set number of rows ... */
+ if( fun->left > 0 ){
+ /* using an index => already filtering, so fewer rows are acceptable */
+ if( (fun->filt == NOFILTER) || (fun->filt->doidx != 1) ){
+ /* otherwise its unexpected */
+ gerror(stderr, "unexpected EOF reading rows\n");
+ }
+ /* ensure we stop after this iteration */
+ fun->left = 0;
+ }
+ /* if we are reading til EOF, we just got there */
+ else{
+ fun->left = 0;
+ }
+ }
+ /* update how many input events we have read */
+ fun->io += got;
+ fun->bytes += (got*fun->rawsize);
+
+ /* filter these events through the co-process */
+ if( fun->dofilt )
+ dofilt = FilterEvents(fun->filt, (void *)ebuf, fun->rawsize, (int)got,
+ (void *)rbuf);
+ else
+ dofilt = 0;
+
+ /* process each event in this batch */
+ for(eptr=ebuf, i=0; i<got; i++, eptr += fun->rawsize){
+ /* if its not a valid event, skip it */
+ if( dofilt && !transparent && !rbuf[i] ){
+ continue;
+ }
+
+ /* get y and x values and convert to int data type for indexing */
+ for(j=0; j<(size_t)fun->dims; j++){
+ switch(fun->header->table->col[fun->bin[j]].type){
+ /* unsigned 8-bit int */
+ case 'B':
+ val[j] = itlp2i((double)eptr[offset[j]], minval[j], binsiz[j], 'B');
+ break;
+ /* 16-bit short */
+ case 'I':
+ ColumnLoad((void *)(eptr+offset[j]), size[j], 1, convert, &sval);
+ val[j] = itlp2i((double)sval, minval[j], binsiz[j], 'I');
+ break;
+ /* unsigned 16-bit int */
+ case 'U':
+ ColumnLoad((void *)(eptr+offset[j]), size[j], 1, convert, &usval);
+ val[j] = itlp2i((double)usval, minval[j], binsiz[j], 'U');
+ break;
+ /* 32-bit int */
+ case 'J':
+ ColumnLoad((void *)(eptr+offset[j]), size[j], 1, convert, &ival);
+ val[j] = itlp2i((double)ival, minval[j], binsiz[j], 'J');
+ break;
+ /* unsigned 32-bit int */
+ case 'K':
+ ColumnLoad((void *)(eptr+offset[j]), size[j], 1, convert, &lval);
+ val[j] = itlp2i((double)ival, minval[j], binsiz[j], 'K');
+ break;
+ /* unsigned 32-bit int */
+ case 'V':
+ ColumnLoad((void *)(eptr+offset[j]), size[j], 1, convert, &uival);
+ val[j] = itlp2i((double)uival, minval[j], binsiz[j], 'V');
+ break;
+ /* 32-bit float */
+ case 'E':
+ ColumnLoad((void *)(eptr+offset[j]), size[j], 1, convert, &fval);
+ val[j] = itlp2i((double)fval, minval[j], binsiz[j], 'E');
+ break;
+ /* 64-bit float */
+ case 'D':
+ ColumnLoad((void *)(eptr+offset[j]), size[j], 1, convert, &dval);
+ val[j] = itlp2i(dval, minval[j], binsiz[j], 'D');
+ break;
+ case 'X':
+ case 'L':
+ case 'C':
+ case 'A':
+ case 'M':
+ case 'P':
+ default:
+ goto done;
+ }
+ }
+
+ /* make sure we are in the desired section, and
+ calculate the offset into the image array */
+ if( fun->dims == 1 ){
+ if( (val[0] < fun->x0) || (val[0] > fun->x1) )
+ continue;
+ index = (val[0] - fun->x0)/fun->block;
+ }
+ else{
+ if( (val[0] < fun->x0) || (val[0] > fun->x1) ||
+ (val[1] < fun->y0) || (val[1] > fun->y1) )
+ continue;
+ index = ((val[1] - fun->y0)/fun->block)*fun->odim1 +
+ (val[0] - fun->x0)/fun->block;
+ }
+
+ /* make sure its valid */
+ if( (index<0) || (index>=imdim) )
+ continue;
+
+ /* use value column or 1 as event value */
+ itsnan = 0;
+ if( dovcol ){
+ /* if vcol is not double, we must do full type convert */
+ if( vtype != 'D' ){
+ ft_acht2('D', &vval, vtype, eptr+voffset, 1, convert, 0);
+ }
+ /* for same types, we might have to swap the bytes of the double */
+ else if( convert ){
+ swap8((char *)&vval, (char *)(eptr+voffset), 8);
+ }
+ /* otherwise its an easy copy */
+ else{
+ memcpy(&vval, eptr+voffset, ft_sizeof(vtype));
+ }
+ /* value operators */
+ switch(fun->vop){
+ case '/':
+ if( vval ){
+ vval = 1.0/vval;
+ }
+ else{
+ itsnan = 1;
+ }
+ break;
+ default:
+ break;
+ }
+ if( vval == getnand() ){
+ itsnan = 1;
+ vval = 0.0;
+ }
+ }
+ else{
+ vval = 1.0;
+ }
+
+ /* add this event to the image array */
+ switch(fun->odtype){
+ case TY_CHAR:
+ if( !dovcol ){
+ if( ((unsigned char *)obuf)[index] >
+ (((unsigned char *)obuf)[index]+(unsigned char)vval) ){
+ fun->overflow++;
+ }
+ else{
+ ((unsigned char *)obuf)[index] += (unsigned char)vval;
+ }
+ }
+ else{
+ if(fun->doblank && (((unsigned char *)obuf)[index] == fun->blank)){
+ ;
+ }
+ else if( itsnan ){
+ ((unsigned char *)obuf)[index] = fun->blank;
+ }
+ else{
+ ((unsigned char *)obuf)[index] += (unsigned char)vval;
+ }
+ }
+ break;
+ /* 32-bit int */
+ case TY_USHORT:
+ if( !dovcol ){
+ if( ((unsigned short *)obuf)[index] >
+ (((unsigned short *)obuf)[index]+(unsigned short)vval) ){
+ fun->overflow++;
+ }
+ else{
+ ((unsigned short *)obuf)[index] += (unsigned short)vval;
+ }
+ }
+ else{
+ if(fun->doblank && (((unsigned short *)obuf)[index] == fun->blank)){
+ ;
+ }
+ else if( itsnan ){
+ ((unsigned short *)buf)[index] = fun->blank;
+ }
+ else{
+ ((unsigned short *)buf)[index] += (unsigned short)vval;
+ }
+ }
+ break;
+ /* 32-bit int */
+ case TY_SHORT:
+ if( !dovcol ){
+ if( ((short *)obuf)[index] > (((short *)obuf)[index]+(short)vval) ){
+ fun->overflow++;
+ }
+ else{
+ ((short *)obuf)[index] += (short)vval;
+ }
+ }
+ else{
+ if( fun->doblank && (((short *)obuf)[index] == fun->blank) ){
+ ;
+ }
+ else if( itsnan ){
+ ((short *)obuf)[index] = fun->blank;
+ }
+ else{
+ ((short *)obuf)[index] += (short)vval;
+ }
+ }
+ break;
+ /* 32-bit int */
+ case TY_INT:
+ if( !dovcol ){
+ if( ((int *)obuf)[index] > (((int *)obuf)[index]+(int)vval) ){
+ fun->overflow++;
+ }
+ else{
+ ((int *)obuf)[index] += (int)vval;
+ }
+ }
+ else{
+ if( fun->doblank && (((int *)obuf)[index] == fun->blank) ){
+ ;
+ }
+ else if( itsnan ){
+ ((int *)obuf)[index] = fun->blank;
+ }
+ else{
+ ((int *)obuf)[index] += (int)vval;
+ }
+ }
+ break;
+ /* 64-bit int */
+ case TY_LONG:
+#if HAVE_LONG_LONG == 0
+ gerror(stderr,
+ "64-bit data support not built (long long not available)\n");
+#endif
+ if( !dovcol ){
+ if(((longlong *)obuf)[index] >
+ (((longlong *)obuf)[index]+(longlong)vval)){
+ fun->overflow++;
+ }
+ else{
+ ((longlong *)obuf)[index] += (longlong)vval;
+ }
+ }
+ else{
+ if( fun->doblank && (((longlong *)obuf)[index] == fun->blank) ){
+ ;
+ }
+ else if( itsnan ){
+ ((longlong *)obuf)[index] = fun->blank;
+ }
+ else{
+ ((longlong *)obuf)[index] += (longlong)vval;
+ }
+ }
+ break;
+ /* 32-bit float */
+ case TY_FLOAT:
+ if( !dovcol ){
+ ((float *)obuf)[index] += (float)vval;
+ }
+ else{
+ if( isnanf(((float *)obuf)[index]) ){
+ ;
+ }
+ else if( itsnan ){
+ ((float *)obuf)[index] = getnanf();
+ }
+ else{
+ ((float *)obuf)[index] += (float)vval;
+ }
+ }
+ break;
+ /* 64-bit float */
+ case TY_DOUBLE:
+ if( !dovcol ){
+ ((double *)obuf)[index] += (double)vval;
+ }
+ else{
+ if( isnand(((double *)obuf)[index]) ){
+ ;
+ }
+ else if( itsnan ){
+ ((double *)obuf)[index] = getnand();
+ }
+ else{
+ ((double *)obuf)[index] += (double)vval;
+ }
+ }
+ break;
+ default:
+ goto done;
+ }
+ }
+ }
+
+ /* binned all events, average each pixel, if necessary */
+ if( fun->btype == FUN_AVG ){
+ for(i=0; i<(size_t)imdim; i++){
+ switch(fun->odtype){
+ case TY_CHAR:
+ ((unsigned char *)obuf)[i] /= block2;
+ break;
+ case TY_USHORT:
+ ((unsigned short *)obuf)[i] /= block2;
+ break;
+ case TY_SHORT:
+ ((short *)obuf)[i] /= block2;
+ break;
+ case TY_INT:
+ ((int *)obuf)[i] /= block2;
+ break;
+ case TY_LONG:
+#if HAVE_LONG_LONG == 0
+ gerror(stderr,
+ "64-bit data support not built (long long not available)\n");
+#endif
+ ((longlong *)obuf)[i] /= block2;
+ break;
+ case TY_FLOAT:
+ ((float *)obuf)[i] /= block2;
+ break;
+ case TY_DOUBLE:
+ ((double *)obuf)[i] /= block2;
+ break;
+ default:
+ break;
+ }
+ }
+ }
+
+done:
+ /* free up space */
+ if( ebuf ) xfree(ebuf);
+ if( rbuf ) xfree(rbuf);
+
+ /* if we have no more rows left, clean up */
+ if( !fun->left ){
+ /* skip to end of extension */
+ if( fun->bytes ){
+ /* skip any unread events */
+ skip = (fun->rawsize*fun->total) - fun->bytes;
+ gskip(fun->gio, (off_t)skip);
+ fun->bytes += skip;
+ /* skip padding */
+ if( (skip = FT_BLOCK - (fun->bytes % FT_BLOCK)) == FT_BLOCK )
+ skip = 0;
+ gskip(fun->gio, (off_t)skip);
+ /* reset io flag */
+ /* fun->io = 0; */
+ }
+ }
+
+ /* close filter if done with the filter process */
+ if( !fun->left && (fun->filt != NOFILTER) ){
+ FilterClose(fun->filt);
+ fun->filt = NULL;
+ }
+
+ /* that's good news */
+ return (void *)obuf;
+}
+
+/*
+ *
+ * _FunImageExtract -- grab a blocked section of data from an array
+ *
+*/
+#ifdef ANSI_FUNC
+static void *
+_FunImageExtract(Fun fun, void *buf, int rstart, int rstop, char *plist)
+#else
+static void *_FunImageExtract(fun, buf, rstart, rstop, plist)
+ Fun fun;
+ void *buf;
+ int rstart;
+ int rstop;
+ char *plist;
+#endif
+{
+ int i, j, k;
+ int ox, oy;
+ int convert;
+ int yoff;
+ int doscale=0;
+ size_t got;
+ off_t newpos;
+ char *obuf;
+ char *pbuf;
+ char **bufs;
+ char tbuf[SZ_LINE];
+ unsigned char cval;
+ unsigned short usval;
+ short sval;
+ int ival;
+ longlong lval;
+ float fval;
+ double dval;
+ double block2;
+ register double val;
+ longlong vall;
+
+ /* make sure we have something */
+ if( !_FunValid(fun) )
+ return NULL;
+
+ /* if we delayed the open before, we have to open now */
+ if( !fun->header && strchr(fun->mode, 'r') )
+ _FunFITSOpen(fun, fun->fname, "r");
+
+ /* just in case ... */
+ if( !fun->header || !fun->header->image ){
+ gerror(stderr, "data is not a table or an image\n");
+ return NULL;
+ }
+
+ /* make sure we have valid dimensions */
+ if( fun->odim1 <=0 ){
+ gerror(stderr, "invalid dim1 (%d) from x0,x1=%d,%d\n",
+ fun->odim1, fun->x0, fun->x1);
+ return NULL;
+ }
+ if( (fun->dims > 1) && (fun->odim2 <=0) ){
+ gerror(stderr, "invalid dim2 (%d) from y0,y1=%d,%d\n",
+ fun->odim2, fun->y0, fun->y1);
+ return NULL;
+ }
+
+ /* we have to convert to native if the data is not the
+ same as the big-endian-ness of the machine */
+ convert = (fun->endian != is_bigendian());
+
+ /* set up some oft-used variables */
+ block2 = (double)fun->block*fun->block;
+
+ /* we cannot handle BLANK==0 */
+ if( (fun->dtype>=-2) && fun->doblank && (fun->blank==0) ){
+ gwarning(stderr, "BLANK==0 detected ... disabling NaN checking\n");
+ fun->doblank = 0;
+ }
+
+ /* see if default data type is to be over-ridden */
+ pbuf = xstrdup(plist);
+ if( _FunKeyword(pbuf, "bitpix", NULL, tbuf, SZ_LINE) ){
+ if( (i = atoi(tbuf)) ){
+ fun->obitpix = i;
+ fun->odtype = fun->obitpix / FT_WORDLEN;
+ }
+ }
+ /* see if scaling is explicitly turned off */
+ if( _FunKeyword(pbuf, "bscale", NULL, tbuf, SZ_LINE) ){
+ if( isfalse(tbuf) )
+ doscale = -1;
+ }
+ if( pbuf ) xfree(pbuf);
+
+ /* actually, we only scale if we have scaled int data */
+ if( (doscale==0) && (fun->scaled&FUN_SCALE_EXISTS) && (fun->dtype>=-2) ){
+ /* flag that we scaled the input data */
+ fun->scaled |= FUN_SCALE_APPLIED;
+ doscale = 1;
+ }
+ else
+ doscale = 0;
+
+ /* make sure we are not out of bounds */
+ if( rstart < 1 ) rstart = 1;
+ if( rstop < 1 ) rstop = 1;
+ if( rstop > MAX(fun->odim2,1) ) rstop = MAX(fun->odim2,1);
+
+ /* allocate output buffer, if necessary */
+ if( buf )
+ obuf = buf;
+ else{
+ obuf = (char *)xcalloc(fun->odim1*(rstop-rstart+1), ABS(fun->odtype));
+ if( !obuf ){
+ gerror(stderr, "can't allocate image buffer (%dx%dx%d)\n",
+ fun->odim1, (rstop-rstart+1), ABS(fun->odtype));
+ return NULL;
+ }
+ }
+
+ /* allocate some row buffers: number of bufs is equal to the block factor */
+ bufs = (char **)xmalloc(fun->block * sizeof(char *));
+ for(i=0; i<fun->block; i++){
+ bufs[i] = (char *)xmalloc((fun->x1-fun->x0+1)*ABS(fun->dtype));
+ }
+
+ /* for all rows in the output image (or image section) ... */
+ for(oy=rstart; oy<=rstop; oy++){
+ yoff = (oy-rstart) * fun->odim1;
+ /* get input rows associated with this output row (by block factor) */
+ for(j=0; j<fun->block; i++, j++){
+ /* calculate specific y line we are processing */
+ i = (fun->y0 - 1) + ((oy-1) * fun->block) + j;
+ /* calculate the new position into the data */
+ newpos = ((i*fun->dim1)+(fun->x0-1)) * ABS(fun->dtype);
+ /* skip to new position */
+ if( gskip(fun->gio, newpos - fun->curpos) < 0 ){
+ gerror(stderr,
+ "_FunImageExtract can't seek %d bytes into image data\n",
+ newpos - fun->curpos);
+ if( !buf && obuf ) xfree(obuf);
+ return(NULL);
+ }
+ /* read the row */
+ _gread(fun->gio, bufs[j], ABS(fun->dtype), (fun->x1-fun->x0+1), &got);
+ /* check for EOF */
+ if( got != (size_t)(fun->x1-fun->x0+1) ){
+ gerror(stderr, "unexpected EOF reading data\n");
+ if( !buf && obuf ) xfree(obuf);
+ return(NULL);
+ }
+ /* convert to native, if necessary */
+ if( convert )
+ swap_data(bufs[j], got, fun->dtype);
+ /* update the current position */
+ fun->curpos = newpos + (got * ABS(fun->dtype));
+ }
+
+ /* process all pixels in this output row */
+ for(ox=0; ox<fun->odim1; ox++){
+ /* reset temp value */
+ val = 0.0;
+ vall = 0;
+ /* for all rows in the input blocks */
+ for(i=0; i<fun->block; i++){
+ /* sum the input pixels which make up this output pixel */
+ for(j=ox*fun->block, k=0; k<fun->block; j++, k++){
+ switch(fun->dtype){
+ case TY_CHAR:
+ cval = ((unsigned char **)bufs)[i][j];
+ if( fun->doblank && (cval == fun->blank) )
+ goto blank;
+ if( doscale )
+ val += ((double)cval*fun->bscale)+fun->bzero;
+ else
+ val += (double)cval;
+ break;
+ case TY_USHORT:
+ usval = ((unsigned short **)bufs)[i][j];
+ if( fun->doblank && (usval == fun->blank) )
+ goto blank;
+ if( doscale )
+ val += ((double)usval*fun->bscale)+fun->bzero;
+ else
+ val += (double)usval;
+ break;
+ case TY_SHORT:
+ sval = ((short **)bufs)[i][j];
+ if( fun->doblank && (sval == fun->blank) )
+ goto blank;
+ if( doscale )
+ val += ((double)sval*fun->bscale)+fun->bzero;
+ else
+ val += (double)sval;
+ break;
+ case TY_INT:
+ ival = ((int **)bufs)[i][j];
+ if( fun->doblank && (ival == fun->blank) )
+ goto blank;
+ if( doscale )
+ val += ((double)ival*fun->bscale)+fun->bzero;
+ else
+ val += (double)ival;
+ break;
+ case TY_LONG:
+#if HAVE_LONG_LONG == 0
+ gerror(stderr,
+ "64-bit data support not built (long long not available)\n");
+#endif
+ /* argh ... 64-bit does not fit into a double! */
+ lval = ((longlong **)bufs)[i][j];
+ if( fun->doblank && (lval == fun->blank) )
+ goto blank;
+ /* probably going to lose precision here */
+ if( doscale )
+ vall += ((double)lval*fun->bscale)+fun->bzero;
+ else
+ vall += lval;
+ break;
+ case TY_FLOAT:
+ fval = ((float **)bufs)[i][j];
+ if( isnanf(fval) )
+ goto blank;
+ val += (double)fval;
+ break;
+ case TY_DOUBLE:
+ dval = ((double **)bufs)[i][j];
+ if( isnand(dval) )
+ goto blank;
+ val += dval;
+ break;
+ default:
+ break;
+ }
+ }
+ } /* x block loop */
+
+ /* take average of this output pixel, if necessary */
+ if( fun->btype == FUN_AVG ){
+ val /= block2;
+ vall /= block2;
+ }
+
+ /* store the resulting value in the output image */
+ /* argh ... 64-bit does not fit into a double! */
+ if( fun->dtype == TY_LONG ){
+ switch(fun->odtype){
+ case TY_CHAR:
+ ((unsigned char *)obuf)[yoff+ox] = (unsigned char)vall;
+ break;
+ case TY_USHORT:
+ ((unsigned short *)obuf)[yoff+ox] = (unsigned short)vall;
+ break;
+ case TY_SHORT:
+ ((short *)obuf)[yoff+ox] = (short)vall;
+ break;
+ case TY_INT:
+ ((int *)obuf)[yoff+ox] = (int)vall;
+ break;
+ case TY_LONG:
+#if HAVE_LONG_LONG == 0
+ gerror(stderr,
+ "64-bit data support not built (long long not available)\n");
+#endif
+ ((longlong *)obuf)[yoff+ox] = (longlong)vall;
+ break;
+ case TY_FLOAT:
+ ((float *)obuf)[yoff+ox] = (float)vall;
+ break;
+ case TY_DOUBLE:
+ ((double *)obuf)[yoff+ox] = (double)vall;
+ break;
+ default:
+ break;
+ }
+ }
+ else{
+ switch(fun->odtype){
+ case TY_CHAR:
+ ((unsigned char *)obuf)[yoff+ox] = (unsigned char)val;
+ break;
+ case TY_USHORT:
+ ((unsigned short *)obuf)[yoff+ox] = (unsigned short)val;
+ break;
+ case TY_SHORT:
+ ((short *)obuf)[yoff+ox] = (short)val;
+ break;
+ case TY_INT:
+ ((int *)obuf)[yoff+ox] = (int)val;
+ break;
+ case TY_LONG:
+#if HAVE_LONG_LONG == 0
+ gerror(stderr,
+ "64-bit data support not built (long long not available)\n");
+#endif
+ ((longlong *)obuf)[yoff+ox] = (longlong)val;
+ break;
+ case TY_FLOAT:
+ ((float *)obuf)[yoff+ox] = (float)val;
+ break;
+ case TY_DOUBLE:
+ ((double *)obuf)[yoff+ox] = (double)val;
+ break;
+ default:
+ break;
+ }
+ }
+
+ /* skip the blank assignment */
+ continue;
+
+blank:
+ /* arrive here is we found a blank or nan */
+ switch(fun->odtype){
+ case TY_CHAR:
+ ((unsigned char *)obuf)[yoff+ox] = fun->blank;
+ break;
+ case TY_USHORT:
+ ((unsigned short *)obuf)[yoff+ox] = fun->blank;
+ break;
+ case TY_SHORT:
+ ((short *)obuf)[yoff+ox] = fun->blank;
+ break;
+ case TY_INT:
+ ((int *)obuf)[yoff+ox] = fun->blank;
+ break;
+ case TY_LONG:
+#if HAVE_LONG_LONG == 0
+ gerror(stderr,
+ "64-bit data support not built (long long not available)\n");
+#endif
+ ((longlong *)obuf)[yoff+ox] = fun->blank;
+ break;
+ case TY_FLOAT:
+ ((float *)obuf)[yoff+ox] = getnanf();
+ break;
+ case TY_DOUBLE:
+ ((double *)obuf)[yoff+ox] = getnand();
+ break;
+ default:
+ break;
+ }
+
+ } /* x loop */
+ } /* y loop */
+
+ /* free allocated space */
+ for(i=0; i<fun->block; i++){
+ if( bufs[i] ) xfree(bufs[i]);
+ }
+ if( bufs ) xfree(bufs);
+
+ /* calculate the position to the end of the data, if we read last line */
+ if( rstop == MAX(fun->odim2,1) ){
+ newpos = ((MAX(fun->dim2,1)*fun->dim1)) * ABS(fun->dtype);
+ /* skip to end of data */
+ gskip(fun->gio, newpos - fun->curpos);
+ /* skip padding as well */
+ gskip(fun->gio, (off_t)fun->dpad);
+ /* reset current position to end of image */
+ fun->curpos = newpos + fun->dpad;
+ }
+
+ return (void *)obuf;
+}
+
+/*
+ *
+ * semi-public routines
+ *
+ */
+
+
+#ifdef ANSI_FUNC
+void *
+_FunImageMask(Fun fun, void *buf, int rstart, int rstop,
+ FilterMask masks, int nmask, char *plist)
+#else
+void *_FunImageMask(fun, buf, rstart, rstop, masks, nmask, plist)
+ Fun fun;
+ void *buf;
+ int rstart;
+ int rstop;
+ FilterMask masks;
+ int nmask;
+ char *plist;
+#endif
+{
+ int len;
+ int i, n;
+ int pixsize;
+ int rowsize;
+ int flag;
+ int bitpix;
+ int dtype;
+ int x=0, y=0;
+ char *pbuf;
+ char *obuf;
+ char *optr;
+ char tbuf[SZ_LINE];
+
+ /* make sure we have something to do */
+ if( !_FunValid(fun) )
+ return buf;
+
+ /* make sure we really want to mask, and figure out which sort of mask */
+ flag = 0;
+ pbuf = xstrdup(plist);
+ if( _FunKeyword(pbuf, "mask", "FUN_MASK", tbuf, SZ_LINE) ){
+ if( isfalse(tbuf) ){
+ if( pbuf ) xfree(pbuf);
+ return buf;
+ }
+ /* 'all' means we change each pixel, setting pixels to 0 outside,
+ and setting pixels to region value inside */
+ else if( !strcasecmp(tbuf, "all") ){
+ /* if we are allocating the buffer, we can change bitpix as well */
+ if( !buf && _FunKeyword(pbuf, "bitpix", NULL, tbuf, SZ_LINE) ){
+ if( (i = atoi(tbuf)) ){
+ fun->obitpix = i;
+ fun->odtype = fun->obitpix / FT_WORDLEN;
+ }
+ }
+ flag = 1;
+ }
+ }
+ if( pbuf ) xfree(pbuf);
+ bitpix = fun->obitpix;
+ dtype = fun->odtype;
+
+ /* make sure we are not out of bounds */
+ if( rstart < 1 ) rstart = 1;
+ if( rstop > MAX(fun->odim2,1) ) rstop = MAX(fun->odim2,1);
+
+ /* allocate output buffer, if necessary */
+ if( buf )
+ obuf = buf;
+ else{
+ obuf = (char *)xcalloc(fun->odim1*(rstop-rstart+1), ABS(dtype));
+ if( !obuf ){
+ gerror(stderr, "can't allocate image buffer (%dx%dx%d)\n",
+ fun->odim1, (rstop-rstart+1), ABS(fun->odtype));
+ return NULL;
+ }
+ }
+
+ /* this is the size of a single pixel in bytes */
+ pixsize = ABS(bitpix/FT_WORDLEN);
+ /* this is the number of bytes per line */
+ rowsize = fun->odim1 * pixsize;
+
+ /* clear y rows before first row containing a mask segment */
+ if( nmask ){
+ for(y=rstart; y<=MIN(masks[0].y-1,rstop); y++){
+ optr = (char *)obuf + ((y-rstart) * rowsize);
+ memset(optr, 0, rowsize);
+ }
+ }
+
+ /* process all mask segments */
+ for(n=0; n<nmask; n++){
+ /* do not process mask values greater than rstop */
+ if( masks[n].y < rstart ) continue;
+ if( masks[n].y > rstop ) break;
+ /* blank out rows before the y row associated with this mask segment */
+ for(; y<masks[n].y; y++){
+ optr = (char *)obuf + ((y-rstart) * rowsize);
+ memset(optr, 0, rowsize);
+ }
+ /* process each mask segment having this same y value */
+ x = 1;
+ while(1){
+ /* clear columns before the start of this segment */
+ optr = (char *)obuf + ((y-rstart) * rowsize) + ((x-1) * pixsize);
+ if( (len = (masks[n].xstart - x) * pixsize) >0 )
+ memset(optr, 0, len);
+ x = masks[n].xstart;
+ /* if flag is set, we set the pixel values to be the region id */
+ if( flag ){
+ optr = (char *)obuf + ((y-rstart) * rowsize) + ((x-1) * pixsize);
+ for(x=masks[n].xstart; x<=masks[n].xstop; x++){
+ switch(dtype){
+ case TY_CHAR:
+ *((char *)optr) = (char)masks[n].region;
+ optr += pixsize;
+ break;
+ case TY_USHORT:
+ *((unsigned short *)optr) = (unsigned short)masks[n].region;
+ optr += pixsize;
+ break;
+ case TY_SHORT:
+ *((short *)optr) = (short)masks[n].region;
+ optr += pixsize;
+ break;
+ case TY_INT:
+ *((int *)optr) = (int)masks[n].region;
+ optr += pixsize;
+ break;
+ case TY_LONG:
+#if HAVE_LONG_LONG == 0
+ gerror(stderr,
+ "64-bit data support not built (long long not available)\n");
+#endif
+ *((longlong *)optr) = (longlong)masks[n].region;
+ optr += pixsize;
+ break;
+ case TY_FLOAT:
+ *((float *)optr) = (float)masks[n].region;
+ optr += pixsize;
+ break;
+ case TY_DOUBLE:
+ *((double *)optr) = (double)masks[n].region;
+ optr += pixsize;
+ break;
+ default:
+ break;
+ }
+ }
+ }
+ /* position after end of segment */
+ x = masks[n].xstop+1;
+ /* if we have another mask segment with the same y,
+ we keep going */
+ if( ((n+1) < nmask) && (masks[n+1].y == y) )
+ n++;
+ else
+ break;
+ }
+ /* clear the rest of the columns for this y row */
+ optr = (char *)obuf + ((y-rstart) * rowsize) + ((x-1) * pixsize);
+ if( (len = (fun->odim1 - x + 1) * pixsize) >0 )
+ memset(optr, 0, len);
+ /* bump to next y row */
+ y++;
+ }
+
+ /* clear y rows after last row containing a mask segment */
+ for(y = (nmask ? MAX(rstart,masks[nmask-1].y+1) : rstart); y<=rstop; y++){
+ optr = (char *)obuf + ((y-rstart) * rowsize);
+ memset(optr, 0, rowsize);
+ }
+
+ /* return result */
+ return obuf;
+}
+
+
+/*
+ *
+ * _FunImagePutHeader -- writes header to file
+ * used by FunImagePut and sometimes FunFlush (if no Put was done)
+ *
+ */
+#ifdef ANSI_FUNC
+int
+_FunImagePutHeader(Fun fun, int dim1, int dim2, int bitpix)
+#else
+int _FunImagePutHeader(fun, dim1, dim2, bitpix)
+ Fun fun;
+ int dim1, dim2, bitpix;
+#endif
+{
+ int i;
+ int extver=1;
+ int domerge=0;
+ char *extname="IMAGE";
+ FITSHead merge;
+
+ /* if no operation has been performed, we have no header */
+ if( !fun->ops ){
+ /* if we do not have header info yet, we should be able to copy it */
+ if( !fun->header ){
+ /* copy header from reference image */
+ if( fun->ifun )
+ _FunCopy2ImageHeader(fun->ifun, fun);
+ /* use passed dims and bitpix, if possible */
+ if( dim1 && bitpix ){
+ /* note that we may not have a header yet -- make one here */
+ if( !fun->header ){
+ fun->header = ft_headinit(NULL, 0);
+ ft_headsetl(fun->header, "SIMPLE", 0, 1, "FITS STANDARD", 1);
+ }
+ /* add the passed dim and bitpix values */
+ ft_headseti(fun->header, "BITPIX", 0, bitpix, "Bits/pixel", 1);
+ if( dim2 ){
+ ft_headseti(fun->header, "NAXIS", 0, 2, "2-D image", 1);
+ ft_headseti(fun->header, "NAXIS", 1, dim1, "Axis 1 dimension", 1);
+ ft_headseti(fun->header, "NAXIS", 2, dim2, "Axis 2 dimension", 1);
+ }
+ else{
+ ft_headseti(fun->header, "NAXIS", 0, 1, "1-D image", 1);
+ ft_headseti(fun->header, "NAXIS", 1, dim1, "Axis 1 dimension", 1);
+ }
+ ft_headsetl(fun->header, "EXTEND", 0, 1,
+ "Standard extensions might follow", 1);
+ /* synchronize the header and the cards after any changes */
+ ft_syncdata(fun->header);
+ /* fill in values */
+ if( dim2 ){
+ fun->dims = 2;
+ fun->dim1 = dim1;
+ fun->dim2 = dim2;
+ }
+ else{
+ fun->dims = 1;
+ fun->dim1 = dim1;
+ fun->dim2 = 1;
+ }
+ fun->bitpix = bitpix;
+ _FunImageSize(fun);
+ }
+ /* hmmm ... could not make a header of any sort */
+ if( !fun->header )
+ return 0;
+ /* flag that we need to merge in user params before writing */
+ domerge = 1;
+ }
+ /* if something is already written, this is an image extension */
+ if( fun->primio || gtell(fun->gio)>0 ){
+ /* required after naxis */
+ ft_headseti(fun->header, "PCOUNT", 0, 0, "Random parameter count", 1);
+ ft_headseti(fun->header, "GCOUNT", 0, 1, "Group count", 1);
+ /* convert first card to a FITS image extension */
+ ft_cardfmt((FITSCard)fun->header->cards,
+ "XTENSION", 0, FT_STRING, "IMAGE", 0, "FITS Image Extension");
+ if( !ft_headfind(fun->header, "EXTNAME", 0, 0) )
+ ft_headsets(fun->header, "EXTNAME", 0, extname, "Ext. name", 1);
+ if( !ft_headfind(fun->header, "EXTVER", 0, 0) )
+ ft_headseti(fun->header, "EXTVER", 0, extver, "Ext. version", 1);
+ /* delete EXTEND keyword, which is for primary only */
+ if( ft_headfind(fun->header, "EXTEND", 0, 0) )
+ ft_headdel(fun->header, "EXTEND", 0);
+ }
+ else{
+ /* flag that we have written the primary header */
+ fun->primio++;
+ }
+ /* if the data are floats, we must delete BSCALE and BZERO */
+ if( fun->bitpix < 0 ){
+ if( ft_headfind(fun->header, "BSCALE", 0, 0) )
+ ft_headdel(fun->header, "BSCALE", 0);
+ if( ft_headfind(fun->header, "BZERO", 0, 0) )
+ ft_headdel(fun->header, "BZERO", 0);
+ }
+ /* merge in user-specified params with overwrite */
+ if( domerge && fun->theader ){
+ merge = ft_headmerge(fun->header, fun->theader, 1);
+ ft_headfree(fun->theader, 1);
+ fun->theader = NULL;
+ ft_headfree(fun->header, 1);
+ fun->header = merge;
+ ft_syncdata(fun->header);
+ }
+ /* for a 1D image, we have to delete 2D header params */
+ if( fun->dims == 1 ){
+ ft_headdel(fun->header, "NAXIS", 2);
+ /* remove 2D WCS header parameters */
+ for(i=0; i<=2; i++){
+ ft_headdel(fun->header, "CTYPE", i);
+ ft_headdel(fun->header, "CRVAL", i);
+ ft_headdel(fun->header, "CRPIX", i);
+ ft_headdel(fun->header, "CDELT", i);
+ ft_headdel(fun->header, "LTV", i);
+ ft_headdel(fun->header, "LTM1_", i);
+ ft_headdel(fun->header, "LTM2_", i);
+ ft_headdel(fun->header, "CD1_", i);
+ ft_headdel(fun->header, "CD2_", i);
+ ft_headdel(fun->header, "CD3_", i);
+ }
+ }
+
+ /* save file position for later updating */
+ fun->headpos = gtell(fun->gio);
+
+ /* write the header */
+ ft_headwrite(fun->gio, fun->header);
+
+ /* we just wrote the header */
+ fun->ops |= OP_WRHEAD;
+ }
+
+ /* good news */
+ return 1;
+}
+
+/*
+ *
+ * public routines
+ *
+ */
+
+#ifdef ANSI_FUNC
+void *
+FunImageGet (Fun fun, void *buf, char *plist)
+#else
+void *FunImageGet (fun, buf, plist)
+ Fun fun;
+ void *buf;
+ char *plist;
+#endif
+{
+ int domask;
+ int doimage;
+ int imagetable;
+ char *pbuf;
+ char *s;
+ void *obuf=NULL;
+ char paint[SZ_LINE];
+ char debug[SZ_LINE];
+ char tbuf[SZ_LINE];
+
+ /* sanity checks */
+ if( !_FunValid(fun) ){
+ gerror(stderr, "data does not contain a valid header\n");
+ return NULL;
+ }
+
+ /* if we delayed the open before, we have to open now */
+ if( !fun->header && strchr(fun->mode, 'r') )
+ _FunFITSOpen(fun, fun->fname, "r");
+
+ /* sanity checks */
+ if( !fun->header ){
+ gerror(stderr, "data does not contain a valid header\n");
+ return NULL;
+ }
+
+ pbuf = xstrdup(plist);
+ domask = 1;
+ doimage = 1;
+ imagetable = 0;
+ if( _FunKeyword(pbuf, "mask", "FUN_MASK", tbuf, SZ_LINE) ){
+ if( isfalse(tbuf) )
+ domask = 0;
+ /* mask=all => no data, just make the mask */
+ else if( !strcasecmp(tbuf, "all") )
+ doimage = 0;
+ /* mask=image => we want to use image masking on a table */
+ else if( !strcasecmp(tbuf, "image") )
+ imagetable = 1;
+ }
+ if( _FunKeyword(pbuf, "paint", "FILTER_PAINT", tbuf, SZ_LINE) ){
+ strncpy(paint, tbuf, SZ_LINE);
+ }
+ else{
+ *paint = '\0';
+ }
+ if( _FunKeyword(pbuf, "debug", "FILTER_DEBUG", tbuf, SZ_LINE) ){
+ strncpy(debug, tbuf, SZ_LINE);
+ }
+ else{
+ *debug = '\0';
+ }
+ if( pbuf ) xfree(pbuf);
+
+ /* pick the appropriate routine for tables and/or images */
+ if( doimage ){
+ if( fun->header->image ){
+ obuf = _FunImageExtract(fun, buf, 1, MAX(fun->odim2,1), plist);
+ }
+ else if( fun->header->table ){
+ obuf = _FunTableBin(fun, buf, plist);
+ /* if mask=image, we remask using image binning, otherwise not */
+ if( !imagetable )
+ domask = 0;
+ }
+ /* don't bother masking if we have no image */
+ if( obuf == NULL )
+ domask = 0;
+ }
+
+ /* filter resulting image, if we have one */
+ if( domask ){
+ if( fun->filter && *fun->filter )
+ s = fun->filter;
+ else if( !doimage )
+ s = "field";
+ else
+ s = NULL;
+ if( s ){
+ /* open the filter, first time through */
+ if( !fun->filt ){
+ /* make up filter plist string */
+ strncpy(tbuf, "type=image", SZ_LINE-1);
+ /* add the binning key */
+ if( fun->bincols ){
+ strncat(tbuf, ",", SZ_LINE-1);
+ strncat(tbuf, fun->bincols, SZ_LINE-1);
+ }
+ /* add paint mode */
+ if( *paint ){
+ strncat(tbuf, ",", SZ_LINE-1);
+ strncat(tbuf, "paint=", SZ_LINE-1);
+ strncat(tbuf, paint, SZ_LINE-1);
+ }
+ /* add debug flag */
+ if( *debug ){
+ strncat(tbuf, ",", SZ_LINE-1);
+ strncat(tbuf, "debug=", SZ_LINE-1);
+ strncat(tbuf, debug, SZ_LINE-1);
+ }
+ fun->filt = FilterOpen(fun->header, s, tbuf);
+ if( fun->filt && (fun->filt != NOFILTER) ){
+ fun->nmask = FilterImage(fun->filt,
+ fun->x0, fun->x1, fun->y0, fun->y1, fun->block,
+ &(fun->masks), NULL);
+ }
+ }
+ /* perform the filter, if necessary */
+ if( fun->filt && (fun->filt != NOFILTER) ){
+ obuf = _FunImageMask(fun, obuf, 1, MAX(fun->odim2,1),
+ fun->masks, fun->nmask, plist);
+ }
+ /* close filters, free masks */
+ if( fun->filt && (fun->filt != NOFILTER) ){
+ FilterClose(fun->filt);
+ fun->filt = NULL;
+ }
+ if( fun->masks ){
+ xfree(fun->masks);
+ fun->masks = NULL;
+ }
+ fun->nmask = 0;
+ }
+ }
+
+ /* we just read an image */
+ fun->ops |= OP_RDIMAGE;
+
+ /* return result */
+ return obuf;
+}
+
+#ifdef ANSI_FUNC
+void *
+FunImageRowGet (Fun fun, void *buf, int rstart, int rstop, char *plist)
+#else
+void *FunImageRowGet (fun, buf, rstart, rstop, plist)
+ Fun fun;
+ void *buf;
+ int rstart;
+ int rstop;
+ char *plist;
+#endif
+{
+ int domask;
+ int doimage;
+ char *pbuf;
+ void *obuf=NULL;
+ char paint[SZ_LINE];
+ char debug[SZ_LINE];
+ char tbuf[SZ_LINE];
+
+ /* sanity checks */
+ if( !_FunValid(fun) ){
+ gerror(stderr, "data does not contain a valid header\n");
+ return NULL;
+ }
+
+ /* if we delayed the open before, we have to open now */
+ if( !fun->header && strchr(fun->mode, 'r') )
+ _FunFITSOpen(fun, fun->fname, "r");
+
+ /* sanity checks */
+ if( !fun->header ){
+ gerror(stderr, "data does not contain a valid header\n");
+ return NULL;
+ }
+
+ pbuf = xstrdup(plist);
+ domask = 1;
+ doimage = 1;
+ if( _FunKeyword(pbuf, "mask", "FUN_MASK", tbuf, SZ_LINE) ){
+ if( isfalse(tbuf) )
+ domask = 0;
+ else if( !strcasecmp(tbuf, "all") )
+ doimage = 0;
+ }
+ if( _FunKeyword(pbuf, "paint", "FILTER_PAINT", tbuf, SZ_LINE) ){
+ strncpy(paint, tbuf, SZ_LINE);
+ }
+ else{
+ *paint = '\0';
+ }
+ if( _FunKeyword(pbuf, "debug", "FILTER_DEBUG", tbuf, SZ_LINE) ){
+ strncpy(debug, tbuf, SZ_LINE);
+ }
+ else{
+ *debug = '\0';
+ }
+
+ if( pbuf ) xfree(pbuf);
+
+ /* pick the appropriate routine for tables and/or images */
+ if( doimage ){
+ if( fun->header->image ){
+ obuf = _FunImageExtract(fun, buf, rstart, rstop, plist);
+ }
+ /* we do not get lines for tables (they are not sorted) */
+ else if( fun->header->table ){
+ obuf = NULL;
+ domask = 0;
+ }
+ else{
+ obuf = NULL;
+ }
+ if( obuf == NULL )
+ domask = 0;
+ }
+
+ /* filter resulting image, if we have one */
+ if( domask ){
+ if( fun->filter && *fun->filter ){
+ /* open the filter, first time through */
+ if( !fun->filt ){
+ /* make up filter plist string */
+ strncpy(tbuf, "type=image", SZ_LINE-1);
+ /* add the binning key */
+ if( fun->bincols ){
+ strncat(tbuf, ",", SZ_LINE-1);
+ strncat(tbuf, fun->bincols, SZ_LINE-1);
+ }
+ /* add paint mode */
+ if( *paint ){
+ strncat(tbuf, ",", SZ_LINE-1);
+ strncat(tbuf, "paint=", SZ_LINE-1);
+ strncat(tbuf, paint, SZ_LINE-1);
+ }
+ /* add debug flag */
+ if( *debug ){
+ strncat(tbuf, ",", SZ_LINE-1);
+ strncat(tbuf, "debug=", SZ_LINE-1);
+ strncat(tbuf, debug, SZ_LINE-1);
+ }
+ fun->filt = FilterOpen(fun->header, fun->filter, tbuf);
+ if( fun->filt && (fun->filt != NOFILTER) ){
+ fun->nmask = FilterImage(fun->filt,
+ fun->x0, fun->x1, fun->y0, fun->y1, fun->block,
+ &(fun->masks), NULL);
+ }
+ }
+ /* perform the filter, if necessary */
+ if( fun->filt && (fun->filt != NOFILTER) ){
+ obuf = _FunImageMask(fun, obuf, rstart, rstop,
+ fun->masks, fun->nmask, plist);
+ }
+ }
+ }
+
+ /* we just read an image row */
+ fun->ops |= OP_RDIMAGE;
+
+ /* return result */
+ return obuf;
+}
+
+#ifdef ANSI_FUNC
+static int
+_FunImagePut (Fun fun, void *buf, int rstart, int rstop,
+ int dim1, int dim2, int bitpix, char *plist)
+#else
+static int _FunImagePut (fun, buf, rstart, rstop, dim1, dim2, bitpix, plist)
+ Fun fun;
+ void *buf;
+ int rstart;
+ int rstop;
+ int dim1;
+ int dim2;
+ int bitpix;
+ char *plist;
+#endif
+{
+ char tbuf[SZ_LINE];
+ char *pbuf;
+ char *zbuf;
+ int convert;
+ int dlen;
+ int dpad;
+ int tdim;
+ off_t newpos;
+
+ /* sanity checks */
+ if( !_FunValid(fun) )
+ return 0;
+
+ /* this extension is an image */
+ fun->type = FUN_IMAGE;
+
+ /* put header */
+ if( !_FunImagePutHeader(fun, dim1, dim2, bitpix) )
+ return 0;
+
+ /* fill in the blanks, if necessary */
+ if( (dim1 == 0) && (dim2 == 0) && (bitpix == 0) ){
+ bitpix = fun->bitpix;
+ dim1 = fun->dim1;
+ dim2 = fun->dim2;
+ }
+
+ /* see if we override the normal data conversion */
+ pbuf = xstrdup(plist);
+ if( _FunKeyword(pbuf, "convert", NULL, tbuf, SZ_LINE) )
+ convert = istrue(tbuf);
+ else
+ convert = !is_bigendian();
+ if( pbuf ) xfree(pbuf);
+
+ /* determine if we are writing one line or the whole image */
+ if( !rstart )
+ tdim = dim1 * MAX(dim2,1);
+ else
+ tdim = dim1 * (rstop-rstart+1);
+
+ /* convert to IEEE, if necessary */
+ if( convert )
+ swap_data(buf, tdim, bitpix/FT_WORDLEN);
+
+ /* determine data size */
+ dlen = tdim * ABS(bitpix/FT_WORDLEN);
+
+ /* we need to skip to the right row */
+ if( rstart ){
+ newpos = ((rstart-1) * dim1) * ABS(fun->dtype);
+ if( gskip(fun->gio, newpos - fun->curpos) < 0 ){
+ gerror(stderr,
+ "FunImagePut can't seek %d bytes into image data\n",
+ newpos - fun->curpos);
+ return 0;
+ }
+ fun->curpos = newpos;
+ }
+
+ /* write the data */
+ gwrite(fun->gio, buf, sizeof(char), dlen);
+
+ /* update the current position */
+ fun->curpos += dlen;
+
+ /* and the current number of bytes read */
+ fun->bytes += dlen;
+
+ /* write the padding, if we just wrote the last row (or full image) */
+ if( !rstart || (rstop == MAX(dim2,1)) ){
+ if( (dpad = FT_BLOCK - (fun->bytes % FT_BLOCK)) != FT_BLOCK ){
+ zbuf = xcalloc(dpad, sizeof(char));
+ gwrite(fun->gio, zbuf, sizeof(char), dpad);
+ if( zbuf ) xfree(zbuf);
+ /* update current position */
+ fun->curpos += dpad;
+ fun->bytes += dpad;
+ }
+ }
+
+ /* we just wrote an image */
+ fun->ops |= OP_WRIMAGE;
+
+ /* success */
+ return 1;
+}
+
+#ifdef ANSI_FUNC
+int
+FunImagePut (Fun fun, void *buf, int dim1, int dim2, int bitpix, char *plist)
+#else
+int FunImagePut (fun, buf, dim1, dim2, bitpix, plist)
+ Fun fun;
+ void *buf;
+ int dim1;
+ int dim2;
+ int bitpix;
+ char *plist;
+#endif
+{
+ return _FunImagePut(fun, buf, 0, 0, dim1, dim2, bitpix, plist);
+}
+
+#ifdef ANSI_FUNC
+int
+FunImageRowPut (Fun fun, void *buf, int rstart, int rstop, int dim1, int dim2,
+ int bitpix, char *plist)
+#else
+int FunImageRowPut (fun, buf, rstart, rstop, dim1, dim2, bitpix, plist)
+ Fun fun;
+ void *buf;
+ int rstart;
+ int rstop;
+ int dim1;
+ int dim2;
+ int bitpix;
+ char *plist;
+#endif
+{
+ return _FunImagePut(fun, buf, rstart, rstop, dim1, dim2, bitpix, plist);
+}
+