Three-Dimensional Surface Reconstruction from Point Clouds Using Euler’s Elastica Regularization

Abstract

Euler’s elastica energy regularizer, initially employed in mathematical and physical systems, has recently garnered much attention in image processing and computer vision tasks. Due to the non-convexity, non-smoothness, and high order of its derivative, however, the term has yet to be effectively applied in 3D reconstruction. To this day, the industry is still searching for 3D reconstruction systems that are robust, accurate, efficient, and easy to use. While implicit surface reconstruction methods generally demonstrate superior robustness and flexibility, the traditional methods rely on initialization and can easily become trapped in local minima. Some low-order variational models are able to overcome these issues, but they still struggle with the reconstruction of object details. Euler’s elastica term, on the other hand, has been found to share the advantages of both the TV regularization term and the curvature regularization term. In this paper, we aim to address the problems of missing details and complex computation in implicit 3D reconstruction by efficiently using Euler’s elastica term. The main contributions of this article can be outlined in three aspects. Firstly, Euler’s elastica is introduced as a regularization term in 3D point cloud reconstruction. Secondly, a new dual algorithm is devised for the proposed model, significantly improving solution efficiency compared to the commonly used TV model. Lastly, numerical experiments conducted in 2D and 3D demonstrate the remarkable performance of Euler’s elastica in enhancing features of curved surfaces during point cloud reconstruction. The reconstructed point cloud surface adheres more closely to the initial point cloud surface when compared to the classical TV model. However, it is worth noting that Euler’s elastica exhibits a lesser capability in handling local extrema compared to the TV model.

1. Introduction

Three-dimensional reconstruction using point clouds has great academic and commercial significance, but the challenge lies in the fact that the point clouds are often disorganized and disconnected, making it difficult to accurately reconstruct surfaces. There are two main methods for the reconstruction of the surface, the data-based method and the model-based method. Data-based methods are typically deep learning methods, such as the Deep level set [1], DeepSDF [2], EdgeConv [3], Point-BERT [4], IterativePFN [5], etc. These method can achieve good results but require a large amount of training data. On the other hand, model-based methods can produce good results without training data [6,7]. Moreover, they can be used as a priori methods [8] to maintain the smoothness of the reconstructed surface, fill in missing surface points, and adjust misalignment [9]. For the above reasons, the classical model-based methods are still being actively researched and improved upon.

There are two general approaches to describe a surface: explicit and implicit. Explicit representation [10,11,12] accurately represents the location of points on the reconstructed surface and is efficient and accurate, but it is less robust and not flexible enough to handle arbitrary and dynamic surface topologies. On the other hand, implicit representation [13,14,15,16] starts with an initial surface and iteratively evolves it through a level set function [17] to achieve an accurate reconstruction. Unlike explicit methods, implicit methods have greater robustness and flexibility, making them better suited for reconstructing noisy and non-uniform point clouds [14,15,18].

The implicit surface reconstruction method proposed by Zhao et al. in [13,19] is a successful algorithm that uses the variational level set method and either the convection method or the energy-minimizing gradient flow method to reconstruct surfaces. This method is efficient in creating an initial surface, but it still has three main challenges. Firstly, the surface evolution process requires reinitialization, which ensures that the evolving level set is equivalent to a signed distance function, but this impacts the calculation efficiency and may lead to the final reconstructed surface shrinking [20]. Secondly, the model is non-convex, and the initial surface greatly affects the reconstruction results. Lastly, if the point cloud data encompasses delicate details or curved inward regions, the initial surface reconstruction may become trapped at a local extreme, despite the accuracy of the initial reconstruction.

To better reconstruct the surface, a method is proposed for reconstructing point cloud surfaces using the L1-sparse method. This method effectively ignores the influence of Gaussian noise, and its global optimization property ensures that the reconstructed surface is not limited to the local extrema. Based on this work, Duan [21,22] proposed a novel variational method using the TV method with a balloon force and a fitting term (TV${}_g$). The method can overcome the issue of getting stuck in local minima but still cause the surface to lose many details. However, when using the TV model for point cloud reconstruction, it can result in the stair-step effect and contrast loss, leading to clearly visible jagged regions in the reconstructed image and the loss of important texture details. To address these issues related to the stair-step effect and contrast loss, He [23] proposed to adopt a curvature regularized term instead of the TV regularized term. Experiments show that the curvature regularized term can be used to better reconstruct the smooth sections of the surface but performs worse around the corners. It remains to be seen whether there is a method that can combine the advantages of the low-order model and the high-order model.

On the other hand, Euler’s elastica term has long been used in image processing, especially in image segmentation [24,25] and denoising [26]. Using Euler’s elastica as a regularization term to address issues such as the staircase effect and contrast loss in low-order models is currently a popular topic [27,28]; however, due to the challenges associated with its application in point cloud reconstruction, there is currently a lack of research in this specific area. Through decades of research, this regularizer has been found to maintain some of the advantages of the TV regularization term and the curvature regularization term at the same time. Utilizing Euler’s elastica regularization as the energy function effectively mitigates the stair-step effect and contrast loss problems inherent in the TV model. It can also address the problem of excessively smooth inflection points in curvature models, resulting in point cloud reconstructions that preserve more details. However, because of its complexity, this term comes with greater computational cost, so many new research works focused on designing a solution for Euler’s elastica term, such as Chan et al.’s ALM method [29] and Duan et al.’s fast method [30]. With the help of these methods, the efficiency of Euler’s elastica term applied to image processing problems is significantly improved. Given the aforementioned advantages, the decision is made to adopt Euler’s elastica term as the regularization term.

To resolve the issues previously identified with the traditional implicit reconstruction technique, we present a new method called Euler’s elastica for creating precise surfaces from point cloud data. This approach transforms the surface reconstruction process into a task of implicit image segmentation. Our approach involves creating a boundary indicator function, a enclosed image, and an initial surface from the input point cloud data. The edge indicator function plays a crucial role in identifying the precise surface area that needs to be reconstructed. The bounded image acts as a marker, indicating if a specific grid point lies inside or outside the object. The level set method used in the proposed model is successfully initialized by the initial surface.

Proposed surface reconstruction approach involves two stages. Stage 1: We use the fast sweeping method to compute a distance function from the point cloud data. This distance function is then utilized: (a) as an edge indicator function, (b) to generate an image (2D or 3D) surrounded by the point clouds, and (c) to determine an initial surface. Stage 2: The initial surface generated in Stage 1 is then advanced through the application of the proposed variational level set model, which effectively combines the edge indicator and the image from Stage 1. Optimizing this model involves utilizing the fast dual method, as detailed in [31], to simplify the computation of intricate curvature terms. The proposed method minimizes the variational model through gradient descent optimization to ultimately produce a precise surface from the point cloud data. Experiments indicate that the Euler elastica method that we propose surpasses other implicit methods in terms of accuracy.

In this work, three main research contributions are presented. Firstly, we introduce Euler’s elastica regularization for the first time in 3D surface reconstruction. By utilizing the Euler elastica function as the regularization term, it addresses the issues of the staircasing effect and contrast loss observed in the TV model. This is mainly because the Euler elastica energy function imposes constraints on the deformation field by considering the elastic properties of objects [32,33], leading to more plausible and realistic deformations. It helps avoid unreasonable local deformations and ensures that the deformation field globally retains the shape and structural features of objects, thus mitigating the occurrence of the staircasing effect. Moreover, the TV model may encounter gradient vanishing problems in certain cases, resulting in the loss of reconstruction details. The Euler elastica energy function partially alleviates this problem, as it is associated with the high-order gradients of the image or deformation field, thereby preserving more detail information. Secondly, we employ a new dual projection algorithm to solve the proposed Euler elastica model. Recently, the alternating direction method of multipliers (ADMMs) algorithm has widely been used for efficiently solving various energy functions incorporating Euler’s elastica regularization term [30,34,35]. However, this method introduces a significant number of additional parameters and operators, which can be a drawback. So the fast dual projection algorithm can be used to instead of the ADMM algorithm. The proposed method introduces only one parameter, significantly speeding up the model computation. We will provide a detailed explanation in the section on the proposed variational method. Thirdly, we validate the effectiveness of proposed method on both 2D and 3D levels. Through numerical experiments, we provide a detailed exploration of the advantages and limitations of the proposed approach. Additionally, we apply the proposed method to face reconstruction, further discussing its potential applications. The advantages include reconstructing surfaces with continuous and nonlinear curvature changes, preserving more details compared to the TV model, and mitigating the stair-step effect and contrast loss. However, a notable drawback is that the Euler elastic model may have limited ability in handling local extrema compared to the TV model.

This paper is structured as follows: Section 2 provides background knowledge. Section 3 offers a detailed description of the Euler elastic energy model proposed in this paper. Section 4 outlines the methods used to evaluate our reconstruction results. Section 5 presents the experimental setup, encompassing both 2D and 3D aspects to illustrate the generality and effectiveness of the proposed approach. Finally, in Section 6, we summarize the findings of this paper and offer prospects for future work.

2. Variational Image 3D Reconstruction Theory

As previously mentioned, the implicit surface reconstruction method proposed by Zhao et al. in [13,19] was one of the first effective methods in this field. Given an unorganized 3D point cloud $\left\{x_i\right\}$, the reconstructed surface is achieved through minimizing the energy of the weighted geodesic active contour (GAC):

$$E(\Gamma )={\int }_{\Gamma }d\left(x\right)d\Gamma ,$$

(1)

The passage discusses the concept of the Snake model, a type of mathematical method known as the variational method. This method can be applied to either a curve in two-dimensional space or a surface in three-dimensional space. The method uses the arc length element or area parametrization as a way of measuring the curve or surface, represented by $\Gamma$. The Snake model was first introduced in a paper by Kass et al. and was later expanded upon by the work of Caselles, Kimmel, and Sapiro in the framework of level sets. The distance function, represented by $d\left(x\right)$ in Equation (1), is an important component of the Snake model and is subject to the Eikonal equation:

$$\left|\nabla d\left(x\right)\right|=w\left(x\right),\forall x\in \Omega \setminus \left\{x_i\right\},w\left(x\right)=1,d\left(x_i\right)=0,\forall x\in \left\{x_i\right\}.$$

(2)

The implementation of Equation (2) using a fast sweeping algorithm was described in detail in [36,37] for both 2D and 3D scenarios. The results of applying the fast sweeping method to two sets of point clouds, (a) and (b), are presented in Figure 1. The output of this method is distance maps, which are displayed in (c) and (d) of the figure. To construct an initial image bounded by the point clouds, the distance $d\left(x\right)$ is used as an edge indicator function.

The authors of [14] utilized the distance function $d\left(x\right)$ calculated from point clouds to determine the boundaries of the point cloud. Given an image $f\left(x\right)$, they suggested using the following as the edge indicator function for $f\left(x\right)$:

$$g\left(x\right)=\frac{1}{\epsilon +|\nabla f\left(x\right){|}^p},$$

where the variables p and $\epsilon$ are both positive values, and g is nearly equal to zero at the edge positions. In order to use the distance function $d\left(x\right)$ as the edge indicator, the authors added a small positive value $\epsilon$ to regularize it, resulting in $g\left(x\right)=d\left(x\right)+\epsilon$. To calculate the initial image $f\left(x\right)$, they then solved the Eikonal equation, which is a mathematical equation used to describe the speed of propagation of waves or fronts:

$$|\nabla f\left(x\right)|=\frac{1}{\epsilon +g{\left(x\right)}^{1/p}}.$$

The authors of [15] utilized the concept of the inner product field to distinguish between the inner and outer regions of a point cloud data set. The basis for computing the image $f\left(x\right)$ is represented by the following equation:

$$f\left(x\right)=\left(x-cp\left(x\right)\right)·n\left(cp\left(x\right)\right).$$

In the equation, x is a point in the defined area $\Omega$. $cp\left(x\right)$ stands for the closest point on the point cloud to x, and $n\left(cp\right(x\left)\right)$ is the normal vector pointing away from $cp\left(x\right)$. The output of $f\left(x\right)$ is negative if x is within the point cloud and positive if x lies outside the point cloud.

In summary, the proposed method starts with a binary image $I\left(x\right)$, which is obtained by thresholding the distance function $d\left(x\right)$. The fast sweeping algorithm is then applied on the binary image to locate it. The right-hand side of Equation (2) is modified by substituting $W\left(x\right)=1$ with $W\left(x\right)=I\left(x\right)$, making $d\left(x\right)$ no longer a distance function. The grid points at the edge of the domain $\Omega$ are assigned a value of zero, while the other internal grid points are given extremely high values. Finally, the rapid sweeping procedure with eight directional sweeping Gauss–Seidel iterations is used to solve $d\left(x\right)$. This method offers a more straightforward way to calculate the initial image bounded by the point cloud data. The proposed method is a simpler way to compute the initial image enclosed by the point clouds compared to other methods. By thresholding the distance function $d\left(x\right)$ obtained from (2), an annular binary image $I\left(x\right)$ is created. This image is then located using the fast sweeping algorithm. On the right-hand side of (2), $W\left(x\right)$ is substituted with $I\left(x\right)$, which makes $d\left(x\right)$ no longer a distance function. The grid points on the six border faces of the domain $\Omega$ are set to zero, and other inner grid points are set to very large values, to solve $d\left(x\right)$ using the rapid sweeping procedure with eight directional sweeping Gauss-Seidel iterations. Once $d\left(x\right)$ is calculated, the new function $f\left(x\right)$ is determined by setting $f\left(x\right)=d\left(x\right)$. The result of thresholding $f\left(x\right)$ to create a new binary image is demonstrated in the final two images of Figure 2. The border points of the object in the binary picture can be found using a simple algorithm such as Marching Cube, and these new points $\left\{{\tilde{x}}_i\right\}$ and the original disorganized point clouds $\left\{x_i\right\}$ are close to each other. To obtain the signed distance function ${\varphi }_0$, the fast sweeping algorithm is performed again on the new point clouds $\left\{{\tilde{x}}_i\right\}$. The inside of the object is set to be negative, and the outside is set to be positive. The signed distance function ${\varphi }_0$ shown in Figure 3 is used as the initial level set function in the proposed method.

The functional of the new method called TV${}_g$ combined balloon force (TV${}_g$-B) was introduced by Duan in [21] and is composed of three separate energy terms, as follows:

$$E(\varphi ,c_1,c_2)=E_R\left(\varphi \right)+\lambda E_I(\varphi ,c_1,c_2)+\beta E_B\left(\varphi \right).$$

(3)

The first term, $E_R$, is responsible for maintaining the level set function $\varphi$ as a signed distance function, which serves as a regularization term to ensure that the reconstructed surface is smooth:

$$E_R\left(\varphi \right)={\int }_{\Omega }d\left(x\right)|\nabla H\left(\varphi \right)|dx+\frac{\mu }{2}{\int }_{\Omega }{\left(\right|\nabla \varphi |-1)}^{2}dx,$$

(4)

The second term $E_R$ in the functional of the TV${}_g$-B method helps keep the level set function $\varphi$ as a signed distance function. The weight $\mu$ ensures that $\varphi$ remains a signed distance function, with larger values of $\mu$ making $\varphi$ closer to a signed distance function. This term corresponds to the GAC model and uses a weighted TV-norm. The third term ensures that the level set function $\varphi$ stays as a signed distance function, eliminating the need for repeated initialization.

The information acquired from the data set is incorporated into the second data fitting term $E_I$, which is:

$$E_I(c_1,c_2,\varphi )={\int }_{\Omega }Q(c_1,c_2)H\left(\varphi \right)dx,$$

(5)

where the values $c_1$ and $c_2$ denote the mean inside and outside the zero level set of $\varphi$, respectively. The function $Q(c_1,c_2)={(c_1-f\left(x\right))}^2-{(c_2-f\left(x\right))}^2$ is established to enhance the precision of the reconstruction by accounting for equal intensity levels of the fine details in the object. The penalty factor $\lambda$ mentioned in Equation (3) has a positive value.

The area information included in the zero level set of $\varphi$ is incorporated into the third balloon force term $E_B$ in order to speed up surface evolution and segment concave objects. The term is expressed as:

$$E_B\left(\varphi \right)={\int }_{\Omega }d\left(x\right)H(-\varphi )dx.$$

(6)

The value of the parameter $\beta$ in $E_B$ depends on whether the zero level set is reached or not. If the initial zero level set function is outside of the object, then the inside of the level set is less than 0 and the outside is greater than 0. In this case, $\beta$ takes on a positive value, which ensures that the contour line shrinks during evolution. Conversely, if the initial zero level set function is inside of the object, then $\beta$ takes on a negative value, which causes the boundary to expand. This mechanism helps to segment concave objects by ensuring that the level set function always contains the object.

An optimization process can improve the performance of the TV${}_g$-B technique (3). Initiating the process involves maintaining the constancy of the value of $\varphi$. This allows us to calculate the values of $c_1$ and $c_2$ as outlined below:

$$c_1=\frac{{\int }_{\Omega }fH\left(\varphi \right)dx}{{\int }_{\Omega }H\left(\varphi \right)dx}\mathrm{and}{c}_2=\frac{{\int }_{\Omega }f(1-H\left(\varphi \right))dx}{{\int }_{\Omega }(1-H\left(\varphi \right))dx}.$$

(7)

After maintaining the constant value of $\varphi$, the optimization process continues by starting with $\varphi ={\varphi }_0$, where ${\varphi }_0$ can be calculated by the fast sweeping method. The minimization of Equation (3) is then carried out by utilizing the gradient descent method as illustrated below:

$$\begin{array}{c}\hfill \frac{\partial \varphi }{\partial t}=\left(\nabla ·\left(d\frac{\nabla \varphi }{\left|\nabla \varphi \right|}\right)-\lambda Q(c_1,c_2)+\beta d\right)\delta \left(\varphi \right)+\mu \left(\Delta \varphi -\nabla ·\frac{\nabla \varphi }{\left|\nabla \varphi \right|}\right).\end{array}$$

(8)

In practical applications, the regularized versions of the Heaviside function $H\left(\varphi \right)$ and Dirac delta function $\delta \left(\varphi \right)$ in Equations (3) and (8) are approximated by using a small positive value $\epsilon$. This approximation helps to avoid numerical issues and ensure stability during the optimization process:

$$H_{\epsilon }\left(\varphi \right)=\frac{1}{2}+\frac{1}{\pi }arctan\left(\frac{\varphi }{\epsilon }\right),{\delta }_{\epsilon }\left(\varphi \right)=\frac{1}{\pi }\frac{\epsilon }{{\epsilon }^2+{\varphi }^2}.$$

3. Proposed Variational Method

When reinitializing the symbol distance function, we consider using the projection method instead of the first term in (3), then the (3) is rewritten as follows:

$$\begin{array}{cc}\hfill E_R\left(\varphi \right)& ={\int }_{\Omega }d\left(x\right)\left(a+b{\left(\nabla ·\frac{\nabla H\left(\varphi \right)}{|\nabla H\left(\varphi \right)|}\right)}^2\right)|\nabla H\left(\varphi \right)|dx,s.t.|\nabla \varphi |=1.\hfill \end{array}$$

(9)

The complex form of Euler’s elastica affects the computational efficiency. If we use the ADMM algorithm proposed by Tai [34] to solve Equation (9), we can obtain the following expression:

$$\begin{array}{cc}\hfill \underset{u,v,\overrightarrow{w},\overrightarrow{n},\overrightarrow{m}}{min}E\left(u,v,\overrightarrow{w},\overrightarrow{n},\overrightarrow{m},\right)=& {\int }_{\Omega }Q\left(c_1,c_2\right)vdx+\gamma {\int }_{\Omega }[a+b|\nabla ·\overrightarrow{n}{|}^2\left]\right|\overrightarrow{w}|dx\hfill \\ \hfill & +\beta {\int }_{\Omega }-d\left(x\right)vdx+\frac{\mu }{2}{\int }_{\Omega }\left(\right|\overrightarrow{w}{|-1)}^{2}dx\hfill \\ \hfill & +{\int }_{\Omega }{\lambda }_1\left(|\overrightarrow{w}|-\overrightarrow{w}·\overrightarrow{m}\right)dx+{\gamma }_1{\int }_{\Omega }\left(|\overrightarrow{w}|-\overrightarrow{w}·\overrightarrow{m}\right)dx\hfill \\ \hfill & +{\int }_{\Omega }\overrightarrow{{\lambda }_2}·\left(\overrightarrow{w}-\nabla u\right)dx+\frac{{\gamma }_2}{2}{\int }_{\Omega }{|\overrightarrow{w}-\nabla u|}^{2}dx\hfill \\ \hfill & +{\int }_{\Omega }{\lambda }_3\left(v-u\right)dx+{\gamma }_3{\int }_{\Omega }{\left(v-u\right)}^{2}dx\hfill \\ \hfill & +{\int }_{\Omega }\overrightarrow{{\lambda }_4}·\left(\overrightarrow{n}-\overrightarrow{m}\right)dx+\frac{{\gamma }_4}{2}{\int }_{\Omega }{|\overrightarrow{n}-\overrightarrow{m}|}^{2}dx\hfill \end{array}$$

(10)

Obviously, the ADMM algorithm introduces four splitting operators $v,\overrightarrow{w},\overrightarrow{n},\overrightarrow{m}$ and four Lagrange multipliers ${\lambda }_1,\overrightarrow{{\lambda }_2},{\lambda }_3,\overrightarrow{{\lambda }_4}$, which undoubtedly impose a significant computational burden on the model. To simplify it, we introduce auxiliary dual variables $\overrightarrow{p}$. This process significantly reduces the computational overhead at each iteration. Furthermore, ADMM requires the tuning of four extra parameters, e.g., ${\gamma }_1$, ${\gamma }_2$, ${\gamma }_3$, and ${\gamma }_4$, whereas our method has only one extra parameter (e.g., $\gamma$) to tune. This introduces a considerable level of complexity and difficulty in utilizing ADMM effectively. With these variables, it is easy to rewrite the energy function in (3) as an Augmented Lagrangian Function, which forms the basis of our proposed dual algorithm for the problem:

$$\begin{array}{cc}\hfill \underset{u}{min}\underset{|\overrightarrow{p}|\leqq g}{max}E(\varphi ,\overrightarrow{p})& =\gamma {\int }_{\Omega }d\left(x\right)\left(\nabla ·\overrightarrow{p}\right)H\left(\varphi \right)dx+\lambda {\int }_{\Omega }Q(c_1,c_2)H\left(\varphi \right)dx\hfill \\ \hfill & +\beta {\int }_{\Omega }d\left(x\right)H(-\varphi )dx+\frac{\mu }{2}{\int }_{\Omega }{\left(\right|\nabla \varphi |-1)}^{2}dx\hfill \\ \hfill & s.t.g=\left(a+b{\left(\nabla ·\frac{\nabla \varphi }{|\nabla \varphi |}\right)}^2\right)\hfill \end{array}$$

(11)

The solution of $c_1$, $c_2$ is the same as (7), but the formula of $\varphi$ is simplified to the following gradient reduction equation:

$$\begin{array}{cc}\hfill \frac{\partial \varphi }{\partial t}=& d\left(x\right)\left(\nabla ·\overrightarrow{p}\right)\delta \left(\varphi \right)+\lambda Q(c_1,c_2)\delta \left(\varphi \right)-\beta d\left(x\right)\delta \left(\varphi \right)+\mu \left(\Delta \varphi -\nabla ·\frac{\nabla \varphi }{\left|\nabla \varphi \right|}\right),\hfill \end{array}$$

(12)

where t is the time step. After $\varphi$ is solved, the other problem $\overrightarrow{p}$ can be optimized as follows:

$$\left\{\begin{array}{cc}\hfill & \frac{\partial \overrightarrow{p}\left(x,t\right)}{\partial t}=\gamma \nabla \left(d\left(x\right)\delta \left(\varphi \right)\right)\hfill \\ \hfill & |\overrightarrow{p}|\leqq g\hfill \end{array}\right.,$$

(13)

where the value of g is directly determined by the $\varphi$ obtained in the previous step. Therefore, the numerical implementations of the function (13) can be calculated as follows:

$$\left\{\begin{array}{cc}\hfill & {\overrightarrow{p}}_{i,j}^{k+1}={\overrightarrow{p}}_{i,j}^k+t\gamma \nabla \left(d\delta \left({\varphi }_{i,j}^{k+1}\right)\right)\hfill \\ \hfill & {\overrightarrow{p}}_{i,j}^{k+1}=\frac{g_{i,j}^{k+1}\overrightarrow{p_{i,j}^k}}{\max \{g_{i,j}^{k+1},|\overrightarrow{p_{i,j}^k}\left|\right\}}\hfill \end{array}\right.,$$

(14)

In each iteration, the following error tolerances should be checked to determine convergence, i.e.,

$${\Sigma }^{k+1}=\frac{\parallel E^{k+1}-E^k\parallel }{\parallel E^k\parallel }\le \mathrm{Tol}.$$

(15)

where $\mathrm{Tol}=0.01$.

The advantage of this method is very obvious; not only is the calculation process greatly simplified but also the convergence speed of the model will increase. We give detailed experimental explanations in the result section. Discretizing the term $\nabla ·\left(\frac{\nabla \varphi }{\left|\nabla \varphi \right|}\right)$ in 3D based on the finite difference scheme is necessary for evolving the level set function $\varphi$ in (8). The 3D grid space of dimension $MNL$ is represented by the symbol $\Omega \to {\mathbb{R}}^{MNL}$. The following formula can be used to calculate the second-order coupled curvature term $\nabla ·\left(\frac{\nabla \varphi }{\left|\nabla \varphi \right|}\right)$ at the voxel $\left(i,j,k\right)$:

$$\begin{array}{c}\hfill \frac{{\nabla }_x^{+}{\varphi }_{i,j,k}}{\sqrt{{\left({\nabla }_x^{+}{\varphi }_{i,j,k}\right)}^2+{\left({\nabla }_y^0{\varphi }_{i+\frac{1}{2},j,k}\right)}^2+{\left({\nabla }_z^0{\varphi }_{i,j,k+\frac{1}{2}}\right)}^2+{\epsilon }^2}}-\frac{{\nabla }_x^{-}{\varphi }_{i,j,k}}{\sqrt{{\left({\nabla }_x^{-}{\varphi }_{i,j,k}\right)}^2+{\left({\nabla }_y^0{\varphi }_{i,j-\frac{1}{2},k}\right)}^2+{\left({\nabla }_z^0{\varphi }_{i,j,k-\frac{1}{2}}\right)}^2+{\epsilon }^2}}\\ \hfill +\frac{{\nabla }_y^{+}{\varphi }_{i,j,k}}{\sqrt{{\left({\nabla }_y^{+}{\varphi }_{i,j,k}\right)}^2+{\left({\nabla }_x^0{\varphi }_{i+\frac{1}{2},j,k}\right)}^2+{\left({\nabla }_z^0{\varphi }_{i,j,k+\frac{1}{2}}\right)}^2+{\epsilon }^2}}-\frac{{\nabla }_y^{-}{\varphi }_{i,j,k}}{\sqrt{{\left({\nabla }_y^{-}{\varphi }_{i,j,k}\right)}^2+{\left({\nabla }_x^0{\varphi }_{i-\frac{1}{2},j,k}\right)}^2+{\left({\nabla }_z^0{\varphi }_{i,j,k-\frac{1}{2}}\right)}^2+{\epsilon }^2}}\\ \hfill +\frac{{\nabla }_z^{+}{\varphi }_{i,j,k}}{\sqrt{{\left({\nabla }_z^{+}{\varphi }_{i,j,k}\right)}^2+{\left({\nabla }_x^0{\varphi }_{i+\frac{1}{2},j,k}\right)}^2+{\left({\nabla }_y^0{\varphi }_{i,j+\frac{1}{2},k}\right)}^2+{\epsilon }^2}}-\frac{{\nabla }_z^{-}{\varphi }_{i,j,k}}{\sqrt{{\left({\nabla }_z^{-}{\varphi }_{i,j,k}\right)}^2+{\left({\nabla }_x^0{\varphi }_{i-\frac{1}{2},j,k}\right)}^2+{\left({\nabla }_y^0{\varphi }_{i,j-\frac{1}{2},k}\right)}^2+{\epsilon }^2}}\end{array},$$

(16)

with the small positive value $\epsilon$ preventing division by zero. To maintain rotation-invariant properties, the half-point difference technique is used in (16). At the voxel $(i,j,k)$, the first-order forward ${\partial }_x^{+}$ and backward ${\partial }_x^{-}$ discrete derivatives can be defined as follows:

$${\nabla }_x^{+}{\varphi }_{i,j,k}={\varphi }_{i+1,j,k}-{\varphi }_{i,j,k},{\nabla }_x^{-}{\varphi }_{i,j,k}={\varphi }_{i,j,k}-{\varphi }_{i-1,j,k}$$

$${\nabla }_y^{+}{\varphi }_{i,j,k}={\varphi }_{i,j+1,k}-{\varphi }_{i,j,k},{\nabla }_y^{-}{\varphi }_{i,j,k}={\varphi }_{i,j,k}-{\varphi }_{i,j-1,k}$$

$${\nabla }_z^{+}{\varphi }_{i,j,k}={\varphi }_{i,j,k+1}-{\varphi }_{i,j,k},{\nabla }_z^{-}{\varphi }_{i,j,k}={\varphi }_{i,j,k}-{\varphi }_{i,j,k-1}$$

4. Evaluation Metrics

Compared with the traditional evaluation metrics, which sweep all the points to calculate the closest distance, we divide the original point clouds into groups using the Delaunay Algorithm. In this way, only the closest tetrahedron $T_i$ needs to be found, and the nearest four points need to be chosen in this tetrahedron. Given a set of discrete points, and four points that are not on the same line can determine a sphere. If all the other points are at the exterior of the sphere, the tetrahedron, which is composed of these four conditional points, is called the Delaunay tetrahedron.

$R{T}_{ij}$ is the distance between $Q_j$ and $T_i$. For expedited determination of the nearest tetrahedron, the calculation of $R{T}_{ij}$ can be performed using the following equation:

$$R{T}_{ij}=\frac{r_{j1}+r_{j2}+r_{j3}+r_{j4}}{4}$$

(17)

while $i=\left\{1,2,\cdots ,m\right\},j=\left\{1,2,\cdots ,n\right\}$,and $\left\{r_{j1},r_{j2},r_{j3},r_{j4}\right\}$ is the distance from $Q_j$ to the point in $\left\{P_{i1},P_{i2},P_{i3},P_{i4}\right\}$. $R{T}_{cj}$ is the shortest distance between $Tc\in \left\{T1,\cdots ,Tm\right\}$ and $Q_j$ in (18):

$$R{T}_{cj}=\underset{j}{min}\left(R{T}_{ij}\right)$$

(18)

We define the $R{P}_j$ as the closest Euclidean Metric in Delaunay tetrahedron $Tc$ as follows:

$$R{P}_j=min(r_{j1},r_{j2},r_{j3},r_{j4})$$

(19)

Conclusively, the summation of $R{P}_c$ is computed to derive the error E for the two point clouds as expressed in the equation:

$$E=\sum _{j=1}^{n}R{P}_j$$

(20)

5. Results

We divide the experimental results into two subsections. The first is the case where the model is applied to 2D images, and then the reconstruction of 3D models is discussed. We mainly compare the different models to analyze the efficacy of the proposed model.

5.1. Experiments for 2D Case

The proposed variational method is applied to obtain 2D results. In Figure 4, a 2D contour of the data set shown in Figure 1a is displayed, while 2–7 lines of Figure 4 show the surface reconstruction results of the Chan–Vese model, GAC model, TV model, TV${}_g$-L1 model, TV${}_g$-B model, and the proposed method, respectively. The initial contour is very similar to the true point clouds, which can hasten the convergence process. We assess the reconstruction quality using the uniform sampling error, which involves converting the reconstructed surface into point cloud data and computing the average of the minimum Euclidean distances between the reconstructed point cloud and the original point cloud data. From

Table 1. Uniform sampling error in different methods.

	Chan–Vese	GAC	TV	TV${}_{\mathit{g}}$-L1	TV${}_{\mathit{g}}$-B	Ours
Point Cloud 1	1.5337	1.5926	1.4644	1.3292	1.1669	1.1278
Point Cloud 2	1.4936	1.7140	2.1126	1.6273	1.2777	1.2256

and Figure 4, we can find that although the reconstruction produced by the Chan–Vese model loses small details of the original data as it evolves, the GAC model yields better results. The TV model and the TV${}_g$-L1 model do not use the level set function, so their results are just rough outlines. However, the term (5) in the TV${}_g$-B method and the proposed variational model are able to preserve the original features, and the level set method has good topology; additionally, the proposed method achieves a slightly higher reconstruction accuracy compared to the TV${}_g$-B method.

The results show that proposed model performs well in 2D reconstruction under various conditions. Next, we present the 3D cases.

5.2. Experiments for 3D Case

In order to demonstrate the superiority of the proposed Euler elastica energy reconstruction model in 3D point cloud reconstruction, we conducted numerical experiments using both synthetic and real data sets. First, in the synthetic experiments, as shown in Figure 5, we aimed to reconstruct a spherical object. We compared proposed Euler’s elastica model with the classical TV reconstruction model. From the experimental results, it is evident that the curves reconstructed using the TV regularization term are affected by the staircase effect, while the surfaces reconstructed using Euler’s elastica as regularization exhibit greater smoothness. This indicates that Euler’s elastica model performs exceptionally well in surface reconstruction.

The results of the 3D surface reconstruction of the bunny point cloud data are shown in Figure 6 using six different methods. The six models all have fitting terms so they all performed fairly well. The contour generated by the Chan–Vese model is not as smooth as that generated by the GAC model, but the latter loses more contour. The TV model and the TV${}_g$-L1 model also lose lots of features and exhibit shrinking behavior. The results of the TV${}_g$-B method and the proposed model are better than the above four models, but the results of the TV${}_g$-B model in reconstruction often suffer from excessive smoothness, leading to significant contrast loss, loss of texture details, and noticeable staircase effects, particularly in regions such as the legs. The proposed method, on the other hand, outperforms the TV${}_g$-B model in preserving fine leg details and texture while minimizing the presence of staircase effects. This shows that the proposed model is the best at preserving details.

In Figure 7, the 2nd–7th lines show that the Chan–Vese model, the GAC model, the TV model, and the TV${}_g$-L1 model all cannot retain some fine details and texture, as those methods all stopped at the local extreme such that the fingers failed to be reconstructed well, and both the TV${}_g$-B method and the proposed methods had no such problems. In that regard, the TV${}_g$-B method and the proposed method performed better than any other method. Nevertheless, our approach demonstrates a lower level of reconstruction accuracy in the connecting region between the index finger and middle finger compared to the TV${}_g$-B method. This discrepancy is primarily attributed to the utilization of non-linear curve interpolation in Euler’s elastica for point cloud surface reconstruction. The non-linear interpolation method exhibits inferior performance when handling curved local extrema in comparison to the piecewise linear TV${}_g$-B method.

Figure 8 shows the convergence curve of the bunny and hand by the Chan–Vese model, the GAC model, the TV model, the TV${}_g$L-1 model, the TV${}_g$-B method, and the proposed model. We plot the normalized energy for clarity. The overall results show that all methods are convergent. In

Table 2. Comparison of the number of iterations and running time required for the energy functions of different methods to reach the convergence condition (15).

	Chan–Vese		GAC		TV		TV${}_{\mathit{g}}$-L1		TV${}_{\mathit{g}}$-B		Ours		ADMM
	Iter	Times	Iter	Times	Iter	Times	Iter	Times	Iter	Times	Iter	Times	Iter	Times
Bunny	43rd	2.45 s	3rd	14.14 s	46th	22.51 s	28th	5.23 s	28th	9.19 s	7th	1.97 s	29th	24.42 s
Hand	43rd	12.04 s	3rd	77.46 s	30th	56.67 s	23rd	22.11 s	13th	22.59 s	17th	24.16 s	20th	103.71 s

, we list the iteration count and runtime required for different models to reach the convergence criterion (15). From the table, it is evident that proposed model performs comparably to the TV${}_g$-B model in both convergence speed and runtime. This can be mainly attributed to the utilization of the fast dual projection algorithm presented in Equation (11). In contrast to the ADMM algorithm shown in (10), which requires computations of five variables and four Laplacian operators at each iteration, our proposed fast method only necessitates the calculation of two variables. Using the ADMM algorithm for reconstructing the bunny and hand requires more running times. Additionally, the ADMM method introduces four new parameters ${\gamma }_1$, ${\gamma }_2$, ${\gamma }_3$, and ${\gamma }_4$, which need to be carefully adjusted during the solving process. In contrast, our proposed method not only exhibits higher efficiency but it also significantly reduces the application complexity. With fast dual projection algorithm, the number of iterations and runtime is reduced, resulting in a more efficient reconstruction process. Moreover, by introducing only one dual variable $\overrightarrow{p}$, we alleviate the burden of adjusting multiple parameters, streamlining the usage of the model in practical applications. As a result, the computational cost per iteration is significantly reduced, allowing the efficient solving of complex Euler’s elastic models to be on par with lower-order TV models in terms of computational efficiency.

Based on the results of Figure 7, we can specifically check the effect of the balloon force by changing the value of $\beta$. As shown in Figure 9, it can be found that it is impossible to cross a local extremum when $\beta$ is too small. With the increase in $\beta$, the result is more accurate, and the number of iterations required is correspondingly reduced. In addition, the number of iterations required by the proposed variational method is always less than that of the TV${}_g$-B method. However, $\beta$ cannot grow indefinitely; once the $\beta$ is too large, the reconstruction will fail.

Table 3. Delaunay tetrahedron error in different iterations.

	10th Error	20th Error	30th Error	40th Error	50th Error
Bunny	41.36	40.90	39.19	39.32	39.34
Hand	67.20	67.15	67.12	67.10	67.09

shows the normalized error corresponding to the different iterations of bunny and hand reconstruction that show the relationship between the iteration and the Delaunay error. It proves that the greater the iteration, the closer the reconstructed hand point clouds are to the original model. The errors for the bunny and the hand model are different. The bunny reconstruction is significantly better than the hand, owing to the fingers’ details and concave regions, which greatly increase the Delaunay error of the hand model.

To better demonstrate the effectiveness of the Euler elastica model compared to the TV model in the 3D reconstruction process, we magnified local images of the reconstruction results as shown in Figure 10. In the reconstruction process of the point clouds “bunny” and “hand”, this figure clearly demonstrates that the TV model exhibits noticeable staircase artifacts in the reconstructed legs of the “bunny” and the pinky of the “hand”. In contrast, the incorporation of Euler’s elastica as a regularization term successfully overcomes this issue. Furthermore, in Figure 10 we also zoomed in on the left ear of the bunny and the index and middle fingers of the hand. As can be seen, TV leads to a less accurate surface reconstruction from the respective point cloud. In contrast, the utilization of Euler’s elastica yields more precise and faithful reconstructions.

Figure 11 shows the Delaunay tetrahedron error for the bunny and hand using six different methods. It is clear that the Chan–Vese model, GAC model, TV model and TV${}_g$-L1 model all move towards stabilization, but their results are not as good as those of the TV${}_g$-B method (3) and the proposed elastica method (11).

When the number of iterations reaches 50, the result of the reconstruction shows a tendency toward stabilization. In

Table 4. Normalized Delaunay tetrahedron error in different methods.

	Chan–Vese	GAC	TV	TV${}_{\mathit{g}}$-L1	TV${}_{\mathit{g}}$-B	Ours
Bunny	0.91	0.92	1.00	0.98	0.46	0.00
Hand	1.00	0.94	0.74	0.83	0.40	0.00
Angel	0.64	0.47	0.81	1.00	0.32	0.00
Armadillo	0.66	0.40	0.14	0.36	1.00	0.00
Buda	0.66	0.75	1.00	0.83	0.47	0.00
Dragon	0.68	0.58	0.79	0.74	1.00	0.00
Head	0.88	0.78	1.00	0.85	0.51	0.00
Horse	0.36	0.31	0.95	1.00	0.00	0.01
Lucy	1.00	0.00	0.39	0.46	0.11	0.40

, it can be seen that different models have different reconstruction errors. In order to better compare the accuracy of the reconstruction, we normalized the results to better present them. By setting the maximum error to 1 and the minimum error to 0, we can represent the errors between different algorithms on a scale of 0 to 1. This allows for a more intuitive observation of the reconstruction results, as the errors are normalized within a standardized range. The Chan–Vese model and GAC model still perform as normal, and in some point clouds with large gradient changes, TV and TV${}_g$-L1 have a lower error value, but these four models do not perform well in the concave regions. However, the TV${}_g$-B method and the proposed elastica model are more stable in that they can achieve a relatively small error in diverse point clouds.

As seen in Figure 12, the surface reconstruction results from point clouds with more intricate shapes are accurately handled by the TV${}_g$-B method and the newly suggested variational method. These two methods produce smooth and visually appealing reconstructed surfaces due to the utilization of a continuous signed distance function. In contrast, when binary-based segmentation models, such as the TV model and TV${}_g$-L1 model, are implemented, the reconstructed shapes are not as smooth as those produced by the proposed method. This difference is due to the discontinuity present in the binary representation compared to the continuous level set function used in the proposed method.

Figure 13 shows the main difference between the results of the TV${}_g$-B method and the proposed variational method. Both have their advantages and disadvantages. The TV${}_g$-B method produced smoother results than our method but led to the loss of details on the curved surface. As shown by the area marked by the red square in Figure 13, the TV${}_g$-B method failed to restore the little details, whereas our method performed better. The better performance in the reparation of details can also be seen in Figure 12 in the 2nd and 6th columns.

In order to further investigate the applicability and effectiveness of our proposed method, Figure 14 displays the reconstructed facial surfaces from face point cloud data. From the experimental results, it is easy to observe that the proposed Euler elastica model can effectively restore the facial surfaces. Moreover, compared to existing TV${}_g$-B algorithms, our method preserves more details in the reconstructed surfaces (including the cheek and chin areas), further emphasizing the superiority of the proposed approach in reconstruction.

6. Conclusions

In this paper, we incorporate the Euler elastica regularization term for 3D point cloud reconstruction, aiming to address issues encountered in existing TV models, such as the stair-step effect, contrast loss, and excessive smoothness. Additionally, our objective is to tackle the computational challenges associated with higher-order models and the inefficiency of the ADMM algorithm. The ADMM algorithm often introduces multiple penalty parameters and augmented Lagrange multipliers, and we seek to improve its efficiency in the context of our proposed approach. The introduction of numerous variables in the model leads to decreased computational efficiency, and the selection of penalty parameters becomes challenging without specific guidelines. Therefore, we propose a fast dual projection algorithm for this model, which introduces only one dual variable and one penalty parameter, significantly reducing the difficulty of parameter tuning compared to the ADMM algorithm and greatly accelerating the model computation speed. Experimental results demonstrate that Euler’s elastica, as a regularization term, outperforms the TV regularization term in preserving contrast and avoiding the stair-step effect. However, it exhibits comparatively less capability in handling local extrema than the TV regularization term. In our future work, we plan to explore the use of landmark-based methods to address the limitations of Euler’s elastica regularization in reconstructing local extremum regions.