include include include include include "ccdred.h" include "icombine.h" # T_COMBINE -- Combine CCD images. # This task is a copy of IMAGES.IMCOMBINE except that it recognizes the # CCD types and can group images by AMP and SUBSET. It also uses header # keyword translation for the exposure times. procedure t_combine () pointer images # Images pointer extns # Image extensions for each subset pointer subsets # Subsets pointer nimages # Number of images in each subset int nsubsets # Number of subsets pointer outroot # Output root image name pointer plroot # Output pixel list root name pointer sigroot # Output root sigma image name pointer logfile # Log filename bool delete # Delete input images? int i pointer sp, output, plfile, sigma bool clgetb() int clgeti(), clgwrd() real clgetr() include "icombine.com" begin call smark (sp) call salloc (outroot, SZ_FNAME, TY_CHAR) call salloc (plroot, SZ_FNAME, TY_CHAR) call salloc (sigroot, SZ_FNAME, TY_CHAR) call salloc (logfile, SZ_FNAME, TY_CHAR) call salloc (output, SZ_FNAME, TY_CHAR) call salloc (plfile, SZ_FNAME, TY_CHAR) call salloc (sigma, SZ_FNAME, TY_CHAR) call salloc (gain, SZ_FNAME, TY_CHAR) call salloc (snoise, SZ_FNAME, TY_CHAR) call salloc (rdnoise, SZ_FNAME, TY_CHAR) call salloc (logfile, SZ_FNAME, TY_CHAR) # Open the header translation which is needed to determine the # amps, subsets and ccdtypes. Get the input images. # There must be a least one image in order to continue. call clgstr ("instrument", Memc[output], SZ_FNAME) if (Memc[output] == EOS) call error (1, "No 'instrument' translation file specified.") call hdmopen (Memc[output]) call cmb_images (images, extns, subsets, nimages, nsubsets) if (nsubsets == 0) call error (0, "No images to combine") # Get task parameters. Some additional parameters are obtained later. call clgstr ("output", Memc[outroot], SZ_FNAME) call clgstr ("plfile", Memc[plroot], SZ_FNAME) call clgstr ("sigma", Memc[sigroot], SZ_FNAME) call clgstr ("logfile", Memc[logfile], SZ_FNAME) call xt_stripwhite (Memc[outroot]) call xt_stripwhite (Memc[sigroot]) call xt_stripwhite (Memc[logfile]) project = clgetb ("project") combine = clgwrd ("combine", Memc[output], SZ_FNAME, COMBINE) reject = clgwrd ("reject", Memc[output], SZ_FNAME, REJECT) blank = clgetr ("blank") call clgstr ("gain", Memc[gain], SZ_FNAME) call clgstr ("rdnoise", Memc[rdnoise], SZ_FNAME) call clgstr ("snoise", Memc[snoise], SZ_FNAME) lthresh = clgetr ("lthreshold") hthresh = clgetr ("hthreshold") lsigma = clgetr ("lsigma") hsigma = clgetr ("hsigma") grow = clgeti ("grow") mclip = clgetb ("mclip") sigscale = clgetr ("sigscale") delete = clgetb ("delete") # Check parameters, map INDEFs, and set threshold flag if (IS_INDEFR (blank)) blank = 0. if (IS_INDEFR (lsigma)) lsigma = MAX_REAL if (IS_INDEFR (hsigma)) hsigma = MAX_REAL if (IS_INDEFI (grow)) grow = 0 if (IS_INDEF (sigscale)) sigscale = 0. if (IS_INDEF(lthresh) && IS_INDEF(hthresh)) dothresh = false else { dothresh = true if (IS_INDEF(lthresh)) lthresh = -MAX_REAL if (IS_INDEF(hthresh)) hthresh = MAX_REAL } # This is here for backward compatibility. if (clgetb ("clobber")) call error (1, "Clobber option is no longer supported") # Combine each input subset. do i = 1, nsubsets { # Set the output, pl, and sigma image names with subset extension. call strcpy (Memc[outroot], Memc[output], SZ_FNAME) call sprintf (Memc[output], SZ_FNAME, "%s%s") call pargstr (Memc[outroot]) call pargstr (Memc[Memi[extns+i-1]]) call strcpy (Memc[plroot], Memc[plfile], SZ_FNAME) if (Memc[plfile] != EOS) { call sprintf (Memc[plfile], SZ_FNAME, "%s%s") call pargstr (Memc[plroot]) # Use this if we can append pl files. #call pargstr (Memc[Memi[extns+i-1]]) call pargstr (Memc[Memi[subsets+i-1]]) } call strcpy (Memc[sigroot], Memc[sigma], SZ_FNAME) if (Memc[sigma] != EOS) { call sprintf (Memc[sigma], SZ_FNAME, "%s%s") call pargstr (Memc[sigroot]) call pargstr (Memc[Memi[extns+i-1]]) } # Combine all images from the (subset) list. iferr (call icombine (Memc[Memi[images+i-1]], Memi[nimages+i-1], Memc[output], Memc[plfile], Memc[sigma], Memc[logfile], NO, delete)) { call erract (EA_WARN) } call mfree (Memi[images+i-1], TY_CHAR) call mfree (Memi[extns+i-1], TY_CHAR) call mfree (Memi[subsets+i-1], TY_CHAR) } # Finish up. call mfree (images, TY_POINTER) call mfree (extns, TY_POINTER) call mfree (subsets, TY_POINTER) call mfree (nimages, TY_INT) call hdmclose () call sfree (sp) end # CMB_IMAGES -- Get images from a list of images. # The images are filtered by ccdtype and sorted by amplifier and subset. # The allocated lists must be freed by the caller. procedure cmb_images (images, extns, subsets, nimages, nsubsets) pointer images # Pointer to lists of subsets (allocated) pointer extns # Image extensions for each subset (allocated) pointer subsets # Subset names (allocated) pointer nimages # Number of images in subset (allocated) int nsubsets # Number of subsets int list # List of input images bool doamps # Divide input into subsets by amplifier? bool dosubsets # Divide input into subsets by subset parameter? bool extend # Add extensions to output image names? int i, nimage, ccdtype pointer sp, type, image, extn, subset, str, ptr, im #int imtopenp(), imtlen(), imtgetim(), ccdtypecl(), ccdtypes() int imtopenp(), imtlen(), imtgetim() pointer immap() bool clgetb(), streq() begin # Get the input image list and check that there is at least one image. nsubsets = 0 list = imtopenp ("input") nimage = imtlen (list) if (nimage == 0) { call imtclose (list) return } # Determine whether to divide images into subsets and append extensions. #doamps = clgetb ("amps") doamps = false dosubsets = clgetb ("subsets") #extend = clgetb ("extensions") extend = true call smark (sp) call salloc (type, SZ_FNAME, TY_CHAR) call salloc (image, SZ_FNAME, TY_CHAR) call salloc (extn, SZ_FNAME, TY_CHAR) call salloc (subset, SZ_FNAME, TY_CHAR) call salloc (str, SZ_FNAME, TY_CHAR) # Go through the input list and eliminate images not satisfying the # CCD image type. Separate into subsets if desired. Create image # and subset lists. #ccdtype = ccdtypecl ("ccdtype", Memc[type], SZ_FNAME) ccdtype = 0 call clgstr ("ccdtype", Memc[type], SZ_FNAME) call xt_stripwhite (Memc[type]) while (imtgetim (list, Memc[image], SZ_FNAME)!=EOF) { iferr (im = immap (Memc[image], READ_ONLY, 0)) { call erract (EA_WARN) next } #ccdtype = ccdtypes (im, Memc[str], SZ_FNAME) call ccdtypes (im, Memc[str], SZ_FNAME) if (Memc[type] != EOS && !streq (Memc[str], Memc[type])) next Memc[extn] = EOS Memc[subset] = EOS if (doamps) { #call ccdamp (im, Memc[str], SZ_FNAME) Memc[str] = EOS if (extend) call strcat (Memc[str], Memc[extn], SZ_FNAME) call strcat (Memc[str], Memc[subset], SZ_FNAME) } if (dosubsets) { call ccdsubset (im, Memc[str], SZ_FNAME) call strcat (Memc[str], Memc[extn], SZ_FNAME) call strcat (Memc[str], Memc[subset], SZ_FNAME) } for (i=1; i <= nsubsets; i=i+1) if (streq (Memc[subset], Memc[Memi[subsets+i-1]])) break if (i > nsubsets) { if (nsubsets == 0) { call malloc (images, nimage, TY_POINTER) call malloc (extns, nimage, TY_POINTER) call malloc (subsets, nimage, TY_POINTER) call malloc (nimages, nimage, TY_INT) } else if (mod (nsubsets, nimage) == 0) { call realloc (images, nsubsets+nimage, TY_POINTER) call realloc (extns, nsubsets+nimage, TY_POINTER) call realloc (subsets, nsubsets+nimage, TY_POINTER) call realloc (nimages, nsubsets+nimage, TY_INT) } nsubsets = i call malloc (ptr, SZ_FNAME, TY_CHAR) call strcpy (Memc[image], Memc[ptr], SZ_FNAME-1) Memi[images+i-1] = ptr call malloc (ptr, SZ_FNAME, TY_CHAR) call strcpy (Memc[extn], Memc[ptr], SZ_FNAME) Memi[extns+i-1] = ptr call malloc (ptr, SZ_FNAME, TY_CHAR) call strcpy (Memc[subset], Memc[ptr], SZ_FNAME) Memi[subsets+i-1] = ptr Memi[nimages+i-1] = 1 } else { ptr = Memi[images+i-1] nimage = Memi[nimages+i-1] + 1 call realloc (ptr, nimage * SZ_FNAME, TY_CHAR) Memi[images+i-1] = ptr Memi[nimages+i-1] = nimage ptr = ptr + (nimage - 1) * SZ_FNAME call strcpy (Memc[image], Memc[ptr], SZ_FNAME-1) } call imunmap (im) } call realloc (images, nsubsets, TY_POINTER) call realloc (extns, nsubsets, TY_POINTER) call realloc (subsets, nsubsets, TY_POINTER) call realloc (nimages, nsubsets, TY_INT) call imtclose (list) call sfree (sp) end # ICOMBINE -- Combine the CCD images in a list. # This procedure maps the images, sets the output dimensions and datatype, # opens the logfile, and sets IMIO parameters. It attempts to adjust # buffer sizes and memory requirements for maximum efficiency. procedure icombine (images, nims, output, plfile, sigma, logfile, stack, delete) char images[SZ_FNAME-1, nims] # Input images int nims # Number of images in list char output[ARB] # Output image char plfile[ARB] # Pixel list file char sigma[ARB] # Output sigma image char logfile[ARB] # Log filename int stack # Stack input images? bool delete # Delete input images? char errstr[SZ_LINE] int i, j, nimages, intype, bufsize, maxsize, memory, oldsize, stack1, err pointer sp, sp1, in, out[3], offsets, temp, key, tmp int getdatatype() real clgetr() char clgetc() int clgeti(), begmem(), errget(), open(), ty_max(), sizeof() pointer immap(), ic_plfile() errchk ic_imstack, immap, ic_plfile, ic_setout, ccddelete include "icombine.com" define retry_ 98 define done_ 99 begin call smark (sp) # Set number of images to combine. if (project) { if (nims > 1) { call sfree (sp) call error (1, "Cannot project combine a list of images") } tmp = immap (images[1,1], READ_ONLY, 0); out[1] = tmp if (IM_NDIM(out[1]) == 1) call error (1, "Can't project one dimensional images") nimages = IM_LEN(out[1],IM_NDIM(out[1])) call imunmap (out[1]) } else nimages = nims # Convert the nkeep parameter if needed. # Convert the pclip parameter to a number of pixels rather than # a fraction. This number stays constant even if pixels are # rejected. The number of low and high pixel rejected, however, # are converted to a fraction of the valid pixels. nkeep = clgeti ("nkeep") if (nkeep < 0) nkeep = max (0, nimages + nkeep) if (reject == PCLIP) { pclip = clgetr ("pclip") if (pclip == 0.) call error (1, "Pclip parameter may not be zero") if (IS_INDEFR (pclip)) pclip = -0.5 i = nimages / 2. if (abs (pclip) < 1.) pclip = pclip * i if (pclip < 0.) pclip = min (-1, max (-i, int (pclip))) else pclip = max (1, min (i, int (pclip))) } if (reject == MINMAX) { flow = clgetr ("nlow") fhigh = clgetr ("nhigh") if (IS_INDEFR (flow)) flow = 0 if (IS_INDEFR (fhigh)) fhigh = 0 if (flow >= 1) flow = flow / nimages if (fhigh >= 1) fhigh = fhigh / nimages i = flow * nimages j = fhigh * nimages if (i + j == 0) reject = NONE else if (i + j >= nimages) { call eprintf ("Bad minmax rejection parameters\n") call sfree (sp) return } } # Map the input images. bufsize = 0 stack1 = stack retry_ iferr { out[1] = NULL out[2] = NULL out[3] = NULL icm = NULL logfd = NULL call smark (sp1) if (stack1 == YES) { call salloc (temp, SZ_FNAME, TY_CHAR) call mktemp ("tmp", Memc[temp], SZ_FNAME) call ic_imstack (images, nimages, Memc[temp]) project = true } # Map the input image(s). if (project) { if (stack1 == YES) { tmp = immap (Memc[temp], READ_ONLY, 0); out[1] = tmp } else { tmp = immap (images[1,1], READ_ONLY, 0); out[1] = tmp } nimages = IM_LEN(out[1],IM_NDIM(out[1])) call calloc (in, nimages, TY_POINTER) call amovki (out[1], Memi[in], nimages) } else { call calloc (in, nimages, TY_POINTER) do i = 1, nimages { tmp = immap (images[1,i], READ_ONLY, 0); Memi[in+i-1] = tmp } } # Map the output image and set dimensions and offsets. tmp = immap (output, NEW_COPY, Memi[in]); out[1] = tmp if (stack1 == YES) { call salloc (key, SZ_FNAME, TY_CHAR) do i = 1, nimages { call sprintf (Memc[key], SZ_FNAME, "stck%04d") call pargi (i) call imdelf (out[1], Memc[key]) } } call salloc (offsets, nimages*IM_NDIM(out[1]), TY_INT) call ic_setout (Memi[in], out, Memi[offsets], nimages) # Determine the highest precedence datatype and set output datatype. intype = IM_PIXTYPE(Memi[in]) do i = 2, nimages intype = ty_max (intype, IM_PIXTYPE(Memi[in+i-1])) IM_PIXTYPE(out[1]) = getdatatype (clgetc ("outtype")) if (IM_PIXTYPE(out[1]) == ERR) IM_PIXTYPE(out[1]) = intype # Open pixel list file if given. if (plfile[1] != EOS) { tmp = ic_plfile (plfile, NEW_COPY, out[1]); out[2] = tmp } else out[2] = NULL # Open the sigma image if given. if (sigma[1] != EOS) { tmp = immap (sigma, NEW_COPY, out[1]); out[3] = tmp IM_PIXTYPE(out[3]) = ty_max (TY_REAL, IM_PIXTYPE(out[1])) call sprintf (IM_TITLE(out[3]), SZ_IMTITLE, "Combine sigma images for %s") call pargstr (output) } else out[3] = NULL # This is done here to work around problem adding a keyword to # an NEW_COPY header and then using that header in a NEW_COPY. # Open masks. call ic_mopen (Memi[in], out, nimages) # Open the log file. logfd = NULL if (logfile[1] != EOS) { iferr (logfd = open (logfile, APPEND, TEXT_FILE)) { logfd = NULL call erract (EA_WARN) } } if (bufsize == 0) { # Set initial IMIO buffer size based on the number of images # and maximum amount of working memory available. The buffer # size may be adjusted later if the task runs out of memory. # The FUDGE factor is used to allow for the size of the # program, memory allocator inefficiencies, and any other # memory requirements besides IMIO. bufsize = 1 do i = 1, IM_NDIM(out[1]) bufsize = bufsize * IM_LEN(out[1],i) bufsize = bufsize * sizeof (intype) bufsize = min (bufsize, DEFBUFSIZE) memory = begmem ((nimages + 1) * bufsize, oldsize, maxsize) memory = min (memory, int (FUDGE * maxsize)) bufsize = memory / (nimages + 1) } # Combine the images. If an out of memory error occurs close all # images and files, divide the IMIO buffer size in half and try # again. switch (intype) { case TY_SHORT: call icombines (Memi[in], out, Memi[offsets], nimages, bufsize) default: call icombiner (Memi[in], out, Memi[offsets], nimages, bufsize) } } then { err = errget (errstr, SZ_LINE) if (icm != NULL) call ic_mclose (nimages) if (!project) { do j = 2, nimages if (Memi[in+j-1] != NULL) call imunmap (Memi[in+j-1]) } if (out[2] != NULL) { call imunmap (out[2]) call imdelete (plfile) } if (out[3] != NULL) { call imunmap (out[3]) call imdelete (sigma) } if (out[1] != NULL) { call imunmap (out[1]) call imdelete (output) } if (Memi[in] != NULL) call imunmap (Memi[in]) if (logfd != NULL) call close (logfd) switch (err) { case SYS_MFULL: bufsize = bufsize / 2 call sfree (sp1) goto retry_ case SYS_FTOOMANYFILES, SYS_IKIOPIX: if (!project) { stack1 = YES call sfree (sp1) goto retry_ } if (stack1 == YES) call imdelete (Memc[temp]) call fixmem (oldsize) call sfree (sp1) call error (err, errstr) default: if (stack1 == YES) call imdelete (Memc[temp]) call fixmem (oldsize) call sfree (sp1) call error (err, errstr) } } # Unmap all the images, close the log file, and restore memory. # The input images must be unmapped first to insure that there # is a FD for the output images since the headers are opened to # update them. However, the order of the NEW_COPY pointers must # be preserved; i.e. the output depends on the first input image, # and the extra output images depend on the output image. if (!project) { do i = 2, nimages { if (Memi[in+i-1] != NULL) { call imunmap (Memi[in+i-1]) if (delete) call ccddelete (images[1,i]) } } } if (out[2] != NULL) call imunmap (out[2]) if (out[3] != NULL) call imunmap (out[3]) if (out[1] != NULL) call imunmap (out[1]) if (Memi[in] != NULL) call imunmap (Memi[in]) if (stack1 == YES) call imdelete (Memc[temp]) if (delete) call ccddelete (images[1,1]) if (logfd != NULL) call close (logfd) if (icm != NULL) call ic_mclose (nimages) call fixmem (oldsize) call sfree (sp) end # TY_MAX -- Return the datatype of highest precedence. int procedure ty_max (type1, type2) int type1, type2 # Datatypes int i, j, type, order[8] data order/TY_SHORT,TY_USHORT,TY_INT,TY_LONG,TY_REAL,TY_DOUBLE,TY_COMPLEX,TY_REAL/ begin for (i=1; (i<=7) && (type1!=order[i]); i=i+1) ; for (j=1; (j<=7) && (type2!=order[j]); j=j+1) ; type = order[max(i,j)] # Special case of mixing short and unsigned short. if (type == TY_USHORT && type1 != type2) type = TY_INT return (type) end # IC_PLFILE -- Map pixel list file # This routine strips any image extensions and then adds .pl. pointer procedure ic_plfile (plfile, mode, refim) char plfile[ARB] # Pixel list file name int mode # Image mode pointer refim # Reference image pointer pl # IMIO pointer (returned) int i, strlen() bool streq pointer sp, str, immap() begin call smark (sp) call salloc (str, SZ_FNAME, TY_CHAR) call imgimage (plfile, Memc[str], SZ_FNAME) # Strip any existing extensions i = strlen(Memc[str]) switch (Memc[str+i-1]) { case 'h': if (i > 3 && Memc[str+i-4] == '.') Memc[str+i-4] = EOS case 'l': if (i > 2 && streq (Memc[str+i-3], ".pl")) Memc[str+i-3] = EOS } call strcat (".pl", Memc[str], SZ_FNAME) pl = immap (Memc[str], NEW_COPY, refim) call sfree (sp) return (pl) end