include	<imhdr.h>
include	<imset.h>
include <pkg/gtools.h>
include	<pkg/xtanswer.h>
include	"ccdred.h"


# SET_OVERSCAN -- Set the overscan vector.
#
#   1.  Return immediately if the overscan correction is not requested or
#	if the image has been previously corrected.
#   2.	Determine the overscan columns or lines.  This may be specifed
#	directly or indirectly through the image header or symbol table.
#   3.	Average the overscan columns or lines.
#   4.	Fit a function with the ICFIT routines to smooth the overscan vector.
#   5.  Set the processing flag.
#   6.  Log the operation (to user, logfile, and output image header).

procedure set_overscan (ccd)

pointer	ccd			# CCD structure pointer

int	i, j, nsec, navg, npts, first, last
int	nc, nl, c1, c2, l1, l2
pointer	sp, str, errstr, buf, overscan, x, y, z

real	asumr()
bool	clgetb(), ccdflag()
pointer	imgl2r()
errchk	imgl2r, fit_overscan

begin
	# Check if the user wants this operation or if it has been done.
	if (!clgetb ("overscan") || ccdflag (IN_IM(ccd), "overscan"))
	    return

	call smark (sp)
	call salloc (str, SZ_LINE, TY_CHAR)
	call salloc (errstr, SZ_LINE, TY_CHAR)
	call imstats (IN_IM(ccd), IM_IMAGENAME, Memc[str], SZ_LINE)

	# Check bias section.
	nc = IM_LEN(IN_IM(ccd),1)
	nl = IM_LEN(IN_IM(ccd),2)

	c1 = BIAS_C1(ccd)
	c2 = BIAS_C2(ccd)
	l1 = BIAS_L1(ccd)
	l2 = BIAS_L2(ccd)
	navg = c2 - c1 + 1
	npts = CCD_L2(ccd) - CCD_L1(ccd) + 1

	nsec = max (1, IN_NSEC(ccd))
	do i = 1, nsec {
	    if (BIAS_SEC(ccd) != NULL) {
		c1 = BIAS_SC1(ccd,i)
		c2 = BIAS_SC2(ccd,i)
		l1 = BIAS_SL1(ccd,i)
		l2 = BIAS_SL2(ccd,i)
	    }
	    if ((c1 < 1) || (c2 > nc) || (l1 < 1) || (l2 > nl)) {
		call sprintf (Memc[errstr], SZ_LINE,
		"Error in bias section: image=%s[%d,%d], biassec=[%d:%d,%d:%d]")
		    call pargstr (Memc[str])
		    call pargi (nc)
		    call pargi (nl)
		    call pargi (c1)
		    call pargi (c2)
		    call pargi (l1)
		    call pargi (l2)
		call error (0, Memc[errstr])
	    }
	    if ((c1 == 1) && (c2 == nc) && (l1 == 1) && (l2 == nl)) {
		call error (0,
		    "Bias section not specified or given as full image")
	    }

	    # If no processing is desired then print overscan strip and return.
	    if (clgetb ("noproc")) {
		call eprintf (
		    "  [TO BE DONE] Overscan section is [%d:%d,%d:%d].\n")
		    call pargi (c1)
		    call pargi (c2)
		    call pargi (l1)
		    call pargi (l2)
		    call sfree (sp)
		    return
	    }
	}

	# Determine the overscan section parameters. The readout axis
	# determines the type of overscan.  The step sizes are ignored.
	# The limits in the long dimension are replaced by the trim limits.

	if (READAXIS(ccd) == 1) {
	    call salloc (buf, nsec*nl, TY_REAL)
	    z = buf
	    do i = 1, nl {
		y = imgl2r (IN_IM(ccd), i)
		do j = 1, nsec {
		    if (BIAS_SEC(ccd) != NULL) {
			l1 = BIAS_SL1(ccd,j)
			l2 = BIAS_SL2(ccd,j)
			if (i < l1 || i > l2)
			    next
			c1 = BIAS_SC1(ccd,j)
			c2 = BIAS_SC2(ccd,j)
			navg = c2 - c1 + 1
			z = buf + (j - 1) * nl
		    }
		    Memr[z+i-1] = asumr (Memr[y+c1-1], navg)
		}
	    }

	    # Trim the overscan vector and set the pixel coordinate.
	    call salloc (x, nl, TY_REAL)
	    call malloc (overscan, nsec*nl, TY_REAL)
	    y = overscan
	    do i = 1, nsec {
		if (BIAS_SEC(ccd) != NULL) {
		    c1 = BIAS_SC1(ccd,i)
		    c2 = BIAS_SC2(ccd,i)
		    l1 = BIAS_SL1(ccd,i)
		    l2 = BIAS_SL2(ccd,i)
		    navg = c2 - c1 + 1
		    npts = l2 - l1 + 1
		    y = overscan + (i - 1) * nl
		    z = buf + (i - 1) * nl
		}
		if (navg > 1)
		    call adivkr (Memr[z+l1-1], real (navg), Memr[z+l1-1],
			npts)
		call trim_overscan (Memr[z], npts, l1, Memr[x], Memr[y])
		call fit_overscan (Memc[str], c1, c2, l1, l2, Memr[x],
		    Memr[y], npts)
	    }

	} else {
	    first = l1
	    last = l2
	    navg = last - first + 1
	    npts = nc
	    call salloc (buf, npts, TY_REAL)
	    call aclrr (Memr[buf], npts)
	    do i = first, last
	        call aaddr (Memr[imgl2r(IN_IM(ccd),i)], Memr[buf], Memr[buf],
		    npts)
	    if (navg > 1)
	        call adivkr (Memr[buf], real (navg), Memr[buf], npts)

	    # Trim the overscan vector and set the pixel coordinate.
	    npts = CCD_C2(ccd) - CCD_C1(ccd) + 1
	    call malloc (overscan, npts, TY_REAL)
	    call salloc (x, npts, TY_REAL)
	    call trim_overscan (Memr[buf], npts, IN_C1(ccd), Memr[x],
		Memr[overscan])

	    call fit_overscan (Memc[str], c1, c2, l1, l2, Memr[x],
		Memr[overscan], npts)
	}
	
	# Set the CCD structure overscan parameters.
	CORS(ccd, OVERSCAN) = O
	COR(ccd) = YES
	OVERSCAN_VEC(ccd) = overscan

	# Log the operation.
	call strcpy ("overscan", Memc[errstr], SZ_LINE)
	y = overscan
	do i = 1, nsec {
	    if (BIAS_SEC(ccd) != NULL) {
		c1 = BIAS_SC1(ccd,i)
		c2 = BIAS_SC2(ccd,i)
		l1 = BIAS_SL1(ccd,i)
		l2 = BIAS_SL2(ccd,i)
		y = overscan + (i - 1) * nl
		npts = c2 - c1 + 1
		if (i > 1) {
		    call sprintf (Memc[errstr], SZ_LINE, "ovrscn%d")
			call pargi (i)
		}
	    }
	    call sprintf (Memc[str], SZ_LINE,
		"Overscan section is [%d:%d,%d:%d] with mean=%g")
		call pargi (c1)
		call pargi (c2)
		call pargi (l1)
		call pargi (l2)
		call pargr (asumr (Memr[y], npts) / npts)
	    call timelog (Memc[str], SZ_LINE)
	    call ccdlog (IN_IM(ccd), Memc[str])
	    call hdmpstr (OUT_IM(ccd), Memc[errstr], Memc[str])
	}

	call sfree (sp)
end


# FIT_OVERSCAN -- Fit a function to smooth the overscan vector.
#   The fitting uses the ICFIT procedures which may be interactive.
#   Changes to these parameters are "learned".  The user is queried with a four
#   valued logical query (XT_ANSWER routine) which may be turned off when
#   multiple images are processed.

procedure fit_overscan (image, c1, c2, l1, l2, x, overscan, npts)

char	image[ARB]		# Image name for query and title
int	c1, c2, l1, l2		# Overscan strip
real	x[npts]			# Pixel coordinates of overscan
real	overscan[npts]		# Input overscan and output fitted overscan
int	npts			# Number of data points

int	interactive, fd
pointer	sp, str, w, ic, cv, gp, gt

int	clgeti(), ic_geti(), open()
real	clgetr(), ic_getr()
pointer	gopen(), gt_init()
errchk	gopen, open

begin
	call smark (sp)
	call salloc (str, SZ_LINE, TY_CHAR)
	call salloc (w, npts, TY_REAL)
	call amovkr (1., Memr[w], npts)

	# Open the ICFIT procedures, get the fitting parameters, and
	# set the fitting limits.

	call ic_open (ic)
	call clgstr ("function", Memc[str], SZ_LINE)
	call ic_pstr (ic, "function", Memc[str])
	call ic_puti (ic, "order", clgeti ("order"))
	call clgstr ("sample", Memc[str], SZ_LINE)
	call ic_pstr (ic, "sample", Memc[str])
	call ic_puti (ic, "naverage", clgeti ("naverage"))
	call ic_puti (ic, "niterate", clgeti ("niterate"))
	call ic_putr (ic, "low", clgetr ("low_reject"))
	call ic_putr (ic, "high", clgetr ("high_reject"))
	call ic_putr (ic, "grow", clgetr ("grow"))
	call ic_putr (ic, "xmin", min (x[1], x[npts]))
	call ic_putr (ic, "xmax", max (x[1], x[npts]))
	call ic_pstr (ic, "xlabel", "Pixel")
	call ic_pstr (ic, "ylabel", "Overscan")

	# If the fitting is done interactively set the GTOOLS and GIO
	# pointers.  Also "learn" the fitting parameters since they may
	# be changed when fitting interactively.

	call sprintf (Memc[str], SZ_LINE,
	    "Fit overscan vector for %s[%d:%d,%d:%d] interactively")
	    call pargstr (image)
	    call pargi (c1)
	    call pargi (c2)
	    call pargi (l1)
	    call pargi (l2)
	call set_interactive (Memc[str], interactive)
	if ((interactive == YES) || (interactive == ALWAYSYES)) {
	    gt = gt_init ()
	    call sprintf (Memc[str], SZ_LINE,
	        "Overscan vector for %s from section [%d:%d,%d:%d]\n")
	        call pargstr (image)
		call pargi (c1)
		call pargi (c2)
		call pargi (l1)
		call pargi (l2)
	    call gt_sets (gt, GTTITLE, Memc[str])
	    call gt_sets (gt, GTTYPE, "line")
	    call gt_setr (gt, GTXMIN, x[1])
	    call gt_setr (gt, GTXMAX, x[npts])
	    call clgstr ("graphics", Memc[str], SZ_FNAME)
	    gp = gopen (Memc[str], NEW_FILE, STDGRAPH)

	    call icg_fit (ic, gp, "cursor", gt, cv, x, overscan, Memr[w], npts)

	    call ic_gstr (ic, "function", Memc[str], SZ_LINE)
	    call clpstr ("function", Memc[str])
	    call clputi ("order", ic_geti (ic, "order"))
	    call ic_gstr (ic, "sample", Memc[str], SZ_LINE)
	    call clpstr ("sample", Memc[str])
	    call clputi ("naverage", ic_geti (ic, "naverage"))
	    call clputi ("niterate", ic_geti (ic, "niterate"))
	    call clputr ("low_reject", ic_getr (ic, "low"))
	    call clputr ("high_reject", ic_getr (ic, "high"))
	    call clputr ("grow", ic_getr (ic, "grow"))

	    call gclose (gp)
	    call gt_free (gt)
	} else
	    call ic_fit (ic, cv, x, overscan, Memr[w], npts, YES, YES, YES, YES)

	# Make a log of the fit in the plot file if given.
	call clgstr ("plotfile", Memc[str], SZ_LINE)
	call xt_stripwhite (Memc[str])
	if (Memc[str] != EOS) {
	    fd = open (Memc[str], APPEND, BINARY_FILE)
	    gp = gopen ("stdvdm", NEW_FILE, fd)
	    gt = gt_init ()
	    call sprintf (Memc[str], SZ_LINE,
	        "Overscan vector for %s from section [%d:%d,%d:%d]\n")
	        call pargstr (image)
		call pargi (c1)
		call pargi (c2)
		call pargi (l1)
		call pargi (l2)
	    call gt_sets (gt, GTTITLE, Memc[str])
	    call gt_sets (gt, GTTYPE, "line")
	    call gt_setr (gt, GTXMIN, 1.)
	    call gt_setr (gt, GTXMAX, real (npts))
	    call icg_graphr (ic, gp, gt, cv, x, overscan, Memr[w], npts)
	    call gclose (gp)
	    call close (fd)
	    call gt_free (gt)
	}

	# Replace the raw overscan vector with the smooth fit.
	call cvvector (cv, x, overscan, npts)

	# Finish up.
	call ic_closer (ic)
	call cvfree (cv)
	call sfree (sp)
end


# TRIM_OVERSCAN -- Trim the overscan vector.

procedure trim_overscan (data, npts, start, x, overscan)

real	data[ARB]		# Full overscan vector
int	npts			# Length of trimmed vector
int	start			# Trim start
real	x[npts]			# Trimmed pixel coordinates (returned)
real	overscan[npts]		# Trimmed overscan vector (returned)

int	i, j

begin
	do i = 1, npts {
	    j = start + i - 1
	    x[i] = j
	    overscan[i] = data[j]
	}
end