############################################################################
#
#	File:     tgrlink.icn
#
#	Subject:  Program to combine TIGER line chains
#
#	Authors:  Gregg M. Townsend and William S. Evans
#
#	Date:     June 23, 2000
#
############################################################################
#
#  Tgrlink connects records from a line chain file to produce a more
#  compact file composed of fewer, longer chains.  Chains having common
#  endpoints and somewhat similar orientations are joined together.
#  Then, wherever three consecutive points are collinear, or nearly so, 
#  the middle point is removed.
#
#  Usage: tgrlink [-e maxerror] [-a maxangle] [file.lch]
#
#  The maxerror parameter, measured in latitude units, sets the maximum
#  distance the middle of three points can deviate from the line connecting
#  its neighbors and still be considered "collinear".  The default value
#  is 4, which is generally large enough to cover quantization errors
#  without introducing visible artifacts.
#
#  The maxangle parameter, defaulting to 30 degrees, limits the change in
#  angle of the chain path due to the removal of a middle point.  This
#  prevents narrow rectangles from turning into pointed triangles.
#
#  The input file must be randomly seekable (a disk file, not a pipe).
#
############################################################################
#
#  The algorithm is effective but not perfect.  It is designed to
#  minimize memory to allow the handling of large input files.
#  Processing the output data a second time may give a little more
#  improvement.
#
#  First, the input file is scanned and each chain is entered in a table.
#  Chains are segregated by feature and boundary code (chains with
#  different codes cannot be combined) and grouped by orientation.
#
#  A table key is formed by concatenating latitude+longitude with
#  latitude (only), using whichever endpoint gives a smaller sum.  The
#  table value for a chain is the chain's offset in the input file.
#  If multiple chains share the same key, a list of offsets is entered
#  in the table.
#
#  Output is generated by iterating through all the codes from the
#  "least important" to "most important" (so that those end up on top
#  when the map is drawn).  Within codes, vertically oriented lines
#  come first, then horizontally oriented lines, followed by others.
#  Within an orientation group, chains are sorted by key, with the
#  effect that they are produced from upper left to lower right
#  along a diagonally oriented wavefront.
#
#  For each generated key, output proceeds as follows, given the file
#  offset o associated with the key.  If offset o has already been
#  processed, as noted in the set "done", then do nothing further.
#  Otherwise, add o to the set and continue.  Seek the input file to
#  offset o and read the chain data into memory.  Calculate the far
#  endpoint of the chain and the key associated with that.  Check the
#  tables for another unprocessed chain of similar orientation beginning
#  there; if successful, append the path and mark that chain as processed.
#  Repeat this as long as a successor chain can be found.
#
#  Now go through the chain in memory and collapse collinear points within
#  the limits permitted by the command options.  Finally, calculate the
#  maximum range of the chain from its starting point, and write it out.
#
#  This all seems to work well in practice.  One possible drawback, at
#  in theory, is that chains heading slightly more north than east will
#  not be connected to chains heading slightly more east than north.
#  The sorting of keys by latitude+longitude means that no matter which
#  chain is processed first, the wrong endpoint of the other one is in
#  the key table and no connection will be seen.
#
############################################################################
#
#  Links:  options
#
############################################################################


link options


$define DefaultError 4		# default max error for removing point
$define DefaultAngle 30		# default max angle for removing point

$define SECTORS 5		# number of different orientations


global ifile			# input file
global maxerr, maxangle		# point removal parameters

global latsin			# scaling factor: sin(latitude)

global chtab			# master chain table (keyed by code)
global done			# set of offsets already output

global xoff, yoff		# lists of deltas in current chain


record crec(code, key, x1, x2, y1, y2, rev, aindex)	# chain record data



#  main procedure

procedure main(args)
   local opts, w, hdr1, hdr2, e, k, l, latmin, latmax

   opts := options(args, "a.e.")	# process command options
   maxangle := \opts["a"] | DefaultAngle
   maxerr := \opts["e"] | DefaultError

   if *args > 1 then
      stop("usage: ", &progname, " file")
   else if *args = 1 then
      ifile := open(args[1]) | stop(&progname, ": can't open ", args[1])
   else
      ifile := &input

   hdr1 := read(ifile) | stop(&progname, ": empty file")
   hdr2 := read(ifile) | stop(&progname, ": file truncated")

   latmin := hdr1[16+:7]
   latmax := hdr2[16+:7]
   latsin := sin(((latmax + latmin) / 2.0) * (&pi / 9999999))

   loadfile()				# load table keys

   write(hdr1)
   write(hdr2)
   every dumpcode(kgen(chtab))		# dump chains in code order
end



#  loadfile() -- load input file keys into tables

procedure loadfile()
   local w, line, alist, t, l, r

   chtab := table()
   repeat {
      w := where(ifile) | stop(&progname, ": input file is not seekable")
      line := read(ifile) | break

      r := crack(line)

      if /(alist := chtab[r.code]) then {
         # first time for this code; make new tables.
         alist := chtab[r.code] := list(SECTORS)
         every !alist := table()
         }

      t := alist[r.aindex]

      ((/t[r.key]) := w) | {
         if type(l := t[r.key]) ~== "list" then
            l := t[r.key] := [t[r.key]]
         put(l, w)
         }

      }
   return
end



# kgen(t) -- generate keys of t in better order, as in the "tgrsort" script

procedure kgen(t)
   local l, k

   l := list()
   every k := key(t) do
      put(l, map(k[1], "FHEABCDX", "ZYXWVUTS") || k)
   l := sort(l)
   while k := pull(l) do
      suspend k[2:0]
   fail
end



#  dumpcode(code) -- output all chains having a particular code

procedure dumpcode(code)
   local h, v, i, l, k, o, alist

   alist := chtab[code]
   done := set()

   every l := sort(alist[aseq()], 3) do 
      while k := get(l) do {
         o := get(l)
         if type(o) == "list" then
            every putchain(code, k, !o)
         else
            putchain(code, k, o)
         }
   return
end



#  aseq() -- generate the orientation table subscripts in proper order

procedure aseq()
   local v, h

   h := 1 + integer(0.25 * SECTORS)
   v := 1 + integer(0.75 * SECTORS)
   suspend h				# sector that includes horizontal lines
   suspend v				# sector that includes vertical lines
   suspend h+1 to v-1			# NW to SE quadrant
   suspend 1 to h-1			# ENE to WSW
   suspend v+1 to SECTORS		# SSW to NNE
   fail
end



#  putchain(code, k, o) -- output chain of given code, key, and offset

procedure putchain(code, k, o)
   local t, r, x, y, x1, y1, xmin, xmax, ymin, ymax, d, w

   if member(done, o) then		# if already processed
      return
   insert(done, o)			# mark as done

   k ? {				# extract (x1, y1) from key
      t := move(8)
      x1 := integer(move(7))
      y1 := t - x1
      }

   xoff := []				# init list of deltas
   yoff := []
   r := putdel(o)			# add this chain's deltas

   while o := successor(r) do {		# while a successor can be found
      insert(done, o)			# mark it as processed
      r := putdel(o)			# append its deltas
      }

   collapse()				# collapse collinear points

   x := xmin := xmax := x1		# find min/max x/y values
   y := ymin := ymax := y1
   every x +:= !xoff do {
      xmin >:= x
      xmax <:= x
      }
   every y +:= !yoff do {
      ymin >:= y
      ymax <:= y
      }

   d := x - xmin			# find max deviation from x1 | y1
   d <:= xmax - x
   d <:= y - ymin
   d <:= ymax - y
   d >:= 9999				# limit to four digits

   # output the resulting chain

   writes(code, right(d, 4), right(x1, 7), right(y1, 7))
   while x := get(xoff) & y := get(yoff) do
      if x ~= 0 | y ~= 0 then
         w := writes(right(5000 + x, 4), right(5000 + y, 4))
   if /w then
      writes("50005000")		# line had degenerated to a point
   write()
   return
end



# putdel(o) -- record deltas (only) for chain at offset o in input file.

procedure putdel(o)
   local line, r, dy, mark

   # read the line located at offset o
   seek(ifile, o) | stop(&progname, ": can't reposition input file")
   line := read(ifile) |
      stop(&progname, ": input file changed during processing")
   # crack its data
   r := crack(line)

   # append deltas
   line ? {
      move(4)
      if ="|" then
         tab(upto('|') + 1)		# skip feature name
      move(18)

      if /r.rev then			# if endpoints were not reversed
         while put(xoff, move(4) - 5000) do
            put(yoff, move(4) - 5000)
      else {
         mark := &pos
         tab(0)				# if must start at far end
         while (mark < &pos) & (put(yoff, 5000 - move(-4))) do {
            put(xoff, 5000 - move(-4))
            }
         }
      }
   return r				# return cracked data
end



#  collapse() -- collapse collinear points in global xoff/yoff lists

procedure collapse()
   local maxsq, maxa, i, x1, y1, a1, x2, y2, a2, da, d, dx, dy

   if maxerr <= 0 then			# if no collapsing allowed
      return
   maxsq := maxerr * maxerr		# square of error (avoid sqrt later)

   maxa := maxangle * &pi / 180
   maxa >:= &pi				# max angle in radians

   x2 := latsin * xoff[1]
   y2 := yoff[1]
   a2 := atan(y2, x2)

   every i := 2 to *xoff do {
      x1 := x2
      y1 := y2
      a1 := a2
      x2 := latsin * xoff[i]
      y2 := yoff[i]
      a2 := atan(y2, x2)

      da := abs(a2 - a1)		# change in angle if removed
      if da > maxa then			# if too big, forget it
         next

      d := abs((x1 * x1 + y1 * y1) * sin(da))	# deviation from straight line
      if d <= maxsq then {			# if close enough
         dx := xoff[i] + xoff[i-1]
         dy := yoff[i] + yoff[i-1]
         if abs(dx) < 5000 & abs(dy) < 5000 then {	# if no overflow
            xoff[i] := dx				# set in curr deltas
            yoff[i] := dy
            xoff[i-1] := yoff[i-1] := 0			# zero previous deltas
            }
         }
      }
   return
end



#  successor(r) -- return offset of successor to chain given by crec record r

procedure successor(r)
   local k, alist, t, i, o, e

   alist := chtab[r.code]			# list, by orientation, for code
   k := right(r.x2 + r.y2, 8) || right(r.x2, 7)	# successor's key would be this
   every i := 0 | 1 | -1 do {			# try same orientation first
      t := alist[r.aindex + i] | next		# table of offsets
      if o := \t[k] then {			# entry can be int or list
         if type(o) ~== "list" then {
            if not member(done, o) then {
               return o
               }
            }
         else if (e := !o) & not member(done, e) then
            return e
         }
      }
   fail
end


# crack(line) -- return crec record giving data about chain

procedure crack(line)
   local angle, x1, y1, x2, y2, a
   static o
   initial o := crec()

   line ? {
      o.code := move(4)
      if o.code ||:= ="|" then		# if feature name present
         o.code ||:= tab(upto('|') + 1)

      move(4)				# skip old dimension measurement

      x1 := x2 := integer(move(7))
      y1 := y2 := integer(move(7))
      while x2 +:= move(4) - 5000 do
         y2 +:= move(4) - 5000

      if x1 + y1 > x2 + y2 then {	# if far endpoint has smaller sum
         o.rev := 1			# chain needs to be reversed
         x1 :=: x2
         y1 :=: y2
         }
      else
         o.rev := &null

      o.key := right(x1 + y1, 8) || right(x1, 7)
      o.x1 := x1
      o.y1 := y1
      o.x2 := x2
      o.y2 := y2
      a := atan(y2 - y1, latsin * (x2 - x1))
      o.aindex := 1 + integer(SECTORS * ((a / &pi) + 2.25)) % SECTORS
      }

   return o
end