A Comparative Study of Methods for Scattered
Data Interpolation Using Minimum Norm
Networks and Quartic Triangular Bézier Surfaces
Krassimira Vlachkova 1 and Krum Radev 1
1
 Faculty of Mathematics and Informatics, Sofia University “St. Kliment Ohridski”, 5
James Bourchier Blvd., Sofia, 1164, Bulgaria


             Abstract
             We consider the problem of scattered data interpolation in R3 using curve
             networks extended to smooth interpolation surfaces. Nielson (1983)
             proposed a solution that constructs smooth interpolation curve network
             with minimal L2-norm of the second derivative. The obtained minimum
             norm network (MNN) is cubic. Vlachkova (2020) generalized Nielson’s
             result to smooth interpolation curve networks with minimal Lp-norm of
             the second derivative for 1<p<∞. Vlachkova and Radev (2020) proposed
             an algorithm that degree elevates the MNN to quartic curve network and
             then extends it to a smooth surface consisting of quartic triangular Bézier
             surfaces. Here we apply this algorithm to the following two curve networks:
             (i) the MNN which is degree elevated to quartic; (ii) the minimum Lp-norm
             network for p=3/2 which is slightly modified to quartic. We evaluate and
             compare the quality and the shape of the obtained surfaces with respect to
             different criteria. We performed a large number of experiments using data
             of increasing complexity. Here we present and comment the results of our
             experiments.

             Keywords
             Scattered data interpolation, curve network, minimum norm network,
             spline, Bézier surface

1. Introduction
     Interpolation of data points in ℝ3 by smooth surface is a fundamental problem
in applied mathematics which finds applications in a variety of fields such as medi-
cine, architecture, archeology, computer graphics and animation, bioinformatics,
scientific visualization, and more. In general the problem can be formulated as
follows: Given a set of points (xi, yi, zi) ∈ ℝ3, i = 1, ..., n, find a bivariate function
Information Systems & Grid Technologies: Fifteenth International Conference ISGT’2022, May 27–28, 2022, Sofia, Bulgaria
EMAIL: krassivl@fmi.uni-sofia.bg (K. Vlachkova); kvradev@uni-sofia.bg (K. Radev)
ORCID: 0000-0002-8832-7336 (K. Vlachkova); 0000-0002-6515-169X (K. Radev)

            © 2022 Copyright for this paper by its authors.
            Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
            CEUR Workshop Proceedings (CEUR-WS.org)
F(x, y) defined in a certain domain D containing points Vi = (xi, yi), such that F pos-
sesses continuous partial derivatives up to a given order, and F(xi, yi) = zi.
     Various methods and approaches for solving this problem were proposed and
discussed, see, e.g., the surveys [1, 2, 3], and also [4, 5]. A standard approach to
solve the problem consists of two steps, see [1]:
     1. Construct a triangulation T = T(V1, ..., VN);
     2. For every triangle in construct a surface which interpolates the data in
     the three vertices.
     The interpolation surface constructed in Step 2 is usually polynomial or
piecewise polynomial. Typically, the patches are computed with a priori pre-
scribed normal vectors at the data points G1 or G2 smoothness of the resulting
surface is achieved either by increasing the degree of the patches, or by the so
called splitting [6] in which for each triangle in T a macro-patch consisting of a
fixed number of Bézier sub-patches is constructed. Splitting allows to keep the
degree of the Bézier patches low by increasing the degrees of freedom. In prac-
tice using patches of least degree and splitting is preferable since it is computa-
tionally simple and efficient.
     Shirman and Séquin [7, 8] construct a G1 smooth surface consisting of quartic
triangular Bézier surfaces. Their method assumes that the normal vectors at the
data points are given as part of the input. Shirman and Séquin construct a smooth
cubic curve network defined on the edges of T, first, and then degree elevate it
to quartic. Next, they apply splitting where for each triangle in T a macro-patch
consisting of three quartic Bézier sub-patches is constructed. To compute the in-
ner Bézier control points closest to the boundary of the macro-patch, Shirman and
Séquin use a method proposed by Chiyokura and Kimura [9, 10]. The interpola-
tion surfaces constructed by Shirman and Séquin’s algorithm often suffer from
unwanted bulges, tilts, and shears as pointed out by the authors in [11] and more
recently by Hettinga and Kosinka in [12].
     Nielson [13] proposed a three-steps method for solving the interpolation
problem as follows:
     Step 1. Triangulation. Construct a triangulation T of Vi, i = 1, ..., n. The do-
main D is the union of all triangles in T.
     Step 2. Minimum norm network. The interpolant F and its first order partial
derivatives are defined on the edges of T to satisfy an extremal property. The
resulting MNN is a cubic curve network, i. e. on every edge of T it is a cubic
polynomial.
     Step 3. Interpolation surface. The MNN obtained is extended to F by an
appropriate blending method based on convex combination schemes. Nielson’s
interpolant F is a rational function on every triangle in T.
     Andersson et al. [14] focused on Step 2 of the above method, namely the
construction of the MNN. The authors gave a new proof of Nielson’s result by

                                         160
using a different approach. They constructed a system of simple linear curve net-
works called basic curve networks and then represented the second derivative
of the MNN as a linear combination of these basic curve networks. The new ap-
proach allows to consider and handle the case where the data are convex and we
seek a convex interpolant. Andersson et al. formulate the corresponding extremal
constrained interpolation problem of finding a minimum norm network that is
convex along the edges of the triangulation. The extremal network is character-
ized as a solution to a nonlinear system of equations.
     Vlachkova and Radev [15] proposed an algorithm for interpolation of data
in ℝ3 which improves on Shirman and Séquin’s approach in the following way.
First, they use Nielson’s MNN which is degree elevated to quartic curve network.
Second, they extend it to a smooth surface consisting of quartic triangular Bézier
patches by applying different strategy for computation of the control points. A
significant advantage of Nielson’s method is that the normal vectors at the data
points are obtained through the computation of the MNN. Moreover, during the
computation of control points closest to the boundary of a macro-patch, Vlach-
kova and Radev [15] adopt additional criteria so that to avoid unwanted distor-
tions and twists which appear in surfaces constructed by Shirman and Séquin’s
method. As a result, the quality of the resulting surfaces is improved.
     Vlachkova [16] extended Nielson’s MNN to minimum Lp-norm networks for
1 < p < ∞. The most important conclusions in [16] are as follows:
     • the results allow the efficient computation of the minimum Lp-norm net-
     works for 1 < p < ∞;
     • the normal vectors at the data points are obtained simultaneously with the
     computation of the minimum Lp-norm networks;
     • the minimum Lp-norm networks are obtained through a global optimiza-
     tion which improves their shape;
     • the results allows the construction of optimized polynomial minimum
     norm networks of a priori given degree.
     In this paper, we consider the following two curve networks:
     (i) Nielson’s MNN which is degree elevated to quartic;
     (ii) the minimum Lp-norm network for p = 3/2 which is slightly modified to
      quartic.
      Using the algorithm proposed in [15], we construct the corresponding two
interpolation surfaces consisting of quartic triangular Bézier patches. Our goal
is to evaluate and compare the quality and the shape of these surfaces. We have
chosen the following criteria for comparison:
      (i) the highlight-line algorithm [17];
       (ii) the color plot of the Gaussian curvature;



                                       161
     ( iii) the maximum distance between the function sampled at the data points
      and the corresponding interpolant.
      Our work and contributions presented here are in the field of experimental
algorithmics. We share and comment on the observations from the experiments
performed, which will help to further optimize the quality and the shape of the
interpolation surfaces generated with the algorithm proposed in [15].

2. Preliminaries and related work
      Let n ≥ 3 be an integer and Pi: = (xi, yi, zi), i = 1, ..., n be different points in
ℝ3. We call this set of points data. We assume that the projections Vi = (xi, yi) of
the data Oxy into the plane are different and non-collinear. Stressing the fact that
Vi are in a general position, although not necessarily irregularly placed, we call
such data scattered.
      Definition 1. A collection of non-overlapping, non-degenerate triangles in is
a triangulation of the points Vi, i = 1, ..., n, if the set of the vertices of the triangles
coincides with the set of the points Vi, i = 1, ..., n.
      Hereafter we assume that a triangulation T of the points Vi, i = 1, ..., n, is
given and fixed. Furthermore, for the sake of simplicity, we assume that the do-
main D formed by the union of the triangles in T is connected. In general D is
a collection of polygons with holes. The set of the edges of the triangles in T is
denoted by E. If there is an edge between Vi and Vj in E, it will be referred to by
eij or simply by e if no ambiguity arises.
      Definition 2. A curve network is a collection of real-valued univariate func-
tions {fe}e∈E defined on the edges in E.
      With any real-valued bivariate function defined on we naturally associate the
curve network defined as the restriction of on the edges in , i. e. for ,

                                                                  (1)

     where 0 ≤ t ≤ || e || and || e || =
     Furthermore, according to the context F will denote either a real-valued bi-
variate function or a curve network defined by (1). For p, such that 1 < p < ∞, we
introduce the following class of smooth interpolants




     where is the class of bivariate continuous functions defined in D, AC is the
class of univariate absolutely continuous functions defined in [0, || e ||], and Lp
is the class of univariate functions defined in [0, || e ||] whose p-th power of the


                                            162
absolute value is Lebesgue integrable. The restrictions on E of the functions in Fp
form the corresponding class of so-called smooth interpolation curve networks
                                                                       (2)
     The smoothness of the interpolation curve network F ∈ Cp(E) geometrically
means that at each point Pi there is a tangent plane to F, where a plane is tangent
to the curve network at a point Pi if it contains the tangent vectors at Pi of the
curves incident to Pi.
     Inner product and Lp-norm are defined in Cp(E) by




    where F ∈ Cp(E) and G: = {ge}e∈E ∈ Cp(E). We denote the networks of the
second derivative of F by F″: = {fe″}e∈E and consider the following extremal
problem:


     For i = 1, ..., n let mi denote the degree of the vertex Vi, i. e. the number of
the edges in E incident to Vi. In [14] for any pair of indices is, such that i = 1, ...,
n and s = 1, ..., mi – 2, a basic curve network Bis is defined on E. The basic curve
networks are linear curve networks of minimal support.
     Let q be the conjugate of p, i.e. 1/p + 1/q = 1. In [16] a full characterization of
the solution F * to the extremal problem (Pp) for 1 < p < ∞ was made. It was shown
that the curve network F * ∈ Cp(E) solves problem (Pp) for 1 < p < ∞ if and only if



     where                                                The coefficients αis are ob-
tained as the unique solution to a system of equations which is nonlinear except
in the case p = 2 where it is linear. We note that Nielson’s MNN is obtained for p
= 2. We also note that in the case where q is an even number then the correspond-
ing minimum Lp-norm network is a polynomial curve network of degree q + 1.
In the case where q is an odd number then on every edge of the triangulation the
corresponding minimum Lp-norm network is either a polynomial of degree q + 1,
or a spline of degree q + 1 with one knot.
     Now we briefly discuss the constrained interpolation problem of finding a
minimum norm network which is convex along the edges of T. We recall that this
problem was set and solved in [14].
     For a given triangulation, there is a unique continuous function L : D → ℝ1
that is linear inside each of the triangles of T and interpolates the data.


                                          163
     Definition 3. Scattered data D in are convex if there exists a triangulation T
of Vi, i = 1, ..., n, such that the corresponding function L is convex. The data are
strictly convex if they are convex and the gradient of L has a jump discontinuity
across each edge inside D.
     We introduce the class of smooth interpolation edge convex curve networks


    and consider the following constrained extremal interpolation problem



    In [14] it was shown that in the case of strictly convex data problem has a
unique solution (the edge convex MNN) such that



     where (x)+ : = max(x, 0), and αis ∈ ℝ. The coefficients αis are obtained as the
unique solution to a nonlinear system of equations. The solution F̂ * on each edge
in T is either a convex cubic polynomial, or a convex cubic spline with one knot
consisting of a linear function and a convex cubic polynomial.
     In this paper we consider the unconstrained problem (Pp) for p = 3/2 and
also the constrained problem (P̂ ). The solution of (Pp) for p = 3/2 (hence q =3) on
each edge in T is a quartic spline with at most one knot. For both problems, we
slightly modify their solutions to obtain a polynomial curve networks (quartic for
(Pp) for p = 3/2, and cubic for (P̂ )) while preserving the same tangent planes at
the data points.




Figure 1: Example 1: (left). The triangulation; (right) The MNN (p=2)




                                        164
Figure 2: Example 1. Example 1. The surface interpolating the MNN: (left)
Highlight lines on the surface; (right) Gaussian curvature: the color scale goes
from blue (low) to red (large)




Figure 3: Example 1. The surface interpolating the modified edge convex
MNN: (left) Highlight lines on the surface; (right) Gaussian curvature: the
color scale goes from blue (low) to red (large)

3. Results from the experiments
     Here we present and comment on two examples from the experiments per-
formed and we share our observations and conclusions. We use BezierView [18]
to visualize the highlight lines and the Gaussian curvature of the interpolation
surfaces.
     Example 1. We consider data obtained from a symmetric triangular pyra-
mid. We have n = 4,
V4 = (0,0), and zi = 0, i = 1,2,3, z4 = –1/2. The triangulation and the correspond-
ing MNN are shown in Figure 1. The highlight lines on the interpolation surface
are visualized in Figure 2 (left). The Gaussian curvature is visualized in Figure
2 (right) where the color scale goes from blue (low) to red (large). The mini-
mum value of the Gaussian curvature is -11.8318, the maximum value is 34.3368.
The surface that interpolates the edge convex MNN is shown in Figure 3. The
highlight lines of this surface are visualized in Figure 3 (left) and its Gaussian
curvature is visualized in Figure 3 (right). The minimum value of the Gaussian
curvature is -0.6046, and the maximum value is 36.1598.



                                       165
Table 1
The data for Example 2




      Example 2. We consider the convex function f = exp((x – 0.5)2 + (y – 0.5)2)
which is sampled at 25 points shown in Table 1. The triangulation, which is
shown in Figure 4 (left), is the Delaunay triangulation. The corresponding edge
convex MNN is shown in Figure 4 (right). The highlight lines on the surface are
visualized in Figure 6.
      For Example 1, see Figure 2, we can see that our interpolant for the uncon-
strained case is visually pleasant: the highlight lines are smooth with a small
number of inflection points (if any), and the Gaussian curvature is evenly dis-
tributed. The highlight lines of the surface interpolating the edge convex MNN
in Figure 3 (left) also look nicely as in the unconstrained case but the Gaussian
curvature in Figure 3 (right) is not evenly distributed.
      For Example 2 where the data are sampled from the convex exponential func-
tion, it is clearly seen that the interpolation surface generated from our algorithm
presented in [15], see Figure 5 (right), significantly improves on the surface gener-
ated from Shirman and Sequin’s algorithm [7, 8], see Figure 5 (left). Although at
first glance the surface in Figure 5 (right) visually even appears convex, when we
examine the highlight lines in Figure 6 in detail, it is clearly seen that it is not convex.
      The above observations shows that the algorithm presented in [15] can be
refined and improved further. Another conclusion obtained from our preliminary
numerical experiments is that using the modified quartic minimum Lp-norm net-
work for p = 3/2 gives better results with respect to the maximum distance than
using the degree elevated MNN.

4. Conclusion and future work
     Given the importance of surface modeling and simulation techniques in
practice, it is important to better understand and create interpolation surfaces with
good approximation properties based on different criteria. In the future we intend
to enlarge our work by analyzing in-depth more data sets, which will support
further optimization of our algorithm [15]. More details will be presented in the
full version of this paper.


                                            166
Figure 4: Example 2: (left) The Delaunay triangulation for n=25; (right) The
edge convex MNN




Figure 5: Comparison of two interpolation surfaces for the data in Example 2:
(left) The surface generated using Shirman and Sequin’s algorithm; (right) The
surface generated using Algorithm 2 in [15]




Figure 6: Example 2: The highlight lines of the surface interpolating the edge
convex MNN and generated using Algorithm 2 in [15]



                                      167
5. Acknowledgments
     This work was supported in part by Sofia University Science Fund Grant
No. 80-10-109/2022, and European Regional Development Fund and the Op-
erational Program “Science and Education for Smart Growth” under contract №
BG05M2OP001-1.001-0004 (2018-2023).

6. References
[1]   S. Mann, C. Loop, M. Lounsbery, D. Meyers, J. Painter, T. DeRose, K.
      Sloan, A survey of parametric scattered data fitting using triangular interpo-
      lants, in: H. Hagen (Ed.), Curve and Surface Design, SIAM, Philadelphia,
      1992, pp. 145–172. doi:10.1137/1.9781611971651.ch8.
[2]   S. Lodha, K. Franke, Scattered data techniques for surfaces, in: Proceed-
      ings of Dagstuhl Conference on Scientific Visualization, IEEE Com-
      puter Society Press, Washington, 1997, pp. 182–222. doi:10.1109/DAG-
      STUHL.1997.1423115.
[3]   M. Berger, A. Tagliasacchi, L. Seversky, P. Alliez, G. Guennebaud, J.
      Levine, A. Sharf, C. Silva, A survey of surface reconstruction from point
      clouds, Comput. Graph. Forum 36 (2017) 301–329. doi:10.1111/cgf.12802.
[4]   T. K. Dey, Curve and surface reconstruction: algorithms with mathematical
      analysis, Cambridge Monographs on Applied and Computational Mathemat-
      ics, Cambridge University Press, 2006. doi:10.1017/CBO9780511546860.
[5]   K. Anjyo, J. Lewis, F. Pighin, Scattered data interpolation for computer
      graphics, SIGGRAPH 2014 course notes, 2014. URL: http://olm.co.jp/rd/
      research_event/scattered-datainterpolation-for-computer-graphics, last ac-
      cessed April 15, 2021.
[6]   R. Clough, J. Tocher, Finite elements stiffness matrices for analysis of plate
      bending, in: Proceedings of the 1st Conference on Matrix Methods in Struc-
      tural Mechanics, volume 66–80, Wright-Patterson A. F. B., Ohio, 1965, pp.
      515–545. URL: http://contrails.iit.edu/reports/8574.
[7]   L. Shirman, C. Séquin, Local surface interpolation with Bézier patch-
      es, Comput. Aided Geom. Des. 4 (1987) 279–295. doi:10.1016/0167-
      8396(87)90003-3.
[8]   L. Shirman, C. Séquin, Local surface interpolation with Bézier patches:
      errata and improvements, Comput. Aided Geom. Des. 8 (1991) 217–221.
      doi:10.1016/0167-8396(91)90005-V.
[9]   H. Chiyokura, F. Kimura, Design of solids with free-form surfaces, in: P. P.
      Tanner (Ed.), SIGGRAPH ’83 Proceedings of the 10th Annual Conference
      on Computer Graphics and Interactive Techniques, volume 17, ACM, New
      York, 1983, pp. 289–298. doi:10.1145/964967.801160.


                                        168
[10] H. Chiyokura, Localized surface interpolation method for irregular meshes,
     in: T. Kunii (Ed.), Advanced Computer Graphics, Proceedings of Comput-
     er Graphics Tokyo’86, volume 66–80, Springer, Tokyo, 1986, pp. 3–19.
     doi:10.1007/978-4-431-68036-9_1.
[11] L. Shirman, C. Séquin, Local surface interpolation with shape parameters
     between adjoining Gregory patches, Comput. Aided Geom. Des. 7 (1990)
     375–388. doi:10.1016/0167-8396(90)90001-8.
[12] G. Hettinga, J. Kosinka, Multisided generalisations of Gregory patch-
     es, Comput. Aided Geom. Des. 62 (2018) 166–180. doi:10.1016/j.
     cagd.2018.03.005.
[13] G. Nielson, A method for interpolating scattered data based upon a mini-
     mum norm network, Math. Comput. 40 (1983) 253–271. doi:10.1090/
     S0025-5718-1983-0679444-7.
[14] L. Andersson, T. Elfving, G. Iliev, K. Vlachkova, Interpolation of convex
     scattered data in R3 based upon an edge convex minimum norm network, J.
     of Approx. Theory 80 (1995) 299 – 320. doi:10.1006/jath.1995.1020.
[15] K. Vlachkova, K. Radev, Interpolation of data in R3 using quartic tri-
     angular Bezier surfaces, in: AIP Conf. Proc., volume 2325, 2021.
     doi:10.1063/5.0040457.
[16] K. Vlachkova, Interpolation of scattered data in R3 using minimum
[17] Lp-norm network, 1 < p < ∞, J. Math. Anal. Appl. 482 (2020) 123824.
     doi:10.1016/j.jmaa.2019.123824.
[18] K.-P. Beier, Y. Chen, Highlight-line algorithm for realtime surface-quality
     assessment, Comput.-Aided Des. 26(4) (1994) 268–277. doi:10.1016/0010-
     4485(94)90073-6.
[19] BezierView, https://www.cise.ufl.edu/research/SurfLab/bview, 2015.
     SurfLab, University of Florida, Last accessed: April 15, 2022.




                                      169