03964.com

文档资料库 文档搜索专家

文档资料库 文档搜索专家

Numerical Algorithms (2005) 39: 289–305

? Springer 2005

Fast surface reconstruction and hole ?lling using positive de?nite radial basis functions

G. Casciola, D. Lazzaro, L.B. Montefusco and S. Morigi

Department of Mathematics, University of Bologna, P.zza di Porta San Donato 5, 40127 Bologna, Italy E-mail: {casciola;montelau;morigi}@dm.unibo.it, lazzaro@csr.unibo.it

Received 29 September 2003; accepted 1 June 2004

Surface reconstruction from large unorganized data sets is very challenging, especially if the data present undesired holes. This is usually the case when the data come from laser scanner 3D acquisitions or if they represent damaged objects to be restored. An attractive ?eld of research focuses on situations in which these holes are too geometrically and topologically complex to ?ll using triangulation algorithms. In this work a local approach to surface reconstruction from point-clouds based on positive de?nite Radial Basis Functions (RBF) is presented that progressively ?lls the holes by expanding the neighbouring information. The method is based on the algorithm introduced in [7] which has been successfully tested for the smooth multivariate interpolation of large scattered data sets. The local nature of the algorithm allows for real time handling of large amounts of data, since the computation is limited to suitable small areas, thus avoiding the critical ef?ciency problem involved in RBF multivariate interpolation. Several tests on simulated and real data sets demonstrate the ef?ciency and the quality of the reconstructions obtained using the proposed algorithm. Keywords: radial basis function, surface reconstruction, local interpolation, hole ?lling AMS subject classi?cation: 65D17, 65D05, 65Y20

1.

Introduction

The digitalization and reconstruction of 3D shapes has numerous applications in areas that include manufacturing, virtual simulation, medicine, entertainment, consumer marketing and archaeology. Methods to digitize and reconstruct the shapes of complex objects have evolved rapidly in recent years. The speed and accuracy of digitizing technologies allows us to take detailed measurements of the objects, but in order to capture the complete shape of an object, many thousands of samples must be acquired. The resulting mass of data requires algorithms that can ef?ciently and reliably generate computer models from these samples. The reconstruction of a 3D model is usually obtained from scattered data points sampled from the surface of the physical 3D object. The acquisition of the data is realized through 3D scanning systems which are able to capture a dense and accurate

290

G. Casciola et al. / Fast surface reconstruction and hole ?lling

sampling usually organized into range images, i.e., sets of distances from the sensor to the object being scanned. Each range image needs to be ?ltered to remove artifacts and noise, and to ?ll the holes. Using laser light scanners, in fact, gives rise to range images which can be noisy and affected by artifacts, depending on their material and on the presence of partially re?ecting surfaces. Moreover, they can present simple or topologically complex holes. Most of them can be naturally ?lled by the overlapping range images process, since the most fundamental cause of holes is occlusion – recesses too deep to be observed using a particular angle. However, holes can also be caused by low re?ectance, constraints on scanner placement, or simply missing views. This is frequently observed in the scanning of works of art, which have a lot of self-occlusions and details. Works of art also impose signi?cant restrictions on the scanner placement. Moreover, the holes can correspond to missing or damaged parts of the object when the scanning process is used, e.g., in a computer aided restoration session. Surface reconstruction from dense range data has been an active area of research for several years. The strategies have proceeded along two basic directions: reconstruction from unorganized points, and reconstruction that preserves the underlying structure of the acquired data. These two strategies can be further subdivided, according to whether they operate by reconstructing piecewise linear or smoother surfaces, or by reconstructing an implicit function. Examples of implicit surface reconstruction include the method of Hoppe [6], Bajaj [1] and, more recently, Beatson et al. [2], which use RBF to ?t a signed distance function followed by an isosurface extraction. Most of these approaches are not particularly concerned with the hole ?lling problem. In fact, they lead to a natural ?lling of the holes, which works well for simple holes in nearly ?at surfaces, but not for convoluted geometric holes. In many applications it is important to be able to reconstruct this missing information for post-processing, as well as for presentation in a more plausible way. One widely used approach to hole ?lling consists in triangulating each connected component of the surface boundary, thereby ?lling the holes with a patch. Recently, in the pioneering work [3], the authors specifically address situations in which the holes are too topologically complex to ?ll using triangulation algorithms, and they propose solving the hole ?lling problem via isotropic diffusion of volumetric data (i.e., iterative Gaussian convolution of some distance function to the known data). The approach proposed by these authors produces watertight digital models, but it does not seem to give intuitive answers to the hole ?lling tests presented. We refer to this paper for an excellent and detailed description of the nature of the holes in scanning statues, and for a review of the literature on this subject. In this work we consider a method which is designed for multivariate interpolation (approximation) of general data sets. It is based on an algorithm presented in an earlier paper [7] which makes local use of radial basis function interpolants. This local approach allows for real time handling of large amounts of data, since the computation is limited to suitable small areas, thus avoiding the critical ef?ciency problem involved in RBF multivariate interpolation (see, e.g., [10], and references herein). Here we adapt the algorithm to the speci?c application of surface reconstruction from 3D scanning range data by including an ef?cient hole ?lling procedure, which progressively ?lls the holes

G. Casciola et al. / Fast surface reconstruction and hole ?lling

291

by expanding the neighbouring information. The locality of the proposed method allows us to ?ll in a selected way only desiderated areas, without forcing a convex watertight model which can be obtained using a unique interpolant function. The paper is organized as follows. In section 2 we discuss some theoretical results on interpolation with RBF, which will be assumptions to derive a fast and robust local reconstruction method given in section 3. In section 4, we present the new hole ?lling method and we discuss some relating theoretical aspects. The corresponding algorithm, integrated in the reconstruction procedure, is described in section 5. Results follow in section 6, where the ef?ciency and the capabilities of our approach are demonstrated by considering real data sets. 2. Interpolation with positive de?nite radial basis functions

The interpolation problem can be stated as follows. Given a set of distinct nodes ? Rd , and a set of function values = {fi }N X = {xi }N i =1 ? i =1 ? R, the radial basis function interpolant s : Rd → R on X is given by

N

s(x) =

j =1

cj ? x ? xj

2

,

(1)

where the function ? : [0, ∞) → R is called radial basis function, and the points xj are referred to as centers of the radial basis functions. The coef?cients c1 , c2 , . . . , cN are determined by the interpolation conditions s(xi ) = fi , which lead to the linear system AX,? c = f (3) with AX,? = (?( xi ? xj 2 ))1 i,j N . In order that the system (3) be solvable, the interpolation matrix AX,? must be nonsingular. In this work, we require that the matrix AX,? is positive de?nite. In this case, following, e.g., the de?nition in [12], the function ?( · ) is said to be positive de?nite in Rd , which we abbreviate by ?( · ) ∈ PDd . Well known PDd radial basis functions are the inverse multiquadrics. However, in this case, the interpolation matrix in (3) is non-sparse, and even modestly large data sets can lead to very ill-conditioned linear systems. Furthermore, the evaluation of the interpolant function (1) can become computationally expensive. In order to overcome these problems, Wendland [13] and Whu [14] have proposed compactly supported radial basis functions, which are positive de?nite on Rd according to their given order of smoothness. Although they represent an effort in the direction of localization, they still suffer from the inherent inability of the radial basis functions to interpolate very large data sets in a numerically stable and accurate way. In fact, the reconstruction may be poor, since it strongly depends on the radius of the support of these functions, which has to be scaled to adapt it to the different densities of the scattered data set. Moreover, i = 1, . . . , N, (2)

292

G. Casciola et al. / Fast surface reconstruction and hole ?lling

improving the smoothness of ? improves the reproduction quality but, at the same time, also blows up the condition of AX,? . For the purpose of this work, we considered inverse multiquadrics and compactly supported Wendland’s functions, thus in the following we investigate the critical aspects of stability and accuracy for these classes of functions in greater details. For this aim, we ?rst introduce the following two measures which characterize the density of the data set X : the separation distance qX , which measures the closest pair of points in X , and the ?ll distance hX , which represents the radius of the largest inner empty sphere, given by qX = 1 min 2 1 i<j

N

xi ? xj

2,

hX = sup min

x∈

1 i N

x ? xi 2 .

If we denote by qX (x) and hX (x) the previous quantities relative to the set X ∪ {x}, for each x ∈ , we get the following relation: qX (x)

1 j

min

N

x ? xj

2

hX (x).

The measures qX (x) and hX (x) are critical to the study of stability and accuracy of radial basis function interpolants. Concerning the conditions for the numerical stability of system (3), these can easily be derived from standard theory of numerical stability. In fact, the sensitivity of the solution c with respect to variations in the data is given by the condition number, which, in our case of symmetric positive de?nitive matrix, is given by the ratio /λ between the biggest and the smallest matrix eigenvalue. Thus the stability analysis consists in ?nding good bounds for and λ. In particular, it is necessary to ?nd lower bounds on λ that are as tight as possible. This depends on the separation distance of the set X . On the other hand, if the data are values of a function f belonging to a function space F, the accuracy of the reconstruction is evaluated by the error of the interpolant de?ned by the solution of (3) with fi = f (xi ). When F is de?ned via ? itself in a natural way [11], its characteristics depend on the d-variate generalized Fourier transform φ of ? . If φ decays at least algebraically at in?nity, then F contains a certain Sobolev space, if φ decays exponentially, then F consists of C ∞ functions. In any case, the space F carries a speci?c seminorm | · |F and the interpolation error can be bounded by f (x) ? s(x) |f |F P (x), (4)

where the power function P (x) is the norm of the error functional on F, evaluated at x. Note that P depends on x, X, ?, and F, but not on f . Theoretical and numerical results have revealed a strict connection between the error and the sensitivity described by P (x) and λ, respectively. Schaback proved in [11] that the accuracy and the condition number of the interpolation matrix cannot both be kept small. This effect is a sort of uncertainty principle, and it takes the form P 2 (x) λ(x), (5)

where λ(x) is the smallest eigenvalue of the matrix AX∪{x},? .

G. Casciola et al. / Fast surface reconstruction and hole ?lling

293

Let us recall some well known results for a better comprehension of the choices we made in the design of our RBF interpolation method. According to [11] the known upper and lower bounds for P (x) and λ(x) take the forms P (x)2 F hX (x) and λ(x) G qX (x) , respectively (6)

where F , G : R+ → R+ are continuous and decreasing functions for small arguments. We refer to [11] for a list of known examples of functions F and G related to some special RBF. Considering the uncertainty principle (5), we then obtain the two-sided bounds G qX (x) λ(x) P 2 (x) F hX (x) . (7)

These bounds can be made more useful for practical applications if conditions can be given under which hX (x) can be bounded from above, in terms of qX (x). This is the case if we deal with a set of centers X ? which is quasi-uniform, i.e., if there are ρ > 0 and h > 0, such that 1. hX (x) 2. qX = h, 1 min 2 1 i<j

N

xi ? xj

2

h · ρ.

(8)

In this case we call h the density and ρ the uniformity √of X with respect to . For 2, while a triangular grid has example, a regular square grid has uniformity ρ = 1 / √ ρ = 3/3 [4]. As proved in [11], for quasi-uniform sets of centers with density h and quasiuniformity ρ , we have two-sided bounds (9) P 2 (x) F (h) √ for all x ∈ with q(x) hε, and G (h) F (h d) for suf?ciently small arguments. For inverse multiquadrics (with constant γ ) these bounds decay exponentially, as we have F (h) = O(e?ρ/ h ), and G (h) = O(h?1 e?γ d/ h ) [8]. For compactly supported radial basis functions ? ∈ PDd ∩ C 2k of minimal degree, the results from [12] state that G (h) = ?d ?2k ?1 ) , and F (h) = O(hk+1/2 ). Thus the condition number of the interpolation O(qX matrix only grows polynomially, in terms of the separation distance, if the latter tends to zero. The previous analysis highlights that, when we interpolate on a set of centers, the best possible situation for balancing the stability and the approximation behaviour is when the uniformity of the set is maximized. This suggests a criterion for selecting suitable subsets of X , on which a good compromise between stability and accuracy can be achieved. This is, in fact, the main feature of the proposed RBF interpolation method, which performs local interpolations on well selected quasi-uniform subsets of X , and obtains a global reconstruction by suitably weighting the local contributions. G (hε) λ(x)

294

G. Casciola et al. / Fast surface reconstruction and hole ?lling

3.

The local approach to the reconstruction

In this section we describe a local method which is designed for the reconstruction of large general data sets. The reconstruction method is based on a technique presented in an earlier paper [7], that is here adapted to the construction of a digital surface model from range data acquired from a physical object. Our local approach to the reconstruction exploits the characteristics of ?exibility and accuracy which have made the radial basis functions a well established tool for multivariate interpolation, overcoming the drawbacks presented by RBF global interpolation methods, such as instability and unacceptable computational cost. The method is a variant of the well-known modi?ed Shepard’s method proposed by Renka [9], and uses RBF interpolants as nodal functions. We consider the set X ? ? Rd , where, for convenience, is assumed to be the convex hull of X , and the corresponding set of function values = {fi }i =1,...,N . For each xk ∈ X, k = 1, . . . , N , we associate a set Xk = {xi ∈ X, i ∈ Ik }, where Ik is the set of indexes of Nq ? 1 suitable neighbours of xk , as well as xk itself, and we determine an RBF nodal function Rk (x) in the form: Rk (x) =

i ∈Ik

ci ? x ? xi

2

,

(10)

by solving on Xk the local interpolation problem: Rk (xi ) = fi , i ∈ Ik .

In order to achieve a real local approach, it is necessary to limit the in?uence of Rk (x) by means of a weight function which decreases with the inverse of the distance from xk . To further control the locality of the reconstruction, Renka in [9] has suggested to associate with each xk a parameter rWk , named radius of in?uence, that truncates the contribution of the k th nodal function outside of the circle with this radius. To choose the radius of in?uence rWk , a certain number NW of points closed to xk is considered, then rWk is taken as the distance from the farthest of these NW points. NW ( Nq ) is thus a parameter which has to be small enough to guarantee the locality of the method. Therefore, we de?ne, for each xk , the k th nodal weight function as Wk (x) = Wk (x)

N i =1

Wi (x)

,

(rWk ? rk (x))+ Wk (x) = rWk rk (x)

p

,

(11)

where rk (x) = x ? xk 2 and rWk is the radius of in?uence around the point xk . Thus the value fk at xk in?uences only the evaluation of points within this radius. It can be shown (see [5]) that the positive weights Wk (x) satisfy the cardinality relations Wk (xi ) = δik , i, k = 1, . . . , N and constitute a partition of unity, namely N k =1 Wk (x) = 1.

G. Casciola et al. / Fast surface reconstruction and hole ?lling

295

Finally, for each evaluation point x, the reconstruction is obtained by a weighted combination of a few local contributions, namely the nodal functions whose radius of in?uence contains x. More precisely, F (x) =

k ∈ Nx

Wk (x)Rk (x),

(12)

where Nx is the set of indexes of all the nodes xk s.t. x ? xk 2 < rWk . Note that the choice of the radial function ?( · ) ∈ C (Rd ) leads to F ∈ C s ( ), where s = min{p, }. Moreover, for each x ∈ , the global approximation error E(x) is given by:

N

E(x) = f (x) ? F (x)

k =1

Wk (x)ek (x),

(13)

ek (x) being the interpolation error (4) of the RBF interpolant Rk (x) on the set Xk . The choice of the sets Xk is therefore crucial for the effectiveness of the reconstruction algorithm and must be made with great care. In fact, as we are interested in the reconstruction of dense range data, a naive choice for Xk , given considering xk and its nearest neighbours, would certainly lead to instability problems. According to the theoretical results of the previous section, we have therefore made a better choice obtained by selecting suitable neighbours of xk , such that the set Xk consists of quasi equi-spread points, thus satisfying the quasi-uniformity conditions (8). This strategy yields an interpolation matrix AXk ,? with minimum condition number still maintaining a good approximation behaviour. Moreover, since each point xj ∈ X is associated with a nodal function Rj (x), the points discarded in the choice of Xk also give a contribution to the ?nal reconstruction, thus avoiding any loss of information. In the case of noisy data, or data which are affected by considerable acquisition errors, a local RBF least squares approximation is a better choice for the nodal function Rk (x). For this aim, we select a subset Sk ? Xk of Nc points, including xk itself, which represents the centers of the radial basis functions considered. Given Sk = {xj ∈ Xk ; j ∈ Jk }, we consider the matrix BXk ,? = (?( xi ? xj 2 ))i ∈Ik ,j ∈Jk , and we solve the following minimization problem min BXk ,? c ? fk 2 ,

c

(14)

where fk is the vector of the function values {fi }, i ∈ Ik , and the entries of c are the coef?cients of the RBF approximant: Rk (x) =

j ∈ Jk

cj ? x ? xj

2

.

(15)

Note that the choice of the centers guarantees that the matrix BXk ,? is full rank, as its columns are in the column space of AXk ,? .

296

G. Casciola et al. / Fast surface reconstruction and hole ?lling

4.

The hole ?lling method

The reconstruction method described in the previous section, for its local nature, preserves the topology of the data sets, that is, roughly speaking, any hole in the data set produces an hole in the reconstruction. This can be considered an advantage both when the original shape presents some empty regions (think, for example, to eye holes in a carnival mask), and when it is necessary to preserve missing regions in an ancient sculpture. But, if we are in the case of lacunary data that we want to restore, the local approach just described can be easily expanded with an ef?cient method for ?lling these holes. This will be the subject of the current section. We now deal with the problem of the reconstruction of an RBF function on a data set (X, ) which presents an empty, lacunary or damaged region. Let us call this region for simplicity the hole. The reconstruction method exploits only the available information on the hole surroundings that we denote by B ? X , thus the processing is limited to the vicinities of the hole. The simplest approach we can consider, is a global reconstruction on B using relation (1); this would lead us to two main drawbacks. First, the gap in the data domain would produce a ?ll distance parameter hB large with respect to qB , thus losing the quasi uniformity conditions. This would give rise to the well known problems regarding accuracy and stability in the reconstruction. Moreover, the global method would ?ll the gap without allowing any kind of control on the shape of the reconstructed region. From the above considerations, we propose a multilayer technique which incrementally expands the information on the hole boundary until the gap has been completely ?lled. This approach reaches the following goals: it maintains the quasi-uniformity conditions for each local interpolant, and enables a local control on the reconstruction by suitably setting the local nodal function parameters. The key idea of the multilayer method is the following. We build a uniform grid, with an estimated grid data density HG , and we overlap it on the boundary data B , thus ?lling completely the missing part in the data domain X . This set of grid points is denoted by G. In order to incrementally extend the original data set (X, ) for including the new set (G, G ), we proceed with a decomposition of B ∪ G into a nested sequence of layers: B = L0 ? L1 ? L2 ? · · · ? LNL ≡ B ∪ G, according to the hole width and to the grid data density. This layer structure, shown in ?gure 1, allows us to ?ll the hole by splitting up the reconstruction into NL steps. Each step m, for m = 1, . . . , NL, consists of the following three phases: ? First we construct the set Lm \Lm?1 by evaluating, using relation (12), an estimate of the reconstruction F (xi ) for each point xi in Lm \Lm?1 , that is: F (xi ) =

k ∈Nxi

Wk (xi )Rk (xi ),

xi ∈ Lm \Lm?1 ,

(16)

G. Casciola et al. / Fast surface reconstruction and hole ?lling

297

Figure 1. Scheme of the layer structure; the white dots represent points in B , while the black are points in G.

where Nxi is the set of indexes of the nodes xk ∈ Lm?1 whose radius of in?uence contains the point xi , and Rk (x) are the corresponding nodal functions. ? In the second phase, the data set is enlarged by including the set (Lm \Lm?1 , Lm \Lm?1 ). ? Then, for each point xk ∈ Lm \Lm?1 , the new nodal function Rk (x) is computed by solving the corresponding local interpolation problem. After the processing of all the layers the original data set has been completed, thus allowing us for the evaluation of the reconstruction everywhere. Determining the layer points represents a critical aspect to the success of the method. As previous mentioned, the local nature of our method allows us to deal always with small sets Xk for each xk ∈ B , characterized by a prede?ned separation distance qXk , and a ?ll distance hXk of the same order of magnitudo. In order to maintain these good properties also for the new points in the layers, the set G has been built with the following grid data density: HG := min hXk .

xk ∈B

Then the subdivision of G into layers is obtained by requiring that each point xk ∈ G is considered into a given layer Lm if there is at least a point in Lm?1 whose distance from xk is less than HG . This choice, in fact, guarantees that the point is in the radius of in?uence of some of its neighbours and that the error given by (13) is limited. The multilayer procedure causes inevitably a propagation of the local interpolation error into the hole region. This is due to the fact that the reconstructed values at each step are not all belonging to the original data set, but they are computed using relation (16) and therefore are affected by a reconstruction error. This error is a convex combination of the errors generated by the local interpolation functions Rk (x). A theoretical study of the propagation of this error is, at ?rst glance, very discouraging. In fact, the upper bound of the difference between the local reconstructions of exact and perturbed data, given in the following proposition, shows that the propagation of the perturbation error can be very large.

298

G. Casciola et al. / Fast surface reconstruction and hole ?lling

Proposition 1. Let (xi , fi ), xi ∈ Xk , fi ∈ , be a set of interpolation points and (xi , f?i ), xi ∈ Xk , fi ∈ , the associated perturbed data set. Using relation (10), the error on the local interpolation function Rk (x), is bounded by Rk (x) ? Rk (x) Proof. Rk (x) ? Rk (x)

i ∈Ik i ∈Ik 1 A? X,? ∞

? x ? xi

i ∈Ik

2

fk ? ? fk

∞

.

(17)

|ci ? c ?i | ? x ? xi max(ci ? c ?i )

i ∈Ik ∞ i ∈Ik

2

? x ? xi ? x ? xi

2

2

? = c?c The result follows replacing relation (3).

.

This theoretical result suggests that, in order to limit the growth factor

1 A? X,? ∞

? x ? xi

i ∈Ik

2

,

the selection of the radial basis functions and of the dimension of the local interpolation problems is really crucial. In our work we are able to maintain a limited growth factor by keeping the system dimension as small as possible, and acting on the inverse multiquadrics constant and on the support of the compactly supported RBF. Moreover this allows us to have a good local control on the shape reconstruction. Our experimental results support our choices revealing a limited error in the reconstruction of test data sets (see section 6) and obtaining a plausible reconstruction by expanding the shape of the boundary data. In the next section, an ef?cient hole ?lling procedure which implements the multilayers method has been integrated in the basic reconstruction algorithm in order to manage any number of holes. 5. The reconstruction algorithm with hole ?lling

The basic algorithm for the reconstruction of a surface F (x), starting from a set (X, ) of N = |X | distinct points, consists of a preprocessing and two main phases. In the preprocessing the data X are scaled in [0, 1] × [0, 1], and they are suitably organized and stored in a grid data structure. Following [9], we used a uniform 2D grid of cells as ef?cient data structure, which allowed us to design a fast procedure for performing a nearest-neighbour search. Algorithm 1 describes the main steps of our approach.

G. Casciola et al. / Fast surface reconstruction and hole ?lling

299

Algorithm 1 (Surface reconstruction). INPUT: (X, ), Nq , NW , q, ? OUTPUT: F (ξ ), ξ ∈ COEFFICIENT PHASE: Compute the coef?cients of the local interpolants Rk (x), k = 1, . . . , N for each xk ∈ X, Construct a set Ik of Nq indexes s.t. Xk = {xk } ∪ {xj ∈ X ; xj ∈ Neigh(xk )} and qXk q Compute the radius of in?uence rWk of xk using NW points; Interpolation step: Build the matrix AXk ,? ∈ RNq ×Nq , AXk ,? = ?( xi ? xj 2 ), i, j ∈ Ik Solve the linear system: AXk ,? c = fk endfor EVALUATION PHASE: Evaluate the reconstructed surface for each evaluation point ξ ∈ : Construct a set Nξ of indexes of all the nodes xk s.t. ξ ? xk 2 < rWk Evaluate Rk (ξ ), ?k ∈ Nξ using (10) Compute the weight Wk (ξ ) , ?k ∈ Nξ , using (11) with p = 2 Compute F (ξ ) = k∈Nξ Wk (ξ )Rk (ξ ) endfor In the ?rst phase (COEFFICIENT PHASE) the most important step is the construction of the sets Xk ? X , as already underlined in the previous section. For this aim, we look for Nq ? 1 points in the area around xk (Neigh(xk )), such that the separation distance of Xk is greater than a given constant q. From an implementation point of view, this means selecting the corresponding set of indexes Ik . In our implementation, q is the same for each Xk . Then we compute the in?uence radii rWk of xk , chosen large enough to include NW points of X . Finally, in order to compute the coef?cients of the nodal functions Rk (x), we solve the linear system AXk ,? c = fk , where the elements of fk are the function values corresponding to the nodes in Xk . The coef?cient matrix AXk ,? is symmetric and positive de?nite. Thus the linear system is solved by applying a Cholesky factorization. The second phase (EVALUATION PHASE) provides the values of the reconstructed function F on a set of points ? R2 . For each evaluation point ξ the algorithm detects the set Nξ of indexes of all the nodes xk , whose radii of in?uence rWk include ξ , and it computes the local contributions (10) in ξ . Thus the reconstructed surface F (ξ ) is evaluated by weighting all the nodal functions according to (12). In order to obtain a reconstructed surface using the least square criterion, a new input parameter Nc is required, which indicates the number of centers to consider for the computation of the local approximant Rk (x). In this case we construct the set Sk ? Xk given by xk and its Nc ? 1 suitable neighbours. This subset Sk can be chosen in several way as, for example, using the farthest star-spread neighbours of xk . We then compute the coef?cients of the nodal function Rk (x) by solving the overdetermined linear system BXk ,? c = fk using

300

G. Casciola et al. / Fast surface reconstruction and hole ?lling

a QR algorithm. Therefore, we follow algorithm 1, in which the interpolation step is replaced with: Approximation step: Construct the set Jk of Nc ( Nq ) indexes of points ∈ Sk Build the matrix BXk ,? ∈ RNq ×Nc , BXk ,? = ?( xi ? xj 2 ), i ∈ Ik , j ∈ Jk Solve the overdetermined linear system: BXk ,? c = fk Note that, in the evaluation phase, the local approximant Rk (x) is given by (15) instead of (10). Now, following section 4, we extend the previous reconstruction algorithm (algorithm 1) by completing it with an ef?cient hole ?lling procedure (named algorithm 2). This has been designed to be able to manage topological complex holes and restore surface shapes, preserving the geometry where it exists, while producing smooth transitions to plausible geometry in unobserved areas. Algorithm 2 (Surface reconstruction with hole ?lling). INPUT: (X, ), Nq , NW , q, ?, Bi , i = 1, . . . , nH OUTPUT: F (ξ ), ξ ∈ PRELIMINARY PHASE for each marked hole i = 1, . . . , nH Overlap a grid of points Gi on Bi ; Label points in Gi in NLi layers, s.t. Bi = L0 ? L1 ? L2 ? · · · ? LNLi ≡ Gi ∪ Bi ; endfor COEFFICIENT PHASE: for the set X HOLE FILLING PHASE: for each marked hole Bi i = 1, . . . , nH for each layer L , = 1, . . . , NLi EVALUATION PHASE: compute { L \L ?1 } X = X ∪ {L \L ?1 }; = ∪ { L \L ?1 }; COEFFICIENT PHASE: compute Rk (x) for xk ∈ {L \L ?1 } endfor endfor EVALUATION PHASE: for the set The critical aspect of algorithms 1 and 2 turns out to be the choice of the free parameters, namely Nq , NW and q , which in?uences both the reconstruction quality and the computational complexity. For the right values of these free parameters we have followed the guidelines given in the theoretical analysis of section 2. Concerning the computational complexity of these algorithms, this is due both to the numerical computation and to the data structure management. The former can be well quanti?ed in the solution of N linear systems of dimension Nq (interpolation case) 3 /6) for the coef?cient using Cholesky factorization, i.e., a computational cost of O(N · Nq

G. Casciola et al. / Fast surface reconstruction and hole ?lling

301

phase, and in the evaluation of the reconstruction in a generic point ξ , whose cost is O(Nq · M), where M = |Nξ |. Regarding the spatial data structure, we can clearly single out a building cost in the preprocessing phase of O(N ), while the complexity of any search inside the ?xed grid depends on the number of cells visited, and therefore on the separation distance q , on the density h and, of course, on the number of searched points Nq for each linear system. 6. Numerical results

The surface reconstruction algorithm described in this work has been implemented and tested with several range images acquired with a PICZA PIX-30 touching probe scanner (Dep. of Math., University of Bologna), and with a VIVID 900 laser scanner (ENEA, Bologna). The scanning technologies considered acquire the shape of an object by generating regular range images. However, since they present outliers and systematic range distorsions, they can be considered quasi-uniform data sets. To illustrate the behaviour of our algorithm, we have reported the results from two data sets. The ?rst data set (angel) has been acquired using PICZA scanner and it does not present any holes. This has been used both to test reconstruction algorithm 1 and to verify the quality of the hole ?lling procedure. For this purpose we have removed parts of the data in different regions to create synthetic holes (angel-h1, angel-h2 and angel-h3) with different topologies. All the regions considered are characterized by convoluted details, because it is well known that, for simple and relatively ?at holes, the reconstruction by radial basis functions naturally works pretty well [2]. The second data set (minerva) presents two holes due to the scanner acquisition, one in the neck and the other under the nose. Moreover, a vertical data zone has intentionally been removed in order to simulate a restoration by hiding the join between the two parts of the statue, while this join has been left on the lips and chin. This data set has been acquired by a VIVID laser scanner and presents inherent noise due to the laser scanner acquisition. The examples have been obtained using the inverse multiquadric radial basis functions, and similar results have been produced by Wendland’s compact support radial basis functions. For each local interpolant/approximant on Xk the constant γk of the inverse multiquadric functions and the support αk of the Wendland’s basis functions, have been automatically adapted to be proportional to the support of the local problems according to the formula: const · maxxi ∈Xk xk ? xi 2 . This yields a common reproducing quality for all the nodal functions Rk (x) and therefore a better quality for the global interpolant/approximant. In our examples, where the data are scaled in [0, 1] × [0, 1], the constant const has been chosen to be 0.35. Figure 2 shows the reconstruction of the angel data set obtained with algorithm 1, using the interpolation criterion. Algorithm 1 has also been applied to the incomplete angel data sets, highlighting the algorithm’s ability to preserve a hole in the data sets when desired. The results of these reconstructions are displayed in ?gures 3(a), 4(a)

302

G. Casciola et al. / Fast surface reconstruction and hole ?lling

Figure 2. Reconstructed angel data set using algorithm 1.

(a)

(b)

Figure 3. Reconstructed angel-h1 data set: (a) algorithm 1; (b) algorithm 2.

(a)

(b)

Figure 4. Reconstructed angel-h2 data set: (a) algorithm 1; (b) algorithm 2.

and 5(a). The corresponding reconstructions, obtained using the hole ?lling procedure (algorithm 2), are given in ?gures 3(b), 4(b) and 5(b). The quality of these reconstructions can be appreciated by comparing them with the reconstructions of the complete data set (?gure 2) and the corresponding reconstruction errors have always been below 5%. Similarly, we applied algorithm 1 to the minerva data set using the interpolation procedure to preserve the inherent noise and the holes of the original data (?gure 6(a)).

G. Casciola et al. / Fast surface reconstruction and hole ?lling

303

(a)

(b)

Figure 5. Reconstructed angel-h3 data set: (a) algorithm 1; (b) algorithm 2.

(a)

(b)

(c)

Figure 6. Reconstructed minerva data set: (a) algorithm 1; (b) algorithm 2 (int.); (c) algorithm 2 (app.). Table 1 Computational results (in seconds) of algorithm 1. DATA SET angel minerva No. points 13956 17050 NW 9 9 Nq 9 9 HX 1.5e?3 1.8e?3 q 0.5e?3 0.5e?3 Tc 4.34 6.19 Te 1.86 1.89 TTOT 6.36 8.28

The surface obtained using algorithm 2 with interpolation is shown in ?gure 6(b), while ?gure 6(c) is an example of removing the noise in the reconstruction using the approximation procedure. Obviously the amount of smoothing applied can easily be modi?ed, considering different choices for the number of centers (Nc ) and the number of points used in the local approximation (Nq ). Note that in ?gure 6(c) we used Nc = 6 and Nq = 20. The illustrations show the high quality of the reconstruction. Another key aspect of the proposed algorithms is their ef?ciency, due to the local nature of the method. Tables 1 and 2 quantify the computational time of the main phases of algorithms 1 and 2, respectively, running on a PC Intel Pentium IV, 2.4 GHz, with 1 Gb of RAM. In the TTOT column we report the execution time of the algorithms, including

304

G. Casciola et al. / Fast surface reconstruction and hole ?lling Table 2 Computational results (in seconds) of algorithm 2.

DATA SET angel-h1 angel-h2 angel-h3 minerva (int) minerva (app)

No. points 13856 13583 13783 17050 17050

NW 9 9 9 9 20

Nq 9 9 9 9 20

HX 1.5e?3 1.5e?3 1.5e?3 1.8e?3 1.8e?3

q 0.5e?3 0.5e?3 0.5e?3 0.5e?3 0.5e?3

Tc 4.30 4.19 4.27 6.18 8.17

Thf 2.78 11.08 5.90 12.16 14.05

Te 2.06 2.56 2.25 5.00 5.02

TTOT 9.34 18.11 12.63 23.66 27.23

the preprocessing phase, while the Tc , Thf and Te report the computational time of the coef?cient, the hole ?lling and the evaluation phases, respectively. The latter strongly depends on the evaluation grid used; all the reconstructed surfaces we considered have been evaluated in a 300 × 300 uniform grid. In the tables we also report the free parameters NW , Nq , the separation distance q and the estimated density HX of the data sets. Note that the latter corresponds to the acquisition step of the scanner and its value represents the ?ll distance of the local interpolation problems and it has been used to estimate the separation distance q . This has allowed us to satisfy the quasi-uniformity conditions (8), guaranteeing a compromise between stability and quality. In fact, the uniformity parameter ρ of the sets Xk turns out to be ρ 0.67 for the angel data set and ρ 0.55 for the minerva data set. Note that, in the minerva data set, we could have considered a slightly bigger value for q yielding a bigger density value, but the algorithm would have been less ef?cient. Acknowledgements This work has been supported by MIUR-Co?n 2002. We are grateful to Ing. S. Petronilli (ENEA, Bologna) for providing the data and motivation for this research. We also thank the referees for many useful comments. References

[1] C.L. Bajaj, F. Bernardini and G. Xu, Automatic reconstruction of surfaces and scalar ?elds from 3D scans, in: Proc. of SIGGRAPH 1995 (ACM Press, New York, 1995). [2] J.C. Carr, R.K. Beatson, J.B. Cherrie, T.J. Mitchell, W.R. Fright, B.C. McCallum and T.R. Evans, Reconstruction and representation of 3D objects with radial basis functions, in: Proc. of SIGGRAPH 2001 (ACM Press, New York, 2001) pp. 67–76. [3] J. Davis, S.R. Marshner, M. Garr and M. Levoy, Filling holes in complex surfaces using volumetric diffusion, in: Proc. of the First Internat. Symposium on 3D Data Processing, Visualization, Transmission (2001). [4] M.S. Floater and A. Iske, Multistep scattered data interpolation using compactly supported radial basis functions, J. Comput. Appl. Math. 73(5) (1996) 65–78. [5] W.J. Gordon and J.A. Wixson, Shepard’s method of metric interpolation to bivariate and multivariate interpolation, Math. Comp. 32(141) (1978) 253–264.

G. Casciola et al. / Fast surface reconstruction and hole ?lling

305

[6] H. Hoppe, T. DeRose, T. Duchamp, J. McDonald and W. Stuetzle, Surface reconstruction from unorganized points, in: Proc. of SIGGRAPH 1992 (ACM Press, New York, 1992). [7] D. Lazzaro and L.B. Montefusco, Radial basis functions for the multivariate interpolation of large scattered data sets, J. Comput. Appl. Math. 140 (2002) 521–536. [8] W.R. Madich and S.A. Nelson, Bounds on multivariate polynomials and exponential error estimates for multiquadric interpolation, J. Approx. Theory 70 (1992) 94–114. [9] R.J. Renka, Multivariate interpolation of large sets of scattared data, ACM Trans. Math. Software 14(2) (1988) 139–148. [10] R. Schaback, Creating surfaces from scattared data using radial basis functions, in: Mathematical Methods for Curves and Surfaces, eds. M. Daehlen, T. Lyche and L. Schumaker (Vanderbilt Univ. Press, Nashville, TN, 1995) 477–496. [11] R. Schaback, Error estimates and condition numbers for radial basis function interpolation, Adv. Comput. Math. 3 (1995) 251–264. [12] H. Wendland, Error estimates for interpolation by compactly supported radial basis functions of minimal degree, J. Approx. Theory 93 (1998) 258–272. [13] H. Wendland, Piecewise polynomial, positive de?nite and compactly supported radial functions of minimal degree, Adv. Comput. Math. 4 (1995) 359–396. [14] Z. Wu, Multivariate compactly supported positive de?nite radial functions, Adv. Comput. Math. 4 (1995) 283–292.

赞助商链接

相关文章:

更多相关标签: