SURFACE APPROXIMATION OF COASTAL REGIONS: LR B-SPLINE FOR DETECTION OF DEFORMATION PATTERN

: Geospatial data acquisition of terrains produces huge, noisy and scattered point clouds. An efficient use of the acquired data requires structured and compact data representations. Working directly in a point cloud is often not appealing. To face this challenge, approximation with tensor product B-spline surfaces is attractive. It reduces the point cloud description to relatively few coefficients compared to the volume of the original point cloud. However, this representation lacks the ability to adapt the resolution of the shape to local variations in the point cloud. The result is frequently that noise is approximated and that surfaces have unwanted oscillations. Locally Refined (LR) B-spline surfaces were introduced to face this challenge and provide a tool for approximating Geographic Information System point clouds. In our LR B-spline based approximation algorithm, iterative least-squares approximation is combined with a Multilevel B-spline Approximation to reduce memory consumption. We apply the approach to data sets from coastal regions in Norway and the Netherlands, and compare the obtained approximation with a raster method. We further highlight the potential of LR B-spline volumes for spatio-temporal visualisation of deformation patterns.


INTRODUCTION
Parametric spline surface approximation is an important step in reverse engineering to convert real object data to a computeraided design model (Raja and Fernandes, 2008), or to perform statistical testing for deformation based on the estimated parameters (Kermarrec et al., 2020). The surface fitting begins by recording a point cloud of an object or scene with, e.g., a contactless sensor such as the Terrestrial Laser Scanner (TLS). Surface approximation techniques of point clouds from Geographic Information System (GIS) data sets can be divided into (i) non-adaptive methods, for which the optimal surface is globally adapted in the approximation process, and (ii) locally adaptive methods. For the latter, each new iteration depends on the result of the previous approximation iteration (Wang 2009). Such procedures save computation time as the number of parameters to estimate is less high than for (i). The widely used adaptive surface fitting strategies using B-splines belongs to (ii), but is limited by the tensor product formulation of the nonuniform rational B-splines (NURBS): it does not allow for local refinement and leads potentially to overfitting where no further refinement would have been necessary. There are three main approaches for extending spline surfaces to support local refinement, see Dokken et al. (2019): • Hierarchical B-splines (HB-splines) were introduced in Forsey and Bartels (19888). Here tensor product Bsplines on the coarser level are removed and B-splines at the finer level added in such a way that linear independence is guaranteed.
• T-splines are described, e.g., in Sederberg et al. (2003). They denote a class of locally refined splines. The starting point for T-spline refinement is a tensor product B-spline surface with control points and meshlines (initial T-mesh) with assigned knot values. The refinement is performed by successively adding new control points in-between two adjacent control points in the T-mesh. They were used in a geodetic context by Kermarrec et al. (2020) to approximate a bridge.
• LR B-splines were developed by Dokken et al. (2009). The refinement is performed successively by inserting axis parallel meshlines in the mesh of knotlines. Each meshline inserted has to split the support of at at least one tensor product B-spline. In this contribution, we will focus on this approximation strategy which was shown to be flexible for approximating various point clouds (Skytt and Dokken, 2021). A toolbox called Gotools (Dokken et al., 2013) is freely available under https://github.com/SINTEF-Geometry/GoTools/wiki/ Module-LRSplines2D. Interoperability with existing GIStechnology is essential for the deployment of LR B-spline within GIS: LR B-spline surfaces can be converted to standard formats used in GIS-systems. They are, thus, particularly attractive within a GIS context where huge, noisy and scattered point clouds are common. Furthermore, LR B-spline volume is a promising tool for visualising and analysing spatio-temporal deformations, without having to manipulate large point clouds. The starting point of iterative surface approximation with LR Bspline is a tensor product B-spline surface. In the LR B-spline surface, extra degrees of freedom are inserted locally, where needed. The desired smoothness is maintained, and the growth in data volume is kept under control. The least square (LS) norm between the parametrised point cloud and the parametric surface is to be minimised. For gridded data, this approach is called surface skinning (Woodward, 1988). The LS method may reach its limits when the number of refinement iterations increases. In that case, it is convenient to use the Multilevel B-spline Approximation (MBA) developed by Lee et al. (1997) after a few initial step with the LS approximation. The MBA is based on a local approximation for which no equation system needs to be solved. The residuals of the data points obtained from the last fitted surface are recursively approximated using finer meshes. This procedure saves storage and was shown to be advantageous in case of scattered and/or non-gridded containing data gaps or outliers occur (Zhang et al. 1998;Skytt et al. 2015). In this contribution, we propose to introduce the LR B-spline and the MBA. To relate LR B-spline to surface representations already in use in GIS, we will compare our approximation results with raster representations and show the potential of LR B-spline volume for spatio-temporal deformation analysis of GIS dataset. We focus on point clouds from a bathymetry data set in Norway and a coastal region in the Netherlands, which are rather smooth and without edges.

MATHEMATICAL BACKGROUND
In this section, we propose to explain shortly what LR B-spline surfaces are starting from the concept of tensor product B-spline surfaces. The iterative approximation and principle of local refinement are shortly adressed. We finally focus on describing the combination between LS and MBA for iterative surface approximation. A bivariate tensor product spline space is made by the tensor product of two univariate spline spaces. Assuming that both univariate spline spaces have a minimal support B-spline basis, the minimal support basis for the tensor product B-spline space is constructed by making all tensor product combinations of the B-splines of the two bases. The minimal support basis for this spline space contains the tensor product B-splines

Surface approximation using LR
As in the univariate case, the basis has useful properties such as non-negativity and partition of unity. Figure 1 illustrates how a tensor product B-spline with the knot vector { } 0,1, 2, 3 in both parameter directions looks like. Spline surfaces are frequently represented using a bivariate minimal support tensor product B-spline basis. , Here , − are the surface coefficients and d is the dimension of the geometry space.

Local refinement
NURBS surfaces can be refined by an axis parallel meshline that crosses the complete parameter domain. LR B-spline can be considered as extension of NURBS, through insertion of shorter meshlines to refine the mesh locally. LR B-spline provides a richer set of refinements than HB-splines and T-splines according to Dokken et al. (2015). When we insert a meshline, at least one local B-spline should be split into two parts. The meshline must start and stop on meshlines in the other direction. In a first step, all tensor product B-splines completely traversed by that meshline are identified. Then these B-splines are split into two B-splines. Duplicates are merged to clean the data structure.
In the volumetric case, meshlines are replaced by mesh rectangles. Figure 2 shows the mesh of an LR B-spline surface. The support of one of the B-splines is illustrated by a yellow box. In a) the meshline goes from the boundary and across the yellow support, and continues on the other side. The support of more than one Bspline is traversed, implying that several B-splines will be split. In b) the support of exactly one B-spline is split. In c) the meshline transvers the edge of the yellow support and will not result in refinement of this tensor product B-spline. However, other tensor product B-splines are affected since the meshline is on top of an existing meshline: the continuity of the surface crossing the line is decreased. The meshline in d) is too short to trigger refinements.
If after this insertion, some tensor product B-spline still do not have minimal support, the refinement process will continue until all the tensor product B-splines have minimal support.

Iterative surface approximation using LS and MBA
The LR B-spline surface is iteratively refined using the principle of 2.1.2. Starting point: The starting point of the iterative surface approximation is a tensor product B-spline surface defined over a coarse mesh. A refinement is performed iteratively in the mesh cells, where one observation is associated with an error term higher than a given threshold TH . With a large threshold, few cells are refined whereas with a low threshold TH , the result will be a global refinement corresponding to the NURBS. The choice of TH depends on the level of accuracy needed, balanced by the computation time and number of surface coefficients A low threshold may lead to a noise fitting and should be avoided. It could result in ripples on the surface that are not the consequence of real pattern.
Here i P are the surface coefficients and i B are tensor product B-spline scaled to form a partition of unity. LS approximation: LS approximation is a global method where the following expression is minimized with respect to the surface coefficients i P over the entire surface domain: are the data points. The expression is differentiated and turned into a linear, sparse equation system in the number of surface coefficients to estimate. A LS approximation method will result in a singular equation system if there exists B-splines with no data points in its support.
As an LR B-spline surface is defined on a rectangular domain and a typical point cloud has a non-rectangular outline and may contain holes, parts of the surface domain will frequently lie outside the domain of the point cloud. To face that challenge, we add ( ) ( ) , J F x y , which is a smoothness term that enables a non-singular equation system even if this situation occurs. In our case, the smoothness term is an approximation to the minimization of curvature and variation of curvature in the surface. These originally intrinsic measures are made parameter dependent to give rise to a linear equation system after differentiation. A more detailed description can be found in Mehlum and Skytt (1997). The weight on the smoothness term is kept low to emphasize the approximation accuracy. In the examples of the next sections we took

MBA:
The multilevel B-spline approximation was introduced in Lee et al. (1997) for data interpolation and extended by Zhang et al. (1998) for approximating scattered data with HB splines. We assume that all observations are of the form ( ) ( ) , , , x y z x y , and that we have an initial surface approximation via the LS method. A mesh refinement is performed, and the coefficients of the new surface approximation are calculated individually using the residuals of the preceding approximation 1 j F − (see Fig.3).
The factor , , , , , x y z h x y z , and the method is directly applicable in the case of LR B-spline volumes. The procedure is summarised in a flow chart in Figure 3. The stop criterion is that the maximum error (locally) is within the threshold. A geospatial point cloud may represent a very rough terrain and contain noise and possible outliers. Thus, it is not necessarily beneficial to continue the iterations until all points are closer to the surface than TH. Typically, the process is stopped by the number of iterations. The accuracy improvement tends to drastically slow down when most points have a distance to the surface less than the threshold (Skytt and Dokken, 2021).

Goodness of fit
The definition of an optimal surface ia a controversial topic and should contain both heuristical (e.g. without "ripples", "smooth") and statistical components. Average and maximum distances are the standard way of measuring the error of a model. They are defined in the z-direction as = ∑ |̂− |⁄

=1
, and = max|̂− |, respectively where n is the number of data values, i z are the given values and ˆi z are the corresponding estimated values. The number of points outside the threshold is defined as the number of points for which the error term exceeds TH for a given iteration level. The degree of freedom or number of coefficients to be estimated for a given iteration k of the refinement is also an important quantity. For the computation time, we used a 64-bit operating system with 8 GB RAM and an Intel(R) Core(TM) i5-63000U CPU @ 2.40 Ghz, 5.8 GHz.

Raster representation
The raster representation is the most frequently used data format in GIS (Li et al., 2005). To make a raster representation of a point cloud, spatial interpolation is used to define values, for instance elevation or precipitation, in a regular grid. The Digital Elevation Model (DEM) is a raster representation and the most commonly used format for processed terrain data. It is a compact, highly structured and efficient representation. The accuracy depends on the resolution of the raster and the selected interpolation method. If the resolution is low compared to the variation of the data then the result will be inaccurate. If it is high, the data volume grows more than necessary. If there are large differences in the local variation of the data in different areas, then a trade-off between accuracy and data volume must be made. A raster represented surface is not completely defined by the sample points. Given a raster, values that are between the sample values have to be estimated. Several approaches are available. A simple method is bivariate evaluation where the estimated value is computed from a bivariate surface interpolating the four surrounding sample values. Inverse distance weighting (Maleika, 2020) can be used also in this context. We will illustrate the LR B-spline approximation by applying the algorithm to a selected data set. The bathymetry data set Figure  4(a) covers an area of slightly less than one square kilometer close to the Norwegian coast. It consists of about 11 millions points with a quite uniform density. The area is shallow, the depth varies from -27.94 m to -0.55 m. The point cloud contains outliers, the most significant one is two points with the same xand y-values and 2.38 m. difference in z. In Figure 4(b) the LR B-spline surface approximating the point cloud is trimmed to fit the point cloud domain. The surface has a rectangular domain in x and y. Since the point cloud has a boundary with large bends and contains holes, the surface will extrapolate the known information in some areas. To face that challenge, a 2-step procedure is applied:

Bathymetry data set Norway
• Restricting the values of the coefficients to an interval In general, both the LS and MBA methods make smooth transitions between existing point groups. Unfortunately, a combination of holes and data with a steep trend can lead to exaggerated heights, in particular with LS. The positivity of the basis functions combined with the partition of unity property implies that B-spline curves are bounded by their coefficients. In non-degenerate cases, this applies also to LR B-spline surfaces. Thus, we ensure well behaved surfaces in areas with missing points by restricting the values of the coefficients to an interval depending on the height range of the data points.
• Trimming We do not know the shape of the terrain in areas without points and use trimming to restrict the surface domain to the point cloud domain. To that aim, we bound the points by curves in theplane. The curves are arranged in one loop for the outer boundary and one curve for each hole and associated to the parameter domain (xy-plane for points parameterized by their x-and yvalues) of the surface. The outer loop is counter clockwise oriented, while eventual inner loops are clockwise oriented. By convention only the areas of the surface situated to the left of such trimming loops are considered valid. The sea floor covered by the data set contains both smooth areas and areas with significant variation. This can be recognized by the element structure of the surface shown in Figure 4(c). We can see that the surface is less refined in areas where there are no points and in areas with little variation in the sea bed. We approximate the point cloud with a LR B-spline surface using an increasing number of iterations. The threshold is set to 0.5 m and an initial tensor-product surface of ten by ten coefficients approximating the points is given. In this contribution, we select an approxition approach that has stable good results with respect to execution time and number of coefficients (Skytt & Dokken, 2021). The results of the approximation are presented in Table 1 and more details are given in Table 2. We compare the approximation accuracy of the LR B-spline surfaces with rasters of varying resolution. The raster surfaces are interpolated using inverse distance weighing with a valid distance of 10, i.e., each grid point is determined on the basis of approximately 8300 data points. The rasters are evaluated by linear interpolation between grid points. The results are presented in Table 1 and 2. The computation time excludes time for reading and writing. Table 1. Results of the approximation with LR B-spline and raster representation. R, is for Raster of a given accuracy in m, LRB is for LR B-spline surface, the number is the iteration step.
The maximum and average distances between the surface and the points are reported in Table 1 along with the number of LR Bspline coefficients or raster size. The number of points with a distance to the surface larger than the threshold and the file sizes are given. Execution time is reported for LR B-spline approximation, but not for computation of the raster surface. This implementation is less refined and the numbers would be misleading. Table 2. Point distribution with respect to the distance to the approximated surface (number of points for each range of distance in m).
-∞,-2.5 -2.5,-1.5 -1.5,-0.5 -0.5,0.5 0.5, The accuracy of the LR B-spline surface is increased with an increasing number of iterations and surface coefficients, but the gain is drastically decreased from the 16 th to the 20 th iteration. The raster surfaces have systematically poorer accuracy than the LR B-spline surfaces, see Table 2. A raster is a very uniform representation while the LR B-spline surfaces are able to adapt to the variation in the data sets. We further note that the LR B-spline surface after 7 iterations and the raster with resolution 1 m have similar file sizes. The raster surface is stored as GeoTiff, a binary format while the LR Bspline is stored in an ascii format. A compressed version of the LR B-spline file after 20 th iteration has size 744KB. Figure 5. (a) The distance field corresponding to the raster representation (1 m). (b) The distance field corresponding to the LR B-spline surface after 7 iterations. The same scaling was used for both representation Figure 5 shows the distance fields between the point cloud and the raster surface with resolution 1 m (a) and the LR B-spline surface after 7 iterations (b). In both cases the majority of the points lies within a short distance to the surface. The largest distances can be found in steep and ragged areas. Much less points have a large distance to the approximated surface even with few iterations to compute the LR B-spline surface, compared to the raster representation. This finding confirms the results presented in Table 1 and highlights the goodness of fit that can be reached with the LR B-spline surfaces.

LR B-spline volume for spatio-temporal analysis
In order to further demonstrate the potential of the LR B-spline to approximate GIS dataset from coastal region method, we made use of a dataset freely available at https://doi.pangaea.de/10.1594/PANGAEA.934058. The dataset was acquired for the CoastScan project (Vos and Kuschnerus, 2017) along the coast of Kijkduin in the Netherlands. The Terrestrial Laser Scanner Riegl VZ2000 laser scanner was used mounted on the roof of a hotel and programmed to perform a scan of the nearby dune and beach area every hour. More information about the pre-processing of the point clouds can be found in the dedicated publication. We selected the same rectangular area of 110 by 140 m from the point clouds for daily acquisitions and removed obvious outliers. The elevation range is [1.56, 7.73]. Within this area the pattern of the point clouds differs between the days as illustrated in Figure 6.   To get an overview of the elevation throughout the period, we add the time as a coordinate in addition to the x-and y-values of the points and interpret the elevation as an associated scalar. A point cloud assembled from all epochs is approximated by a LR B-spline volume. We apply a threshold of 0.2 m and 8 iteration steps. This is sufficient to get a very good approximation of the smooth component (see also Tab.3). The maximum distance between the points and the volume is 2.05 m and the average distance is 0.02 m. 1583 points out of 2.4 million have a larger distance than the threshold (0.066 %) and the average distance in these points is 0.41 m. Figure 7 shows the corresponding result. The elevation is visualized as colours: blue close to the seaside represents the minimum height. The elevation is increased as we move further from the sea and is relatively constant during time with the exception of the area marked with A in the figure. Here the elevation is increased for a period of a couple of days. At the position marked with B, there are two spots where the elevation is significantly higher than in the surroundings.  Figure  7 through one of the sand piles. Bottom: a horizontal cut in the period 5th to 6th of January.
We made some cuts in the volumetric approximation of the elevation during time to reveal information in the inner, see Figure 8. The direction from sea to land is now left to right while the time direction is vertical.
• We can recognize sand mounds at B. From Figure 8 (top), one is present during the entire period (B1) while the other appears mid-way (B2). • A height difference at A is identified together with a second less pronounced difference towards the end of the month at C. • Some objects seem to have appeared on the surface in a shorter period for later to be removed, e.g., at D, not visible in Figure 8 (top) but in Figure 8 (bottom).

Surface approximation to study deformation pattern
To reveal even more detail, we turn to the individual point clouds and approximate them with LR B-spline. In 16 cases the approximation converges totally with all points within the threshold. The approximation accuracy is reported in Table 3 for some data sets. In addition, the range of the reported numbers is presented. Distances are given in m. Their small values highlight the high confidence that can be put on the approximation. The number of coefficients varies with the data density as illustrated in Figure 6. The average distance is below 0.025 m. We note that the maximum distance is slightly higher for days 170105 and 170106. Fortunately, the average distance is similar to the other days: we link this discrepancy with a worthier approximation in a very specific domain. Indeed, we were able to identify local height variations related to some items on the surface. Figure 9 shows the positions of these items at January 6 th . The point cloud of this day is compared to the surface approximating the points of January 8 th when the items have been removed. Figure 9 shows the corresponding residuals. Table 3. Results of the surface approximation. The number of points, the maximum and average distance in m to the surface, the number of unresolved points and the number of surface coefficients are given.
To go further into details, we investigate specifically the aforementioned domains A, B, and C ( Fig.8) with the aim to analyse the deformation and find out the "story" beyond. From Figure 7 and 8 we know that there was a temporary increase in elevation in a part of the domain (indicated with the flag A). This happens between point cloud 170110_210054 and 170113_230038 as illustrated in Figure 10, which shows the distance fields between the point cloud from 170111 and the surface approximating day 170113. A surface to surface approximation comparison is shown in Figure 12 but not for other examples: it could indicate misleading differences in areas with missing information, i.e., with sparse and varying coverage. Kuschnerus et al. (2021) stated that on a small area of the beach at the bottom of a footpath, sand accumulated after a storm. Afterwards, the beach was cleared, and the superfluous sand forms a pile. This domain is indicated by B2 in Figure 8 and further analysed in Figure 11 where the points from 170117_030135 are compared with the surface approximating the cloud from 170116_010114. The pile marked B1 has been there during the entire period but has increased (see Fig.12). The pile B2 is new from 170117_030135. A surface-to-surface comparison for the same two point clouds is shown in Figure 12 for the sake of completeness. The point clouds are approximated by LR B-spline functional surfaces, 17016 and 17017 . The difference between them can be represented as a new surface, ( , ) = 170117 ( , ) − 170116 ( , ). In Figure 12, this surface is shown together with contour lines at intervals of 0.2 m. reaching from 0.0 m. to 1.6 m. In most of the area, the height of is less than 0.2 m. The zero level is show in black, the remaining contours are red. The elevation is higher at the 17 th , with some exceptions. The major differences are at the sand piles: B2 is new while B1 is extended.

Dealing with holes and ragged boundaries
The coverage of the selected point clouds varies, and several of the clouds have holes and/or ragged boundaries. The structures of the clouds related to the two surfaces 17016 and 17017 are relatively similar although 170117_030135 has one more hole than 170116_010114, and the size of the other hole is larger. To avoid comparing surface information without corresponding point cloud information, is trimmed with respect to 170117_030135. Finally, we investigate the feature marked with C in Figure 8. Figure 13 highlights that the difference is linked with a local alteration in sand height.

CONCLUSION
The LR B-spline surface allows adaptive surface refinement, and is a promising alternative to more standard NURBS representation for approximating GIS point clouds. In this contribution, we have set up the mathematical background to understand how local refinement is performed by means of Bsplines splitting and insertion of new meshlines in the LR-mesh. We have presented an effective method called MBA that allows to face the drawbacks associated with the LS method for approximating large point clouds. The LR B-spline representation has some advantages over the effective regularity of the raster format. We have highlighted how an effective surface approximation after only seven iterations could be performed with a much lower storage size. With respect to the raster approximation, the maximum distance was improved by 1 dm compared with the raster representation of resolution 0.5 m B1 A for our case study. Particularly in steep and ragged areas, the LR B-spline approximation performed better than the raster one. We have shown that increasing the number of iterations was not linked with an improved surface fitting. We have further shown that LR B-spline surfaces can be used to identify patterns of change in point clouds from TLS recorded at different days. A sand dune was taken as an example. We could reconstruct the spatio-temporal story of the dune before and after a sturm and highlighted specific height variations due to objects that were removed. We have shown that the proposed approximation method combining LS and MBA with LR Bspline could be used to approximate point clouds even in challenging domains with a variable data density and holes. LR B-spline surfaces are smooth and can, due to adaptivity, represent local details without a drastic increase in data size. Being a novel format, it is not supported in GIS systems. Fortunately, the surfaces can be exported as rasters in various resolutions as well as collections of tensor product spline surfaces. Outliers may obstruct the surface computation. Single points do not influence the surface significantly, but collection of outliers points and outliers in areas with sparse data points will drag the surface in their direction. The outlier problem is an important topic for further work using robust estimation method. Ongoing research focuses on the approximation threshold and a stop criterion for the number of iterations in the approximation algorithm using statistical quantity such as information criterion. The validation of the obtained surface is also an impotant topic with the intent to find a balance between the accuracy of the approximation and the number of coefficients representing the surface.