DIGITAL ELEVATION MODEL INTERPOLATION BY FUSION OF MORPHOLOGICAL RECONSTRUCTION AND DISTANCE TRANSFORMATION

Interpolation methods have significant impacts on the accuracy of the digital elevation model (DEM) from contours which are one of frequently employed data sources. In this paper, an interpolation method is presented to build DEM from contour lines by fusion/integration of morphological reconstruction and distance transformation with obstacles. Particularly, morphological reconstruction is used to get the elevation values of the higher contour lines and the lower contour lines of any a spatial point between two contour lines, and distance transformation with obstacles is used to get the geodesic distances of the spatial point to the higher contour lines and the lower contour lines respectively. At last, linear interpolation along water flow line is used to get the elevation values of the pixels to be interpolated. The experiment demonstrates that feasibility of our proposed method. * Corresponding author


INTRODUCTION
Digital elevation model (DEM) is a discrete digital expression of earth surface terrain (Li and Zhu, 2003).Since the late 1950s, DEM has drawn great attention and been widely used in many fields, such as mapping, civil engineering, geology, mine engineering, landscape architecture, road design, flood control, agriculture, plan, military engineering, aircraft and battlefield simulation, etc.Generally, different data sources need different interpolation methods for DEM generation.At present, the data of generating DEM mainly come from topographic map, remote sensing data (including image data of aerospace and synthetic aperture radar interferometry and LiDAR data), ground survey and the existing DEM, etc., among which the method from topographic map is the most commonly-used one (Tang et al., 2005).Survey and mapping department in China generates DEM with different resolutions by applying digital line graphics of 1:10 000, 1:50 000 and 1:250 000 scale.Usually, there are three ways to generate DEM from topographic map due to the distribution characteristics of contour lines (Li and Zhu, 2003): contour line discretization, contour line interpolation and contour line constructing Delaunay triangulated irregular network (TIN).The essence of contour line discretization is to regard the contour lines as irregularly distributed data without considering its own topographic characteristics (Li and Zhu, 2003), which causes anomaly to DEM generation.The interpolation of generating DEM by most fast line (water flow line) from contour lines is one method with simple principle.But there are many problems in realizing this method.Due to the problem of "flat triangle" (namely horizontal triangle) (Hu et al., 2007a) in building TIN directly by contour lines, now contour lines as well as the attached "feature data" (such as terrain structure lines and feature data points like surface peak, concave points and saddle points) are commonly used to construct TIN in practical engineering production.In recent years, many new interpolation methods have been proposed, and the research achievements by Hu et al. (2006) and Hu et al. (2007b) are of representativeness."Feature data" is the dual pattern of contour lines in essence, which is nonessential.Particularly, it is hard to control the density of feature data to balance the accuracy and workload of DEM during the production.Therefore, map algebra can be used in interpolating DEM by contour lines directly, which is MADEM.Map algebra is an image manipulation built upon the distance transformation operation.When using it to interpolate DEM, extra supplementary feature data is not necessary.Besides, DEM generated in this way has higher accuracy, which satisfies DEM accuracy evaluation criteria of "elevation order isomorphism" (Hu et al., 2009;Hu et al., 2011).However, the interpolation method based on map algebra is urgent to be improved.This method is based on solving half-interval contours (which are lines of equal distance to the adjoining contour lines) C l /2, C l /4, C l /8, C l /16, C l /32… (C 1 is the base contour line interval in topographic map) through iteration to generate DEM, which means iterating the edges of Voronoi diagram of the adjoining contour lines and assigning the mean value of two adjoining contour lines to the corresponding edge until it is not necessary to divide.Then bounded by 1/2n+1 Voronoi diagram (n is the largest iteration times), the elevation is valued for layers (Fan, 2007), which is a linear interpolation method in essence.But by this method, the largest number of iteration times has to be input in advance.Once the actual time exceeds the input number, the points failed to be interpolated are all assigned to the increment of last half-interval contour.It is clear that the largest number of iteration time is a critical index of interpolating DEM.However, this index adopts manual intervention without self-adaption, which limits the engineering application to a certain extent.In this paper, according to the conception of geodesic distance in mathematical morphology, a contour line interpolation method is proposed based on morphological reconstruction (Lantuejoul and Maisonneuve, 1984) and distance transformation with obstacles (Hu et al., 2007c).This method needs no manual intervention.

FUNDAMENTAL PRINCIPLES
The interpolation method proposed in this paper refers to two classical methods, including geodesic transformation (Lantuejoul and Maisonneuve, 1981) and geodesic measurements (Soille, 2008).

Geodesic transformation
The fundamental operation of geodesic transformation includes two operators, geodesic dilation and geodesic erosion (Lantuejoul and Maisonneuve, 1981;Soille, 2008).Geodesic dilation revolves mark image and mask image.The size and domain of definition of these two images are the same.But every pixel value of mask image has to be no less than the corresponding pixel value of mark image.The process of realizing geodesic dilation is dilation operation for the mark image by using basic isotropic structural elements.During the process, the resulting image has to be maintained under the mask image, which means the mask image restrains the dilation and extension of mark image.Similarly, geodesic erosion demands the mask image to be no greater than the mark image.The process of realizing geodesic erosion is erosion operation for the mark image by applying basic isotropic structural elements.During the process, the resulting image has to be maintained above the mask image, which means the mask image limits the shrink of mark image.The transformation of geodesic dilation or geodesic erosion of bounded images will converge after a certain cycles until the expansion or shrink of mark image is stopped by mask image thoroughly.After that, the cycle starts again.The value of any a pixel of mark image does not change this time, which is the morphological reconstruction principle of mask image from mark image. .In this paper, morphological reconstruction is used to generate two plateau images to identify the elevation value of the contour lines (the upper contour line and the lower contour line) adjoining any a spatial point.

Geodesic distance
Compared to the Euclidean distance, geodesic distance of mathematical morphology refers to the shortest path linking two pixels of a subset in the image, in which the subset area is called geodesic mask (Lantuejoul and Maisonneuve, 1981;Soille, 2008).Moreover, geodesic distance and geodesic path between two independent pixels in geodesic mask are mainly determined by mask shape.If the shape of mask is convex, then geodesic distance is equivalent to Euclidean distance and geodesic path is a straight line.Otherwise, geodesic distance is affected by pixel location and geodesic mask shape.The calculation of geodesic distance in this paper is realized by Euclidean distance transformation with obstacle space (Hu et al., 2007c).

Euclidean distance transformation with obstacles
Distance transformation is the transformation or process of calculating and marking the distance between each point of spatial point set and reference (Hu et al., 2007c), which is divided into Euclidean distance transformation and taxi distance transformation.The former marks the Euclidean distance, which is appropriate for natural-shaped graphic.The latter marks Manhattan distance, which is appropriate for regular-shaped graphic.Distance transformation in this paper refers in particular to Euclidean distance transformation.Moreover, the detailed technological process of Euclidean distance transformation with obstacles can be referred to Hu et al. (2007c).In this paper, Euclidean distance transformation with obstacles is applied to identify geodesic distance of any a spatial point to the adjoining upper contour line and lower contour line.

THE PROPOSED METHOD
Based on the above principles, the following interpolation method of DEM generation by contour lines is designed in this paper, and it is composed of three main parts.1) Acquirement of vector contour lines If the original topographic map is papery, the contour lines data in the form of vector can be obtained through digitalization.If the supplied original data are in the form of digital line graph, its contour lines can be extracted.Here are the following requirements of vector contour lines: ① one contour line within a map sheet is not allowed to break off without reason, otherwise the complete data of contour line is obtained through join operation; ② each contour line has to be assigned with corresponding elevation value.The elevation values of all contour lines have to be sorted.The type partition of odd or even contour lines starts from the minimum value.In principle, the adjoining elevation values have to be an odd and an even alternately.
2) Generation of two elevation images First, vector contour line is rasterized.In the raster image, each pixel of contour line is filled with the corresponding elevation value.For the pixel out of any contour line is filled with an invalid value.As shown in Fig. 1, the image is noted as contour line image CL.In Fig. 1, Pixel point p in the white area is surrounded by two contour lines.The contour line with greater elevation value is expressed as Cu(p).The contour line with lower elevation value is expressed as Cl(p).Meanwhile, p can be surrounded only by Cu(p) or Cl(p) when p is surface peak or concave point.At this point, Cu(p) = Cl(p).

EXPERIMENTS AND ANALYSIS
Experimental data adopts one 1:10 000 contour line data, and the contour line interval is 10.0 m.The contour line elevation of the scene distributes between 630.0-920.0m.The length and width are 1400.0m × 1100.0 m.DEM, whose resolution is 1 m, generated by the method in this paper is showed in Fig. 2(f).In addition, as is showed in Fig. 2

Figure 1 .
Figure 1.The upper contour line and the lower contour line associated with some a pixel p As for each pixel point, if the following two elevation images Pl (lower plateau image) and Pu (upper plateau image) can be (g), the contour line of 10.0 m interval generated by DEM in Fig. 2(g) can approximate the original contour line precisely, which proves the correctness of our method.(a) The contour lines (b) The lower plateau image superimposed by the contour lines (c) The upper plateau image superimposed by the contour lines (d) The geodesic distance from the odd contourr liens to the even contour liens (e) The geodesic distance from the even contourr liens to the odd contour liens (f) The interpolated DEM by our method (g) The generated contour lines (yellow) from the DEM in (f) and the original contour lines (black) Figure 2 DEM interpolation from contour lines by integration of morphological reconstruction and Euclidean distance transformation with obstacles 5. CONCLUSIONS Interpolation method is an important factor on influencing the accuracy of digital elevation model.In the process of obtaining DEM from contour line, the selection of interpolation method is an important technological part, which has a significant effect on the accuracy of DEM generation.In this paper, an interpolation method of DEM generation based on the integration of morphological reconstruction and Euclidean distance transformation with obstacles is proposed.Morphological reconstruction is used to get elevation values of the corresponding upper and lower contour lines of the pixel.Euclidean distance transformation with obstacles is used to get geodesic distance of the point associated with the upper and lower contour lines.Linear interpolation is used to get the elevation value of the point.The experiment suggests that the proposed method is feasible.