procedure crzap(inlist,outlist)

#  Cosmic ray cleaning routine using median filtering.
#
#  970227: created from dimsum.xzap (original courtesy Mark Dickinson)
#  980622: iterstat calculation not done if iterstat header keywords extant
#  980709: hadd parameter: ensure iterstat params not wrongly added to header
#  990516: use imcalc instead of imreplace where appropriate
#  990520: fix inadvertent imcalc bug
#  990530: add CRMASK to output image header only
#  010730: added explicit .fits suffix to _peaks to force new image creation
#
#  Calls scripts fileroot.cl, iterstat.cl; possibly makeobjmask.cl, minv.cl

string	inlist	    {prompt="Image(s) for cosmic ray cleaning"}
string	outlist	    {prompt="Output image(s)"}
int	zboxsz	    {7,min=3,prompt="Box size for zapping"}
int	oboxsz {7,min=3,prompt="Box size for smoothing before object flagging"}
real	nsigma	  {5.,min=0.,prompt="Number of sky sigma for zapping threshold"}
real	nnegsigma {0.,min=0.,prompt="Number of sky sigma for negative zapping"}
int	nrings	{1,min=0,prompt="Number of pixels to flag as buffer around CRs"}
real    nobjsigma   {2.,min=0.,prompt="Number of sky sigma for object masking"}
int skyfiltsize {15,min=0,prompt="Median filter size for sky in object masking"}
int skysubsample {1,min=0,prompt="Block averaging factor for object masking"}
int ngrowobj {2,min=0,prompt="Number of buffer pixels to grow around objects"}
string	statsec	    {"",prompt="Image section to use for computing sky sigma"}
bool	deletemask  {yes, prompt="Delete CR mask after execution?"}
bool	cleanpl     {yes,prompt="Delete other working masks after execution?"}
bool	cleanimh    {yes,prompt="Delete working images after execution?"}
bool	verbose     {yes,prompt="Verbose output?"}
bool  	unzap       {no,prompt="Unzap using objects from objmask files?"}
bool   checklimits {yes,prompt="Check min and max pix values before filtering?"}
bool  	hadd	    {no,prompt="Add iterstat results to header if missing?"}
real	zmin        {-99999.,prompt="Minimum data value for fmedian"}
real	zmax        {99999.,prompt="Minimum data value for fmedian"}
bool	echog	    {no,prompt="Beep when task is done?"}

struct	*inimglist
struct	*outimglist
struct	*statlist

begin

	real	skysig,skymode,cr_thresh,objthresh,negthresh,dmin,dmax
	int	nbox,nbox2,nin
	string 	infile, outfile, img, outimg, statfile, objmask
	string  itmean, itsig, origimg
	bool	maskobj, maskneg

# Check input parameters.
	if (nnegsigma > 0.) 
		maskneg = yes 
	   else 
		maskneg = no
	if (nobjsigma == 0.) {
		print ("No object masking will be used during zapping.")
		maskobj = no
	   } else {
		maskobj = yes
	}

# Expand file lists into temporary files.
	infile =  mktemp("tmp$zap")
	outfile = mktemp("tmp$zap")
	sections (inlist,option="fullname",>infile)
	nin = sections.nimages
	sections (outlist,option="fullname",>outfile)
	if (nin != sections.nimages) {
	print ("ERROR:   Numbers of input and output images do not match.")
		return
	}
	inimglist = infile
	outimglist = outfile

# Loop through input files:
	while (fscan(inimglist,origimg) != EOF) {

# Strip extension off input file name.
		fileroot (origimg,validim+)
		img = fileroot.root

# Read name of output file.
		if (fscan(outimglist,outimg) == EOF) {
		print ("ERROR:   Problem scanning output file name in crzap.")
			return
		}
		fileroot (outimg,validim+)
		outimg = fileroot.root

		print ("Working on ",img," ---ZAP!--> ",outimg)

# Calulate sky mean and RMS noise using iterative sigma rejection.
		if (verbose) print ("     Computing statistics.")
		imgets(img, "itermean", >& "dev$null")
		itmean=str(imgets.value)
        	if(itmean!="0") {
			skymode = real(imgets.value)
			imgets(img, "itersig", >& "dev$null")
			skysig = real(imgets.value)
		} else {
			iterstat (img//statsec,nsigrej=5.,maxiter=10,
				verbose=verbose,print=verbose,addh=hadd)
			skymode = iterstat.mean
			skysig  = iterstat.sigma
		}

# Median filter image to produce file _fmed_//img
		if (verbose) print ("     Median filtering.")

# First, check limits of data range.  Wildly large (positive or negative) data
# values will screw up fmedian calculation unless checklimits = yes and zmin 
# and zmax are set to appropriate values, e.g. zmin=-32768, zmax=32767.   
# Note that the old fmedian bug which occurred if zmin=hmin and zmax=hmax has 
# been fixed in IRAF version 2.10.1.  The pipe to dev$null, however, is 
# included because V2.10.1 and .2 have debugging print statements accidentally
# left in the fmedian code.  

		if (checklimits) {
			minmax (img,force+,update+,ve-)
if (verbose) print("Data minimum = ",minmax.minval, " maximum = ",minmax.maxval)
		}
	if (checklimits && (minmax.minval < zmin || minmax.maxval > zmax)) {
			if (minmax.minval < zmin) {
				dmin = zmin
			   } else { 
				dmin = minmax.minval
			}
			if (minmax.maxval > zmax) {
				dmax = zmax
			   } else { 
				dmax = minmax.maxval
			}
	if (verbose) print ("     Truncating data range ",dmin," to ",dmax)
			if(!access("_fmed_"//origimg)) {
				 fmedian (img,"_fmed_"//img,xw=zboxsz,yw=zboxsz,
				 boundary="reflect",hmin=-32768,hmax=32767,
				 zmin=dmin,zmax=dmax,>&"dev$null")
		   	}
		   } else {
			if(!access("_fmed_"//origimg)) {
	fmedian (img,"_fmed_"//img,xw=zboxsz,yw=zboxsz,boundary="reflect",
		 hmin=-32768,hmax=32767,zmin=INDEF,zmax=INDEF,>&"dev$null")
		   	}
		}

# Take difference to produce "unsharp masked" image _cr_//img
		if(!access("_cr_"//origimg)) {
			imarith (img,"-","_fmed_"//img,"_cr_"//img)
		}
# Threshold _cr_//img at nsigma*skysig to make CR mask _peaks_//img.  
#    Potential CRs --> 1
#    "Blank sky"   --> 0
# Note that cr_thresh will be positive by definition.
		if (verbose) print ("     Masking potential CR events.")
		if (!access("_peaks_"//img//".pl")) {
		    imcopy ("_cr_"//img,"_peaks_"//img//".fits",ve-)
		    cr_thresh = nsigma*skysig   
		    #imreplace ("_peaks_"//img,0.,upper=cr_thresh,lower=INDEF)
		    #imreplace ("_peaks_"//img,1.,lower=cr_thresh,upper=INDEF)
		    imcalc ("_peaks_"//img//".fits", "_peaks_"//img//".fits", "if im1 .le. "//cr_thresh//" then 0 else 1", ver-)
# Make PL format version of _peaks_//img and delete old format version.
		    imcopy ("_peaks_"//img//".fits","_tmpeaks_"//img//".pl",ve-)
		    imdelete ("_peaks_"//img//".fits",ve-)
		    imren    ("_tmpeaks_"//img//".pl","_peaks_"//img//".pl")
		}

# Object masking:  create mask identifying where objects might be,
# inverting mask after creation to make  objects --> 0 and "sky" --> 1.
		if (maskobj) {
			if (verbose) print ("     Creating object mask.")
		        objthresh = nobjsigma*skysig  

# If smoothing scale for object detection is different, create another image
			if (!access("_obj_fmed_"//img//".pl")) {
			  if (oboxsz==zboxsz) {
			    makeobjmask ("_fmed_"//img,suffix="_obj",headlist="",subs=skysubsample,filts=skyfiltsize,nsmoo=0,threshty="constant",constth=objthresh,ngrow=ngrowobj,statsec=statsec,checklimits=checklimits,verbose=verbose,inv+)
		            imren ("_fmed_"//img//"_obj.pl", "_obj_fmed_"//img//".pl",ver-)
			  } else {
			    fmedian (img,"_omed_"//img,xw=oboxsz,yw=oboxsz,boundary="reflect",hmin=-32768,hmax=32767,zmin=INDEF,zmax=INDEF,>&"dev$null")
			    makeobjmask ("_omed_"//img,suffix="_obj",headlist="",subs=skysubsample,filts=skyfiltsize,nsmoo=0,threshty="constant",constth=objthresh,ngrow=ngrowobj,statsec=statsec,checklimits=checklimits,verbose=verbose,inv+)
			    imren ("_omed_"//img//"_obj.pl", "_obj_fmed_"//img//".pl",ver-)
			  }
			}
		}

# If not masking objects, final CR mask is just _peaks_//img.  If we are
# masking objects, take product of object and CR masks to make _crmask_//img
		if (maskobj) {
		imarith ("_peaks_"//img//".pl","*","_obj_fmed_"//img//".pl",
				 "_crmask_"//img//".pl")
		   } else {
		       imcopy ("_peaks_"//img//".pl","_crmask_"//img//".pl",ve-)
		}

# Grow additional buffer region around identified CRs.
		if (nrings > 0) {
		if (verbose) print ("     Growing mask rings around CR hits.")
			nbox = 2 * nrings + 1
			nbox2 = nbox * nbox
	       imarith ("_crmask_"//img//".pl","*",nbox2,"_crmask_"//img//".pl")
			boxcar ("_crmask_"//img//".pl","_crmask_"//img//".pl",
				xw=nbox,yw=nbox)
			imreplace ("_crmask_"//img//".pl",1,lower=1,upper=INDEF)
		}

# Identify negative pixels if desired.  No "rings" are grown around them.
		if (maskneg) {
		if (verbose) print ("     Masking deviant negative pixels.")
			imcopy ("_cr_"//img,"_neg_"//img,ve-)
			negthresh = -1.*nnegsigma*skysig
			#imreplace ("_neg_"//img,0.,lower=negthresh,upper=INDEF)
			#imreplace ("_neg_"//img,1.,upper=negthresh,lower=INDEF)
		        imcalc ("_neg_"//img, "_neg_"//img, "if im1 .ge. "//negthresh//" then 0 else 1", ver-)
			imcopy ("_neg_"//img,"_neg_"//img//".pl",ve-)
			imdelete ("_neg_"//img,ve-)
			imarith ("_crmask_"//img//".pl","+","_neg_"//img//".pl",
				 "_crmask_"//img//".pl")
			imreplace ("_crmask_"//img//".pl",1,lower=1,upper=INDEF)
		}

# unzap pixels which are where objects are as defined by objmask files
               if  (unzap) {
			imgets(img, "objmask", >& "dev$null")
        		objmask=str(imgets.value)
        		if(objmask!="0") {
                     		imarith(objmask,"*","_crmask_"//img//".pl",
                                       "_crmask_"//img//".pl",ver-)
        		} else {
      print ("No object mask defined in image header-- cosmic rays left as is.")
               		}
               }

		if (verbose) print ("     Replacing CR hits with local median.")
# Multiply CR mask by _crmask_//img to make "comic rays only" image _cronly_//img
		imarith ("_cr_"//img,"*","_crmask_"//img//".pl",
			 "_cronly_"//img)

# Subtract _cronly_//img from data to produce clean image "outimg". 
#      Note that this effectively replaces the masked regions with the local
#      median, since _cronly_//img = img - _fmed_//img.
		imarith (img,"-","_cronly_"//img,outimg)

## Record CR mask name in headers of input and output images
#		hedit (img//","//outimg,"CRMASK","_crmask_"//img//".pl",
# Record CR mask name in header of output image
		hedit (outimg,"CRMASK","_crmask_"//img//".pl",
		       add+,ver-,show+,update+,>&"dev$null")

# Clean up.
		if (deletemask) delete ("_crmask_"//img//".pl",ve-)
		if (cleanpl) {
			delete ("_peaks_"//img//".pl",ve-)
			if (maskobj) delete ("_obj_fmed_"//img//".pl",ve-)
			if (maskneg) delete ("_neg_"//img//".pl",ve-)
		}
		if (cleanimh) {
			imdelete ("_fmed_"//img,ver-,>&"dev$null")
			imdelete ("_omed_"//img,ver-,>&"dev$null")
			imdelete ("_cr_"//img,ver-,>&"dev$null")
			imdelete ("_cronly_"//img,ver-,>&"dev$null")
		}
	if (verbose) print ("     "//outimg//" created.")
	}

	inimglist = ""; delete (infile, verify-)
	outimglist = ""; delete (outfile, verify-)
	statlist = ""

if(echog)beep

end