58 - 
In the case of the bicubic polynomials, the total surface is constrained 
to be smooth. The sum which expresses its flatness is here the sum of squares of 
the first derivatives of the surface with respect to the two planimetric coordinates 
z and y at all grid nodes. The unknowns solved from the equations are the height 
h, the two first derivatives h, and h,,, and the mixed second derivative h…, at each 
grid node. Therefore, the number of {Unknowns is four times the number of Grid 
nodes. The 16 parameters of each local polynomial are then computed from the 16 
. values at the four surrounding nodes. As will be shown later, this ensures con- 
tinuity and smoothness of the total surface. 
In the case of the linear polynomials, the sum which expresses the smooth- 
ness is the sum of squares of the slopes along the grid lines. The unknowns are 
here the heights at the grid nodes and, therefore, they are equal in number to the 
number of nodes. 
Junkins and Jancaitis [36] make use of local polynomials which similarly 
must satisfy continuity and smoothness conditions along their boundaries and which 
are Submitted to further constraints "which could be altitude or slope requirements". 
They mention as a typical local polynomial the 15-term fourth-degree polynomial. 
In an experiment, not the terrain height but its reciprocal is expressed as such 
a polynomial. There is one equation for each constraint and the weighted sum of 
the height residuals at the reference points is minimized, subject to the 
satisfaction of the constraints. The method of least squares is here used with 
Lagrangian multipliers. These multipliers must be computed together with the 
coefficients of all local polynomials. For details of the method, reference is 
made to the doctoral thesis of Junkins. 
Without going into details, Schult [64] mentions that such a method has 
been programmed also at the IAGB of Stuttgart University. Presumably, as in a 
program described earlier in his paper, bilinear local polynomials are used. Only 
continuity conditions are mentioned. 
De Masson d'Autume [16] presents an interpolation method which is derived 
from the approximation of a function of one variable by a spline function. The 
approximating function has the property that the integral of the square of its 
second derivative over the range of the variable is a minimum. It consists of a 
set of third-degree polynomials, one in each interval between reference points, and 
having the same first and second-degree derivatives at their joins. Interpolating 
along a grid line y = constant with n equidistant reference points, the interpo- 
lation formula for each interval between two reference points can be written 
h = [1 x x? x3]Dz (5a) 
Here, z is the vector of the n reference heights and D is a 4xa matrix. 
De Masson d'Autume has listed the matrix D for the cases of from two to seven 
reference points and in each case for each interval between two points. 
This function can be used to derive a local polynomial for each grid 
element as follows. The interpolated height of a point with coordinates x,y can 
be written as a function of the heights at the grid nodes by first interpolating 
with eq. (5a) on all grid lines y=constant for the required value of x and then 
Injerpiiating in the y-direction along the line x=constant. The result can be 
written: 
h » x D4HD]y (5b) 
Here, x is the column vector whose components are the four powers of x, y is the 
corresponding column vector for the y-coordinate, D, and D, are the D-matrices for 
the x- and the y-interpolation, respectively, and H is the matrix of grid heights, 
arranged in each column and in each row according to increasing x- and y-coordinates, 
respectively. If one interpolates first in y-direction and then in x-direction,