nircombine.cl0

# NIRCOMBINE: 17JUL00 KMM expects IRAF 2.11Export or later
# NIRCOMBINE - produce IMCOMBINED IR image of heavily overlapped region
# NIRCOMBINE: 13OCT92 KMM
# NIRCOMBINE: incorporates Valdes' new IMCOMBINE 
# NIRCOMBINE: incorporates frame_num list to combine subset of image database
# NIRCOMBINE: 09JUN93 KMM
# NIRCOMBINE: incoporates imstack to get around 102 image limit
# NIRCOMBINE: 16FEB94 KMM
# NIRCOMBINE: fix so that all input images to imstack have the same dimension
#  02MAR94 KMM  shift internal files to avoid internal use of "type" command
#  06APR94 KMM Replace "type" with "concatenate" or "copy"
#  08APR94 KMM Remove prior imcombine version parts so that "nircomine" is the
#                current version and onircombine is the prior version.
#  27APR94 KMM Modify stack database.
#              Add apply_zero flag to allow choice of 
#                (1) applying zero-point shifts via imarith prior to imcombine,
#                    so that zeropoint adjusted results will be saved via
#                    make_stack+ or save_images+.
#                (2) submitting intensity offsets file directly to imcombine,
#                    so that zeropoint unadjusted results will be saved via
#                    make_stack+ or save_images+.
#                NB: IMCOMBINE conserves flux.  The mean of zero-point offsets
#                    is computed, the adjustments are made relative to this mean
#                    value.  The net result is that the level of the resultant
#		     image shifts by the mean value of the zero-point offsets.
#              Fix so that frame_nums=""," ","all" works to include all frames.
#              Fix bug in recovering intensity offsets from imcombine log;
#                 When less than 10 images are combined, the format is
#                     imname[n] median adjust xoff yoff
#                 However when 10 or more images are combined, the  format is:
#                     imname[ n] median adjust xoff yoff ;  n < 10
#                     imname[nn] median adjust xoff yoff ;  n >= 10
#                 Since the field count varies for ten or more images, the
#                    assumption that field 3 contained the "adjust" value
#                    fails on images 1-9 (median is extracted)
#              Replace fscan with scan from a pipe to get database parameters
#              Note: IRAF V2.10EXPORT April 1992 imcombine has an inconsistency
#                 in the sign convention on applied zeropoint.  The sign of the
#                 file based zeropoint assumes the value to be subtracted,
#                 whereas the reported number is the value to be added.  These
#                 have been reconciled to "value to be added" in the next
#                 release and current IRAFX.  zero_invert flag toggles these
#                 situations.
#  23JUN94 KMM Identify and use smallest and largest COM# when all frames are
#              requested
#  22JUL94 KMM replace fscan with scan from file a key spots
#  08AUG94 KMM utilize printf for formatted output (instead of AWK)
#  18APR96 KMM allow specification of statsec to determine zero offsets
#  17MAR98 KMM typo fixes
#  25JUL98 KMM add global image extension
#              replace access with imaccess where appropriate
#              incorporate imexpr into setpix option
#  17AUG98 KMM modified imexpr to produce explicitly real image for setpix
#              implement option to scale images to a common effective
#                 integration time (to_scale) according to image header
#                 keyword (to_name)
#  15MAY00 KMM replace default for reject_value with -1e7 (was -1e8) to
#                 get setpix to work with floating point images
#              modify geotran paramerization to match new geotran lexicon
#  17JUL00 KMM modify parsing of IMCOMBINE output in response to new IMCOMBINE
#                 output format
             
procedure nircombine (match_name,infofile)

string match_name   {prompt="name of resultant composite image"}
file   infofile     {prompt="file produced by XYGET|XYLAP|XYADOPT"}

string frame_nums   {"",
                   prompt="Include only selected frame numbers in MOSAIC|@list"}
file   outfile      {"", prompt="Output information file name"}
# trim values applied to final image
string trimlimits   {"[2:2,2:2]",
                      prompt="added trim limits for input subrasters"}
bool   register     {yes,prompt="Maintain input image origin and size"}
string common       {"median",enum="none|adopt|mean|median|mode",
              prompt="Pre-combine common offset: |none|adopt|mean|median|mode|"}
string comb_opt     {"median", enum="average|median",
                       prompt="Type of combine operation: |average|median|"}
string reject_opt   {"none", prompt="Type of pixel rejection operation",
                    enum="none|minmax|ccdclip|crreject|sigclip|avsigclip|pclip"}
string stat_sec     {"overlap",
                       prompt="Image section for calculating statistics"}
real   lthreshold   {-200.,prompt="Rejection floor for imcombine and stats"}
real   hthreshold   {65000.,prompt="Rejection ceiling for imcombine and stats"}
real   blank        {-1000.,prompt="Value if there are no pixels"}
bool   setpix       {no,prompt="Run SETPIX on data?"}
string maskimage    {"badmask", prompt="bad pixel image mask"}
real   svalue       {-1.0e7, prompt="Setpix value (far below lthreshold)"}
bool   fixpix       {no,prompt="Run FIXPIX on data?"}
string fixfile      {"badpix", prompt="badpix file in FIXPIX format"}
real   to_scale     {0.0,
                   prompt="Scale to value (0 = noscale) via to_name"}
string to_name      {"INT_S", prompt="Image header to_name keyword"}
# bool   project      {no,prompt="Project highest dimension of input images?"}
bool   make_stack   {no,prompt="Imstack images prior to imcombine?"}
bool   apply_zero   {yes,prompt="Apply zeropoint to images prior to imcombine?"}
bool   invert_zero  {no, 
                     prompt="Invert database zeropoints prior to application?"}
# prior to IRAF 2.10 one must invert to database zeropoints when using imcombine
bool   save_images  {no,prompt="Save images which were IMCOMBINED via imstack?"}
bool   do_combine   {yes,prompt="IMCOMBINE (else IRMATCH) into final image?"}
bool   verbose      {yes, prompt="Verbose output?"}

bool   size         {no, prompt="Compute image size?"}
string interp_shift {"linear",enum="nearest|linear|poly3|poly5|spline3",
              prompt="IMSHIFT interpolant (nearest,linear,poly3,poly5,spline3)"}
string bound_shift  {"nearest",enum="constant|nearest|reflect|wrap",
                      prompt="IMSHIFT boundary (constant,nearest,reflect,wrap)"}
real   const_shift  {0.0,prompt="IMSHIFT Constant for boundary extension"}
real   gxshift      {0.0, prompt="global xshift for final image"}
real   gyshift      {0.0, prompt="global yshift for final image"}
string goverlap     {"", prompt="Global overlap section (overrides all else)"}
string gsize        {"", prompt="Global image size (overrides all else)"}
string weight       {"none",prompt="Image weights"}
bool   mclip        {no, prompt="Use median, not mean, in clip algorithms"}
real   pclip        {-0.5, prompt="pclip: Percentile clipping parameter"}
int    nlow         {1, prompt="minmax: Number of low pixels to reject"}
int    nhigh        {1, prompt="minmax: Number of high pixels to reject"}
real   lsigma       {3., prompt="Lower sigma clipping factor"}
real   hsigma       {3., prompt="Upper sigma clipping factor"}
real   sigscale     {0.1,
                     prompt="Tolerance for sigma clipping scaling correction"}
int    grow         {0, prompt="Radius (pixels) for 1D neighbor rejection"}
string expname      {"", prompt="Image header exposure time keyword"}
string rdnoise      {"0.", prompt="ccdclip: CCD readout noise (electrons)"}
string gain         {"1.", prompt="ccdclip: CCD gain (electrons/DN)"}

bool   answer       {yes, prompt="Do you want to continue?", mode="q"}
bool   compute_size {yes, 
                       prompt="Do you want to [re]compute image size?",mode="q"}

struct  *list1, *list2

begin

   int    i,stat,nim,maxnim,slen,slenmax,njunk,pos1b,pos1e,ncolsout,nrowsout,
          nxoff0,nyoff0,im_xoffset,im_yoffset,c_opt,maxsizex,maxsizey,
          ncols, nrows, nxhiref, nxloref, nyhiref, nyloref,mincom,maxcom,
          nxhi, nxlo, nyhi, nylo, nxhisrc, nxlosrc, nyhisrc, nylosrc,
          nxlotrim,nxhitrim,nylotrim,nyhitrim,
          nxlolap,nxhilap,nylolap,nyhilap,
          nxhimat,nxlomat,nyhimat,nylomat,
          nxlonew,nxhinew,nylonew,nyhinew,
          nxmat0, nymat0, nxmos0, nymos0, ixs, iys, ndim, max_nim
   real   zoff, rjunk, xin, yin, xmat, ymat, fxs, fys,
          xs, ys, xoff, yoff, g_xshift, g_yshift, xmax, xmin, ymax, ymin,
          xofftran,yofftran,gtxmax, gtxmin, gtymax, gtymin,
          oxmax,oxmin,oymax,oymin,xoffset,yoffset,
          xlo, xhi, ylo, yhi, xshift, yshift, reject_value
   real   roff, delta, rmedian, rmode, avestat, avemean, avemode, avemedian,
          rmean, stddev_mean,stddev_median,stddev_mode
   bool   do_tran, tran, max_tran, prior_tran,fluxtran,found,new_origin,update,
          project, zero_invert, scale_it
   string uniq,imname,tmpname,matchname,slist,sjunk,soffset,encsub,sopt,
          src,srcsub,mos,mossub,mat,matsub,ref,refsub,stasub,lapsec,outsec,
          vigsec,sformat,traninfo,sxrot,syrot,sxmag,symag,sxshift,syshift,
          dbtran,cotran,consttran,geomtran,interptran,boundtran,
          zero, vcheck, reject, combopt, statsec, dbpath,image_nims
   file   info,out,dbinfo,tmp1,tmp2,im1,im2,pre_zero,comblist,
          tmpimg,do1info,do2info,matinfo,newinfo,misc,
          comblog,im_offsets,ncomblist,nimlist,stacklist,combinfo,im_zero
   int    nex
   string gimextn, imextn, imroot
   bool   debug=no
   real   from_scale, mult_by
   struct line = ""

   matchname   = match_name
   info        = infofile
   g_xshift    = gxshift
   g_yshift    = gyshift
   
# get IRAF global image extension
   show("imtype") | translit ("",","," ",delete-) | scan (gimextn)
   nex     = strlen(gimextn)
   
   uniq        = mktemp ("_Ticb")
   newinfo     = mktemp("tmp$icb")
   dbinfo      = mktemp("tmp$icb")
   matinfo     = mktemp("tmp$icb")
   comblist    = mktemp("tmp$icb")
   stacklist   = mktemp("tmp$icb")
   ncomblist   = mktemp("tmp$icb")
   nimlist     = mktemp("tmp$icb")
   do1info     = mktemp("tmp$icb")
   do2info     = mktemp("tmp$icb")
   tmp1        = mktemp("tmp$icb")
   tmp2        = mktemp("tmp$icb")
   combinfo    = mktemp("tmp$icb")
   misc        = mktemp("tmp$icb")
   im_offsets  = mktemp("tmp$icb")
   im_zero     = mktemp("tmp$icb")
   pre_zero    = mktemp("tmp$icb")
   comblog     = mktemp("tmp$icb")
   im1         = uniq // "_im1"
   im2         = uniq // "_im2"
   tmpimg      = uniq // "_tim"

   reject  = reject_opt
   combopt = comb_opt
   
   if (to_scale != 0)
      scale_it = yes
   else
      scale_it = no
      
   i = strlen(matchname)
   if (imaccess(matchname)) {
      print("Image ",matchname," already exists!")
      goto err
   } else if (substr(matchname,i-nex,i) == "."//gimextn)  { # Strip off imextn
      matchname = substr(matchname,1,i-nex-1)
   }  
   if (! access(info)) { 		# Exit if can't find info
      print ("Cannot access info_file: ",info)
      goto err
   }
# establish ID of output info file
   if (outfile == "" || outfile == " " || outfile == "default")
      out = matchname//".src"
   else
      out = outfile
   if (out != "STDOUT" && imaccess(out)) {
      print("Will overwrite output_file ",out,"!")
      if (!answer) goto err
   } else
      print ("Output_file= ",out)

   if (setpix && !imaccess(maskimage)) {	# Exit if can't find mask    
      print ("Cannot access maskimage ",maskimage)
      goto err
   }

   if (svalue >= lthreshold) {
      reject_value = lthreshold - 1.0
      print ("Resetting reject_value from ",svalue," to ",reject_value)
   } else
      reject_value = svalue

# select zero and statsec options
   c_opt = 0				# Establish common offset statistic
   if (common == "none" || common == "" || common == " ") {
      c_opt = 0
      zero="none"
   } else if (common == "adopt") {
      c_opt = 1
      if (apply_zero) {
         zero = "none"
         print ("Note: zero-points will be applied prior to IMCOMBINE.")
      } else {
         zero = "@"//im_zero 
         print ("Note: zero-points will be applied by IMCOMBINE.")
      }
   } else {
      if (common == "mean")   c_opt = 2
      if (common == "median") c_opt = 3
      if (common == "mode")   c_opt = 4
      zero = common
   }

   if (zero != "none") {
      statsec = stat_sec
      if (statsec == "" || statsec == " " || statsec == "overlap")
         statsec = "overlap"
   } else
      statsec = ""

   print("#NOTE: ",common," statistics for data within ",
      lthreshold," to ",hthreshold)

   if (setpix) {
      print ("#SETPIX to ",reject_value," using ", maskimage," mask")
   }

   if (scale_it) {
      print ("#SCALING images to ",to_scale," using ", to_name,
      " header keyword")
   }
   
   print (trimlimits) | translit ("", "[:,]", "    ") |
      scan(nxlotrim,nxhitrim,nylotrim,nyhitrim)

# Transfer appropriate information from reference to output file
   match ("^\#DB",info,meta+,stop-,print-, > dbinfo)

# count number of images, including reference images
   match ("^COM",info,meta+,stop-,print-,> tmp1)
   count(tmp1) | scan(maxcom)     # here maxcom = total number in COMfile
# test for uniqueness of COM lines
   fields (tmp1,1,lines="",quit-,print-) | sort ("",col=0,num-,rev-) |
      unique (>> tmp2)
   count(tmp1) | scan(i)
   if (i < maxcom) {
      i = maxcom - i
      print ("Warning: ",i," COM lines are not unique!")
      goto err
   }
   match ("^COM_000",tmp2,meta+,stop+,print-) | head (nl=1) | scan(sjunk)
   pos1b = stridx("_",sjunk) + 1
   mincom = int(substr(sjunk,pos1b,strlen(sjunk))) # mincom = smallest COM
   tail (tmp2,nl=1) | scan(sjunk)
   pos1b = stridx("_",sjunk) + 1
   maxcom = int(substr(sjunk,pos1b,strlen(sjunk))) # maxcom = largest COM
   delete (tmp2, ver-, >& "dev$null")

# Expand the range of images and extract COMlines from database

   if (frame_nums != "" && frame_nums != " " && frame_nums != "all") {
      print("Selecting image subset: ",frame_nums)
      print (frame_nums) | translit ("", "|", " ") | scan(image_nims)
      expandnim(image_nims,ref_nim=0,max_nim=maxcom,>> nimlist)
      list1 = nimlist
      for (i = 0; fscan(list1,nim) != EOF; i += 1) {
         sjunk = "COM_000" + nim
         match ("^"//sjunk,tmp1,meta+,stop-,print-,>> matinfo)
         list2 = matinfo
         for (slen = 0; fscan(list2,imname) != EOF; slen += 1) {
            pos1b = stridx("_",imname)+1
            stat = int(substr(imname,pos1b,strlen(imname)))
            if (stat == nim) {
               found = yes
               break
            }
            found = no
         }
         if (!found) {
            print("Image path position missing for ",sjunk)
            next
         }
      }
      list1 = ""; list2 = ""
   } else {
      sjunk = mincom//"-"//maxcom
      print("Selecting entire image set: ",sjunk)
      match ("^COM_",info,meta+,stop-,print-, > matinfo)
   }
   delete (tmp1, ver-, >& "dev$null")

# count number of images, excluding reference images
   match ("COM_000",matinfo,meta-,stop+,print-) | count() | scan(maxnim)
   if ((! make_stack) && (maxnim > 102)) {
      print ("Number of images: ",maxnim," exceeds 102 limit.  Use make_stack+")
      goto err
   }
# make list of temporary images
   for (nim = 1; nim <= maxnim; nim += 1) {
      tmpname = uniq//"_000" + nim
      print (tmpname//"."//gimextn, >> comblist)
   }  
   print ("#NOTE: will combine ",maxnim," images...")

# Size of output image
   list1 = tmp1; delete (tmp1, ver-, >& "dev$null")
   match ("out_sec",dbinfo,meta-,stop-,print-) |		# MIG order
      sort (col=1,ignore+,numeric-,reverse+, > tmp1)
   if (fscan(list1, sjunk, sjunk, outsec) != 3) {
      outsec = ""; ncolsout = 256; nrowsout = 256
#         outsec = "[1:256,1:256]"; ncolsout = 256; nrowsout = 256
   } else { 
      print (outsec) | translit ("", "[:,]", "    ") |
         scan(nxlosrc,nxhisrc,nylosrc,nyhisrc)
      ncolsout = nxhisrc - nxlosrc + 1
      nrowsout = nyhisrc - nylosrc + 1
   }
   list1=""; delete (tmp1, ver-, >& "dev$null")
   print ("#NOTE: prior_outsec: ",outsec)

# Get GEOTRAN info
   match ("do_tran",dbinfo,meta-,stop-,print-) | scan(sjunk, sjunk, do_tran)
   if (do_tran) {		# Fetch GEOTRAN info from database file
      print ("Note: images will be GEOTRANed prior to IMSHIFT and COMBINE")
      match ("db_tran",dbinfo,meta-,stop-,print-) | scan(sjunk, sjunk, dbtran)
      if (!access(dbtran)) {
         print ("GEOTRAN database file db_tran ",dbtran," not found!")
         goto err
      }
      match ("begin",dbtran,meta-,stop-,print-) | scan (sjunk,cotran)
      if (!access(cotran)) {
         print ("GEOTRAN database coord file co_tran ",cotran," not found!")
         goto err
      }
      match ("geom_tran",dbinfo,meta-,stop-,print-) |
         scan(sjunk, sjunk, geomtran)
	 
   # handle new geotran lexicon	 
      if (geomtran != "general" && geomtran != "geometric")
         geomtran = "linear"
      else
         geomtran = "geometric"
	 
      match ("interp_tran",dbinfo,meta-,stop-,print-) |
         scan(sjunk, sjunk, interptran)
      match ("bound_tran",dbinfo,meta-,stop-,print-) |
         scan(sjunk, sjunk, boundtran)
      match ("const_tran",dbinfo,meta-,stop-,print-) |
         scan(sjunk, sjunk, consttran)
      if (real(consttran) > reject_value)
         print ("Note: overriding const_tran with ",reject_value)
      match ("fluxconserve",dbinfo,meta-,stop-,print-) |
         scan(sjunk, sjunk, fluxtran)
      match ("xoffset_tran",dbinfo,meta-,stop-,print-) |
         scan(sjunk, sjunk, xofftran)
      match ("yoffset_tran",dbinfo,meta-,stop-,print-) |
         scan(sjunk, sjunk, yofftran)
      match ("max_tran",dbinfo,meta-,stop-,print-) |
         scan(sjunk, sjunk, max_tran)
   }

   time() | scan(line)
# Log parameters 
   print("#DBI ",line," IRCOMBINE:",>> dbinfo)
   print("#DBI    info_file       ",info,>> dbinfo)
   print("#DBI    global_xshift   ",g_xshift ,>> dbinfo)
   print("#DBI    global_yshift   ",g_yshift ,>> dbinfo)
   print("#DBI    global_overlap  ",goverlap,>> dbinfo)
   print("#DBI    global_size     ",gsize,>> dbinfo)
   print("#DBI    more_trimlimits ",trimlimits,>> dbinfo)
   print("#DBI    common_opt      ",common,>> dbinfo)
   print("#DBI    lthreshold      ",lthreshold,>> dbinfo)
   print("#DBI    hthreshold      ",hthreshold,>> dbinfo)
   print("#DBI    fixpix          ",fixpix,>> dbinfo)
   print("#DBI    fixfile         ",fixfile,>> dbinfo)
   print("#DBI    setpix          ",setpix,>> dbinfo)
   print("#DBI    set_mask        ",maskimage,>> dbinfo)
   print("#DBI    set_value       ",reject_value,>> dbinfo)
   print("#DBI    interp_shift    ",interp_shift,>> dbinfo)
   print("#DBI    bound_shift     ",bound_shift,>> dbinfo)
   print("#DBI    const_shift     ",const_shift,>> dbinfo)
   print("#DBI    imcomb_floor    ",lthreshold,>> dbinfo)
   print("#DBI    imcomb_ceiling  ",hthreshold,>> dbinfo)
   print("#DBI    imcomb_blank    ",blank,>> dbinfo)
   print("#DBI    apply_zero      ",apply_zero,>> dbinfo)
   if (do_tran) {
      print("#Note: Images GEOTRANed prior to IMSHIFT and COMBINE",>> dbinfo)
      if (real(consttran) > reject_value)
         print ("#Note: overriding const_tran with ",reject_value,>> dbinfo)
   }

# Put in global shifts
   xoffset = g_xshift; yoffset = g_yshift

   if (outsec == "" || size) {
      if (compute_size || outsec == "")
# Determine corners of minimum rectangle enclosing region
         new_origin = yes
      else
         new_origin = no
   } else
      new_origin = no
   if (new_origin) print ("Will optimize origin.")

# Compute minimum rectangle enclosing region and overlap region
   closure (matinfo,xoffset,yoffset,trimlimits=trimlimits,
      interp_shift=interp_shift,origin=new_origin,verbose+,>> do1info)
   match ("^COM",do1info,meta+,stop-,print-, > newinfo)
   if (verbose) type (newinfo)
   match ("^ENCLOSED_SIZE",do1info,meta+,stop-,print-) |
      scan(sjunk,encsub)
   print (encsub) | translit ("", "[:,]", "    ") |
      scan(nxlomat,nxhimat,nylomat,nyhimat)
   match ("^UNAPPLIED_OFFSET",do1info,meta+,stop-,print-) |
      scan(sjunk,xoffset,yoffset)
   match ("^OVERLAP",do1info,meta+,stop-,print-) |
      scan(sjunk,lapsec,vigsec)
   print (lapsec) | translit ("", "[:,]", "    ") |
      scan(nxlolap,nxhilap,nylolap,nyhilap)
   if (lapsec == "[0:0,0:0]") lapsec = ""

   if (new_origin) {
# Establishes origin at (0,0)
      ncolsout = nxhimat - nxlomat + 1
      nrowsout = nyhimat - nylomat + 1
# Override minimum rectangle
      outsec  = "[1:"// ncolsout //",1:"// nrowsout //"]"
   } else {	# null unapplied offsets since we don't want to apply them
      xoffset = 0
      yoffset = 0
   }

   if (gsize != "") {			# Report image composite size
      print ("#Note: computed composite image size: ", outsec)
      print ("#Note: adopted composite image size: ", gsize)
      print ("#Note: computed composite image size: ", outsec,>> dbinfo)
      outsec = gsize
      print ("#Note: adopted composite image size: ", outsec,>> dbinfo)
   } else {
      print ("#Note: computed composite image size: ", outsec)
   }

   if (goverlap != "") {
      print ("#Note: computed composite overlap region: ", lapsec)
      print ("#Note: adopted composite overlap region: ", goverlap)
      print ("#Note: computed composite overlap region: ", lapsec,>> dbinfo)
      lapsec = goverlap
      print ("#Note: adopted composite overlap region: ", lapsec,>> dbinfo)
   } else {
      print ("Computed composite overlap region: ", lapsec)
      if (lapsec == "" && c_opt > 1) { 
         print ("Null overlap: can't compute statistics")
         goto err
      }
   }

   print("#DBI    overlap_sec     ",lapsec,>> dbinfo)
   print("#DBI    out_sec         ",outsec,>> dbinfo)
        
# Mark offsets as absolute or relative
   if (register)
      print ("# Absolute",>> im_offsets)
   else
      print ("# Relative",>> im_offsets)

   maxsizex = 0; maxsizey = 0
   slenmax = 0; list1 = newinfo; list2 = comblist
# Loop through frames: fixpix|setpix|geotran|imshift|zoffset
   while (fscan(list1,imname,src,nxlosrc,nxhisrc,nylosrc,nyhisrc,
      nxmat0,nymat0,fxs,fys,soffset) != EOF) {
      pos1b = stridx("_",imname)+1
      if (pos1b < 2) {
         print("Image path position (COM_pos) missing at ",imname)
         goto err
      } else
         nim = int(substr(imname,pos1b,strlen(imname)))

      if (nim == 0 ) {	 			# SKIP COM_000 REF IMAGE
         print("#NOTE: skipping COM_000 ",src)
         next
      }

      if (soffset == "INDEF") soffset = "0.0"
      zoff = real (soffset)
      slenmax = max(slenmax,strlen(src))
      srcsub = "["//nxlosrc//":"//nxhisrc//","//nylosrc//":"//nyhisrc//"]"
      nxlomat = nxmat0 + nxlosrc ; nxhimat = nxmat0 + nxhisrc
      nylomat = nymat0 + nylosrc ; nyhimat = nymat0 + nyhisrc

# SETPIX
      if (setpix) {
# replace bad pixels with selected value  
         imexpr("a==0?b:c",im1,maskimage,reject_value,src,dims="auto",
	    outtype="real",verbose-) 
      }  else
         imcopy(src,im1,verbose-)
 
      if (fixpix) fixpix(im1, badpix, verbose-) 	# FIXPIX

      if (do_tran) {
         print (src) | translit ("", "[:,]", "    ") |
            scan(sjunk,nxlonew,nxhinew,nylonew,nyhinew)
         ncols = nxhinew - nxlonew + 1
         nrows = nyhinew - nylonew + 1
         gtxmin = 1 + xofftran; gtxmax = gtxmin + ncols
         gtymin = 1 + yofftran; gtymax = gtymin + nrows
         geotran(im1,im2,dbtran,cotran,geometry=geomtran,
            xin=INDEF,yin=INDEF,xshift=INDEF,yshift=INDEF,
            xout=INDEF,yout=INDEF,xmag=INDEF,ymag=INDEF,
            xrot=INDEF,yrot=INDEF,xmin=gtxmin,xmax=gtxmax,
            ymin=gtymin,ymax=gtymax,xscale=1.,yscale=1.,
            ncols=INDEF,nlines=INDEF,xsample=1.,ysample=1.,
            interpolant=interptran,boundary=boundtran,
            constant=reject_value,flux=fluxtran,nxblock=256,nyblock=256)
         imdelete(im1,verify-,>& "dev$null")
         imrename(im2,im1,verbose-)
      }
      
# shift image appropriate fractional pixel value
      stat = fscan(list2,tmpname)
      imshift(im1,tmpname,fxs,fys,shifts_file="",interp_type=interp_shift,
         boundary_type=bound_shift,constant=const_shift)	# IMSHIFT
      imdelete(im1,verify-,>& "dev$null")

# scale images to common integration time
      if (scale_it) {						# SCALE
         imgets (src, to_name); from_scale = int (imgets.value)
	 if (from_scale > 0.)   
	    mult_by = to_scale/from_scale
	 else
	    mult_by = 1.0
	 if (mult_by != 1.0) {
            imarith (tmpname,"*",mult_by,tmpname,pixtype="real",calctype="real",
	       hparams="")
            if (verbose) {
	       print ("Scaling ", imname, " from ", to_name, "=", from_scale,
	       "to ", to_scale)
	    }       
	 }
      }
	           
# Note: "pre_zero" has zeropoints prior to, and "im_zero" applied by, IMCOMBINE
      if (apply_zero && c_opt == 1) {	 			# OFFSET
         imarith (tmpname,"+",zoff,tmpname,pixtype="",calctype="",hparams="")
         print (zoff,>> pre_zero)
         print ("0.0",>> im_zero)
      } else if (c_opt == 0) {
         print ("0.0",>> pre_zero)
         print ("0.0",>> im_zero)
      } else {
         print ("0.0",>> pre_zero)
         if (invert_zero)
            print (-zoff,>> im_zero)
         else
            print (zoff,>> im_zero)
      }

# Redetermine the source and destination sections
      nxlomat += xoffset; nxhimat += xoffset
      nylomat += yoffset; nyhimat += yoffset
      nxmat0  += xoffset; nymat0  += yoffset
      matsub = "["//nxlomat//":"//nxhimat//","//nylomat//":"//nyhimat//"]"
# Test for out-of-range and fix as needed
      if ((nxmat0 < 0) || (nymat0 < 0) || (nxhimat > ncolsout) ||
         (nyhimat > nrowsout)) {
         print ("#NOTE: ",src," destination ",matsub," out of range!")
         print ("#NOTE: ",src," destination ",matsub," out of range!",
            >> dbinfo)
         nxlosrc += max(0, (1 - nxlomat)); nxlomat = max(1,nxlomat)
         nylosrc += max(0, (1 - nylomat)); nylomat = max(1,nylomat)
         nxhisrc -= max(0, (nxhimat - ncolsout))
         nxhimat = min(ncolsout,nxhimat)
         nyhisrc -= max(0, (nyhimat - nrowsout))
         nyhimat = min(nrowsout,nyhimat)
      }
      matsub = "["//nxlomat//":"//nxhimat//","//nylomat//":"//nyhimat//"]"
      srcsub = "["//nxlosrc//":"//nxhisrc//","//nylosrc//":"//nyhisrc//"]"
      nxlo = nxlolap - nxlomat + nxlosrc
      nxhi = nxhilap - nxhimat + nxhisrc
      nylo = nylolap - nylomat + nylosrc
      nyhi = nyhilap - nyhimat + nyhisrc
# Determine overlap section within imshifted source image
      stasub = "["//nxlo//":"//nxhi//","//nylo//":"//nyhi//"]"
      print (tmpname," ",srcsub," ",zoff," ",matsub," ",stasub, >> do2info)
      if (! make_stack) {
         print (tmpname//srcsub,>> ncomblist)
      } else {
  mossub = "["//1//":"//(nxhisrc-nxlosrc+1)//","//1//":"//nyhisrc-nylosrc+1//"]"
         ixs = nxhisrc - nxlosrc + 1
         iys = nyhisrc - nylosrc + 1
         print (tmpname," ",srcsub," ",mossub," ",ixs,iys,>> ncomblist)
      }
# Establish offset of origin for each frame
      im_xoffset = nxlomat - 1
      im_yoffset = nylomat - 1
      print (im_xoffset,im_yoffset,tmpname//srcsub,>> im_offsets)
      if (verbose) {
         print (tmpname," ",srcsub," ",zoff," ",matsub," ",stasub)
      }
      if (((nxhisrc-nxlosrc) != (nxhimat-nxlomat)) ||
          ((nyhisrc-nylosrc) != (nyhimat-nylomat))) {
         print (tmpname," dimension mismatch!")
         print (tmpname," dimension mismatch!", >> dbinfo)
         beep
      }
      maxsizex = max((nxhisrc-nxlosrc+1),maxsizex)
      maxsizey = max((nyhisrc-nylosrc+1),maxsizey)
   }
if(debug){##DEBUG
print(maxsizex,maxsizey)
goto err
}
# Extract intensity offsets
#      fields(do2info,"3",lines="1-",print-,quit-,> im_zero)
   if (do_combine) {
      if (make_stack) {
         imdelete(tmpimg//","//im1,verify-,>& "dev$null")
#         delete(tmp2,verify-,>& "dev$null")
#         print ("repl("//reject_value//","//maxsizex//")",>> tmp2)
#         imexpr("@"//,tmpimg,outtype=real,dim=maxsizex//","//maxsizey)
         mkpattern (tmpimg,output="",pattern="constant",v1=reject_value,
           v2=0.,title="",pixtype="real",ndim=2,ncols=maxsizex,
           nlines=maxsizey)
         print ("Making uniform size images for IMSTACK ",maxsizex,maxsizey)
         list1 = ncomblist
         stasub = "[1:"//maxsizex//",1:"//maxsizey//"]"
         while (fscan(list1,src,srcsub,mossub,ixs,iys) != EOF) {
           if ((ixs < maxsizex) || (iys < maxsizey)) {
              imcopy(tmpimg,im1,verbose-)
              imcopy(src//srcsub,im1//mossub,verbose+)
              imdelete(src,verify-,>& "dev$null")
              imrename(im1,src,verbose-)
              print(src//stasub,>> stacklist)
           } else {
              print(src//srcsub,>> stacklist)
           }
         }
         imdelete(tmpimg,verify-,>& "dev$null")
         print ("Creating stack_"//matchname//"...")
         imstack ("@"//stacklist,"stack_"//matchname,title="*",pixtype="*")
         delete (ncomblist, ver-, >& "dev$null")
         print ("stack_"//matchname, > ncomblist)
         if (access(comblist) && ! save_images) {
            print ("Deleting excess images...")
            imdelete("@"//comblist,verify-,>& "dev$null")
         }
         project=yes
      } else {
         project=no
      }
      print ("Performing IMCOMBINE: reject= ",reject," combine= ",
         combopt," zero= ",zero," output= ", matchname)
      print ("#Performing IMCOMBINE: reject= ",reject," combine= ",
         combopt," zero= ",zero," output= ", matchname,>> out)
      imcombine("@"//ncomblist,tmpimg,plfil="",sigma="",logfile=comblog,
         combine=combopt,reject=reject,project=project,outtype="real",
         offsets=im_offsets,masktype="none",maskvalue=0,blank=blank,
         scale="none",zero=zero,weight=weight,statsec=statsec,
         lthreshold=lthreshold,hthreshold=hthreshold,
         nlow=nlow,nhigh=nhigh,pclip=pclip,lsigma=lsigma,hsigma=hsigma,
         mclip=mclip,sigscale=sigscale,expname=expname,
         rdnoise=rdnoise,gain=gain,grow=grow)
      if (register) {
         imgets (tmpimg, "i_naxis1"); ncols = int (imgets.value)
         imgets (tmpimg, "i_naxis2"); nrows = int (imgets.value)
         if (ncols != ncolsout || nrows != nrowsout) {
            mkpattern (matchname,output="",pattern="constant",
               v1=blank,v2=0.,title="",pixtype="real",ndim=2,
               ncols=ncolsout,nlines=nrowsout,header=tmpimg)
            nxlosrc = 1; nxhisrc = min(ncols,ncolsout)
            nylosrc = 1; nyhisrc = min(nrows,nrowsout)
            srcsub = "["//nxlosrc//":"//nxhisrc//","//
               nylosrc//":"//nyhisrc//"]"
            print ("# NOTE: Restoring image origin and size: ",
               tmpimg//srcsub," -> ",matchname//srcsub,>> out)
            imcopy(tmpimg//srcsub,matchname//srcsub,verbose+)
         }
      } else
         imrename(tmpimg,matchname,verbose-)
   } else {
      print ("Performing OVERLAY:")
      imdelete(tmpimg,verify-,>& "dev$null")
      mkpattern (tmpimg,output="",pattern="constant",v1=reject_value,v2=0.,
         title="",pixtype="real",ndim=2,ncols=ncolsout,nlines=nrowsout,
         header=tmpimg)
      list1 = do2info
      for (i = 1; fscan(list1,src,srcsub,zoff,matsub) != EOF ; i += 1) {
            imcopy(src//srcsub,tmpimg//matsub,verbose+)
      }
      imrename(tmpimg,matchname,verbose-)
   }

# output information
   delete (out, ver-,>& "dev$null") # delete prior version (we said OK above)
   copy (dbinfo,out,verbose-)
# update intensity offset data
   delete(tmp2//","//matinfo//","//dbinfo,verify-,>& "dev$null")
   match ("^COM_000",newinfo,meta+,stop-,print-,> dbinfo)
# Skip reference frame
   match ("^COM_000",newinfo,meta+,stop+,print-,> matinfo)
   delete (newinfo,verify-,>& "dev$null")
   copy (dbinfo,newinfo,verbose+)

   if (do_combine) {	# extract per_image info from IMCOMBINE log
      match ("_",comblog,meta-,stop-,print-) |
         match ("combine",,meta-,stop+,print-, >> combinfo)  
# Extract number of output categories in IMCOMBINE log to stat
	 match ("Images",comblog,meta-,stop-,print-) | count() | scan(i,stat)
   }
#   
#  Images                              Offsets
#  _Ticb2469gj_001.imh[3:1022,3:1022]  314  330
#  Images                             Median    Zero   Offsets
#  _Ticb2469gj_001.imh[3:1022,3:1022]  691.71  -599.1  314  330

# Extract new database info to "newinfo" and applied zero_points to "tmp1"
   if (c_opt <= 1) {	# extract zoff from input data to avoid roundoff
      fields(matinfo,"1-11",lines="1-",print-,quit-,>> newinfo)
      fields(matinfo,"11",lines="1-",print-,quit-,>> tmp1)
   } else {		# extract zoff from IMCOMBINE log
      delete (tmp1,verify-,>& "dev$null")
      if (do_combine) {	# extract zoff to tmp1 and offsets to combinfo
         list1 = combinfo
         nim = fscan (list1,imname,src,mat,mos,ref,sjunk)
         if (stat == 4) {		# Images Median Zero Offsets 
            if (nim == 5) 		# imname[nn] median zero xoff yoff
               fields (combinfo,"3",print-,quit-,>> tmp1)
            else if (nim == 6)	# imname[ nn] median zero xoff yoff
               fields (combinfo,"4",print-,quit-,>> tmp1)
            else {
               print ("Warning: incorrect field count in IMCOMBINE log! ",
	          imname," ",c_opt,nim,stat)
               print ("Warning: incorrect field count in IMCOMBINE log! ",
	          imname," ",c_opt,nim,stat,>> out)
               fields (combinfo,"3",print-,quit-,>> tmp1)
            }
         } else if (stat == 3) {	# Images Zero Offsets 
            if (nim == 4) 		# imname[nn] zero xoff yoff
               fields (combinfo,"2",print-,quit-,>> tmp1)
            else if (nim == 5)		# imname[ nn] zero xoff yoff
               fields (combinfo,"3",print-,quit-,>> tmp1)
            else {
               print ("Warning: incorrect field count in IMCOMBINE log! ",
	          imname," ",c_opt,nim,stat)
               print ("Warning: incorrect field count in IMCOMBINE log! ",
	          imname," ",c_opt,nim,stat,>> out)
               fields (combinfo,"2",print-,quit-,>> tmp1)
            }
         } else if (stat == 2) {         # Images Offsets
 	    print ("0.00", >> tmp1)
            print ("Warning: incorrect field count in IMCOMBINE log! ",
	       imname," ",c_opt,nim,stat)
            print ("Warning: incorrect field count in IMCOMBINE log! ",
	       imname," ",c_opt,nim,stat,>> out)	    
	 } else {
               print ("Warning: incorrect field count in IMCOMBINE log! ",
	          imname," ",c_opt,nim,stat)
               print ("Warning: incorrect field count in IMCOMBINE log! ",
	          imname," ",c_opt,nim,stat,>> out)
         }
      } else {
         fields (do2info,"3",print-,quit-,>> tmp1)
      }
      
# check for inconsistency in number of image lines      
      count (matinfo) | scan(ixs)
      count (tmp1)    | scan(iys)
      if (ixs != iys) {
         print("Warning: lines in INFO: ",ixs," not matched to ZOFF: ",iys)
         print("Warning: lines in INFO: ",ixs," not matched to ZOFF: ",iys,
	    >> out)	 
      } 
      fields(matinfo,"1-10",lines="1-",print-,quit-) |
         join("STDIN",tmp1,out=newinfo,missing="Missing",delim=" ",
             maxchars=161,shortest-,verbose+)
   }

# update database with new COM_nnn
   sformat = '%-7s %'//-slenmax//'s %3d %3d %3d %3d %4d %4d %5.2f %5.2f %9.3f\n'
   list1 = ""; list1 = newinfo
   for (i = 0; fscan(list1,imname,src,nxlosrc,nxhisrc,nylosrc,nyhisrc,
         nxmat0,nymat0,xs,ys,soffset) != EOF; i += 1) {
      printf(sformat,imname,src,nxlosrc,nxhisrc,nylosrc,nyhisrc,
            nxmat0,nymat0,xs,ys,real(soffset),>>out)
   }

   if (do_combine) {	# output IMCOMBINE per_image log
      concatenate (comblog,out,append+)
   }
   if (do_combine && make_stack) {
# output STA_nnn to database for use by RECOMBINE
# format STA_nnn xoff yoff imcombine_applied_zero imstack_applied_zero
#        COM_nnn original_source_image total_applied_zero
# Assumes stat contains the number of output categories in IMCOMBINE log
# Assumes IMCOMBINE applied zero_points are in "tmp1" and in "im_zero"
# Assumes pre-IMCOMBINE applied zero_points are in "pre_zero"
      delete (tmp2//","//matinfo//","//do1info,verify-,>& "dev$null")
# extract COM_nnn and original source for later ID
      match ("^COM_000",newinfo,meta+,stop+,print-) |
         fields("STDIN","1,2,11",lines="1-",print-,quit-,>> do1info)
      list1 = ""; list1 = combinfo
      for (i = 1; fscan(list1,imname,src,mat,mos,ref,sjunk) != EOF; i += 1) {
         tmpname = "STA_000" + i
         nim = nscan ()
# format STA_nnn xoff yoff imcombine_applied_zero imstack_applied_zero
         if (c_opt == 0) {		# common = none 
            if (nim == 3)  		# imname[nn] xoff yoff
               print (tmpname," ",src," ",mat,>> tmp2)
            else if (nim == 4) 	# imname[ nn] xoff yoff
               print (tmpname," ",mat," ",mos,>> tmp2)
            else {
               print ("Warning: incorrect field count in IMCOMBINE log! ",
	          imname," ",tmpname," ",c_opt,nim,stat)
               print ("Warning: incorrect field count in IMCOMBINE log! ",
	          imname," ",tmpname," ",c_opt,nim,stat,>> out) 
            }
         } else if (c_opt == 1) {	# common = adopt
            if (nim == 4) 		# imname[nn] zoff xoff yoff
               print (tmpname," ",mat," ",mos,>> tmp2)
            else if (nim == 5)	# imname[ nn] zero xoff yoff
               print (tmpname," ",mos," ",ref,>> tmp2)
            else {
               print ("Warning: incorrect field count in IMCOMBINE log! ",
	          imname," ",tmpname," ",c_opt,nim,stat)
               print ("Warning: incorrect field count in IMCOMBINE log! ",
	          imname," ",tmpname," ",c_opt,nim,stat,>> out) 
            }
         } else {			# common = median|mean|mode
            if (nim == 5)		# imname[nn] median zero xoff yoff
               print (tmpname," ",mos," ",ref,>> tmp2)
            else if (nim == 6)	# imname[ nn] median zero xoff yoff
               print (tmpname," ",ref," ",sjunk,>> tmp2)
            else {
               print ("Warning: incorrect field count in IMCOMBINE log! ",
	          imname," ",tmpname," ",c_opt,nim,stat)
               print ("Warning: incorrect field count in IMCOMBINE log! ",
	          imname," ",tmpname," ",c_opt,nim,stat,>> out)
            }
         }
      }
      delete (newinfo,verify-,>& "dev$null")
      if (invert_zero) {
         rename (im_zero,newinfo)
         list1 = newinfo
         while (fscan(list1,sjunk) != EOF) {
            print (-real(sjunk),>> im_zero)
         }
         delete (newinfo,verify-,>& "dev$null")
      }
      join(tmp2,im_zero,out="STDOUT",miss="Missing",delim=" ",maxch=161,short-,
         ver+) | join("STDIN",pre_zero,out="STDOUT",miss="Missing",delim=" ",
         maxch=161,short-,ver+) | join("STDIN",do1info,out=newinfo,
         miss="Missing",delim=" ",maxch=161,short-,verb+)
      print ("STA_000 stack_"//matchname," ",out," ",ncolsout,nrowsout,
         slenmax,>> out)
   # Fancy format
      sformat = '%-7s %3d %3d %9.3f %9.3f %-7s %'//slenmax//'s'//' %9.3f\\n'
      list1 = ""; list1 = newinfo
      for (i = 0; fscan(list1,src,nxlo,nxhi,xs,ys,mat,ref,fxs) != EOF; i += 1) {
         printf(sformat,src,nxlo,nxhi,xs,ys,mat,ref,fxs,>>out)
      }
   }
 
   if (verbose)    concatenate (do2info,out,append+)
   
# Finish up

err:  list1 = ""; list2 = ""
   imdelete(tmpimg//","//im1//","//im2, verify-,>& "dev$null")
   if (access(comblist) && ! save_images)
      imdelete("@"//comblist,verify-,>& "dev$null")
   delete (comblist//","//tmp1//","//tmp2//","//nimlist,ver-,>& "dev$null")
   delete (combinfo//","//do1info//","//do2info,ver-,>& "dev$null")
   delete (ncomblist//","//comblog,ver-,>& "dev$null")
   delete (misc//","//dbinfo//","//matinfo,ver-,>& "dev$null")
   delete (stacklist//","//pre_zero,ver-,>& "dev$null")
   delete (im_zero//","//im_offsets//","//uniq//"*",ver-,>& "dev$null")
   
end