# version 1.0.5
# task overscan = path$overscan.cl

procedure overscan (inimage, outimage)

string inimage	{prompt="Input images"}
string outimage	{prompt="Output images"}

struct *list0
struct *list1
struct *list2

begin
	int  imgmin2,imgmax2
	
	int  limgmin1,limgmax1,lovsmin1,lovsmax1
	real lgain
	int  rimgmin1,rimgmax1,rovsmin1,rovsmax1
	real rgain

	int  oimgmin1,oimgmax1
	int  olimgmax1,orimgmin1

	real lovslevel,rovslevel
	real lovssigma,rovssigma

	string inimg,outimg
	string tmp0,tmp1,tmp2,tmp3,tmp4

	tmp3 = mktemp("tmp$tmp3")
	sections(inimage, option="fullname", > tmp3)
	list0 = tmp3
	tmp4 = mktemp("tmp$tmp4")
	sections(outimage, option="fullname", > tmp4)
	list1 = tmp4

	while( (fscan(list0,inimg) != EOF) && (fscan(list1,outimg) != EOF))
	{
	# image starts from 1
	imgmin2  = 1
	limgmin1 = 1
	oimgmin  = 1

	# get parameters from image header
	imgets(image=inimg, param="N2XIS2")
	imgmax2  = int(imgets.value)

	imgets(image=inimg, param="H_OSMIN1")
	lovsmin1 = int(imgets.value)
	imgets(image=inimg, param="H_GAIN1")
	lgain    = real(imgets.valu)

	imgets(image=inimg, param="N2XIS1")
	rimgmax1 = int(imgets.value)
	imgets(image=inimg, param="H_OSMAX1")
	rovsmax1 = int(imgets.value)
	imgets(image=inimg, param="H_GAIN2")
	rgain    = real(imgets.value)

	# calculate the rest parameters
	limgmax1 = lovsmin1 - 1 
	rimgmin1 = rovsmax1 + 1
	lovsmax1 = rimgmax1 / 2
	rovsmin1 = lovsmax1 + 1

	oimgmax1 = limgmax1 + rimgmax1 - rovsmax1
	olimgmax1= limgmax
	orimgmin1= olimgmax + 1

	# not to use the next line ... to avoid negative tail of strong pixel
	lovsmin1 = lovsmin1 + 1
	rovsmax1 = rovsmax1 - 1

	# calculate average pixel value of overscan region
	tmp0 = mktemp("tmp$tmp0")
	imstat(images=inimg//"["//lovsmin1//":"//lovsmax1//","//imgmin2//":"//imgmax2//"]",fields="mean,stddev",format-, > tmp0)
	list2=tmp0
	while(fscan(list2,lovslevel,lovssigma)!=EOF){}
	del(tmp0)

	tmp0 = mktemp("tmp$tmp0")
	imstat(images=inimg//"["//rovsmin1//":"//rovsmax1//","//imgmin2//":"//imgmax2//"]",fields="mean,stddev",format-, > tmp0)
	list2 = tmp0
	while(fscan(list2,rovslevel,rovssigma)!=EOF){}
	del(tmp0)

	# print overscan level / stddev
	print("Overscan mean / stddev of ",inimg," = ",lovslevel,"/ ",lovssigma,rovslevel,"/ ",rovssigma)

	# copy the real image region
	tmp1 = mktemp("home$tmp1")
	tmp2 = mktemp("home$tmp2")

	imcopy(inimg//"["//limgmin1//":"//limgmax1//","//imgmin2//":"//imgmax2//"]", tmp1, verbose-)
	imcopy(inimg//"["//rimgmin1//":"//rimgmax1//","//imgmin2//":"//imgmax2//"]", tmp2, verbose-)

	# overscan subtraction	
	imarith(operand1=tmp1,op="-",operand2=lovslevel,result=tmp1,pixtype="r")
	imarith(operand1=tmp2,op="-",operand2=rovslevel,result=tmp2,pixtype="r")

	# gain correction
	imarith(operand1=tmp1,op="*",operand2=lgain,result=tmp1,pixtype="r")
	imarith(operand1=tmp2,op="*",operand2=rgain,result=tmp2,pixtype="r")

	# join images
	tmp0 = mktemp("tmp$tmp0")
	print(tmp1, >> tmp0)
	print(tmp2, >> tmp0)
	imjoin(input="@"//tmp0, output=outimg, join_dim=1, ver-)

	# correct header
	hedit(outimg, fields="N2XIS1", value=oimgmax1, ver-, show-)
	hedit(outimg, fields="PRD-RNG1", value=oimgmax1, ver-, show-)
	hedit(outimg, fields="GAIN", value="1.00", ver-, show-)

	# clean up
	del(tmp0)
	imdel(tmp1, ver-, >& "dev$null")
	imdel(tmp2, ver-, >& "dev$null")
	}
	# clean up
	del(tmp3)
	del(tmp4)
end