A Comprehensive Study of the

Rational Function Model for

Photogrammetric Processing
C. Vincent Tao and Yong Hu

Abstract uncorrelated because each parameter has a physical signifi-

The rational function model (RFM) has gained considerable cance. Further refinement is also possible by extending the
interest recently mainly due to the fact that Space Imaging model with the addition of calibration parameters to describe
Inc. (Thornton, Colorado) has adopted the RFM' as an alter- effects known or effects suspected to be present.
native sensor model for image exploitation. The RFM has also Although rigorous physical sensor models are more accu-
been implemented in some digital photogrammetric systems rate, the development of generalized sensor models indepen-
to replace the physical sensor mode for photogrammetric dent of sensor platforms and sensor types becomes attractive.
processing. However, there have been few publications addres- In a generalized sensor model, the transformation between the
sing the theoretical properties and practical aspects of the image and the object space is represented as some general func-
RFM until recently. In this paper a comprehensive study of the tion without modeling the physical imaging process. The func-
RFM is reported upon. Technical issues such as the solutions, tion can be of several different forms, such as polynomials or
feasibility, accuracy, numerical stability, and requirements for rational functions.
control information are addressed. Both the direct and itera-
tive least-squares solutions to the RFM are derived, and the Replacement Sensor Models
solutions under terrain-dependent and terrain-independent Sensor models are used extensively in photogrammetric reduc-
computation scenarios are discussed. Finally, evaluations of tion (e.g., stereo reconstruction, ortho-rectification, DEM gener-
the numerous tests with different data sets are analyzed. The ation, feature editing, quality checking, etc.). Many of the above
objective of this study is to provide readers with a comprehen- applications require real-time implementations. Therefore,
sive understanding of the issues pertaining to applications of use of the generalized sensor model to replace the physical sen-
the RFM. sor model for real-time photogrammetric processing is a natural
choice provided that the accuracy degradation is negligible
Background (Paderes et al., 1989).A further advantage is that the photogram-
metric processing software can be kept unchanged when deal-
Physical versus Genefallzed Sensor Models ing with different sensor data because the generalized sensor
Sensor models are required to establish the functional relation- model is sensor independent. For new sensors, only the values
ship between the image space and the object space. They are of of coefficients in the generalized sensor model need to be
particular importance to stereo reconstruction and image updated. The concept of the replacement sensor model is not
ortho-rectification.Sensor models are typically classified into new and has been implemented in some commercial photo-
two categories: physical and generalized models. The choice of grammetric systems (Madani, 1999; Dowman and Dolloff,
a sensor model depends primarily on the performance and 2000; Tao and I-Iu, 2001a).
accuracy required, and the camera and control information The key to the use of a replacement sensor model is that the
available WcGlone, 1996). replacement sensor model must fit the physical sensor model
A physical sensor model represents the physical imaging very well. Normally, the unknown coefficients in the replace-
process. The parameters involved describe the position and ori- ment sensor model are solved using a 3D object grid and the
entation of a sensor with respect to an object-space coordinate corresponding image grid. The parameters of a replacement
system. Physical models, such as the collinearity equations, sensor model are then solved for by fitting the model to the
are rigorous, very suitable for adjustment by analytical triangu- object grid and corresponding image grid.
lation, and normally yield high modeling accuracy (a fraction There are three replacement sensor models that have been
of one pixel). In physical models, parameters are normaIly used (OGC, 1999):the grid interpolation model, the RFM, and
the universal real-time model. In the grid interpolation model,
a 3D grid in ground space is generated and the image coordinates
for each grid point are computed using a physical sensor model.
'The terms Rational Polynomial Camera (RPC)model and Image Geome- To find the image coordinates corresponding to specified
try Model (IGM) are used by Space Imaging in its product line. They
are the same a s the RFM when used in this context.
Department of Geomatics Engineering, The University of Cal-
gary, 2500 University Dr., NW, Calgary, Alberta T2N 1N4, Can-
ada ([email protected]). Photogrammetric Engineering & Remote Sensing
Vol. 67, No. 12, December 2001, pp. 1347-1357.
C.V. Tao is presently with the Geospatial Information and Com-
munication (Geo-ICT) Lab, Dept. of Earth and Atmospheric 0099-lll2/01/6712-1347$3.00/0
Science, York University, 4700 Keele Street, Toronto, Ontario O 2001 American Society for Photogrammetry
M3J 1P3, Canada ([email protected]). and Remote Sensing

I PHOTOGRAMMETRIC ENGINEERING 81REMOTE SENSING December 2001 1347

-
ground coordinates, the surrounding grid points are found. Tri- (X, Y , Z). In order to minimize the introduction of errors during
linear interpolation is then used between these eight points. the computations and improve the numerical stability of equa-
The grid interpolation model does not produce adequate accu- tions, the two image coordinates and three ground coordinates
racy. The RFM uses ratios of polynomials to establish the rela- are each offset and scaled to fit the range from -1.0 to +1.0
tionship between the image coordinates and the object (NIMA, 2000). For the ground-to-image transformation, the
coordinates. The universal real-time model is in fact an exten- defined ratios of polynomials have the following form for each
sion to the RFM. It employs interpolation of high-order correc- image section (OGC, 1999):
tion functions. Because the RFM is the most popular model in
use, the emphasis in this study is placed on the investigation of
the RFM.
Chatacterlstlcs of the RFM
Mathematically, the disadvantage of using polynomials for
approximation is their tendency for oscillation. This often
causes error bounds in polynomial approximation to signifi-
cantly exceed the average approximation error. The RFM has where rn and cn are the normalized row and column indices,
better interpolation properties. It is typically smoother and can respectively, of pixels in image space, and Xn,Yn,and Znare
spread the approximation error more evenly between exact fit the normalized coordinate values of object points in ground
points. The RFM has the added advantage of permitting efficient space. The normalization of the coordinates is computed using
approximation of functions that have infinite discontinuities the following equations (OGC, 1999):
near, but outside, the interval of fitting, while a polynomial
approximation is generally unacceptable in this situation (Bur-
den and Faires, 1997).
The RFM is independent of sensors and platforms. It also
has coordinate system flexibility. It can accommodate object
coordinates in any system such as geocentric, geographic, or
any map projection coordinate system (Paderes et al., 1989).
The RFM resembles projective equations very well
(Madani, 1999). With adequate control information, the RFM where r, and c, are offset values for the two image coordinates,
can achieve a very high fitting accuracy with sufficient speed to and r, and c, are scale values for the two image coordinates.
support real time implementations. This is the primary reason Similarly, X,, Y,, and 4 are offset values for the three ground
why the RFM has been used as a replacement sensor model. coordinates, and Xs, Y,,and Z, are scale values for the three
Because there are no functional relationships between the ground coordinates.
parameters of the physical sensor model and those of the RFM, The maximum power of each ground coordinate is typi-
the physical parameters can hardly be recovered from the W M cally limited to 3; and the total power of all ground coordinates
and the sensor information can be kept confidential. In order is also limited to 3. In such a case, each polynomial is of 20-term
to protect the undisclosed sensor information, some commer- cubic form (the subscripts are omitted for convenience): i.e.,
cial satellite data vendors, such as Space Imaging Inc., only
provide users with the RF'M instead of the physical sensor mod-
els. As a result, without knowing the physical sensor models,
users are still able to perform photogrammetric processing such
as ortho-rectification, stereo reconstruction, and DEM genera-
tion at no discernible loss of accuracy (Grodecki, 2001).
Introduction
The RFM is already available within some digital photogram-
metric software packages (Madani, 1988; Paderes et al., 1989; where aijkare polynomial coefficients called rational function
Greve et al., 1992;Toutin and Cheng, 2000; Yang, 2000). Except coefficients (RFCS). The order of the terms is trivial and may
for some experimental results reported, the research on the differ in different literature.
various aspects of the RFM has not been published. For this rea- The distortions caused by the optical projection can gener-
son, the University of Calgary initiated a research project in ally be represented by the ratios of first-order terms, while cor-
1999, "A Comprehensive Study on the RFM for Photogrammet- rections such as Earth curvature, atmospheric refraction, lens
ric Processing" (Tao and Hu, Zoola). This paper addresses the distortion, etc., can be well approximated by the second-order
important practical issues pertaining to the use of the RFM. terms. Some other unknown distortions with high-order com-
The derivation of the iterative and direct least-squares ponents, such as camera vibration, can be modeled with the
solution to the RFM is described in the next section. The two third-order terms.
computational scenarios, i.e., terrain-independent and terrain- The RFM is essentially a generic form of polynomials.
dependent, are then discussed. For the terrain-dependent sce- When the denominator is equal to 1, Equations l a and l b
nario, where the physical sensor model is not available, a become reguIar 3D polynomials. The RFM resembles the projec-
robust bucketing technique is proposed and developed for tive equations, such as the Direct Linear Transformation (repre-
automatic selection of evenly distributed control points in sented as the ratios of two first-order polynomials). A com-
ground space. Finally, several tests using both aerial photograph parative study of the with other imaging geometric models
and SPOT data sets are described, and evaluations, findings, and is given in Tao and Hu (2000b).
conclusions are given.
lteratlve and Direct LeastSquaresSolutions
Solutions to the Rational Function Model Two methods have been developed to solve for the RFM,direct
Bask Equations and iterative least-squares solutions (Tao and Hu, 2000a). The
In the rational function model, image pixel coordinates (r, c) are derivation of these two solutions is given below.
expressed as the ratios of polynomials of ground coordinates First, we can rewrite Equations l a and l b as

1348 December 2001 PHOTOGRAMMETRICENGINEERING & REMOTE SENSING

 (I Z YX Y3 X3) (a, a, where
r= (44
(1 Z Y X ... Y3X3)' (1 bl "' b1gIT

The observation error equations can then be formed as

1 ZYX --
----... Y3X3-- rY ...
r z -- r
BBBB B B B B

1 Z Y X ... --
----
DDDD D D D D
c z --
Y3X3-- cY ..* --D --D
cy3 cx31 . K --
C

(5b)
W, can be considered as the weight matrix for the residuals on
the left side ofEquation 6a. Consequently,the obtained normal
or
equation is
vi = Bv, = [lZ YX ... Y3 X3 -rZ -rY .-.-rY3 -rX3] . J - r MW,2MJ - MTWfR= 0 (8)
If W, is set to be the identity matrix, the direct solution of
RFCS can be represented as

For the iterative solution, the initial values J(O) of the coeffi-
where cients are first solved for using the direction solution method,
i.e., Equation 9. Then W,(')and Jfi) are calculated by solving the
B = (1 Z YX ... Y3 X3) ' (1 bl '.' bl9IT normal Equation 8 iteratively until the absolute difference of
the residuals between two consecutive iterations is below some
J = (a, al alg b, b, ... b1g)T
threshold. A check for zero crossing in the denominator B is
D = (1 Z Y X . . . Y3 X3) ' (1 dl "' dlJT performed during the computation in order to protect the
weight matrix from division by zero.
K = (c, c1 ... clg dl d2 d19)T Comparing the direct solution, Equation 9, and the itera-
tive solution, Equation 8, the latter is theoretically more rigor-
Given n, the number of ground control points (GCPS) and ous because the weights have been accounted for in the
the corresponding image points, the matrix form of Equation solution. The weight matrix multiplication is introduced to
5a can be written as form the new design matrix of the normal equations in the itera-
tive solution. By making this simplification of weighting, the
direct solution is obtained. Therefore, the computational bur-
den of each iteration is close to the time needed for the direct
solution.
The solution of column-wise Equation 5b is the same as
that of row-wise Equation 5a, provided that the symbols r, B, J,
a, b, R, and M in Equations 5 through 8 are replaced by c, D, K,
c, d, C, and N, respectively. The row-wise and column-wise
equations, Equations 5a and 5b, can be adjusted together. The
following error equations can be used:

V = WTI - WG.

The normal equation is

T W ~ T-ITTW2G= 0. (11)
Similarly, after solving for the the initial values I(O1of the

PHOTOGRAMMETRIC ENGINEERING 81REMOTE SENSING December 2001 1349

coefficientsby the direct solution, W('1 and the iterative solu- range of 0.0002 to 0.004. The result is also not sensitive to the
tion I(fican be computed by solving the normal Equation 11 particular h value as long as h is within that range (see Figure 1).
iteratively. Therefore, h = 0.001 was chosen for the rest of the tests.
The combined form (Equation 11) is used in our experi-
ments, but the results computed using Equation 8 and Equa- Terrain-Independent versus Dependent Computational Scenarios
tion 11 have only a trivial difference. The RFCs of the RFM can be solved for with or without knowing
the physical sensor models. If the physical sensor model is
ReguiarlraUon of the NomlEquatlon available, the terrain-independent solution can be developed.
It has been found that the values of the denominators Biand Di Otherwise, the RFM solution will be highly dependent on the
(i= 1, n) vary widely when the control points are not evenly input control points from the terrain surface.
distributed. As a result, the design matrix T becomes ill condi-
tioned and the matrix T w 2 T in normal Equation 11 could TerrainJndependentComputatlonal Scenario
become singular. It happens often when the high order (i.e., With the physical sensor model available, the RFM can be
more than second-order) polynomials in the RFM are used. The solved using an object grid with its grid-point coordinates
negative impact of this is that the iterative solution cannot be determined using a physical sensor model. This method
converged. involves the following steps:
In order to improve the condition number of the matrix
T W q , we applied the regularization technique in which a Determination of an Image Grid. The image grid contains m by
n image points. These points are evenly distributed across the
small multiplication of the identity matrix E is added. Because full extent of the image. The number of rows m and columns
T W 2is usually~ symmetric and positive semi-definite, the n should be adequate (e.g., 10 by 10).
matrix TrW2T + h2Ehas its eigenvalues in [h2,h2 + 11 Tw2TII] Establishment of a 30 Object Grid in Ground Space. The physi-
and, hence, a condition number is not greater than (h2 + cal sensor model is used to compute the point positions of the
I v T I I ) / h 2 ,which becomes small as h2increases. With this object grid (shown in Figure 2). The dimension of this 3D grid
modification, the normal Equation 11 turns into a regularized is based on the full extent of the image and the range of the
one: i.e., estimated terrain relief, i.e., the dimension of the grid covers
the range of the 3D terrain surface. The grid contains several
elevation layers. To prevent the design matrix from becoming
ill-conditioned, the number of layers should be greater than
three (based on our experience).
To obtain a smoother solution, the above equation is solved RFM Fitting. The unknown WCs are then solved for using the
iteratively using the Tikhonov method (Neumaier, 1998):i.e., corresponding image and object grid points.
Accumcy Checking. A check object grid is generated in a similar
way but with double density in each dimension. The obtained
RFM is then used to calculate the image positions of the check
object grid points. By calculating the difference between the
coordinates of the original image grid points and those calcu-
where W(,) = W(~,I), V(,I = G - TI(:]. lated from the WM, the accuracy of the RFM can be evaluated.
Determination of the regularization parameter h is non- The above method can be used to determine the RFCs of
trivial. Typically, solutions are computed for a large number of each image and provide the estimated accuracy for the RFM fit-
different h values, and the best one is selected by suitable heu- ting solutions, without knowledge of detailed terrain informa-
ristics, e.g., the L-curve based method (Neumaier, 1998). tion. It is thus a solution that is terrain-independent. The
In order to determine the best value of h, different values of
h were tried in Equation 12. An experimental result is shown
in Figure 1, where the third-order RFM was used. In this experi-
ment, 60 GCPs and their corresponding image points were used
to solve for the RFCs, and 249 checking points (CKPS) in an image
were employed for accuracy checking. Based on the tests, it
was found that a reasonably good convergence could be
reached in six iterations or less in most cases when h is in the

-
error (pixel)
7-

6-

.'.......... ' ."

5

4 -
3 - m
m
m
2- m
, ---__--__----------.--------
- mwm

O i , . , , , , , , , , , T , , , , , q , , , , a q ! , , , l h,
0 2 4 6 8 10 20 30 40 50 80 70 BO 90 lOOXlO

Error versus h at CKPs

Figure 2. Use of a 30 object grid and image grid to solve
Figure 1. Determination of regularization parameter h. the RFM (an illustration for the aerial frame sensor).

-
1350 December 2001 PHOTOGRAMM€rRIC ENGINEERING & REMOTE SENSING

-
solved for RFCS can be provided to users for image ortho-rectifi- buckets shown in Figure 3 (from Zhang (1995)).Each bucket
cation, or for stereo reconstruction if a pair of stereo images may contain a number of points, and the buckets having no
each with a set of RFCS is available (Yang, 2000; Tao and Hu, points are excluded. To generate a subset of n GCPS, we select
2001b). This characteristic of the RFM is very attractive because the two points with maximum and minimum elevation values,
users can make use of the images without knowing the physical then randomly select n -2 mutually different buckets and
sensor model that might be treated as confidential information choose one point in each selected bucket randomly.
by commercial vendors. On the other hand, the users do not However, the number of points in a bucket may be quite
need to change the ortho-rectification and stereo reconstruc- different from one to another. There is, consequently, a higher
tion software if the software was developed based on the RFM. probability of selecting a point pertaining to a bucket having
The software can also be used to deal with images acquired by fewer points. In order for each point to have almost the same
different sensors as long as the RFCS of each image are provided. probability of being selected, it is thus preferred that a bucket
having more points should have a higher probability of being
Terrain-DependentComputational Scenarios selected than a bucket having fewer points. This is to ensure
With no physical sensor models at hand, the 3D object grid can- that points selected are evenly distributed. This concept can be
not be established. Therefore, the GCPS and checkpoints have realized using the following procedure (Buckley, 1994;Zhang,
to be collected in a conventional manner (e.g.,from maps or a 1995):if there are in total L buckets, we divide the range [O, 11
DEM).In this case, the solution is highly dependent on the into L intervals such that the length of the i~ interval is equal to
actual terrain relief, the number of GCPS, and their distribution lilSili,where li is the number of points attached to the i~ bucket
across the scene. This method is obviously terrain-dependent. shown in Figure 3. During the bucket selection procedure, a
This is a popular approach for image rectification of remote random number from [O, 11, generated by a uniform random
sensing imagery when the rigorous physical sensor model is generator and falling in the id interval, implies that the id
not available, and when the accuracy requirement is not bucket is selected. Then a point in the i~ bucket is randomly
stringent. selected as a GCP and is excluded from the next selection.
Because the RFM has a very high fitting capability, we are
interested in comparing these two computational scenarios in Robust Selection of GCPs
terms of accuracy and numerical stability. The testing results The bucketing technique only accounts for the planimetric dis-
will be discussed in a later section. tribution of points and not for the vertical (elevation) distribu-
tion. In addition, the generated subset is not unique. Because
Robust Bucketing for Automatic Selection of GCPs the program executes several times, the produced RFM results
For terrain-dependent scenarios, different results can be may be not comparable. In order to determine the best subset,
obtained with different input control information, such as the we could, theoretically, complete trials for all the possible sub-
number of GCPS and their distribution. Given a set of known sets. The best subset is the one in which the RMS (root-mean-
control points, the question is how to select a group of control square) error of the solution at checkpoints is minimum. For
points by which the RFM can reach the best overall fitting accu-
racy. In other words, how do we select a group of points that are
evenly distributed in each dimension? To address this ques-
-
example, given 71 input points, the number of subsets each
with 64 points is C!f 1.33 X lo9. This is computationally
impossible in practice, especially when the number of input
tion, we developed a robust bucketing method for automatic points is large. Thus, we choose t(s1)trials in such a way that
selection of GCPS from a given set of control points. the ~robabilitvof at least one subset with evenlv distributed
poiits will beg large value A close to 1. ~ e c a u s i t h entire
e set
Bucketing Technique has been divided into L buckets, the expression for this proba-
This method is based on a regular random selection bucketing bility is
technique (Zhanget al., 1995).The idea of selecting a subset of
points h a t are evenly distributed is not to select more than one
point for a given neighborhood. The bucketing technique is A = 1 - (1 - d'lL)* (14)
given as follows:
First, the minimum and maximum of the planimetric coor- where a is the probability measuring the extent to which the
dinates of the entire point set in the ground space are calcu- entire set distributes evenly; then cu"JL is the probability of
lated. Then the test region is divided into w X w (e.g., w = 8) obtaining a perfect subset by bucketing. In our experiments, a

I
I
I
I
f
I
I
I
I ,.
',"
, ',* m; ; I uniform
0 1 2 ... w-l 0 1 variable

Figure 3. Automatic selection of GCPS based on the bucketing technique.

PHOTOGRAMMETRIC ENGINEERING 81REMOTE SENSING December 2001 1351

 Test Data Sets

Aerial Photograph Data [provided by ERDAS Inc., USA)

The original stereo image pair at a scale of 1:40,000was taken
over Colorado Springs, Colorado, and was scanned at 100
micrometers per pixel. Within the overlap area of the stereo
pair, 7499 well distributed points in the left image and their
corresponding points in the right image were determined using
an automatic image matching algorithm. The ground coordi-
nates of these 7499 points were intersected using the rigorous
collinearity equations established after a photogrammetric
bundle block adjustment with ERDAS OrthoBASE. The average
is defined as Ll(w x w) and is usually avalue close to 1.Then t standard deviations after adjustment in ground space are (mx,
can be computed using my, mZ)= (1.7008,2.1577,0.2957) meters at five control points,
and (mx, my, mZ) = (4.2964,0.7726, 3.8165) meters at one
checkpoint. The above processing was done by ERDAS Inc. The
left image is used in the experiments. Figure 4 shows a 3D view
of the terrain surface generated using these ground points.
Other relevant information is given in Table 3.
This method is a simple variation to the random sample
consensus paradigm (Fischler and Bolles, 1981).As shown in SPOT Data Set I (provided by ERDAS Inc., USA)
Table 1, the number of bucketing trials for given A and a, when
n and ware known (e.g., n = 64; w = a), is much less than the The original image pair is part of the example data in ERDAS
possible combinations even if a very high reliability is ensured. IMAGINE OrthoBASE (Yang,2000). They were taken over the
For example, when a is 80 percent, only four trials are needed Palm Springs, California region. The relief range is from - 75.0
to reach a probability as high as 0.99. meters to +3414.0 meters. The overlap between the two images
is about 60 percent. The image sizes are 6000 by 6000 pixels
with 13 micrometers for a pixel, and the ground pixel size is 10
Test Results and Evaluation meters. The physical sensor parameters of both images were
derived using the OrthoBASE block triangulation software.
Design of Tests The unit-weight RMS errors are 0.3 pixels. The grid consists
In order to provide a comprehensive evaluation of the RFM, the of five layers each with 11by 11points. These coordinates
tests have been designed for these purposes:
were calculated using a rigorous physical sensor model in
to evaluate the overall fitting accuracy of the RFM under nine OrthoBASE by considering Earth curvature. An independent
different cases, i.e., three different polynomial orders against set of 21 by 2 1 checkpoints on the actual terrain surface was
three different denominator configurations (see Table 2); generated based on the inverse form of the corresponding rig-
to compare the RFM performance resulting &om the two scenar- orous model of the left image (Yang, 2000). Part of the informa-
ios, terrain-independent scenario versus terrain-dependent; tion is shown in Table 3, and a 3D view of the distribution of the
to compare the RFM results derived using the direct least-squares
441 checkpoints is shown in Figure 5.
solution and the iterative least-squares solution;
to evaluate the numerical stability of the RFM solutions as well
as the effects of regularization method; SPOT Data Set I1 (provided by PCI Geomatics, Canada)
to evaluate the RFM results based on the GCPS selected using the
developed automatic method in the case of terrain-dependent In total, 71 ground control points and the corresponding image
scenario; and points were collected manually. Figure 6 shows a 3D view of
to evaluate the sensitivity of the RFM to the input control infor- the distribution of these points on the terrain surface. Most of
mation (number and distribution of GCPS).
Table 2 shows the number of unknown WCs and the
required minimum number of m s under these nine cases. It is
.....6.<. ..,
. ..
worth noting that the case of p2 = p4 = 1is a regular 3D polyno-
. .
, , "'i.... . :
mial model, and the case of p2 = p4 with the first order is the .. . ' .. , : ' ....
.. .... .
Direct Linear Transformation (DLT) model. For the regular 3D
polynomial model, only the direct least-squares solution is
applied.

TABLE2. NINECASESFOR THE RFM TESTS

Order of Number of Min. Number
Polynomials Cases RFCs of GCPs

Figure 4. 3D view of the terrain surface from the aerial phote

graph data.

1352 December 2001 PHOTOGRAMMETRIC ENGINEERING 81REMOTE SENSING

 TABLE3. I N ~ R M A T ~ O N ON THE THREEDATASETS
Image Size Ground Pixel Size Film Pixel Size Elevation Range
Data Set (pixel) (meters) (mfi) Location (meters)
Aerial Frame 2313 X 2309 4.5 Colorado Springs, CO 1846.6 - 2205.1
SPOT I 6000 X 6000 10 Palm Springs, CA -75.0 - 3414.0
SPOT I1 6000 X 6000 10 N/ A -
1.8 1128.3

using the rigorous sensor model established. The RFCSwere

. .
., ........: ... . then solved for using the control grid coordinates and the cor-
responding image coordinates. The accuracy of the solution
was examined by comparing the differences in the image coor-
dinates of the check grid computed using the generated RFM
and using the rigorous sensor mode!. Tests under nine different
cases (Table 2) were performed. Both the direct and the iterative
.. . .. .. ..
solutions were also used. A portion of the experimental results
using the iterative solution without regularization are given in
Table 4. Both the RMS errors and the maximum errors in the
image at the control points (CNPS) and the checkpoints (CKPS)
are listed. This result is based on the established control grid
consisting of five elevation layers each with 10 by 10 grid
points, and the check grid consisting of ten elevation layers
each with 20 by 20 grid points.
Based on this test, the following findings were obtained:

-
The approximating accuracy of the RFM with denominators
(cases p2 # p4 and p2 = p4) is extremely high and is signifi-
cantly better than the regular polynomial cases (p2 = p4 1).
Figure 5. 3D view of the distribution of checkpoints from the The accuracy reaches the level of lo-' to pixels in compar-
SPOT data I. ison to the regular polynomial cases where the accuracy level
is at 10-I to pixels for the second- and third-order polyno-
mials. This shows that the RFM with rational components
(denominators) resembles the collinearity equations well and
has no accuracy loss for the aerial photograph data.
. .
..' ..-...
. .
... ..... . ..
The cases of first-order RFM with denominators gain the best
results at both CNPS and CKPS (shaded area in Table 4). It indicates
..'. . ..
,

. ,.1'
. .
. . .. .... that the higher order RFM may not be necessary when dealing
with aerial frame data sets.

Accuracy Aspect-SPOT Data I

In this test, the provided 3D object grid was used to solve the
RFM, and 441 checkpoints from an actual DEM were used for
accuracy checking. The test results obtained using the direct
solution without regularization are provided in Table 5. The
following findings can be summarized:
The approximation accuracy is also extremely high for the cases
when the second- or the third-order polynomial is used (shaded
area in Table 5). The accuracy can reach at CKPs. The
maximum errors are also controlled under 0.06 pixels.
None of the first-order cases perform well compared to the aerial
photography test in which the cases (first-order RFM with
Figure 6. 3~ view of the distribution of points from the SPOT denominator) gain the best results (see shaded areas in Table 4).
data 11. The RFM cases with denominator always perform better than
the regular polynomials but the third-order regular polynomial
model can reach the same level of accuracy at CKPS.
The RFM cases with unequal denominator (p2 # p4) perform
better than the cases with equal denominator (p2 = p4).
the points are situated on a flat foothill, while only a few
appear in the mountain area. Other information is given in
Table 3. The accuracy of point measurements is not known. Computation Aspect of Terrain-Independen t Scenario

Experiments and Evaluation : Tenaiblndependent Scenario

In the above two tests, an object grid was used to solve for the
RFM. It was found that the condition of the design matrix in
Accuracy Aspect-Aerial Photograph Data normal Equation 11is always good and regularization of the
normal eauation is not necessarv. The accuracv obtained with
A rigorous physical sensor model based on the collinearity the iterathe solution is only sliihtly better than that obtained
equations was developed using the imaging parameters pro- with the direct solution. Therefore, the direct solution method
vided by ERDAS Inc. An image grid as well as its corresponding is recommended due to its simplicity and no requirements for
control grid and check grid in object space were generated iterations.

PHOTOGRAMMETRIC ENGINEERING 81REMOTE SENSING December 2001 1353

 TABLE
4. RMS (MAX). ERRORSIN IMAGE WITH THE AERIALPHOTOGRAPH
DATA(UNIT:
PIXELS)
-

First-Order Second-Order Third-Order

Case at CNPs at CKPs at CNPs at CKPs at CNPs at CKPs

TABLE5. RMS (MAX.) ERRORS

IN IMAGE WITH SPOT DATAI (UNIT: PIXELS)
First-Order Second-Order Third-Order
Case at CNPs at CKPs at CNPs at CKPs at CNPs at CKPs

The number of layers used in the object grid affects the RMS C5.0e - 02 at CKPS),
if the control points are adequate and
numerical stability of the RFM. The design matrix becomes sin- evenly distributed for the aerial frame data. However, the level
gular and the normal equation is ill-conditioned if only two or of approximation is not as high as that in the terrain-indepen-
three elevation layers in the object grid are used. However, the dent scenario.
The WM cases with denominator are better than the regular
approximation accuracies are very close when four or more lay- polynomial models overall, but the accuracy of the third-order
ers are used. regular polynomial model is still very good.
It is interesting that, from Table 6, the case p2 f p4 always has
Experiments and Evaluation :Terraln-Dependent Scenario a better fitting accuracy than the case p2 = p4 at CNPS for
different orders, but the case p2 = p4 works better than the
Accuracy Aspect-Aerial Photograph Data case p2 # p4 at cKPs. However, there is no significant difference
between both cases.
In the terrain-dependent scenario, the GCPs are from the real ter- The first-order RFM cases (shaded areas in Table 6) gain higher
rain surface. For a comparison, we used the same data set, accuracy at CKPS,while the third-order RFM cases have higher
aerial photograph data from ERDAS, in which 7499 ground accuracy at CNPS. None of these cases can achieve the best accu-
points onthe terrain surface are available. In this test, 100points racy at both CNPS and CKPS. This implies that the RFM approxi-
out of 7499 points were first selected as GCPs using the auto- mation has not achieved the best numerical status under the
matic bucketing selection method, and the remaining points terrain-dependent scenario. Therefore, for terrain-dependent
solutions, one may need to try different orders in order to deter-
are treated as checkpoints. Due to the good distribution of these mine the best RFM case or the best trade-off.
points, only one trial was needed in the bucketing selection
process. As described earlier, this selection ensures that these
100 GCPS are evenly distributed from the given 7499 points. Computation Aspect-Aerial Photograph Data
In order to provide a reliable assessment on accuracy, the
rigorous collinearity equations were used to transform these Due to the fact that the design matrix is almost rank-deficient
7499 points to both the left and right images. The obtained under the terrain-dependent scenario, regularization is particu-
accuracy of the RFM solution based on the 100 GePs and the larly important especially for the cases with the second- and
7399 checkpoints, as well as their corresponding image points, third-order polynomials. Table 7 shows a comparison of results
is shown in Table 6. The results are computed from the left solved for with and without the use of regularization (Table 6
image using the direct least-squares solution with regulariza- is the result with regularization). The improvements are very
tion ( h = 0.001). significant in terms of accuracy. It is worth mentioning that
regularization also makes the solutions converge. This result
This result shows that, under the terrain-dependent scenario, indicates that the numerical stability of the RFM may be poor
the RFM can also achieve a high approximation accuracy (i.e., under the terrain-dependent scenario, although it is able to

TABLE6. RMS (MAX.) ERRORSIN IMAGE WITH FRAME

DEM DATA
(TERRAINDEPENDENT
SCENARIO,
UNIT: PIXELS)
First-Order Second-Order Third-Order
Case at CNPs at CKPs at CNPs at CKPs At CNPs At CKPs

PHOTOGRAMMETRIC ENGINEERING & REMOTE SENSING

 TABLE7. EFFECTSOF REGULARIZATION
FOR ACCURACY IMPROVEMENTS (UNIT:
PIXELS)
Case Order Regularization RMS at CNPs Max at CNPs RMS at CKPs Max at CKPs
1 h =0 4.2924e-02 1.4294e-01 4.9006e-02 1.8375e-01
h = 0.001 4.2924e-02 1.4295e-01 4.8997e-02 1.8368e-01
2 h=O 9.2097e-02 5.5913e-01 1.5101e+00 7.5617e+01
h = 0.001 3.9456e-02 1.3810e-01 5.1097e-02 2.2157e-01
3 h=O 7.0171e-02 2.5211e-01 6.3360e+00 5.3717e+02
h = 0.001 3.4249e-02 9.8376e-02 6.1061e-02 4.1028e-01
1 h=O 4.3339e-02 1.3759e-01 4.8654e-02 1.6661e-01
h = 0.001 4.3340e-02 1.3756e-01 4.8646e-02 1.6657e-01
7 h=O 7.2629e-02 3.3755e-01 2.5887e+00 1.8341e+02

.
l
h = 0.001
h 0
h = 0.001
- 4.0481e-02
8.2221e-02
3.5615e-02
1.3900e-01
4.4792e-01
1.0335e-01
5.0317e-02
3.1974ef00
5.5134e-02
2.0361e-01
2.4258e+02
2.8847e-01

achieve a high approximation accuracy provided the ccps are The third-order case p2 # p4 gains the best approximation
adequate and well distributed. accuracy at the CWS,while the regular 3D polynomial model
For this data set, due to the adequate number of well-distributed reaches the same level of accuracy at the CWS.
ccps, the accuracy obtained by the direct solution is comparable The cases with denominator are better than the regular polyno-
to that from the iterative solution. Therefore, only the results mial cases (p2 = p4 = 1).
from the direct solution are given in Table 7. The cases with non-equal denominator (p2 # p4) are always
Figure 7 shows the relationship between the number of GCPS better than the cases with equal denominator (p2 = p4). This
used (from 40 to 1600) and the accuracy obtained at CKPS. In conclusion is the same as that from the SPOT data I test under
Figure 7, the computation result is based on the third-order RFM the terrain-independent scenario.
with unequal denominator case where the required minimum For this data set, the higher the order, the better the approxima-
number of W s is 39 (see Table 2). It can be observed from tion accuracy.
Figure 7 that the accuracy will not be improved much once
the number of ccps used is larger than 80 (which doubles the
required minimum number, 39). This is a useful finding
because the use of too many ~s would not be beneficial to Computation Aspect-SPOT Data ZI
the accuracy. It also indicates that the proposed automatic
selection method of ccps is efficient. The iterative least-squares solution with regularization pro-
vides better results than does the direct soIution. This result
shows that the iterative solution with regularization is much
Accuracy Aspect-SPOT Data I1 more robust. It performs best for data sets in which the normal
Sixty-four out of the total given 71 GCPS were selected automati- equations are not well-conditioned or GCPs are not well distrib-
cally using the robust bucketing method, and the remaining uted. More comparative experiments have been reported in Tao
points were used as checkpoints. It is a fact that the distribution and Hu (2000a; 2000b).
of these GCPS (see Figure 5) is not good. The obtained a is only
56.3 percent (see Equation 14) and the number of trials is seven Conclusions and Discussions
when the probability A is assumed to be 95 percent. The results In this paper, the iterative and direct least-squares solutions to
shown in Table 8 were computed using the iterative least- the WMare derived. The regularization technique is proposed
squares method with regularization (h = 0.001). The following to improve the solutions when the condition of the normal
findings were obtained: equations is poor. The terrain-independent and terrain-depen-
dent scenarios are analyzed and compared. Based on the
Because the available ccps are not evenly distributed, the
obtained accuracy is only at the one-pixel level. This is under- numerous tests, the following conclusions can be drawn:
standable because the RFM solution is sensitive to the distribu- With the physical sensor model available-terrain-inde-
tion of GCPS under the terrain-dependent scenario. pendent scenario:
The RFM can approximate the rigorous physical sensor model
extremely well for aerial frame data and SPOT pushbroom data

-
used in the tests. Thus, it can be used as a replacement sensor
model to communicate to end users for photogrammetric proc-
0.3 -error- (pixel)
- -5 essing such as ortho-rectification and stereo reconstruction.
RMSE The high order RFM may not be necessary when dealing with
-MA)( -- 4 aerial frame data sets, because the RFM resembles the projective
............................... -...-.. ........ model very well. However, the high order RFM is preferable
-- 3 when SPOT data are processed.
..... .-. The RFM cases with unequal denominator achieve a better accu-
-- 2 racy than the cases with equal denominator at cws overall.
... .. The RFM cases with denominator perform better than the regular
polynomials, but the regular third-order polynomial model can
. =. ......w:"
1 achieve an accuracy that is also reasonably good.
0. r . .
*---*-*-*-*-*-*-*-4
. . . . . . .
7 r v o
The conditions of normal equations are very good. The iterative
least-squares solution is only slightly better than the direct
40 50 60 80 100 120 160 200 300 403 800 1M)O 1800 . least-squares solution. There is no need to apply regularization.
no.ofaFs The direct solution without regularization is recommended for
solving the RFM due to its simplicity and no iterations involved.
Figure 7. The relationship between the accuracy and the In the establishment of the 3D object grid for the RFM solutions,
number of GCPs. at least four or more elevation layers are needed, otherwise,
the design matrix would be singular.

I PHOTOGRAMMETRIC ENGINEERING & REMOTE SENSING December 2001 1355

 TABLE
8. RMS (MAX.)ERRORSIN IMAGE WITH SPOT DATAII (TERRAIN-DEPENDENT
SCENARIO,
UNIT:
PIXELS)
-- - - -

First-order Second-order ~hird-order

Case at CNPs at CKPs at CNPs at CKPs at CNPs At CKPs

With the physical sensor model unknown-temin-depen- number of coefficients. This causes instability i n the least-
dent scenario: squares solution. Madani (1999) proposed selecting the signifi-
cant coefficients for a particular sensor by the trial-and-error
The RFM can also achieve a good approximation accuracy pro- method. In this paper w e used the regularization technique
vided the GCPS are adequate and are evenly distributed. How- with some success to overcome this drawback. Further investi-
ever, the accuracy of the approximation is not as high as that gations on the accuracy analysis of high resolution images
in the terrain-independent scenario. Moreover, in practice, it
is not easy to collect an "adequate number" of GCPS. (such as IKONOS-2 imagery) using the RFM approach are of inter-
The RFM solution is sensitive to the distribution of GCPS as well est (Tao and Hu, 2001b).
as to the number of GCPs. The result cannot be expected to be
good once the distribution of the GCPS is not good (see tests on Acknowledgments
SPOT Data 11). Discussions with Drs. Xinghe Yang and Yongnia Wang from
It is difficult to determine which case of the RFM can provide ERDAS Inc., Dr. Philip Cheng from PC1 Geomatics, Dr. Thierry
the best result, because the solution is not numerically stable. Toutin from the Canada Centre for Remote S e ~ s i n gMr.
, Brian
Therefore, one may need to try different cases to determine the Robertson from MacDonald Dettwiler and Associates Ltd., Dr.
best RFM case (i.e., order) or to find the best trade-off. Clive Fraser from the University of Melbourne, and Prof. Gor-
The design matrix is almost rank-deficient;therefore, regulariza- don Petrie from the University of Glasgow were extremely val-
tion is extremely important, especially for the high order RFM uable. We also appreciate the comments made by Dr. Mostafa
cases. Use of regularization to improve the conditions of the
design matrix will result in a much better accuracy and make Madani from Z/I Imaging, Dr. Ian Dowman from University
the solutions converged. College London, Mr. Kris Morin from LH systems, as well as
The iterative least-squares solution always offers better results anonymous reviewers. Special thanks go to ERDAS Inc. (Drs.
than does the direct solution, particularly for the terrain-depen- X. Yang and Y. Wang), PC1 Geomatics (Dr. P. Cheng), and
dent scenario.Therefore, the iterative least-squaresmethod with Intermap Technologies (Dr. J. Bryan Mercer) for providing test
regularization should be adopted for solving the RFM under the data sets to support the research.
terrain-dependent scenario.
The proposed method on automatic selection of GCPS is useful
for improving the numerical stability of RFM solutions because References
the method can produce from the given control points a set of Buckley, J.J., 1994. Fuzzy genetic algorithm and applications, Fuzzy
GCPS that are evenly distributed. Sets and Systems, 61(2):129-136.
Burden, R.L., and J.D. Faires, 1997. Numerical Analysis, Sixth Edition,
Consequently, the RFM determined under the terrain-inde- BrooksICole Publishing Company, Pacific Grove, California,
pendent scenario can offer a robust a n d very accurate approxi- 811 p.
mation to the rigorous physical sensor model. The RFM Dowman, I., and J.T. Dolloff, 2000. An evaluation of rational functions
soIutions derived from the terrain-dependent scenario are sen- for photogrammetric restitution, Int'l Achieve of Photogrammetry
sitive to the number of GCPS and their distribution, and are not and Remote Sensing, 33(Part B3):254-266.
numerically stable. Based o n the concept that the RFM can be Fischler M.A., and R.C. Bolles, 1981. Random sample consensus: A
used to deal with different sensor imagery, a prototype system, paradigm for model fitting with applications to image analysis
Rational Mapper, has been developed for the RFM-basedortho- and automated cartography, Communications of the ACM,
rectification and stereo reconstruction using images whose 24(6):381-395.
WCs are provided (Tao and Hu, 2001a; Tao and Hu, 2001b). The Greve, C.W., C.W. Molander, and D.K. Gordon, 1992. Image processing
advantage of the system is that the software is sensor indepen- on open systems, PhotogrammetricEngineering & Remote Sens-
dent and the photogrammetric processing becomes versatile. ing, 58(1):85-89.
It is a fact that the RFM solutions are usually determined by Grodecki, J., 2001. IKONOS stereo feature extraction-RPC approach.
the data vendor using a proprietary physical sensor model. The Proceedings of 2001 ASPRS Annual Convention (CD ROM), 23-27
accuracy of the RFM solutions is dependent on the availability April, St. Louis, Missouri, unpaginated.
and the usage of the GCPS (Grodecki, 2001). If accurate RFM Hu, Y., and C.V. Tao, 2001. Updating solutions of the rational function
solutions are required, GCPs are needed and are incorporated model using additional control points for enhanced photogram-
into the RFM solution process. In many cases, GCPs are not metric processing, Proceedings of Joint ISPRS Workshop on "High
available at the time of processing or cannot be supplied due to Resolution Mapping from Space 2001 " (CD ROM), 19-21 Septem-
some reasons (e.g., politics or confidentiality). H u and Tao ber, Hannover, Germany, unpaginated.
(2001) have investigated a method to update and improve the Madani, M., 1999. Real-time sensor-independent positioning by
existing RFM solutions (provided, for example, by the vendor) rational functions. Proceedings of ISPRS Workshop on "Direct
versus Indirect Methods of Sensor Orientation", 25-26 November,
with additional GCPs using the Kalman-filter-based incremen- Barcelona, Spain, pp. 64-75.
tal technique and covariance propagation. McGlone, C., 1996. Sensor modeling in image registration, Digital Pho-
Because low-order RFM cases may obtain better approxi- togrammetry: An Addendum (C. W. Greve, editor), American Soci-
mating accuracy, a remaining problem with the high-order RFM ety for Photogrammetry and Remote Sensing, Bethesda, Maryland,
cases is their over-parameterization because of an excessive pp. 115-123.

PHOTOGRAMMETRIC ENGINEERING 81REMOTE SENSING

 Neumaier, A,, 1998. Solving ill-conditioned and singular linear system, 2001a. The rational function model: A tool for processing high-
SIAM Review, 40(3):636-666. resolution imagery, Earth Observation Magazine (EOM),
NIMA (National Imaging and Mapping Agency), 2000. The Compen- 10(1):13-16.
dium of Controlled Extensions (CE) for the National Imagery 2001b. 3-D reconstruction algorithms with the rational function
7kansmission Format (NITF), Version 2.1, URL: http:/l model and their applications for IKONOS stereo imagery, Proceed-
www,ismc.nima.mil/ntb/superceded/STDI-0002~v2.l.pdf ings of Joint ISPRS Workshop on "High Resolution Mapping from
OGC (OpenGIS Consortium), 1999. The OpenGIS Abstract Specifica- Space 2001" (CD ROM), 19-21 September, Hannover, Ger-
tion-Topic 7: The Earth Imagery Case, URL: http:// many, unpaginated.
www,opengis.org/public/abstract/99-107.pdf. Toutin, T., and P. Cheng, 2000. Demystification of IKONOS, Earth
Paderes, Jr, F.C., E.M. Mikhail, and J.A. Fagerman, 1989. Batch and Observation Magazine (EOM), 9(7):17-21.
on-line evaluation of stereo SPOT imagery, Proceedings of the Yang, X., 2000. Accuracy of rational function approximation in photo-
ASPRS-ACSM Convention, 02-07 April, Baltimore, Maryland, grammetry, Proceedings of 2000 ASPRS Annual Convention (CD
pp. 31-40. ROM), 22-26 May, Washington D.C., unpaginated.
Tao, C.V., and Y. Hu, 2000a. Investigation on the Rational Function Zhang, Z., R. Deriche, 0. Faugeras, and Q.T. Luong, 1995. A robust
Model, Proceedings of 2000 ASPRS Annual Convention (CD technique for matching two uncalibrated images through the
ROM), 24-26 May, Washington, D.C., unpaginated. recovery of the unknown epipolar geometry, Artificial Intelli-
gence, 78(1-2):87-119.
-, 2000b. Study of the rational function model for image rectifica-
tion, Proceedings of Canadian Conference on Remote Sensing, (Received 07 February 2001; accepted 10 May 2001; revised 17
22-25 August, Victoria, BC, Canada, pp. 55-64. August 2001)

M A X I M I Z E YOUR COMPANY'S EXPOSURE

BY HAVING YOUR IMAGES FEATURED I N ASPRS PROMOTIONAL MATERIALS!

As technology is constantly changing, it 1s important for the society to provide

the most current representations of this industry to our members, readers, and po-
tential clients or employees.
We need your help with locating applicable images that convey exceptional ex-
amples of the work done by our profession so that we can expand the variety of im-
ages in our photo image bank. We are asking all ASPS Sustaining Members to par-
ticipate by donating five to ten of their company's images illustrating good
examples of GIs, Photogrammetry or Remote Sensing.
The images would be featured throughout ASPRS promotional materials - prima-
rily conference brochures and advertisements, membership pieces, and the ASPS
Career Brochure and accompanying materials.
In exchange for your donation, your company will receive photo recognition
wherever the image is used, thereby increasing your company's visibility. We en-
courage you to take advantage of this opportunity-please participate.
To donate images, please follow the submission instructions below. For more in-
formation, please contact ASPRS Production Manager, Rae Kelley by phone: 301-
493-0290 ext. 107 or email: [email protected].

Submission Instructions
Color
300 DPI, tiff or eps format
40-50 word description including image type, GIS, Photogrammetry or Remote Sensing
Submit on a disk or to the ASPS FTP site (Contact Rae Kelley for logon instructions.)

Thank you in advance for your support!

Donate an image today and it may be featured in the ASPRS Career Poster!
The poster will accompany the ASPRS Career Brochure and serve as an additional marketing piece to help
promote the education initiative. The poster will be given to high school teachers and college professors
teaching subjects relative to our industry.

PHOTOGRAMMETRIC ENGINEERING & REMOTE SENSING December 2 0 0 1 1357