David N. Fogel

[IMAGE] Rocky Flats Sub-Area Map [IMAGE] Rocky Flats Sub-Area Image

Image Rectification with Radial Basis Functions: Application to RS/GIS Data Integration

Back to top


Remotely sensed digital imagery is inherently in need of geometric correction to maximize its usefulness within a GIS for subsequent analyses. This has long been problematic for large scale imagery such as acquired from airborne scanners. Radial basis functions, such as Hardy's multiquadrics and reciprocal multiquadrics, thin plate splines and variations of these methods offer ready alternatives to conventional polynomials and finite element methods for correction of these data. The nonstochastic methods of multiquadrics, thin plate splines and polynomials are reviewed in the context of scattered data interpolation and approximation. Applications to actual remotely sensed data acquired from airborne and satellite remote sensing systems (passive and active sensors) indicate that multiquadrics, reciprocal multiquadrics and thin plate splines produce excellent results while maintaining control point correspondence. Implementation specifics of these analytical functions, e.g. the development of the multiquadric algorithm for image warping, are reviewed based on the development of the software program MQReg. MQReg was conceived and implemented as a research tool during NCGIA Initiative 12: Integration of Remote Sensing and GIS. The status of MQReg as a research tool, its availability and extensions to current research are also discussed.

Back to top


Remotely sensed data preprocessing generally refers to geometric and radiometric corrections by data providers, to varying degrees, as well as subsequent image registration, rectification or orthorectification and further radiometric correction. Preprocessing by data providers (e.g. SPOT, EOSAT, USGS) infers a sense of data fidelity to the user who may not in fact be aware of the nature and the extent of the operations on the "raw" data. For example, the Landsat-5 imagery in Figures 1 and 2 illustrate different degrees of data preprocessing (Wivell 1995). However, far less stringent requirements apply to data fidelity with respect to subsequent geometric "preprocessing" by the users in applications such as change detection and classification. Specifically, this paper addresses this subsequent geometric preprocessing by reviewing conventional practice, offering new directions and providing an illustrative examples of the potential in this areas as an impetus for change. In comparison, radiometric preprocessing is a rather more dynamic research area, however it is beyond the scope of this paper.

Integrating remotely sensed data into a geographic information system provides valuable input for environmental modeling much the same way ancillary data has been used within digital image processing systems for classification and change detection. Common to both processing flows is remotely sensed data preprocessing to a format useful for subsequent analyses. This preprocessing requirement is increasingly important given developments in sensor technology and the increasing data stream of high resolution remotely sensed data. These new data are being acquired from both satellite and airborne sensors.

Issues surrounding the integration of remotely sensed data into a geographic information system (GIS) warranted the establishment of Initiative-12 by the National Center for Geographic Information and Analysis (NCGIA) in response to interests in both the public and private sectors. Information about the Initiative research agenda and substantive issues may be found in the proceedings from a special session at the 1991 ACSM-ASPRS Annual Convention (Star 1991) and NCGIA Technical Paper 93-4 (Estes and Star 1993).

[IMAGE] Landsat-5 Image with Level 0R Correction

Figure 1: Landsat-5 Image with Level 0R Correction.

Remote Sensing/Geographic Information System (RS/GIS) data integration is not a new topic. A raster-based information system for georeferenced data was coupled to an image processing system developed at the Jet Propulsion Laboratory (JPL) at the California Institute of Technology. This system is often referred to colloquially as VICAR/IBIS (Video Information Communication and Retrieval/Image Based Information System). Conceived as a subset of VICAR, IBIS included the ability to manipulate vector data. Marble and Peuquet (1983) provide an overview of GIS, RS/GIS, and an application using VICAR/IBIS.

[IMAGE] Landsat-5 Image with Level 1G Correction

Figure 2: Landsat-5 Image with Level 1G Correction.

This paper will review digital image warping, i.e. image registration and image rectification, in the context of change detection. Trotter (1991) provides an overview of the use of remotely-sensed data as an input to a GIS. Brown (1992) provides a good useful overview on image registration in general. Warping, as used here, refers to geometric transformations of digital imagery. An assessment of image registration issues and requirements follows in the form of a critique of conventional practice and a review of scattered data interpolation methods within a RS/GIS framework. Radial basis function methods are introduced and applied to airborne scanner data at the Department of Energy's Rocky Flats facility in Colorado.

Back to top


Three broad categories of geocorrection techniques may be identified for post-preprocessing of remotely sensed data based on inputs to the procedure: (1) correction models based on ephemeris platform and sensor information (scene/sensor models); (2) correction models relating corresponding ground control points or features; and (3) hybrid models. In (2) above, ground control may include vertical (height) information as well as relative location and features may be extended to include collections of points. The focus here is on the second category, or more specifically, on mathematically modeling the geometric distortions based on control point correspondence. This research is biased toward airborne scanner data which is intrinsically more difficult than satellite data to geometrically correct. It is assumed that, in some sense, the imagery is reasonable and not given to excessive sensor or platform variations.

Bivariate mapping polynomials are the most commonly used method for the registration and rectification of digital imagery. An early application of bivariate mapping polynomials to image warping is documented in Markarian, Bernstein et al. (1973) for the correction of "high-resolution" images from the ERTS series RBV (Landsat Remote-Beam Vidicon sensor system). Computing resources were comparatively scarce and expensive and these polynomials of the trend surface form are simple and can be reasonably effective. Regardless, these bivariate polynomials were still too computationally intensive for all subsequent preprocessing tasks.

One solution to this problem involved generating a grid from the ground control points based on an underlying geometric model (e.g. bivariate polynomials, piecewise linear methods) and approximating the warp. Thormodsgard (1987) provides information on this gridding approach which offers a trade-off between geometric fidelity and computational intensiveness. With satellite data, these methods are widely applied and most often used when registration or rectification is required. Grid density may be adjusted depending on the nature of the geometric distortions such that the amount of information lost during the gridding approximation is minimized.

The use of bivariate polynomial mapping functions in remote sensing is the most widely known approach to image warping despite its well-known shortcomings. Christensen, Jensen et al. (1988) describe common frustrations using bivariate mapping polynomials to warp airborne scanner data (5.6 meters-squared). In this case, geometric correction is both labor and computationally intensive. The lack of attention to geometric correction is evident in publications where additional geometric preprocessing has taken place with no mention of the method (e.g. Curran and Pedley 1990). It must be conjectured that bivariate polynomial mapping functions were applied to warp the data since the approach is most prevalent and rarely contrasted against alternatives.

Seriousness consequences of misregistration are underscored in the change detection literature. In a review of change detection methodology, Singh (p. 990, 1989) states, "It may be mentioned here that accurate spatial registration of the two images is essential for most change detection methods." On a global scale, Townshend, Justice et al. (p. 1054, 1992) conclude that more comprehensive satellite preprocessing is required for change detection stating, "The registration of data sets to a common spatial framework is an essential precursor to the use of remotely sensed data for monitoring change." Townshend et al. conduct their research using third-order bivariate polynomials seemingly because ".… geometric correction is [polynomials are] …. an established technique." (Heard, Mather and Higgins 1992).

A well-reasoned statement by Billingsley (p. 429, 1982) provides some perspective, "Given a certain loss in accuracy due to misregistration, how does that damage the ability to use the data analysis results? These evaluations will be discipline dependent, and must be sought separately." Billingsley's observation does not permit geometric distortions to be assumed away (out of convenience). It is of little consequence that there exists the common practice of reporting error figures loosely summarizing the warping accuracy. These root mean squared error (RMSE) figures usually underestimate misregistration. Interestingly, proprietary systems may not even calculate the same values for a metric as benign as a RMSE (cf. ERDAS Imagine 8.1 and PCI EASI/PACE 5.3).

Back to top

RS/GIS Data Integration

A basic knowledge of digital image processing and remote sensing fundamentals is requisite to understanding RS/GIS data integration issues. Reviews of remote sensing concepts may be found in Duggin and Robinove (1990) and Forshaw, Haskell et al. (1983). A digital image is basically a two-dimensional array of pixels having both radiometric characteristics and spatial components. The radiometric resolution is a function of spectral wavelength and bandwidth (among other things) and operationalized as a digital number (DN). Spatial resolution refers to a measurement or observation's local nature most often represented as a raster data type with relative location and size attributes. These attributes are not always sufficient to define spatial resolution since spatial resolution is not necessarily equivalent to pixel size (Forshaw et al. 1983). Remotely sensed data are often resampled to an arbitrary or convenient size bearing little relationship to the original spatial resolution.

Geometric distortions are generally classified as either internal sensor-related effects and external effects (Bernstein 1983). We will assume that most of the internal effects have been corrected in a prior preprocessing operation by the data provider (e.g. EOSAT). The principal distortions remaining are limited to internal effects (low-frequency) as well as high-frequency attitude and altitude effects (sensor-platform variations), and terrain effects (horizontal displacements). Airborne scanner distortions may be complicated by pronounced panoramic effects as well as S-Bend scanning effects.

A clear distinction between image registration and image rectification may be made based on whether the output domain is an image, a georeferenced or geocoded image, or a map (projection). In the most basic sense, registration refers to an image-to-image geometric transformation. Rectification may be either an image-to-geocoded image transformation or an image-to-map transformation. Significantly, there is not explicit modeling and removal of terrain effects. Orthorectification involves image rectification with explicit terrain effects removal. The distinction between registration and rectification is important with respect to discussing the ground control points (GCP) accuracy.

Image registration is perhaps the most straightforward procedure. It is likely that sufficient corresponding points can be identified given images with similar spatial and spectral characteristics. Extracting GCPs becomes more difficult when registering imagery from different data sources and/or at different resolutions. Even with proper attention to detail, it is reasonable to assume that there will be some error in both sets of GCPs, input and output (see, for example, the Total Least Squares Problem in Van Huffel and Zha 1993).

Image rectification is often more complicated. Recognizable features on a map are not always found on an image or may be difficult to recognize and vice-versa. Additionally, existing map accuracy standards may inhibit selection of representative GCPs. This situation arises when points in map space are less certain than the corresponding image points. This will become more problematic with the coming generation of high spatial resolution satellites.

In change detection applications, image-to-image registration or image-to-image rectification may have less geometric error than rectifying multiple images to a common map projection (see Niblack, 1986). On a case-specific basis, it may be determined that only the output product(s) be integrated into a GIS. If the output image is not orthorectified, then geometric distortions between images may be compounded during the correction process.

Digital orthophotos are not always well-suited for GCP extraction either. Clearly identifiable corresponding points are often lacking (e.g. cultural features) especially in the case of natural resource management. Rectification suffers from the same problems inherent in registration where errors exist in both the input and output GCPs. It is usually the case that these positional errors, horizontal and vertical, are ignored. This is true of digital elevation modeling (DEM) modeling as well where errors exist not only in the elevation data but in the locations of the elevation data too.

Back to top

Scattered Data Modeling

Geometric correction of imagery may be approached using bivariate functions where no assumptions are made with respect to the distribution of the locations of the control points in either map or image space. In fact, bivariate mapping polynomials are but one such method under this definition. Most literature on scattered data is in reference to interpolation. An interpolating function maintains the correspondence of the control points between images, whereas an approximating function will not honor the control point locations. Accordingly, it may be argued that the warping problem may be better approached as an approximation problem given the previously discussed difficulties in accurately locating corresponding GCPs. Bivariate mapping polynomials are most often approximating functions, though this is not justification for their (seemingly) universal application.

The expression "scattered data interpolation" was is used to distinguish applicable methods from otherwise structured-dependent techniques. Often, the objective is to obtain a contoured representation of the data. Survey literature may be found in reports, conference proceedings and monographs, as well as journal articles. These include, but are not limited to Franke (1982, 1987) Franke and Nielson (1991), Lam (1983), Lancaster and Salkauskas (1986), Watson (1992) and for 3D problems, Alfeld (1989). The RS/GIS community will recognize the potential applicability of contouring methods, terrain modeling schemes and data visualization techniques to image warping since these too are interpolation/approximation problems.

A distinction may be made between local and global models based on the whether or not the function varies with respect to location, e.g. piecewise linear functions vis-à-vis bivariate polynomials. The former involves a tessellation of the data points, usually a triangulation, and a series of local functions for each element. Bivariate polynomials involve the solution of a single basis, or set of equations, to characterize the entire surface. To summarize, a local transformation is dependent on a subset of nearby data close to the point of evaluation whereas a global method is dependent on all of the data. If the removal of a single datum affects the entire surface, then the method must be considered global in nature.

The nomenclature can be confusing with respect to interpolation and approximation across disciplines. Much of the scattered data literature originates in mathematics and computer science. In a cartographical context, Lam (1983) classifies point (vis-à-vis areal) interpolation as either exact or approximate. This classification carries over into the GIS literature where Burrough (1986) uses the expression "exact interpolator" in the context of interpolation as defined above. The NCGIA Core Curriculum in GIS also uses this terminology (Goodchild and Kemp 1990). If these cases, if a method is not an exact interpolator, then it is an approximating technique.

Back to top

Image Registration with Scattered Data Methods

At a sufficiently large scale, as in the case of typical airborne scanned imagery, the Earth's curvature may be ignored greatly simplifying the distortion modeling. In real Euclidean space, using rectangular Cartesian coordinates, the image registration problem may be separated into two components. Each component consists of fitting a surface to a plane based on a set of control points. The general surface fitting problem may be represented in explicit form as:

Equation 1: z=f(x,y)

In this case, the surface value z is expressed as a function of its location in the plane (u,v). The general image warping or geometric transformation problem may be represented in an equation of the form:

Equation 2: (x,y)=f(u,v)

This states that the (x,y) input coordinates depend on some function of the (u,v) output coordinates. This is expressed as in inverse transformation as it greatly simplifies the resampling process. A forward transformation is more difficult to implement due to problems in the reconstruction and sampling of the DNs (Wolberg 1990). In inverse warping, the input locations are determined based on functions of the output coordinates and the resampling takes place in the input image. It should be clear that resampling—image reconstruction and sampling—does not constitute a form of geometric correction. For example, the following description of "geometric correction" yields little practical information, "the September data were geometrically rectified to the June data using the nearest neighbor algorithm… " (Curran and Pedley 1990). The complete image transformation or warping consists of both a geometric correction and a resampling.

The geometric modeling of the input pixel locations may be calculated separately as follows:

Equation 3: x=f(u,v) and y=f(u,v)

In this form, a wide range of surface fitting methods may be used to model the geometric distortion. Three approaches will be described, each having previously been applied to image warping in recent published literature: (1) inverse distance weighted methods; (2) finite element methods; and (3) bivariate polynomials. Radial basis functions will be discussed in the next section.

In its simplest form, inverse distance weighting is often referred to as Shepard's method (see Franke and Nielson 1991). Its deficiencies are well-known and many modifications have been suggested and reviewed (see Alfeld 1989; Watson 1992). It is a global function that is almost always localized based on the assumption that far away data have little effect on the interpolation. Ironically, the basic methodology survives within the RS/GIS community with a persistence akin to the polynomial approach to image warping. For example, Hodgson, Cheng et al. (1995) apply heuristics and supercomputing power to hurry the construction of this generally unacceptable interpolant. In modified forms, Shepard's method may be suitable for surface fitting and image warping due to its use of ad hoc parameterization (Ruprecht and Müller 1995).

Finite element methods were developed in the engineering literature and have been used for image warping for almost two decades. The simplest form, piecewise linear patches consisting of affine transformations, are constructed based on a triangulation of the control points in the output image. There exist many different triangulation criteria. However, for image warping, no data-dependent or surface-dependent triangulations were located in an extensive literature review.

Two triangulation schemes representing different criteria are the "greedy" (Manacher and Zobrist 1979) and Delaunay (Aurenhammer 1985) triangulations are in The RAND Corporation's CAGIS (Cartographic Analysis and Geographic Information System; Zobrist, Marcellino and Daniels 1991). The Delaunay triangulation has many desirable properties for surface modeling including the avoidance of long narrow triangles and a minimal roughness property in the piecewise linear case (Rippa 1990). In particular, the addition or deletion of a point has only a local effect preventing a rippling throughout the dependent surface as with any other triangulation criteria. The constrained Delaunay triangulation (not strictly a Delaunay triangulation) allows an edge or edges to be fixed (Chew 1989). This makes the triangulation more suitable for applications with "discontinuous" surfaces, e.g. cliffs or cultural features. Adding interior points to a Delaunay triangulation without information on the surface or on nearby points only distorts the rectification of airborne data in an arbitrary manner contrary to research by Devereux, Fuller et al. (1990).

Piecewise linear methods generate surfaces that are continuous but are not differentiable across boundaries, i.e. C0 surfaces. This type of warp has been referred to as tin-sheeting as opposed to rubber-sheeting. C1 surfaces with at least continuous first derivatives across the surface can be more desirable due to a smoother appearance. In remote sensing applications, it will not always be clear if a smoother surface is more accurate as well as being more visually pleasing.

Cubic (Goshtasby 1987) and quintic (Parr and Comer 1990) surface patches are two approaches to improving interpolation quality by including derivative information based on nearby points. The quality of the derivative information greatly influences the surface fit. Data-dependent interpolation is discussed by Dyn and Rippa (1993) and Schumaker (1993) and is becoming increasingly visible.

The methods above interpolate the data. The final method discussed in this section are bivariate polynomials. The following expression is equivalent to an affine transformation when N=1, i.e. a first-order polynomial. The general bivariate formulation may be given as:

Equation 4: Bivariate Mapping Polynomial Formula

where N is the degree or order of the polynomial. These polynomials are equivalent to the trend surface form, although no assumptions are made regarding the existence of any underlying stochastic processes. These polynomials are ill-conditioned, especially at higher orders in which case the solutions must be viewed with skepticism (Ripley 1981). The problem is exacerbated in image-to-map rectification when the output units are quite large, such as might be found in a UTM projection or State Plane Coordinate System.

Modern algorithms are available for solving bivariate polynomial equations such as the Modified Gram-Schmidt, and QR and SVD decompositions (Golub and Van Loan 1989). The pseudoinverse and normal equation approaches should be strictly avoided. Numerical methods are important, especially if RS/GIS investigations are to be conducted under the rubric of science. However, least-squares solutions for bivariate polynomials are well-known and should not require repetition as in Mather (1995). The discussion in Mather on bivariate polynomials perpetuates the almost single-minded focus on bivariate mapping polynomials in remote sensing.

An important and often over-looked characteristic of the least-squares bivariate modeling is its affine invariant form. The coordinate system is arbitrary; the same surface results regardless of scale changes, rotation and translation of the output coordinates. However, the coefficients may be different (Ripley 1981). This information may be used to scale map coordinates, e.g. UTM, and improve the results of the numerical computations. Tensor-product forms as given in Hall (1979) do not possess this property. Similarly deficient is the stepwise selection of parameters for inclusion in the basis (system of equations) if the full trend surface form is not specified.

If the coordinate system arbitrarily affects the surface fit, then any arbitrary origin (translation) and rotation of the output space would effectively change the model. Following Tobler (1994) and Ripley (1981) it would seem prudent to fit coordinate invariant models to the problem or process, deterministic or stochastic. Shepard's method and the greedy and Delaunay triangulations are translation, rotation and scalar invariant, but depending on the triangulation criteria, data-dependent triangulations may be only translation and scalar invariant. In fact, without attention to the physical units in the warp specification, it is possible to generate a myriad of different solutions to the geometric correction problem. This includes local methods such as triangulations as well as global methods. Principally, this refers to expressing coordinates in image space versus map space.

Based on the previous discussion and the survey papers by Franke (1982), Alfeld (1989) and book by Lancaster and Salkauskas (1986), it is possible to identify a set of characteristics that are desirable for scattered data interpolation methods and image warping of remotely sensed data. In particular, the methods must be general—no specific structure to exploit or peculiar configurations to confound—and the data assumed to be relatively accurate. The complexity of the method and its computational requirements must be reasonable. The surface should be smooth—C0 at a minimum and preferably C1 (continuous first derivatives and visually smooth). The method must be coordinate invariant in the sense that the origin shall not influence the surface construction and the implementation of any given scheme must be reproducible. Ideally, the surface should depend continuously on the data. The existence and uniqueness of the solution is presumed. Finally, the effects of free parameters should be understood well-enough to provide acceptable defaults. These criteria are for image warping. The merging of coincident vector coverages, or vector data from a distorted image are different, though related problems. In particular, a forward transformation may be appropriate.

Back to top

Radial Basis Function Methods

Hardy's multiquadrics, and thin plate splines, first appeared in publications on image warping by Göpfert (1982) and Goshtasby (1988), respectively. Both Hardy's multiquadric (MQ) and the thin plate spline (TPS) are global methods and both have been the subject of considerable theoretical and applied research. Other radial basis function methods have been applied to image warping as well, but are not investigate here (see Arad et al. 1994; and Arad and Reisfeld 1995).

In most theoretical work, the functional formulation is very general. For clarity, this presentation is limited to the two-dimensional case. As in the previous section, the surface fitting process is solved separately for image warping in x and y. The radial basis functions may be constructed as a linear combination of the following equations for x and y (presented once for z, an arbitrary variable):

Equation 5: Radial Basis Function Formula

subject to the following constraints:

Equation 6: Radial Basis Function Constraints

where N is number of GCPs and M is the number of terms in a bivariate polynomial as formulated above (referred to as the polynomial precision of the method). The basis for the MQ and reciprocal multiquadrics (RMQ), may be expressed as follows:

Equation 7: Hardy's Multiquadrics and Reciprocal Multiquadrics

where d is the distance between (x,y) to every other point and c is the multiquadric parameter. The effect from a point ui to all equidistant points is the same, expressed as a translate of the radial function h. Originally, Hardy (1971) included no polynomial precision in the formulation. The existence and uniqueness of a solution to the multiquadric equations with or without polynomial precision was proven by Micchelli (1986; see Franke and Nielson 1991). The basis functions are formulations of the upper sheet of a hyperboloid of revolution.

The thin plate spline (TPS) is so-called due to its approximation the least bent infinite, infinitely thin plate or lamina under point loads z at locations (u,v). Its formulation may be given as:

Equation 8: Thin Plate Spline

where the latter form simplifies computations (the a coefficients in (8) are scaled by one-half). This method was introduced by Harder and Desmarais (1972); background references are given in Franke and Nielson (1991). The TPS has linear polynomial precision based on its theoretical development. Matrix formulations of the MQ and TPS are given in Ehlers and Fogel (1994) and Bookstein (1989), respectively. The separation of the radial and monomial (polynomial) bases is given in Lancaster and Salkauskas (1986).

Polynomial precision refers to a method's ability to reproduce a polynomial of a given order. Important in other fields, such as computer aided graphics design, and approximation theory, polynomial precision has less significance in remote sensing. Conceivably, polynomial precision could be important for mosaicking or in localized, adaptive variations of radial basis functions. Carlson and Foley (1991) note that unless the surface may be closely approximated by a low-order polynomial, the MQ basis should not be augmented except perhaps by a constant. Franke (1987) and Fogel and Tinney (1996) found that even a constant may lessen the quality of the MQ fit based on the RMSE on known functions or at known check points. Though not studied here, the completely regularized spline of Mitásová and Mitás (1993) has a constant trend or constant polynomial precision.

It is well-known that the solutions to radial basis functions can be problematic for large numbers of points—large being application specific. The condition numbers for these systems rise with N and inevitably threaten the computations due to rounding error. For image warping with the TPS, Barrodale, Kuwahara et al. (1993) follow Franke's (1982) comments and scale the data to the unit square [0,1]2 greatly reducing the condition number. Franke and Nielson (1991) report the MQ method has a milder condition number. The scaling strategy is known to allow for the solution of problems involving more than 1000 and less than 2000 points. The computational cost is rather small, and unlike properly scaled bivariate polynomials, the coefficients to the scaled problem may be rescaled back to the original problem space.

Hardy's Multiquadrics

In an extensive study by reviewed in Franke (1982), the multiquadric method proved to be one of the best methods for interpolating over a set of different "known" surfaces from (not too) scattered observations. A review of the method with more than 100 references is given by Hardy (1990). The method has only recently been applied to digital image warping although it finds mention in Wolberg (1990). In computer graphics, criteria for good interpolation are different than in remote sensing applications, however the work by Ruprecht and Müller (1993, 1995) is instructive.

Ruprecht and Müller (1993) place a premium on visual and computational aspects in a comparison of a modified version of Shepard's method, the MQ method and the TPS. The radial basis functions were preferred although the parameter adjustment of the modified Shepard's method was more intuitive. In this paper they also use a variable multiquadric parameter. The existence and uniqueness of a solution is no longer guaranteed.

In subsequent work Ruprecht and Müller (1995) implemented a locally-bounded MQ and implemented another free parameter. This made the resulting surface more "tunable" while increasing the computational efficiency and limited the influence of far away points. Notably, the applications are ad hoc, a characteristic of both Shepard's method and its (improved) variants and the MQ method (see Alfeld, 1989). However, the lack of ready alternatives encourages further study. The effects of the multiquadric parameter on interpolation are reviewed in Hardy (1990) and Carlson and Foley (1991).

Göpfert (1982) specified a two-stage approach to image warping which involves estimating the trend using bivariate polynomials and then applying the MQ method to the residuals (residuals are not errors, per se) thus interpolating the GCPs (see Ehlers and Fogel 1994). Fogel and Tinney (1996) found this approach performed better than if low-order polynomial precision had been included in the system of equations, however they provide no results. Fogel and Tinney also noted that the algorithm may prove less useful at extrapolating beyond the convex hull without at least linear (first-order) precision.

Thin Plate Splines

The thin plate spline has received more attention for image warping than the MQ method. In part, this may reflect the well-understood (hence desirable) minimum-curvature properties. It is certainly more intuitive than the MQ method. Image warping using the TPS in computer graphics (Arad et al. 1994; Arad and Reisfeld 1995; and Ruprecht and Müller 1993); and remote sensing (Goshtasby 1988; Barrodale, Kuwahara et al. 1993); Flusser 1993; and Ehlers and Fogel 1994) has generally found favor. Bookstein (1989) provides a decomposition of the TPS mapping from a morphological viewpoint that complements the image warping research.

Flusser (1993) comments on the computational cost of this method and proposes an adaptive approach that depends on the nature of the distortions. However, properties of the thin plate spline have been exploited by Powell (1992a) and implemented by Barrodale, Skea et al. (1993) dramatically reducing the computational requirements with little loss in accuracy. Depending on the nature of the distortions, it appears that the method requires many control points to be effective.

Barrodale's (1995) implementation of Powell's (1992a) scheme transforms entire skeletons or arbitrary shapes as well. This effectively extends the control point-based TPS warp to what may loosely referred to as a feature-based method. However, it differs greatly from mesh-based warping (see Wolberg, 1990) and feature-based schemes (Beier and Neely 1992) in computer graphics. The utility for RS/GIS data integration is obvious given the difficulties in locating coincident points distributed in such a manner as to accurately characterize the distortions.

Bookstein and Green (1993) and Bookstein (1995) build on properties of the thin plate spline to characterize biological shape (biometrics; morphometrics) to include derivative information, called edgels, that may be further used to further specify the deformations. Arad (1994) and Arad and Reisfeld (1995) provide some indication that grid-based patches may be used to refine the warp by coupling the deformation calculations of (u,v) and (x,y). It is not clear to what extent these methods will affect fast evaluation of radial basis function methods such as given in Powell (1992a).

Back to top

Rocky Flats Environmental Technology Site

The motivation for alternative scattered data warping methods follows from the operational shortcomings of the two most commonly used warping methods: bivariate mapping polynomials and piecewise linear finite elements. For many applications, these methods are adequate. However, in a general sense, high spatial resolution airborne scanner remotely sensed data over rapidly varying terrain exposes these models to criticism. The environmental restoration of Rocky Flats provides some context for discussion. The scene and data are described in general terms since the principal goal in this section is to communicate the nature of the rectification problem.

The Rocky Flats facility in Colorado was formerly engaged in nuclear components production. The site is characterized by various contamination problems including, but not limited too ground water contamination. The integration of high spatial resolution remotely sensed data into a GIS for environmental monitoring has a sense of urgency about it. Site specific information and the location of Rocky Flats is shown in Figure 1. More information may be found at the web site for the U.S. Department of Energy Office of Environmental Management (1995).

[IMAGE] Rocky Flats is located approximately 16 miles northwest of Denver, Colorado

Figure 3: Rocky Flats is located approximately 16 miles northwest of Denver, Colorado (North is toward the top of the page).

The site occupies eleven square miles including a buffer surrounding the operating facilities. The elevation varies from over 6200 feet in the southwest to approximately 5600 feet in the southeast corner. Outside the plant inner boundaries, it is difficult to find adequate ground control for registration or rectification. Rectifying data sampled to less than 5 meter-squared to the local Colorado State Plane Coordinate System using bivariate mapping polynomials has a lack of fit at both control points and check points. The ground control is this case is a 2-feet-squared digital orthophoto of the facility and surrounding area. Piecewise finite element methods improve upon the polynomials provided that the GCPs are relatively accurate and representative of scene variations and the tessellation is similarly representative.

In an application using 98 control points and 24 check points, the radial basis function methods provided a better fit compared to polynomials based on two criteria: RMSE and visual inspection. Finite element methods were not included in this comparison due to time and space constraints. Furthermore, it is thought that implementing data dependent triangulations will provide better results and research should be focused in this direction. The polynomial RMSE results for control and check points are presented in Tables 1 and 2, respectively. The bivariate polynomials were applied only to the third order for two reasons: (1) polynomials are prone to oscillate widely between observations for higher degrees, and (2) polynomials are a distraction from research into better alternative techniques.

[IMAGE] Table 1: Bivariate mapping polynomials.  Control Points

Table 1: Bivariate mapping polynomials. Control Points

[IMAGE] Table 2: Bivariate mapping polynomials.  Check Points

Table 2: Bivariate mapping polynomials. Check Points

The implementation of Hardy's multiquadric (Ehlers and Fogel 1994) is a two-stage procedure usually involving the removal of a low order trend and the specification of the multiquadric parameter. Currently, an approximation of the best multiquadric parameter is based on an approximation by Göpfert (1982). The results for the two radial basis functions described previously are presented in Table 3. There are no residuals for these two methods since they are interpolating functions.

[IMAGE] Table 3: Radial Basis Function Methods

Table 3: Radial basis function methods. Hardy's multiquadric (first order polynomial, Göpfert's G=1.10) and the thin plate spline.

The grouping of observations into control and check points is useful for estimating the relative accuracy of the geometric transformations. This measure of accuracy is relative to how well the distortions have been sampled while it also reflects the adequacy of the model for these points. A useful exercise in evaluating both the adequacy of the sampled distortions as well as the overall model fit is ordinary cross validation. A conscience effort is made to avoid the use of the terms "bias" and "predictors" so that this is not confused with statistical methodologies (see Goodall, p. 494). In this case, the "projection" matrices for the polynomials (Hadi 1996) and the matrices for radial basis functions (see Golub, Heath and Wahba 1979) may be readily used to invoke "leave one out" cross validation. In the former case, the calculations are trivial. For radial basis functions, the computations are more intensive and problematic for large numbers of control points. Generalized cross validation is an alternative yet to be explored in this image warping context (for information on generalized cross-validation, see Craven and Wahba 1979 and Bates, Lindstrom, Wahba and Yandell 1987).

The results of ordinary cross validation based a single calculation of the "projection" or "influence" matrices are shown in Table 4. This is the same result one would get by deleting each observation in turn and estimating the model fit without that observation. This cross validated RMSE still suffers as a metric given that, among other things, some areas may have little ground control and consequently reflect poorly on overall model fitness. In general, it is recommended that the ordinary cross validation method be used primarily to identify problems with geometric models and not solely to establish accuracy summaries that are all too often misleading.

[IMAGE] Table 4: Ordinary Cross Validation

Table 4: Ordinary Cross Validation. Third degree polynomial, Hardy's multiquadric (first order polynomial, Göpfert's G=1.10) and the thin plate spline.

An example of the multiquadric method applied to this data set using all points for control is shown in Figure 2. The corrected image is shown with some ancillary data included for visualization (i.e. no comments are made here with respect to the data quality). The image was acquired using a Daedalus 1268 MSS scanner. It is an RGB composite of the thermal band (inverted), an infrared band and a portion of the red part of the electromagnetic spectrum. PC-based ArcView 2.1 (ESRI 1995) is used to display the image and other features.

[IMAGE] Rocky Flats Environmental Technology Site

Figure 4: Rocky Flats Environmental Technology Site. 25 foot contours and roads from GIS coverages are presented over airborne scanned image data (see text).

It is reasonable to question the effects of improved registration on subsequent analyses though this is a much tougher task than the analysis of misregistration effects such as in the work cited previously. RS/GIS integrated analyses will improve measurably in amount and quality. Figures 5 and 6 show two rectified areas from the region using PCIs Image Handler to illustrate the potential of new methodologies.

[IMAGE] Potential Misregistration Example #1

Figure 5: Left-to-right. Third degree polynomial, Hardy's multiquadric (first order polynomial, Göpfert's G=1.10) and the thin plate spline using nearest neighbor resampling.

[IMAGE] Potential Misregistration Example #2

Figure 6: Left-to-right. Third degree polynomial, Hardy's multiquadric (first order polynomial, Göpfert's G=1.10) and the thin plate spline using nearest neighbor resampling.

Back to top

A Software Implementation

The radial basis function methods described here were implemented in several stages. The current software is based on source code for a program called MQE and written by Flottemesch (1993). The program was continually refined over the next sixteen months. Stable and reliable numerical procedures were implemented and some conceptual shortcomings were corrected. At this point, the software was more widely tested. In its current manifestation, MQReg (Multiquadric Registration) uses Flottemesch's code for the user-interface and file-handling. Otherwise it is quite different hence the new name.

The software and a manual will soon be made available via anonymous ftp from the NCGIA (Fogel 1996). Additionally, a revision of a report written is available as an NCGIA Technical Report 96-01 (Fogel and Tinney 1996). The technical report and Ehlers and Fogel (1994) provide a more thorough presentation of radial basis functions and include many more citations.

The program is research code. It is operational, but hardly equivalent to a commercial application. Most notably, it has been structured to be maintainable by others than the programmers. Of course, this means that it must have been written in the C programming language.

At this time, PCI has implemented the thin plate spline (v. 6.0) and ERDAS has been also working on implementing these methods in their commercial packages. Barrodale (1995) has implemented and developed a working demonstration of Powell's (1992a) fast evaluation of thin plate splines. The advanced image rectification software, Spider Warp, can match and warp vector-based objects for highly accurate RS/GIS data integration (Barrodale 1995).

It is expected that more computational details and applications information will be presented in the very near future (Ehlers and Fogel 1996, and McGwire and Fogel 1996). In any case, additional technical information will be deferred until after the general release of this program in the Spring of 1996.

Back to top


The geometric correction of remotely sensed imagery in the absence of ephemeris platform data has been largely ignored by researchers and application specialists. During this period, significant progress has been made in approximation theory, scattered data modeling and surface fitting. The potential utility of these methods is rather obvious, especially given increasingly high spatial resolution data being acquired now and in the future.

Purposely, this presentation has been quite critical of two areas of research: (1) the unfortunate and unwise practice of all but dismissing geometric distortions and their implications in research, and (2) the almost singular focus on bivariate mapping polynomials. This may be regarded as a disservice to those outside the academic and research communities.

The radial basis function methods are intriguing and clearly worthy of increased attention. As Powell (1992b) states, "…the author believes radial basis function methods are becoming as important as piecewise polynomials although good software is not yet available for general computer calculations." It is hoped that this paper will stimulate interest and perhaps provoke debate on this subject within the remote sensing and GIS communities.

Back to top


This research was supported in part by an EG&G Energy Measurements Contract (P.O. 27381) and a NASA Office of University Affairs Grant NAGW-1743 (J.E. Estes, Principal Investigator). Larry Tinney, Michael Howard and Lynne Christel of Bechtel Nevada provided valuable input and data for this paper and Ken McGwire (Desert Research Institute, Reno, NV), and Phil Townsend and Prof. Aaron Moody (University of North Carolina exercised the software used in this paper for their own applications helping to validate these results and making the program more robust. Kim Rauenzahn of TASC deserves all the credit for helping arrange this paper in HTML format (and none of the blame for its current state).

Back to top


Alfeld, P. (1989) Scattered data interpolation in three or more variables. Mathematical methods in computer aided geometric design, T. Lyche and L. L. Schumaker, eds., San Diego: Academic Press, pp. 1-34.
Arad, N., Dyn, N., Reisfeld, D., and Yeshurin, Y. (1994) Image warping by radial basis functions: application to facial expressions. CVGIP: Graphical Models and Image Processing 56(2): 161-172.
Arad, N. (1994) Designing and implementing a grid-distortion mapping based on variational principles. Computer Graphics Forum 13(3): C-259-C-270.
Arad, N., and Reisfeld, D. (1995) Image warping using few anchor points and radial functions. Computer Graphics Forum 14(1): 35-46.
Aurenhammer, F. (1991) Voronoi diagrams--A survey of a fundamental geometric structure. ACM Computing Surveys 23(3): 344-405.
Barrodale, I., Kuwahara, R., Poeckert, R., and Skea, D. (1993) Side-scan sonar image processing using thin plate splines and control point matching. Numerical Algorithms 5: 85-98.
Barrodale, I., Skea, D., Berkley, M., Kuwahara, R., and Poeckert, R. (1993) Warping digital images using thin plate splines. Pattern Recognition 26(2): 375-376.
Barrodale, I. (1995). SpiderWarp. http://www.barrodale.com/~Ubarro. Victoria, British Columbia: Barrodale Computing Services.
Bates, D.M., Lindstrom, M.J., Wahba, G., and Yandell, B. (1987) GCVPACK - Routines for generalized cross-validation. Communications in Statistics. B. Simulation and Computation 16(1): 263-297.
Beier, T., and Neely, S. (1992) Feature-based image metamorphosis. Computer Graphics 26(2): 35-42.
Bernstein, R. (1983) Image geometry and rectification. Manual of Remote Sensing, Second Edition, R. N. Colwell, ed., Falls Church, Virginia: American Society of Photogrammetry, pp. (Chapter 21) 873-922.
Billingsley, F.C. (1982) Modeling misregistration and related effects on multispectral classification. Photogrammetric Engineering and Remote Sensing 48(3): 421-430.
Bookstein, F.L. (1989) Principal warps: thin plate splines and the decomposition of deformations. IEEE Transactions on Pattern Analysis and Machine Intelligence 11(6): 567-585.
Bookstein, F.L., and Green, W.D.K. (1993) A feature space for edgels in images with landmarks. Journal of Mathematical Imaging and Vision 3: 231-261.
Bookstein, F.L., and Green, W.D. (1993) A thin-plate spline for deformations with specified derivatives. Mathematical Methods in Medical Imaging II, San Diego, CA, USA, 07/11 - 07/16/93, SPIE, vol. 2035, pp. 14-28.
Bookstein, F.L. (1995) How to produce a landmark point: the statistical geometry of incompletely registered images. Vision Geometry IV, San Diego, CA, USA, 07/09 - 07/14/95, SPIE -- The International Society for Optical Engineering, vol. 2573, pp. 266-277.
Brown, L.G. (1992) A survey of image registration techniques. ACM Computing Surveys 24(4): 325-376.
Burrough, P.A. (1986) Principles of geographical information systems for land resources assessment. Oxford: Clarendon Press.
Carlson, R.E., and Foley, T.A. (1991) The parameter R2 in multiquadric interpolation. Computers & Mathematics with Applications 21(9): 29-42.
Chew, L.P. (1989) Constrained Delaunay triangulations. Algorithmica 4: 97-108.
Christensen, E.J., Jensen, J.R., Ramsey, E.W., and Mackey, H.E., Jr. (1988) Aircraft MSS data registration and vegetation classification for wetland change detection. International Journal of Remote Sensing 9(1): 23-38.
Craven, P., and Wahba, G. (1979) Smoothing noisy data with spline functions. Numerische Mathematik 31: 377-403.
Curran, P.J., and Pedley, M.I. (1990) Airborne MSS for Land Cover Classification II. Geocarto International 2: 15-26.
Devereux, B.J., Fuller, R.M., Carter, L., and Parsell, R.J. (1990) Geometric correction of airborne scanner imagery by matching Delaunay triangles. International Journal of Remote Sensing 11(12): 2237-2251.
Duggin, M.J., and Robinove, C.J. (1990) Assumptions implicit in remote sensing data acquisition and analysis. International Journal of Remote Sensing 11(10): 1669-1694.
Dyn, N., and Rippa, S. (1993) Data-dependent triangulations for scattered data interpolation and finite element approximation. Applied Numerical Mathematics 12(1-3): 89-105.
Ehlers, M., and Fogel, D.N. (1994) High-precision geometric correction of airborne remote sensing revisited: the multiquadric method. Image and Signal Processing for Remote Sensing, Rome, Italy, 26-30 September 1994, SPIE, vol. 2314, pp. 814-824.
Ehlers, M., and Fogel, D.N. (1996) Radial basis function methods for image rectification: background and implementation (in preparation). ERDAS. (1995) ERDAS/Imagine 8.1. , ERDAS, Inc., Atlanta, Georgia.
ESRI. (1995) ESRI Arcview 2.1. , ESRI, Redlands, California.
Estes, J.E., and Star, J.L. (1993) Remote Sensing and GIS Integration: Towards a Prioritized Research Agenda., Technical Report Number 93-04, National Center for Geographic Information and Analysis, University of California, Santa Barbara.
Fogel, D.N., and Tinney, L. (1996) Image registration using multiquadric functions., Technical Report Number 96-01, National Center for Geographic Information and Analysis, University of California, Santa Barbara.
Forshaw, M.R.B., Haskell, A., Miller, P.F., Stanley, D.J., and Townshend, J.R.G. (1983) Spatial resolution of remotely sensed imagery: a review paper. International Journal of Remote Sensing 4: 497-520.
Franke, R., and Nielson, G. (1980) Smooth interpolation of large sets of scattered data. International Journal for Numerical Methods in Engineering 15: 1691-1704.
Franke, R. (1982) Scattered data interpolation: tests of some methods. Mathematics of Computation 38(157): 181-200.
Franke, R. (1987) Recent advances in the approximation of surfaces from scattered data. Topics in Multivariate Approximation, C. K. Chui, L. L. Schumaker, and F. I. Utreras, eds., Orlando, Florida: Academic Press, Inc., pp. 79-98.
Franke, R., and Nielson, G. (1991) Scattered data interpolation: a tutorial and survey. Geometric modelling: methods and applications, H. Hagen and D. Roller, eds., Berlin: Springer-Verlag, pp. 131-160.
Golub, G.H., and Loan, C.F.V. (1989) Matrix computations, second edition. Baltimore, Md.: Johns Hopkins University Press.
Goodall, C.R. (1993) Computation using the QR decomposition. Computational statistics, C. R. Rao, ed., New York: North-Holland, pp. 467-508.
Goodchild, M.F., and Kemp, K.K., eds. (1990) NCGIA Core Curriculum: Technical Issues in GIS. Santa Barbara, California, 93106: National Center for Geographic Information and Analysis University of California, Santa Barbara.
Goshtasby, A. (1987) Piecewise cubic mapping functions for image registration. Pattern Recognition 20(5): 525-533.
Goshtasby, A. (1988) Registration of images with geometric distortion. IEEE Transactions on Geoscience and Remote Sensing 26(1): 60-64.
Göpfert, W. (1982) Methodology for thematic image processing using thematic and topographic data bases and base-integrated multi-sensor imagery. Proceedings, ISPRS Commission VII Symposium, Toulouse, France, vol. I, pp. 13-19.
Hadi, A.S. (1996) Matrix algebra as a tool. Belmont, California, Duxbury Press.
Hall, E.L. (1979) Computer image processing and recognition. New York: Academic Press.
Harder, R.L., and Desmarais, R.N. (1972) Interpolation using surface splines. Journal of Aircraft 9: 189-191.
Hardy, R.L. (1971) Multiquadric equations of topography and other irregular surfaces. Journal of Geophysical Research 76(8): 1905-1915.
Hardy, R.L. (1990 ) Theory and applications of the multiquadric-biharmonic method . Computers and Mathematical Applications 19 (8/9 ): 163-208 .
Heard, M.I., Mather, P.M., and Higgins, C. (1992) GERES: a prototype expert system for the geometric correction of remotely-sensed images. International Journal of Remote Sensing 13(17): 3381-3385.
Hodgson, M.E., Cheng, Y., Coleman, P., and Durfee, R. (1995) Computational GIS burdens. Geo Info Systems 5(4): 28-37.
Lam, N.S.-N. (1983) Spatial interpolation methods: a review. The American Cartographer 10(2): 129-149.
Lancaster, P., and Salkauskas, K. (1986) Curve and surface fitting: an introduction. New York: Academic Press.
Manacher, G., and Zobrist, A.L. (1979) Neither the greedy nor the Delaunay triangulation of a planar point set approximates the optimal triangulation. Information Processing Letters 9(1): 31-34.
Marble, D.F., and Peuquet, D.J. (1983) Geographic information systems and remote sensing. Manual of Remote Sensing, Second Edition, R. N. Colwell, ed., Falls Church, Virginia: American Society of Photogrammetry, pp. (Chapter 22) 923-958.
Markarian, H., Bernstein, R., Ferneyhough, D.G., Gregg, L.E., and Sharp, F.S. (1973) Digital correction for high resolution images. Photogrammetric Engineering 39: 1311-1320.
Mather, P.M. (1995) Map-image registration accuracy using least-squares polynomials. International Journal of Geographical Information Systems 9(5): 543-554.
McGwire, K. and Fogel, D.N. (1996) Radial basis function methods for image rectification: methdology and applications (in preparation). Micchelli, C.A. (1986) Interpolation of scattered data: distance matrices and conditionally positive definite functions. Constructive Approximation 2(1): 11-22.
Mitásová, H., and Mitás, L. (1993 ) Interpolation by regularized spline with tension: I. Theory and implementation . Mathematical Geology 25 (6 ): 641-655 .
NCGIA. (1991) Initiative 12: Integration of Remote Sensing and Geographic Information Systems, Report of the Specialist Meeting, Technical Report Number 91-16, National Center for Geographic Information and Analysis, University of California, Santa Barbara.
Niblack, W. (1986) An introduction to digital image processing. Englewood Cliffs, N.J.: Prentice-Hall International.
Parr, J.T., and Comer, R.P. (1990) A local area geometric warping algorithm. ACSM/ASPRS Annual Convention, Denver, Colorado, 18-23 March 1990, ASPRS, vol. 4, pp. 316-320.
Powell, M.J.D. (1992a) Tabulation of thin plate splines on a very fine two-dimensional grid. Numerical methods of approximation theory, D. Braess and L. L. Schumaker, eds., Basel: Birkhäuser-Verlag, pp. 105-210.
Powell, M.J.D. (1992b) The theory of radial basis function approximation in 1990. Advances in Numerical Analysis, W. Light, ed., Oxford: Oxford Science Publications, pp. 105-210.
Ripley, B.D. (1981) Spatial Statistics. New York: John Wiley.
Ruprecht, D., and Müller, H. (1993) Free form deformations with scattered data interpolation methods. Geometric Modelling (Computing Supplement 8), G. Farin, H. Hagen, and H. Noltemeier, eds., Wien, Germany: Springer-Verlag, pp. 267-281.
Ruprecht, D., and Müller, H. (1995) Image warping with scattered data interpolation. IEEE Computer Graphics and Applications 15(2): 37-43.
Schumaker, L.L. (1993) Computing optimal triangulations using simulated annealing. Computer Aided Geometric Design 10: 329-345.
Singh, A. (1989) Digital change detection techniques using remotely-sensed data. International Journal of Remote Sensing 10(6): 989-1003.
Star, J.L., ed. (1991) The integration of remote sensing and geographic information systems. Bethesda, Maryland: American Society of Photogrammetry and Remote Sensing.
Thormodsgard, J.M., and Lillesand, T. (1987) Comparison of the gridded finite element and the polynomial interpolations for geometric rectification and mosaicking of Landsat data. Proceedings of the ASPRS-ACSM Annual Convention, Baltimore, Maryland, vol. 2, pp. 139-151.
Townshend, J.R.G., Justice, C.O., Gurney, C., and McManus, J. (1992) The impact of misregistration on change detection. IEEE Transactions on Geoscience and Remote Sensing 30(5): 1054-1060.
Trotter, C.M. (1991) Remotely-sensed data as an information source for geographical information systems in natural resource management. International Journal of Geographical Information Systems 5(2): 225-239.
U.S. Department of Energy. (1995) http://www.em.doe.gov/. U.S. Department of Energy Office of Environmental Management: U.S. Department of Energy Office of Environmental Management.
Van Huffel, S., and Zha, H. (1993) The total least squares problem. Computational statistics, C. R. Rao, ed., New York: North-Holland, pp. 377-408.
Wivell, C. (1995). http://geo.arc.nasa.gov/esd/esdstaff/landsat/landsat.html. Sioux Falls, South Dakota: EROS Data Center.
Wolberg, G. (1990) Digital Image Warping. Los Alamitos, California: IEEE Computer Society Press.
Zobrist, A.L., Marcellino, L.J., and Daniels, G.S. (1991) RAND's Cartographic Analysis and Geographic Information System (RAND-CAGIS): A Guide to System Use. Santa Monica, CA: The Rand Corporation.

Back to top

Author Information

David N. Fogel
Research Fellow
National Center for Geographic Information and Analysis
University of California, Santa Barbara, CA 93106-4060

Back to top