procedure imregister (input,output,refimage)

string	input	=""		{prompt="Input image"}
string	output	=""		{prompt="Output image"}
string	refimage=""		{prompt="Reference image"}
string	coords	=""		{prompt="Reference coordinates file"}
int	cntrbox	=7		{prompt="Size of the fine centering box"}
int	srchbox	=11		{prompt="Size of the coarse centering box"}
int	niterate=3		{prompt="Maximum number of iterations"}
bool	refrecen=no		{prompt="Recenter reference coordinates?"}
bool	execute	=yes		{prompt="Perform the image transformation?"}
string	interpol="spline3"	{prompt="Interpolation type (near|lin|poly3|poly5|spline3|sinc)",
				enum="nearest|linear|poly3|poly5|spline3|sinc"}
string	boundary="wrap"		{prompt="Boundary type (near|reflect|wrap)",
				enum="nearest|reflect|wrap"}
bool	verbose	=yes		{prompt="Print centers, and shifts?"}
string	logfile ="register.log"	{prompt="Text logfile name"}

# @(#)task imregister                 Author: R.A. Jansen -- Jun 19 2000
# @(#)
# @(#)Task to register (shift,scale and rotate) an image using a list of
# @(#)matched input coordinates.  The coordinate file need not be given;
# @(#)if absent, the reference image is displayed and reference objects
# @(#)can be marked interactively.
# @(#)Uses: tvmark, geomap and geotran, and several user written tasks.

begin
	string	img, refimg, coofil, outimg, intpl, bndry, txtlog
	string	tmptvmark, tmpgeo, tmpsol, cstr
	int	cbox, sbox, ncol, nrow, NN, N1, N2
	real	x, y, x1, y1, x2, y2
	bool	oldcoo

	img	= input
	refimg	= refimage
	outimg	= output
	coofil	= coords
	cbox	= cntrbox
	sbox	= srchbox
	intpl	= interpol
	bndry	= boundary
	txtlog	= logfile

	# Check whether the required packages are loaded...
	if (!defpac("tv")) { error (1,"Please load \"image.tv\" first!") }
	if (!defpac("immatch")) { error (1,"Please load \"images.immatch\" first!") }

	# Check existence of input, reference and output images...
	unlearn chkimg
	chkimg (img, "access", imtype="", verbose=yes)
	if (!chkimg.ok) { return }
	img = chkimg.imroot//chkimg.imtype//chkimg.imrpstr

	chkimg (refimg, "access", imtype="", verbose=yes)
	if (!chkimg.ok) { return }
	refimg = chkimg.imroot//chkimg.imtype//chkimg.imrpstr

	chkimg (outimg, "noaccess", imtype="", verbose=yes)
	if (!chkimg.ok) { return }
	outimg = chkimg.imroot//chkimg.imtype//chkimg.imrpstr

	# Determine the size of reference image ...
	imgets (refimg, "naxis1") ; ncol = int(imgets.value)
	imgets (refimg, "naxis1") ; ncol = int(imgets.value)
	imgets (refimg, "naxis2") ; nrow = int(imgets.value)

	# Opening message...
	gdate() ; sysinfo()
	if (verbose) { 
	   print ("IMREGISTER: input     = \""//img//"\"")
	   print ("            reference = \""//refimg//"\"")
	}
	print (" ", >> txtlog)
	print ("IMREGISTER:  "//sysinfo.version//"  "//sysinfo.user//"@"//sysinfo.host//"  "//gdate.adate, >> txtlog)
	print ("    reference = "//refimg, >> txtlog)

	# Check whether coordinate file is supplied ...
	chkimg (coofil, "access", imtype="file", verbose=no)
	oldcoo = yes

	if (!chkimg.ok) {
	   # No file was supplied; create one using 'tvmark'.
	   # Display the reference image ...
	   xdisplay (refimg, z1=0., z2=0., zmode="linear", verbose=no)

	   oldcoo = no
	   if (access(coofil)) { delete (coofil, yes, verify=no) }
	   if (coofil=="") { coofil = mktemp("coo") }

	   # Let user mark reference objects interactively and write 'coofil'...
	   print ("Wait for the cursor cross to appear in the image display area, then")
	   print ("mark at least 3 objects (stars) for coordinate referencing using")
	   print ("the 'c' key (press q to quit).")
	   sleep(1)
	   tmptvmark = mktemp("tvmark")
	   tvmark (frame=1, coor="", logfile=tmptvmark, autolog=yes,
		outimage="", mark="circle", radii=str(sbox), lengths="0",
		label=no, number=no, nxoffset=9, nyoffset=0, txsize=1,
		pointsize=3, font="raster", color=204, commands="",
		interactive=yes, > "dev$null", >& "STDERR" )

	   # Write the coordinate file
	   list = tmptvmark
	   while( fscan(list,x,y) != EOF ) {
	      if ( x < 1. )  x = 1.
	      if ( x > ncol) x = ncol
	      if ( y < 1. )  y = 1.
	      if ( y > nrow) y = nrow
	      print (x,y, >> coofil)
	   }
	   delete (tmptvmark, yes, verify=no)

	}

	if (refrecen) {
	   print ("Recentering objects in the reference coordinate list ...")
	   ctrcoo (coofil, coofil, refimg, cntrbox=cbox, redispl=no,
		   verbose=no)
	}
	iwc (coofil) ; N1 = iwc.nlines
	rdlist (coofil, 1)
	cparse (rdlist.cline, delim=" ")
	x1 = real(cparse.field1)
	y1 = real(cparse.field2)


	# If a new reference coordinate file was created, save it...
	if (!oldcoo && coofil != chkimg.imroot ) {
	   chkimg (refimg, "access", verbose=no)
	   if (substr(chkimg.imrpstr,3,3) == "]") {
	      copy (coofil, chkimg.imroot//"."//substr(chkimg.imrpstr,2,2)//".coo", verbose=no)
	      print ("    refcoords = "//chkimg.imroot//"."//substr(chkimg.imrpstr,2,2)//".coo", >> txtlog)
	   } else {
	      copy (coofil, chkimg.imroot//".coo", verbose=no)
	      print ("    refcoords = "//chkimg.imroot//".coo", >> txtlog)
	   }
	} else {
	   print ("    refcoords = "//coofil, >> txtlog)
	}
	print ("    input     = "//img, >> txtlog)
	   
	# Display the image to be registered, and overlay the positions
	# of the stars in the reference frame ...
	xdisplay (img, z1=0., z2=0., zmode="linear", verbose=no)
	tvmark (frame=1, coor=coofil, logfile="", autolog=yes,
                outimage="", mark="plus", radii="0", lengths="0",
                label=no, number=yes, nxoffset=9, nyoffset=0, txsize=1,
                pointsize=3, font="raster", color=207, commands="",
                interactive=no, > "dev$null", >& "STDERR" )

	# Ask the user to mark all the stars ...
	print ("Wait for the cursor cross to appear in the image display area, then")
	print ("mark all reference stars in the same order using the 'c' key.")
        sleep(1)
	tmptvmark = mktemp("tvmark")
	tvmark (frame=1, coor="", logfile=tmptvmark, autolog=yes,
                outimage="", mark="circle", radii=str(cbox), lengths="0",
                label=no, number=no, nxoffset=9, nyoffset=0, txsize=1,
                pointsize=3, font="raster", color=204, commands="",
                interactive=yes, > "dev$null", >& "STDERR" )

	# Check whether all objects in the list were marked or whether
	# the positions in the reference image may serve as an initial
	# guess after taking into account a possible shift.
	iwc (tmptvmark) ; N2 = iwc.nlines
	rdlist (tmptvmark, 1)
	cparse (rdlist.cline, delim=" ")
	x2 = real(cparse.field1)
	y2 = real(cparse.field2)
	if ( N2 > N1 ) { 
	   print ("ERROR: More objects marked than read from reference list!", >> txtlog)
	   error (1,"More objects marked than read from reference list!")
	}
	if ( N2 < N1 ) {
	   print ("Guessing the positions of the unmarked reference objects...")
	   for ( NN = (N2+1) ; NN <= N1; NN = NN + 1 ) {
	      rdlist (coofil, NN)
	      cparse (rdlist.cline, delim=" ")
	      print ((real(cparse.field1)+x2-x1),(real(cparse.field2)+y2-y1), >> tmptvmark)
	   }
	}
	print ("TVMARK:  N_obj = "//str(N1)//",  N_marked = "//str(N2)//",  N_guessed = "//str(N1-N2), >> txtlog)

	# Show what we're doing...
	print ("Recentering objects in the coordinate list ...")
	ctrcoo (tmptvmark, tmptvmark, img, cntrbox=cbox, redispl=no,
		verbose=no)

	# Preparing the input coordinate file for 'geomap' ...
	tmpgeo = mktemp("tmpgeo_in")
	for ( NN = 1 ; NN < N1 ; NN = NN + 1 ) {
	    rdlist (coofil, NN) ; cstr = rdlist.cline
	    rdlist (tmptvmark, NN)
	    print (cstr//" "//rdlist.cline, >> tmpgeo)
	}

	# Derive the transformation (involving shifts, scaling and/or
	# image rotation.
	print ("GEOMAP:  Computing the spatial transformation function...")
	print ("GEOMAP:  Computing the spatial transformation function...", >> txtlog)
	tmpsol = mktemp("tmpgeo_sol")
	geomap (tmpgeo, "geomap.sol", 1., real(ncol), 1., real(nrow),
		transforms="", results="", fitgeometry="rscale",
		function="polynomial", xxorder=2, xyorder=2, xxterms="none",
		yxorder=2, yyorder=2, yxterms="none", calctype="real",
		reject=4., verbose=yes, interactive=no, graphics="",
		cursor="", > tmpsol)
	tail (tmpsol, nlines=8) | match ("Coordinate", stop=yes, >> txtlog)
	if (verbose) { tail (tmpsol, nlines=8) | match ("Coordinate", stop=yes) }
	delete (tmpsol, yes, verify=no)

	if (execute) {
	   print ("GEOTRAN:  "//img//" --> "//outimg)
	   print ("GEOTRAN:  interpol  = "//intpl//"   boundary="//bndry, >> txtlog)
	   print ("    "//img//" --> "//outimg, >> txtlog)
	   geotran (img, outimg, "geomap.sol", tmpgeo, geometry="linear",
		xin=INDEF, yin=INDEF, xshift=INDEF, yshift=INDEF, xout=INDEF,
		yout=INDEF, xmag=INDEF, ymag=INDEF, xrotation=INDEF,
		yrotation=INDEF, xmin=INDEF, xmax=INDEF, ymin=INDEF, ymax=INDEF,
		xscale=1., yscale=1., ncols=INDEF, nlines=INDEF, xsample=1.,
		ysample=1., interpolant=intpl, boundary=bndry, constant=0.,
		fluxconserve=yes, nxblock=512, nyblock=512, verbose=verbose)

	}

	if (access(tmptvmark))  { delete (tmptvmark, yes, verify=no) }
	if (access(tmpgeo))      { delete (tmpgeo, yes, verify=no) }
	if (access("geomap.sol")) { delete ("geomap.sol", verify=no) }

	if (verbose) { print ("Finished.") }
	print ("IMREGISTER: Finished.", >> txtlog)

end