include <imhdr.h>
include <error.h>		# for EA_ERROR and EA_WARN
include	<foc.h>
include "resfind.h"

# resfind -- reseau finding routine
# This task calls subroutines from rfindx and Dave Baxter's gridextrap
# to find reseau positions.
#
# Phil Hodge,  9-Jun-1993  Task created.
# Phil Hodge, 15-Feb-1994  Assume Xo,Yo from resinit is central reseau mark.
# Phil Hodge,  6-Feb-1998  Three calls to strcpy had only two arguments.

procedure resfind()

pointer input		# scratch for names of input images
pointer outres		# scratch for name of output reseau table
pointer entry		# scratch for entry name in reseau table
pointer rmodel		# scratch for name of reseau model text file
real	Xo, Yo		# estimated position of central reseau
real	Dxr		# X-offset to next row    \  Depends on position 
real	Dyr		# Y-offset to next row     | of (Irx,Iry) relative 
real	Dxc		# X-offset to next column  | to reseau (Px,Py)
real	Dyc		# Y-offset to next column /
int	searchx, searchy	# size of search box for finding reseau
real	depth		# correlation depth required for finding reseau
int	ncols, nrows	# size of reseau grid
int	verbosity	# how much info should we print?
bool	noindefs	# save approximate position if reseau not found?
#--
pointer sp
pointer optcrly		# scratch for value of optcrly from header
pointer geocorr		# scratch for value of geocorr from header
pointer image		# scratch for name of one image from list
pointer im_list		# image template list
pointer im		# for input image
pointer rm		# for reseau model structure
pointer rp, rentp	# pointers to reseau structures
char	an_entry[SZ_R_ENTNAME]	# entry name for current image
int	ix1, ix2	# indexes in IRp where the reseaux are located
int	iy1, iy2	# as ix1,ix2 but for Y coordinate
int	nfound		# number of reseaux actually found
int	samppln, linepfm		# size of current image
int	samppln_ref, linepfm_ref	# size of first (reference) image
real	IRp[NRX,NRY,3]	# reseau positions
real	xoff, yoff	# X & Y offsets of image from beginning of photocathode
bool	center_mark	# true if Xo,Yo are from resinit, presumed at center
bool	new_table	# true if output table doesn't exist yet
bool	first_image	# true for first input image
bool	verbose		# true if verbosity > 1; this is for rfindx
pointer immap()
pointer imtopen()
int	imtlen(), imtgetim()
real	clgetr()
int	clgeti()
bool	clgetb()
pointer rs_open()
int	imaccf()
int	tbtacc()

begin
	call smark (sp)
	call salloc (input, SZ_FNAME, TY_CHAR)
	call salloc (outres, SZ_FNAME, TY_CHAR)
	call salloc (entry, SZ_FNAME, TY_CHAR)
	call salloc (rmodel, SZ_FNAME, TY_CHAR)
	call salloc (optcrly, SZ_FNAME, TY_CHAR)
	call salloc (geocorr, SZ_FNAME, TY_CHAR)

	call salloc (image, SZ_FNAME, TY_CHAR)

	call clgstr ("input", Memc[input], SZ_FNAME)
	call clgstr ("outres", Memc[outres], SZ_FNAME)
	call clgstr ("entry", Memc[entry], SZ_FNAME)
	call clgstr ("rmodel", Memc[rmodel], SZ_FNAME)

	# ensure search area dimensions are odd integers
	searchx = 2 * (clgeti ("search1") / 2) + 1
	searchy = 2 * (clgeti ("search2") / 2) + 1
	depth = clgetr ("depth")

	verbosity = clgeti ("verbosity")
	verbose = (verbosity > 1)		# print tons of stuff?
	noindefs = clgetb ("noindefs")

	# Open image template list.
	im_list = imtopen (Memc[input])
	if (imtlen (im_list) < 1)
	    call error (1, "no input image specified")

	# Get the model data.
	call xrfgmodel (Memc[rmodel], rm)

	# Open output reseau table.
	if (tbtacc (Memc[outres]) == YES) {
	    # Write to an existing reseau table.
	    rp = rs_open (Memc[outres], READ_WRITE, NULL)
	    nrows = RES_NROWS(rp)
	    ncols = RES_NCOLS(rp)

	    # get space for output coords
	    call rs_alloc (rp, TY_REAL, rentp)

	    new_table = false
	} else {
	    # Create a new reseau table.
	    rp = rs_open (Memc[outres], NEW_FILE, NULL)
	    ncols = clgeti ("ncols")
	    nrows = clgeti ("nrows")

	    # This will be reset to false after we call rs_new.
	    new_table = true
	}


	# Do for each input image.
	first_image = true
	while (imtgetim (im_list, Memc[image], SZ_FNAME) != EOF) {

	    im = immap (Memc[image], READ_ONLY, NULL)
	    if (IM_NDIM(im) < 2) {
		call imunmap (im)
		call error (1, "image must be 2-D")
	    }

	    if (verbosity > 0) {
		if (!first_image)
		    call printf ("\n")
		call printf ("image `%s'\n")
		    call pargstr (Memc[image])
		call flush (STDOUT)
	    }

	    # Get initial estimates of reseau location and spacing.
	    iferr {
		call rf_get_xy (im, Xo, Yo, Dxr, Dyr, Dxc, Dyc,
			xoff, yoff, samppln, linepfm, center_mark)
	    } then {
		call rs_close (rp)		# close the reseau table
		call imunmap (im)
		call erract (EA_ERROR)
	    }

	    if (first_image) {
		samppln_ref = samppln
		linepfm_ref = linepfm
		first_image = false
	    } else if (samppln != samppln_ref || linepfm != linepfm_ref) {
		call eprintf (
	"Warning:  This image is not the same size as the first image:\n")
		call eprintf ("    `%s'\n")
		    call pargstr (Memc[image])
	    }

	    # If not given explicitly then use image ROOTNAME if present.
	    call xt_stripwhite (Memc[entry])
	    if (Memc[entry] == EOS) {
		if (imaccf (im, "ROOTNAME") == YES) {
		    call imgstr (im, "ROOTNAME", an_entry, SZ_R_ENTNAME)
		    call strlwr (an_entry)
		} else {
		    call strcpy (Memc[image], an_entry, SZ_R_ENTNAME)
		}
	    } else {
		call strcpy (Memc[entry], an_entry, SZ_R_ENTNAME)
	    }

	    # Set size of output reseau table, if it has just been created.
	    if (new_table) {
		# If the user didn't specify the number of columns or rows
		# find the expected number based on the image size.
		if (IS_INDEFI(ncols) || ncols < 1) {
		    ncols = IM_LEN(im,1) / RES_SPACING_X + 1
		    ncols = min (ncols, NRESX)
		}
		if (IS_INDEFI(nrows) || nrows < 1) {
		    nrows = IM_LEN(im,2) / RES_SPACING_Y + 1
		    nrows = min (nrows, NRESY)
		}
		call rs_new (1, 1, nrows, ncols, TY_REAL, rp)

		# get space for output coords
		call rs_alloc (rp, TY_REAL, rentp)

		new_table = false
	    }

	    # copy entry name to entry structure
	    call strcpy (an_entry, RES_ENT_ENTRY(rentp), SZ_R_ENTRY)

	    # Find the reseau positions, starting with (Xo, Yo).
	    iferr {
		call gridextrap (im, rm, Xo, Yo, Dxr, Dyr, Dxc, Dyc,
			xoff, yoff, searchx, searchy, depth, verbose, noindefs,
			IRp, nfound)
	    } then {
		call imunmap (im)
		call erract (EA_WARN)
		call eprintf ("skipping image `%s'\n")
		    call pargstr (Memc[image])
		next				# skip this image
	    }

	    # Find the subset of IRp that contains reseaux positions.
	    call rf_subset (im, IRp, ncols, nrows, center_mark, verbosity,
			ix1, ix2, iy1, iy2)

	    call imunmap (im)

	    if (verbosity > 0) {
		call printf ("%d reseaux found out of %d total\n")
		    call pargi (nfound)
		    call pargi ((ix2-ix1+1) * (iy2-iy1+1))
		call flush (STDOUT)
	    }

	    # Write the reseau positions.
	    call rf_save_res (IRp, ix1, ix2, iy1, iy2, noindefs,
		    Memr[RES_XPT(rentp)], Memr[RES_YPT(rentp)],
		    Memb[RES_FPT(rentp)])

	    # Insert the new entry; set date, tracking, and copy number of
	    # columns & rows.
	    call rf_set_info (rp, rentp)
	    call rs_insert (NULL, rp, rentp)
	}

	call rs_close (rp)		# close the reseau table

	call imtclose (im_list)		# close image template list

	call sfree (sp)
end

# rf_get_xy -- get initial reseau position and spacing
# Get initial estimates of the location of one reseau mark and the
# spacing between reseaux.  These depend on the relay and whether
# geometric correction has been done.  In addition, the location of
# the central reseau depends on image format (samppln, linepfm).

procedure rf_get_xy (im, Xo, Yo, Dxr, Dyr, Dxc, Dyc,
		xoff, yoff, samppln, linepfm, center_mark)

pointer im		# i: imhdr pointer
real	Xo, Yo		# o: estimated position of central reseau
real	Dxr		# o: X-offset to next row
real	Dyr		# o: Y-offset to next row
real	Dxc		# o: X-offset to next column
real	Dyc		# o: Y-offset to next column
real	xoff, yoff	# o: offset of image origin from photocathode origin
int	samppln		# o: width of image
int	linepfm		# o: height of image
bool	center_mark	# o: true if Xo,Yo are from resinit, presumed at center
#--
pointer sp
pointer optcrly		# scratch for value of optcrly from header
pointer geocorr		# scratch for value of geocorr from header
pointer pxformt		# scratch for value of pxformt from header
pointer pxlcorr		# scratch for value of pxlcorr from header
pointer x0par		# scratch for name of parameter for getting Xo
pointer y0par		# scratch for name of parameter for getting Yo
bool	got_XY		# did we get non-INDEF values from parameters x0,y0?
bool	subset		# true if image is a subset of original
real	clgetr()
int	imaccf(), imgeti()
bool	streq()

begin
	call smark (sp)
	call salloc (optcrly, SZ_FNAME, TY_CHAR)
	call salloc (geocorr, SZ_FNAME, TY_CHAR)
	call salloc (pxformt, SZ_FNAME, TY_CHAR)
	call salloc (pxlcorr, SZ_FNAME, TY_CHAR)
	call salloc (x0par, SZ_FNAME, TY_CHAR)
	call salloc (y0par, SZ_FNAME, TY_CHAR)

	# Get info from image header.
	if (imaccf (im, "geocorr") == YES) {
	    call imgstr (im, "geocorr", Memc[geocorr], SZ_FNAME)
	    call strlwr (Memc[geocorr])
	} else {
	    call strcpy ("perform", Memc[geocorr], SZ_FNAME)
	}
	if (imaccf (im, "optcrly") == YES) {
	    call imgstr (im, "optcrly", Memc[optcrly], SZ_FNAME)
	    call strlwr (Memc[optcrly])
	} else {
	    call strcpy ("f96", Memc[optcrly], SZ_FNAME)
	}
	if (imaccf (im, "pxformt") == YES) {
	    call imgstr (im, "pxformt", Memc[pxformt], SZ_FNAME)
	    call strlwr (Memc[pxformt])
	} else {
	    call strcpy ("normal", Memc[pxformt], SZ_FNAME)
	}

	# Get the offset between the current image and the (flipped)
	# photocathode origin.  Also get the flag that says whether
	# we have just a subset of the original image.
	call res_offset (im, xoff, yoff, subset)

	# Normally these will be INDEF, meaning that we should get the
	# values from the relay-dependent value in the resinit pset.
	Xo = clgetr ("x0")
	Yo = clgetr ("y0")
	got_XY = (!IS_INDEFR(Xo) && !IS_INDEFR(Yo))

	# If the user did not specify Xo,Yo with keywords "x0", "y0",
	# then we'll use the photocathode coordinates in the resinit pset.
	# These coordinates are for the central reseau mark (for the
	# format of our image).  But if we have a subset of the original
	# image, this position may not be the center of the subset.
	center_mark = !got_XY && !subset

	# Get the image size.
	if (imaccf (im, "samppln") == YES)
	    samppln = imgeti (im, "samppln")
	else
	    samppln = IM_LEN(im,1)
	if (imaccf (im, "linepfm") == YES)
	    linepfm = imgeti (im, "linepfm")
	else
	    linepfm = IM_LEN(im,2)

	if (streq (Memc[pxformt], "zoom")) {
	    if (imaccf (im, "pxlcorr") == YES) {
		call imgstr (im, "pxlcorr", Memc[pxlcorr], SZ_FNAME)
		call strlwr (Memc[pxlcorr])
	    } else {
		call error (1, "can't tell if image has been dezoomed")
	    }
	    if (streq (Memc[pxlcorr], "complete")) {
		samppln = 2 * samppln
	    } else {
		call error (1, "zoom-format image has not been dezoomed")
	    }
	}

	# Have we got Xo,Yo yet?
	if (!got_XY) {

	    if (streq (Memc[geocorr], "complete")) {
		call strcpy ("gx_", Memc[x0par], SZ_FNAME)
		call strcpy ("gy_", Memc[y0par], SZ_FNAME)
	    } else {			# raw, i.e. not geometrically corrected
		call strcpy ("rx_", Memc[x0par], SZ_FNAME)
		call strcpy ("ry_", Memc[y0par], SZ_FNAME)
	    }

	    if (streq (Memc[optcrly], "f48")) {
		call strcat ("48", Memc[x0par], SZ_FNAME)
		call strcat ("48", Memc[y0par], SZ_FNAME)
	    } else if (streq (Memc[optcrly], "f96")) {
		call strcat ("96", Memc[x0par], SZ_FNAME)
		call strcat ("96", Memc[y0par], SZ_FNAME)
	    } else {
		call error (1, "can't interpret OPTCRLY in image header")
	    }

	    if (samppln == 1024 && linepfm == 1024) {
		call strcat ("zlrg", Memc[x0par], SZ_FNAME)
		call strcat ("zlrg", Memc[y0par], SZ_FNAME)
	    } else if (samppln == 512 && linepfm == 1024) {
		call strcat ("nlrg", Memc[x0par], SZ_FNAME)
		call strcat ("nlrg", Memc[y0par], SZ_FNAME)
	    } else if (samppln == 512 && linepfm == 512) {
		call strcat ("n512", Memc[x0par], SZ_FNAME)
		call strcat ("n512", Memc[y0par], SZ_FNAME)
	    } else if (samppln == 256 && linepfm == 256) {
		call strcat ("n256", Memc[x0par], SZ_FNAME)
		call strcat ("n256", Memc[y0par], SZ_FNAME)
	    } else {
		call error (1,
		"no parameters for this format; specify x0,y0 explicitly")
	    }
	    Xo = clgetr (Memc[x0par])
	    Yo = clgetr (Memc[y0par])

	    # Subtract offsets so that Xo,Yo will be the position
	    # in the image, rather than in photocathode coordinates.
	    Xo = Xo - xoff
	    Yo = Yo - yoff
	}

	# Use different parameter names depending on whether the image has
	# been geometrically corrected and whether it's f/48 or f/96.
	if (streq (Memc[geocorr], "complete")) {
	    if (streq (Memc[optcrly], "f48")) {
		Dxr = clgetr ("gdxr_48")
		Dyr = clgetr ("gdyr_48")
		Dxc = clgetr ("gdxc_48")
		Dyc = clgetr ("gdyc_48")
	    } else if (streq (Memc[optcrly], "f96")) {
		Dxr = clgetr ("gdxr_96")
		Dyr = clgetr ("gdyr_96")
		Dxc = clgetr ("gdxc_96")
		Dyc = clgetr ("gdyc_96")
	    } else {
		call error (1, "can't interpret OPTCRLY in image header")
	    }
	} else {				# not geometrically corrected
	    if (streq (Memc[optcrly], "f48")) {
		Dxr = clgetr ("rdxr_48")
		Dyr = clgetr ("rdyr_48")
		Dxc = clgetr ("rdxc_48")
		Dyc = clgetr ("rdyc_48")
	    } else if (streq (Memc[optcrly], "f96")) {
		Dxr = clgetr ("rdxr_96")
		Dyr = clgetr ("rdyr_96")
		Dxc = clgetr ("rdxc_96")
		Dyc = clgetr ("rdyc_96")
	    } else {
		call error (1, "can't interpret OPTCRLY in image header")
	    }
	}

	call sfree (sp)
end

# res_offset -- determine image offset
# This routine gets the offset in X and Y of the image with respect to
# the full-frame photocathode in "flipped" coordinates.  The first step
# is to get header keywords sampoff & lineoff, which give the offset in
# X and Y of the FOC image with respect to the unflipped photocathode
# origin.  Then sampoff is modified to refer to the flipped photocathode.
# Finally, any offset due to an image section is included.
# If the image was taken in zoom mode, we assume that it was dezoomed.
# The subset variable will be true if 

procedure res_offset (im, xoff, yoff, subset)

pointer im		# i: imhdr struct for input image
real	xoff		# o: X offset of image from beginning of photocathode
real	yoff		# o: Y offset of image from beginning of photocathode
bool	subset		# o: true if an image section was taken
#--
pointer sp
pointer pxformt		# scratch for keyword value
pointer mw, ct		# mwcs pointers
real	lx, ly, px, py	# logical & physical coords of a pixel
real	sampoff		# offset not including image section
real	lineoff		# offset not including image section
int	samppln		# width of FOC image (not including image section)
int	linepfm		# height of image (not including section)
pointer mw_openim(), mw_sctran()
int	imaccf()
int	imgeti()
real	imgetr()
bool	streq()


begin
	call smark (sp)
	call salloc (pxformt, SZ_FNAME, TY_CHAR)

	if (imaccf (im, "samppln") == YES) {
	    samppln = imgeti (im, "samppln")
	    if (imaccf (im, "pxformt") == YES) {
		call imgstr (im, "pxformt", Memc[pxformt], SZ_FNAME)
		call strlwr (Memc[pxformt])
		if (streq (Memc[pxformt], "zoom"))
		    samppln = samppln * 2		# assumed it's dezoomed
	    }
	} else {
	    samppln = IM_LEN(im,1)
	}

	if (imaccf (im, "linepfm") == YES)
	    linepfm = imgeti (im, "linepfm")
	else
	    linepfm = IM_LEN(im,2)

	if (imaccf (im, "sampoff") == YES) {
	    sampoff = imgetr (im, "sampoff")
	    sampoff = MAX_NSAMPS - samppln - sampoff	# "flip" sampoff
	} else {
	    sampoff = 0.
	}

	if (imaccf (im, "lineoff") == YES)
	    lineoff = imgetr (im, "lineoff")
	else
	    lineoff = 0.

	# Get the offset due to the image section, if any.
	mw = mw_openim (im)
	ct = mw_sctran (mw, "logical", "physical", 3b)
	lx = 1.
	ly = 1.
	call mw_c2tranr (ct, lx, ly, px, py)
        call mw_ctfree (ct)
        call mw_close (mw)

	# Include the offset from the image section.
	xoff = sampoff + px - lx
	yoff = lineoff + py - ly

	# Do we have an image section, rather than the full image?
	subset = (IM_LEN(im,1) != samppln || IM_LEN(im,2) != linepfm)

	call sfree (sp)
end

# rf_subset -- locate subset
# Locate the relevant subset of the reseau positions, and return
# the first and last indices in X and Y for the subset.

procedure rf_subset (im, IRp, ncols, nrows, center_mark, verbosity,
		ix1, ix2, iy1, iy2)

pointer im		# i: imhdr struct for input image
real	IRp[NRX,NRY,3]	# i: reseau positions
int	ncols, nrows	# i: maximum number of reseau marks in X and Y
bool	center_mark	# i: true if Xo,Yo are from resinit, presumed at center
int	verbosity	# i: print pattern of reseaux found if verbosity > 0
int	ix1, ix2	# o: indexes in IRp containing relevant subset
int	iy1, iy2	# o: as ix1,ix2 but for Y coordinate
#--
int	sumx[NRX]	# sum of reseau found along each column
int	sumy[NRY]	# sum of reseau found along each row
int	i, j
int	i_l, i_r, j_b, j_t	# min & max indices, used for printing


begin
	# Initial values.
	i_l = NRX
	i_r = 1
	j_b = NRY
	j_t = 1

	do i = 1, NRX
	    sumx[i] = 0
	do i = 1, NRY
	    sumy[i] = 0

	do j = 1, NRY {
	    do i = 1, NRX {
		if (IRp(i,j,3) == YES) {	# position was found
		    sumx[i] = sumx[i] + 1
		    sumy[j] = sumy[j] + 1
		    i_l = min (i_l, i)
		    i_r = max (i_r, i)
		    j_b = min (j_b, j)
		    j_t = max (j_t, j)
		}
	    }
	}

	if (center_mark) {
	    # If we started with Xo,Yo as the central reseau mark for
	    # our image format, we can use that information to extract
	    # the appropriate subset.
	    call rf_center (ncols, nrows, ix1, ix2, iy1, iy2)
	} else {
	    # Find the appropriate range based on how many reseau were found.
	    call rf_range (im, sumx, NRX, ncols, ix1, ix2)
	    call rf_range (im, sumy, NRY, nrows, iy1, iy2)
	}

	# Print a pattern showing the locations of reseaux that were
	# found; use "x" if the reseau is within the range that we're
	# keeping and "o" if it's outside.
	if (verbosity > 0) {
	    do j = j_t, j_b, -1 {
		call printf ("  ")
		do i = i_l, i_r {
		    if (IRp(i,j,3) == YES) {
			if (i >= ix1 && i <= ix2 && j >= iy1 && j <= iy2)
			    call printf (" x")
			else
			    call printf (" o")
		    } else {
			call printf ("  ")
		    }
		}
		call printf ("\n")
	    }
	    call printf ("\n")
	    call flush (STDOUT)
	}
end

# rf_center -- take central subset in X and Y

procedure rf_center (ncols, nrows, ix1, ix2, iy1, iy2)

int	ncols, nrows	# i: number of reseau marks in subset X and Y
int	ix1, ix2	# o: indexes in IRp containing relevant subset
int	iy1, iy2	# o: as ix1,ix2 but for Y coordinate
#--

begin
	ix1 = NRX / 2 + 1 - ncols / 2
	ix2 = ix1 + ncols - 1

	iy1 = NRY / 2 + 1 - nrows / 2
	iy2 = iy1 + nrows - 1

	if (ix1 < 1 || ix2 > NRX || iy1 < 1 || iy2 > NRY)
	    call error (1,
		"reseau position gotten from resinit is not in the center")
end

# rf_range -- locate subset in X or Y
# Find the first and last indices where reseaux were found.

procedure rf_range (im, sum, nvals, maxres, i1, i2)

pointer im		# i: imhdr struct for input image
int	sum[ARB]	# i: sum of reseau found along each column (or row)
int	nvals		# i: size of array (number of columns or rows)
int	maxres		# i: max number of reseau for this size image
int	i1, i2		# o: indexes in IRp containing relevant subset
#--
int	sum_s		# sum of values in SUM array
int	half_sum	# half of the number of reseau found
int	midpt		# index of middle of SUM array
int	i

begin
	sum_s = 0
	do i = 1, nvals
	    sum_s = sum_s + sum[i]
	half_sum = (sum_s + 1) / 2

	# Find index that splits the number of reseau in half.
	sum_s = 0
	midpt = 1
	do i = 1, nvals {
	    sum_s = sum_s + sum[i]
	    if (sum_s >= half_sum)
		break
	    midpt = midpt + 1
	}

	i1 = midpt - maxres / 2
	i1 = max (i1, 1)
	i2 = i1 + maxres - 1
	if (i2 > nvals) {
	    i2 = nvals
	    i1 = i2 - maxres + 1
	    if (i1 < 1)
		call error (1, "grid size is too large")
	}

	# Was the range truncated by the limits (1,nvals)?
	if (i2 - i1 + 1 < maxres)
	    return

	# Consider shifting left or right by one.

	if (i1 <= 1) {

	    # No room to shift left, but should we shift right?
	    if (i2 < nvals) {
		if (sum[i1] < sum[i2+1]) {
		    i1 = i1 + 1			# shift right
		    i2 = i2 + 1
		}
	    }

	} else if (i2 >= nvals) {

	    # No room to shift right, but should we shift left?
	    if (i1 > 1) {
		if (sum[i1-1] > sum[i2]) {
		    i1 = i1 - 1			# shift left
		    i2 = i2 - 1
		}
	    }

	# Else consider shifting either left or right.

	} else if (sum[i1] < sum[i2+1]) {
	    i1 = i1 + 1				# shift right
	    i2 = i2 + 1
	} else if (sum[i1-1] > sum[i2]) {
	    i1 = i1 - 1				# shift left
	    i2 = i2 - 1
	}
end

# rf_save_res -- write coordinates to reseau table

procedure rf_save_res (IRp, ix1, ix2, iy1, iy2, noindefs,
		x, y, f)

real	IRp[NRX,NRY,3]	# i: reseau positions
int	ix1, ix2	# i: indexes in IRp containing relevant subset
int	iy1, iy2	# i: as ix1,ix2 but for Y coordinate
bool	noindefs	# i: save approximate position if reseau not found?
real	x[ARB]		# o: X coordinates of reseaux
real	y[ARB]		# o: Y coordinates of reseaux
bool	f[ARB]		# o: true if the reseau was not found
#--
int	n		# index in x, y, fout arrays
int	i, j

begin
	n = 1
	if (noindefs) {
	    # Copy position to output, even if it's only approximate.
	    do j = iy1, iy2 {
		do i = ix1, ix2 {
		    x[n] = IRp[i,j,1]
		    y[n] = IRp[i,j,2]
		    # true if position was not found
		    f[n] = (IRp(i,j,3) == NO)
		    n = n + 1
		}
	    }
	} else {
	    # Only copy position to output if it was actually found.
	    do j = iy1, iy2 {
		do i = ix1, ix2 {
		    if (IRp(i,j,3) == YES) {	# position was found
			x[n] = IRp[i,j,1]
			y[n] = IRp[i,j,2]
			f[n] = false
		    } else {
			x[n] = INDEFR
			y[n] = INDEFR
			f[n] = true		# yes, it was not found
		    }
		    n = n + 1
		}
	    }
	}
end

# rf_set_info -- assign stuff for entry in reseau table

procedure rf_set_info (rp, rentp)

pointer	rp		# i: res pointer
pointer	rentp		# i: entry pointer
#--
long	clktime()

begin
	call strcpy ("F", RES_ENT_TRACKING(rentp), SZ_R_TRACKING)
	RES_ENT_DATE( rentp) = clktime(0)

	RES_ENT_NROWS(rentp) = RES_NROWS(rp)
	RES_ENT_NCOLS(rentp) = RES_NCOLS(rp)

	RES_ENT_SFIELD1(rentp) = 1
	RES_ENT_SFIELD2(rentp) = 1
	RES_ENT_SFIELD3(rentp) = 1024
	RES_ENT_SFIELD4(rentp) = 1024
end