On the Constrained Solution of RBF Surface Approximation [original]

Citation: Pasioti, A. On the

Constrained Solution of RBF Surface

Approximation. Mathematics 2022,10,

2582. https://doi.org/10.3390/

math10152582

Academic Editors: Janusz Brzd˛ek

and Alexander Felshtyn

Received: 14 June 2022

Accepted: 20 July 2022

Published: 25 July 2022

Publisher’s Note: MDPI stays neutral

with regard to jurisdictional claims in

published maps and institutional affil-

iations.

Licensee MDPI, Basel, Switzerland.

This article is an open access article

distributed under the terms and

conditions of the Creative Commons

Attribution (CC BY) license (https://

creativecommons.org/licenses/by/

4.0/).

mathematics

Article

On the Constrained Solution of RBF Surface Approximation

Anastasia Pasioti

Institute of Geodesy and Geoinformation Science, Technische Universität Berlin, Kaiserin–Augusta–Allee 104–106,

10553 Berlin, Germany; [email protected]

Abstract:

In this contribution, a scattered data approximation problem, chosen from the literature

and using the Radial Basis Function (RBF) approach, is considered for the application of point cloud

modelling. Three solutions are investigated for the approximation problem. First, a technique

known from the literature is investigated using a linear combination of thin-plate splines and linear

polynomials, with additional constraint equations. Then, using the same approximation function as

before, a technique is developed for a rigorous consideration of the constraint equations. Finally, a

technique is presented in which the approximation function consists only of a linear combination

of thin-plate splines, without the introduction of linear polynomials and constraint equations. In

addition, some interpolation problems with the RBF approach are discussed to present the differences

between an interpolation with thin-plate splines only and an interpolation with thin-plate splines

together with linear polynomials and constraint equations. Numerical examples are given to illustrate

and discuss the solutions from the different techniques.

Keywords:

Radial Basis Function (RBF); thin-plate spline; surface approximation; surface interpolation;

point cloud

MSC: 65D12

1. Introduction

Nowadays, the use of 3D laser scanning technology is widespread in various appli-

cations for digitizing real world objects. The acquired point cloud consists of measured

coordinates in

IR3

of numerous points. The number of points can be up to several millions,

resulting in a huge data set. Moreover, the points are considered unorganized and scat-

tered, which means that there is no specific structure or topological relation between them.

Hence, the points and their coordinates are not sufficient to provide additional informa-

tion about the object. As Bureick et al.

[1]

mentioned, only by modelling the geometric

information can further analysis be performed, and several questions which could not be

addressed before can be answered. For instance, in Computer-Aided Design (CAD) or in

deformation analysis, a surface approximation with a continuous mathematical function is

required, as, e.g., pointed out by Bureick et al.

[1]

. Therefore, modelling the point cloud is a

necessity and surface approximation is taken into account.

There are numerous approaches to approximate a surface involving functions such

as polynomials, splines, B–splines, etc. For instance, Bureick et al.

[1]

presented surface

approximations based on polynomials, Bézier, B-splines and Non-Uniform Rational B-

splines (NURBS). Each of these approaches has its advantages and some disadvantages. In

this article, the Radial Basis Function (RBF) approach is discussed. This choice is motivated

by typical tasks in Geodesy and Geoinformation Science. There, the challenge is often to

model point clouds that include a very large number of points in the 3D space dimension

with unorganized points (scattered points), which point clouds have been captured with

laser scanners or photogrammetric methods. As Buhmann

[2]

(p. 1ff.) explained, when the

functions to be approximated (a) depend on many variables or parameters, (b) are defined

Mathematics 2022,10, 2582. https://doi.org/10.3390/math10152582 https://www.mdpi.com/journal/mathematics

Mathematics 2022,10, 2582 2 of 26

by possibly very many data, and (c) the data are scattered in their domain, then the Radial

Basis Function approach is especially well-suited.

Another advantage of the Radial Basis Function approach is that it is a meshfree

approach. In contrast to some other approaches that require a mesh such as those using

wavelets, multivariate splines, finite elements, etc., as Wendland

[3]

(p. ix) explained.

According to Fasshauer

[4]

(p. 1), mesh generation can be the most time-consuming part

of any mesh-based numerical approach. Moreover, he continued, a meshfree approach is

often better suited to cope with changes in the geometry of the domain of interest (e.g., free

surfaces and large deformations) than approaches such as finite elements.

Mainly, in this article, the approximation of point clouds is considered, since with high

point density and the fact that measured coordinates are subject to random measurement

errors, interpolation is no longer appropriate. In the case of interpolation, measurement er-

rors and other abrupt changes in the data points would be modelled, resulting in a strongly

oscillating modelled surface. Nevertheless, some numerical examples of interpolation, with

a few data points, will be presented to reveal the effect of the applied technique.

For the scattered data approximation problem with the RBF approach, an approxima-

tion function is built that consists of a linear combination of radial basis functions. Each

radial “basis” function

Φ(||·||2)

is centered at each fixed center point, and by composi-

tion with the Euclidean norm

||·||2

is radially (or spherically) symmetric about this

center point. A formal definition of a radial function is given by Fasshauer

[4]

(p. 17).

Furthermore, each radial “basis” function is a multivariate function which is based on a

univariate continuous “basic” function

φ(r)

. The basic function is used to derive all the

radial basis functions in the approximation function. The terminology “basis” and “basic”

function is adopted from the textbook by Fasshauer

[4]

(p. 18). Based on the contribution by

Flyer et al.

[5]

, Table 1shows some of the well-known basic functions, where

is a shape

parameter defined by the user.

Table 1. Basic functions.

Name φ(r)

Gaussian e−(er)2

Multiquadric p1+ (er)2

Inverse multiquadric 1

√1+(er)2

Thin-plate spline r2log(r)

In this article, the basic function used, namely thin-plate spline, is chosen for the impor-

tant advantage that no shape parameter is required. Based on the analysis of

Fasshauer [4]

(p. 23ff.), the value of the shape parameter has a great influence on the numerical sta-

bility of the interpolation problem, and if not chosen appropriately, can result in severe

ill-conditioning of the interpolation matrix. A sophisticated method is, e.g., presented in the

paper by Zheng et al.

[6]

, who used a test differential equation for the selection of a good

shape parameter. In the following, a shape parameter-free basic function is preferred. In ad-

dition to the previous advantage, the thin-plate spline function, according to

Fasshauer [4]

(p. 170), has the tendency to produce “visually pleasing” smooth and tight surfaces.

The first appearance of the thin-plate spline function in an interpolation problem was

presented by Harder and Desmarais

[7]

under the name “surface spline” for applications

related to aircraft design. Duchon

[8]

studied their method and extended the mathematical

representation based on Hilbert kernel theory, which had a major impact on subsequent

research on the interpolation of spline functions with radial basis functions. In his article,

Duchon

[8]

was interested in an interpolating surface with minimal bending energy, which

geometrically means that he was aiming for an interpolation that was as flat as possible, that

is, as close to a plane as possible. He composed, for this purpose, the interpolating function

as a linear combination of thin-plate spline functions with an addition of linear polynomials.

Furthermore, condition equations were taken into account to ensure a unique solution.

Mathematics 2022,10, 2582 3 of 26

As Harder and Desmarais [7] explained, the thin-plate spline function has a physical

interpretation which is a plate of infinite extend that deforms in bending only. Since the

bending energy is represented by the integral of the partial derivatives, the minimal bend-

ing energy of the surface that Duchon

[8]

was interested in is approximately equivalent to

the minimal integral curvature, and hence with the minimum of the functional of second

derivatives. For this minimisation, Duchon

[8]

defined a semi-norm in the infinitely dimen-

sional vector subspace of functions of finite energy. Furthermore, since linear polynomials

defining a plane are equal to zero in their second derivatives, they are also included in the

interpolation function. Finally, Duchon

[8]

proved that the functions, thin-plate splines

and polynomials can be used as a basis for surface interpolation, and the space of this

basis is a subspace of the infinite dimensional Hilbert vector space of all real functions

of two real variables. Later on, this form of interpolation presented by Duchon

[8]

was

widely used in the literature, for instance, by Buhmann

[9]

, by Wendland

[3]

(p. 116ff.),

Fasshauer [4]

(p. 70ff.), and by Flyer et al.

[5]

, to name a few, when interpolating with

thin-plate spline functions.

Here, it should be clarified that this form of interpolation that is adding a series of

polynomials to the interpolation function and also considering condition equations is an

interpolation technique that was first introduced by Hardy

[10]

, but with different basis

functions, under the title “osculating mode”. Specifically, in his article, the interpolation

function is constructed as a linear combination of multiquadric functions plus a series

of polynomials. Additionally, condition equations were used. Hardy

[10]

applied this

technique for the problem of representing a topographic surface by a continuous function

while a set of few discrete data on the surface was given. In his application, he was

seeking to minimize horizontal and vertical displacements of the maximum (hilltops)

and minimum (depressions) points in the resulting surface of topography. Therefore, he

chose an osculating, as he named it, interpolation technique, where the surface coordinates

collocate with the data point coordinates, but also surface tangents coincide at specified

points. The latter was achieved by the usage of low degree polynomials in the interpolation

function. Later on, this interpolation technique was used in various applications, as were

described collectively by Hardy

[11]

. In this article, Hardy

[11]

named this interpolation

technique the multiquadric method.

An additional clarification follows here. In the literature on RBF, the term “condition”

equations is usually used, while in the geodetic literature these equations are referred

to as “constraint” equations since only unknown parameters (and possibly error-free

values) are considered and related to each other. Mikhail

[12]

(p. 213ff.) explained that

constraints for the parameters occur when some or all of the parameters must conform to

some relationships arising from either geometric or physical characteristics of the problem.

Therefore, in the remainder of this article, these equations will only be referred to as

constraint equations.

In this article, it is shown that the approximation problem with the RBF approach,

consisting of a linear combination of thin-plate splines, has a unique solution without

modifying the approximation function, i.e., without additional polynomials and without

constraint equations. The unknown parameters are estimated with the use of the least

squares method. Additionally, it is shown that when constraint equations must be intro-

duced, a least squares solution with constraint equations for the unknown parameters

can be computed. Moreover, some interpolation problems of a few points, with thin-plate

splines, without additional polynomial terms and constraint equations, revealed that they

are well-posed. This article is organized as follows.

In Section 2, the theoretical foundations for the solution of the approximation problem

with thin-plate splines are shown. The result is a linear least squares adjustment problem,

which is well-known, e.g., in Geodesy.

In Section 3, a numerical example from the textbook by Fasshauer

[4]

is presented.

The solution technique used by Fasshauer [4] (p. 170) is described and discussed.

Mathematics 2022,10, 2582 4 of 26

In Section 4, the numerical example of Section 3is used again. However, the solution

is now performed using the adjustment technique presented in Section 2, which allows for

the rigorous consideration of the constraint equations.

In Section 5, the numerical example of Section 3is solved again, but without the

addition of the polynomial terms and the constraint equations for the unknown parameters.

The solution is obtained based on the theory presented in Section 2.

Section 6discusses the investigation results of the surface approximations from

Sections 3–5.

Section 7contains a comparison of eight interpolation problems with thin-plate splines

of four data sets. Finally, in Section 8, conclusions are drawn and some considerations for

further investigations are made.

2. The Scattered Data Approximation Problem

In the following, an overdetermined configuration for the surface approximation

problem with scattered data using RBF is considered. In this case, the parameters can

be determined via least squares adjustment (

–norm minimisation). First, the ordinary

adjustment technique is shown, then the constrained adjustment technique is presented

where the solution is obtained with additional consideration of constraint equations for

the unknowns. Detailed derivations for the solution of adjustment problems with or

without constraint equations can be found in many places in the geodetic literature, e.g., in

the textbooks by Dermanis

[13]

(pp. 4ff., 43ff.), Ghilani and Wolf

[14]

(pp. 173ff., 388ff.),

Mikhail [12] (pp. 159ff., 213ff.) and Niemeier [15] (pp. 112ff., 188ff.), to name a few.

2.1. Ordinary Least Squares Solution

The problem of approximating scattered data is considered, in which a part of the data

is regarded as measurements (observations) subject to random errors. If

n>m

, which is

that the number of measurements

is larger than the number of unknowns

, the following

adjustment problem can be formulated.

A set

X={x1

, ...,

xn}

, with

xi= (xi

yi)∈IR2

, named “data sites”, is given. In

addition, the corresponding measurements are given as scalar-valued data, denoted by

li∈IR

with

1, ...,

. Additionally, a set

Xk={xk1

, ...,

xkm}

, with

xkj= (xkj

ykj)∈IR2

is given and the points are named “center points”. For this problem, the data sites and

the center points are considered error-free while the measurements are subject to random

errors. The function that describes the problem is composed of a linear combination of

radial basis functions

si(xi) =

∑

i=1

ckjΦ(||xi−xkj||2)(1)

where

ckj

1, ...,

are the unknown coefficients to be determined and

Φ(||·||2)

is a

radial basis function which is centered at each center point

. The norm

||·||2

is the

Euclidean norm in IR2. The chosen basic function is the thin-plate spline function

φ(r) = r2log(r)(2)

which is used to derive all the radial basis functions in (1) with r=||xi−xkj||2.

According to the function given in (1), the functional model of this problem can be

written as

li=

∑

i=1

ckjΦ(||xi−xkj||2). (3)

Since an overdetermined system of equations is considered using erroneous measurements,

this system can only be satisfied by (theoretical) true measurements

and true unknowns

ckj

. For real measurements

(where errors are involved), the functional model must be

Mathematics 2022,10, 2582 5 of 26

fulfilled by estimated parameters

ckj

. Therefore, residuals must be introduced, symbolized

by vi,i=1, ..., n, and the observation equations

li+vi=

∑

i=1

ckjΦ(||xi−xkj||2)(4)

can be formed. The measurements

and the residuals

can be represented by the vectors

L= [l1l2. . . ln]T(5)

respectively

v= [v1v2. . . vn]T. (6)

The coefficients of the linear observation equations given in (4), i.e., the radial basis

functions, can be placed in the design matrix

An×m=





Φ(||x1−xk1||2)Φ(||x1−xk2||2). . . Φ(||x1−xkm||2)

Φ(||x2−xk1||2)Φ(||x2−xk2||2). . . Φ(||x2−xkm||2)

.....

Φ(||xn−xk1||2)Φ(||xn−xk2||2). . . Φ(||xn−xkm||2)







(7)

while the estimated unknown parameters in the vector

X= [ˆ

ck1ˆ

ck2. . . ˆ

ckm]T(8)

and the observation equations given in (4) can be expressed as

L+v=Aˆ

X. (9)

To find the least squares solution, the criterion or the objective function of the least

squares method

Ω=vTPv (10)

has to be minimized, where

denotes the weight matrix of the observations. This matrix

results from

P=Q−1

LL (11)

with the cofactor matrix

QLL =1

σ2

ΣLL. (12)

The matrix

ΣLL

is the variance–covariance matrix of the observations and

σ2

is the

arbitrary theoretical (or a priori) variance factor. In the following, it is assumed that the

observations are equally weighted and uncorrelated, so that the weight matrix results as an

identity matrix, P=I. Thus, in this case, (10) can also be written in the form

Ω=

∑

i=1

i. (13)

Rearranging (9) as

v=Aˆ

X−L(14)

Loading more pages...