# quickfringe.cl -- interactively adjust defringing ...and fast!

# 950322: created from tpl.cl
# 970227: added iterstat determination of RMS of each defringing test
# 010727: added option to query for first scaling

procedure quickfringe(fiximage,fringe,outname)

string fiximage   {prompt="Name of (original) image to be defringed"}
string fringe     {prompt="Fringe image to scale and subtract"}
string outname    {prompt="Name of defringed image to create"}

bool   dis=yes   {prompt="Display fringe and original image at start?",mode="q"}
bool   fit=yes    {prompt="Adjust parameters and refit image?",mode="q"}
real   scale      {prompt="Scaling to use on fringe image",mode="q"}
int    z          {0,prompt="Number of image to save?",mode="q"}
bool   echog        {yes,prompt="Beep when task is done?"}

begin

string logimex, logsurf, countfile, rrootname, iimagelist, ooutname, ffunc
string wweighting, fiximg, ffringe, titlestring
real max, nncols, nnlines, sscale, indiv_exptime, avg_exptime, sum_exptime
real skyfring, ll, ul, mn, sig, nx, npx
int x, y, zz, xxord, yyord, xx, m
bool go, ffit, xxterms, ddis

fiximg=fiximage
ffringe=fringe
ooutname=outname
ddis=dis
zz=0

# check for addcomment; if not defined abort task.
go=yes
if(!deftask("addcomment")) {go=no}
if(go) {

# Get initial scaling
    	if(ddis){
		print("===> Displaying fringe image and original image...")
		display(ffringe,1)
		display(fiximg,2)
    	}
        imgets(fiximg,"FRINGMED",>&"dev$null")
	if(imgets.value!="0"){
		sscale=real(imgets.value)
	}
	else {
		print(" ")
		print("===> Calculating initial fringe image scaling...")
#     iterstat(fiximg,nsigrej=5.,maxiter=10,lower=INDEF,upper=INDEF,print-,ver-)
# new scaling = new median from iterstat * old scaling * exptime * avg_exptime
#	sscale=iterstat.median*(real(imgets.value))*indiv_exptime/avg_exptime
        	imgets(fiximg,"SKYFRING")
		if(imgets.value!="0"){
			skyfring=real(imgets.value)
        		imgets(ffringe,"exptime")
			avg_exptime=real(imgets.value)
        		imgets(fiximg,"exptime")
			indiv_exptime=real(imgets.value)
        		imgets(fiximg,"itermed")
        	  sscale=(real(imgets.value))*skyfring*indiv_exptime/avg_exptime
		} else {
			sscale=scale
		}
	}

#need to create countfile
x=1
y=2
ffit=yes
# Fitting loop
while(ffit) {
	if(x!=1){
# get new scaling
		sscale=scale
#		vi(countfile//"_"//str(x))
	}

# scale and subtract the fringe image 
  print ("===> Scaling the fringe image by "//sscale//" and subtracting...") 
#  imar (ffringe, "/", sscale, ffringe//"_"//str(x), ver-)
#  imar (ooutname//"_0", "-", ffringe//"_"//str(x), ooutname//"_"//str(x), ver-)
  imcalc(fiximg//","//ffringe,ooutname//"_"//str(x),"im1-im2*"//sscale,ver-)
  hedit(ooutname//"_"//str(x),"title","Defringed image using scale="//sscale,ver-,add+,update+,show-)

# inspect the subtracted image 
	print ("===> Displaying defringed image #"//str(x)//"...")
	y=y+1
	if(y>4)y=1
	displ (ooutname//"_"//str(x), y)
	print ("===> Computing RMS of defringed image #"//str(x)//"...")
#	iterstat (ooutname//"_"//str(x),nsigrej=5.,maxiter=10,lower=INDEF,upper=INDEF,print-,ver-)
	imstat (ooutname//"_"//str(x),fields="mean,stddev,npix",lower=INDEF,upper=INDEF,for-) | scan(mn,sig,npx)
	m = 1
        while (m <= 10)  {
		ll = mn-5*sig
		ul = mn+5*sig
		imstat (ooutname//"_"//str(x),fields="mean,stddev,npix",lower=ll,upper=ul,for-) | scan(mn,sig,nx)
		if (nx == npx)
                        break
                npx = nx
                m = m + 1
	}
	beep
	print ("===> RMS of defringed image #"//str(x)//" is "//sig)
	beep

# abort loop?
	ffit=fit
	x+=1
}

# clean up
zz=z
if(zz!=0){
  imren(ooutname//"_"//str(zz),ooutname,ver-) 
  hedit(ooutname,"FRINGMED",sscale,add+,ver-,show-,update+)
  addcomment(ooutname,"Defringed using FRINGMED-scaled copy of "//ffringe,ver-)
  imgets(fiximg,"title")
  titlestring=str(imgets.value)
  hedit(ooutname,"title",titlestring,add+,ver-,update+,show-) 
}
else{
  print("No images saved.")
}

y=0
zz=0
while(y<x){
	imdel(ooutname//"_"//str(y),ver-, >& "dev$null")
	y+=1
}

}
else{
print("The tasks `addcomment' are not loaded,")
print("and this task needs them to run.  Please load it and start again.")
}

if(echog)beep

end