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

define	SZ_KEYVAL	11	# size of string for a keyword value
define	MAX_NSAMPS    1024	# width of photocathode

# Evaluate a polynomial fit file at each pixel of an image, writing the
# x,y coordinates to an output image which will serve as input to calfoc
# using the newgeom algorithm.  This geometric correction file will have
# width 2 * (npix1 + 1) and height (npix2 + 1), where npix1 & npix2
# are the axis lengths of the image to be calibrated.
#
# As of 9-Feb-1994, the polynomial fit should be computed in photocathode
# coordinates; this is the coordinate system that is used by focgeom tasks
# that read or write reseau positions.  The GEO file written by this task
# still uses image pixel coordinates rather than photocathode coordinates.
# Prior to this change, the observed and reference positions had to be
# shifted by subtracting (flipped) sampoff and lineoff before running rfitx
# to do the fit.  
#
# Phil Hodge,    Jan-1991  Extracted from newgeom task.
# Phil Hodge,  8-Aug-1991  Get & put coordinate parameters.
# Phil Hodge, 22-Jan-1993  Use IS_INDEF instead IS_INDEFR (of).
# Phil Hodge,  9-Feb-1994  Add and then subtract sampoff & lineoff;
#			also get pedigree from cl and write to output header.
# Phil Hodge, 12-Apr-1994  Get imscale from cl and write to output header;
#			don't get (or put) c_valid; change length of pedigree
#			string from SZ_FNAME to SZ_PATHNAME.

procedure t_p2geo()

char	inpoly[SZ_FNAME]	# polynomial database file
char	output[SZ_FNAME]	# output image (geo file)
char	entry[SZ_R_PATTERN]	# transform entry to use
char	pedigree[SZ_PATHNAME]	# value for pedigree keyword to add to header
char	template[SZ_FNAME]	# template image
#--
pointer	tim, oim		# pointers to imhdr struct
pointer	out			# pointer to image data
pointer dt
pointer	fit
pointer	sx1, sx2, sy1, sy2
pointer	xout, yout
pointer immap()
pointer	impl2r()
pointer dtmap()
pointer sp

real	x, y
real	xoff, yoff		# sampoff, lineoff (if centered format)
int	i, j, npix1, npix2
int	dirn			# Direction number,1=f/w,2=b/w
int	axis
int	rec			# record number of the entry in database
int	tenth			# print message after this many lines
bool	verbose
int	db_readr()
double	clgetd()
real	clgetr()
int	clgeti()
bool	clgetb()

char	optcrly[SZ_KEYVAL]		# F48 or F96
char	pxformt[SZ_KEYVAL]		# NORMAL or ZOOM
char	cammode[SZ_KEYVAL]		# INBEAM or NOTUSED
char	smmmode[SZ_KEYVAL]		# INBEAM or NOTUSED
real	sampoff, lineoff
int	samppln, linepfm
int	optelt[2]
double	imscale
char	ctype1[SZ_KEYVAL]
char	ctype2[SZ_KEYVAL]
double	crval[2]
real	crpix[2]
real	cd[2,2]

begin
	call smark (sp)

	# 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")

	# Get the string to add to the output header with keyword 'pedigree'.
	call clgstr ("pedigree", pedigree, SZ_PATHNAME)

	# Get name of template image.  This defines the output header.
	call clgstr ("template", template, SZ_FNAME)

	# This is the size of the image which will be created when
	# newgeom is run using the geo file that we're creating.
	npix1 = clgeti ("npix1")
	npix2 = clgeti ("npix2")
	if (IS_INDEFI (npix1) || IS_INDEFI (npix2))
	    call error (1, "size of image must be specified")

	# Get the other parameters that will be included in the header
	# of the output image.
	# Note that pxformt, samppln, sampoff and lineoff, are not just
	# copied to the output header, they are used to find the appropriate
	# offset between (flipped) photocathode coordinates and pixel coords.
	call clgstr ("optcrly", optcrly, SZ_KEYVAL)
	call clgstr ("pxformt", pxformt, SZ_KEYVAL)
	call clgstr ("cammode", cammode, SZ_KEYVAL)
	call clgstr ("smmmode", smmmode, SZ_KEYVAL)
	sampoff = clgetr ("sampoff")
	lineoff = clgetr ("lineoff")
	samppln = clgeti ("samppln")
	linepfm = clgeti ("linepfm")
	optelt[1] = clgeti ("optelt1")
	optelt[2] = clgeti ("optelt2")

	imscale = clgetd ("imscale")
	crval[1] = clgetd ("crval1")
	crval[2] = clgetd ("crval2")
	crpix[1] = clgetr ("crpix1")
	crpix[2] = clgetr ("crpix2")
	cd[1,1] = clgetr ("cd1_1")
	cd[1,2] = clgetr ("cd1_2")
	cd[2,1] = clgetr ("cd2_1")
	cd[2,2] = clgetr ("cd2_2")
	call clgstr ("ctype1", ctype1, SZ_KEYVAL)
	call clgstr ("ctype2", ctype2, SZ_KEYVAL)

	verbose = clgetb ("verbose")

	call strupr (optcrly)
	call strupr (pxformt)
	call strupr (cammode)
	call strupr (smmmode)

	# These are the offsets from image pixel coordinates to
	# photocathode coordinates.  Well, sort of.  The difference
	# is that we're using coordinates that are flipped around
	# the Y axis.  xoff and sampoff will be the same unless the
	# format is not centered.
	if (pxformt[1] == 'Z')				# zoom format
	    xoff = MAX_NSAMPS - 2*samppln - sampoff	# "flip" sampoff
	else						# normal format
	    xoff = MAX_NSAMPS - samppln - sampoff
	yoff = lineoff

	# Allocate space for fit pointers
	call salloc (fit, LEN_FIT, TY_STRUCT)

	# Allocate space for x and y coordinates of corners
	call salloc (xout, npix1+1, TY_REAL)
	call salloc (yout, npix1+1, TY_REAL)

	# Set the direction of the transformation
	dirn = FORWARD

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

	# Open template image.
	tim = immap (template, READ_ONLY, NULL)

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

	# Set parameters for output file
	IM_LEN (oim, 1) = 2 * (npix1 + 1)
	IM_LEN (oim, 2) = npix2 + 1
	IM_PIXTYPE (oim) = TY_REAL

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

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

	# Evaluate the polynomial at the corner of each pixel.
	# Note that we offset to get (flipped) photocathode coordinates.
	do j = 0, npix2 {
	    y = real (j) + yoff
	    do i = 0, npix1 {
		x = real (i) + xoff
		call polyeval (x, y, sx1, sx2, sy1, sy2,
			Memr[xout+i], Memr[yout+i])
	    }
	    # Subtract xoff & yoff to get back to image pixel coordinates.
	    call asubkr (Memr[xout], xoff, Memr[xout], npix1+1)
	    call asubkr (Memr[yout], yoff, Memr[yout], npix1+1)

	    # Put the values into the output image.
	    out = impl2r (oim, j+1)		# line number is one indexed
	    call amovr (Memr[xout], Memr[out], npix1+1)
	    call amovr (Memr[yout], Memr[out+npix1+1], npix1+1)

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

	# Set header parameters.
	call add_headname (oim, output, INDEF, INDEF)

	call imastr (oim, "optcrly", optcrly)
	call imastr (oim, "pxformt", pxformt)
	call imastr (oim, "cammode", cammode)
	call imastr (oim, "smmmode", smmmode)
	call imaddr (oim, "sampoff", sampoff)
	call imaddr (oim, "lineoff", lineoff)
	call imaddi (oim, "samppln", samppln)
	call imaddi (oim, "linepfm", linepfm)
	call imaddi (oim, "optelt1", optelt[1])
	call imaddi (oim, "optelt2", optelt[2])

	call imaddi (oim, "npix1", npix1)
	call imaddi (oim, "npix2", npix2)

	call imaddd (oim, "imscale", imscale)
	call imaddd (oim, "crval1", crval[1])
	call imaddd (oim, "crval2", crval[2])
	call imaddr (oim, "crpix1", crpix[1])
	call imaddr (oim, "crpix2", crpix[2])
	call imaddr (oim, "cd1_1", cd[1,1])
	call imaddr (oim, "cd1_2", cd[1,2])
	call imaddr (oim, "cd2_1", cd[2,1])
	call imaddr (oim, "cd2_2", cd[2,2])
	call imastr (oim, "ctype1", ctype1)
	call imastr (oim, "ctype2", ctype2)

	call imastr (oim, "pedigree", pedigree)

	call imputh (oim, "history", inpoly)

	call imunmap (oim)
	call imunmap (tim)

	# Close the database file.
	call dtunmap (dt)

	call sfree (sp)
end

# POLYEVAL -- evaluate the polynomial

procedure polyeval (x, y, sx1, sx2, sy1, sy2, x_sol, y_sol)

pointer sx1, sx2, sy1, sy2
real	x_sol, y_sol
real	x, y
real	gseval()

begin
	x_sol = gseval (sx1, x, y)
	y_sol = gseval (sy1, x, y)
	if (sx2 != NULL)
	    x_sol = x_sol + gseval (sx2, x, y)
	if (sy2 != NULL)
	    y_sol = y_sol + gseval (sy2, x, y)
end