include <imhdr.h>

define	SZ_KEYW		11	# size of buffer for header keyword
define	FEX_NONE	1	# no background subtraction
define	FEX_MEAN	2	# subtract mean
define	FEX_SIGCLIP	3	# subtract mean, clipping by n sigma
define	NUM_STARS	20	# size & increment for arrays of star coords

# fextract -- extract stars to another image
# Given a list of images and a list of star coordinates, 
# Given an input image, an output image, and a list of star coordinates
# in pixels, this task copies a region around each star on the input image
# to the output image.  If the output image does not exist it will be
# created.  The stars for the current input image are copied to a
# horizontal strip of the output; the strips for different input images
# are stacked vertically.  The input images are at different focus positions,
# so the focus changes vertically in the output image.  A global value for
# the background may be subtracted, and the current output strip will be
# divided by the exposure time for the current input image.
#
# Phil Hodge, 29-Jan-1990  Task created.
# Phil Hodge, 23-Jan-1992  Divide by exposure time.

procedure fextract()

pointer input			# scr for name of the input image
pointer output			# scr for name of the output image
pointer coords			# scr for name of file containing star coords
pointer option			# scr for background-subtraction option
int	refoc			# position of focussing mirror
int	maximages		# allow for this many input images
int	rownum			# position of input image in output image
int	box[2]			# size of box to extract around each star
int	b_option		# background-subtraction option
int	lowcut, highcut		# low & high cutoffs for computing mean
real	nsigma			# clip values outside nsigma * sigma
#--
pointer sp
pointer x, y			# scr for arrays of coordinates of stars
pointer iim, oim		# input, output imio pointers
char	keyword[SZ_KEYW]	# keyword for writing to output header
real	background		# background level in an input image
real	exptime			# exposure time for input image
int	numstars		# number of stars per image
int	j			# loop index for star
pointer immap()
real	clgetr(), imgetr()
int	imaccess(), imgeti(), clgwrd(), clgeti()

begin
	call smark (sp)
	call salloc (input, SZ_FNAME, TY_CHAR)
	call salloc (output, SZ_FNAME, TY_CHAR)
	call salloc (coords, SZ_FNAME, TY_CHAR)
	call salloc (option, SZ_FNAME, TY_CHAR)

	call clgstr ("input", Memc[input], SZ_FNAME)
	call clgstr ("output", Memc[output], SZ_FNAME)
	refoc = clgeti ("refoc")
	call clgstr ("coords", Memc[coords], SZ_FNAME)

	rownum = clgeti ("rownum")

	# Open the input image.
	iim = immap (Memc[input], READ_ONLY, NULL)

	# Get exposure time for this image.
	exptime = imgetr (iim, "exptime")
	if (exptime <= 0) {
	    call eprintf ("exptime = %f; give new value\n")
		call pargr (exptime)
	    exptime = clgetr ("exptime")
	}

	# Get star coordinates; x & y are allocated by this routine.
	call fex_xypairs (Memc[coords], numstars, x, y)

	# Check whether the output image already exists.
	if (imaccess (Memc[output], 0) == YES) {

	    oim = immap (Memc[output], READ_WRITE, NULL)
	    box[1] = imgeti (oim, "box1")
	    box[2] = imgeti (oim, "box2")

	} else {		# we must create the output image

	    maximages = clgeti ("maximages")
	    box[1] = clgeti ("box1")
	    box[2] = clgeti ("box2")

	    # Create the output image.
	    call fex_output (iim, oim, Memc[output],
			box, numstars, maximages)
	}

	# Replace rownum if it was not specified, and add header parameters.
	call fex_header (Memc[input], oim, rownum, refoc,
			Memi[x], Memi[y], numstars)

	b_option = clgwrd ("b_option", Memc[option], SZ_LINE,
			"|none|mean|sigclip|")
	lowcut = 0			# just to have definite values
	highcut = 10
	nsigma = 3.
	if (b_option == FEX_MEAN) {
	    lowcut = clgeti ("lowcut")
	    highcut = clgeti ("highcut")
	} else if (b_option == FEX_SIGCLIP) {
	    nsigma = clgetr ("nsigma")
	}

	# Find out what is the background level.
	call fex_background (iim, b_option, lowcut, highcut, nsigma,
			background)

	# Add background level to output header.
	call sprintf (keyword, SZ_KEYW, "backg%d")
		    call pargi (rownum)
	call imaddr (oim, keyword, background)

	# Copy a region around each star to the output image, subtracting
	# the background.
	do j = 1, numstars {
	    call fex_copy (iim, oim, j, rownum, Memi[x+j-1], Memi[y+j-1],
			box, background, exptime)
	}

	call imunmap (iim)
	call imunmap (oim)

	if (x != NULL)
	    call mfree (x, TY_INT)
	if (y != NULL)
	    call mfree (y, TY_INT)

	call sfree (sp)
end

# fex_xypairs -- get star coordinates
# This routine opens an ascii file, allocates pointers to int arrays,
# reads coordinate pairs into those arrays, and closes the file.

procedure fex_xypairs (coords, numstars, x, y)

char	coords[ARB]	# i: name of file of coordinates
int	numstars	# o: number of stars
pointer x		# o: pointer to x coordinates of stars
pointer y		# o: pointer to y coordinates of stars
#--
char	linebuf[SZ_FNAME]	# input buffer
int	fd			# fd for file of coordinates
int	j
int	maxstars		# size of x & y arrays
int	ncoords			# number of coords read from line
int	dummy
int	open(), getline(), sscan(), nscan()

begin
	# Open the file of coordinates.
	fd = open (coords, READ_ONLY, TEXT_FILE)

	# Allocate x & y for coordinates.
	maxstars = NUM_STARS			# initial value
	call malloc (x, maxstars, TY_INT)
	call malloc (y, maxstars, TY_INT)

	j = 0
	while (getline (fd, linebuf) != EOF) {
	    if (linebuf[1] == '#')
		next
	    dummy = sscan (linebuf)
		call gargi (Memi[x+j])		# zero indexed
		call gargi (Memi[y+j])
	    ncoords = nscan()

	    if (ncoords == 1)
		call error (1, "must give both x & y coordinates of each star")
	    else if (ncoords < 1)
		next				# blank line or comment?

	    j = j + 1
	    if (j >= maxstars) {
		maxstars = maxstars + NUM_STARS
		call realloc (x, maxstars, TY_INT)
		call realloc (y, maxstars, TY_INT)
	    }
	}

	if (j < 1)
	    call error (1, "no coordinates specified")

	numstars = j

	call close (fd)
end


# fex_output -- create output image
# This routine creates the output image as a new copy of the (first) input
# image except that the size is specified depending on the box size and
# the number of stars and total number of input images.  Keywords are added
# to the header giving the box size and star coordinates.  The data values
# are initialized to zero.

procedure fex_output (iim, oim, output,
		box, numstars, maximages)

pointer iim		# i: pointer to imio struct for input image
pointer oim		# o: pointer to imio struct for output image
char	output[ARB]	# i: name of output image
int	box[2]		# i: size of box for each star
int	numstars	# i: needed for length of horizontal axis
int	maximages	# i: needed for length of vertical axis
#--
pointer x		# points to data values
int	naxis1, naxis2	# size of output image
int	k		# loop index
pointer immap()
pointer impl2r()

begin
	oim = immap (output, NEW_COPY, iim)		# open output

	# The output will have different stars in the horizontal direction
	# and different focus positions (input images) vertically.  Set the
	# data type to real.
	naxis1 = numstars * box[1]
	naxis2 = maximages * box[2]
	IM_NDIM(oim) = 2
	IM_LEN(oim,1) = naxis1
	IM_LEN(oim,2) = naxis2
	IM_PIXTYPE(oim) = TY_REAL

	# Add box size to output image.
	call imaddi (oim, "box1", box[1])
	call imaddi (oim, "box2", box[2])

	# Initialize the image to zero.
	do k = 1, naxis2 {
	    x = impl2r (oim, k)
	    call amovkr (0., Memr[x], naxis1)
	}
end


# fex_header -- update rownum, add header parameters
# If the row number was not specified, this routine sets it to one more than
# the largest rownum that has been used so far.  The value of rownum and the
# coordinates of the stars are added to the output header.

procedure fex_header (input, oim, rownum, refoc,
		x, y, numstars)

char	input[ARB]	# i: name of input image
pointer oim		# i: pointer to imio struct for output image
int	rownum		# io: position of input stars in output image
int	refoc		# i: current refoc mirror position
int	x[numstars]	# i: x coordinates of stars
int	y[numstars]	# i: y coordinates of stars
int	numstars	# i: number of stars
#--
char	rowstr[SZ_FNAME]	# string for row numbers
char	keyword[SZ_KEYW]	# keyword for writing to output header
int	j			# loop index
int	maxrow			# largest row number used so far
int	ip
int	imaccf(), ctoi(), clgeti()

begin
	# Update the rownum string in the output header, and possibly
	# update the value of rownum.

	if (imaccf (oim, "rownum") == NO) {

	    if (rownum < 1)
		rownum = 1		# set the value of rownum
	    call sprintf (rowstr, SZ_FNAME, "%d")
		call pargi (rownum)
	    call imastr (oim, "rownum", rowstr)

	} else {

	    call imgstr (oim, "rownum", rowstr, SZ_FNAME)

	    if (rownum >= 1) {
		# Check whether this rownum has been used before.  If so,
		# copy current value to newrow parameter, and get newrow.
		ip = 1
		while (ctoi (rowstr, ip, j) > 0) {
		    if (j == rownum) {
			call eprintf (
	"Warning:  This rownum has been used before.  Please verify value.\n")
			call eprintf ("    (set to -1 for default)\n")
			call clputi ("newrow", rownum)
			rownum = clgeti ("newrow")
			break
		    }
		}
	    }

	    if (rownum < 1) {
		# Set rownum to maxrow + 1.
		maxrow = 0
		ip = 1
		while (ctoi (rowstr, ip, j) > 0) {
		    if (j > maxrow)
			maxrow = j
		}
		rownum = maxrow + 1
	    }

	    # Add rownum to the rownum string, and put it in the header.
	    # (keyword is not really a keyword here, just a short string)
	    call sprintf (keyword, SZ_KEYW, " %d")
		call pargi (rownum)
	    call strcat (keyword, rowstr, SZ_FNAME)
	    call imastr (oim, "rownum", rowstr)
	}

	# Add the name of the input image to the output header.
	call sprintf (keyword, SZ_KEYW, "input%d")
	    call pargi (rownum)
	call imastr (oim, keyword, input)

	# Add the current refoc position to the header.
	call sprintf (keyword, SZ_KEYW, "refoc%d")
	    call pargi (rownum)
	call imaddi (oim, keyword, refoc)

	# Keyword x3_5 means the x coordinate of the third star extracted
	# from the fifth input image.  3 & 5 are the (x,y) box numbers in
	# the output image.
	do j = 1, numstars {
	    call sprintf (keyword, SZ_KEYW, "x%d_%d")
		call pargi (j)
		call pargi (rownum)
	    call imaddi (oim, keyword, x[j])
	    call sprintf (keyword, SZ_KEYW, "y%d_%d")
		call pargi (j)
		call pargi (rownum)
	    call imaddi (oim, keyword, y[j])
	}
end


# fex_background -- get the background level in an image
# The lowcut & highcut parameters are used for the FEX_MEAN option,
# while nsigma is used for FEX_SIGCLIP.

procedure fex_background (iim, b_option, lowcut, highcut, nsigma,
		background)

pointer iim		# i: pointer to imio struct for input image
int	b_option	# i: background-subtraction option
int	lowcut		# i: values below this are ignored
int	highcut		# i: values above this are ignored
real	nsigma		# i: clip values farther than nsigma * sigma from mean
real	background	# o: background level in image
#--
real	mean		# mean of pixel values
real	sigma		# sqrt of mean used as approximation of std dev
int	low, high	# local values of lowcut, highcut based on sigma

begin
	if (b_option == FEX_NONE) {

	    background = 0.		# don't subtract any background

	} else if (b_option == FEX_MEAN) {

	    # background is the mean of all the values within [lowcut, highcut].
	    call fex_gmean (iim, lowcut, highcut, background)

	} else if (b_option == FEX_SIGCLIP) {

	    # Get the mean of all values.
	    call fex_gmean (iim, INDEFI, INDEFI, mean)

	    # Use the approximation that the mean is equal to the variance.
	    sigma = sqrt (mean)
	    low = nint (mean - nsigma * sigma)
	    high = nint (mean + nsigma * sigma)

	    # Get the mean again, clipping beyond [low, high].
	    call fex_gmean (iim, low, high, background)
	}
end


# fex_gmean -- get the mean value in an image
# If either lowcut or highcut is INDEF, all the values in the image will
# be averaged (i.e. no clipping).  If both have definite values, the mean
# will be computed using only values within [lowcut, highcut] inclusive.

procedure fex_gmean (iim, lowcut, highcut, mean)

pointer iim		# i: pointer to imio struct for input image
int	lowcut		# i: values below this are ignored
int	highcut		# i: values above this are ignored
real	mean		# o: mean value in image
#--
pointer x		# pointer to input data
int	naxis1, naxis2	# size of input image
int	j, k		# loop indexes
int	n		# number of points within lowcut, highcut
real	value		# a pixel value
real	sum		# sum of pixel values
pointer imgl2r()

begin
	naxis1 = IM_LEN(iim,1)
	naxis2 = IM_LEN(iim,2)

	n = 0
	sum = 0.

	if (IS_INDEFI(lowcut) || IS_INDEFI(highcut)) {

	    # Average all values.
	    do k = 1, naxis2 {
		x = imgl2r (iim, k)
		do j = 0, naxis1-1 {
		    value = Memr[x+j]
		    n = n + 1
		    sum = sum + value
		}
	    }

	    if (n <= 0)
		call error (1, "fex_gmean:  image section is null")

	} else {

	    # Include only values within range when computing average.
	    do k = 1, naxis2 {
		x = imgl2r (iim, k)
		do j = 0, naxis1-1 {
		    value = Memr[x+j]
		    if (value >= lowcut && value <= highcut) {
			n = n + 1
			sum = sum + value
		    }
		}
	    }

	    if (n <= 0)
		call error (1, "fex_gmean:  no values within lowcut, highcut")
	}

	mean = sum / n
end


# fex_copy -- copy a part of an image to a part of another image
# A rectangular section is copied from one image to another.  The location
# of the section in the input may be changed in order to keep the section
# entirely within the input image.  If lowcut is not indef, pixel values
# below lowcut will be set to lowcut.  Background will be subtracted, and
# the image section will be normalized by dividing by the exposure time.

procedure fex_copy (iim, oim, j, k, x, y,
		box, background, exptime)

pointer iim			# i: pointer to imio struct for input image
pointer oim			# i: pointer to imio struct for output image
int	j, k			# i: location of output (unit = boxes)
int	x, y			# i: pixel location of star in input
int	box[2]			# i: size of box to copy (unit = pixels)
real	background		# i: background level to subtract
real	exptime			# i: exposure time for current image
#--
pointer ix			# pointer to input section
pointer ox			# pointer to output section
int	hbox[2]			# half of box
int	naxis[2]		# size of input image
int	ij1, ij2, ik1, ik2	# image section in input image
int	oj1, oj2, ok1, ok2	# image section in output image
pointer imgs2r(), imps2r()

begin
	hbox[1] = box[1] / 2
	hbox[2] = box[2] / 2
	naxis[1] = IM_LEN(iim,1)
	naxis[2] = IM_LEN(iim,2)

	# Get the limits of the image section in the input image.
	ij1 = x - hbox[1]
	ik1 = y - hbox[2]

	ij1 = max (1, ij1)
	ik1 = max (1, ik1)

	ij2 = ij1 + box[1] - 1
	ik2 = ik1 + box[2] - 1

	ij2 = min (ij2, naxis[1])
	ik2 = min (ik2, naxis[2])

	ij1 = ij2 - box[1] + 1
	ik1 = ik2 - box[1] + 1
	if (ij1 < 1)
	    call error (1, "box is wider than image")
	if (ik1 < 1)
	    call error (1, "box is taller than image")

	# Get the limits of the image section in the output image.
	oj1 = (j-1) * box[1] + 1
	oj2 = j * box[1]
	ok1 = (k-1) * box[2] + 1
	ok2 = k * box[2]

	ix = imgs2r (iim, ij1, ij2, ik1, ik2)
	ox = imps2r (oim, oj1, oj2, ok1, ok2)

	# Copy the image section.
	call amovr (Memr[ix], Memr[ox], box[1]*box[2])

	# Subtract background.
	if (background != 0.)
	    call asubkr (Memr[ox], background, Memr[ox], box[1]*box[2])

	# Normalize by exposure time.
	call adivkr (Memr[ox], exptime, Memr[ox], box[1]*box[2])
end