IMPLEMENTING SIFT AND BI-TRIANGULAR PLANE TRANSFORMATION FOR INTEGRATING DIGITAL TERRAIN MODELS

Since their inception in the middle of the twentieth century, Digital Terrain Models (DTMs) have played an important role in many fields and applications that are used by geospatial professionals, ranging from commercial companies to government agencies. Thus, both the scientific community and the industry have introduced many methods and technologies for DTM generation and data handling. These resulted in a high volume and variety of DTM databases, each having different coverage and data-characteristics, such as accuracy, resolution, level-of-detail amongst others. These various factors can cause a dilemma for scientists, mappers, and engineers that now have to choose a DTM to work with, let alone if several of these representations exist for a specified area. Traditionally, researchers tackled this problem by using only one DTM (e.g., the most accurate or detailed one), and only rarely tried to implement data fusion approaches, combining several DTMs into one cohesive unit. Although to some extent this was successful in reducing errors and improving the overall integrated DTM accuracy, two prominent problems are still scarcely addressed. The first is that the horizontal datum distortions and discrepancies between the DTMs are mostly ignored, with only the height dimension taken into account, even though in most cases these are evident. The second is that most approaches operate on a global scale, and thus do not address the more localized variations and discrepancies that are presented in the different DTMs. Both problems affect the resulting integrated DTM quality, which retains these unresolved distortions and discrepancies, resulting in a representation that is to some extent inferior and ambiguous. In order to tackle this, we propose an image based fusion approach: using the SIFT algorithm for matching and registration of the different representations, alongside localized morphing. Implementing the proposed approach and algorithms on various DTMs, the results are promising, with the capacity correctly geospatially align the DTMs, thus reducing the mean height difference variance between the databases to close to zero, as well as reducing the standard deviation between them by more than 30%.


BACKGROUND
Digital Terrain Model (DTM) is a digital representation of the bare earth used for a variety of processes and applications, such as GIS, city modelling, land use studies, drainage control, geologyto name a few.Since it is a digital model, a discrete one, one of the main goals is that it will be "realistic" (in terms of similarity, precision, and accuracy), i.e., serve as a reducednonetheless reliable and accurate -representation of our reality (Meyer, 2012).
In an earlier inception, DTMs were generated using scanning and digitization of contour lines in existing topographic maps.As a result, the accuracies of these models were tied to both the source data and the digitization process.Nowadays, modern data acquisition technologies have made a vast evolution, with techniques like laser scanners, radar, and photogrammetry, which can provide an improvement in terms of accuracy, resolution, and coverage, i.e., significantly improve the modelling stage.Today, large coverage models exist, having high quality, sometimes covering the entire earth, such as SRTM (Shuttle Radar Topography Mission), ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) and WorldDEM.
The fact that several different models exist for the same (coverage) area, possesses a dilemma for experts and authorities making use of DTMs, whereas mostly they will have to choose a single DTM to work with, overlooking the other sources.In most cases, the decision is made based on a direct comparison between the databases and their characteristics, while choosing the most up-to-date or accurate source.This, however, will not always be the appropriate and suitable decision, as each DTM has its own characteristics, thus important and relevant information might be overlooked by choosing a single DTM, both on a global and local level.This fact intensifies in cases where it is not an easy task to choose the most 'appropriate' DTM to use since all DTMs are relevant in terms of up-to-date, accuracy, and resolution.Figure 1 (right).SRTM presents a more rugged terrain, having more detailed features; however, it also has some data-holes.
Alternatively, one can suggest an optimized solution that makes use of the various available models (sources).However, as each source presents different types of detailing, scale, resolution, and accuracy (both horizontal and vertical, as well as local)such a solution is not a simple task.An optimized solution can suggest the use of all the models (data) simultaneously in the form of DTM manipulation, such as data fusion, merging, and rubber shitting.Such a solution aims to take the various existing DTMs, and use the cumulative data to present it as one cohesive continuous model (e.g., Katzil and Doytsher, 2005;Dalyot, 2010), with the aim to improve the overall quality and representation of the resulting fused DTM (when compared to any of the DTMs) by correctly using and modeling the existing data from the multiple sources.
Still, none of the above techniques consider the horizontal differences between the DTMs, thus assuming that there are no planar data distortions between the databases.In most situations, however, this is not the case, as a DTM acquired from one source will not only show discrepancies in the height values to a DTM from another source but also show horizontal discrepancies.Accordingly, a planar (2D) registration of the different sources should take place prior to fusion, introducing more complex challenges to the problem, as distortions between the DTMs are evident not only on a global level but also on a local one, which are needed to be addressed in order to acquire reliable and robust integration implementation and accurate results.
The more classical approach to registration is a vector-based approach, also known as geometric registration.This approach aligns two point datasets using their geometric characteristics.Some methods in this family include the Hausdorff distance (Huttenlocher et al., 1993;Dubuisson and Jain, 1994), The Iterative Closest Point (ICP) algorithm (Besl and McKay, 1992;Eggert and Dalyot, 2012), and The Coherent Point Drift (CPD) algorithm (Myronenko and Song, 2010).All these approaches, however, are very susceptible to noise and outliers (Ben Haim et al., 2015).
A different approach is a raster-based registration, which applies image-processing techniques for solving dataset registration problems.In this approach, a grayscale image represents the DTM where each pixel value is the height data.For point matching, the most prominent method is the Scale Invariant Feature Transform (SIFT) (Lowe, 1999;Lowe, 2004) and Speeded-Up Robust Feature (SURF) (Bay et al., 2006).The main problem with the image transformation approach is that the transformation occurs from one image to another, while in the case of DTMs, we aim to find an in-between state that represents both datasets simultaneously.One interesting approach to this is image morphing, where an initial image is transformed seamlessly into another, where each step in the process is defined in advance and calculated.Techniques in this field include mesh warping (Wolberg, 1990), field morphing (Beier and Neely, 1992), radial basis functions (Arad et al., 1994), thin plate splines (Lee et al., 1994;Litwinowicz and Williams 1994), energy minimization (Lee et al., 1996), and multilevel free-form deformations (Lee et al., 1995).
In this research, we propose a method for the transformation of one DTM to the other as a first stage to data integration, using image processing tools, namely SIFT and image morphing.
Registration and matching are achieved on a local level to monitor and quantify local discrepancies between the DTMs, while still ensuring that the models remain continuous and intact when fused.

METHODOLOGY
As this research focuses on grid type datasets (i.e., raster format), image processing approaches are investigated and analyzed.As such, the following four phase workflow is conducted: (1) interest point detection, (2) point matching, (3) TIN construction, and (4) DTM morphing.

Interest Points Detection
A DTM grid data is a datatype constructed from points with regular intervals between them, with each point having a certain height data attached to it.This type of data is equivalent to a greyscale image data, with a certain pixel size that represents the intervals and a pixel color that represents the height.Therefore, the tool chosen for interest point detection is the Scale-Invariant Feature Transform (SIFT).The SIFT algorithm uses a Difference of Gaussian (DoG) calculation to detect key points within an image.The DoG method takes an image and resamples it to multiple scales, called octaves.For each of those octaves two or more Gaussian filters are calculated using different σ values.Following the DoG process, a local extremum is detected by comparing a pixel in an image to its 8 neighbors, as well as to the 9 closest pixels on both the scale below and the scale above.In addition, a final test is performed on the detected key point -edge removal using a 2x2 Hessian matrix.This edge removal is important as the DoG process is very sensitive to edges, and can lead to a high noise in the process if not preliminary dealt with.

Point Matching
Following the extraction of interest points from two datasets, we aim to match homologous interest points that exist in both datasets; henceforth, a point matching process needs to take place.In order to find the corresponding point-pairs, the point descriptor result from the SIFT algorithm is used.This descriptor is calculated using two steps.In the first step, an orientation is assigned to each key point using the gradient calculation from the neighboring points; this gives each point a description for both scale and direction.The second step creates a more robust description that takes into account not only the point itself, but also its neighborhood.A 16x16 window is used around the point, with each point in that window also assigned an orientation.In addition, a 16 sub-blocks of 4x4 are picked, with an orientation calculation of each sub-block.The final result is a 128 bin descriptor.
It is important to note that even after the SIFT matching, outliers might still exist.As such, an outlier elimination process is implemented.This process takes into account the fact that different DTMs of the same area should have almost the same orientation value, and that the Euclidean distance between corresponding point-pair should be similar to the one calculated for all pairs.Thus, based on a statistical calculation (average and standard deviation), in case a Euclidean distance of a certain point-pair is not in the acceptable statistical range, this point-pair is considered an outlier and is removed.After removal, a new average value is calculated, and the process continues until no outliers are detected.

TIN Construction
After finding corresponding point-pairs in the two DTM databases, a Triangular Irregular Network (TIN) is constructed for both datasets.These triangles represent a local area in the data, which will allow the handling of localized areas within the DTMs.Thus, instead of handling the entire surface by a single averaged matching, we divide it to multiple intendent matched surfaces.The TIN itself is constructed using the Delaunay Triangulation (DT) method.
It is important that both TINs in the two DTMs are identical, i.e., each triangle-surface is constructed relying on the same interest points (required for the morphing stage in Section 2.4), such that a single TIN is constructed for one of the DTMs, and its topological structure is copied to the other one.In addition, in order to have a TIN structure that covers the entire coverage area represented in the DTMs, the corner points in both DTMs are artificially chosen as point-pairs.This, however, will affect the accuracy of the process alongside the databases' edges, but is important to keep to handle the complete data represented in the DTMs.
Another important aspect of this process that needs to be taken into consideration is the area of a single triangle in the TIN structure.In case a small triangle exists in the TIN structure, it means that relatively close interest points exist, such that a small triangle can be considered as redundant, introducing noise to the overall process.Henceforth, a supplementary process is implemented designed to eliminate small triangles, where a minimum threshold area value is set.This value is a function of two parameters, the resolution of the DTM and the density of the interest points found; the higher the resolution is the smaller the area value allowed will be.The process itself inspect each triangle in the set, when a triangle is found that does not meet the threshold, its vertices are removed from the set and are replaced with a new point calculated based on the vertices' centroid, and a new DT calculation takes place, as depicted in Figure 2.This process will continue until no more triangles are found that meet the requirement.

DTM Morphing
The morphing process consists of creating a local transformation between the two databases that hinges on the TIN of one database, to the TIN on the second database.Basically meaning that a point that exists in a certain triangle in one DTM is transformed to its new location in the corresponding triangle in the second DTM database using a local transformation between the triangles, as depicted in Figure 3.To define the new location, an affine transformation is calculated for each triangle that relies on the vertices' coordinate values, as depicted in Equation 1.
where XD, YD is the destination coordinates, XS, YS is the source coordinates, and a, b, c, d, e, f are the local transformation parameters.The values of the six unknown parameters are calculated using an equation system consisting of six equations, two for each point-pairs that define the corresponding triangles.This triangle transformation process is akin to a rubber sheeting process found in different GIS tools.The major difference here is that the height value that is different for both triangles' location is also taken into consideration, and is calculated using local interpolation within the triangle, depicted in Equation 2.
where HP1, HP2, HP3 are the heights of the triangle vertices, and d1, d2, d3 are the Euclidean distances from the corresponding vertices to point Pi.After calculating the heights of the point in both triangles, the bias is calculated as the difference between them.Thus, the final point height in the destination will be the height of the source point plus the bias, as depicted in equation 3.This is calculated for all DTM points, achieving a more accurate local discrepancies monitoring, as opposed to a global one.The main idea behind this local morphing process is the creation of a seamless transition between two DTMs, which is a first stage before integration, while maintaining the continuous structure of both DTMs involved in the process.

EXPERIMENTAL RESULTS
An area covering close to 15 sq.km in a mountainous area in northern Israel was analyzed.To evaluate our methodology, two DTMs were used for that area, one from the SoI (Survey of Israel), and one from NASA's SRTM.Both DTMs are depicted in Figure 4. SIFT automatically detected 48 interest points in the SoI DTM, and 50 in the SRTM DTM, with 38 point-pairs.The SIFT results are depicted for both DTMs in Figure 5.After eliminating outliers and duplicates, a total of 31 point-pairs remained for the TIN process, which calculated 64 triangles.A 100 sq.m triangle area threshold was used to eliminate redundant data.
Based on the above, the comprehensive geospatial local morphing process was implemented, transforming points from one DTM to their corresponding and accurate location in the other.Table 1 and Table 2 depict the statistical results for the height difference values per location for the entire DTMs.The Mean value of the height difference calculated for more than 20,000 points was reduced to almost zero, where both the RMSE and STDEV values were reduced by close to half.This means that the proposed methodology is able to accurately align the DTMs, and thus can be used to integrate both DTMs.Observing the values in the tables, it is clear that the process correctly identified the local geospatial correspondence(s) of both DTMs, precisely modeling any local discrepancies exists between them, accurately positioning them geospatially.This is presented in Figure 6, showing that in general the large height discrepancies (in white) are removed after transformationthe right image (after morphing) depicts much less white areas, when compared to the left image (before transformation).The white areas that still remain are on the borders, and are the result of the corner points that are artificial point-pairs.Another area from both databases was analyzed with the same coverage area, depicted in Figure 7. SIFT has detected 58 interest points in the SoI database and 80 in the SRTM database, and out of those 43 matched point-pairs were found.After eliminating outliers and duplicates, a total of 38 pairs remained for the TIN process.A 50 sq.m triangle area threshold was used.Table 3 and Table 4 depict the statistical results for the height difference values per location for the entire DTMs.These results are not as good as the previous area due to higher signal-to-noise ratio.For the RMSE and STDEV values we can see that the initial values are lower than the previous area, but the final improvement was not as effective, which is a result of the transformation not being as accurate due to the existing noise.The SRTM database in particular was noisier than the SOI database, as can be seen by the amount of interest points found in the SIFT process, as opposed to the SOI one.This means that the algorithm is less effective when noisy or rugged surfaces are involved, but can still produce reliable results, with a meaningful statistical improvement.

DISCUSSION
Based on the preliminary experiments presented here, the results are very promising.The most notable improvement is the reduction of the mean height difference values between the DTMs, which have been reduced to almost zero, implying that the process manages to monitor local discrepancies, and in doing that manages to co-model and align both databases.Another important issue is that the noisier the data is the more effective the triangle removal, as shown in the second experiment, where the mean difference reduction was 56% as opposed to 36% when no triangles were removed.However, more experiments are needed, with other DTM sources and topographic representations, to find the optimal adaptive threshold value for triangle removal.Overall, the proposed methodology and algorithms are robust and effective.
Future work will try to extend this methodology and improve the algorithms to higher resolution DTMs, as well as to DSM type models.In addition, other problems related to DTM fusion are planned for investigation, such as filling missing data.The overall objective is to extend this process into an automated fusion process, where the local transformation (geospatial morphing) quantifications (values) will be used to transform both DTMs simultaneously 'towards' each-other to create a seamless fused DTM product.
depicts an example of an area represented by two DTMs, showing a different representation and characterization of the topography.

Figure 1 .
Figure 1.Shaded relief representation of the same area by two distinctive DTM sources: SRTM (left) and Survey of Israel(right).SRTM presents a more rugged terrain, having more detailed features; however, it also has some data-holes.

Figure 2 :
Figure 2: An example of a small (redundant) triangle removal: the triangle in red (top) is smaller than the predefined threshold, and thus is removed, where its vertices' centroid is defined as a new point in the new DT process (bottom).

Figure 3 .
Figure 3.A transformation example between two corresponding triangles: a point in the left triangle is locally transformed to its corresponding position on the right triangle.

Figure 6 .
Figure 6.TIN construction and shift vector for each point in area #1.Grey values depict height difference value per point: before morphing (left) and after (right).

Table 1 .
Height difference values for area #1 before and after morphing, with no triangle threshold.

Table 2 .
Height difference values for area #1 before and after morphing, with triangle threshold.

Table 3 .
Height difference values for area #2 before and after morphing, with no triangle threshold.

Table 4 .
Height difference values for area #2 before and after morphing, with triangle threshold.