Merge commit '339420dd5dd874c41f6bab5808291fb4036dd022' as 'funtools'

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);
+}
+