include <imhdr.h>

# fl_corner -- merge a corner of external flat into internal flat
# This task takes an internal flat and an external flat and splices the
# upper left corner of the external flat to the internal, writing the
# output to a new image.  The splice region is defined by a pair of
# parallel diagonal lines.
# If the output name is null or equal to input1, then the output image
# will be renamed to the input1 image.
#
# Phil Hodge,  8-Mar-1990  Task created.
# Phil Hodge, 18-Sep-1990  Replace upper left corner instead of upper right.

procedure flcorner()

pointer input1			# scr for name of internal flat
pointer input2			# scr for name of external flat
pointer output			# scr for name of merged output image
real	slope			# slope of splice region
real	y1			# y intercept for lower limit of splice region
real	y2			# y intercept for upper limit of splice region
#--
pointer sp
pointer iim1, iim2		# imhdr ptrs for two input flat fields
pointer oim			# imhdr ptr for output image
pointer ix, ox			# pointers to data
pointer in1, in2, out		# scratch space for data
real	y1m, y2m		# y1, y2 modified to apply to corner box
real	x, y			# lower right corner of box containing splice
real	sum1, sum2		# sum of values within overlap region
real	ratio			# normalization factor:  internal / external
int	m1, m2			# size of corner including splice (m for mask)
int	naxis1, naxis2		# full size of flat fields
int	k
bool	inplace
pointer immap(), imgl2r(), impl2r(), imgs2r(), imps2r()
real	clgetr()
bool	streq()

begin
	call smark (sp)
	call salloc (input1, SZ_FNAME, TY_CHAR)
	call salloc (input2, SZ_FNAME, TY_CHAR)
	call salloc (output, SZ_FNAME, TY_CHAR)

	call clgstr ("input1", Memc[input1], SZ_FNAME)
	call clgstr ("input2", Memc[input2], SZ_FNAME)
	call clgstr ("output", Memc[output], SZ_FNAME)
	y1 = clgetr ("y1")
	y2 = clgetr ("y2")
	if (y1 > y2) {			# y1 should be smaller, so swap
	    slope = y1
	    y1 = y2
	    y2 = slope
	}
	slope = clgetr ("slope")
	if (slope == 0.)
	    call error (1, "slope can't be zero")

	if (streq (Memc[input1], Memc[output]) || (Memc[output] == EOS)) {
	    inplace = true
	    call mktemp ("flattemp", Memc[output], SZ_FNAME)
	} else {
	    inplace = false
	}

	# Open the input image and get its size.
	iim1 = immap (Memc[input1], READ_ONLY, NULL)
	naxis1 = IM_LEN(iim1,1)
	naxis2 = IM_LEN(iim1,2)

	# Find out how large a region is needed to contain the splice.
	x = (naxis2 - y1) / slope
	y = y1
	if (x < 0) {
	    call eprintf ("warning in flcorner:  x is negative\n")
	    x = naxis1
	}
	m1 = x + 1.			# add 1 to round up
	m2 = naxis2 - y + 4.		# extra space in y
	m1 = min (m1, naxis1)
	m1 = max (m1, 1)
	m2 = min (m2, naxis2)
	m2 = max (m2, 1)

	# Compute the y intercepts relative to the image subset.
	y1m = y1 - (naxis2 - m2)
	y2m = y2 - (naxis2 - m2)

	# Allocate memory for an image section for each of the input images
	# (internal & external flats) and for the output image.
	call salloc (in1, m1 * m2, TY_REAL)
	call salloc (in2, m1 * m2, TY_REAL)
	call salloc (out, m1 * m2, TY_REAL)

	# Copy the internal flat to the output image.
	oim = immap (Memc[output], NEW_COPY, iim1)
	do k = 1, naxis2 {
	    ix = imgl2r (iim1, k)
	    ox = impl2r (oim, k)
	    call amovr (Memr[ix], Memr[ox], naxis1)
	}
	# Close output.  We'll open it again later to replace the corner.
	call imunmap (oim)

	# Read the corner of the internal flat into memory.
	ix = imgs2r (iim1, 1, m1, naxis2-m2+1, naxis2)
	call amovr (Memr[ix], Memr[in1], m1 * m2)

	call imunmap (iim1)			# done with internal flat

	# Read the corner of the external flat into memory.
	iim2 = immap (Memc[input2], READ_ONLY, NULL)
	if (naxis1 != IM_LEN(iim2,1) || naxis2 != IM_LEN(iim2,2))
	    call error (1, "input1 and input2 are not the same size")
	ix = imgs2r (iim2, 1, m1, naxis2-m2+1, naxis2)
	call amovr (Memr[ix], Memr[in2], m1 * m2)
	call imunmap (iim2)			# done with external flat

	# Get the sum of pixel values in each input image within the overlap
	# region, as defined by the pair of lines.
	call fl_sum_inmask (Memr[in1], m1, m2, slope, y1m, y2m, sum1)
	call fl_sum_inmask (Memr[in2], m1, m2, slope, y1m, y2m, sum2)

	# Multiply external flat by ratio to normalize with internal flat.
	ratio = sum1 / sum2
	call amulkr (Memr[in2], ratio, Memr[in2], m1 * m2)

	# Splice in1 & in2 together, writing to Memr[out].
	call fl_splice (Memr[in1], Memr[in2], Memr[out], m1, m2,
			slope, y1m, y2m)
	
	# Insert output data into corner of output image.
	oim = immap (Memc[output], READ_WRITE, NULL)
	ox = imps2r (oim, 1, m1, naxis2-m2+1, naxis2)
	call amovr (Memr[out], Memr[ox], m1 * m2)
	call imunmap (oim)

	if (inplace) {
	    call imdelete (Memc[input1])
	    call imrename (Memc[output], Memc[input1])
	}

	call sfree (sp)
end

# fl_sum_inmask -- add pixel values within splice region

procedure fl_sum_inmask (input, m1, m2, slope, y1m, y2m, sum)

real	input[m1,m2]	# i: input data
int	m1, m2		# i: size of input array
real	slope		# i: slope of lines defining splice region
real	y1m		# i: y intercept for lower limit of splice region
real	y2m		# i: y intercept for upper limit of splice region
real	sum		# o: sum of values in splice region
#--
real	s		# local variable for sum
int	x, y		# loop indexes for x & coordinates within input data
int	ylow		# for a given x, lower limit of splice region
int	yhigh		# for a given x, upper limit of splice region

begin
	s = 0.
	do x = 1, m1 {
	    ylow = nint (y1m + slope * x)
	    yhigh = nint (y2m + slope * x)
	    ylow = max (ylow, 1)
	    yhigh = min (yhigh, m2)

	    if (ylow <= m2 && yhigh > 0) {
		do y = ylow, yhigh
		    s = s + input[x,y]
	    }
	}

	sum = s
end


# fl_splice -- splice two input into output
# This routine takes two rectangular input regions and splices them together
# to form the output.  The splice region is defined by two lines with the
# same slope but different y intercepts.  Below the lower line, in1 is copied
# to the output; above the upper line, in2 is copied.  The lower and upper
# lines cross a given x column at ylow & yhigh.  The nearest integers to
# those points are the limits (inclusive) of the splice region.  The weights
# for combining in1 and in2 are linear functions of the vertical distance
# of the pixel from ylow-1; the weight for in1 is one at ylow-1 and zero
# at yhigh+1.

procedure fl_splice (in1, in2, out, m1, m2,
			slope, y1m, y2m)

real	in1[m1,m2]	# i: input data for internal flat
real	in2[m1,m2]	# i: input data for external flat
real	out[m1,m2]	# o: in1 & in2 spliced together
int	m1, m2		# i: size of arrays
real	slope		# i: slope of lines defining splice region
real	y1m		# i: y intercept for lower limit of splice region
real	y2m		# i: y intercept for upper limit of splice region
#--
real	width		# vertical width of splice region
real	w		# weight factor for splicing
int	x, y		# loop indexes for x & coordinates within input data
int	ylow		# for a given x, lower limit of splice region
int	yhigh		# for a given x, upper limit of splice region

begin
	width = (y2m - y1m) + 1.

	do x = 1, m1 {

	    ylow = nint (y1m + slope * x)
	    yhigh = nint (y2m + slope * x)

	    ylow = max (ylow, 1)
	    yhigh = min (yhigh, m2)

	    # Fill lower region with in1.
	    do y = 1, min (ylow-1, m2)
		out[x,y] = in1[x,y]

	    # Fill upper region with in2.
	    do y = max (yhigh+1, 1), m2
		out[x,y] = in2[x,y]

	    # Fill splice region.
	    if (ylow <= m2 && yhigh > 0) {
		do y = ylow, yhigh {
		    w = real (y - (ylow-1)) / (width + 1.)
		    out[x,y] = w * in1[x,y] + (1.-w) * in2[x,y]
		}
	    }
	}
end