Local Stretch Theory

Unlike aerial imagery made with a lens, old maps and other pictorial sources such as satellite scanner data do not obey optical laws which can be described by a global mathematical transformation  applied to the entire image. For well over a quarter of a century, satellite imagery made with line ("push-broom") scanners has been treated using bi-variate polynomials (Bernstein, 1976), usually of order 2. This method is frequently incorporated in the "rubber stretch" plugins available for GIS programmes. It is adequate for the application for which it was designed, but it may produce very large errors when the object to be transformed is an ancient map with irregular distortion to a modern large scale digital map. 

Old maps may contain valuable information about archaeological sites which have been destroyed by modern agriculture or construction. The problem of fitting old property boundary maps to modern ones is well known to cartographers who call the process "conflation", but many solutions are considered to be unsatisfactory.

Transformation of such data requires a localised approach which restricts operations to small regions around corresponding points. It is desirable that errors at these points be minimal or zero, and that the distribution of errors in the intervening areas be a linear function of the distance between the nearest points. It is also desirable that straight lines crossing the boundaries of local regions not be broken. It was shown that transforming small corresponding triangles drawn from the control points with a linear (affine) transformation would meet these conditions (Goshtasby, 1986) . However, practical implementation of this method for real-world imagery is rather more complex. The process of drawing non-intersecting lines between a set of points is called triangulation. An optimal triangulation was described by the Russian mathematician Delaunay and carries his name (Delaunay, 1934). The Delaunay triangulation has the desirable property that it maximises the smallest angle at any vertex over all triangles.  No vertex in the vertex set falls in the interior of the circumcircle (circle that passes through all three vertices) of any triangle in the triangulation:

For example, there are three ways of triangulating 5 points:  

but only one of them obeys the rule cited above:

A related problem is the construction of a polygon around each point which delimits its surrounding area relative to all near neighbour points. This problem was discussed earlier by another Russian mathematician, Voronoi (Voronoi, 1908).  The partitioning of a plane with n points into n convex polygons such that each polygon contains exactly one point and every point in a given polygon is closer to its central point than to any other is a Voronoi diagram. It is sometimes also known as a Dirichlet or Theissen tessellation. A Delaunay triangulation is said to be the dual of a Voronoi tessellation, because the lines of the triangulation intersect those of the tessellation at right angles and one diagram can easily be constructed from the other. 

This is a Voronoi tessellation of a collection of points:

and this is the dual Delaunay triangulation constructed by drawing lines at right angles to the Voronoi lines to connect the points:

Geometric construction of such figures for large numbers of points was greatly speeded up by the development of "sweep-line" or "divide and conquer" algorithms which reduce the number of operations to a linear multiple of the number of points times the logarithm of that number (Fortune, 1987). This improved on earlier less efficient methods (Sibson, 1981). 

Once the triangulation has been carried out, transformation of the corresponding triangles may proceed one pair at a time along with suitable interpolation from the modern map to the old map, pixel by pixel. Goshtasby recognised that a problem arises about what to do with the area outside the lines connecting the outermost points. This polygon is called the convex hull of the control points. He suggested a technique for dividing up the remaining area into triangles outside the hull which might be ambiguous for two images. A Voronoi tessellation of the whole space which is consistent right out to the corners might be used followed by affine transformation on the triangles contained within the Voronoi polygons, but there is an easier way.

Unlike the artificial examples given in Goshtasby's publication, it is likely that the corresponding control point set of an old map or non-optical image is randomly oriented with respect to those on the modern map. A preliminary projective transformation to bring these into approximate correspondence is needed or it won't be possible to be sure which pairs of triangles correspond. If that is done first, then one might as well use the corners of the whole image of the old map as artificial control points with a transformation to project them on to the new map which would then use these exterior points as the outer boundary. That puts all real control points inside the convex hulls of the four artificial points in image and map.  The resulting artificial triangles on the outer boundaries should not be not too narrow, hence a constant should be added or subtracted from the coordinates of the corner points all around the original image. 

Local error in the projective transformation due to hilly terrain may be partially compensated with the secondary bi-variate polynomial correction already available in the AirPhotoSE Multi-Point method due to Weidner. The degree of the polynomial may be chosen automatically as a function of the available number of control points. Since transformation with a polynomial greater than degree five fails due to rounding errors caused by hardware limitations, that is the most that is needed regardless of how many additional control points may be available.

In addition, there is a problem of stability of the independent triangulations in old and new images, with the possibility of degeneracy. Triangles become degenerate if three adjacent points lie on a straight line or when four points lie exactly on a circle, since there are two ways of drawing the pair of triangles depending on the orientation of the diagonal. These may have opposite orientations in old and new images. Degeneracy can be prevented by using a small random linear transformation on each control point before triangulation as suggested by (Beichl, 2002), but triangles which are inconsistent between old and new images may still result. Occasionally,  a triangle may be created in one image which does not exist in the other. A scan through the triangle lists prior to transformation must therefore reorder all entries for correspondence where needed, since the Delaunay algorithm may not create triangles in the same sequence for both sets of points. Furthermore, the vertex coordinates of corresponding triangles may come in any order. By comparing the ratios of the areas to the perimeters of each pair of triangles, remaining inconsistency may be detected. Offering a user-chosen threshold to select those which differ by more than a chosen percentage for separate treatment is then appropriate. A dual pass approach, flagging the consistent triangles in the first pass as 'don't touch', then using only the orientation of the new map triangles for the second pass is a simple solution.

In AirPhotoSE, control points are first checked for consistency. Then a Multi-Point Weidner transformation is used to produce approximate correspondence between source and target using constants determined by Singular Value Decomposition with all control points. If there are between 4 and 6 points, a pure projective transformation) is computed, If there are more, the additional bi-variate polynomial correction described above is introduced, the order of the correction polynomial depending on the available number of points. With 7 to 9 points, order 2 is used, 10 to 12, order 3, 13 to 17 order 4, 18 or more, order 5. This completes the first correction phase. 

The set of points in the preliminary image now become a new source and those of the modern map are used to compute Delaunay triangulations for each. Then for corresponding triangles, an affine transformation is computed and applied to transfer the data from the previously rectified source image to the original modern map. This produces the desired additional small local stretch  around each of the control points. Any errors in the other image picture element positions are confined locally to the immediate area around the point. 

Several different options for displaying the result are offered, including merging of the lines of the old with the new map where the new has priority over the old, overwriting the new map with the old one leaving other areas unchanged, creating an undistorted version of the old map or mixing the lines both in user-selectable proportions. A diagnostic option is also included to aid with placement of control points. It shows either all the triangles in both images or only those which are inconsistent. 

The method may also be used on oblique aerial images if there are a fairly large number (7 or more) of visible control points. This may be used to improve the accuracy of single image rectification if desired.