29- 
the same formula is obtained. It is obvious from eq. (5b) that in each grid 
element the interpolated surface is represented by a bicubic polynomial. It also 
follows from the formulation that the local surfaces agree at their boundaries in 
height and in tilt. 
To compute the as yet unknown heights at the grid nodes, each reference 
point provides one equation (5b). These equations, if too small in number, are 
augmented by condition equations which express regularity of the total surface and 
minimum curvature. The latter are used with a generally very small weight and the 
heights at the grid nodes are solved from the equations by the method of least 
Squares. 
Because in the methods of this group the number of parameters that must 
be computed during the first step is equal to the number of grid nodes or is a 
multiple thereof, these methods are most attractive in the case where the number 
is smaller than the number of given reference heights. Bosman et al give one 
application in which a region with 234 reference values [50] or maybe even 312 
values [22] is subdivided into only five rectangular units. In one of Junkins 
and Jancaitis' applications, a region is first subdivided into four units. Sub- 
sequently, each unit in which the residuals at the reference points prove to be 
too large is itself subdivided. In the one example given by Schult, the inter- 
polation region is divided into eight units. De Masson d'Autume lists his D- 
matrices for no more than seven grid nodes per row and per column and he suggests 
to use only the nearest four or six in each row and column. The interpolation in 
each unit is then treated as an interpolation in the central interval of a 4x4 or 
a 6x6 grid and only one D-matrix is needed. This is a simplified solution in which 
the total surface is still continuous. However, it is no longer smooth at the 
boundaries. 
6. Interpolation in a rectangular grid 
Heights at the nodes of a rectangular grid may be obtained either by one 
of the preceding interpolation methods or by measuring along a set of equidistant 
parallel lines in a photogrammetric model. These heights may be augmented by other 
values such as slopes of the terrain surface. In each case, the density of the 
grid is made sufficiently great to allow interpolation in each grid element by 
means of a separate polynomial that is of the first-, second-, or third-degree with 
respect to each of the planimetric coordinates. The interpolation, therefore, will 
make use of some or all of the terms of the bicubic polynomial 
h = agg * ajox + agiy * a2ox? + ajay + ag? + agex’® 
* a»1z?y * ajomy? * aggy? * a31x%y + a22x2y? 
*ajgmy? * agom)y? * a23z?7y? * a3323y? 
In matrix notation, this equation can be written 
h* [1222 23]A[1 y y? y?]! 
or simply 
h = xTAy (6c) 
Here, x and y are the vectors used already in eqs. (5) and A is the matrix which 
has aij as the element in row Z+1 and column j+1. 
The parameters of each local polynomial must be computed from the given 
heights at the nodes of its grid element and from any other provided values at 
the boundaries of that element. It follows from eq. (6) that the parameters are 
linear functions of those values. By shifting the origin of the coordinate system