[GRASSLIST:5238] General purpose rubbersheeting, i.points, i.vpoints, and i.rectify.

[Tom Russo]

I've recently been trying to georeference scanned images of old (i.e. 100
year old) US Geological Survey topo maps.  I'm using GRASS 5.7.  I've run into 
some problems, naturally.

The maps all have Lat/Lon grids, and are all in "polyconic projection" --- 
with no specified central meridian, so it is pointless to try to georeference 
them into a GRASS polyconic projection location.  Trying to use a first order
affine transformation in i.rectify to georeference the map into a UTM location 
does not work, because we need more than a simple scaling, rotation, and 
translation.

The USGS itself doesn't even do that when it georeferences its DRGs:
they create a set of control points based on the lat/lon grid intersections,
and do a piecewise linear rubbersheeting into UTM.  

Since GRASS doesn't have a piecewise linear rubbersheeting module, I tried
to use i.group,i.target,i.points, and i.rectify with a transformation of 
higher order than the simple linear affine transformation, in the hopes that
it would be a reasonable approximation.  I ran into some issues:

1) i.points has an "analyze" function that is supposed to tell the RMS error
   of the transformation based on the points you supply.  However, it only
   uses the linear affine transformation, and cannot tell you the quality of
   any higher order transformation you might use in i.rectify

2) i.rectify's documentation implies that it reports RMS error, but in fact
   it does not --- in order to get any output at all that tells how well
   higher order transformations are doing, you must recompile i.rectify with
   a "#define BDEBUG" uncommented in imagery/i.rectify/crs.c to get any such
   output, and it reports error in the image coordinates, not in meters in
   projected coordinates.
 
3) there is a module called i.vpoints that apparently does have code that would
   allow you to specify 1st, 2nd, or 3rd order transformations in the analys
   function, but when I attempt to use that menu it enters some sort of 
   infinte loop waiting for mouse input, and no matter what I do it doesn't
   respond to mouse clicks or keyboard input.  It has to be killed from another
   window.

By using the debugging output of i.rectify with the debug #define included
I can have a tiny bit of confidence that choosing a 3rd order transformation
using the 16 lat/lon grid intersections is somewhat reasonable, but it's not
as good as having the simple "is the RMS error below the pixel size?" test
that the first order affine transformation can have.  

Does anybody use i.vpoints?  Has anyone else seen this strange behavior when
choosing the transformation order menu?

Is anyone aware of free (or cheap) software that can perform general
piecewise linear rubbersheeting of rasters to georeference them into a 
projection that doesn't match the ones they were created in?  USGS refers
to some in-house software that is unavailable.  There's apparently a
set of Matlab image processing functions that can do this (the method 
creates a Delaunay triangulation of the net of control points, and does
a local linear affine transformation of points within each triangle).  But
I can't find anything that could be used easily with GRASS.  Does anyone have
suggestions?

-- 
Tom Russo    KM5VY     SAR502  DM64ux         http://www.swcp.com/~russo/
Tijeras, NM  QRPL#1592 K2#398  SOC#236 AHTB#1 http://www.qsl.net/~km5vy/
 "When life gives you lemons, find someone with a paper cut."