David N. Fogel
![[IMAGE] Rocky Flats Sub-Area Map](rf-mapt.gif)
![[IMAGE] Rocky Flats Sub-Area Image](rf-mqrt.gif)
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](sf-0r-th.gif)
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](sf-1g-th.gif)
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.
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).
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.
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.
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:
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:
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:
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.
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):
subject to the following 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:
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:
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).
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](rf-site.gif)
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](tab-1.gif)
Table 1: Bivariate mapping polynomials. Control Points
![[IMAGE] Table 2: Bivariate mapping polynomials. Check Points](tab-2.gif)
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](tab-3.gif)
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](tab-4.gif)
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](rf-mqr.gif)
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](dns-ex1.gif)
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](dns-ex2.gif)
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.
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.
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.