include <time.h>
include <pkg/dttext.h>
include	<imhdr.h>
include	<foc.h>
include "../dbfit.h"

define	SZ_KEYVAL	11	# size of string for a keyword value

# geoevalfit -- apply polynomial fit to geo correction file
# This task reads a geometric correction file (image), evaluates a
# polynomial fit at each pixel, and writes the results to an output
# geometric correction file.
#
# Phil Hodge,  7-Apr-1992  Task created based on p2geo.
# Phil Hodge,  9-Apr-1993  Remove call to flush after printing "percent done:"

procedure geoevalfit()

char	input[SZ_FNAME]		# input image (geo file)
char	inpoly[SZ_FNAME]	# polynomial database file
char	output[SZ_FNAME]	# output image (geo file)
char	entry[SZ_R_PATTERN]	# transform entry to use
bool	offset			# apply sampoff & lineoff offsets?
bool	verbose			# print percent done?
#--
pointer	iim, oim		# pointers to imhdr struct for input & output
pointer	iv, ov			# pointers to input & output data
pointer dt
pointer	sx1, sx2, sy1, sy2
pointer immap()
pointer	imgl2r(), impl2r()
pointer dtmap()
pointer sp
pointer	fit			# scratch for fit struct
pointer message			# scratch for history info for output image

char	timestr[SZ_TIME]	# date & time
long	zero			# used by clktime
long	ctime			# returned by clktime
long	clktime()

real	sampoff, lineoff	# offsets between X & Y and photocathode coords
real	x, y			# X & Y coords from input image at (i,j)
real	xout, yout		# x, y after applying polynomial function
int	i, j			# loop indexes
int	nx, ny			# lengths of X & Y axes of images
int	halfnx			# nx / 2
int	direction		# direction number for evaluating fit
int	axis
int	tenth			# print message after this many lines
real	imgetr()
pointer db_template(), dlist
int	db_list(), db_readr()
int	strlen()
bool	clgetb()

begin
	# Allocate space for fit pointers
	call smark (sp)
	call salloc (fit, LEN_FIT, TY_STRUCT)
	call salloc (message, SZ_FNAME, TY_CHAR)

	# Get name of input image (geo file).
	call clgstr ("input", input, SZ_FNAME)

	# Get name of polynomial coefficients file, and map it.
	call clgstr ("inpoly", inpoly, SZ_FNAME)
	dt = dtmap (inpoly, READ_ONLY)

	# Get name of entry to process
	call clgstr ("entry", entry, SZ_R_PATTERN)

	# Get name of geo file to create.
	call clgstr ("output", output, SZ_FNAME)
	if (output[1] == EOS)
	    call error (1, "output file name must be given")

	# This can be set to either FORWARD or BACKWARD with the same
	# results because the evaluation is only done forward.
	direction = FORWARD

	offset = clgetb ("offset")
	verbose = clgetb ("verbose")

	# Open input image.
	iim = immap (input, READ_ONLY, NULL)

	# Create output image.
	oim = immap (output, NEW_COPY, iim)

	# Get the size of the input image.
	nx = IM_LEN(iim,1)
	ny = IM_LEN(iim,2)
	halfnx = nx / 2				# should divide exactly

	if (offset) {
	    sampoff = imgetr (iim, "sampoff")
	    lineoff = imgetr (iim, "lineoff")
	} else {
	    sampoff = 0.
	    lineoff = 0.
	}

	# Read the polynomial parameters into the fit struct.
	dlist = db_template (dt, entry)
	if (db_readr (dt, direction, "", db_list(dlist), axis, fit) ==
			RES_F_NONEXTENTRY)
	    call error (1, "entry not found")

	sx1 = FIT_SX1(fit)
	sy1 = FIT_SY1(fit)
	sx2 = FIT_SX2(fit)
	sy2 = FIT_SY2(fit)

	if (verbose)
	    call printf ("percent done: ")
	tenth = (ny + 9) / 10

	# Evaluate the polynomial at each pixel.
	do j = 0, ny-1 {			# zero indexed
	    iv = imgl2r (iim, j+1)		# line number is one indexed
	    ov = impl2r (oim, j+1)
	    do i = 0, halfnx-1 {		# zero indexed
		x = Memr[iv+i] + sampoff
		y = Memr[iv+halfnx+i] + lineoff

		call polyeval (x, y, sx1, sx2, sy1, sy2, xout, yout)

		Memr[ov+i] = xout - sampoff
		Memr[ov+halfnx+i] = yout - lineoff
	    }

	    if (verbose) {
		# Print a message every ten percent of the way through.
		if (mod (j+1, tenth) == 0 && j < ny-2) {
		    call printf (" %d")
			call pargi ((j+1) / tenth * 10)
		    call flush (STDOUT)
		}
	    }
	}
	if (verbose) {
	    call printf (" 100\n")
	    call flush (STDOUT)
	}

	# Add time & date and names of input and poly files as history info.
	zero = 0
	ctime = clktime (zero)
	call cnvtime (ctime, timestr, SZ_TIME)
	call strcpy ("Transformed by geoeval on ", Memc[message], SZ_FNAME)
	call strcat (timestr, Memc[message], SZ_FNAME)
	call strcat (".", Memc[message], SZ_FNAME)
	call imputh (oim, "history", Memc[message])

	# 27 is the length of the string "geoeval.input (geo file) = ".
	if (strlen (input) > SZ_FNAME-27) {
	    call strcpy (input, Memc[message], SZ_FNAME)
	} else {
	    call strcpy ("geoeval.input (geo file) = ",
			Memc[message], SZ_FNAME)
	    call strcat (input, Memc[message], SZ_FNAME)
	}
	call imputh (oim, "history", Memc[message])

	# 28 is the length of the string "geoeval.inpoly (fit file) = ".
	if (strlen (inpoly) > SZ_FNAME-28) {
	    call strcpy (inpoly, Memc[message], SZ_FNAME)
	} else {
	    call strcpy ("geoeval.inpoly (fit file) = ",
			Memc[message], SZ_FNAME)
	    call strcat (inpoly, Memc[message], SZ_FNAME)
	}
	call imputh (oim, "history", Memc[message])

	call imunmap (oim)
	call imunmap (iim)

	# Close the database file.
	call dtunmap (dt)

	call sfree (sp)
end