An algorithm and implementation for displaying 3D prism maps showing data depending on geographic regions by raising a prism above each region. The program preprocessed the 2D map in advance so that new data could be displayed, with hidden surfaces removed, quickly.

1.1. Papers

[wrf-cga-82] WR Franklin. Program translates statistics into 3-D color map of europe. IEEE Computer and Applications, 2(5):front cover and p. 4, July 1982. (invited).
[wrf-prism-hu] WR Franklin. Prism -- a prism plotting program. In A. H. Schmidt, editor, Mapping Software and Cartographic Data Bases, Harvard Library of Computer Mapping, 1979 collection, pages 75-79. 1979.
[fl-3gdds-78] WR Franklin and H. R. Lewis. 3-D graphic display of discrete spatial data by prism maps. In Proc. SIGGRAPH'78, volume 12, pages 70-75, August 1978.

2. Overlaying Two Maps

An algorithm for calculating the areas of overlaid polygons without calculating the overlay itself, is presented. OVERPROP is useful when the sole purpose of overlaying two maps is to find some mass property of the resulting polygons, or for an areal interpolation of data from one map to the other. Finding the areas of all the output polygons is both simpler and more robust than finding the polygons themselves. OVERPROP works from a reduced representation of each map as a set of `half-edges'; with no global topology. It uses the uniform grid to find edge intersections. The method is not statistical, but is exact within the arithmetic precision of the machine. It is well suited to a parallel machine, and could be extended to overlaying more than two maps simultaneously, and to determining other properties of the output polygons, such as perimeter or center of mass. OVERPROP has been implemented as a C program, and is very fast.

One subtask is to locate all the points of each map in a specific polygon of the other, which is a useful application in its own right.

The largest example intersected counties of the coterminous US against hydrography polygons. They have these stats:

Map #0: 55068 vertices, 46116 edges, 2985 polygons
Map #1: 76215 vertices, 69835 edges, 2075 polygons

The total CPU time, in seconds, on an IBM T30 Thinkpad with a 1600MHz Pentium is about:

Execution times:
   Read map:                      0.87
   Scale vertices:                0.02
   Extract edges:                 0.08
   Calculate areas:               0.08
   Make grid:                     0.00
   Add map to grid:               0.10
   Intersect edges:               0.13
   Locate map 0 points in map 1:  0.40
   Locate map 1 points in map 0:  0.48
   Accumulate areas:              0.28
   Print areas:                   0.20

TOTAL TIME=                       2.65

2.1. Papers

[fs-ocamo-90] WR Franklin and Venkatesh Sivaswami. OVERPROP -- calculating areas of map overlay polygons without calculating the overlay. In Second National Conference on Geographic Information Systems, pages 1646-1654, Ottawa, 5-8 March 1990.
[fsskn-ofmop-93] WR Franklin, V. Sivaswami, D. Sun, M. Kankanhalli, and C. Narayanaswami. Calculating the area of overlaid polygons without constructing the overlay, Cartography and Geographic Information Systems, pages 81-89, April 1994. (paper), (code),
An algorithm and implementation for calculating the areas of overlaid polygons without calculating the overlay itself, is presented. OVERPROP is useful when the sole purpose of overlaying two maps is to find some mass property of the resulting polygons, or for an areal interpolation of data from one map to the other. Finding the areas of all the output polygons is both simpler and more robust than finding the polygons themselves. OVERPROP works from a reduced representation of each map as a set of `half-edges'; with no global topology. It uses the uniform grid to find edge intersections. The method is not statistical, but is exact within the arithmetic precision of the machine. It is well suited to a parallel machine, and could be extended to overlaying more than two maps simultaneously, and to determining other properties of the output polygons, such as perimeter or center of mass. OVERPROP has been implemented as a C program, and is very fast. The execution time on a Sun Sparcstation 10/30 to combine US counties, with 55068 vertices and 2985 polygons with hydrography polygons, with 76215 vertices and 2075 polygons, was 19 CPU seconds, excluding I/O time. This included the subtask of locating all the points of each map in a specific polygon of the other, which is a useful application in its own right.
[f-moaap-92] WR Franklin. Map overlay area animation and parallel simulation. In D. H. Douglas, editor, Proceedings, SORSA'92 Symposium and Workshop, pages 200-203, July 28-August 2 1992.
[f-cmopa-90] WR Franklin. Calculating map overlay polygon' areas without explicitly calculating the polygons -- implementation. In 4th International Symposium on Spatial Data Handling, pages 151-160, Zürich, 23-27 July 1990.
[fsmoa-83] WR Franklin. A simplified map overlay algorithm. In Harvard Computer Graphics Conference, Cambridge, Mass, USA, 31 July - 4 August 1983.

Code.

3. Logic Programming for Map Overlay

3.1. Paper

[wf-lpacm-90] P. YF Wu and WR Franklin. A logic programming approach to cartographic map overlay. Canadian Computational Intelligence Journal, 6(2):61-70, 1990.

4. Compressing Terrain Data

Very efficient, wavelet-based, lossy and lossless compression of gridded terrain elevation data.

4.1. Papers

[wrflossy] WR Franklin and Amir Said. Lossy compression of elevation data. In Seventh International Symposium on Spatial Data Handling, Delft, August 1996.
Generic lossy image compression algorithms, such as Said & Pearlman's, perform excellently on gridded terrain elevation data. We used 24 random 1201x1201 USGS DEMs as test data. Lossless compression required 2.12 bits per point on average, with a range from 0.52 to 4.09. With lossy compression at a rate of 0.1 bits per point, the RMS elevation error averaged only 3 meters, and ranged from 0.25 to 15 meters. Even compressing down to 0.005 bits per point, or 902 bytes per cell, gave an indication of the surface's characteristics. This performance compares favorably to compressing with Triangulated Irregular Networks. If it is desired to partition the cell into smaller blocks before compressing, then even 32x32 blocks work well for lossless compression. However, for lossy compression, specifying the same bit rate for each block is inadequate, since different blocks may have different elevation ranges. Finally, preliminary tests suggest that the visibility indices of the points are robust with respect to even a quite aggressive lossy compression.
Keywords: lossy compression, Triangulated Irregular Network, terrain compression, elevation compression, visibility index.
[wrfelev95] WR Franklin. Compressing elevation data. In Fourth International Symposium on Large Spatial Databases - SSD '95. Portland, Maine, USA, 6-9 Aug 1995.
This paper compares several, text and image, lossless and lossy, compression techniques for regular gridded elevation data, such as DEMs. Sp_compress and progcode, the best lossless methods average 2.0 bits per point on USGS DEMs, about half the size of gzipped files, and 6.2 bits per point on ETOPO5 samples. Lossy compression produces even smaller files at moderate error rates. Finally, some techniques for compressing TINs will be introduced.

This paper is largely obsoleted by Lossy compression of elevation data. However, this paper also gives examples where partitioning the data before compressing generally does not increase the total compressed size, and suggests that that generally holds.

5. CONTOUR->GRID

Conversion of Elevation Contours to a Grid: Mike Gousie and I have new techniques for this classic problem, which reduce the terracing effect seen in many other methods, and produce a smoother, more realistic, result.

5.1. Papers

[gousie05] MB Gousie and WR Franklin. Augmenting Grid-based Contours to Improve Thin Plate DEM Generation Photogrammetric Engineering & Remote Sensing 71(1), pp. 69-79, 2005.
[gousie97] MB Gousie and WR Franklin. Converting elevation contours to a grid In Eighth International Symposium on Spatial Data Handling (SDH), Vancouver BC Canada, July 1998. Dept of Geography, Simon Fraser University, Burnaby, BC, Canada.
We present two new methods for approximating elevation data from contours to a grid. The first repeatedly interpolates new contour lines between the original ones. The second starts with any interpolated or approximated surface, determines its gradient lines, and does a Catmull-Rom spline interpolation along them to improve the elevations. We compare the new methods to a more classical thin-plate approximation on various data sets. The new methods appear visually smoother, with the undesirable terracing effect much reduced.

Talk slides

An old version of the Matlab code for the preliminary overdetermined Laplacian method is here
[acmgis2003] M Gousie & WR Franklin, Constructing a DEM from Grid-Based Data by Computing Intermediate Contours, ACM-GIS 2003, New Orleans, 7-8 Nov 2003

5.2. PhD

Michael B Gousie, Contours to Digital Elevation Models: Grid-based Surface Reconstruction Methods, 1998.

5.3. Masters

John Childs, Development of A Two-Level Iterative Computational Method for Solution of the Franklin Approximation Algorithm for the Interpolation of Large Contour Line Data Sets, 2003.

6. TIN-73

This is the Triangulated Irregular Network (TIN) program that I wrote at Simon Fraser University in the summer of 1973, at Tom Poiker's and Dave Douglas's suggestion. This was possibly the first TIN program written anywhere in a GIS context. The code is in PL/1 There are also some related programs, such as a simple convex hull code.

6.1. Code

6.2. Interview

Jeff Thurston, Looking Back and Ahead: The Triangulated Irregular Network (TIN) GEOinformatics, Vol 6, oct/nov 2003.

7. TIN

TIN is the current version of my efficient Triangulated Irregular Network computation program, for input points on a grid. Due to its compact and simple data structures, TIN can easily process all 1,442,401 points in a 1201x1201 level-1 DEM. Degeneracies, such as included points falling on the edge of the square, or on an existing edge, are handled. TINning the Lake Champlain West DEM until all points were included used about 70 CPU secs on a 1600-MHz Pentium linux box.

This is the 2001 version. It fixes an error in the 1999 version, where edges were not always swapped properly, and cleans things up. It replaces my 1973 TIN program.

7.1. Code (compressed tar file)

Alternative archive formats: uncompressed tar file, zip file

8. Visibility and Observer Siting

This work produced new algorithms and implementations relating to terrain visibility. These involve fast approximations to the Line of Sight (LOS) problem, where we wish to compute the visibility index of every point in a Digital Elevation Model (DEM) array describing some terrain. The object is to find a set of site which, together, can cover the whole area of interest. This work has produced some counterintuitive results. For example, for some terrain, there is no significant correlation between elevation and visibility index; on the average, higher is no better for observation.

8.1. Papers and contract reports

[siting-isprs-2004] WR Franklin & C Vogt Efficient Multiple Observer Siting on Large Terrain Cells, GIScience 2004 Third International Conference on Geographic Information Science University of Maryland Conference Center October 20-23, 2004. Extended abstract, the unpublished paper that it abstracts.
[siting-isprs-2004] WR Franklin & C Vogt Multiple Observer Siting on Terrain with Intervisibility or Lo-Res Data, XXth Congress, International Society for Photogrammetry and Remote Sensing, Istanbul, 12-23 July 2004. paper poster
[siting-teclostwg-2004] WR Franklin Siting Observers on Terrain, US Army Topographic Engineering Center, 28 Jan 2004.
[siting-2002] WR Franklin, Siting Observers on Terrain, Symposium on Spatial Data Handling, Ottawa, July 2002. (paper), (talk). Once upon a time, the conference had a website at http://www.geomatics2002.org/index_e.html .
This paper presents an experimental study of a new algorithm that synthesizes separate programs, for fast viewshed, and for fast approximate visibility index determination, into a working testbed for siting multiple observers jointly to cover terrain from a full level-1 DEM, and to do it so quickly that multiple experiments are easily possible. Both the observer and target may be at a given fixed distance above the terrain. The process operates as follows. (1) An approximate visibility index is calculated for each point in the cell under consideration. (2) A set of tentative observers is selected from the highly visible points. (3) The actual observers are selected from them, so as to cover as much of the cell of possibe, using a greedy algorithm. Various experiments with varying parameters were performed on the Lake Champlain West cell, with observations such as the following. (1) Forcing tentative observers to be well spaced was more important than using the most visible tentative observers. (2) Most of the new observers added (because they covered the most unseen points) were already visible to an existing observer. (3) Randomly deleting many tentative observers before final selection didn't reduce the final area covered.

KEYWORDS: terrain visibility, viewshed, line of sight, multiple observers

This is an expanded colorized version of conference paper.
[battelle] WR Franklin, CK Ray and S Mehta. Geometric algorithms for siting of air defense missile batteries, March 1994.
Geometric aspects of the visibility problem in the siting of air defense missile batteries were studied in this project. It theoretically analyzed the problem, produced several new and efficient algorithms, implemented them, and tested them on many cells of data. The theoretical analysis studied the terrain characteristics, formally defined the Hawk siting problem, formalized the optimal placement of observers, and considered the minimum elevation of an airplane flying over varying terrain. Three algorithms for finding the viewshed around a particular observer were studied. R3 is slow but accurate. R2 is much faster, yet almost as accurate. Xdraw computes an approximate viewshed with error bounds. Four visibility index algorithms were studied. VL runs a fixed number of lines of sight out from every possible observer. WeightF weights the points by distance. Weight approximates WeightF, and DEM/LOS skips points along each line of sight for increased speed. The visibility indices of 20 DTED level I cells of the Korean peninsula, totaling over 28 million points, were processed with VL. Experimentally, visibility index is not strongly correlated with elevation. A PC-based demonstration program with this data was exhibited to the user community. The accuracy was compared on data ranging from Korea to Dijon and Pripyat', showing that these methods do select the highest-visibility points.
Keywords: line of sight, visibility index, viewshed, siting, air defense, geometry, Digital Terrain Elevation Data.
[fr-hinbv-94] WR Franklin and C Ray. Higher isn't necessarily better: visibility algorithms and experiments. In Thomas C. Waugh and Richard G. Healey, editors, Advances in GIS Research: Sixth International Symposium on Spatial Data Handling, pages 751-770, Edinburgh, 5-9 Sept 1994. Taylor & Francis.
We describe several fast programs to compute viewsheds and weighted visibility indices for observation points in a raster terrain. These programs explore various tradeoffs between speed and accuracy. We have analyzed many cells of data; there is no strong correlation between a point's elevation and its weighted visibility index. However, the, very few, high visibility points tend to characterize features of the terrain. This work can form a basis for automatically locating observers jointly to cover a terrain region of interest.
[wrf-appvis-00] WR Franklin, Approximating visibility, GIScience 2000 plenary talk, 30 October 2000, Savannah, Georgia.
Determining the visibility indices of points in a terrain is an important application in cartography, but is very compute-intensive. This paper studies the relative importance of various factors in visibility computation, with a view to making the process more efficient. Preliminary results include the following. Scaling down elevations to 8-bits precision doesn't change the general visibility indices, but can cause artifacts, where the visibility index changes in steps. The height of the observers and targets above the terrain have little effect on the result. Changing the radius of interest completely changes the result. A larger ROI leads to finer details in the visibility index map. It also leads to a smaller average visibility index. In every test case but one the visibility indices appear Poisson distributed. We approximate the visibility index by running various numbers of lines of sight from each possible observer to many random targets. Ten targets per observer is insufficient. 30 targets per point appears to lead to a visibility map that resembles many more points. Nevertheless, there is still a visible improvement on increasing from 100 to 300.

Our goal is to use these properties in a visibility index program that will use industrial engineering sampling techniques to estimate the most visible points with the minimum number of observer-target tests.

9. GENERAL

Here are broader papers and talks covering more than one topic, or reviewing the field, or that don't fit into one of the above categories.

9.1. Papers and talks

[ica99] WR Franklin and MK. Gousie. Terrain elevation data structure operations. In 19th International Cartographic Conference & 11th General Assembly of the International Cartographic Association (ICA), Ottawa Canada, August 1999. The conference's former website was here: http://www.ccrs.nrcan.gc.ca/ica1999/.
We describe several past and proposed projects for converting elevation data to different formats, such as gridded, TIN, and spline, lossily compressing, applying to visibility and drainage to test the compression quality, and integrating everything into a test suite.

Our predominant elevation data structure is a grid (array) of elevations. We have also worked with a Triangulated Irregular Network (TIN), which the first author first implemented (under the direction of Poiker and Douglas) in 1973. This work has several themes, such as leveraging from recent results in applied math, skillfully implementing our algorithms, and testing them on realisticly large data sets.
Efficient calculation of visibility information was one project. This includes determination of the the viewshed polygon of a specific point, or the visibility index of all points. The data processed includes several 1201x1201 level-1 DEMs and DTEDs. One surprise was that often there was no correlation between elevation and visibility index.
Lossless and lossy compression of elevation data was another project. We used a wavelet-based technique on 24 random level-1 USGS DEMs. Lossless compression required 2.12 bits per point on average, with a range from 0.52 to 4.09. With lossy compression at a rate of 0.1 bits per point, the RMS elevation error averaged only 3 meters, and ranged from 0.25 to 15 meters. Even compressing down to 0.005 bits per point, or 902 bytes per cell, gave an indication of the surface's characteristics. This performance compares favorably to compressing with TINs. If it is desired to partition the cell into smaller blocks before compressing, then even 32 x 32 blocks work well for lossless compression.
Approximating elevation data from contours to a grid is the latest project. Our first new method repeatedly interpolates new contour lines between the original ones. The second starts with any interpolated or approximated surface, determines its gradient lines, and does a Catmull-Rom spline interpolation along them to improve the elevations. We compare the new methods to a more classical thin-plate approximation on various data sets. The new methods appear visually smoother, with the undesirable terracing effect much reduced. The third new method is a realization of the classic proposals to use a partial differential equation (PDE) to interpolate or approximate the surface. Previously, this was possible only on small data sets. However, with recent sparse matrix and multigrid techniques, it is now feasible to solve a Lagrangian (heat flow) PDE on a 1000x1000 grid. It is also possible to solve an overdetermined system of equations for a best fit. Here, even the known points are constrained to be the average of their neighbors; so there are more equations than unknowns. This causes the slope to be continuous across the contours, and smooths the result.
Our proposed future work is to combine these various projects into a test suite, and to test how well whether other terrain properties such as drainage patterns are preserved by these transformations. The goal is better to understand terrain with a view to formalizing a model of elevation.
Keywords: contours, interpolation, grids, elevation, visibility, compression, triangulated irregular network

Talk slides
[bu-apr2003] WR Franklin Computational and Geometric Cartography, Boston University, 30 Apr 2003.
[siena-apr2003] WR Franklin Computational and Geometric Cartography, Siena College, 22 Apr 2003.
Updated version of my GIScience 2002 talk.
[gpr2002] WR Franklin Observations in Support of Automation with GPR, (invited talk) The Use of Ground Penetrating Radar in Assessing the Condition of Transportation Infrastructure (Workshop) CenSSIS, RPI, 29-30 Oct, 2002.
[giscience2002] WR Franklin Computational and Geometric Cartography, (invited keynote talk) GIScience 2002, Boulder, Colorado, USA, 26-29 Sept, 2002.
[geospatial] WR Franklin. Geospatial terrain: algorithms and representations. unpublished paper, and related talk, April 2002.
Geospatial terrain elevation contains long-range, nonlinear relationships, such as drainage networks. Representations and operations, such as compression, that ignore this will be inefficient and possibly erroneous. For example, multiple correlated layers may become internally inconsistent. This paper presents various representations, and terrain operations (such as visibility), discusses criteria for a good representation, and attempts a start towards a formal mathematical foundation.
[opaal102] Wm Randolph Franklin Mathematical Issues in Terrain Representation, (talk), DARPA/NSF OPAAL Workshop and Final Program Review, 10 Jan 2002.
[ICIP] H Pedrini, WR Schwartz, and WR Franklin. Automatic Extraction of Topographic Features Using Adaptive Triangular Meshes, 2001 International Conference on Image Processing (ICIP-2001), Thessaloniki, Greece, 7-10 Oct 2001.
[ucgis-ancart] Harold Moellering, Keith Clarke, Robert Cromley, Wm. Randolph Franklin, Alan Saalfeld, Jon Kimerling, and Marc Armstrong, Analytical Cartography, UCGIS Emerging Research Themes in GIScience , white paper, (UCGIS = University Consortium for Geographic Information Science). Once upon a time it was online at www.ucgis.org/emerging/ac.pdf. Dec 2000.
[wrf-cagis-00] WR Franklin. Applications of analytical cartography, Cartography and Geographic Information Systems 27(3), 2000, pp. 225-237. This was also presented as a talk on 20 April 2000.
Several applications of analytical cartography are presented. They include terrain visibility (including visibility indices, viewsheds, and intervisibility), map overlay (including solving roundoff errors with C++ class libraries and computing polygon areas from incomplete information), mobility, and interpolation and approximation of curves and of terrain (including curves and surfaces in CAD/CAM, smoothing terrains with overdetermined systems of equations, and drainage patterns). General themes become apparent, such as simplicity, robustness, and the tradeoff between different data types. Finally several future applications are discussed, such as the lossy compression of correlated layers, and just good enough computation when high precision is not justified.
[dagstuhl] WR Franklin. Elevation data operations, Dagstuhl-Workshop on Computational Cartography, (by invitation), Nov. 1996.
We describe several past and proposed projects for converting elevation data to different formats, such as gridded, TIN, and spline, lossily compressing, applying to visibility and drainage to test the compression quality, and integrating everything into a test suite.
Keywords: TIN, Triangulated Irregular Network, gridded elevation data, lossy compression, visibility indices.
[wrf-ucgia] UCGIA Steering Committee. On the possible role(s) of a ``university consortium for geographic information and analysis'' (UCGIA). In Proceedings, American Congress on Surveying and Mapping / American Society for Photogrammetry and Remote Sensing '93, New Orleans, 1993.
[f-tcfg-92] WR Franklin. Tutorial on curve fitting for GIS. In D. H. Douglas, editor, Proceedings, SORSA'92 Symposium and Workshop, July 28-August 2 1992.
[f-cslld-91] WR Franklin. Computer systems and low level data structures for GIS. In D. Maguire, David Rhind, and M. Goodchild, editors, GIS: Principles and Practice, volume 1, pages 215-225. Longman Higher Education and Reference, London UK, 1991.