# An Improved Skewness Balancing Filtering Algorithm Based on Thin Plate Spline Interpolation
## 2 Materials and Methods
### Principle of the Improved Method  
Although skewness balancing can obtain a promising filtering result in flat areas, it will cause a larger omission error in abrupt terrain. This effect is because there is an implicit assumption for skewness balancing that all ground points are located above the objects. However, this assumption is obviously unreasonable in undulating areas, such as mountains, cliffs, and ridges. As shown in Figure 4, some ground points are higher than some nonground points (building roof) in abrupt terrain environments. Thus, according to the principle of the skewness balancing, those ground points with higher elevations will be misclassified as nonground points. Meanwhile, due to the undulations of the ground, some nonground points with lower elevations will also be accepted as ground points. Therefore, the classical skewness balancing algorithm cannot perform well in these kinds of areas. As shown in Figure 5, although both ground points and nonground points have elevation jumps, that effect will be changed if we build an interpolated surface based on the ground seeds. The ground seeds are points with the lowest elevations in local areas. Thus, the interpolated surface can reflect the fluctuation of the terrain. Although the elevations of some ground points are higher than those of nonground points, their corresponding elevation differences to the interpolated surface are totally different. Figure 4: Undulation terrains with ground points located above objects. Figure 5: Principle of the improved method. As shown in Figure 5, it is easy to find that the elevation differences between the ground points and the interpolated surface are obviously smaller than the elevation differences between the nonground points and the interpolated surface. This finding can avoid the misclassification of ground points and nonground points in undulating terrains. Thus, the improved method can be realized by modifying Equation (3) as:  
\[sk=\frac{1}{N\sigma^{3}}\sum_{i=1}^{N}(\Delta_{i}-\mu)^{3}, \tag{7}\]  
\[\Delta_{i}=Z_{i}-Z_{iter}, \tag{8}\]  
where \(\Delta_{i}\) is the elevation difference of each point to the interpolated surface, and \(Z_{iter}\) is the interpolated elevation of each point. The key procedures for realizing the improved method are ground seed selection and interpolated surface construction. Ground seeds are the points with the lowest elevations and can be selected using pseudo-grids as presented in the Section 2.1. Note that the cell size used here should be larger than the largest size of an object. As such, we can guarantee that the seeds are true ground points. A thin plate spline (TPS) was adopted to build the interpolated surface, since it can produce smooth, oscillation-free surface [40]. Moreover, in comparison with the other three interpolation methods, namely, ordinary Kriging, inverse distance weighting, and triangular irregular network, the TPS performs the best for fitting the LiDAR data and is the only interpolator that does not erode the ground surface [41]. For a set of ground seeds \((X_{i},Y_{i},Z_{i})\), \(i=1,2,\cdots,n\), the TPS can conduct interpolation to build a smooth surface \(f(X,Y)\), which is denoted as Equation (9):  
\[f(X,Y)=a_{0}+a_{1}X+a_{2}Y+\sum_{i=1}^{N}w_{i}R(\|(X_{i},Y_{i})-(X,Y)\|), \tag{9}\]  
where \(R\) is a kernel function that is defined as \(R(r)=r^{2}\log r\), \(r\) is the Euclidean distance, and column vector \(a(a_{0},a_{1},a_{2})\) and \(w\) can be estimated according to Equation (10):  
\[\left[ \begin{array}{c}KP\\ P^{T}0\end{array} \right]\left[ \begin{array}{c}w\\ a\end{array} \right]=\left[ \begin{array}{c}T\\ 0\end{array} \right], \tag{10}\]  
where \(K\) is a square matrix with a size of \(N\times N\), \(K_{ij}\) is equal to \(R\left(\|(X_{i},Y_{i})-\left(X_{j},Y_{j}\right)\|\right)\), \(P\) is an \(N\times 3\) matrix with the \(i\)-th row defined as \((1,X_{i},Y_{i})\), and \(T\) is a tensor of the ground seeds. The interpolated surface is obtained by minimizing the bending energy given by Equation (11) [42].