Comparison of interpolation methods and image correlation in the collection of data for digital elevation model S

Original scientific paper Digital Relief Model (DRM) is the name used officially in Croatia, and we distinguish between the original and derived DRM. Referring to the regular raster structure of elevations stored in the derived DRM, the term digital elevation model (DEM) was used for this model. Three geomorphometrically different landscapes have been tested. The elevation points were created using two methods, interpolation from the data obtained in the photogrammetric survey and image correlation, in two resolutions (5×5 and 25×25 meters). The paper presents the differences in the application of two methods of calculating DEM. The causes of differences discussed and the relative estimation of the accuracy of each method are elaborated by comparing the height of points on new models with the heights on original DRM. The recommendation for the usage of the acceptable method is given related to the accuracy, the speed of data collection and a higher degree of automation.


Introduction
The definition of DRM as a set of positionally (and vertically) determined points and geometric elements (break lines, form lines and exclusion zones) that present the surface of the Earth and provide the calculation of the mathematical model of that surface [1] includes both the original and derived DRM.Due to the regular raster structure of elevations stored in the derived DRM, and for the purpose of the simpler distinction between the original and the derived DRM, the term digital relief model is used in the article for the original DRM, and digital elevation model (DEM) for the derived DRM.
The digital relief model is the basis for the spatial analyses.The parameters related to the creation of a digital model (resolution, generalisation degree,...) are dictated by various applications and user needs (urbanism, geomorphometric analyses, hydrological analyses,...).However, the mutual need of all users is the samedigital model that presents the true relief most reliably within the given limits.
A few authors have been dealing with the influence of the selection of optimal DRM resolution related to the geomorphometric characteristics of the true terrain.The samples having representative morphometric characteristics of the test area that we want to present using DRM are selected in the field.These samples are observed with much larger spatial density (resolution) than it is the expected resolution of DRM, and then the discretization density is determined by thus obtained small representative samples using the spectral analysis [2].In order to avoid the application of extremely complicated procedure of determining the density of data capturing, it is possible to make the analyses by the method of acquiring the characteristic profiles from the terrain and to obtain optimal result with allowed elevation deviation due to discretization of the continuous variable [3].This method is far more simple than the previously mentioned method of spectral analysis and convenient for practical use.
The creation of the digital elevation model with the resolution of 25×25 meters from the data of digital relief model (original DRM) using the method of interpolation has been done for the State Geodetic Administration (SGA).The so far made research has shown that the application of the method of image correlation (i.e. the application of automatic correlation needed in the production of digital elevation model) with the resolution of 5×5 and 25×25 meters has not been researched in relation to the results obtained by means of interpolation.Since the image correlation automates the manual procedures of data collection, the question of the optimal method of data collection for geomorphometrically characteristic areas of the Republic of Croatia has been raised.The paper presents the comparison of DEM obtained using the method of interpolation from original DRM with the DEM obtained using image correlation from aerial images and aerotriangulation data.

Research area and data
The comparison was made in three geomorphometrically different areas of the Republic of Croatia.The first area is a part of the very steep mountain range Biokovo (characterized by sinkholes and karren).The second area is a part of gently to significantly sloping terrain of the mountain Papuk (covered mostly by forests).The third area is the lowland of Eastern Croatia Technical Gazette 24, Suppl.2(2017), 419-425 in the municipality Nijemci (cut by drainage and irrigation channels).
The original DRM25 (Fig. 1), as well as the data for the reconstruction of the project of the aerial survey with included aerial photographs, the parameters of external orientation and the data about cameras, were used as input information.Digital relief model (DRM25) obtained from the Central Office of the State Geodetic Administration is the collection of individual characteristic points, rasters of height points, break lines and form lines needed for the height presentation of the Earth's surface, and a part of topographic and cartographic database.It is produced for TK25 in the reference scale of 1:25000 according to the data model of CROTIS [4].
The aerial photographs were made with analogous aerial cameras at the scale of 1:20000.They were scanned with the resolution of 21 µm in TIFF format, and the size of individual photographs (format 23×23 cm) is 11200×11200 pixels.
Since the entire area of the Republic of Croatia was not surveyed at the same time (various surveying projects), the data for various areas have different characteristics.In the area of Biokovo, the Cyclic aerial photogrammetric survey "Dubrovnik and the Surrounding area" was applied that was made in 2005 with analogous camera ZEISS RMK TOP 15.Those were colour photographs, and 16 of them were used for the analysis.For the area of Papuk, the Cyclic aerial photogrammetric survey "Western Slavonia" was applied that was made in 2001 with analogous camera WILD RC20 15/4.Those were black and white photographs, and 14 of them were used.For the area of Nijemci, the Cyclic aerial photogrammetric survey "Eastern Slavonia" was applied that was made in 2002 with analogous camera ZEISS RMK TOP 15.Those were black and white photographs, and 11 of them were used.
Two sets of data were made for each area to be used in the research.From the original DRM data obtained from SGA, the DEM with the resolution of 5x5 meters and DEM with the resolution of 25×25 meters (INT data) were made by means of interpolation.DEM with the resolution of 5×5 meters and DEM with the resolution 25×25 meters (CORR data) were made by means of image correlation from the aerial photographs.

Data capturing and processing
The evenly distributed elevation points (raster) were captured by INPHO digital photogrammetric system.The raster of elevation points was captured by means of two different methods: interpolation and image correlation.The interpolation method is the procedure that detects the function value between two known values [1].The image correlation is the procedure of determining the position of some detail on the image on the basis of its correspondence with the sample that was taken from another image [5].DRM25 was used for the purpose of interpolation, and aerial photographs and the data of aerotriangulation were used for capturing the raster of elevation points.
Both the method of interpolation and image correlation were made using the software MATCH-T.It combines linear method, the method of least squares and the Finite Element Method for the interpolation procedure.For the procedure of image correlation, it combines the area based data and the method based on texture properties.
The basic idea of the linear method of interpolation relies on the spline concept, using the first-degree polynomials instead of one high degree polynomial [6].The interpolation with the method of least squares (Least Squares Interpolation) approximates the default function by another of a certain type, so in a certain sense their mutual distance is as small as possible, regardless of what functions may not coincide at any point [7].The interpolation by means of the Finite Element Method FEM is used for the approximation per parts [8].
In the correlation of the area-based data, the entities for the comparison are tone values where it is not sufficient to compare only one pixel, but the classification of tone values in individual parts of the image should be compared with the sample taken from the other image or computer base.The feature based matching is based on the comparison of the object features like points, edges or areas presuming that the relative orientation has been done.The method shows the best results in the areas where the texture is clearly expressed.The procedure itself is based on the calculation of the parameters of interest (depending on grey tone values) within the analysed part of the digital image (window) for each window on the image, and on the comparison of the values with minimal, i.e. maximal values given in advance.As the points of interest, only those windows are accepted having the parameter values larger or smaller than the minimal, i.e. maximal values [9].
Two data sets were created for each area for the purpose of the research.DEM with the resolution 5×5 and 25×25 meters (INT data) was produced from the original DRM data by means of the interpolation, and the method of image correlation from aerial photographs and aerotriangulation data was used to produce DEM with the resolution 5×5 and 25×25 meters (CORR data).After being captured, the data were checked and corrected in DTMaster module.For the purpose of analysing the differences between DEMs obtained by the interpolation and correlation methods, the gathered data were transformed from the vector into raster form.Due to different methods of capturing the points for the same areas, smaller differences appeared in the coverage with DEM, and all presentations for individual area were reduced to the same coordinates and the same number of pixels.DEM5 was reduced to 614 400 (960 × 640) pixels, and DEM25 to 24 576 (192 × 128) pixels.
To notice the differences more clearly, the parts of the areas Biokovo, Papuk and Nijemci were separated covering 16577 (137 × 121) pixels for DEM5 and 621 (27 × 23) pixels for DEM25.The data obtained by means of the interpolation method and the method of image correlation were presented for 5 and 25 meter grid (Figure 2, 3 and 4) with matched height ranges in legends.

Analysis of the data obtained using interpolation and correlation methods
The noticing of differences in the application of both methods of creating DEM (interpolation/correlation) is an important factor in making a decision about the correctness of using individual methods on certain terrain.The elevation differences between the model created using the interpolation method and the model created by means of the method of image correlation vary depending on the test area (Figs. 5, 6 and 7).For the purpose of the analysis using DEM created by means of the interpolation method and DEM created by means of the method of correlation, only those pixels were separated that are arranged on the same coordinates like the points at which the heights were measured in the original DRM.The problem arising in the separation of the points on the same coordinates was caused by the resolution of the digital elevation model and by defining the position of the pixel starting point.Not all software programmes have the pixel starting point at the same place.Hence, the export of raster data into vector form causes a movement related to the DRM points created in the vector software.The problem was solved by converting the point of the original DRM into raster form using the software IDRISI with the movements at coordinates (within the pixel sizes) appearing in the process.
Hence, the estimation of accuracy was made to define the acceptable method (Tab.1).In the case of significantly sloping terrain with a large number of sinkholes and cracks (Biokovo) and of gently sloping terrain covered with forests (Papuk), break lines and form lines (Fig. 1) were added to the raster of elevation points, so the interpolation method has yielded better results.Relatively high value of standard deviation is the consequence of raster resolution and of the error in converting the vector into raster.Smaller differences in the values of standard deviations for the area of Nijemci are the consequence of the plain relief where both methods yield equal results.The error resulting from the raster resolution and the conversion of the vector into raster is caused by the inability to define the position of each point within the raster accurately.Two or more points of the original DRM with various elevations fall within one raster pixel that can be defined by only one value (Fig. 8).The pixels with the value 0 are the pixels where no point of the original DRM is located.The programme calculated their value by means of interpolation from the neighbouring pixels.
The IDRISI module calculated the pixel elevation from the points within a pixel for the purpose of converting vectors into raster.Hence, the deviations at individual points of the original DRM are larger, which finally results in the standard deviation larger than the allowed one.The method of Inverse Distance Weighing (IDW) is used for the determination of pixel elevation out of more points within IDRISI [11].The values at the Technical Gazette 24, Suppl.2(2017), 419-425 selected points are estimated on the basis of the following formula , 1 where: Z ˆ is the interpolated value of elevation; d the distance to the known point i; Z the value of the elevation of the known point i; p the distance weighting exponent; n the number of points included into the search, and i is the number of known points included into the calculation.The influence of each point is inversely proportional to its distance from the position where the value is estimated.The number of points included into the estimation is determined with the radius of the circle set around the mentioned position.The result of the method depends greatly on the value of distance exponent that is most often set to 2 because the map form, i.e. the rounding of lines is the simplest in this way [10].The larger the distance exponent, the larger is the influence on the closest points.
Due to the raster resolution, the estimation of accuracy is not absolute accuracy of an individual model.However, the method of converting the vector data into the raster system of IDRISI is the same for both methods of data creation, and the accuracy estimation obtained by calculating the standards deviation of elevations on the pixels created from the set of points yields unanimous and comparable results.
The terrain of the Biokovo area is rather rugged, and large elevation differences on a small area create various perspectives on the images of one stereo model that makes the identification of identical points uncertain, while in the area of Papuk, the surface covered with forests prevents the correlation on the ground.In the case of uncertain detail identification, the operator can move the cursor manually to the identical detail on the model for the area of Biokovo because of the terrain configuration in order to improve the accuracy of the obtained data, which would improve the accuracy of data obtained by means of the method of image correlation, but a large number of sinkholes and cracks would make this procedure rather difficult and time consuming especially for DEM5.
On the forested hilly terrain in the area of Papuk, the image correlation is a rather demanding task.Due to dense forests, the terrain is invisible at many places, which makes the procedures of stereoscopic point measurement in the field rather demanding and long lasting.The operator must check every correlation made in the forest area and "bring down" the selected point or the area to the true terrain elevation.
The lowlands of the Municipality Nijemci show small deviations between the methods of interpolation and correlation on the created models with both methods yielding satisfactory standard deviation.The participation of the operator who would correct the inaccurate identifications is a demanding procedure on the plain terrain without forests.
The automatic correlation with additional stereoscopic vectorising of channel network would result in better accuracy for the area of the Municipality Nijemci although the values of standard deviation are equal for both methods.Apart from that, special attention should be paid to the selection of adequate interpolation function referring to morphometric characteristics of the terrain.In this way, the appearance of unreal oscillations would be avoided in digital terrain model that are not present in reality (Fig. 3).
Tab. 2 shows the recommendations for the selection of methods in individual areas needed for the creation of DEM5 and DEM25 from DTM25 and from photogrammetric images and aerotriangulation data.

Conclusion
The method of image correlation is an automated procedure in which the programme identifies the identical details on images, which largely accelerates the procedure of capturing the rasters of elevation points.Although the presented results are obtained by means of the method of image correlation without the participation of an operator, the possibility of including the operator into the procedure of correlation is allowed in the programme settings for image correlation.In this case, it is possible for the operator to correct every "uncertain" correlation by setting the cursor onto the same detail on both images.
The values of standard deviation for the area of Biokovo are high for both resolutions and both methods of data capturing.The area of Biokovo is very rugged with a large number of sinkholes and karren, and it can be presumed that the programme could hardly identify the identical details by means of the image correlation procedure either due to the similarity of details or due to bad perspective of the stereo model.Although the DRM made for DEM25 was used for the interpolation, the values of the standard deviation are significantly worse for the resolution 25×25.The geomorphometric forms of the researched area vary in their size and are often significantly smaller than the defined resolution, which resulted in higher values of standard deviation.Completely automated and unselective data capturing (CORR data) for DRM in the areas of extremely rugged geomorphometry cannot yield reliable results, but the data obtained in such a way should be supplemented by other elements of relief modelling (structural lines and characteristic elevations).INT-data contain exactly such elements and yield, therefore, more reliable and more accurate results in all situations.
The values of standard deviation for the area of Papuk are smaller when using the interpolation method than the values when using the method of image correlation.In the area of Papuk, the interpolation method yields higher accuracy due to the terrain configuration and high-quality input data (the areas that were invisible because of forests were supplemented with the contour lines from the map TK25).The area of Papuk is covered with forests that make the automatic correlation difficult in the field where the terrain is invisible due to the density of forests.The procedure of bringing down the elevation to the terrain is rather demanding and long-lasting because the operator must check each correlation made in the forest area (most of the area), and "bring down" each individual point or selected area to the real elevation, which makes the interpolation method more acceptable.
The plain area of the Municipality Nijemci shows on the created models DEM5 and DEM25 small deviation between interpolation method and correlation method with both methods resulting in satisfactory standard deviation, and considering thereby the speed of data capturing, the image correlation has been evaluated as more acceptable.

Figure 2 Figure 3 Figure 4
Figure 2 DEM5 (above) and DEM25 (down), created by means of interpolation method (a) and the method of correlation (b), Biokovo a) b)

Figure 5 Figure 6 Figure 7 Figure 8
Figure 5 Differences in the application of interpolation methods (a) and the method of correlation (b) in meters, Biokovo, DEM5 (a) and DEM25 (b) a) b)

Table 1
Standard deviations for all models in meters

Table 2
Recommendations for the selection of methods for creating DEMs