Smoothing

Both Laplacian smoothing and Taubin smoothing^{1 2} are smoothing operations that adjust the positions of the nodes in a finite element mesh.

Laplacian smoothing, based on the Laplacian operator, computes the average position of a point's neighbors and moves the point toward the average. This reduces high-frequency noise, but can result in a loss of shape and detail, with overall shrinkage.

Taubin smoothing is an extension of Laplacian smoothing that seeks to overcome the shrinkage drawback associated with the Laplacian approach. Taubin is a two-pass approach. The first pass smooths the mesh. The second pass re-expands the mesh.

Laplacian Smoothing

Consider a subject node with position $\mathbf{p}=\mathbf{p}(x,y,z)$. The subject node connects to $n$ neighbor points $\mathbf{q}\_i$ for $i \in [1, n]$ through $n$ edges.

For concereteness, consider a node with four neighbors, shown in the figure below.

Figure: The subject node $\mathbf{p}$ with edge connections (dotted lines) to neighbor nodes $\mathbf{q}\_i$ with $i \in [1, n]$ (withouth loss of generality, the specific example of $n=4$ is shown). The average position of all neighbors of $\mathbf{p}$ is denoted $\bar{\mathbf{q}}$, and the gap $\Delta \mathbf{p}$ (dashed line) originates at $\mathbf{p}$ and terminates at $\bar{\mathbf{q}}$.

Define $\bar{\mathbf{q}}$ as the average position of all $\mathbf{q}\_i$ neighbors of $\mathbf{p}$,

$$\bar{\mathbf{q}}:=\frac{1}{n}\sum {i=1}^n\mathbf{q}i.$$

Define the gap vector $\Delta \mathbf{p}$ as originating at $\mathbf{p}$ and terminating at $\bar{\mathbf{q}}$ (viz., $\mathbf{p}+\Delta \mathbf{p}=\bar{\mathbf{q}}$),

$$\Delta \mathbf{p}:=\bar{\mathbf{q}}-\mathbf{p}.$$

Let $\lambda \in {\mathbb{R}}^{+}$ be the positive scaling factor for the gap $\Delta \mathbf{p}$.

Since

$$\bar{\mathbf{q}}=\mathbf{p}+\lambda \Delta \mathbf{p}\mathrm{when}\lambda =1,$$

subdivision of this relationship into several substeps gives rise to an iterative approach. We typically select $\lambda < 1$ to avoid overshoot of the update, $\lambda \in {\mathbb{R}}^{+}\subset (0,1)$.

At iteration $k$, we update the position of ${\mathbf{p}}^{(k)}$ by an amount $\lambda \Delta {\mathbf{p}}^{(k)}$ to ${\mathbf{p}}^{(k+1)}$ as

$${\mathbf{p}}^{(k+1)}:={\mathbf{p}}^{(k)}+\lambda \Delta {\mathbf{p}}^{(k)},$$

with

$$\Delta {\mathbf{p}}^{(k)}={\bar{\mathbf{q}}}^{(k)}-{\mathbf{p}}^{(k)}.$$

Thus

$${\mathbf{p}}^{(k+1)}:={\mathbf{p}}^{(k)}+\lambda \left(\Delta {\mathbf{p}}^{(k)}\right),$$

$${\mathbf{p}}^{(k+1)}:={\mathbf{p}}^{(k)}+\lambda \left({\bar{\mathbf{q}}}^{(k)}-{\mathbf{p}}^{(k)}\right),$$

and finally

$${\mathbf{p}}^{(k+1)}:={\mathbf{p}}^{(k)}+\lambda (\frac{1}{n}\sum {i=1}^n\mathbf{q}{i}^{(k)}-{\mathbf{p}}^{(k)}).$$

The formulation above, based on the average position of the neighbors, is a special case of the more generalized presentation of Laplace smoothing, wherein a normalized weighting factor, $w_i$, is used:

$${\mathbf{p}}^{(k+1)}:={\mathbf{p}}^{(k)}+\lambda \sum {i=1}^n{w}_i(\mathbf{q}{i}^{(k)}-{\mathbf{p}}^{(k)})$$

When all weights are equal and normalized by the number of neighbors, $w_i=\frac{1}{n}$, the special case presented in the box above is recovered.

Example

For a 1D configuration, consider a node with initial position $\mathbf{p}=1.5$ with two neighbors (that never move) with positions $\mathbf{q}\_1=0.0$ and $\mathbf{q}\_2=1.0$ ($\bar{\mathbf{q}}=0.5$). With $\lambda =0.3$, the table below shows updates for for position $\mathbf{p}$.

Table: Iteration updates of a 1D example.

Figure: Convergence of position $\mathbf{p}$ toward $0.5$ as a function of iteration $k$.

Taubin Smoothing

Taubin smoothing is a two-parameter, two-pass iterative variation of Laplace smoothing. Specifically with the definitions used in Laplacian smoothing, a second negative parameter $\mu$ is used, where

$$ \lambda \in \mathbb{R}^+ \subset (0, 1) \hspace{0.5cm} \rm{and} \hspace{0.5cm} \lambda < -\mu. $$

The first parameter, $\lambda$, tends to smooth (and shrink) the domain. The second parameter, $\mu$, tends to expand the domain.

Taubin smoothing is written as, for $k=0$, $k < k_{\rm{max}}$, $k=k+1$, with $k\_\max$ an even number,

• First pass (if $k$ is even):

$${\mathbf{p}}^{(k+1)}:={\mathbf{p}}^{(k)}+\lambda (\frac{1}{n}\sum {i=1}^n\mathbf{q}{i}^{(k)}-{\mathbf{p}}^{(k)}),$$

• Second pass (if $k$ is odd):

$${\mathbf{p}}^{(k+1)}:={\mathbf{p}}^{(k)}+\mu (\frac{1}{n}\sum {i=1}^n\mathbf{q}{i}^{(k)}-{\mathbf{p}}^{(k)}),$$

In any second pass (any pass with $k$ odd), the algorithm uses the updated positions from the previous (even) iteration to compute the new positions. So, the average is taken from the updated neighbor positions rather than the original neighbor positions. Some presentation of Taubin smoothing do not carefully state the second pass update, and so we emphasize it here.

Hierarchical Control

Hierarchical control classifies all nodes in a mesh as belonging to a surface $\mathbb{A}$, interface $\mathbb{B}$, or interior $\mathbb{C}$. These categories are mutually exclusive. Any and all nodes must belong to one, and only one, of these three categories. For a given node $\mathbf{p}$, let

• the set of surface neighbors be denoted $\mathbf{q}\_\mathbb{A}$,
• the set of interface neighbors be denoted $\mathbf{q}\_\mathbb{B}$, and
• the set of interior neighbors be denoted $\mathbf{q}\_\mathbb{C}$.

Hierarchical control states that

• for any surface node $\mathbf{p}\in \mathbb{A}$, only surface nodes $\mathbf{q}\_\mathbb{A}$ connected via edges to $\mathbf{p}$ are considered neighbors,
• for any interface node $\mathbf{p}\in \mathbb{B}$, only surface nodes $\mathbf{q}\ast \mathbb{A}$ and interface nodes $\mathbf{q}\ast \mathbb{B}$ connected via edges to $\mathbf{p}$ are neighbors, and
• for any interior node $\mathbf{p}\in \mathbb{C}$, surface nodes $\mathbf{q}\ast \mathbb{A}$, interface nodes $\mathbf{q}\ast \mathbb{B}$, and interior nodes $\mathbf{q}\_\mathbb{C}$ connecting via edges to $\mathbf{p}$ are neighbors.

The following figure shows this concept:

Figure: Classification of nodes into categories of surface nodes $\mathbb{A}$, interface nodes $\mathbb{B}$, and interior nodes $\mathbb{C}$. Hierarchical relationship: a surface node's neighbors are other other edge-connected surface nodes, an interface node's neighbors are other edge-connected interface nodes or surface nodes, and an interior node's neighbors are edge-connected nodes of any category.

Chen Example

Chen³ used medical image voxel data to create a structured hexahedral mesh. They noded that the approach generated a mesh with "jagged edges on mesh surface and material interfaces," which can cause numerical artifacts.

Chen used hierarchical Taubin mesh smoothing for eight (8) iterations, with $\lambda =0.6307$ and $\mu =-0.6732$ to smooth the outer and inner surfaces of the mesh.

References