Estimating Horizontal Displacement between DEMs by Means of Particle Image Velocimetry Techniques

To date, digital terrain model (DTM) accuracy has been studied almost exclusively by computing its height variable. However, the largely ignored horizontal component bears a great influence on the positional accuracy of certain linear features, e.g., in hydrological features. In an effort to fill this gap, we propose a means of measurement different from the geomatic approach, involving fluid mechanics (water and air flows) or aerodynamics. The particle image velocimetry (PIV) algorithm is proposed as an estimator of horizontal differences between digital elevation models (DEM) in grid format. After applying a scale factor to the displacement estimated by the PIV algorithm, the mean error predicted is around one-seventh of the cell size of the DEM with the greatest spatial resolution, and around one-nineteenth of the cell size of the DEM with the least spatial resolution. Our methodology allows all kinds of DTMs to be compared once they are transformed into DEM format, while also allowing comparison of data from diverse capture methods, i.e., LiDAR versus photogrammetric data sources.


Introduction
A lot of geomatic applications use DTMs as the basis to derive other products such as cartography, hydrological features, engineering solutions, urban planning, environmental studies, time series analyses, reports, statistics, etc.The precision and accuracy of DTMs, required for deriving the above-mentioned products, however, may vary substantially.The data capture methodology is one key to explain the precision and accuracy of DTMs and how they influence the derived product.There is a need for reliable measures of accuracy allowing comparative evaluations of DTMs from different sources, e.g., LiDAR versus photogrammetry.Most studies thus far focus on the height variable [1,2].As the easiest variable to statistically compute and interpret, such studies focus on vertical accuracy.Some authors [3,4] studied the DEM accuracy as a function of some features (topographic index, hydrological objects).There are no reports to date of horizontal accuracy of DTMs or DEMs linked with the derived features, most likely because any such measure, if it exists, would hardly be applicable to other features [5].More recently, algorithms were defined to link the differences between horizontal accuracy and the derived parameters of features; such algorithms could be based on a triangulated irregular network (TIN) [6,7] or DEM grid models [8].The first approach generalizes the separation between homologous contours as the measure of horizontal discrepancy for the whole DEM, whereas the second compares several cells in two or more DEMs to estimate the horizontal discrepancy.
We can define the vertical accuracy for a DEM (grid structure composed of rows and columns) as the magnitude indicating how close the height of a point is, located by its horizontal position (row and column), on a DEM with respect its true value.If we have a reference DEM (RefDEM) as the ground truth we could estimate the vertical accuracy from any other DEM (CompDEM) just by computing the differences between the corresponding cell's height from both DEMs.The vertical accuracy can be seen as the matrix: accMat = RefDEM ´ComDEM.Horizontal accuracy is a little bit more complicated to understand because it may be difficult, given a point (RefP) on a RefDEM, to find its homologous point (HomP) on a CompDEM, except, for example, when RefP and HomP are a building corner, but that is not usually the case.We could understand the horizontal accuracy of a CompDEM with respect to its RefDEM; imagining the case where the CompDEM is identical to the RefDEM, then their contours would be the same.However, if we move the identical CompDEM horizontally, then the homologous contours from the CompDEM and RefDEM would not match and the discrepancy between the homologous contours indicates the horizontal discrepancy, which, consequently, is a measure of the horizontal accuracy.It is possible for each cell in the CompDEM to be horizontally displaced in any direction, so that homologous contours cross each other as shown in Figure 1, where the area between homologous contours has been filled with gray color.Both vertical and horizontal accuracy are important features to keep in every DEM because error in vertical and horizontal positioning could have bad consequences when, for example, computing flood risk areas.If we miscompute flood areas as a consequence of inaccurate DEMs, it would be possible to authorize, by error, the construction of some houses in flood risk areas.An illustration for both errors (vertical and horizontal inaccuracy) is shown in Figure 2: on the left the real flood area has been computed, and the image on the upper right shows an error due to vertical inaccuracy (the CompDEM altitude is higher than the RefDEM altitude and an urban area is flooded); the image on the bottom right shows the error due to horizontal inaccuracy (the CompRef is moved to the left and up direction and an urban area is flooded).Another example of horizontal accuracy importance comes from the ability of LiDAR data to produce a highly detailed drainage network [9] and it would be interesting to know the variation in the streams if computed with different DEM sources.
Landslide is another field where horizontal displacement has to be computed.There are several approaches to get such calculations, and one of them is based on series of images taken in different epochs, e.g., [10,11].Those studies base their computation on minimizing a correlation function between a local window in the reference image and its homologous window in the compared image, according to the method explained by [12].
In our present study we use an approach similar to that of [8], but apply the particle image velocimetry (PIV) algorithm.PIV comes from a distant research field, more commonly used in studies involving aerodynamics [13], transonic flows [14], wind tunnels [15], or aeroacoustic predictions [16].PIV can be used to determine the velocity field generated by a flow, tracking the particles inside that flow by repeated photography of the same area; an algorithm then computes the displacement of a particle (or a set of particles) from one photograph to the next.Dividing the particle displacement (number of pixels multiplied by the pixel size on the ground) by the time between photos gives us the velocity.Figure 3 shows how two sets of particles are displaced from Photo 1 to Photo 2.  A DEM resembles a digital photograph in two ways: -Both of them are made up of a regular grid of cells.
-Each cell is represented by a numerical value.
For these reasons, the PIV algorithm used on photos can almost directly be applied to DEMs.The only difference is that RGB (Red, Green, Blue) color photos must be transformed to a gray scale in the case of PIV, while DEM does not require transforming the cell values.

Materials and Methodology
In order to test the ability of PIV in estimating the horizontal displacement between two DEMs, we adopted the dataset used in [8].It consists of four samples (corresponding to four different places) pertaining to terrains with different slopes; they are labeled as upland, hill and mountain according to the classification from [17].Each sample is composed of two DEMs from two different sources: "Instituto Geográfico Nacional" (IGN) [18] and "Instituto Cartográfico de Andalucía" (ICA) [19], Spanish cartographic agencies at the national and regional level, respectively.The data captured correspond to the years 2005 (IGN) and 2004 (ICA).The cell size for the ICA DEMs is 10 m, and for the IGN ones it is 25 m.This means that the IGN DEMs must be resized in order to have the same dimensions in homologous DEMs (we will refer to homologous DEMs to the sample composed for the ICA and the IGN DEMs from the same place).Once we resize an IGN DEM, its cell size is the same as the ICA one, 10 m.There are several interpolating methods that could be used for resizing the IGN DEM, taking into account the cell neighborhood, e.g., the nearest, the bilinear or the bicubic.We decided to use the bicubic one, because the interpolating method influence on the estimated horizontal displacement is small, according to the study shown in Section 3.1.
The names and the areas of each sample DEM are as follows: Seville (63.98 km 2 ), Huétor-Tájar (34.39 km 2 ), Colmenar (70.64 km 2 ) and Quéntar (51.85 km 2 ).Table 1 shows the classification and the slope value for each place.To compare the PIV-estimated horizontal displacement (PIVEHD) with a reference displacement, we compute what we called the "real horizontal displacement" (RHD), which will be considered as the ground truth.
The RHD was computed by digitizing, as a polygonal shape (linear polyline), the same hydrological feature (creek or river) from the two DEM sources (ICA and IGN).Because both polygonals are in the same reference system (ETRS89 Datum) and they do not coincide exactly, the horizontal displacement can be computed as the area comprised between the two homologous polygonals and divided by their mean length.Each hydrological feature was digitized by four different people to enhance RHD accuracy.For the sake of redundancy we compute the RHD for four different hydrological features in each DEM (denoted as hydrological features 1, 2, 3 and 4 in Table 2, distributed as shown in Figure 4).Table 2 shows the RHD between homologous hydrological features and their corresponding length; the last columns give the RHD mean value for each place, taking into account the four hydrological features.It may be argued that applying the PIV to the whole DEM-comparing the PIVEHD of the entire DEM with the RHD from linear hydrological features-is an over-generalization.Hence, we computed the PIVEHD inside a buffer 600 m wide surrounding each hydrological feature.This was compared to the corresponding RHD.The size and shape of the buffers can be seen in Figure 4.
Depending on the type of analysis needed, the PIV algorithm can be set with different parameters.A complete overview of PIV setting variability is found in [20].We opted to run the algorithm based on the cross-correlation function performed under Fast Fourier Transform (FFT), taking into account that displacement could entail translation (case A in Figure 3) or a mixture of translation and rotation (case B in Figure 3), as explained in [20].
We used PIVlab software [21] to arrive at parameters suitable for the above requirements.
The PIV method consists of defining an interrogation window in the first DEM (D) and displacing it over the second DEM (D 1 ) around the original position of the window in D (Figure 5).The PIV algorithm attempts to find the best match between images (in our case DEMs D and D 1 ) using the discrete cross-correlation function defined in Equation (1).
where H and H 1 respectively represent the cell height in D and D 1 .For each displacement (x, y) selected, we obtain a R I I px, yq cross-correlation value.If we select a shift range p´Mq ď x ď M, ´N ď y ď N we obtain a correlation plane of size p2M `1q ˆp2N `1q.The cell value for the correlation plane is the sum of the products of all overlapping cells in D and D 1 .A graphic illustration is shown in Figure 6.
The best match would be the highest cell value in the correlation plane, and therefore the displacement is the corresponding (x, y) shifts that generated such a high cell value.The cross-correlation function as expressed in Equation 1 has two weak points: 1. High computational time.

2.
It may produce different maximum correlation values for the same degree of matching because the function is not normalized.For example, areas with many (or brighter) particle images will produce much higher correlation values than interrogation windows with fewer (or weaker) particle images.
The solution to the first weak point is to compute the cross-correlation function as a complex conjugate multiplication of their Fourier transforms: R I I " p H x H 1 ˚, where p H and x H 1 are the Fourier transforms from H and H 1 .
The second weak point can be overcome using the cross-correlation coefficient, which normalizes the cross-correlation function and its definition as expressed by Equation ( 2): where The value µ H is the average of the template and is computed only once while µ H 1 px,yq is the average of H 1 coincident with the template H at position px, yq.
PIVlab [21] software offers the option of applying a multiple pass interrogation window that improves the PIV results in the sense that more particle matching is achieved in the successive passes.Each pass decreases the interrogation window size.We tried three hierarchy passes, where in the next pass the interrogation window was one half of the previous one.

Effect of the Interpolation Method for Resampling IGN DEM to ICA DEM on the PIVEHD
In order to study the influence of the resampling method to resize the IGN DEM to the size of the ICA DEM, we compute the differences in height for the two extreme cases, i.e., Quéntar (the steepest DEM) and Seville (the flattest DEM).If we denote a resampled DEM as R m x where x represents the DEM (Quéntar (Q) or Seville (S)) and m represents the interpolation method (nearest (near), bilinear (bil) or bicubic (bic)), we can compute the mean difference in height for a particular DEM and two different methods, although it is most informative if we compute the mean of the absolute value of the height difference.For example considering Quéntar and the bicubic and nearest method, such a mean value will be computed as follows: Table 3 shows the possible combination represented by Equation ( 6).A histogram representation referring to Quéntar for the three combinations mentioned in Table 3 is shown in Figure 7. From Table 3 and Figure 7 we observe than the differences between bicubic and bilinear methods are small and, when the terrain is flat like Seville, the differences between all the methods are also small.For this reason we start to compute PIV displacement for a steep area like Quéntar because if the differences between the methods are small, then in flat areas the differences will be smaller.We compute the PIVEHD using the three different IGN resampled DEMs (nearest, bilinear and bicubic) and compute the mean discrepancy in the absolute value for the module of the vector, i.e., if U m is the matrix containing the horizontal axis component for the displacement vectors, V m is the vertical axis component for the displacement vectors and m is the interpolation method (nearest, bilinear an bicubic), we can compute the effect in computing the PIVEHD using Equation ( 7), which considers the computation made between the bicubic and nearest interpolation methods: The parameters used for computing PIVEHD where pass 1:32 ˆ32; pass 2:16 ˆ16; pass 3:8 ˆ8.
The Table 4 shows the values for the differences between methods when computing PIVEHD and the corresponding standard deviation.
Quéntar 0.04 0.11 0.16 0.17 0.16 0.18 From Table 4 we observe that the effect on PIVEHD is not significant either computed by the bicubic or bilinear method (4 cm), and just the near method differs with respect to bilinear and bicubic method in the same magnitude (16 cm), which is not a big value considering the cell size for ICA (10 m) and IGN (25 m).So we decide to use the bicubic interpolation method in order to resample the IGN DEM to the ICA DEM size.

Results Analysis
In order to analyze several interrogation window sizes, we tested three settings classified as High, Medium and Small.The values for such a classification are as follows: -High: (pass 1:64 ˆ64; pass 2:32 ˆ32; pass 3:16 ˆ16).
A graphic example of applying the Medium interrogation window to the Colmenar region is shown in Figure 8.The displacement vectors are shown in green.
Table 5 shows the PIVEHD expressed in meters (Medium interrogation window setting) for each place (Quéntar, Colmenar, Huétor-Tajar and Seville), evaluated over each of their hydrological features (1, 2, 3 and 4).Although the PIVlab algorithm outputs are vectors, we compute their modules in order to know the magnitude of the displacement; the comparison shown in Table 5 refers to those modules.Remember, the horizontal displacement has been computed inside the buffer area surrounding each hydrological feature (creeks or river), as shown in Figure 4. Given in Table 5 are the respective RHDs in order to compute the error.The two last columns show the mean absolute (ABS) error for each place and the total mean ABS (error).The mean ABS (error) shown in Table 5 seems somewhat high, particularly if compared with the estimations obtained using the optical flow algorithm [8] which had a mean error of 1.13 m.Nevertheless, the optical flow algorithm features a parameter (α) that is similar to a scale factor.For this reason, if we find a suitable scale factor that reduces the errors then PIVEHD can be considered as scalable and a good estimator of displacement.We thus propose to compute the best scale factor using the least squares method, minimizing the function v T v obtained from Ax ´B " v, where A is the estimation from PIV, B is the RHD, and x is the scale factor.For example, considering Mountain terrain (Quéntar), the corresponding values are as follows: A " r0.92 0.94 1.16 0.95s T , B " r4.45 4.75 4.55 4.01s T and v " r3.53 3.81 3.39 3.06s T Equation ( 8) computes x and the standard deviation (s) x " ´AT A s " where dof means degrees of freedom.The scale factors for the High, Medium and Small interrogation windows after applying the least squares methods are shown in Table 6.The values are similar, ranging from 2.937 to 2.388.After applying the scale factor to the PIVEHD, the adjusted PIVEHD (APIVEHD) and the corresponding errors are shown (Tables 7-9), representing the High, Medium and Small interrogation windows.Comparing Tables 5 and 8 both corresponding to the Medium interrogation window, we note an improvement (Impv) of 2.22 m through application of the scale factor.This amounts to a substantial improvement (62.9%): Impv " errorWithoutScaFactor ´errorA f terScaFactor " 3.53 ´1.31 " 2.22 (10) Given that the cell size of the ICA DEMs is 10 m, and the best error corresponds to the Medium interrogation window (1.31 m), we may conclude that the PIV algorithm is a sound estimator of horizontal displacement between homologous DEMs: the error is around one-seventh of the corresponding cell size.On the other hand, if we consider the IGN cell size, the proportion decreases to one-nineteenth.
In general, no great differences are observed among the errors of the three categories (High, Medium and Small); the dispersion appears to decrease when the interrogation window decreases (see the values for the Total Standard Deviation ABS (error) in Tables 7-9).However, if we were to recommend an interrogation window, it would be the Small one, in view of its reduced Total Standard Deviation ABS (error) value.Notwithstanding, the Medium category offers the advantage of a shorter processing time.

Performance in Comparison to Previous Approaches
As mentioned in the abstract and the Introduction section, there are barely any algorithms analyzing horizontal displacement on homologous DEMs as a whole.Just the contour-based [7] and optical flow [8] approaches study the DEM horizontal variable.So we consider it necessary to compare the performance between both algorithms, and such a performance is studied considering the absolute error between the ground truth and the estimation achieved by every approach (contour, optical flow and PIV).An additional comparison concerning computation time has been added in order to assess the global performance.Absolute error and computation time are shown in Table 10.The computer used to test the mentioned algorithms had the following features: Intel i5 processor, 8 Gb RAM, 256 SSD hard drive, Windows 7 operating system.From Table 10, Mountain and Hill terrains have been similarly estimated by contour and optical flow algorithms, and both get better estimations than the APIVEHD one.On the computation time side, optical flow had the best performance; it is even better than APIVEHD although the differences are not very large (around 14% better than APIVEHD).Although the best horizontal displacement comes from the contour-based algorithm, its computation time is 4.28ˆlarger than optical flow and 3.71ˆlarger than APIVEHD.
Likewise, the error values obtained for Seville, a nearly flat region (2.6% slope), are noteworthy.On such a terrain type, the horizontal displacement errors increase quickly [8].In our case, the APIVEHD algorithm improves results for nearly flat terrains with respect to the optical flow approach, namely 1.9 m for the APIVEHD algorithm versus 3.18 m for the optical flow algorithm, although the contour algorithm keeps the best performance.We can see from Table 10 that the steeper the terrain the better the optical flow algorithm performs; nevertheless, in near-flat areas the APIVEHD approach achieves the best estimation.Taking into account both absolute error and computation time, we think that our new proposed algorithm APIVEHD gets a better tradeoff when we estimate the horizontal displacement in near-flat areas such as Seville.

Assessing the Influence of DEM Cell Size on the APIVEHD Algorithm Results
After computing the scale factor, the results obtained from the APIVEHD applied to our dataset seem to be suitable for estimating the horizontal displacement (Tables 7-9), but there is a variable that could be modified for suitable estimation yield by APIVEHD; this variable is the cell size ratio between the DEMs compared.Our study, resumed in Tables 7-9 compared IGN (25 m cell size) and ICA (10 m cell size).In order to evaluate the ability of APIVEHD to get a suitable estimation of horizontal displacement across different cell size ratios, we modified the cell size for ICA to 15 m and 25 m and then we applied the APIVEHD algorithm again, comparing those cell sizes (15 m and 25 m) against the IGN DEM (25 m).According to Section 3.1, the bicubic and bilinear resampling methods are similar and practically do not influence the results yielded by APIVEHD; for this reason, we select the bicubic method in order to compute the derived ICA 15 m and 25 m cell sizes in order to assess the effect of using different cell size ratios between DEMs when computing the horizontal displacement with APIVEHD.
For the sake of readability we just present in Table 11 the horizontal displacement between the following homologous DEMs whose cell sizes are inserted in their names: (1) IGN25m-ICA15m; (2) IGN25m-ICA25m.Note that IGN25m-ICA10m corresponds to Table 8.We selected Table 8 as the reference because the Medium interrogation window has the smaller mean absolute error.From the comparison of Tables 8 and 11 we do not observe large variations in the TOTAL mean ABS (error).The error estimating the horizontal displacement when we increase the ICA cell size by 5 m is 1.47, just 0.16 m if compared with original ICA cell size, i.e., the influence of the cell size ratio (IGN25m-ICA10m versus IGN25m-ICA10m) in the horizontal displacement estimation by APIVEHD is just 16 cm.The influence in the case of IGN25m-ICA10m versus IGN25m-ICA25m is 0.26 m.Both effects, 0.16 m and 0.26 m, seem to indicate that the APIVEHD algorithm remains robust to the variation in the cell size ratio.

Conclusions
In our study we compute the horizontal displacement between two homologous DEMs coming from different sources.This is important because the horizontal discrepancy in both DEMs puts some derived features in non-corresponding places (e.g., hydrological linear features).Such discrepancies may be estimated by using the PIV algorithm as proposed in this paper.However, a larger number of cases must be computed in order to completely generalize our results.Our findings are a novelty in that an algorithm from the distant disciplines of aerodynamics or water flows is here applied to the DTM field.
The PIV algorithm underestimates the real horizontal displacement, but it does so in a scalar sense, which means that it as a scale factor can be used to arrive at a suitable estimation.
Using PIV and applying the scale factor, the error is around one-seventh of the cell size of the smaller DEMs (the ICA DEMs) and around one-nineteenth in the case of the larger cell size DEMs (the IGN DEMs).We consider these to be acceptable error values.
One remarkable finding is that the PIV approach provides reliable horizontal displacement values for nearly flat terrain types, which are the most difficult ones to compute due to their horizontal displacement.
Future works along these lines will involve exploring PIV behavior in the horizontal displacement estimation between DEMs derived from LiDAR (very small cell sizes) and photogrammetry (cell sizes similar to the one used here).It would be interesting to determine whether the three interrogation window categories produce results that agree with the results presented here, or if new dimensions for the interrogation windows are necessary.In the case of LiDAR, because the cell size could be around a few centimeters, it would likely be better to use smaller features than hydrological ones, in order to compute the ground truth for the horizontal displacement, e.g., houses, buildings, fences, small dams, small pieces of roads.Another important point will be the computation time, because if we want to analyze a significant area, the file size will increase considerably, as will the PIV computation.Such a situation could need to parallelize some functions from the PIV software.

Figure 1 .
Figure 1.Homologous contours (coming from homologous DEM: RefDEM and CompDEM) crossing each other and filling the area with gray color.The black color is a river piece.

Figure 2 .
Figure 2. Importance of vertical and horizontal accuracy.

Figure 3 .
Figure 3. Particle displacement derived from the position captured in two different photographs.

Figure 4 .
Figure 4. Sample DEMs (four different places) and distribution of the digitized linear hydrological features.Surrounding each hydrological feature is a buffer 600 m wide.

Figure 5 .
Figure 5. Interrogation window on the first DEM (D) and a displacement over the homologous DEM (D 1 ).

Figure 6 .
Figure 6.Cross-correlation plane derived from the cross-correlation function using an interrogation window 4 ˆ4 in size.

Table 1 .
Slopes (%), heights (m) and category for each DEM from the dataset.

Table 2 .
Real displacement (RD) between homologous DEMs computed comparing linear hydrological features.

Table 3 .
Mean absolute differences between interpolation methods to resample the ICA DEM.

Table 4 .
Mean differences for PIVEHD modules and standard deviations comparing pairs of interpolation methods.

Table 6 .
Scale factor minimizing the PIVED errors.

Table 10 .
Horizontal displacement absolute errors and computation time obtained using contours, optical flow and PIV approaches.

Table 11 .
APIVEHD obtained by comparison between the homologous pairs DEMs IGN25m-ICA15m and between the homologous IGN25m-ICA25m.