Implicit progressive-iterative approximation for curve and surface reconstruction

Computer Aided Geometric Design 77 (2020) 101817

Contents lists available at ScienceDirect

Computer Aided Geometric Design www.elsevier.com/locate/cagd

Implicit progressive-iterative approximation for curve and surface reconstruction ✩ Yusuf Fatihu Hamza a , Hongwei Lin a,b,∗ , Zihao Li b a b

School of Mathematical Science, Zhejiang University, Hangzhou, 310027, China State Key Laboratory of CAD&CG, Zhejiang University, Hangzhou, 310058, China

a r t i c l e

i n f o

Article history: Received 21 June 2019 Received in revised form 29 December 2019 Accepted 3 January 2020 Available online 23 January 2020 Keywords: Implicit curve and surface Curve and surface fitting Progressive-iterative approximation

a b s t r a c t Implicit curve and surface reconstruction attracts the attention of many researchers and gains a wide range of applications, due to its ability to describe objects with complicated geometry and topology. However, extra zero-level sets or spurious sheets arising in the reconstruction process make the reconstruction result challenging to be interpreted and damage the final result. In this paper, we propose an implicit curve and surface reconstruction method based on the progressive-iterative approximation, named implicit progressive-iterative approximation (I-PIA). The proposed method elegantly and naturally eliminates the spurious sheets without requiring any explicit minimization procedure, thus significantly reducing the computing cost and providing high-quality reconstruction results. Numerical examples are presented to demonstrate the efficiency and effectiveness of the proposed method. © 2020 Elsevier B.V. All rights reserved.

1. Introduction Implicit representation and parametric representation are two common representation techniques in geometric design. For parametric representation, it is difficult to fit a data set with complicated geometry, and parameterization is always a challenging problem. Without requiring any parameterization, implicit functions can describe an object with complicated geometry and supply a flexible and smooth surface representation. Thus, implicit surface reconstruction receives great attention due to its capability to create an object with complicated topology and geometry. However, the extra zero-level sets generated in the procedure of implicit curve and surface reconstruction make the reconstruction results challenging to be interpreted, and damage the resulting curve and surface. To eliminate the extra zero-level sets, regularization terms are usually required to be added in the objective functions of the minimization problems for implicit curve and surface reconstruction. For example, Liu et al. (2017) incorporated the total variation of implicit representation to reduce the appearance of the extra zero-level sets as minimum as possible. In Refs. (Jüttler and Felis, 2002; Rouhani and Sappa, 2011; Rouhani et al., 2014; Yang et al., 2005), tension terms are added in the minimization problem to get rid of extra zero-level sets, and thus a singular system of equations is avoided. Because the forms of the regularization terms are usually complicated, their addition to the objective functions seriously affects the efficiency of the implicit curve and surface reconstruction algorithms.

✩

*

Editor: Thomas A. Grandine. Corresponding author at: School of Mathematical Science, Zhejiang University, Hangzhou, 310027, China. E-mail address: [email protected] (H. Lin).

https://doi.org/10.1016/j.cagd.2020.101817 0167-8396/© 2020 Elsevier B.V. All rights reserved.

2

Y.F. Hamza et al. / Computer Aided Geometric Design 77 (2020) 101817

Progressive-iterative approximation (PIA) is a series of efficient data fitting methods with intuitive geometric meaning. They have been extensively employed in parametric curve and surface fitting, and subdivision curve and surface fitting (Chen et al., 2008; Cheng et al., 2009; Deng and Lin, 2014; Deng and Ma, 2012; Lin, 2010; Lin et al., 2018a, 2004; Lin and Zhang, 2011; Lin et al., 2005; Lu, 2010), but never been used in implicit curve and surface reconstruction. In this paper, we develop a progressive-iterative approximation method for implicit curve and surface reconstruction, named implicit progressive-iterative approximation (I-PIA). We prove the convergence of I-PIA, and show that the I-PIA method itself naturally solves a minimization problem with a regularization term, without any extra computation effort. Therefore, the I-PIA method can not only eliminate the extra level sets, but also, more importantly, improve the reconstruction efficiency of the implicit curves and surfaces. Numerous of numerical examples illustrated in this paper show that, by I-PIA, the implicit curve and surface reconstruction time is reduced one to three orders of magnitude, compared with the state-of-the-art implicit curve and surface reconstruction methods. In conclusion, the main contributions of this paper include:

• I-PIA naturally solves a minimization problem with a regularization term in its iteration procedure, without any extra computation effort.

• No extra zero-level set exists in the iteration of I-PIA, so the reconstruction result generated by I-PIA is very clear. • The implicit curve and surface reconstruction time by I-PIA is reduced at least one to three orders of magnitude, compared with the state-of-the-art methods. The structure of this paper is as follows. In Section 1.1, we review the related work. Preliminary definitions and the statement of the problems are given in Section 2. In Section 3, we present the I-PIA method for implicit curve and surface reconstruction. Experimental results and discussions are given in Section 4. Finally, we conclude the paper in Section 5. 1.1. Related work In this section, some work related to implicit curve and surface reconstruction and progressive iterative approximation (PIA) will be briefly reviewed. Implicit curve and surface reconstruction: Most implicit surface reconstruction algorithms blend local implicit primitives to represent surfaces based on the idea developed by Blinn (1982). Muraki (1991) presented the Blobby model to fit complicated data set by blending implicit primitives. Hoppe et al. (1992) proposed a surface reconstruction algorithm based on the locally defined signed distance function. Curless and Levoy (1996) used the volumetric representation consisting of a cumulative weighted signed distance function. Carr et al. (2001) proposed a fast method for fitting and evaluating radial basis functions (RBFs) to model a large data set. Recently, Huang et al. (2019) presented the variational implicit point set surfaces (VIPSS) for constructing an implicit surface, which particularly suits for sparse and non-uniform sampling. So far, many schemata had been considered to reduce the amount of storage cost in implicit reconstruction. Morse et al. (2005) proposed using compactly supported radial basis functions to reduce the computational cost and memory requirement, which allows reconstructing surface from large data sets that are impractical in the previous method like (Turk and O’brien, 1999). Furthermore, the authors in Refs. (Kojekine et al., 2003; Ohtake et al., 2005) employed compactly supported RBF to reduce the computational cost and improve the efficiency of the reconstruction process. Pan et al. (2016) reduced the storage requirement of implicit reconstruction by incorporating a low-rank tensor approximation technique in the reconstruction process. Ohtake et al. (2003) proposed a multi-level partition of unity (MPU) representation to reconstruct surface models from a huge set of points. Wang et al. (2011) presented a surface reconstruction algorithm based on the implicit PHT-spline, which reconstructs a surface from a large point cloud efficiently. Surface reconstruction was restated to be a Poisson problem by Kazhdan et al. (2006), and the proposed method of Poisson surface reconstruction approximates the indicator function of the surface. Moreover, screened Poisson surface reconstruction was presented to avoid the over-smoothing by incorporating the positional constraints in the optimization problem (Kazhdan and Hoppe, 2013). Progressive iterative approximation (PIA): PIA is widely used due to its ability to fit data efficiently. PIA method elegantly generates a sequence of curves/surfaces by refining the control points of blending curves/surfaces, and the data points are interpolated by the limit of the sequence. Initially, Qi et al. (1975) and de Boor (1979) discovered the property of PIA for uniform cubic B-spline curve. Later, Lin et al. (2004) showed that non-uniform B-spline curve and surface have PIA property. Moreover, PIA property holds for curves and surfaces with normalized totally positive (NTP) basis (Lin et al., 2005). Also rational B-spline curve and surface possess this property (Shi and Wang, 2006). Lu (2010) proposed a weighted PIA technique to increase the convergence rate of the PIA method. Both the numbers of control points and data points are required to be equal in the classical PIA. With the recent advancement of the big data era, it is infeasible to fit large-scale data points by PIA. To overcome this drawback, Lin and Zhang (2011) developed an extended PIA (EPIA) method, which allows the number of data points to be higher than that of control points. Progressive and iterative approximation for least square fitting (LSPIA) (Deng and Lin, 2014) is another elegant PIA method, which allows the number of data points to be higher than that of control points, and its limit is the

Y.F. Hamza et al. / Computer Aided Geometric Design 77 (2020) 101817

3

least square fitting result to a given data set. For more details, a comprehensive overview of PIA was provided in Lin et al. (2018b). As stated above, some PIA methods have been developed for both parametric curve and surface fitting, as well as subdivision curve and surface fitting. However, they have never been used in implicit curve and surface reconstruction. 2. Definitions and preliminaries In this section, the statement of the problems and the definition of implicit B-spline curve and surface are given. Specifically, the following problems will be handled in this paper. 2.1. Implicit curve reconstruction problem Given a collection of unordered data points in the two-dimensional space,

{ p i = (xi , y i ), i = 1, 2, ..., n},

(1)

with a set of associated oriented unit normals {ni , i = 1, 2, ..., n}, find a function f (x, y ) so that the zero-level sets of f (x, y ), i.e., f (x, y ) = 0, fit the unordered point set (1). Let f (x, y ) be a bivariate tensor-product B-spline function of bi-degree (d1 , d2 ) defined over some domain  ⊆ R2 (Liu et al., 2017; Yang et al., 2005):

f (x, y ) =

Nu  Nv 

C i j B i (x) B j ( y ),

(2)

i =1 j =1

where C i j are the control coefficients and B i (x), B j ( y ) are the B-spline basis functions with some given knot sequences. We consider the bi-cubic B-spline functions, i.e. d1 = d2 = 3 with uniform knot sequences. The implicit B-spline curve reconstructed by fitting the data point set (1) is given by

z f = {(x, y ) ∈  ⊆ R2 : f (x, y ) = 0}.

(3)

Actually, the implicit curve is reconstructed by minimizing the sum of the squared algebraic distances (Jüttler and Felis, 2002), i.e.,

minimize E (C 11 , C 12 , · · · , C N u , N v ) =

n 

f 2 ( p i ),

(4)

i =1

where the sum E is a quadratic form of the control coefficients C i j , i = 1, 2, · · · , N u , j = 1, 2, · · · , N v of the tensor product B-spline function (2). However, in the implicit curve and surface reconstruction problem, the minimization problem (4) is usually underdetermined, i.e., the number of unknowns is larger than that of the data points, thus leading to extra zero-level sets in the reconstruction results. To make the undetermined problem (4) determined, and to eliminate the extra zero-level sets, global regularization terms are usually added in the minimization problem (4) (Jüttler and Felis, 2002; Rouhani and Sappa, 2011; Rouhani et al., 2014; Yang et al., 2005). However, the addition of the regularization terms seriously affects the efficiency of implicit curve and surface reconstruction. In this paper, we develop the I-PIA method for implicit curve and surface reconstruction, which naturally solves the minimization problem with regularization terms (refer to Eq. (21)), and significantly improves the efficiency of implicit curve and surface reconstruction, while eliminating the extra zero-level sets. 2.2. Implicit surface reconstruction problem Given a collection of unordered data points in the three-dimensional space,

{ p i = (xi , y i , zi ), i = 1, 2, ..., n},

(5)

with a set of associated oriented unit normals {n i , i = 1, 2, ..., n}, find a function f (x, y , z) so that the zero-level sets of f (x, y , z), i.e., f (x, y , z) = 0, fit the unordered point set (5). Define the trivariate tensor product B-spline function f (x, y , z) on some domain  ⊆ R3 as (Liu et al., 2017; Yang et al., 2005):

f (x, y , z) =

Nu  Nv  Nw  i =1 j =1 k =1

C i jk B i (x) B j ( y ) B k ( z),

(6)

4

Y.F. Hamza et al. / Computer Aided Geometric Design 77 (2020) 101817

where C i jk are the control coefficients, and B i (x), B j ( y ), B k ( z) are the B-spline basis functions. In our implementation, the uniform tri-cubic B-spline functions are employed. Analogously, the implicit B-spline surface reconstructed by fitting the data point set (5) is,

z f = {(x, y , z) ∈  ⊆ R3 : f (x, y , z) = 0}.

(7)

And, it is reconstructed by minimizing the sum of the squared algebraic distances (Jüttler and Felis, 2002), i.e.,

minimize E (C 111 , C 112 , · · · , C N u , N v , N w ) =

n 

f 2 ( p i ),

(8)

i =1

where the sum E is a quadratic form of the control coefficients C i jk , i = 0, 1, · · · , N u , j = 0, 1, · · · , N v , k = 1, 2, · · · , N w of the tensor product B-spline function (6). 3. Implicit progressive iterative approximation In this section, we will develop the I-PIA iteration method for the implicit curve and surface reconstruction, and prove its convergence. 3.1. I-PIA for implicit curve reconstruction Given an unordered data point set { p i = (xi , y i ), i = 1, 2, ..., n} (1), with a set of associated oriented unit normals {ni , i = 1, 2, ..., n}, we want to reconstruct an implicit curve from the data point set. To avoid the trivial solution f = 0, we need to add extra offset points { pl = (xl , yl ), l = n + 1, n + 2, ..., 2n} to the data point set (1). The offset points are generated along the normal vector n at a small distance σ (Carr et al., 2001; Fasshauer, 2007; Rouhani et al., 2014), i.e.,

pl = p i + σ ni , Let

l = n + i,

i = 1, 2, ..., n.

 be the value of the implicit function at the offset points, i.e., f ( pl ) =  ,

l = n + 1, n + 2, ..., 2n.

Define the initial implicit B-spline function as follows:

f (0) (x, y ) =

Nu  Nv 

( 0)

C i j B i (x) B j ( y ),

(9)

i =1 j =1

(0)

where the initial control coefficients are taken as C i j = 0. (0)

Let δk , k = 1, 2, · · · , 2n be the difference vectors for the data points, which can be calculated as, ( 0)

δk = 0 − f (0) (xk , yk ), ( 0)

δl

=− f

( 0)

(xl , yl ),

k = 1, 2, ..., n, l = n + 1, n + 2, ..., 2n.

(0)

Moreover, let i j , i = 1, 2, ..., N u , j = 1, 2, ..., N v be the difference vectors for the control coefficients, defined as,

(i j0) =

2n 

( 0)

B i (xk ) B j ( yk )δk

k =1

=



( 0)

B i (xk ) B j ( yk )δk ,

k∈ I i j

where I i j is the index set of the data points pk = (xk , yk ) which satisfy B i (xk ) B j ( yk ) = 0. The new coefficients are obtained by, (1 )

( 0)

( 0)

C i j = C i j + μi j ,

where,

μ is a weight, and its selection method is explained in Remark 3.1. Thus, the new implicit B-spline curve is, f (1) (x, y ) =

Nu  Nv  i =1 j =1

(1 )

C i j B i (x) B j ( y ).

Y.F. Hamza et al. / Computer Aided Geometric Design 77 (2020) 101817

Likewise, suppose that we have obtained the (α )

δk

= 0 − f (α ) (xk , yk ),

5

α -th implicit curves f (α ) (x, y ) after the α -th iteration, and let,

k = 1, 2, ..., n,

(α )

=  − f (α ) (xl , yl ), l = n + 1, n + 2, ..., 2n,  (α ) (i jα ) = B i (xk ) B j ( yk )δk , δl

k∈ I i j

( α +1 )

Cij

= C i(jα ) + μ(i jα ) ,

f (α +1) (x, y ) =

Nu  Nv 

(10) ( α +1 )

Cij

B i (x) B j ( y ).

(11)

i =1 j =1

From the above iterative procedure, a series of implicit B-spline functions { f (α ) (x, y ), α = 0, 1, 2, ...} can be generated. Now, arrange the control coefficients C i j , i = 1, 2, · · · , N u , j = 1, 2, · · · , N v in the lexicographic order, i.e., (α )

(α )

(α )

(α )

C (α ) = [C 11 , C 12 , ..., C 1N v , ..., C N u , N v ]T ,

(12)

and let,

b = [b1 , b2 , · · · , b2n ]T = [0, 0, · · · , 0,  ,  , · · · ,  ]T .





 



n



n

According to (10), we have ( α +1 )

Cij

= C i(jα ) + μ

2n 

(α )

B i (xk ) B j ( yk )δk ,

k =1

(α )

= Cij + μ

2n 

⎡ B i (xk ) B j ( yk ) ⎣bk −

Nu  Nv 

⎤ C i j B i (xk ) B j ( yk )⎦ . (α )

i =1 j =1

k =1

Then, the I-PIA (10) for implicit curve reconstruction can be represented in the matrix form,





C ( α +1 ) = C ( α ) + μ B T b − B C ( α ) ,



= I − μ B T B C (α ) + μ B T b ,

α = 0, 1, 2, ...

(13)

where, B is the collocation matrix of the basis functions (arranged in the lexicographical order),

{ B 1 (x) B 1 ( y ), B 1 (x) B 2 ( y ), · · · , B 1 (x) B N v ( y ), · · · , B N u (x) B 1 ( y ), · · · , B N u (x) B N v ( y )} on the point sequence {(xk , yk ), k = 1, 2, · · · , 2n}, i.e.,

⎡

··· ··· B =⎣ ··· ··· ··· B 1 (x2n ) B 1 ( y 2n ) B 1 (x2n ) B 2 ( y 2n ) · · · B 1 (x1 ) B 1 ( y 1 ) ⎢ B 1 (x2 ) B 1 ( y 2 ) ⎢

B 1 (x1 ) B 2 ( y 1 ) B 1 (x2 ) B 2 ( y 2 )

⎤

B N u (x1 ) B N v ( y 1 ) B N u (x2 ) B N v ( y 2 ) ⎥ ⎥

⎦ ··· B N u (x2n ) B N v ( y 2n )

Remark 3.1. For the convergence of the I-PIA iterative method (13), the weight T

T

(14)

.

μ (13) should satisfy 0 < μ <

2

λmax ( B T B )

,

where λmax ( B B ) is the largest eigenvalue of B B (13). Moreover, for the fast convergence of I-PIA, a practical selection of μ in our implementation is (refer to Deng and Lin, 2014),





μ = , where, C =  B T B  = max 2

C

∞

ij



B i (xk ) B j ( yk ).

k

It can be shown that, the I-PIA iterative procedure is a kind of gradient descent method with μ as the step size. The convergence analysis of the I-PIA for the implicit curve reconstruction (13) will be presented in Section 3.3.

6

Y.F. Hamza et al. / Computer Aided Geometric Design 77 (2020) 101817

3.2. I-PIA for implicit surface reconstruction The I-PIA iterative method for the implicit curve reconstruction can be easily extended to the implicit surface reconstruction. In the following, the I-PIA for implicit surface reconstruction is highlighted. Given the unordered data point set { p i = (xi , y i , zi ), i = 1, 2, ..., n} (5), with a set of associated oriented unit normals {ni , i = 1, 2, ..., n}, we want to reconstruct an implicit surface from the data point set. To avoid the trivial solution f = 0, extra offset points { pl = (xl , yl , zl ), l = n + 1, n + 2, ..., 2n} should be added to the data point set (5). The offset points are generated along the normal vector n at a small distance σ (Carr et al., 2001; Fasshauer, 2007; Rouhani et al., 2014), i.e.,

pl = p i + σ ni , Moreover, let

l = n + i,

i = 1, 2, ..., n.

 be the value of the implicit function at the offset points, i.e.,

f ( pl ) =  ,

l = n + 1, n + 2, ..., 2n.

Similar to the implicit curve reconstruction case, suppose that we have obtained the f (α ) (x, y , z) after the α -th iteration, and let, (α )

δr

(α )

δl

= 0 − f (α ) (xr , yr , zr ), =− f

α) (i jk =

2n 

(α )

(xl , yl , zl ),

r = 1, 2, · · · , n , l = n + 1, n + 2, · · · , 2n, (α )

B i (xr ) B j ( y r ) B k ( zr )δr

=

r =1

( α +1 )

C i jk

α -th implicit B-spline surface



(α )

B i (xr ) B j ( y r ) B k ( zr )δr

,

r ∈ I i jk

(α )

(α )

= C i jk + μi jk ,

f (α +1) (x, y , z) =

(15)

Nu  Nv  Nw 

( α +1 )

C i jk

B i (x) B j ( y ) B k ( z),

(16)

i =1 j =1 k =1

where I i jk is the index set of the data points p r = (xr , y r , zr ), which satisfy B i (xr ) B j ( y r ) B k ( zr ) = 0. In this way, a series of implicit B-spline functions { f (α ) (x, y , z), α = 0, 1, 2, ...} is generated. Let,

b = [b1 , b2 , · · · , b2n ]T = [0, 0, · · · , 0,  ,  , · · · ,  ]T .





 

n





n

The I-PIA (15) for implicit surface reconstruction can be represented in the matrix form,





C ( α +1 ) = C ( α ) + μ B T b − B C ( α ) ,



= I − μ B T B C (α ) + μ B T b ,

α = 0, 1, 2, ...

(17)

where, C (α ) is the vector of the control coefficients arranged in the lexicographical order, (α )

(α )

(α )

(α )

C (α ) = [C 111 , C 112 , ..., C 1,1, N w , ..., C N u , N v , N w ]T , and B is the collocation matrix of the basis functions (arranged in the lexicographical order),

{ B 1 (x) B 1 ( y ) B 1 ( z), B 1 (x) B 1 ( y ) B 2 ( z), · · · , B 1 (x) B 1 ( y ) B N v ( z), · · · , B N u (x) B N v ( y ) B 1 ( z), · · · , B N u (x) B N v ( y ) B N w ( z)} on the point sequence {(xk , yk , zk ), k = 1, 2, · · · , 2n}, i.e.,

⎡

··· ··· B =⎣ ··· ··· ··· B 1 (x2n ) B 1 ( y 2n ) B 1 ( z2n ) B 1 (x2n ) B 1 ( y 2n ) B 2 ( z2n ) · · · B 1 (x1 ) B 1 ( y 1 ) B 1 ( z1 ) ⎢ B 1 (x2 ) B 1 ( y 2 ) B 1 ( z2 ) ⎢

B 1 (x1 ) B 1 ( y 1 ) B 2 ( z1 ) B 1 (x2 ) B 1 ( y 2 ) B 2 ( z2 )

⎤

B N u (x1 ) B N v ( y 1 ) B N w ( z1 ) B N u (x2 ) B N v ( y 2 ) B N w ( z2 ) ⎥ ⎥

⎦ ··· B N u (x2n ) B N v ( y 2n ) B N w ( z2n )

.

(18)

Remark 3.2. Similar as Remark 3.1, for the convergence of the I-PIA iterative method (17), the weight μ (17) should satisfy 0 < μ < λ (2B T B ) , where λmax ( B T B ) is the largest eigenvalue of B T B (17). Moreover, a practical selection of μ for the fast max

convergence of I-PIA is (refer to Deng and Lin, 2014),

Y.F. Hamza et al. / Computer Aided Geometric Design 77 (2020) 101817

 

2

 

μ = , where, C =  B T B  = max ∞

C



i jk

7

B i (xl ) B j ( yl ) B k ( zl ).

l

μ is obtained by bounding the spectral norm in terms of the infinity norm.

The value of

The convergence analysis of the I-PIA for the implicit surface reconstruction (17) will be presented in Section 3.3. 3.3. Convergence analysis As shown above, the I-PIA iteration method for implicit curve reconstruction (13) and that for implicit surface reconstruction (17) can be represented in a unified form, i.e.,



C ( α +1 ) = I − μ B T B C ( α ) + μ B T b , The matrices

α = 0, 1, 2, ...

(19)

T

μ B B in Eqs. (13) and (17) holds some common properties:

Property 1: The matrix B T B is positive semi-definite. So, its eigenvalues are nonnegative real numbers. Property 2: The matrix B T B is singular. In the implicit B-spline curve and surface reconstruction, the number of data points (supposed to be n) is less than that of the control coefficients of implicit curve and surface (supposed to be m), i.e., m > n. Because the order of the matrix B is n × m, its rank is at most n. So the rank of the m × m matrix B T B is at most n (n < m). It means that the matrix B T B is singular. Property 3: The eigenvalues of μ B T B satisfy 0 < λ(μ B T B ) < 2. This is because of the choice of the weight μ, as well as Property 1 and 2. Remark 3.3. According to Property 2, the m × m matrix B T B is singular. Suppose the dimension of its zero eigenspace is m0 . Then, the rank of the matrices B T B and μ B T B is the same, i.e.,

rank(μ B T B ) = rank( B T B ) = m − m0 . Moreover, the implicit curve and surface reconstruction problem can be formulated as solving the following normal equation of a least-squares fitting problem,

B T B X = B T b,

(20)

where X is an unknown vector, B and b are the same as in Eq. (19). As pointed out by Property 2, the coefficient matrix B T B is singular. So, if the solution of the least-square fitting system (20) exists, it has infinite solutions, which usually leads to extra zero-level sets. Therefore, to eliminate the extra zero-level sets and get desirable results, it is required to solve the constrained minimization problem,

min  X  E X

s.t .

(21)

B T B X = B T b,

where, · E is the Euclidean norm. In the following, we show that, the I-PIA iterative method (19) converges to the solution of the constrained minimization problem (21).



Theorem 3.4. When the initial values C (0) equals 0 C (0) = 0 , the I-PIA iterative method (19) converges to the solution of the constrained minimization problem (21), i.e., ( B T B )+ B T b, where ( B T B )+ is the Moore-Penrose (M-P) pseudo-inverse of the matrix B T B. Proof. By Remark 3.3, rank(μ B T B ) = rank( B T B ) = m − m0 . Because the matrix B T B is both a normal matrix and a positive semi-definite matrix, its eigen decomposition and singular value decomposition are the same,

B T B = V diag (λ1 , λ2 , ..., λm−m0 , 0, 0, ..., 0) V T ,

(22)

   m0

where V is an orthogonal matrix, and λi > 0, i = 1, 2, · · · , m − m0 are both the eigen values and singular values of the matrix B T B. Then, the M-P inverse of B T B is,

⎛

( B T B )+ = V diag ⎝

⎞

1

1

1

, , ..., , 0, 0, ..., 0⎠ V T . λ1 λ2 λm−m0    m0

Therefore, we have,

8

Y.F. Hamza et al. / Computer Aided Geometric Design 77 (2020) 101817

( B T B )+ ( B T B ) = V diag (1, 1, · · · , 1, 0, 0, ..., 0) V T .       m−m0

m0

Due to Eq. (22), it holds,

μ B T B = V diag (μλ1 , μλ2 , ..., μλm−m0 , 0, 0, ..., 0) V T ,    m0

where μλi , i = 1, 2, · · · , m − m0 are the eigen values of the matrix 1, 2, · · · , m − m0 . Therefore,

lim

α →∞

I − μBT B

α

μ B T B. Based on Property 3, they satisfy 0 < μλi < 2, i =

= lim V diag ((1 − μλ1 )α , (1 − μλ2 )α , · · · , (1 − μλm−m0 )α , 1, 1, · · · , 1) V T    α →∞ m0

= V diag (0, ..., 0, 1, ..., 1) V       m−m0

T

m0

(23)

= I − V V T ( B T B )+ ( B T B ) V V T = I − ( B T B )+ ( B T B ). Note that the linear system B T B X = B T b has solutions, if and only if (James, 1978),

( B T B )( B T B )+ ( B T b) = B T b.

(24)

Then, subtracting ( B T B )+ B T b from both sides of Eq. (19), leads to,

C (α +1) − ( B T B )+ B T b = ( I − μ B T B )C (α ) + μ B T b − ( B T B )+ B T b

= ( I − μ B T B )C (α ) + μ( B T B )( B T B )+ ( B T b) − ( B T B )+ B T b = ( I − μ B T B )C (α ) − ( I − μ B T B )( B T B )+ B T b = ( I − μ B T B )(C (α ) − ( B T B )+ B T b) = ( I − μ B T B )α +1 (C (0) − ( B T B )+ B T b). So, together with Eqs. (23) and (24), we have,

C (∞) − ( B T B )+ B T b = lim ( I − μ B T B )(α +1) (C (0) − ( B T B )+ B T b) α →∞

= ( I − ( B T B )+ ( B T B ))(C (0) − ( B T B )+ B T b) = ( I − ( B T B )+ ( B T B ))C (0) . Therefore,

C (∞) = ( B T B )+ B T b + ( I − ( B T B )+ ( B T B ))C (0) , which are the solutions of the singular linear system B T B X = B T b (note that C (0) can take an arbitrary value). Among them, ( B T B )+ B T b is the one with minimum Euclidean norm. So, when the initial value C (0) = 0, the I-PIA iterative format converges to C ∞ = ( B T B )+ B T b, the solution of the singular linear system B T B X = B T b with the minimum Euclidean norm. It is the solution of the constrained minimization problem (21). 2 As aforementioned, there are infinite solutions to the singular normal equation (20), and each one corresponds to a B-spline function (i.e., the solutions of the equation (20) are the control coefficients of the B-spline functions). Suppose all B-spline functions so generated contain the true zero set, and there exists one B-spline function which contains just the true zero set (without extra zero sets). Then, the solution X ∗ of Eq. (20), generated by the I-PIA iterations with zero initial value, corresponds to the B-spline function without extra zero sets. Otherwise, suppose the B-spline function taking X ∗ as its control coefficient contains extra zero sets, and the control coefficient of the B-spline function without extra zero sets is X t . So, X ∗ and X t are both the solutions of Eq. (20), and  X t  E <  X ∗  E . This is a contradiction, because X ∗ is the solution of the constrained minimization problem (21), i.e., the solution of Eq. (20) with the minimum Euclidean norm, according to Theorem 3.4. Therefore, we have the following corollary. Corollary 3.5. If the implicit B-spline functions which take the solutions of Eq. (20) as control coefficients all contain the true zero sets, and there exists one implicit B-spline function without extra zero sets, then the solution generated by I-PIA iterations with zero initial value is the control coefficients of the implicit B-spline function without extra zero sets.

Y.F. Hamza et al. / Computer Aided Geometric Design 77 (2020) 101817

9

Fig. 1. Iterations in the reconstruction of 2D data sets: First row, flower model; second row, Coons curve model. Blue points are the given data sets, and the red line is the reconstructed curve. From left to right: the 1st, 5th, 10th, and 15th iteration step. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

4. Experiments and discussion Several experiments have been carried out to evaluate the performance of I-PIA. All the experiments were performed in MATLAB on a PC with an Intel-core i7 @ 3.6 GHz processor and 16 GB of RAM. The results and discussions specifically focus on the following areas: effectiveness, robustness to inaccurate distance field, hole filling, noisy data, non-uniform sampling, grid resolution, porous surface (open surface), and fine details. We take the fitting error as error = max δk . The performance of I-PIA, as well as the comparison with the state-of-the-art methods (Huang et al., 2019; Liu et al., 2017), is described in Table 1 and 2. 4.1. Effectiveness As stated above, the extra zero-level sets generated in the implicit curve and surface reconstruction procedure make the reconstruction results challenging to be interpreted, and the elimination of extra zero-level sets is the main problem in designing implicit curve and surface reconstruction method. With the I-PIA developed in this paper, no extra zero-level set appears in the reconstruction procedure. To show the effectiveness of I-PIA in the reconstruction of implicit curves and surfaces without the appearance of extra level sets, we tested our algorithm on planar curves and 3D surfaces. The initial control coefficients are taken as zero. After every iteration, the control coefficients are updated and refined by the difference vectors for the control coefficients, and the resulting implicit curve/surface approach the given data sets closer than the previous implicit curve/surface. Fig. 1 shows the reconstruction process of two 2D data sets. Similarly, Fig. 2 illustrates the reconstruction process of two 3D data sets with sparse sampling. From the results presented in Fig. 1 and Fig. 2, we can see that no extra zero-level set exists in the reconstruction process of the 2D and 3D data sets, even the 3D data sets are sparsely sampled. 4.2. Robustness to inaccurate distance field The auxiliary offset points appended to the data points help orient the surface and avoid the appearance of artifacts in curve and surface reconstruction (Carr et al., 2001; Liu et al., 2017). The I-PIA algorithm is insensitive to the distance values assigned at the offset points, thus robust to the inaccurate distance field. The first row of Fig. 3 illustrates the reconstruction of a 2D data set with synthesis offset points. The function values on the offset points are assigned in a random strategy. That is, the preassigned values are selected as uniformly distributed random numbers in [0.5 − ω, 0.5 + ω] for the outside offset points and [−0.5 − ω, −0.5 + ω] for the inside offset points, respectively, where ω = 0, 0.02, 0.05, 0.1, as demonstrated in Figs. 3(b)-3(e). The results demonstrate that I-PIA can still reconstruct the given data sets in a corrupted distance field. Moreover, given a data set, we generate two sets of offset points inside and outside the data set, respectively. However, the distance values are assigned in the reverse order, i.e., larger distance values are given to closer offset points, and smaller distance values are given to farther offset points. As illustrated in Figs. 3(f) and 3(h), from inside to outside, the distance values of the four sets of offset points are −0.1, −0.2, 0.2, 0.1. The implicit curves are robustly reconstructed using I-PIA, and demonstrated in Figs. 3(g) and 3(i).

10

Y.F. Hamza et al. / Computer Aided Geometric Design 77 (2020) 101817

Fig. 2. Iterations in the reconstruction of 3D data sets with sparse sampling, Torus and double-Torus. The red points are the given data sets. From left to right: the 1st, 5th, 10th, and 15th iteration step(s).

Fig. 3. Reconstruction of 2D data sets with inaccurate distance fields. First row: (a) the input data of dolphin (blue) and offset points (outside offset in green and inside offset in magenta). The preassigned function values are selected as uniformly distributed random numbers in [±0.5 − ω, ±0.5 + ω] for the outside and inside offsets, respectively, (b)–(e): ω = 0, 0.02, 0.05, 0.1. Second row: (f) and (h): the input data of butterfly and dolphin (blue), respectively, and offset points (outside offset in green and inside offset in magenta). Larger values are assigned to closer offset points and vise versa. From inside to outside, the distance values of the four sets of offset points are −0.1, −0.2, 0.2, 0.1. (g) and (i): the reconstructed curves from (f) and (h), respectively.

4.3. Hole filling Holes and missing data often occur in the point cloud generated by scanning devices. Fig. 4 (a)–(d) depict the reconstructed surfaces of human face (Wullbert, 2019), bunny, elephant, and fertility, respectively, in which some portions of the data points are missing or artificially removed to create holes and gaps. Specifically, while the missing data are notable

Y.F. Hamza et al. / Computer Aided Geometric Design 77 (2020) 101817

11

Fig. 4. Hole and gap filling. (a, e) Human face; (b, f) Bunny; (c, g) Elephant; (d, h) Fertility. First row: reconstructed surfaces. Second row: the reconstructed surfaces with the data points superimposed.

Fig. 5. Surface reconstruction from noisy data and non-uniform sampling. (a) Noisy data. (b) Reconstructed surface of the noisy data by I-PIA. (c) Nonuniform sampling: 98% of the middle part of the torus model are removed. (d) Reconstructed surface of non-uniform sampling by I-PIA with data point superimposed.

around the eyebrows and naris of the human face (Fig. 4 (e)), and the bottom of the bunny (Fig. 4 (f)), the remaining holes in the middle of the bunny, elephant and fertility models are artificially created. In Fig. 4, the reconstructed surfaces are illustrated in the first row, and the second row shows the reconstructed surfaces with the data points superimposed. As can be seen, I-PIA successfully fills the holes and gaps. Thus, I-PIA is robust to holes and missing data. 4.4. Noisy data, non-uniform sampling, and grid resolution In this section, we demonstrate the capability of I-PIA in handling noisy data and non-uniform sampling. Fig. 5(a) displays the noisy data. In Fig. 5(c), the middle part of the torus model is down-sampled, and 98% of the original data points are removed. As can be seen in Fig. 5(b, d), I-PIA is robust to noisy data and non-uniform sampling. Moreover, the features and details of the reconstructed surface by I-PIA depend on the resolution of the grid. Fig. 6 shows the reconstructed surfaces from the point clouds of the model armadillo with different densities, where Fig. 6(a) is the original point cloud, and Fig. 6(e) is the low-density point cloud with 30% points of the original point cloud. From left to right of Fig. 6, they are the input point clouds, and reconstructed surfaces based on coarse to fine grids (50 × 50 × 50, 100 × 100 × 100, 190 × 190 × 190), respectively. As the grid is made denser and denser, the details and features of the armadillo become more and more pronounced. As can be seen, though the sampling densities of the two point clouds in Figs. 6(a) and 6(e) are different, the implicit surfaces reconstructed from them with the same grid resolution are nearly the same (Fig. 6).

12

Y.F. Hamza et al. / Computer Aided Geometric Design 77 (2020) 101817

Fig. 6. Surfaces reconstructed from the point clouds with different sampling densities of the armadillo model. (a) and (e): Input point clouds with different sampling densities. (b) and (f): Reconstructed surfaces based on a coarse grid (50 × 50 × 50). (c) and (g): Reconstructed surfaces based on a grid (100 × 100 × 100). (d) and (h): Reconstructed surfaces based on a fine grid (190 × 190 × 190).

Fig. 7. Surface reconstruction of Gyroid surface (a triply periodic minimal surface). (a) Gyroid surface. (b) The point cloud sampled from the Gyroid surface. (c) Reconstructed surface by I-PIA. (d) Reconstructed surface with data points superimposed.

4.5. Open surface (porous surface) Gyroid surface is a triply periodic minimal surface, which has zero mean curvature. An example of a Gyroid surface is given in Fig. 7(a). To show the ability of I-PIA in processing the surface reconstruction for the open boundary, we sampled the data points in Fig. 7(b) from Fig. 7(a) (gyroid surface). Fig. 7(c) depicts the reconstructed surface from the sampled data points (Fig. 7(b)). The reconstructed surface with the given data points superimposed is displayed in Fig. 7(d). It can be seen that the implicit B-spline surface reconstructed by I-PIA well conforms to the given point cloud. 4.6. Fine details To show the level of quality of the surfaces reconstructed by our algorithm, we enlarge some parts of BU and Buddha in Fig. 8. The features of the models are correctly recovered using I-PIA. 4.7. Performance and discussion Table 1 summarizes the running time performance of I-PIA on different planar curves. Table 2 shows the comparison of running time performance of I-PIA and the state-of-the-art methods in Huang et al. (2019); Liu et al. (2017). On the one hand, the computational cost of implicit surface reconstruction with total variation regularization (ISRTV) (Liu et al., 2017) is higher than that of I-PIA, because in each iteration of the ISRTV method (Liu et al., 2017), a linear system needs to be solved, while I-PIA avoids solving the linear system. The surface reconstructions of Armadillo, Buddha, and BU are not completed within 24 hours by the ISRTV method (Liu et al., 2017). In summary, the implicit curve and surface reconstruction time

Y.F. Hamza et al. / Computer Aided Geometric Design 77 (2020) 101817

13

Fig. 8. Surface reconstruction of BU and Buddha with sectional enlargements.

Table 1 Performance of our method on different planar curves. Model

Number of data point

Grid size

Maximum Error

Time (seconds)

Flower Coons curve Dolphin Butterfly

509 780 371 508

30 × 30 30 × 30 30 × 30 28 × 28

2.924e–3 1.2524e–3 3.2457e–3 2.8881e–3

1.88 2.05 1.67 1.79

Table 2 Performance of I-PIA and the methods in Huang et al. (2019) and Liu et al. (2017) on 3D point clouds (all timings are measured in seconds). Model

Double-torus Max Plank Torus Gyroid Elephant Bunny Fertility Armadillo Buddha Bu

Data point

766 4000 10100 36627 52100 35947 241607 51892 80000 70000

Grid size 30 × 30 × 30 50 × 50 × 50 30 × 30 × 30 50 × 50 × 50 50 × 50 × 50 70 × 70 × 70 70 × 70 × 70 190 × 190 × 190 150 × 150 × 150 160 × 160 × 160

Max. Error

1.0574e–3 9.5715e–3 3.0412e–6 4.8708e–6 9.6839e–4 9.5983e–2 9.5989e–5 1.9831e–2 7.4594e–2 1.1899e–3

Method in Huang et al. (2019)

Liu et al. (2017)

I-PIA

Time

Time

Time

333.78 36364.6 out of memory out of memory out of memory out of memory out of memory out of memory out of memory out of memory

90.07 5525 189.13 11149.28 2773.65 26428.69 25298.49 >24hr >24hr >24hr

3.61 4.75 8.22 34.26 44.55 88.24 902.35 2732.32 1835.73 2028.01

by I-PIA is reduced at least one to three orders of magnitude, compared with the state-of-the-art method ISRTV (Liu et al., 2017). On the other hand, VIPSS method (Huang et al., 2019) is suited for sparse data points due to its high computational complexity and high space requirement. Thus, the sparse data points are provided (Double-torus and Max Plank) to compare the I-PIA method with the VIPSS method. Referring to Table 2, when the number of points is larger than 10000, the program of the VIPSS method is out of memory, where the code downloaded from https://github.com/adshhzy/VIPSS is employed. Moreover, Fig. 9 illustrates the reconstructed surfaces of Max Plank using I-PIA, VIPSS (Huang et al., 2019), and ISRTV (Liu et al., 2017), respectively. Twenty iterations were performed in the reconstruction process of the Max Plank using I-PIA and ISRTV, and the time consumed for the reconstruction by I-PIA and ISRTV are 4.75 seconds and 5525 seconds, respectively (refer to Table 2). Furthermore, the maximum fitting error for I-PIA and ISRTV are 9.5715e–3 and 1.0230e–1, respectively. With I-PIA, the construction of the iterative matrix B (Eq. (14) (18)) consumes more than 90% of the running time. For example, 88.24 seconds were spent in the reconstruction of the bunny model, but the construction of the iterative matrix B cost 81.92 seconds. In general, there are two approaches for accelerating the computation. In the first approach, tabulating the univariate B-spline values can speed up the computations. In the second approach, the calculation of matrix B can be accelerated by parallel computing in virtue of the local property of the B-spline basis. The main limitation of I-PIA is the weak sharp feature preservation capability. Feature preservation is a traditionally difficult problem in the implicit curve and surface reconstruction. In future work, it is worth developing the sharp feature preservation techniques in the I-PIA iterations, and the parallelization methods for acceleration.

14

Y.F. Hamza et al. / Computer Aided Geometric Design 77 (2020) 101817

Fig. 9. Surface reconstruction of Max Plank: (a) The Max Plank data point. (b) Reconstructed surface by I-PIA method. (c) Reconstructed surface by VIPSS method. (d) Reconstructed surface by method in Liu et al. (2017).

5. Conclusions This paper presents a novel approach for implicit curve and surface reconstruction, i.e., I-PIA. I-PIA naturally solves the minimization problem with regularization terms in its iterations. Thus, it not only avoids the spurious sheets and artifacts but also reduces the computational cost effectively. Several kinds of experiments demonstrate that I-PIA is robust to inaccurate distance fields, holes and gaps, non-uniform sampling and noisy data, and thus produces high-quality reconstruction results. Declaration of competing interest The authors wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome. Acknowledgements This work is supported by the National Natural Science Foundation of China under Grant Nos. 61872316, 61932018, and the National Key R&D Plan of China under Grant No.2016YFB1001501. References Blinn, J.F., 1982. A generalization of algebraic surface drawing. ACM Trans. Graph. 1 (3), 235–256. Carr, J.C., Beatson, R.K., Cherrie, J.B., Mitchell, T.J., Fright, W.R., McCallum, B.C., Evans, T.R., 2001. Reconstruction and representation of 3D objects with radial basis functions. In: Proceedings of the 28th Annual Conference on Computer Graphics and Interactive Techniques. SIGGRAPH ’01, pp. 67–76. Chen, Z., Luo, X., Tan, L., Ye, B., Chen, J., 2008. Progressive interpolation based on Catmull-Clark subdivision surfaces. Comput. Graph. Forum 27 (7), 1823–1827. Cheng, F.-H.F., Fan, F.-T., Lai, S.-H., Huang, C.-L., Wang, J.-X., Yong, J.-H., 2009. Loop subdivision surface based progressive interpolation. J. Comput. Sci. Technol. 24 (1), 39–46. Curless, B., Levoy, M., 1996. A volumetric method for building complex models from range images. In: Proceedings of the 23rd Annual Conference on Computer Graphics and Interactive Techniques. ACM, pp. 303–312. de Boor, C., 1979. How does Agee’s smoothing method work. In: Proceedings of the 1979 Army Numerical Analysis and Computers Conference. ARO Report. 79-3. Deng, C., Lin, H., 2014. Progressive and iterative approximation for least squares B-spline curve and surface fitting. Comput. Aided Des. 47, 32–44. Deng, C., Ma, W., 2012. Weighted progressive interpolation of Loop subdivision surfaces. Comput. Aided Des. 44 (5), 424–431. Fasshauer, G.E., 2007. Meshfree Approximation Methods with MATLAB, vol. 6. World Scientific. Hoppe, H., DeRose, T., Duchamp, T., McDonald, J., Stuetzle, W., 1992. Surface reconstruction from unorganized points. Comput. Graph. 26 (2), 71–78. Huang, Z., Carr, N., Ju, T., 2019. Variational implicit point set surfaces. ACM Trans. Graph. 38 (4), 124:1–124:13. James, M., 1978. The generalised inverse. Math. Gaz. 62 (420), 109–114. Jüttler, B., Felis, A., 2002. Least-squares fitting of algebraic spline surfaces. Adv. Comput. Math. 17 (1–2), 135–152. Kazhdan, M., Bolitho, M., Hoppe, H., 2006. Poisson surface reconstruction. In: Proceedings of the Fourth Eurographics Symposium on Geometry Processing, vol. 7, pp. 61–70. Kazhdan, M., Hoppe, H., 2013. Screened Poisson surface reconstruction. ACM Trans. Graph. 32 (3), 29. Kojekine, N., Hagiwara, I., Savchenko, V., 2003. Software tools using CSRBFs for processing scattered data. Comput. Graph. 27 (2), 311–319. Lin, H., 2010. Local progressive-iterative approximation format for blending curves and patches. Comput. Aided Geom. Des. 27 (4), 322–339. Lin, H.-W., Bao, H.-J., Wang, G.-J., 2005. Totally positive bases and progressive iteration approximation. Comput. Math. Appl. 50 (3–4), 575–586. Lin, H., Cao, Q., Zhang, X., 2018a. The convergence of least-squares progressive iterative approximation for singular least-squares fitting system. J. Syst. Sci. Complex. 31 (6), 1618–1632. Lin, H., Maekawa, T., Deng, C., 2018b. Survey on geometric iterative methods and their applications. Comput. Aided Des. 95, 40–51. Lin, H., Wang, G., Dong, C., 2004. Constructing iterative non-uniform B-spline curve and surface to fit data points. Sci. China, Ser. F, Inf. Sci. 47 (3), 315–331. Lin, H., Zhang, Z., 2011. An extended iterative format for the progressive-iteration approximation. Comput. Graph. 35 (5), 967–975. Liu, Y., Song, Y., Yang, Z., Deng, J., 2017. Implicit surface reconstruction with total variation regularization. Comput. Aided Geom. Des. 52, 135–153. Lu, L., 2010. Weighted progressive iteration approximation and convergence analysis. Comput. Aided Geom. Des. 27 (2), 129–137.

Y.F. Hamza et al. / Computer Aided Geometric Design 77 (2020) 101817

15

Morse, B.S., Yoo, T.S., Rheingans, P., Chen, D.T., Subramanian, K.R., 2005. Interpolating implicit surfaces from scattered surface data using compactly supported radial basis functions. In: ACM SIGGRAPH 2005 Courses. ACM, p. 78. Muraki, S., 1991. Volumetric shape description of range data using “blobby model”. Comput. Graph. 25 (4), 227–235. Ohtake, Y., Belyaev, A., Alexa, M., Turk, G., Seidel, H.-P., 2003. Multi-level partition of unity implicits. ACM Trans. Graph. 22 (3), 463–470. Ohtake, Y., Belyaev, A., Seidel, H.-P., 2005. Multi-scale and adaptive CS-RBFS for shape reconstruction from clouds of points. In: Advances in Multiresolution for Geometric Modelling. Springer, pp. 143–154. Pan, M., Tong, W., Chen, F., 2016. Compact implicit surface reconstruction via low-rank tensor approximation. Comput. Aided Des. 78, 158–167. Qi, D., Tian, Z., Zhang, Y., Feng, J., 1975. The method of numeric Polish in curve fitting. Acta Math. Sin. 18 (3), 173–184. Rouhani, M., Sappa, A.D., 2011. Implicit B-spline fitting using the 3L algorithm. In: 2011 18th IEEE International Conference on Image Processing. IEEE, pp. 893–896. Rouhani, M., Sappa, A.D., Boyer, E., 2014. Implicit B-spline surface reconstruction. IEEE Trans. Image Process. 24 (1), 22–32. Shi, L., Wang, R., 2006. An iterative algorithm of NURBS interpolation and approximation. J. Math. Res. Expo. 26 (4), 735–743. Turk, G., O’brien, J.F., 1999. Variational implicit surfaces. Technical Report GIT-GUV-99-15. Georgia Institute of Technology. Wang, J., Yang, Z., Jin, L., Deng, J., Chen, F., 2011. Parallel and adaptive surface reconstruction based on implicit PHT-splines. Comput. Aided Geom. Des. 28 (8), 463–474. Wullbert, 2019. Grabcad. https://grabcad.com/library/jaffar-face-1. (Accessed 24 November 2019). Yang, Z., Deng, J., Chen, F., 2005. Fitting unorganized point clouds with active implicit B-spline curves. Vis. Comput. 21 (8–10), 831–839.