Curvature Regularized Surface Reconstruction from Point Cloud

2020·Arxiv

Abstract

Abstract

We propose a variational functional with a curvature constraint to reconstruct implicit surfaces from point cloud data. In the point cloud data, only locations are assumed to be given, without any normal direction nor any curvature estimation. The minimizing functional balances two terms, the distance function from the point cloud to the surface and the mean curvature of the surface itself. We explore both norm for the curvature constraint. With the added curvature constraint, the computation becomes particularly challenging. We propose two efficient algorithms. The first algorithm is a novel Operator Splitting Method (OSM). It replaces the original high-order PDEs by a decoupled PDE system, which is solved by a semi-implicit method. We also discuss an approach based on an Augmented Lagrangian Method (ALM). The proposed model shows robustness against noise and recovers concave features and corners better compared to models without curvature constraint. Numerical experiments in two and three dimensional data sets, noisy, and sparse data are presented to validate the model. Experiments show that the operator splitting semi-implicit method is flexible and robust.

1 Introduction

In industrial and scientific fields [15, 21], surface reconstruction from point cloud data is a critical step in informative data visualization and successful high-level data processing. Effective methods to reconstruct a continuous surface from finitely many points can reduce the burden in data transmission and facilitate shape manipulations. One of the main goals of surface reconstruction is to render a meaningful and reliable surface which captures the geometrical features of the point cloud. In this paper, we use the celebrated level set function [30] to represent the surface. For -dimensional implicit surface is represented by the set

for a level set function Implicit surfaces enjoy the flexibility in topological changes, and via , one can easily derive and express geometric features of Γ, such as normals, mean curvature, and Gaussian curvature.

In [47, 48], the authors proposed the minimal surface model by interpreting the reconstructed surface as an elastic membrane attached to the given point cloud. It finds the zero-level set surface Γ that minimizes the following energy:

Here, is the area element, s > 0 is an exponent coefficient, the point-to-point-cloud distance, and D is the set of point cloud data. The energy (1) minimizes the surface area weighted by the distance from surface to point cloud.

There are a number of related works using implicit surface reconstruction: a data-driven logarithmic prior for noisy data was considered in [37], surface tension was used to enrich the Euler-Lagrange equations in [18], and principal component analysis was used to reconstruct curves which are embedded in sub-manifolds in [26]. Fast algorithms solving (1) was proposed in [19]. In [24], convexified image segmentation model with a fast algorithm was proposed for implicit surface reconstruction for point clouds. An efficient algorithm incorporating heat kernel convolution and thresholding was introdcued in [43]. In [8], an efficient algorithm for level set method which preserves distance function was proposed. Open surface reconstruction using graph-cuts was proposed in [42], where reconstruction of open surface based on domain decomposition was also proposed. In [25], the authors proposed a variational model consisting of the distance, the normal direction, and the smoothness terms. In [22], a ridge and corner preserving model based on vectorial TV regularization for surface restoration was introduced. In [24], the authors defined the surface via a collection of anisotropic Gaussians centered at each entry of the input point cloud, and used TVG-L1 model for minimization. A similar strategy addressing an gradient regularization model can be found in [23].

We propose a variational functional with a curvature constraint to improve surface reconstruction. One such curvature constraint is the squared mean curvature, , such as Euler’s elastica minimization model [35]. In addition to image inpainting, it has been applied to de-noising [6], segmentation problem [50], and others. For any closed surface Γ in , the bending energy,

where denotes the surface area element, is a conformal invariant [2], and it has a universal lower bound [44]: Another curvature constraint we consider is the absolute mean curvature , which preserves sharp edges and corners in various cases, e.g., denoising [6, 49], and segmentation [1]. As a related work, in [36], graph cuts algorithm was explored for a functional with the absolute mean curvature term. In [7], a variation using a function which is sensitive to large curvature was considered. Other works used weighted mean curvature [16], principle curvature [32], Gaussian curvature [17], Menger curvature [14], and other high-order geometrical information, e.g., conformal factor [27] and elastic ratio [34].

Optimizing a curvature regularized functional is a non-convex and non-linear problem. Computation of such functional is particularly challenging. There are a number of different approaches to design a fast and efficient algorithm, e.g., multigrid method [4], graph-cut algorithm [36], homotopy method [45] and convex relaxation [3, 33], just to name a few. A semi-implicit scheme introduced in [38] simulates the curvature and surface diffusion motion of the interface. One of the major class of methods is based on splitting [10]. The common spirit of these methods is to cast the complicated primal problem into a series of more tamable subproblems, then to find the minimizer using alternative direction method. There are various strategies to obtain such decompositions from the optimization problem. One can derive the associated Euler-Lagrange equations, then apply operator splitting methods on the differential equations, e.g., Lie-Trotter method [41]. A new operator splitting algorithm was proposed for Euler’s elastica model for image smoothing in [6]. One can also introduce auxiliary variables and transform the primal problem into a constrained one, then obtain a series of subproblems by alternatively optimizing one variable at a time while keeping the others fixed; e.g., augmented Lagrangian method (ALM) [1, 39, 40, 46]. We refer the readers to [13, 9, 5, 11] for a detailed discussion on the splitting method in image processing.

In this paper, we propose a variational functional with a curvature constraint to reconstruct implicit surfaces from point cloud data, and explore fast algorithms to solve the associated non-convex, non-linear optimization problem. The minimizing functional balances two terms, the Euclidean distance from the point cloud to the surface and the mean curvature of the surface. We show that the curvature term improves corner reconstruction and recovers non-convex features of the underlying shape of the point cloud data. To avoid dealing with the high-order PDEs resulting from the gradient descent approach, we introduce a semi-implicit method to solve an easier, but equivalent, problem derived by the operator splitting method (OSM). We also explore an ALM method recently proposed by Bae et al. [1], which reduces the number of parameters compared to other curvature regularized models. Our approaches work effectively for 2D/3D cases, as well as for noisy and sparse point cloud. The contributions of this paper are as follows:

1. We propose a new minimal surface model with a curvature constraint, and develop efficient algorithms. It is based on operator splitting and solved by a semi-implicit approach. We also explore an augmented Lagrangian type algorithm and discuss the effects of the parameters involved.

2. We explore and numerically compare results using norms of the mean curvature as the regularization term. We compare the performance of the OSM approach with that of the ALM approach.

This paper is organized as follows. In Section 2, we introduce our curvature regularized minimal surface model and present fast algorithms for the proposed model. In Section 2.1, we introduce a new operator splitting scheme. This approach transforms the original high-order PDE into an easier differential system, which is solved by a semi-implicit method. In Section 2.2, we describe an ALM method for the model, followed by the numerical details in Section 2.3. Various numerical experiments are presented in Section 3, and analytical discussion is presented in Section 4. We conclude this paper in Section 5.

2 Curvature Regularized Surface Reconstruction Model

Let D be the set of given point cloud data, measures the point-to-point-cloud distance, and is the area element. To reconstruct a surface from a given point cloud data D, we propose the following curvature-constrained minimal surface energy:

Here, is the mean curvature of Γ, and the exponent coefficient s > 0 is a constant integer. We explore the cases when s = 1 and s = 2. The first term, the surface integral of the distance from point cloud to the surface, signifies the fidelity of reconstruction. It moves the surface Γ towards the point cloud. The second term, which is the integral of the surface mean curvature along the reconstructed surface, is the regularization. This induces regularized geometric features for Γ independent to the point cloud location. Geometric features can include sharp corners, smooth corners, or straight segments, depending on the choice of s. The parameter 0 controls the influence of the curvature regularization. When minimal surface model proposed in [47].

Remark 2.1. The 1/s power in (2) comes from the original model (1). If one takes the distance function as the potential function of the point cloud, then the energy is an the potential on Γ. Based on this, we add the regularization term related to the mean curvature of Γ, and the 1/s power is used to keep the two terms in the same format. One may remove this power to get a simpler energy and still get the same minimizer when is chosen appropriately.

We employ the implicit surface representation [29] to rewrite the energy (2), and the level set function is defined such that Γ is its zero level set:

Using , the functional (Γ) restricted to the surface Γ can be expressed as

where ) is the Dirac Delta function with H being the Heaviside step function: 0, and 0 otherwise. We use the smooth approximation [1] for practical computation

with 0, which is a constant controlling the smoothness. For any point x on Γ, its mean curvature can be computed as, . Putting these together, we focus on the smoothed energy

and the reconstructed surface is defined as the minimizer of the energy (5), i.e.,

Here represents the optimal level set function.

2.1 Operator Splitting Method

One of our main challenges is that in (3) is highly nonlinear in terms of , and the corresponding Euler-Lagrange equation is a high-order nonlinear PDE. See (30) and (31) in Section 4. To circumvent this difficulty, we propose a new operator splitting strategy, which leads to an equivalent differential equation system that is much easier to solve.

We follow the direction of gradient flow; however, we first decouple the data fidelity term and the curvature regularization term, then minimize the simplified functional via its gradient flow. For s = 2, using (5), we rewrite the energy (2) in the following equivalent form:

with the notation ) as in (4). We then compute the variation of ) with respect to [47]. For denoting the Sobolev space, we have

with

The steady state of (8) is a minimizer of in (6). On the right hand side of the first equation in (8), the two terms are of the same form, only differing by . The first term is the driving velocity to minimize the squared distance from the surface to the given data. The second term is the driving velocity to minimize the squared curvature along the reconstructed surface. The parameter controls the trade-off between these two terms.

We adopt the Lie type of operator splitting and refer the readers to [12] for a complete discussion of different splitting schemes. Given -th step, we update in two fractional steps. In particular, for k > 0, we update the variables through as follows:

and set

and set

We have two subproblems (9) and (10) to address. There is no difficulty to solve (10), since we have the closed form solution

To solve (9) for , the simplest way is to use the explicit scheme as the following:

However, due to the stability consideration, one needs to choose a very small time step of order is the spatial step size. To relax the time step constraint, for some add on both sides of (9), as in [38], to get

We discretize (11) in time semi-implicitly:

φ

We fix and can be solved efficiently by fast Fourier transformation (FFT). The updating formula is summarized as

To solve (12), we need the initial condition (). We choose to be a signed distance function whose zero level set encloses all data. is assigned as ). The algorithm of OSM with s = 2 is stated in Algorithm 1. In Algorithm 1, the reinitialization is used to keep to be a signed-distance function near its zero level set. The details of reinitialization are discussed in Section 2.3.

In the following, we give a brief derivation of OSM for has the same minimizer as

With the Lie type splitting in time and introducing the term on both sides in the first

We note that (6) can be solved similarly by TV method proposed in [46] for image inpainting problem. One advantage of OSM in this paper is its simplicity and having less parameters: once the model parameter, i.e., , is fixed, there is only one parameter, the artificial time step, to tune. Numerical experiments show that there is a wide range of the time step we can choose.

2.2 Augmented Lagrangian Method

We present another efficient algorithm to find the minimizer of (5) with s = 1. The reason for only focusing on s = 1 in this case is that we can take advantage of the shrinkage operator. We first introduce three new variables: ). Finding the minimizer of is equivalent to solving

This can be addressed via alternating direction method of multipliers by introducing Lagrange multipliers . The associated Augmented Lagrangian functional is

where are vectors, are scalars, are fixed constants. To find the saddle point of L, we update each variable in an alternative manner. In each iteration, for each variable, we minimize the corresponding functional while keeping other variables fixed. After all variables are updated, we update Lagrange multipliers. This procedure is repeated until we achieve a steady state. In each iteration, we have four subproblems to minimize:

After those four variables being updated correspondingly, Lagrange multipliers are updated as

These subproblems can be solved efficiently as described in the following. ) in (16), the corresponding Euler-Lagrange equation is:

where 0 is a frozen coefficient. We discretize the time as follows

This is the Laplacian equation of , and we efficiently solve it by FFT.

) can be written as:

where C is independent of q. Then, the minimizer can be found via the shrinkage operator

with

) can be rewritten as

where is independent of ) can be simplified as

Following the idea of Theorem 2 in [1], we can minimize this energy efficiently.

Theorem 2.1. Assume that be the angle between a and the minimum vector of is the angle between . Then the following arguments hold:

•

•

3. if 4. if , the angles satisfy the equation:

and denotes the 2 matrix with the vector being the first and second column respectively.

Note that the condition in Theorem 2.1 is always satisfied since and . From (22), is solved by Newton’s method.

) in (19), the Euler-Lagrange equation is:

where is a small positive number. We discretize in time as

2.3 Numerical Implementation Details

being positive integers, we discretize it by a Cartesian grid with ∆defined on Ω, we use ) to denote for 0 . Denote the standard forward and backward difference as

The gradient, divergence and the Laplacian operators are approximated as follows:

Denote the discrete Fourier transform and its inverse by , respectively. For a function

which gives rise to

and ) can be computed similarly. Both OSM and ALM use FFT to enhance the computational efficiency. The first equation of (14) and the first equation of (20) belong to the same class:

for a function u with constants a, b > 0 and constant c. Using the Fourier transform, we have

Then (25) can be solved by

The equation (23) is in the form of

for some vector valued function are constant positive scalars and is a vector valued constant. After distretization, (26) can be written as

Applying the discrete Fourier transform on both sides of (27), we get

with

Hence, v can be computed by first solving (28) for (and then apply inverse Fourier transform.

For any ) is the distance from x to the collection of the point cloud D, and it can be computed by solving the Eikonal equation

The simplest monotonic scheme to discretize (29) is the Lax-Friedrich scheme which leads to the updating formula [20]:

Figure 1: Effect of in ALM. For fixed induces better reconstruction on the concave part.

This is updated with a fast sweeping method.

To make our algorithm robust, when updating , we reinitialize the level set to be a signed distance function via solving

In practice, after each iteration, we only solve this PDE for a few iterations. For three dimensional space, we use a simple extension of the two dimensional case.

Remark 2.2. In three-dimensional surface reconstruction problems, the point cloud can be large and the computer memory is limited. One can consider a narrow tube which encloses the point cloud, and assume the reconstructed surface lies inside this tube during the evolution. Adopting the local level set method such as [31], only the values of the level set function on grid points inside the tube need to be stored. ALM and OSM can be applied under the local level set method framework except for solving the Laplace equations. Under this framework, one can derive a corresponding linear system for each Laplace equation which can be solved efficiently by the conjugate gradient method.

3 Numerical Results and Comparisons

We present numerical results of the proposed model (2). Without specification, when OSM is used, we refer to Algorithm 1 for model (6) with we refer to Algorithm 2 for model (15) with 1. For a fixed , OSM only has one (the time step) parameter, whereas ALM has three parameters (). We use domain [0two dimensional problems and [0for three dimensional problems. In all examples, used.

In this section, we first consider the effect of parameters for ALM. Second, we compare the performance of ALM and OSM. We find that using OSM with s = 2 gives the best results. We conclude this section by several examples to further explore the performance of OSM when s = 2.

3.1 Choice of Parameters for ALM Method

In the case of ALM, the choice and combinations of the parameters are delicate. When is increased, the reconstruction becomes closer to the point cloud. In Figure 1, we fix vary. With increased , ALM renders the curve closer to the point cloud. In Figure 2, we fix vary; larger better reconstruction. For this example, has little influence, yet, with a large , the results may become unstable or divergent.

Figure 2: Effect of in ALM. For fixed induces better reconstruction on the concave part.

In Figure 3, we fix . With more influence on the curvature, as is increased from 0 to 1, the indent on the rectangle is better reconstructed. When we increase from 1 to 5, we see that, the indent is preserved, yet the tip of the wedge does not extend inward as much as in the case where of encourages both small mean curvature and short curve length. We find that increasing also helps to avoid oscillation during the iteration. In Figure 3, we plot the energy curves corresponding to the these cases. Before the 100-th iteration, these curves are indistinguishable; however, after the 100-th iteration, larger values of suppress the oscillation of the energy curve, which gives a more stable convergence.

Figure 4 shows the robustness against noise. The point cloud is sampled from a circle in , and Gaussian noise with standard deviation 2 is added to the data. Using ALM with 1 and 10. Figure 4 (a) shows similar performances, while the zoomed-in results in Figure 4 (b) show that the larger the the smoother the result becomes.

3.2 Comparison between OSM and ALM

Figure 5 (a) shows the energy convergence comparison between OSM (with s = 2) and ALM (for OSM. The reconstructed curves are shown in Figure 5 (b) and (c). ALM with s = 1 prefers to shorten the length, since it allows sharp corners. There is a balance between the distance term and the regularization term in the functional. OSM performs better in preserving corners, while the results extrude out a little bit at all corners. In this case of s = 2, the reconstruction is more circular, since sharp corners are not allowed.

Using OSM, we fix s = 2 for the distance term in (2) and explore the difference between using (I) no curvature term norm of the mean curvature, and (III) the mean curvature. Figure 6 (a) shows the comparison between (red curve), and s = 2 (blue curve). The blue curve (OSM with s = 2) in (a) is presented separately in Figure 6 (b). OSM using s = 2 gives the best result capturing the structure of the underlying surface more accurately. The rest of the numerical experiments use OSM with s = 2, which gives more stable results with less number of parameters and faster convergence with smaller minimized energies.

3.3 Effect of curvature constraint: OSM with s = 2

Figure 7 shows the comparison between the algorithm from [8] (the first row), row), and OSM with 0 (the third row) for different surfaces. In the algorithm from [8], For the boomerang shape in column (d), the surface constructed by [8] first shrinks to a point

Figure 3: Effect of in ALM. Here reconstruction of the concave wedge. Although in cases, the energy curves are identical before the 100-th iteration, larger suppresses the oscillation of the energy curve: yellow line (is more stable compared to red (

Figure 4: (a) Results by ALM with noisy data and zoom-in of the right-bottom of (a). The noise is additive Gaussian with standard deviation 2. As increases, the curve becomes less oscillatory.

Figure 5: With 10are shown in (b) for ALM and in (c) for OSM: ALM may shorten the curve, while OSM can extrude a corner to make a circular reconstruction.

Figure 6: (a) By OSM with s = 2 for the distance term in (2), the comparison between (I) (b) OSM with s = 2 gives the best result in terms of capturing the structure of the underlying surface more accurately.

Figure 7: Comparison between the algorithm from [8] and OSM with or without curvature constraints: The first row results are by the algorithm proposed in [8] with results are by OSM with curvature constraint (model with the curvature constraint.

and then disappears. From this comparison, OSM with the curvature regularization provides the best results. The performance of the algorithm from [8] is similar to that of OSM without curvature regularization.

With the curvature term, the shape of the underlying surface is better captured in the third row. The interior triangle shape of Figure 7 (a), the first column, is better captured with the curvature term. Without the curvature term, the reconstructed curve does not move further inside, since the prominent part has attained the balance between the curve length and its distance to the two sides of the triangle. Notice that the prominent part in the second row has non-zero curvature. With a positive , this balance is broken and the prominent part will further move towards the upper vertex of the triangle. Also, for the cases with sharp corners in Figure 7 (c)-(d), our model with the curvature term improves the results and recover the underlying shape better.

As increases, different effect can be shown, see Figure 8. As one increases sharp corners are recovered better. However, if is too large, like 3 and 4 in this example, the corners get more circular. As , the weight of curvature term, gets larger, the reconstructed curve becomes even more circular to avoid large curvature.

Figure 9 shows results when the given point cloud is sparse. The point cloud for (a) and (b) is a boomerang shape and that for (c) and (d) is a sparse square with an indent at the bottom. With sparse boomerang data, just using the distance term (part of the given point cloud; see Figure 9 (a). With

Figure 8: Effect of in OSM. (a) increases, the two sharp corners are recovered better. As gets larger, the corners get more circular to avoid large curvature.

Figure 9: By OSM, sparse data results with or without curvature constraints. (a) (b) sparse square shape where only eight corners and one point on each side are given. For both examples, with curvature constraint, the recovery is more accurate and sharper.

of boomerang in Figure 9 (b). For the sparse square shape, only eight corners and one point on each side are given in Figure 9 (c)-(d). While for shape is smooth. With bottom area in Figure 9 (d). Figure 10 shows the case where even less number of points are given. See Figure 10 (a). Only two points around each corner are given. Figure 10 (b) and (c) show results with 5, respectively. Even with extremely sparse data, curvature constraint model can reconstruct the corners well. We observe that when the given point cloud is non-uniform, or data are missing in some region, our algorithm yields results with straight edges between distant points and smooth corners. These are due to the fact that the curvature along a straight edge is zero, and smooth corners have smaller curvature than sharp corners for s = 2 in the discrete setting.

The next experiment is for the noisy boomerang data, where Gaussian noise with standard deviation 1 is added to the locations of the point cloud. The results with 2 are shown in Figure 11. As gets larger, the two lower corners get recovered better. Even with noisy data, OSM shows a strong competence of recovering the sharp corners.

Figure 10: By OSM, extremely sparse data: (a) Given data. (b) The recovered result with 5. Even with extremely sparse data, curvature constraint model can reconstruct the square corners well.

Figure 11: By OSM, reconstruction with noisy data: (a) noise is Gaussian with standard deviation 1. As gets larger, the two lower corners are better recovered.

Figure 12: Examples of three dimensional point cloud data. (a) A pyramid. (b) A yoyo. (c) An ice cream cone.

Figure 13: Reconstruction of the pyramid by OSM with Result with

3.4 Three Dimensional Examples

We conclude this section with experiments of reconstruction of surfaces in three dimensional space. We use OSM with s = 2 to reconstruct the pyramid, the yoyo, and the ice cream cone, whose point clouds are shown in Figure 12. The data in these examples are concentrated within a cube [0. The pyramid has a relatively simple geometry structure: it is convex and its surface only consists of five planes. We can use a large time step ∆t = 500. For the yoyo and the ice-cream cone we use ∆t = 100, since the underlying surfaces have more details, e.g., the neck of the yoyo and the upper concave part of the ice cream cone.

For the pyramid, the reconstructed surfaces with 10 and the comparison of cross sections along y = 25 (a middle section) are shown in Figure 13. In this case, we see limited improvements of capturing the vertices when the curvature constraint is included.

For the yoyo, the reconstructed surface with 5 and the comparison of cross sections along y = 25 are shown in Figure 14. The advantage of the curvature term is obvious. With at some location away from the middle neck part. Since the curvature at that part is non-zero, given a positive , the surface further evolves to capture the neck. For comparison, the reconstructed surface by the algorithm in [8] is shown in Figure 14 (a). Similar to the result by OSM with

The ice cream cone surface consists of two layers and its cross section looks like a boomerang. For this example, if we use reconstructed surfaces with 10 and the comparison of cross sections along y = 25 are

Figure 14: Reconstruction of the yoyo by the algorithm from [8] and OSM with s = 2: (a) Result by the algorithm proposed in [8] with along y = 25.

Figure 15: Reconstruction of the ice cream cone by the algorithm from [8] and OSM with s = 2: (a) Result by the algorithm proposed in [8] with with OSM along y = 25.

shown in Figure 15 (b)-(d). The effect of the value of on this ice cream cone is similar to that on the boomerang. Results with larger values of capture better the features of the underlying surface such as corners. The reconstructed surface by the algorithm in [8] is shown in Figure 15 (a). [8] recovers the bottom corner better but fails to reconstruct the upper concave part of the surface.

4 Analytical Aspects of the Model

We consider the first variation of each term of the functional (2). The first variation of (2) when

which shows the interaction between the data-dependent driving force, d, and the shape geometric feature, . When Γ is close to the point cloud, i.e., d is small, the shape of Γ becomes flexible, i.e., can be large. When Γ is away from the point cloud, i.e., d is large, the shape of Γ becomes rigid, and must be small. For the regularization term, the effect of the meancurvature of Γ is adjusted by the surface area. Notice that this term only focuses on the geometry of Γ, and the point cloud data has no influence. The first variation of the functional can be derived as [7]:

Here is the tangent component of the gradient, and divis its dual operator; W is the Weingarten map of Γ, and |W| equals the Gaussian curvature if Γ is a 2D surface. Compared to (30), the first variation related to the regularization term (31) is more complicated. While the first variation (30) connects the distance and curvature, (31) only depends on the geometric feature of the surface Γ.

When of the minimizer can not be constantly zero, since there is no compact minimal surface. When condition becomes the minimizer where 0, this condition becomes ], while on the region where ]. Hence, the curvature regularization modifies the distance weighted area of the minimizing surface depending on the local concavity/convexity and the Gaussian curvature. These modifications introduce more flexibility when fitting the point cloud, and our experiments show that they can help to improve reconstruction results.

When

where ∆is the Laplace-Beltrami operator, and G is the Gaussian curvature. Since (32) contains ∆, we expect to see that our model with s = 2 will be influenced by the locally averaged mean curvature, which leads to smoothing effects. Here we show this model’s behavior in the following special case.

Proposition 4.1. Suppose the point cloud is sampled from a smooth closed surface Γ with mean curvature and Gaussian curvature G satisfying:

If the point cloud is sufficiently dense, i.e. the computed d is very close to the exact distance function, then Γ is a minimizer of (2) only if it is a sphere of radius

Proof. Since Γ passes through the point cloud, (30) degenerates to 0. Moreover, by (33) together with (31), the necessary condition for Γ being a minimizer is that ∆is a non-zero constant. By (33), this implies that G is also constant; hence, we know that Γ can only be a sphere. Finally, since the radius of the sphere, we can solve for the radius of Γ, which is

The curvature regularization terminflates the membrane supported by the point cloud. Mylar balloon [28], which resembles slightly flattened sphere, satisfies the condition (33). Proposition 4.1 claims that, if the point cloud lies on a Mylar balloon, the minimizer of the functional (2) deviates from the underlying surface.

The following result shows a two dimensional example where the object is a circle, denoted by , which is centered at the origin with radius . It is fair to assume that a local minimizer of (Γ) is a circle, denoted by C, with the same center and radius close to . Denote the radius of C by r. Then we have the following proposition:

Proposition 4.2. Under the setting of the above example, for is a local minimizer of is a local minimizer of (is a local minimizer if

Proof. Note that for a circle with radius r and the same center as ) can be written as

For is a local minimizer for any For whose subdifferential is

For , it can be easily shown that if 12. In other words, if sufficiently close to

For 0 and thus is a local minimizer of then Thus 6 is a local minimizer of

These properties show that the minimizer of the model (2) is not easy to be analyzed even in a simple case such as a circle, and the results heavily depend on the combination of

5 Conclusion

In this paper, we explored the surface reconstruction models based on point cloud data with curvature constraints. The proposed model combines the distance from surface to the point cloud with a global curvature regularization. Introducing this high-order geometric information allows us to impose geometric features around the corners and to reconstruct concave features of the point cloud better. We find that the interactions between two terms of the functional are subtle. For example, for s = 1, since it allows sharp corners, it may be more relaxed around the corners and give shorter length reconstruction. For the curvature term in model (2), larger s gives more weight to the part of the surface which has large curvature, i.e., corners. As a result, the corners of the reconstructed surface are smoothed and extruded out a little bit. For a fast computation, instead of directly solving the complicated terms in its Euler-Lagrange equations, we use a new operator splitting strategy and minimize the energy by a semi-implicit scheme. We also explore an augmented Lagrangian method, which has the advantage of having less parameters compared to other ADMM approaches. Both methods produce reliable results in many cases, including those where the point cloud is noisy or sparse. Comparison between OSM and ALM shows advances of OSM in flexibility, stability, and efficiency. Comparing the results of OSM using s = 1 and s = 2, we find that OSM with s = 2 provides better results when reconstructing features of the point clouds. There are a number of extensions to be considered, including different curvature constraints and using additional information such as surface normal directions to facilitate the reconstruction. Applications to segmentation and image inpainting can also be considered.

Acknowledgment

The authors would like to thank Dr.Martin Huska at University of Bologna, Italy for the valuable discussions exploring different approaches of fast algorithm for high order functionals.

References

[1] E. Bae, X.-C. Tai, and W. Zhu. Augmented Lagrangian method for an Euler’s elastica based segmentation model that promotes convex contours. Inverse Problems & Imaging, 11(1):1–23, 2017.

[2] W. Blaschke. ¨Uber topologische fragen der differentialgeometrie. Jahresbericht der Deutschen Mathematiker-Vereinigung, 38:193–205, 1929.

[3] K. Bredies, T. Pock, and B. Wirth. A convex, lower semicontinuous approximation of Euler’s elastica energy. SIAM Journal on Mathematical Analysis, 47(1):566–613, 2015.

[4] C. Brito-Loeza and K. Chen. Multigrid method for a modified curvature driven diffusion model for image inpainting. Journal of Computational Mathematics, pages 856–875, 2008.

[5] M. Burger, A. Sawatzky, and G. Steidl. First order algorithms in variational image pro- cessing. In Splitting Methods in Communication, Imaging, Science, and Engineering, pages 345–407. Springer, 2016.

[6] L.-J. Deng, R. Glowinski, and X.-C. Tai. A new operator splitting method for the Euler elastica model for image smoothing. SIAM Journal on Imaging Sciences, 2019.

[7] M. Droske and A. Bertozzi. Higher-order feature-preserving geometric regularization. SIAM Journal on Imaging Sciences, 3(1):21–51, 2010.

[8] V. Estellers, D. Zosso, R. Lai, S. Osher, J.-P. Thiran, and X. Bresson. Efficient algorithm for level set method preserving distance function. IEEE Transactions on Image Processing, 21(12):4722–4734, 2012.

[9] R. Glowinski. Admm and non-convex variational problems. In Splitting Methods in Communication, Imaging, Science, and Engineering, pages 251–299. Springer, 2016.

[10] R. Glowinski and P. Le Tallec. Augmented Lagrangian and Operator-Splitting Methods in Nonlinear Mechanics, volume 9. SIAM, 1989.

[11] R. Glowinski, S. Luo, and X.-C. Tai. Fast operator-splitting algorithms for variational imag- ing models: Some recent developments. In Handbook of Numerical Analysis, volume 20, pages 191–232. Elsevier, 2019.

[12] R. Glowinski, S. J. Osher, and W. Yin. Splitting Methods in Communication, Imaging, Science, and Engineering. Springer, 2017.

[13] R. Glowinski, T.-W. Pan, and X.-C. Tai. Some facts about operator-splitting and alter- nating direction methods. In Splitting Methods in Communication, Imaging, Science, and Engineering, pages 19–94. Springer, 2016.

[14] B. Goldluecke and D. Cremers. Introducing total curvature for image processing. In 2011 International Conference on Computer Vision, pages 1267–1274. IEEE, 2011.

[15] L. Gomes, O. R. P. Bellon, and L. Silva. 3D reconstruction methods for digital preservation of cultural heritage: A survey. Pattern Recognition Letters, 50:3–14, 2014.

[16] Y. Gong and O. Goksel. Weighted mean curvature. Signal Processing, 164:329 – 339, 2019.

[17] Y. Gong and I. F. Sbalzarini. Local weighted gaussian curvature for image processing. In 2013 IEEE International Conference on Image Processing, pages 534–538. IEEE, 2013.

[18] J. Haliˇckov´a and K. Mikula. Level set method for surface reconstruction and its application in surveying. Journal of Surveying Engineering, 142(3):04016007, 2016.

[19] Y. He, M. Huska, S. H. Kang, and H. Liu. Fast algorithms for surface reconstruction from point cloud. arXiv preprint arXiv:1907.01142, 2019.

[20] C. Y. Kao, S. Osher, and J. Qian. Lax–Friedrichs sweeping scheme for static Hamilton– Jacobi equations. Journal of Computational Physics, 196(1):367–391, 2004.

[21] D. Khan, M. A. Shirazi, and M. Y. Kim. Single shot laser speckle based 3D acquisition system for medical applications. Optics and Lasers in Engineering, 105:43–53, 2018.

[22] R. Lai, X.-C. Tai, and T. F. Chan. A ridge and corner preserving model for surface restoration. SIAM Journal on Scientific Computing, 35(2):A675–A695, 2013.

[23] H. Li, Y. Li, R. Yu, J. Sun, and J. Kim. Surface reconstruction from unorganized points with gradient minimization. Computer Vision and Image Understanding, 169:108–118, 2018.

[24] J. Liang, F. Park, and H.-K. Zhao. Robust and efficient implicit surface reconstruction for point clouds based on convexified image segmentation. Journal of Scientific Computing, 54(2-3):577–602, 2013.

[25] H. Liu, X. Wang, and W. Qiang. Implicit surface reconstruction from 3D scattered points based on variational level set method. In 2008 2nd International Symposium on Systems and Control in Aerospace and Astronautics, pages 1–5. IEEE, 2008.

[26] H. Liu, Z. Yao, S. Leung, and T. F. Chan. A level set based variational principal flow method for nonparametric dimension reduction on Riemannian manifolds. SIAM Journal on Scientific Computing, 39(4):A1616–A1646, 2017.

[27] L. M. Lui, C. Wen, and X. Gu. A conformal approach for surface inpainting. Inverse Problems and Imaging, 7(3):863–884, 2013.

[28] I. M. Mladenov and J. Oprea. The mylar ballon: New viewpoints and generalizations. In Proceedings of the Eighth International Conference on Geometry, Integrability and Quantization, pages 246–263. Institute of Biophysics and Biomedical Engineering, Bulgarian Academy of Sciences, 2007.

[29] S. Osher, R. Fedkiw, and K. Piechor. Level set methods and dynamic implicit surfaces. Appl. Mech. Rev., 57(3):B15–B15, 2004.

[30] S. Osher and J. A. Sethian. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations. Journal of Computational Physics, 79(1):12–49, 1988.

[31] D. Peng, B. Merriman, S. Osher, H. Zhao, and M. Kang. A PDE-based fast local level set method. Journal of computational physics, 155(2):410–438, 1999.

[32] X. Qiao. The principle curvature-driven diffusion model for image de-noising. In 3rd International Conference on Multimedia Technology (ICMT-13). Atlantis Press, 2013.

[33] T. Schoenemann, F. Kahl, S. Masnou, and D. Cremers. A linear framework for region- based image segmentation and inpainting involving curvature penalization. International Journal of Computer Vision, 99(1):53–68, 2012.

[34] T. Schoenemann, S. Masnou, and D. Cremers. The elastic ratio: Introducing curvature into ratio-based image segmentation. IEEE Transactions on Image Processing, 20(9):2565– 2581, 2011.

[35] J. Shen, S. H. Kang, and T. F. Chan. Euler’s elastica and curvature-based inpainting. SIAM Journal on Applied Mathematics, 63(2):564–592, 2003.

[36] J. Shi, M. Wan, X.-C. Tai, and D. Wang. Curvature minimization for surface reconstruction with features. In International Conference on Scale Space and Variational Methods in Computer Vision, pages 495–507. Springer, 2011.

[37] Y. Shi and W. C. Karl. Shape reconstruction from unorganized points with a data-driven level set method. In 2004 IEEE International Conference on Acoustics, Speech, and Signal Processing, volume 3, pages iii–13. IEEE, 2004.

[38] P. Smereka. Semi-implicit level set methods for curvature and surface diffusion motion. Journal of Scientific Computing, 19(1):439–456, 2003.

[39] X.-C. Tai, J. Hahn, and G. J. Chung. A fast algorithm for Euler’s elastica model using augmented Lagrangian method. SIAM Journal on Imaging Sciences, 4(1):313–344, 2011.

[40] X.-C. Tai and C. Wu. Augmented Lagrangian method, dual methods and split Bregman iteration for ROF model. In International Conference on Scale Space and Variational Methods in Computer Vision, pages 502–513. Springer, 2009.

[41] H. F. Trotter. On the product of semi-groups of operators. Proceedings of the American Mathematical Society, 10(4):545–551, 1959.

[42] M. Wan, Y. Wang, E. Bae, X.-C. Tai, and D. Wang. Reconstructing open surfaces via graph-cuts. IEEE transactions on visualization and computer graphics, 19(2):306–318, 2012.

[43] D. Wang. An efficient iterative method for reconstructing surface from point clouds. arXiv preprint arXiv:2005.11864, 2020.

[44] T. Willmore. Mean curvature of immersed surfaces. An. Sti. Univ.Al. I. Cuza Iasi, Sec. I. a Mat.(NS), 14:99–103, 1968.

[45] F. Yang, K. Chen, and B. Yu. Homotopy method for a mean curvature-based denoising model. Applied Numerical Mathematics, 62(3):185–200, 2012.

[46] M. Yashtini and S. H. Kang. A fast relaxed normal two split method and an effective weighted TV approach for Euler’s elastica image inpainting. SIAM Journal on Imaging Sciences, 9(4):1552–1581, 2016.

[47] H.-K. Zhao, S. Osher, and R. Fedkiw. Fast surface reconstruction using the level set method. In Proceedings IEEE Workshop on Variational and Level Set Methods in Computer Vision, pages 194–201. IEEE, 2001.

[48] H.-K. Zhao, S. Osher, B. Merriman, and M. Kang. Implicit and nonparametric shape reconstruction from unorganized data using a variational level set method. Computer Vision and Image Understanding, 80(3):295–314, 2000.

[49] W. Zhu and T. Chan. Image denoising using mean curvature of image surface. SIAM Journal on Imaging Sciences, 5(1):1–32, 2012.

[50] W. Zhu, X.-C. Tai, and T. Chan. Image segmentation using Euler’s elastica as the regu- larization. Journal of Scientific Computing, 57(2):414–438, 2013.