LINEAR AND NONLINEAR APPROACH FOR DEM SMOOTHENING

One of the biggest problems faced while analyzing digital elevation models (DEMs), particularly DEMs that are produced using photogrammetry, is to avoid pits and peaks in DEMs. Peaks and pits, which are errors, are generated during the surface generation process. DEM smoothening is an important preprocessing step meant for removing these errors. This paper discusses two linear DEM smoothening methods, Gaussian blurring and mean smoothening, and two nonlinear DEM smoothening methods, morphological smoothening and morphological smoothening by reconstruction. The four methods are implemented on a photogrammetrically generated DEM. The drainage network of the resultant DEM is obtained using skeletonization by morphological thinning, and the fractal dimension of the extracted network is computed using the box dimension method. The fractal dimensions are then compared to study the effects of the four smoothening methods. The advantages of nonlinear DEM smoothening over linear DEM smoothening are discussed. This study is useful in landscape descriptions.


Introduction
A digital elevation model (DEM) is a set of points defined in a three-dimensional Cartesian space (X,Y ,Z) that approximates a real surface.The Xand Y -axes may be expressed as geographic coordinates (i.e., longitude and latitude), whereas the Z-axis usually represents the altitude above sea level [3].It is a digital file consisting of the terrain elevations for ground positions at regularly spaced horizontal intervals [12].The data may be presented in ASCII or 2-byte integer binary formats.The most commonly used data structure for DEMs is the regular square grid in which available elevations are equally spaced in two orthogonal directions [15].
DEMs are a popular source for digital terrain modeling and watershed characterization because of their simple data structure; widespread availability [12] and they lend themselves to many GIS processes and operations [15].The most important use of DEMs is in terrain characterization.Terrain characterization refers to the determination of attributes of terrain, such as elevation at any point, slope and aspect and finding features on the terrain, such as drainage basins and watersheds, drainage networks and channels, peaks and pits [12].The nongraphic applications of DEMs include modeling terrain and gravity data for use in the search for energy resources, calculating the volume of proposed reservoirs, and determining landslide probability [9].A comprehensive review of the scientific, commercial, industrial, operational, and military applications of DEMs is proposed in [11].
DEMs can be derived indirectly from digitized topographic maps or directly through photogrammetric processing of medium scale, black and white, metric arrays photographs [9].The generation process of DEMs from digitized topographic maps is discussed in [1], while the generation process of DEMs through photogrammetric processing is discussed in [2].
One of the biggest problems faced while analyzing DEMs, particularly DEMs that are produced using photogrammetry, is a number of pits and peaks in DEMs are errors caused by poor correlation between images and buildings and trees on the terrain surface in the DEM surface generation process [8].These errors are also generated during the surface generation process [1].The process of removing these errors, DEM smoothening, is an important preprocessing step in most DEM analysis algorithms.
In this paper, two linear DEM smoothening methods, Gaussian blurring and mean smoothening and two nonlinear DEM smoothening methods, morphological smoothening and morphological smoothening by reconstruction.A linear function is a function f , which satisfies: for all x and y in the domain and scalars a [16].Nonlinear functions are functions which do not satisfy the conditions in (1.1).
In Section 2 of this paper, mathematical morphology and the relevant morphological operators are introduced.In Section 3, the two linear DEM smoothening methods, Gaussian blurring and mean smoothening are explained.In Section 4, the two nonlinear DEM smoothening methods are discussed.In Section 5, the results obtained from the described procedures are discussed.

Mathematical morphology
Mathematical morphology deals with extracting image components that are useful in representation and description of region shape, such as boundaries, skeletons, and the convex hull [7].Mathematical morphology is well suited to the processing of elevation data because in morphology, any image is viewed as a topographic surface, the gray level of a pixel standing for its elevation [15].Morphological operators generally require two pieces of data as input.The first input is the input image (referred to as M), which may be either binary or grayscale for most of the operators.The other input is the kernel (referred to as S).The kernel is used to determine the precise details of the effect of the operator on the image [6].The two basic morphological operators are erosion ((M S) ⊂ M) and dilation ((M ⊕ S) ⊃ M).The basic effect of erosion on a binary image is to erode away the boundaries of regions of foreground pixels (resulting in the image being shrunk).The basic effect of dilation on a binary image is to gradually enlarge the boundaries of regions of foreground pixels (resulting in the image being expanded).Erosion and dilation are used to form opening and closing [5].Opening (M • S = (M S) ⊕ S) is done by performing erosion followed by dilation.The effect of the operator is to preserve foreground regions that have a similar shape to this kernel, or that can completely contain the kernel, while eliminating all other regions of foreground pixels.Closing is done by performing dilation followed by erosion (M•S = (M ⊕ S) S).Closing is similar in some ways to dilation in that it tends to enlarge the boundaries of foreground regions in an image and shrink background color holes in such regions, but it is less destructive of the original boundary shape.

Morphological grayscale reconstruction.
Morphological reconstruction allows for the isolation of certain features within an image based on the manipulation of a mask image, M and a marker image, Y .It is founded on the concept of geodesic transformations, where dilations or erosion of a marker image are performed until stability is achieved (represented by a mask image).The marker must be the same size as the mask and its elements must be less than or equal to the corresponding elements of mask [17].
The geodesic dilation used in the reconstruction process is performed through iteration of elementary geodesic dilations, until stability is achieved. ( The elementary dilation process is performed using standard dilation of size one followed by an intersection. (2. 2) The operation in (2.2) is used for elementary dilation in binary reconstruction.In grayscale reconstruction, the intersection in the equation is replaced with a pointwise minimum [13].

Threshold decomposition.
An image M is a mapping from a finite rectangular subset D M of discrete plane Z 2 into a discrete set {0, 1,...,N − 1} of gray levels.Gray level images are often regarded as functions or mappings [14].Let ψ be an increasing binary morphological operator, that is, (2.3) In order to extend the binary transformation ψ to grayscale images M, taking their values in {0, 1,...,N − 1}, it suffices to consider successive thresholds T k (M), for k = 0 to N − 1.
(2.4) Equation (2.4) constitutes the threshold decomposition of M [5].As shown in Figure 2.1, these sets satisfy the following inclusion relationship.
(2.5) When applying the increasing operation ψ to each of these sets, their inclusion relationships are preserved [14].Hence, ψ can be extended to grayscale images as follows: Grayscale erosion with a flat disc shaped structuring element will generally darken the image.Bright regions surrounded by dark regions shrink in size, and dark regions surrounded by bright regions grow in size.Small bright spots in images will disappear as they get eroded away down to the surrounding intensity value, and small dark spots will become larger spots.Figure 2.2 shows a vertical cross section through a grayscale image and the effect of erosion using a disc shaped kernel.Grayscale dilation with a flat disc shaped structuring element will generally brighten the image.Bright regions surrounded by dark regions grow in size, and dark regions surrounded by bright regions shrink in size.Small dark spots in images will disappear as they get filled in to the surrounding intensity value.Small bright spots will become larger spots.Figure 2.3 shows a vertical cross section through a grayscale image and the effect of dilation using a disc shaped S. Dinesh  structuring element.Grayscale opening can be used to select and preserve particular intensity patterns while attenuating others.Bright features in an image that is smaller than the kernel have been greatly reduced in intensity, while the larger features are unchanged in intensity.On grayscale functions, opening reduces the sharp peaks and may remove them.Grayscale closing can be used to select and preserve particular intensity patterns while attenuating others.Dark features in an image that is smaller than the structuring element are largely filled in to the same colour as the surrounding bright features.On grayscale functions, closing fills valleys and may completely fill the shallow ones [6].

Linear DEM smoothening
The two Linear DEM smoothening methods discussed in this section are Gaussian blurring and Mean Smoothening.The Gaussian blurring operator is a 2D convolution operator that is used to remove detail and noise from images (Fisher et al. [6]).The kernel used to implement the operator represents the shape of a Gaussian ("bell shaped") hump.In 2D, an isotropic Gaussian has the following form: This distribution is shown in Figure 3.1.The idea of Gaussian blurring is to use this 2D distribution as a point-spread function, and this is achieved by convolution.Since the image is stored as a collection of discrete pixels we need to produce a discrete approximation to the Gaussian function before we can perform the convolution.In theory, the Gaussian distribution is nonzero everywhere, which would require an infinitely large convolution mask, but in practice it is effectively zero more than about three standard deviations from the mean, and so we can truncate the mask at this point.Figure 2.3 shows a suitable integer valued convolution mask that approximates a Gaussian with a σ of 1.4.Using the kernel shown in Figure 3.2, DEM smoothening using Gaussian blurring can be performed using standard convolution methods.1/9 1/9 1/9 1/9 1/9 1/9 1/9 1/9 1/9 The idea of mean smoothening is to replace each pixel value in an image with the mean ("average") value of its neighbours, including itself.This has the effect of eliminating pixel values, which are unrepresentative of their surroundings.Mean smoothening is usually thought of as a convolution filter.Like other convolutions it is based around a kernel, which represents the shape and size of the neighbourhood to be sampled when calculating the mean.Often a 3×3 square kernel is used, as shown in Figure 3.3, although larger kernels can be used for more severe smoothing [6].

Nonlinear DEM smoothening methods
The two nonlinear DEM smoothening methods discussed in this section are morphological smoothening and morphological smoothening by reconstruction.
Morphological smoothening is done by performing opening followed by closing.Grayscale morphological opening is used to darken small bright areas and to reduce sharp peaks in images.Grayscale morphological closing is used to brighten small dark areas and to fill valleys in images [5].In DEM processing, opening is used to remove the peaks of a DEM and closing is used to remove the pits of a DEM.The net result is the peaks and pits in the DEM are removed.

Morphological smoothening by reconstruction.
The biggest disadvantage of morphological smoothening is the closing operation can create big pits by filling the narrow valleys while the opening operations can create big peaks by enlarging valleys [6].Morphological reconstruction can be used to maintain the peak and pit removal effect while avoiding the valley filling and enlargement.However, morphological reconstruction cannot be used on a closed image, as the closed image (the marker) can contain elements that are larger than the corresponding elements in the original image (the mask).The steps involved in morphological smoothening by reconstruction are: Step 1: Perform opening on the original DEM This step removes the peaks on the DEM.
Step 2: Perform reconstruction on the opened DEM using the opened DEM as the marker and the original DEM as the mask This is to avoid the valley enlargement caused by the opening process.
Step 3: The DEM obtained from Step 2 is inverted The peaks in the inverted DEM are pits in DEM obtained from Step 2.
Step 4: Performing greyscale opening on the inverted DEM This step removes the peaks in the inverted DEM, which is equivalent to pit removal in the DEM from Step 2. Since the elements in the opened inverted DEM (the marker) are less than or equal to the corresponding elements in the inverted DEM (the mask), morphological reconstruction can be used.
Step 5: Implementing morphological reconstruction on the opened inverted image using the DEM from Step 3 as the marker and the DEM from Step 4 as the mask This is to avoid the valley enlargement caused by the opening process.
Step 6: Inverting the reconstructed DEM This process performs pit and peak removal, but avoids the valley filling and enlargement effects of opening and closing.

Case study
The DEM in Figure 5.1(a) is a 250K : 1 DEM of Salt Lake City, USA.The DEM was obtained from http://edc.usgs.gov/geodata.Gaussian blurring, mean smoothening, morphological smoothening and morphological smoothening by reconstruction are performed on the DEM at levels of 5.The channel networks of the resultant DEMs are extracted using the algorithm developed in [4] and their fractal dimensions are computed using the box dimension method, which is explained in [10].As observed from the results in Figure 5.2, Gaussian blurring and mean smoothening produce results that are almost similar.Morphological smoothening initially produces results that are similar to Gaussian blurring and mean smoothening.However, when the opening process in morphological smoothening starts removing objects in the DEM, its fractal dimension value drops significantly.For real DEMs, morphological smoothening by reconstruction is able to produce results that are similar to mean smoothening and Gaussian blurring until the opening process in the morphological smoothening removes all the objects in the DEM.

Conclusion
Morphological smoothening by reconstruction is able to produce results that are almost similar to Gaussian blurring and mean smoothening, until the opening process in the S. Dinesh and P. Radhakrishnan 9 morphological smoothening removes all objects in the DEM.Since for most real DEMs, the opening step only removes all objects in the DEM after a large number of operation levels (190 and above), morphological smoothening by reconstruction can be considered a valid DEM smoothening tool for real DEMs.
One of the advantages of morphological smoothening by reconstruction is that it provides the user with more flexibility in choosing the threshold level of peak and pit removal by varying the number of opening and closing steps in each smoothening step and type step.
Our future work plan is to develop an algorithm to avoid the drop in the fractal dimension value caused by the removal of all objects in the DEM by the opening process in morphological smoothening.We also plan to design a more robust and efficient DEM smoothening tool, by making use of some of the more advanced mathematical morphological operators, including, area morphology and morphological filtering.This work will also be extended to study the fractal geomorphometric properties of simulated and real DEMs.

Figure 5 . 1 .Figure 5 . 2 .
Figure 5.1.Implementation of the DEM smoothening algorithms on the DEM of Lake Mary (a) The DEM of Lake Mary, (b) The channel network of the DEM.
and P. Radhakrishnan 5