AN EXPLICIT NUMERICAL ALGORITHM FOR SURFACE RECONSTRUCTION FROM UNORGANIZED POINTS USING GAUSSIAN FILTER

Abstract

We present an explicit numerical algorithm for surface reconstruction from unorganized points using the Gaussian filter. We construct a surface from unorganized points and solve the modified heat equation coupled with a fidelity term which keeps the given points. We apply the operator splitting method. First, instead of solving the diffusion term, we use the Gaussian filter which has the effect of diffusion. Next, we solve the fidelity term by using the fully implicit scheme. To investigate the proposed algorithm, we perform computational experiments and observe good results.

1. INTRODUCTION

An important work in the surface reconstruction process is to preserve the original image and achieve the suitable smoothing effect. There are a number of papers that have performed the surface smoothing, the image processing, the surface or volume reconstruction by applying various techniques such as Gaussian smoothing [1], spherical diffusion [2], anisotropic diffusion [3], non-iterative feature preserving filtering [4], minimizing total variation [5, 6], level-set method [7], embedding narrow band [8, 9], and phase-field method [10].

The following model is used for removing noise from given images, which based on minimizing total variation of image [11]:

\(\frac{\partial \phi}{\partial t}=\nabla \cdot\left(\frac{\nabla \phi}{|\nabla \phi|}\right)+\lambda(f(\mathbf{x})-\phi), \quad \mathbf{x} \in \Omega, t>0,\)       (1.1)

where ϕ(x, t) is an image to be restored in a domain Ω, f(x) is the original image, and λ is a nonzero Lagrange multiplier.

In Eq. (1.1), since the term ∇ ·(∇ϕ/|∇ϕ|) describes the effect of diffusion, we equivalently consider the following evolutionary governing equation:

\(\frac{\partial \phi(\mathbf{x}, t)}{\partial t}=\Delta \phi(\mathbf{x}, t)+\lambda(f(\mathbf{x})-\phi(\mathbf{x}, t)), \quad \mathbf{x} \in \Omega\)       (1.2)

where ϕ(x, t) is the phase-field, f(x) is the original data set, and λ is the fidelity parameter. We define ϕ = 1 on the points in the domain Ω and ϕ = −1 otherwise. Assuming that the reconstructed surface is sufficiently away from the boundary, the Dirichlet boundary condition, ϕ = −1, can be used [12].

The main purpose of this article is to present an explicit numerical algorithm for surface reconstruction from unorganized points by using the Gaussian filter.

This article is organized the following manner. In Section 2, we describe a mathematical configuration and numerical solution algorithm using the Gaussian filter. To investigate the proposed algorithm, we perform computational experiments in Section 3. Conclusions are derived in Section 4.

2. MATHEMATICAL CONFIGURATION AND NUMERICAL SOLUTION ALGORITHM

In this section, we present an explicit numerical algorithm for surface reconstruction from unorganized points using the Gaussian filter. Let the computational domain Ω be (a, b) × (c, d)×(e, f), the number of grid points N_x, N_y, and N_z be positive integers, and the uniform grid size h be (b − a)/(Nx − 1) = (d − c)/(N_y − 1) = (f − d)/(N_z − 1). We define a discrete domain Ω_h = {x = (x_i , y_j , z_k) : x_i = a + (i − 0.5)h, yj = c + (j − 0.5)h, z_k = e + (k − 0.5)h, 1 ≤ i ≤ N_x, 1 ≤ j ≤ N_y, 1 ≤ k ≤ N_z}. The notation ϕ(x_i , y_j , z_k, n∆t) is simply denoted by \(ϕ^n_{ijk}\),, where ∆t is the temporal step size. By using the operator splitting method, we split the diffusion term and fidelity term of Eq. (1.2) into the following equations:

\(\frac{\partial \phi}{\partial t}=\Delta \phi,\)       (2.1)

\(\frac{\partial \phi}{\partial t}=\lambda(f(\mathbf{x})-\phi).\)       (2.2)

When solving Eq. (1.2), we can replace the diffusion Eq. (2.1) using a Gaussian function [1]. For the sake of clarity, we first explain the function in two-dimensional space:

\(G(x, y)=\frac{1}{2 \pi \sigma^{2}} e^{-\frac{x^{2}+y^{2}}{2 \sigma^{2}}}\)

where σ is the standard deviation of distribution. Figure 1 shows a graphical representation of two-dimensional Gaussian function with σ = 1.0, which resembles a bell shape. As the first step to find the numerical solution, instead of solving the diffusion Eq. (2.1), we use the numerical convolution [13] G ∗ \(\phi_{i j}^{0}\) with an initial condition \(\phi_{i j}^{0}\).

 

FIGURE 1. Two-dimensional Gaussian function with σ = 1.0.

Now, we consider (m × m) smoothing kernel. Here, m is an odd integer greater than 2. We define W_pq = G((p − (m + 1)/2)h,(q − (m + 1)/2)h) for 1 ≤ p, q ≤ m. In general, the value of sum of weighting coefficients holds \(\sum_{p=1}^{m} \sum_{q=1}^{m}\) Wpq ̸= 1. Therefore, we normalize the weighting coefficients to make the value of sum equal to be 1, i.e., wpq = Wpq/\(\sum_{i=1}^{m} \sum_{j=1}^{m}\)W_ij. Figure 2 shows that the shape of wpq is smoother as the value of σ increases when m = 3.

 

FIGURE 2. When m = 3, (a), (b), and (c) are the graphical representations of Gaussian filter with σ = 0.5, 1.0, and 1.5, respectively.

The solution of Eq. (2.1) obtained by using the Gaussian filter is as follows:

\(\phi_{i j}^{*}=\sum_{p=1}^{m} \sum_{q=1}^{m} w_{p q} \phi_{i-\frac{m+1}{2}+p, j-\frac{m+1}{2}+q}^{n}, \quad 1 \leq i \leq N_{x}, 1 \leq j \leq N_{y},\)       (2.3)

which has the effect of diffusion.

When m = 3, Fig. 3 illustrates the stencil of ϕ_{i−(m+1)/2+p,j−(m+1)/2+q} and normalized weighting coefficient w_pq for 1 ≤ p, q ≤ m. Similarly to the two-dimensional Gaussian function, the three-dimensional Gaussian function is expressed as

 

FIGURE 3. Schematic illustration of ϕ_i-2+p, j-2+q and normalized weighting coefficient w_pq for 1 ≤ p, q ≤ 3.

\(G(x, y, z)=\frac{1}{(\sqrt{2 \pi} \sigma)^{3}} e^{-\frac{x^{2}+y^{2}+z^{2}}{2 \sigma^{2}}}.\)

For 1 ≤ i ≤ N_x, 1 ≤ j ≤ N_y, 1 ≤ k ≤ N_z, we can extend the two-dimensional numerical solution \(\phi_{ij}^{∗}\) in Eq. (2.3) to the following three-dimensional numerical solution \(\phi_{ijk}^{∗}\):

\(\phi_{i j k}^{*}=\sum_{p=1}^{m} \sum_{q=1}^{m} \sum_{r=1}^{m} w_{p q r} \phi_{i-\frac{m+1}{2}}^{n}+p, j-\frac{m+1}{2}+q, k-\frac{m+1}{2}+r.\)

Figure 4 illustrates stencil of ϕ_i−2+p,j−2+q,k−2+r and normalized weighting coefficient w_pqr.

 

FIGURE 4. Schematic illustration of ϕ_i-2+p, j-2+q, k-2+r and normalized weighting coefficient w_pqr for 1 ≤ p, q, r ≤ 3.

As the second step, we solve Eq. (2.2) by using the fully implicit scheme:

\(\frac{\phi_{i j k}^{n+1}-\phi_{i j k}^{*}}{\Delta t}=\lambda\left(f-\phi_{i j k}^{n+1}\right).\)       (2.5)

To simplify the representation, we take θ = 1/(1 + ∆tλ). Then Eq. (2.5) is given as

\(\phi_{i j k}^{n+1}\) = \(\theta \phi_{i j k}^{*}\)+ (1 − θ)f, 0 < θ < 1

Since f is the initial data \(\phi_{ijk}^{0}\), we can straightforwardly find the numerical solution as follows:

\(\phi_{i j k}^{n+1}=\theta \phi_{i j k}^{*}+(1-\theta) \phi_{i j k}^{0}.\)       (2.6)

As a consequence, an explicit numerical algorithm consists of Eqs. (2.4) and (2.6). Therefore, the implicit numerical solver such as multigrid solver is not needed to find the numerical solution. Also, the execution of the algorithm is fast.

3. COMPUTATIONAL EXPERIMENTS

In this section, we perform numerical tests using the proposed numerical method. Figure 5 shows surface reconstruction of the Stanford bunny [14] on the domain Ω = (0, 1) × (0, 1) × (0, 1). We take the following computational parameters: N_x = N_y = N_z = 60, h = 1/60, ∆t = h², σ = h, and λ = 700. From left to right, the given unorganized points set, initial reconstruction, and the numerical results after 5 iterations.

 

FIGURE 5. Surface reconstruction of bunny: (a) the given points set, (b) initial reconstruction, and (c) the numerical results after 5 iterations, respectively.

3.1. Effect of smoothing kernel. On the three-dimensional space Ω = (0, 1) × (0, 1.5) × (0, 1), we simulate the effect of smoothing kernel. The initial data is set to the Stanford dragon model [14]. We take the following computational parameters: N_x = 120, N_y = 180, N_z = 120, h = 1/120, ∆t = 0.05h, σ = 2h/\(\sqrt{3}\), and λ = 3000. Figure 6(a) shows the initial condition. At the final time T = 5∆t, the results of 3 × 3 × 3 and 5 × 5 × 5 smoothing kernel are shown in Fig. 6(b) and (c), respectively. We can see that 5 × 5 × 5 smoothing kernel are smoother than 3 × 3 × 3 smoothing kernel. It means that the initial data is lost. Therefore, all computational experiments shall use the Gaussian filter computed by 3 × 3 × 3 smoothing kernel, from now on.

 

FIGURE 6. Right and left views of Stanford dragon: (a) initial point set, (b) initial condition (c) 3×3×3, and (d) 5×5×5 smoothing kernel, respectively.

3.2. Effect of standard deviation σ. In this section, we confirm the effects of the standard deviation σ on the three-variable Gaussian filter dynamics. The initial data is taken from the Stanford bunny model. We can observe the shape of the model at time T = 5∆t with respect to different standard deviations. In Fig. 7, the parameters are given as N_x = 120, N_y = 180, N_z = 120, h = 0.0147, ∆t = 0.05h, λ = 0, and Ω = (0, 1) × (0, 1.5) × (0, 1). As shown in Fig. 7, when the standard deviation of the given Gaussian filter is small, it is seen that the resulting surface closely resembles the initial model. On the other hand, when the standard deviation of the filter becomes larger, the resulting surface is smoother. Thus, the standard deviation of the filter must be set to an appropriate value for each set of initial data.

 

FIGURE 7. The graphical views of Stanford bunny: (a) shows the initial condition and (b), (c), and (d) are the numerical results with σ = 0:01h; h, and 100h, respectively.

3.3. Effect of parameter λ. From now on, we investigate the effect of parameter λ on the three-variable Gaussian filter dynamics. The initial condition is used as Stanford bunny and dragon model. We can observe the shape of model at the total time T = 5∆t with respect to different values of λ. The common performance parameters are taken as following: N_x = 120, N_y = 180, N_z = 120, ∆t = 0.05h, σ = 2h/\(\sqrt{3}\), and Ω = (0, 1) × (0, 1.5) × (0, 1), where h = 0.01 in Fig. 8 and h = 0.0147 in Fig. 9. As shown in Fig. 8, 9, when the value of parameter λ is small, we can see that the result is oversmoothing by the Gaussian filter with a sufficiently large standard deviation σ. Meanwhile, if the value of parameter λ is large, then the results are rougher, which seems like an initial condition. Therefore, λ must be set to a suitable value for each given initial data.

 

FIGURE 8. The graphical views of Stanford dragon: (a) shows the initial condition and (b), (c), and (d) are the numerical results with λ = 1; 1000, and 10000, respectively.

 

FIGURE 9. The graphical views of Stanford bunny: (a) shows the initial condition and (b), (c), and (d) are the numerical results with λ = 1; 1000, and 10000, respectively.

4. CONCLUSION

In this article, we presented an explicit numerical algorithm for surface reconstruction from unorganized points by using the Gaussian filter. We have investigated the diffusion effect of the Gaussian filter through the computational experiments and observed the numerical results. The proposed algorithm can be utilized in many aspects of industry, such as the three-dimensional structure printing from scattered scanned data.