On the fast approximation of point clouds using Chebyshev polynomials

: Suppose a large and dense point cloud of an object with complex geometry is available that can be approximated by a smooth univariate function. In general, for such point clouds the “best” approximation using the method of least squares is usually hard or sometimes even impossible to compute. In most cases, however, a “near-best” approximation is just as good as the “best”, but usually much easier and faster to calculate. Therefore, a fast approach for the approximation of point clouds using Chebyshev polynomials is described, which is based on an interpolation in the Chebyshev points of the second kind. This allows to calculate the unknown coefficients of the polynomial by means of the Fast Fourier transform (FFT), which can be extremely efficient, especially for high-order polynomials. Thus, the focus of the presented approach is not on sparse point clouds or point clouds which can be approximated by functions with few parameters, but rather on large dense point clouds for whose approximation perhaps even millions of unknown coefficients have to be determined.


Introduction
Modern measuring instruments can record information about the geometry of physical objects and provide the user with discrete points in a Cartesian coordinate system.In geodetic practice, such measurements are mainly carried out with terrestrial laser scanners or photogrammet-ric methods, which are capable of capturing millions of points for the observed object.Often an appropriate curve or surface approximation of the measured object with a continuous function is required, either for further analysis of the measurement results (e. g. deformation analysis) or for the representation of the data in a CAD (computeraided design) model.For objects with complex shapes, freeform curves or surfaces are mainly used for data modelling, usually represented by a linear combination of basis functions, which can be freely chosen by the user depending on the problem at hand.Some sets of basis functions that are widely used and can be easily found in the literature are, for example, B-Splines by Wang [43] or Harmening [21], Chebyshev polynomials by Chang [11], Clenshaw [13] or Smitka [40], Fourier series by Wang [42], radial basis functions by Majdisova [26] or Ohtake [33], spherical harmonics by Gannon [18,30] or Mousa [30], wavelets by Manson [27] or Zernike polynomials by Ibanez [24].The variety of different sets of basis functions in combination with the different methods available for calculating the unknown parameters, provides the user with an almost infinite number of possibilities for an approximation of point clouds using freeform curves or surfaces.
In geodetic science, B-Splines are mainly used as basis functions for the approximation of point clouds, given by the recurrence relation introduced by de Boor [15] and Cox [14], while the unknown coefficients of the resulting Spline are usually estimated using the method of least squares.An overview on the approximation of 3D point clouds with freeform curves and surfaces, including Bézier, B-Spline and NURBS is given in the work of Bureick et al. [10] and different mathematical representations of 2D splines and their least squares solution for data approximation are discussed by Ezhov et al. [16].A general review of point cloud modelling methods for deformation analysis has been presented by Neuner et al. [31] and the quantification of object deformation using B-spline approximation is thoroughly discussed by Harmening [20].
An aspect that has a significant influence on the result of the spline approximation is the optimal knot placement, while different strategies exist.The location can be predefined as proposed in the study of de Boor [15] or Yanagihara and Ohtaki [46], or treated as unknown parameters leading to a nonlinear optimisation problem, which can be solved either directly as a non-linear adjustment problem, as e. g. shown by Schwetlick and Schütze [39], or by using an additional optimisation algorithm, as e. g. addressed by Galvez et al. [17], Bureick et al. [8,9] or Schmitt and Neuner [37].Additionally, several criteria have been presented in the literature for choosing the optimal number of control points, which is known as model selection.In this regard, there are many studies that use the Akaike information criterion [2,3] or the Bayesian information criterion [38].A detailed overview of model selection strategies is given by Harmening [20, p. 25 ff.].Both, the problems of knot placement and model selection can be considered as separate research topics and have a significant impact on the quality of the approximation using B-Splines.
Motivated by these challenges and since these approaches are usually used for small point clouds with only a few unknown parameters to be determined, we present a strategy for the approximation of point clouds using Chebyshev polynomials, which is based on an interpolation in the Chebyshev points of the second kind.Thus, we deliberately avoid the computation of the "best" approximation using the method of least squares in favour of a "near-best" solution, which can be computed much easier and extremely fast using the Fast Fourier transform (FFT).The main focus in this paper is on the curve approximation of large and dense point clouds of objects with complex geometry, consisting of hundreds of millions of points with maybe even millions of unknown parameters to be determined.
Although the basic idea of our presented approach is not new and has been already applied, see e. g. the contributions by Clenshaw [13], Chang [11] or McKinley and Ishihara [29], it has unfortunately received little attention in the current literature.Recent contributions in geodetic science that are known to the authors in this regard are those of Hu et al. [23] and Xi et al. [44], without however exploiting or at least addressing the full potential of this approach.Since we approximate the whole point cloud by a single polynomial, a split of the entire domain, as e. g. required for an approximation with B-splines by choosing an appropriate knot vector, is not needed.In addition, the determination of the optimal polynomial degree is solely based on the unknown coefficients and not, as usually, on the residuals.Therefore, an additional model selection criterion is also not needed.
In this contribution we consider at first only point clouds that represent curves without any kinks, jumps or data gaps, which consequently can be approximated by continuous smooth functions.Since for the determination of the unknown parameters, as well as the optimal polynomial degree, the residuals no longer play an important role, a requirement for orthogonal residuals is not considered or rather explicitly omitted in this paper.
The paper is organized as follows: In Section 2 we present the least squares approximation of functions by a polynomial in the Chebyshev basis and its numerical solution, as it as known in the literature.Afterwards, we adapt the basic methodology of this approach and show in Section 3 how it can be used to efficiently approximate point clouds, which is the main contribution of our work.An example is given in Section 4, while Section 5 briefly summarises the investigations that were carried out and critically reviews the results.An outlook on current research regarding the presented approach concludes this paper.

Approximation of functions
At first, we present the approximation of functions by a polynomial in the Chebyshev basis using the method of least squares.Since usually the calculation of the coefficients of the polynomial is very time-consuming and inaccurate, we also present how the coefficients can be calculated fast and efficiently by an interpolation in the Chebyshev points using the Fast Fourier transform (FFT), as it is known in the literature.

Least squares approximation
The Chebyshev polynomial of the first kind T j (x) of degree j ≥ 0 is defined by se e. g. Mason [28,p. 2] or Rivlin [36, p. 2].These polynomials satisfy the three term recurrence relation with the initial conditions Let us consider a function f (x) ∈ L 2 ([−1, 1]), the Hilbert space of square integrable functions, which we want to approximate by a polynomial P p (x) of degree p in the Chebyshev basis.P p (x) lies in the vector space V spanned by the p + 1 basis functions and is given by the linear combination with the unknown coefficients c j .In general, f (x) is not within V and is therefore different from P p (x).The difference between f (x) and P(x) is known as residual function v(x) and can be expressed by To determine the unknown coefficients c j of a least squares approximation of f (x), we minimise the length of v(x), given by the with ⟨⋅, ⋅⟩ being the inner product of two functions, see e. g.Bronshtein et al. [7, p. 920].Thus, the objective function reads with Minimising the length of v(x) is equivalent to demand that v(x) is orthogonal to vector space V, or as Langtangen [25, p. 10] has shown, that v(x) is orthogonal to the p + 1 basis functions.This leads to the inner products with i = 0, 1, . . ., p. Rearranging (6) and inserting into (10) yields and introducing (5) leads to the system of normal equations Since the inner product is bilinear, this equation system results in for i = 0, 1, . . ., p, which can be written in matrix notation as follows As Chebyshev polynomials are orthogonal with respect to a weighted inner product the normal equation system ( 14) simplifies to The normal matrix N is diagonal and the unknown coefficients can be directly determined by with for j = 0, 1, . . ., p.
The derivation of an analytical solution of the integrals in (17) for arbitrary functions f (x) is hardly feasible, so that a numerical solution is usually used.But, due to singularities of the weighting function in ( 18) at x = ±1, the numerical integration of the right hand side is slow and inaccurate for higher polynomial degrees.To derive a very fast and accurate solution for the numerical integration of (18), we present in the following the relationship between the Chebyshev and Fourier series.

Connection to the Fourier series
Let us introduce the following change of variables and with we obtain dx = − sin(θ) dθ.
Inserting ( 19) and ( 21) in the right hand side of ( 18) and changing the limits according to cos(π) = −1 and cos(0) = 1 yields Hence, the unknown coefficients of a least squares approximation of f (x) in the Chebyshev basis can be determined by The presented change of variables can also be interpreted as a mapping of f (x) onto the upper half of the unit circle, as depicted in Figure 1.An arbitrary point z on the unit cir- cle is uniquely defined by its coordinates (x, y), satisfying or by an angle θ, with Re(z) = x = cos θ, while the upper equation in (25) is the introduced change of variables (19).Due to it follows and the integral in (23) can also be written as which corresponds to a scaled Fourier cosine transform.While mapping f (x) onto the unit circle we obtain a transplanted function F(θ) for θ ∈ [−π, π] which is an infinitely differentiable, even and periodic function, as discussed by Orszag [34].Therefore, F(θ) has a Fourier cosine series expansion where the coefficients a j are given by and which are equivalent to the coefficients c j in equation (28) except for a scaling factor, see the handbook by Bronshtein et al. [7, p. 927].Based on this fact, a fast and accurate numerical solution for the integrals in ( 28) is presented in the next section and the resulting coefficients of a polynomial approximation in the Chebyshev basis are discussed in more detail on an example.

Numerical solution by an interpolation in the Chebyshev points
In case F(θ) in ( 30) is known at N + 1 equally spaced points θ k , with k = 0, 1, 2, . . ., N, the integral in ( 30) can be determined approximately by for j = 0, 1, . . ., N and, hence, can be efficiently computed by the Fast Fourier Transform (FFT) or Fast Cosine Transform, see the handbook by Bronshtein et al. [7, p. 928].According to Trefethen [41, p. 15] this has been probably observed around 1970, see for example the investigations of Ahmed [1] or Orszag [34].
Let us consider some equally spaced points on the upper half of the unit circle, including the boundary points on the intersection with the real axis at ±1, as depicted in Figure 2. The projection of these equally spaced points onto the real axis are known as Chebyshev points of the second kind, Chebyshev extreme points or Chebyshev-Lobatto points and are given by In the following we just call them Chebyshev points.These points are the extrema of the Chebyshev polynomials of the first kind, in contrast to the Chebyshev points of the first kind, which are the roots.How the Chebyshev points of the first kind can be used for the approximation of functions will not be considered here any further, but is discussed e. g. by Xu [45].
If we discretize f (x) at the Chebyshev points we obtain which can directly be used to calculate the unknown coefficients c j using FFT.Thereby the solution of the integrals in equation ( 18) are only approximated and the calculated coefficients c j correspond to those of an interpolation in the Chebyshev points.
Although in theory we approximate functions by the method of least squares, this approach is unsuitable in practice as the integrals in equation ( 17) can be in general solved only by numerical integration, which is in this case slow and inaccurate.It is rather more efficient to approximate functions by interpolation in the Chebyshev points, where the coefficients calculated by FFT converge to the values of the exact solutions for equation (17) as the number of points increases.
To demonstrate such an approximation, we consider the following function for x ∈ [−1, 1], which is depicted in Figure 3.
To show a very important property of an approximation in the Chebyshev basis, we approximate f (x) by a polynomial of degree p = 200.The coefficients c j are depicted in Figure 4 and decrease geometrically until a plateau of rounding errors near level of machine precision is reached for j ≥ 117 and a floating-point relative accuracy of ϵ = 2.22 ⋅ 10 −16 (vertical red line).As long as f (x) is continuous and satisfies the Dini-Lipschitz condition, the Chebyshev series expansion will converge to f (x), see the textbook by Mason and Handscomb [28,Chapter 5.3.2].As a rule of thumb, the truncation error of this approximation is in the same order-of-magnitude as the absolute value of the last coefficient, as mentioned by Boyd [6,p. 51].Aurentz [4] designed a chopping algorithm that detects and chops the series at the beginning of this plateau.Consequently, for this example a polynomial of degree p = 116 is found as the best approximation of f (x).The resulting residual function v(x) is shown in Figure 5 and reveals, as expected, differences near machine precision.As illustrated by this  17) can be calculated very fast and accurately by an interpolation in the Chebyshev points using FFT.Furthermore, the optimal polynomial degree of the approximation can be identified by an iterative increase until the coefficients reach the plateau of rounding errors near level of machine precision, which can then be detected and truncated.This provides an efficient and accurate approach for the approximation of functions and is the core of the open-source package Chebfun which was in 2004 by Battles and Trefethen [5].

Approximation of point clouds
Since we can consider functions as a special case of a point cloud, i. e. a point cloud with an infinite number of errorfree points, we can directly adopt the presented procedure for the approximation of point clouds.At first, let us generate a point cloud by discretising function (34) at n + 1 linearly spaced points (x i , y i ) for i = 0, 1, . . ., n.Now, the actual problem concerning the discretised function arises from the fact, that there are usually no more functional values f (x k ) at the Chebyshev points x k available.Therefore, we have to extract or interpolate the corresponding values for the Chebyshev points from the nearest surrounding points or a cluster of neighboured points.
As long as the true values for f (x k ) at the Chebyshev points are not given any more, the coefficients will not decrease until they reach a plateau of rounding errors near level of machine precision.In this case, a plateau will be formed at a level ≫ ϵ that depends essentially on the density of the point cloud and the amount of noise.Consequently, the task is to identify this plateau and its level and to truncate the coefficients at its beginning using the chopping algorithm of Aurentz and Trefethen [4].These are in principle the main differences to the approximation of functions using Chebfun [12], whereas the functional values f (x k ) and the level of the plateau are directly given.Hence, Algorithm 1 describes a simple and fast approach for the approximation of point clouds using Chebyshev polynomials.
Algorithm 1 Pseudo code for the calculation of the coefficients c j for a polynomial approximation of the point cloud x i , y i .
⊳ Initialisation for p. return c j 17: end function Since we approximate the point cloud by an interpolation in the Chebyshev points x k , the accuracy of the approximation depends directly on the accuracy of their functional values f (x k ).Besides the shape of the underlying unknown function, the density of the point cloud and the amount of noise mainly influences the accuracy of f (x k ) and consequently of the approximation itself.Both aspects will be addressed in more detail in the following sections.

Impact of point density
At first, we would like to demonstrate the impact of the point density on the level of the plateau of the coefficients and, thus, on the accuracy of the approximation itself.Therefore, we generate point clouds using function (34) at n + 1 = 10 2 10 4 10 6 linearly spaced points and without taking any measurement noise into account.Each point cloud will be approximated by a polynomial of degree p = 200, while the functional values f (x k ) at the Chebyshev points x k will be determined by a simple linear interpolation.The resulting coefficients for each point cloud are depicted in Figure 6.As already mentioned, and also illustrated by Figure 6, the coefficients decrease until they reach a plateau at a certain level, which is directly related to the point density.This is essential for the presented approach, because the higher the level of the plateau the lower the resulting polynomial degree and consequently, the less accurate the approximation will be.
For this example, the mean level of the plateau for 10 2 (blue), 10 4 (red) and 10 6 (yellow) points is located at approximately 10 −6 , 10 −10 and 10 −14 , respectively.It can be observed that an increase in the number of points by a factor of 10 2 reduces the level of the plateau by ≈ 10 4 .For a point cloud with 10 6 points, the resulting coefficients differ to the true values given in Figure 4 by only about ±5 ⋅ 10 −13 .It should also be noted, that the values for the coefficients forming the plateau do not have the same random characteristics as those in Figure 4, but rather show some systematic variations, which probably result from the calculation of the functional values f (x k ) by a linear interpolation.In general we can conclude that the accuracy of the approximation of point clouds, using the proposed algorithm, increases proportionally with the point density.

Impact of noise
In this section we illustrate the impact of measurement noise on the level of the plateau.Therefore, we generate a point cloud with 10 4 linearly spaced points and add normally distributed noise with a standard deviation of σ = 10 % 1 % 0.1 % .The resulting coefficients for each point cloud are depicted in Figure 7. First of all, it is apparent in Figure 7 that the three plateaus are quite close to each other and their levels only differ by a factor of 10 and, therefore, share the same difference as the corresponding standard deviation σ to which they refer.Furthermore, the values of the coefficients forming each plateau have, in this example, a rather random characteristic.In a direct comparison with the coefficients of the point cloud without noise (purple), it is also noticeable that noise has a considerable impact on the level of the plateau.Even if only low noise with σ = 0.1 % is present, the level of the plateau already increases from 10 −10 to ≈ 10 −4 .A part of this increase is certainly due to the choice of calculating the functional values f (x k ) by a simple linear interpolation.Thus, the result could be significantly improved by using better methods than a linear interpolation.However, the measurement noise will always have a considerable impact on the level of the plateau and consequently on the accuracy of the approximation.
As we have shown, the accuracy of an approximation of point clouds by an interpolation in the Chebyshev points is directly related to the point density and measurement noise.In general we can state, that especially for sparse point clouds or point clouds of simple geometries, which can be approximated by low-order polynomials, the presented approach is not the best choice.In such cases it is preferable to estimate the polynomial coefficients us-ing the method of least squares and determine the optimal polynomial degree iteratively using a suitable model selection criterion, as for instance discussed by Harmening and Neuner [22].
However, the presented approach is particularly suitable for the approximation of dense point clouds of complex geometries, where even polynomials of millionth degree can be usually calculated within a few seconds.So the focus is on point clouds where the number of points is much larger than the degree of the polynomial approximation.This will be illustrated in the following section.

Example
To demonstrate the suitability of the presented approach we consider the following example of a profile scan of a bust of C. F. Gauss with 100 000 points, as depicted in Figure 8.The profile scan was extracted from a surface scan of the bust and the number of points was increased afterwards.In addition, the resulting point cloud was smoothed in order to reduce the impact of measurement noise as much as possible.The derived profile scan is considered as a point cloud without measurement noise and is used as a reference for the approximation in this example.
To analyse the accuracy of the approximation using our approach, we generate a synthetic point cloud by adding normally distributed noise with a standard deviation of σ x = σ y = 0.5 mm to the true coordinates of the points of the reference profile scan.To approximate such a point cloud, we consider the parametric representation of the x and y coordinates in the form with the new parameter u ∈ [−1, 1].In general, there is no corresponding value u i for each point (x i , y i ) available.Therefore, three different methods can be found in the literature for calculating u i for each point, namely "equally spaced", "chord length" and "centripetal", see e. g. the textbook by Piegl and Tiller [35, p. 364 ff.].In this example, we choose equally spaced values for u i .Since in this case we determine the functional values at the Chebyshev points by a simple linear interpolation, the stochastic properties of the observations (x i , y i ) will have no influence on the determination of these values.Therefore, we can assume equally weighted and uncorrelated observations (x i , y i ) and calculate the unknown parameters c x j and c y j independently.This results in the two point clouds with measurement noise depicted in Figure 9, which will be approximated in the following sections using both our proposed approach and the method of least squares.As already mentioned, we have decided, for reasons of simplicity, to determine the functional values at the Chebyshev points by a linear interpolation.To obtain more precise results and also to take the stochastic properties of the observations (x i , y i ) into account, the function values at the Chebyshev points could be estimated by e. g. local regression.However, this is out of the scope of this contribution.

Approximation by interpolation
At first, we approximate the two point clouds from Figure 9 by an interpolation in the Chebyshev points according to the proposed approach described in Algorithm 1.As a solution we obtain a curve which is described by two polynomials of degree p x = 528, p y = 422 and which can not be distinguished visually from the two point clouds in Figure 9. Therefore, we first of all will have a look at the resulting residuals v x i = P x,528 (u i ) − x i and v y i = P y,422 (u i ) − y i , which in this case are shown as a relative frequency histogram in Figure 10.The two histograms are just slightly different with maximum values of max |v x i | ≈ 2.07 mm and max |v y i | ≈ 2.26 mm.The distribution of both sets of residuals correspond very well to a normal distribution with an expectation nearly zero and a standard deviation of s x = 0.5021 mm and s y = 0.5022 mm for v x i and v y i , respectively.The residuals of both approximations therefore correspond to the previously added normally distributed noise with a standard deviation of σ x = σ y = 0.5 mm.
To analyse the accuracy of the approximation in more detail, the deviations Δ x i and Δ y i , which result as the difference between the polynomials and the point cloud of the reference profile scan without measurement noise from Figure 8, are shown in Figure 11.Due to the calculation Figure 11: Deviations Δ x i and Δ y i between the resulting polynomials and the point cloud of the reference profile scan without measurement noise from Figure 8.For a better visualisation of the deviations, the value range was limited to ±0.3 mm.Therefore, the largest deviations at the limits at u ≈ ±1 are not visible.
of the functional values f (x k ) at the Chebyshev points by a simple linear interpolation, the largest deviations occur very close to the limits at u ≈ ±1, with maximum values of max |Δ x i | ≈ 0.59 mm and max |Δ y i | ≈ 1.21 mm.However, for larger polynomial degrees, these areas are negligibly small in relation to the total point cloud and can even be removed afterwards if desired.Alternatively, these deviations at the limits can also be reduced by choosing more precise methods for determining the functional values f (x k ).In addition, larger deviations occur in the range of u = −0.1 in both point clouds, which corresponds to a spot on the forehead of the bust.There is a small dent which could not be approximated properly, due to the large measurement noise.Nevertheless, these deviations are still much smaller than the introduced measurement noise of σ x = σ y = 0.5 mm.Apart from some few exceptions, the deviations for the x and y coordinates are mainly below 0.1 mm, with a resulting standard deviation of s x ≈ 0.049 mm and s Δ y ≈ 0.051 mm for Δ x i Δ y i , respectively.

Least squares approximation
In this section we will derive the basic formulas to compute the solution for a polynomial approximation using the method of least squares, see e. g. the textbooks by Ghilani and Wolf [19, p. 173 ff.] or Niemeier [32, p. 129 ff.].Since the derivation of the solution for a least squares approximation is identical for x i and y i , we will present it only for x i .The functional model for x i for the n + 1 points according to equation (35) results in with the unknown parameters c x j and for i = 0, 1, . . ., n.If we consider x i as observations and u i as error-free values and introduce a residual v x i for each observation, we obtain the following observation equations which can be written in matrix notation as x 0 x 1 . . .
) Assuming equally weighted and uncorrelated observations, the weight matrix P results in an identity matrix which represents the stochastic model.To determine the unknown parameters c x j of a least squares approximation, we minimise the sum of weighted squared residuals, resulting in the so-called normal equation system with the normal matrix and the vector of the right hand side Solving (40) for c x yields the unknown coefficients c x j of a least squares approximation.In contrast to the resulting normal matrix for the least squares approximation of functions in equation ( 16), the normal matrix in ( 41) is a dense matrix.Therefore, the normal equation system (40) can not solved as efficiently as (16).
For the two point clouds in Figure 9 we determine the polynomial approximation by solving (40) with a polynomial degree of p x = 528 and p y = 422.In principle, the same distribution of the residuals results for both polynomials with almost the same properties as already obtained for the approximation by interpolation from Figure 10.Also the resulting deviations hardly differ visually from those in Figure 11, with just slightly smaller values at the limits of max |Δ x i | ≈ 0.44 mm and max |Δ y i | ≈ 0.97 mm.
Accordingly, the polynomial approximation by the method of least squares seems to correspond to the one by an interpolation in the Chebyshev points from the previous section, which will be analysed in more detail in the following.

Comparison
At first, we will compare the solution of the polynomial coefficients obtained by an interpolation in the Chebyshev points from Section 4.1 and by the method of least squares from Section 4.2.The absolute differences between these two sets of coefficients for the polynomials P x,528 (u) and P y,422 (u) are shown in Figure 12.In general, small differences have been observed for most of the coefficients, that range below 10 −3 mm, whereby the values for P y,422 (u) (red) are a bit smaller than for P x,528 (u) (blue).Moreover, both graphs show a similar trend, with an exponential increase of the difference for larger j.Therefore, the maximum absolute difference of ≈ 7 ⋅ 10 −3 mm occurs at the last coefficient for both polynomials, which is still nearly 100 times smaller than the added measurement noise of Figure 12: Logarithmic plot of the absolute difference between the polynomial coefficients obtained by an interpolation in the Chebyshev points and by the method of least squares for P x,528 (u) in blue and P y,422 (u) in red.σ x = σ y = 0.5 mm.Consequently, the two sets of polynomials only differ slightly from each other, as can be seen in Figure 13.
As in the previous sections, the largest deviations of up to 0.25 mm occur at the limits of the two polynomials.However, since both polynomials approximate the true function rather poorly at the limits, compare Figure 11, no conclusion can be drawn at this point about which approximation is better or worse.Additionally, ∼ 99.95 % of all deviations are smaller than 5 ⋅ 10 −2 mm and even ∼ 85 % are smaller than 5 ⋅ 10 −3 mm, thus, they are 100 times smaller than the added measurement noise.Whether these very small deviations play a significant role depends, of course, on the problem itself and has to be answered individually.In the present example, however, these deviations are not relevant.
Figure 13: Difference between the polynomials obtained by an interpolation in the Chebyshev points and by the method of least squares.For a better visualisation of the differences, the value range was limited to ±0.015 mm.Therefore, the largest differences at the limits at u ≈ ±1 are not visible.
Except for the limits, a polynomial approximation of dense point clouds by interpolation in the Chebyshev points or by the method of least squares yields almost the same result.But with the important difference, that the approximation by interpolation can be calculated much faster and also for polynomials of millionth degree, which is hardly possible for an approximation using the method of least squares according to Section 4.2.To clearly emphasise this fact, Table 1 shows the computation time required for the approximation of different sized point clouds and polynomial degrees.For small point clouds and low polynomial degrees, there is hardly any significant difference in the computation time between these two approaches, but this changes quickly with increasing number of points and polynomial degrees.For case 3, the computation time using the method of least squares is already ≈ 10 min and, thus, nearly 150 000 times longer than for interpolation.Most of the time is thereby needed to set up the normal matrix (41) and, based on extrapolations, this takes already about 2 weeks for case 4.This is also the reason why in Table 1 the computation time for an approximation using the method of least squares for the last 3 cases could no longer be given.In contrast, even a point cloud with 100 million points can be approximated by a polynomial of 10 millionth degree in less than 2 s using the presented approach of an interpolation in the Chebyshev points.

Conclusion and outlook
By slightly modifying the general methodology for an approximation of functions by an interpolation in the Chebyshev points, as applied by Chebfun [12], the basic idea of an approach was shown that allows even large point clouds with complex geometry to be approximated very fast.The identification of the optimal polynomial degree is based on the determination and truncation of the plateau and is therefore solely based on the computed polynomial coefficients.The impact of the point density and measurement noise on the level of the plateau and consequently on the accuracy of the approximation has been shown in more detail.It has been demonstrated that the measurement noise has the largest impact.The main advantages of the presented approach for the univariate approximation of point clouds by interpolation in the Chebyshev points in contrast to an approximation using B-Splines can be summarised as follows: -No knot placement strategy needed.
-No additional model selection criteria needed.
-Extremely fast and efficient.
So far, the main challenge of the presented approach is the determination of the functional values f (x k ) at the Chebyshev points x k and also the determination of the level of the plateau.Both aspects have a significant impact on the approximation itself.Since we wanted to demonstrate the basic concept of an approximation in the Chebyshev basis by an interpolation in the Chebyshev points in this contribution, we have limited ourselves to the simplest approaches for both aspects.How both aspects can be solved in a smart way is currently being researched.Moreover, we have neglected the stochastics of the unknown parameters so far, which is, among other aspects, one of the essential result of a least squares approximation.In case the approximation is used e. g. for a geodetic deformation analysis, the stochastic properties of the unknown parameters play an essential role and must be taken into account.Therefore, the determination of the stochastic properties as well as the identification of outliers is currently being investigated in more detail.
Until now, only point clouds that can be described by univariate functions, for example, a trajectory, can be approximated by the presented approach.However, the basic methodology presented in this article, could be directly extended to the approximation of surfaces and is currently being developed.

Figure 1 :
Figure 1: Unit circle in the complex plane with the imaginary unit = −1.

Figure 2 :
Figure 2: Equispaced points on the upper half on the unit circle in blue and Chebyshev points in red.

Figure 4 :
Figure 4: Logarithmic plot of the coefficients |c j |.
Furthermore, we only consider point clouds on the interval [−1, 1] since any interval [a, b] can be easily scaled to [−1, 1].

Figure 8 :
Figure 8: Reference profile scan of a bust of C. F. Gauss with 100 000 points.

Figure 9 :
Figure 9: Resulting two point clouds for the x and y coordinates with measurement noise of σ x = σ y = 0.5 mm and 100 000 points.

Figure 10 :
Figure 10: Relative frequency histogram of the resulting residuals v x i and v y i for the approximation of the two point clouds with measurement noise of σ x = σ y = 0.5 mm and 100 000 points.

Table 1 :
Computation time required for an approximation by interpolation in the Chebyshev points (IP) and by the method of least squares (LS) of different sized point clouds and polynomial degrees using a budget PC with an Intel Core i5-6400 @ 2.70 GHz.All computation have been performed in Matlab and the normal equation system (40) has been solved using the backslash operator.