The goal of the paper is to introduce an efficient algorithm for computation of derivatives of bicubic spline surfaces over equispaced grids with C2 class continuity. The algorithm is based on a recently proposed approach using a special approximation property between quartic and cubic polynomials. The proposed solution replaces the classical de Boor's systems of equations with systems of reduced size and simple remainder explicit formulas. We will show that the proposed new algorithm is numerically equivalent to de Boor's algorithm and the former is more than 50% faster.
Modeling of surfaces plays a key role in a wide variety of computer science fields such as graphics or CAD applications. The paper proposes an efficient computational algorithm for interpolating uniform bicubic spline surfaces of class C2. The standard way to accomplish that is an almost 60 years old classic algorithm designed by Carl de Boor [2] that uses tridiagonal systems of equations. We will refer to it as the full algorithm.
Even though the evaluation of such systems is quite straightforward in general, running in linear time, it can be still viewed as slow especially in real time computing scenarios. For this reason we present a new reduced algorithm for uniform bicubic spline surfaces that, while being in the same complexity class, is faster and needs lower memory requirements.
The structure of the paper is as follows. The remaining part of this section is devoted to the problem statement and introducing terminology. Section 2 presents the new reduced algorithm along with proof of its equality to the full algorithm. The third section provides experimental measurements of actual time savings of the new approach along with some words about its efficient implementation. To be self contained, we provide the equations of de Boor’s full algorithm for surfaces and the reduced algorithm for curves in Appendix, for easier comparison with ours.
Now let’s jump into the problem statement. Consider a
uniform grid
[x0, x1, x2, . . . , xI−1] × [y0, y1, y2, . . . , yJ−1],
(
According to [2], a spline surface is defined by given values zi, j, i = 0, 1, 2, . . . , I − 1, j = 0, 1, 2, . . . , J − 1 at the grid-points, and given first directional derivatives dix, j, i = 0, I − 1, j = 0, 1, 2, . . . , J − 1 at the boundary verticals,
y di, j, i = 0, 1, 2, . . . , I − 1, j = 0, J − 1 at the boundary horizontals and cross derivatives dix,,jy, i = 0, I − 1, j = 0, J − 1 at the four corners of the grid.
The task is to define a quadruple [zi, j, dix, j, diy, j, dix,,jy] at
every grid-point [xi, y j], based on which a bicubic clamped
spline surface S of class C2 is constructed with properties
(
S(xi, y j) = zi, j, ∂ S(xi, y j) ∂ x = dix, j, ∂ S(xi, y j)
∂ y ∂ 2S(xi, y j) ∂ x∂ y
y = di, j, = dix,,jy.
For I = 7 and J = 5, we provide visual illustration of
the input situation in Figure 1, where bold values represent
(
y y y y y y y z0,4 d0,4 z1,4 d1,4 z2,4 d2,4 z3,4 d3,4 z4,4 d4,4 z5,4 d5,4 z6,4 d6,4 x,y d1x,4 d1x,,4y d2x,4 d2x,,4y d3x,4 d3x,,4y d4x,4 d4x,,4y d5x,4 d5x,,4y d6x,4 d6,4 d0x,4 d0,4 x,y y y y y y y y z0,3 d0,3 z1,3 d1,3 z2,3 d2,3 z3,3 d3,3 z4,3 d4,3 z5,3 d5,3 z6,3 d6,3 dx0,3 d0x,,3y d1x,3 d1x,,3y d2x,3 d2x,,3y d3x,3 d3x,,3y d4x,3 d4x,,3y d5x,3 d5x,,3y d6x,3 d6x,,3y y y y y y y y z0,2 d0,2 z1,2 d1,2 z2,2 d2,2 z3,2 d3,2 z4,2 d4,2 z5,2 d5,2 z6,2 d6,2 dx0,2 d0x,,2y d1x,2 d1x,,2y d2x,2 d2x,,2y d3x,2 d3x,,2y d4x,2 d4x,,2y d5x,2 d5x,,2y d6x,2 d6x,,2y y y y y y y y z0,1 d0,1 z1,1 d1,1 z2,1 d2,1 z3,1 d3,1 z4,1 d4,1 z5,1 d5,1 z6,1 d6,1 dx0,1 d0x,,1y d1x,1 d1x,,1y d2x,1 d2x,,1y d3x,1 d3x,,1y d4x,1 d4x,,1y d5x,1 d5x,,1y d6x,1 d6x,,1y y y y y y y y z0,0 d0,0 z1,0 d1,0 z2,0 d2,0 z3,0 d3,0 z4,0 d4,0 z5,0 d5,0 z6,0 d6,0 x,y d1x,0 d1x,,0y d2x,0 d2x,,0y d3x,0 d3x,,0y d4x,0 d4x,,0y d5x,0 d5x,,0y d6x,0 d6,0 d0x,0 d0,0 x,y 2
The reduced algorithm for uniform cubic spline curves was proposed in [11], where the model equations were obtained using a special approximation property between cubic splines and quartic polynomials [10]. This property was extended for spline surfaces proving that a biquartic polynomial fully determines a 2 × 2-component uniform bicubic spline surface of class C2 [9], [5]. This approach was generalized in [3], [6], [4], where a new algorithm for solving the unknown derivatives was proposed and new model equations and explicit formulas were derived. The proposed algorithm introduced equation systems with reduced dimensionality while the remaining unknown derivatives were computed from simple explicit formulas. Measured results showed, that the proposed algorithm is faster than the full one [4]. Thorough analysis has shown that we can further increase the performance of the reduced approach. In the next subsection every model equation will be expressed by simple expressions that correspond to the model equations of [11].
This section consists of two parts. Firstly we show that one of the four systems of de Boor’s algorithm can be equivalently solved by a reduced system and additional explicit formulas. In the second part a reduced algorithm is proposed for the solution of the unknown derivatives of a spline surface of class C2 that is equivalent to de Boor’s full algorithm. 2.1
Consider the uniform grid (
Consider a system of linear equations for the j-th horizontal ( j = 0, 1, 2, . . . , J − 1), see (20) in Appendix, that takes a form 4 1 1 4 0 1 0 1
Lemma 1. Let hx R+, I ≥ 4,
{z0, j, z1, j, z2, j, . . . zI−1, j, d0x, j, dIx−1∈,j} ⊂ R be some known
values and {d1x, j, d2x, j, d3x, j, . . . , dIx−2, j} ⊂ R be unknown
values. Then the full equation system (
Based on Lemma 1 the full systems (21) – (23) can be solved using analogical reduced systems and explicit formulas. Hence we get immediately the reduced algorithm for solving the unknown derivatives of a C2 spline surface that appears to provide faster computations than de Boor’s algorithm.
Theorem 1 (Reduced algorithm for surfaces). If the
values (
3 where i = 2, 4, 6, . . . , I − 3 =
12 hx (zi+2, j − zi−2, j) − hx (zi+1, j − zi−1, j), (12) Step 1b. Computation of dx along the horizontals from explicit formulas for inner odd-indexed grid-points. For each i = 1, 3, 5, . . . , I − 2 and j = 0, 1, 2, . . . , J − 1, Step 2a. Computation of dy along the verticals from equation systems for inner even-indexed grid-points. For each i = 0, 1, 2, . . . , I − 1, 3 12 hy (zi, j+2 − zi, j−2) − hy (zi, j+1 − zi, j−1), (14) where j = 2, 4, 6, . . . , J − 3 Step 2b. Computation of dy along the verticals from explicit formulas for inner odd-indexed grid-points. For each i = 0, 1, 2, . . . , I − 1 and j = 1, 3, 5, . . . , J − 2, diy, j = 3 4hy
1 (zi, j+1 − zi, j−1) − 4 (diy, j+1 + diy, j−1). (15) Step 3a. Computation of dx,y along the horizontals from equation systems for inner even-indexed grid-points. For each j = 0, J − 1 , solve_system( dix+,y2, j − 14dix,,jy + dix−,y2, j = =
h3x (dix+2, j − dix−2, j) − 1h2x (dix+1, j − dix−1, j), (16) where i = 2, 4, 6, . . . , I − 3 Step 3b. Computation of dx,y along the horizontals from explicit formulas for inner odd-indexed grid-points. For each i = 1, 3, 5, . . . , I − 2 and j = 0, J − 1, dix,,jy = 3 4hx
x 1 (dix+1, j − di−1, j) − 4 (dix+,y1, j + dix−,y1, j). (17) Step 4a. Computation of dx,y along the verticals from equation systems for inner even-indexed grid-points. For each i = 0, 1, 2, . . . , I − 1, solve_system( dix,,jy+2 − 14dix,,jy + dix,,jy−2 = =
h3y (diy, j+2 − diy, j−2) − 1h2y (diy, j+1 − diy, j−1), (18) where j = 2, 4, 6, . . . , J − 3 Step 4b. Computation of dx,y along the verticals from explicit formulas for inner odd-indexed grid-points. for each i = 0, 1, 2, . . . , I − 1 and i = 1, 3, 5, . . . , I − 2, dix,,jy = 3 4hy
y 1 (diy, j+1 − di, j−1) − 4 (dix,,jy+1 + dix,,jy−1). (19)
If I is even, then the last model equation in steps (12)
and (16) needs to be accordingly replaced by the model
equation derived according to (
We underline that the systems (12), (14), (16) and (18) are diagonally dominant and use the same tridiagonal matrix with constant main diagonal and upper, lower diagonals.
In the next section it will be shown that the reduced algorithm is faster than the full one and additionally it needs lower memory requirements.
The main difference between original full algorithm and the reduced one is that the latter computes only half of the unknown values using LU factorization where the remaining half of the unknowns are solved by simple remainder formulas. Because of this, the reduced algorithm is supposed to be more efficient in terms of execution speed and memory requirements as well. We will discuss and measure it’s speed in Section 3. 3
To be worth of actual implementation, the reduced algorithm should be significantly faster than full one. The elementary task in both algorithms is computation of tridiagonal systems of equations. As the full and the reduced systems are diagonally dominant we use a modified LU factorization method as the internal solver for such systems [4]. Our optimizations of the LU factorization relies on a fact, that the left-hand side matrices comprise only constant values.
While the LU factorization for tridiagonal linear systems is a quite efficient way to solve the equations which are interdependent, it has some drawbacks in the sense of modern CPU’s capabilities like usage of vectorized SIMD instructions especially due to the structure of tridiagonal matrices. By breaking down the computational task of the full approach into reduced systems and simple explicit formulas, half of the equations are independent and therefore more versatile to manual or automatic compiler optimizations.
Memory requirements For the sake of completeness some words about applied data structures and memory requirements should be given.
The input grid (
The needs of the full and reduced algorithms are quite low regarding the size of the grid. The full tridiagonal systems in Lemma 2 require 2 · max(I, J) space to store the right-hand side vector and an auxiliary buffer vector used for the LU factorization. In case of the reduced algorithm, only half number of the equations form the tridiagonal system, therefore they require only max(I, J) space for the LU part.
Data structures Since the grid may contain tens of thousands or more grid-points, the most effective representation is the jagged array structure for each of the z, dx, dy, dx,y values. Each equation system from any of the two algorithms always depends on one row of a jagged array, therefore entire rows of the jagged structure can be effectively cached under the assumption, that the size of the row is not very large. Notice that the iterations for computing the diy, j and most of the dix,,jy values in both algorithms have interchanged indices compared to the iteration throughout the dix, j values. We mention that an efficient implementation needs to setup the jagged arrays in accordance with how we want to iterate the data [7]. 3.1
Now let us compare the implementations of both algorithms. We implemented a benchmark in C++ 17 compiled with MSVC 2017 using -O2 optimization level and individual code generation for each tested CPU, i.e. CPUs with AVX2 support were running binaries compiled to AVX2 instruction set whereas older CPUs received for instance only SSE2 compiled binaries. Testing environments comprised several computers with various CPUs ranging from rather obsolete Nehalem to recent Skylake microarchitecture. All systems had 8 – 32 GB of RAM, SSD and Windows 10 installed. The tests were conducted on freshly booted PCs after 2 minutes of idle time without running any non-essential services or processes like web browsers, database engines, etc.
On all testing computers the tests were conducted on
two datasets, a small one on a grid of size I, J = 100,
and a large one, where grid dimensions were I, J =
1000. Both datasets comprised the grid [x0, x1, . . . , xI−1] ×
[y0, y1, . . . , yJ−1], where x0 = −20, xyI−1 = 20, y0 = −20,
yJ−1 = 20 and the values zi, j, dix, j, di, j, dix,,jy, see (
Tables 1 and 2 contain results for both datasets and consist of four columns. The first column contains the tested CPUs ordered by their release date. Columns two through four contain measured execution times in microseconds for both algorithms and their speed ratios.
To provide better comparison for more datasets, Table 3 provides measurements on more different input sizes, however for the sake of readability it contains measurements only from Intel i7 6700K as the fastest testing CPU.
Let us review the measured performance improvement of the reduced algorithm. Results of Table 1 say that the reduced algorithm is more than fifty percent faster than the full one as long as the grid dimensions remains relatively small to enable the entire rows of the grid as well as the auxiliary LU vectors to be cached. CPU
The speedup of the reduced algorithm in Table 2 with a larger dataset is not so high due to difficulties of keeping data in fast but small L1 cache, resulting in higher ratio of cache misses. In the case of even larger dataset, let’s say in the order of billions of grid-points, the performances of both algorithms are similar with the reduced one gaining only small advantage.
In our experiments we observed that the reduced algorithm was always faster than the full one, however exact speedup depends on the size of the input grid. The general rule is: the larger the grid the smaller the speedup.
The highest speedup does not depend on the grid size if it is from range of 50 to 200.
CPU I, J = 50 I, J = 100 I, J = 200 I, J = 300 I, J = 400 I, J = 1000 I, J = 1500 I, J = 2000
Full
Reduced
Speedup Let us briefly discuss the new algorithm from the numerical and experimental point of view. While the classic full algorithm is composed of four series of tridiagonal systems of equations, the reduced algorithm breaks down each equation system into a reduced one approximately half the size of the original one and simple mutually independent explicit formulas.
In addition, from the numerical point of view, the
reduced tridiagonal subsystems are diagonally dominant and
therefore computationally stable [1], similarly to the full
systems. The remainder explicit formulas (
The maximal numerical difference in our C++ implementation was in the order of 10−16 on several different datasets so the reduced algorithm yields numerically accurate results.
Since the algorithm consists of many independent systems of linear equations, it can be also effectively parallelized for both CPU and GPU architectures.
There is also a future work for further reduction of the number of equations in the tridiagonal systems. 5
The paper introduced a new algorithm to compute the unknown derivatives for uniform bicubic spline surfaces of class C2. The algorithm reduces the size of the tridiagonal systems of equations by half and computes the remaining unknown derivatives using simple explicit formulas. A substantial decrease of execution time of derivatives at grid-points has been achieved with lower memory space requirements at the cost of a slightly more complex implementation where the measured speedup ranges from 1.3 to 1.9 depending on the grid size and CPU architecture.
This work was partially supported by projects Technicom ITMS 26220220182 and APVV-15-0091 Effective algorithms, automata and data structures.
Carl de Boor in [2] proposed an algorithm for computation of the unknown derivatives of spline surface of class C2 over any grid with four types of full systems of linear equations. The following lemma reformulates this algorithm for uniform grid.
Lemma 2 (Full algorithm for surfaces). If the values (
For each i = 0, 1, 2, . . . , I − 1, solve_system( diy, j+1 + 4 · diy, j + dix, j−1 = h3y · (zi, j+1 − zi, j−1), where j = 1, 2, 3, . . . , J − 2 Step 3. Computation of dx,y along the horizontals from equation systems.
For each j = 0, J − 1, solve_system( dix+,y1, j + 4 · dix,,jy + dix−,y1, j = h3x · (dix+1, j − dix−1, j), (22) where i = 1, 2, 3, . . . , I − 2 Step 4. Computation of dx,y along the verticals from equation systems.
For each i = 0, 1, 2, . . . , I − 1, solve_system( diz,,jy+1 + 4 · diz,,jy + diz,,jy−1 = h3y · (diy, j+1 − diy, j−1), (23) where j = 1, 2, 3, . . . , J − 2
A new algorithm was proposed in [11] for computation of the unknown derivatives of a spline curve of class C2 over uniform grids with a reduced system of linear equations instead of the well known full algorithm with parameters 1,4,1.
Lemma 3 (Reduced algorithm for curves). Consider a cubic clamped spline of class C2 over a uniform grid [u0, u1, . . . , uN−1], where ui = u0 + ih, i = 1, 2, 3, . . . , N − 1, −14 1 0 defined by given values y 0, . . . , yN−1 and d0, dN−1 at gridpoints. The tridiagonal system
1