include <foc.h>

# rsubset -- extract a subset of a reseau grid
# This task extracts a rectangular subset of a reseau entry.
# The region to extract is specified by starting and ending column
# and row numbers.  These are not pixel numbers but rather the index
# numbers of reseau marks.
# If the output 'outres' is "STDOUT" or null or blank, this task prints
# the X and Y positions within the subset on the standard output.  The
# X and Y positions for a reseau position are separated by a comma,
# and the lines are printed last row first (i.e. 'rlast' at the top).
#
# Phil Hodge, 19-Mar-1992  Task created

procedure rsubset()

char	inres_list[SZ_FNAME]	# list of input reseau table names
char	entry[SZ_FNAME]		# entry name
char	outres_list[SZ_FNAME]	# list of output reseau table names
int	cfirst, clast		# first and last columns to be extracted
int	rfirst, rlast		# first and last rows to be extracted
bool	verbose			# print table names?
#--
pointer sp
pointer irp, orp		# pointers to reseau table structs
pointer rlist			# list of reseau entries
pointer ient, oent		# reseau entry pointers
char	inres[SZ_FNAME]		# name of one input reseau table
char	outres[SZ_FNAME]	# name of one output reseau table
int	ilist			# for list of input tables
int	olist			# for list of output tables
int	ent_num			# entry number in input
int	inx, iny		# size of input reseau grid
int	onx, ony		# size of output reseau grid
int	junk
bool	just_print		# print positions rather than extract?
int	fntopnb(), fntlenb(), fntgfnb()
int	clgeti()
bool	clgetb()
bool	streq()
pointer rs_open(), rs_template()
int	rs_readr(), rs_list()
int	tbtacc()
bool	rs_check()
define	next_table_	91

begin
	call clgstr ("inres", inres_list, SZ_FNAME)
	call clgstr ("entry", entry, SZ_FNAME)
	call clgstr ("outres", outres_list, SZ_FNAME)

	ilist = fntopnb (inres_list, NO)		# don't sort

	if (outres_list[1] == EOS || outres_list[1] == ' ' ||
		streq (outres_list, "STDOUT")) {
	    just_print = true
	} else {
	    just_print = false
	    olist = fntopnb (outres_list, NO)		# don't sort
	    if (fntlenb (ilist) != fntlenb (olist))
		call error (1, "input & output lists not the same length")
	}

	cfirst = clgeti ("cfirst")
	clast  = clgeti ("clast")
	rfirst = clgeti ("rfirst")
	rlast  = clgeti ("rlast")
	verbose = clgetb ("verbose")

	if (cfirst < 1 || rfirst < 1 || clast < cfirst || rlast < rfirst)
	    call error (1, "invalid subset")

	# Loop over each table in list.
	while (fntgfnb (ilist, inres, SZ_FNAME) != EOF) {
	    if (!just_print)
		junk = fntgfnb (olist, outres, SZ_FNAME)

	    call smark (sp)

	    # Open input reseau table.
	    iferr {
		irp = rs_open (inres, READ_ONLY, NULL)
	    } then {
		call eprintf ("%s is not a reseau table\n")
		    call pargstr (inres)
		goto next_table_
	    }

	    if (verbose && !just_print) {
		call printf ("%s --> %s\n")
		    call pargstr (inres)
		    call pargstr (outres)
		call flush (STDOUT)
	    }

	    # Allocate space for an input entry.
	    call rs_alloc (irp, TY_REAL, ient)

	    # Get the sizes of the input & output grids.
	    inx = RES_NCOLS(irp)
	    iny = RES_NROWS(irp)
	    onx = clast - cfirst + 1
	    ony = rlast - rfirst + 1
	    if (clast > inx || rlast > iny) {
		call rs_close (irp)
		call eprintf ("input table `%s' is only %d by %d\n")
		    call pargstr (inres)
		    call pargi (inx)
		    call pargi (iny)
		goto next_table_
	    }

	    if (!just_print) {
		# Open output reseau table and specify size.
		if (tbtacc (outres) == YES) {		# does it exist?
		    orp = rs_open (outres, READ_WRITE, NULL)
		    if (!rs_check (orp, ony, onx)) {
			call rs_close (orp)
			call rs_close (irp)
			call eprintf (
			"existing output table `%s' is the wrong size\n")
			    call pargstr (outres)
			goto next_table_
		    }
		} else {
		    orp = rs_open (outres, NEW_FILE, NULL)
		    call rs_new (1, 1, ony, onx, TY_REAL, orp)
		}

		# Allocate space for an output entry.
		call rs_alloc (orp, TY_REAL, oent)
	    }

	    rlist = rs_template (irp, entry)		# process these entries

	    # Do for each reseau entry in current table.
	    ent_num = rs_list (rlist)
	    while (rs_readr (irp, ient, "", ent_num) != RES_F_NONEXTENTRY) {

		if (just_print) {
		    # Print table and entry names.
		    call printf ("\n# table `%s', entry `%s'\n")
			call pargstr (inres)
			call pargstr (RES_ENT_ENTRY(ient))
		    # Print reseau positions.
		    call rsub_print (Memr[RES_XPT(ient)], Memr[RES_YPT(ient)],
			    inx, iny, cfirst, clast, rfirst, rlast)

		} else {

		    # Fill in values for info columns.
		    call rsub_upd (orp, ient, oent)

		    # Copy a subset of the input grid into the output grid.
		    call rsub_copy (Memr[RES_XPT(ient)], Memr[RES_YPT(ient)],
			    Memb[RES_FPT(ient)],
			    Memr[RES_XPT(oent)], Memr[RES_YPT(oent)],
			    Memb[RES_FPT(oent)],
			    inx, iny, cfirst, rfirst, onx, ony)

		    # Save entry in output reseau table.
		    call rs_write (orp, oent, 0)
		}

		ent_num = rs_list (rlist)		# get next number
	    }

	    if (!just_print)
		call rs_close (orp)
	    call rs_close (irp)

next_table_
	    call sfree (sp)
	}

	if (!just_print)
	    call fntclsb (olist)
	call fntclsb (ilist)
end


# rsub_upd -- update info columns

procedure rsub_upd (orp, ient, oent)

pointer	orp		# i: res pointer for output table
pointer ient		# i: entry pointer for input table
pointer	oent		# i: entry pointer for output table
#--

begin
	call strcpy (RES_ENT_ENTRY(ient), RES_ENT_ENTRY(oent), SZ_R_ENTRY)

	# Copy tracking info and append an S for subset.
	call strcpy (RES_ENT_TRACKING(ient), RES_ENT_TRACKING(oent),
		SZ_R_TRACKING)
	call rs_t_update (oent, "S", SZ_R_TRACKING)

	RES_ENT_NROWS(oent) = RES_NROWS(orp)
	RES_ENT_NCOLS(oent) = RES_NCOLS(orp)

	RES_ENT_SFIELD1(oent) = RES_ENT_SFIELD1(ient)
	RES_ENT_SFIELD2(oent) = RES_ENT_SFIELD2(ient)
	RES_ENT_SFIELD3(oent) = RES_ENT_SFIELD3(ient)
	RES_ENT_SFIELD4(oent) = RES_ENT_SFIELD4(ient)
end


# rsub_print -- print X & Y positions
# This routine prints the X and Y positions within the subset.

procedure rsub_print (xin, yin,
		inx, iny, cfirst, clast, rfirst, rlast)

real	xin[inx,iny]		# i: x coordinates of input reseau marks
real	yin[inx,iny]		# i: y coordinates of input reseau marks
int	inx, iny		# i: size of input reseau grid
int	cfirst, clast		# i: first and last columns
int	rfirst, rlast		# i: first and last rows
#--
int	i, j			# indexes

begin
	# Print the values.
	do j = rfirst, rlast {
	    do i = cfirst, clast {
		call printf ("%10.3f %10.3f\n")
		    call pargr (xin[i,j])
		    call pargr (yin[i,j])
	    }
	}
	call flush (STDOUT)
end

# rsub_copy -- copy a subset to output
# This routine copies a subset of the input grid to the output grid.

procedure rsub_copy (xin, yin, fin, xout, yout, fout,
		inx, iny, cfirst, rfirst, onx, ony)

real	xin[inx,iny]		# i: x coordinates of input reseau marks
real	yin[inx,iny]		# i: y coordinates of input reseau marks
bool	fin[inx,iny]		# i: flag -- is reseau indef?
real	xout[onx,ony]		# o: x coordinates of output reseau marks
real	yout[onx,ony]		# o: y coordinates of output reseau marks
bool	fout[onx,ony]		# o: flag -- is reseau indef?
int	inx, iny		# i: size of input reseau grid
int	cfirst, rfirst		# i: column & row of beginning reseau number
int	onx, ony		# i: size of output reseau grid
#--
int	ii, ji			# indexes in input arrays
int	io, jo			# indexes in output arrays

begin
	ji = rfirst		# row number of first reseau to copy
	do jo = 1, ony {
	    ii = cfirst		# column number of first reseau to copy
	    do io = 1, onx {
		xout[io,jo] = xin[ii,ji]
		yout[io,jo] = yin[ii,ji]
		fout[io,jo] = fin[ii,ji]
		ii = ii + 1
	    }
	    ji = ji + 1
	}
end