# autodefringe.cl -- automatically defringe using RMS-minimizing scalings

# 980730: created from testfringe.cl
# 010729: speeded up a bit

procedure autodefringe(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 final defringed image to create"}
string logfile    {"STDOUT",prompt="Log file (either STDOUT or a filename)"}
real   startscale {1.,min=0.,prompt="Scaling to start with on fringe image"}
real   stepscale  {0.01,min=0.,prompt="Step between scaling values"}
string imsect     {"",prompt="[Image section] to use for statistics & display"}
bool   addhead    {no,prompt="Add iterstat output to image header?"}
bool   display    {no,prompt="Display original & final images when done?"}
bool   echog        {yes,prompt="Beep when task is done?"}

begin

string llogfile, rrootname, ooutname, ffunc
string wweighting, fiximg, ffringe, titlestring, statsec
real max, nncols, nnlines, sscale, minrms
real sc_start, sc_step, itermed, oldrms, newrms, sc_final, ll, ul, mn, sig
int x, y, zz, xxord, yyord, xx, steps, w, nx, npx, m, v
bool go, ffit, xxterms, ddisplay, hadd, disfinal, getsig

fiximg=fiximage
ffringe=fringe
ooutname=outname

llogfile=logfile
sc_start=startscale
if(sc_start==0.)sc_start=1
sc_step=stepscale
if(sc_step==0.)sc_step=0.01
statsec=imsect
hadd=addhead
ddisplay=display

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

# Get initial RMS
getsig=yes
print("", >> llogfile)
print("===> Getting stats on original image section...", >> llogfile)
if(statsec==""){
	imgets(fiximg,"ITERSIG",>&"dev$null")
	if(imgets.value!="0"){
		sig=real(imgets.value)
		getsig=no
		oldrms=sig
		minrms=sig
	print("===> Using sig="//sig//" from image header...", >> llogfile)
	}
}
if(getsig){
	iterstat (fiximg//statsec,nsigrej=5.,maxiter=10,lower=INDEF,upper=INDEF,print+,ver-,addh=hadd, >> llogfile)
	oldrms=iterstat.sigma
	minrms=iterstat.sigma
#	imstat (fiximg//statsec,fields="mean,stddev,npix",lower=INDEF,upper=INDEF,for-) | scan(mn,sig,npx)
#	m = 1
#        while (m <= 10)  {
#		ll = mn-5*sig
#		ul = mn+10*sig
#		imstat (fiximg//statsec,fields="mean,stddev,npix",lower=ll,upper=ul,for-) | scan(mn,sig,nx)
#		if (nx == npx)
#                        break
#                npx = nx
#                m = m + 1
#	}
#oldrms=sig
#minrms=sig
}
sc_final=0.

# First, increase scalings in positive steps until RMS starts going up;
# then decrease them in negative steps until same thing happens
for (v=1; v>=-1; v+=-2) {
    sscale=sc_start
    sc_step=sc_step*v
    print ("===> Subtracting fringe images scaled from "//sscale//" in "//sc_step//" steps.", >> llogfile) 
# Loop over requested scalings, 100 steps maximum
    for (w=1; w<=100; w+=1) {
# if new image...
	if(w!=1){
# ...get new scaling
		sscale=sscale+sc_step
	} else {
# ...skip redo of initial scaling at start of negative step phase
		if(v==-1){ sscale=sscale+sc_step }
	}
# subtract & iterstat the scaled fringe image 
	imcalc(fiximg//statsec//","//ffringe//statsec,ooutname//"_"//str(sscale),"im1-im2*"//sscale,ver-)
  	hedit(ooutname//"_"//str(sscale),"title","Defringed image using scale="//sscale,ver-,add+,update+,show-)
	imstat (ooutname//"_"//str(sscale),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(sscale),fields="mean,stddev,npix",lower=ll,upper=ul,for-) | scan(mn,sig,nx)
		if (nx == npx)
			break
                npx = nx
		m = m + 1
        }
	print ("RMS of "//ooutname//"_"//str(sscale)//" is "//sig, >> llogfile)
	imdel (ooutname//"_"//str(sscale),ver-)
	newrms=sig
	if(newrms<=minrms){
		minrms=newrms
		sc_final=sscale
	}
	if(newrms>=oldrms){
		w=101
	}
	oldrms=newrms
    }
}

# create full-size final image
if(sc_final!=0){
    print("Creating final image with scale "//sc_final, >> llogfile)
    imcalc(fiximg//","//ffringe,ooutname,"im1-im2*"//sc_final,ver-)
    hedit(ooutname,"FRINGMED",sc_final,add+,ver-,show-,update+)
    addcomment(ooutname,"Defringed using FRINGMED-scaled copy of "//ffringe,ve-)
    if(ddisplay){
	print("===> Displaying original and defringed images...")
	displ(fiximg, 1)
	displ(ooutname, 2)
    }
} else {
    beep
    print ("No defringed image had RMS lower than original!  Wrong fringe image used?", >> llogfile)
}

}
else{
print("Tasks stsdas.toolbox.imgtools.imcalc & phiirs.addcomment aren't loaded,")
print("and this task needs them to run.  Please load them and start again.")
}

if(echog)beep

end