«COMMISSION IV, WORKING GROUP 1 PRESENTED PAPER HEINRICH EBNER, BERNHARD HOFMANN-WELLENHOF, PETER REISS, FRANZ STEIDLER LEHRSTUHL FOR PHOTOGRAMMETRIE ...»
14. CONGRESS OF THE INTERNATIONAL SOCIETY FOR PHOTOGRAMMETRY
COMMISSION IV, WORKING GROUP 1
HEINRICH EBNER, BERNHARD HOFMANN-WELLENHOF,
PETER REISS, FRANZ STEIDLER
LEHRSTUHL FOR PHOTOGRAMMETRIE
TECHNISCHE UNIVERSIT~T MONCHENARCISSTR. 21 POSTFACH 202420 0- 8000 MONCHEN 2
HIFI - A ~1INICOMPUTER PROGRAM PACKAGE FOR
HEIGHT INTERPOLATION BY FINITE ELEMENTS *)
ABSTRACTA general minicomputer program package for height interpola- tion is presented. It is written in FORTRAN and interpolates a digital he i ght model (OHM) from arbitrarily distributed re - ference points and breaklines, using the Finite Element Me- thod. From the OHM digital profiles can be derived for ortho- photo production and digital contour lines can be interpola- ted and plotted. The paper describes the principles of inter- polation and the structure of the programs. Finally the effi - ciency of the package is demonstrated by the results of pro- cessed practical examples.
*)A German version of the paper is published in Ze i tschrift fUr Vermessungswesen, May 1980 202.
1. INTROOUDUCTION AND REV! EWThe paper presents a general FORTRAN program packag~ for height interpolation. It is developed for the use w1th the minicomputer Hewlett-Packard HP 1000, but it can also be trans- ferred to other minicomputers and to large computers.
The program package interpolates the heights of grid points from arbitrarily distributed reference points and points along break lines. The meshes of the grid are squares, the size of which can be chosen by the user. From the digital height model (OHM), formed by the grid points, digital contour lines and digital profiles can be derived optionally.
For interpolation of the OHM the Finite Element ·Method is used. This general method is successfully applied in numeri- cal mathematics and several fields of engineering over more than two decades /1/, /2/. During the last few years the me- thod has also been applied in geodesy /3/, /4/. Similar to these applications the OHM interpolation is based on rela- tively simple and regular structures. Local surfaces are de- fined for all meshes of the grid of interpolation points.
These FINITE ELEMENTS are tied together accordingly and form the total interpolation surface. The HEIGHT INTERPOLATION BY FINITE ELEMENTS -leading to the name-HIFI or the program package - results in an interpolation surface of minimum cur- vature, which approximates the given reference points with optional filtering and represents the break lines adequately.
Two types of finite surface elements are used for the indivi- dual grid meshes, bilinear or bicubic ones. The bilinear va- riant was presented already in 1978 at the International Geodetic Week in Obergurgl /5/ and at the ISP Commission III Symposium in Moscow /6/. If the bilinear elements are replaced by bicubic elements the mesh width of the grid can be extended and the computational effort decreases. On the other hand, consideration of break lines is much more complicated with bicubic elements. The program package HIFI therefore applies the bilinear variant for terrain with break lines and the bicubic variant if no break lines are to be considered.
The development of the program package is supported by Carl Zeiss. HIFI can be used in combination with the analytical plotter PLANICOMP C 100 and the new analytical orthoprojector ORTHOCOMP Z2 /7/.
Two program versions are available. The smaller one HIFI-P interpolates a OHM from terrain data without break 1 ines and uses bicubic finite elements. Reference points may be arranged along contour lines, along profiles, in a grid or may even be arbitrarily distributed. From the OHM, interpolated by HIFI-P digital PROFILES can be derived and used for orthophoto pro~ duction with the ORTHOCOMP Z2.
The extended program version HIFI-PC interpolates a OHM from arbitrarily distributed reference points and break lines, using bilinear finite elements. If no break lines are given, automatically the bicubic variant is applied. Further on 203.
HIFI-PC offers the options to derive digital PROFILES (as above) and d1gital CONTOUR LINES from the OHM. The contour interval can be chosen by the user and the computed contour lines can be plotted by a digital drawing table of Zeiss.
Both program versions supply absolute orientation of the model data onto given control points. Besides profiles and contour lines further informations can be derived from the OHM princi pally. An according example for production of slope maps is given in /8/.
A more detailed description of the principles of interpolat ion of OHM, profiles and contour lines is given in the following chapter. Then the program versions HIFI - P and HIFI-PC are presented in the form of block diagrams. Finally the efficiency of the program package is demonstrated by the results of processed practical examples.
2. PRINCIPLES OF OHM INTERPOLATION AND DERIVATION OF DIGITAL
PROFILES AND DIGITAL CONTOUR LINES
2.1 OHM interpolation by bilinear finite elements In this case for each grid mesh a separate bilinear polynomial is used and the interpolation surface, formed by these finite elements is assumed to be continuous. The mathematical formu lation is founded on the use of bilinear base sp·finc:s Sij which are defined for all grid points Pij · The base of these splines is locally bounded to the four surrounding grid meshes with a maximum value of 1 at the grid point Pi j· For a more detailed respresentati.on see for instance /J/. The interpolation surface then is the result of a superposition of local surfaces a.; j · S; j· This means, that each base splineS; j makes a contribution to the interpolation surface. In the case of bilinear finite elements the values a.; j are identical with the functional values or heights of the'1nterpolation surface at the grid points Pi,j.
If a point Pk of coordinates xk, Yk is located inside a mesh, as shown in figure 1, the corresponding height hk of the in terpolation surface follows as a linear function of the heights of the four surrounding grid points P; j• P;~ 1 j' P; j+1• Pi+1 j+l and the respective coefficients are' ' bil;near functions of the local coordinates 6xk, 6Yk· P.. p i+l, j-+2 pi-1,)+2 I, J+2
The requested grid heights of the OHM are defined as unknowns and estimated from the heights hk of the available reference points and a general curvature minimization of the interpo lation surface with filtering at the reference points. This problem can be interpreted as an adjustment of indirect observa tions and be solved by 1 east squares.
The height hk of the respective reference point Pk (xk,Ybhk) is treated as an observation of the unknown height of the interpolation surface at position Xk,Yk · This unknown height however, is a linear function of the four surrounding grid heights, as shown above. So, reference point Pk leads to an observation equation, containing four unknowns, the coeffici~nts of which are linear functions of 6xk and 6Yk· Equation (1) represents this relation in detail, vk is the least squares correction of hk.
vk = (1-6xk) (1-6yk)hi,j + 6xk(1-6yk)h. + 1 'J.
1 + (1 - 6xk)6yk h.. 1 + 6 Xk 6 Yk hi+1,j+1 - hk ( 1 ) 'J + 1 All equations of type (1) can be weighted individually, ac cording to the accuracy of the observations hk and to the desired amount of filtering.
With the bilinear variant curvature minimization can not be realized rigorously, because all second derivatives o.f the interpolation surface are zero. Therefore second differences of the heiqhts are used instead. For each grid point such second differences are formulated in x- and y-direction and are interpreted as fictitious observations of amount zero and certain stochastic properties. The corresponding observation
equations for grid point Pi,j• shown in figure 1, then read :
v.. = h.. - 2 h.. +h.. 0 1+ 1,J 1- 1,J (2) 1,1,J 1,J o v.. = h.. - 2 h.. +h..
1 1,J - 1 2, 1,J 1,J 1,J+ T he m s t s i mp l e way i s, to t rea t a 1 1 f i c t i t i o u s o b.s e r vat i o n s a s o uncorrelated and of equal accuracy, which is adequate for homogeneous and isotrop terrain. The suggested concept howeve~ is not restricted to this assumption, but allows for more general stochastic properties of the fictitious observations, by which the local shape of the terrain surface can be modelled too. A relevant example for one dimensional interpolation is given in /6/. The general realization of the concept is still in work.
Break lines, given as a sequence of single points in space, are considered in a direct and simple way, as shown in fig. 2.
For this purpose the respective breakline firstly is represented as a polygon, eventually as a filtered one. Then the intersections of the break line with the vertical planes, passing through the grid lines of the OHM, are computed and the observation equations (2) are altered in such a way, that connections across the break line are excluded and sharp bends may appear at the intersection points.
205. Fig. 2
The equations of type (1), put up for all reference points and the total of equations (2), altered accordingly in the neighbourhood of break lines, form a system of observation equations, from which, in case of a diagonal weight matrix, a banded structured system of normal equations is obtained.
Assuming a rectangular area of n grid lines a:-rrd m grid points per line, the number of unknowns is m·n and the band width is proportional to m. A direct solution of the normal equations, as used in the program package HIFI, then leads to a total computing time of const.m3n, whilst the corresponding computing time per grid point is proportional to m2. From there it follows, that the computational effort can be reduced con siderably by a limitation of the number m of grid points per line. This is possible by subdividing the total area into se veral strips with full number n of lines, but less grid points per line. These strips have to overlap adequately and are treated successively.
In connection with the presented OHM interpolation by bilinear finite elements two related publications shall be mentioned /9/, /10/. The concepts of interpolat i on suggested in the~e articles, show close similarities to the method used here.
Both publications and the present one are independent of each other and confirm the topicality of the suggested method and its efficiency.
2.2 OHM INTERPOLATION BY BICUBIC FINITE ELEMENTSIn this case for each grid mesh a bicubic polynomial is used and the interpolation surface, formed by these more general finite elements, is assumed to be continuous. Continuity here however, is not restricted to the functional values, but ex tended also to the first and second derivatives. The resul ting smoothness of the interpolation surface allows for an ex tension of the mesh width of the grid in comparison to the bi linear variant. This is shown in figure 3.
Fig. 3 Interpolation surfaces built of bilinear (left) and bibubic (right) finite elements Similar as in 2.1 the mathematical formulation is founded on the use of base splines Si j. which are defined for all grid points Pi j. In the present case however, bicubic base splines are used. 'The base of these splines is extended to the 16 surrounding grid meshes, with a maximum value at the grid point Pi j· The interpolation surface again is the result of a superposition of local surfaces a; j'Si,j, but the individual values a; j. specifying the contr1 butions of the base splines Si j to the interpolation surface, are not identical with the he1ghts hi,j at the grid points.
With the assumptions of figure 1 the height hk of the interpolation surface at·position Xk,Yk then follows as a linear function of the parameters a of the 16 surrounding grid points.
The respective coefficients are bicubic functions of the local coordinates 6xk, 6Yk· In the adjustment the parameters a are defined as unknowns and estimated from observation equations for the reference points and for the curvatures at the grid points. These equations are much more complicated than equations (1) and (2) and can not be represented here in detail. It shall be stated however, that the observation equations for the reference points typi cally contain 16 unknown prameters a, whilst the curvature equations contain only 9. The curvatures themselves result from the second derivatives of the interpolation surface. As mentioned previously no consideration of break lines takes place with the bicubic variant.
The informations from 2. 1, concerning the structure of the system of normal equations, the band width and the subdivision of the total area into overlapping strips are also valid here.
Result of the adjustment are the parameters a, by which the interpolation surface is completely determined. The heights at grid points or at arbitrary positions xk, Yk can directly be computed from the parameters a of the surrounding 9 or 16 grid points. This allows for a rigorous densification of the OHM in a simple way.
A OHM interpolation, based on piecewise bicubic polynomials was suggested by Kubik already in 1971 /11/. The bicubic variant of the program package HIFI is similar to Kubik 1 s concept.
The chosen formulation by base splines however, 1 eads to a considerable reduction of the computational effort.
2.3 Derivation of digital profiles from the OHM From the OHM, interpolated according to 2.2, digital profiles can e as i l y be d e r i v ed and us e d for control l i n g an or tho pro j e c tor. To exhaust the height information of the bicubic OHM, it is recommended to choose the distance between adjacent profiles and between points along the profiles smaller than the mesh width of the OHM. The heights of the individual profile points are computed from the parameters a, as described in 2.2.