Implicitly smooth a triangular mesh — vcg

Applies smoothing algorithms on a triangular mesh. Vertices that belong to no face carry no connectivity, so they are excluded from the computation and returned at their input positions with zero normal vectors; the remaining vertices are unaffected by their presence.

vcg_smooth_implicit solves, separately for each coordinate, the sparse linear system $(M + \lambda L^k) x = M v$, where $v$ holds the input positions, $M$ the per-vertex areas scaled to a maximum of one, $L$ the 'Laplacian' and $k = 2^{degree - 1}$. Degrees 1 and 2 are solved iteratively (conjugate gradients), with memory growing linearly with the mesh; higher degrees are solved by a sparse factorization. The memory the solve needs is estimated before anything large is allocated, and the function stops with an explanation when it exceeds max_memory; very large meshes can be simplified first with vcg_decimate.

Usage

vcg_smooth_implicit(
  mesh,
  lambda = 0.2,
  use_mass_matrix = TRUE,
  fix_border = FALSE,
  use_cot_weight = FALSE,
  degree = 1L,
  laplacian_weight = 1,
  max_memory = 2
)

vcg_smooth_explicit(
  mesh,
  type = c("taubin", "laplace", "HClaplace", "fujiLaplace", "angWeight",
    "surfPreserveLaplace"),
  iteration = 10,
  lambda = 0.5,
  mu = -0.53,
  delta = 0.1
)

Arguments

mesh: triangular mesh stored as object of class 'mesh3d'.
lambda: In vcg_smooth_implicit, the amount of smoothness, which has no effect when use_mass_matrix is FALSE; default is 0.2. In vcg_smooth_explicit, parameter for 'taubin' smoothing.
use_mass_matrix: logical: whether to keep the mesh close to its input position with an area-weighted (mass matrix) term; default is TRUE. With FALSE there is no such term and the fixed border vertices alone determine the result, the smoothest surface spanning the border; this needs fix_border = TRUE and a border in every connected part of the mesh, and lambda then has no effect
fix_border: logical: whether border vertices (vertices on an edge that belongs to exactly one face) keep exactly their input positions; default is FALSE. Closed meshes have no border
use_cot_weight: logical: whether to use cotangent weight; default is FALSE (using uniform 'Laplacian')
degree: integer: degree of the 'Laplacian' operator; the system uses $L^k$ with $k = 2^{degree - 1}$ (the 'Laplacian' squared degree - 1 times), so degree = 2 uses $L^2$ and degree = 3 uses $L^4$; default is 1
laplacian_weight: numeric: weight when use_cot_weight is FALSE; default is 1.0
max_memory: maximum memory in GiB that vcg_smooth_implicit may use for the solve; default is 2. The need is estimated from the mesh before any large allocation, and exceeding it raises an error that suggests alternatives
type: method name of explicit smooth, choices are 'taubin', 'laplace', 'HClaplace', 'fujiLaplace', 'angWeight', 'surfPreserveLaplace'.
iteration: number of iterations
mu: parameter for 'taubin' explicit smoothing.
delta: parameter for scale-dependent 'Laplacian' smoothing or maximum allowed angle (in 'Radian') for deviation between surface preserving 'Laplacian'.

Value

An object of class "mesh3d" with:

vb: vertex coordinates
normals: vertex normal vectors
it: triangular face index

Coercing Surface Inputs

The surface objects are converted to 'mesh3d' object before applying further calculations.

When surface is a surface ieegio object, the returned mesh3d$vb contains vertices that have been left-multiplied by surface$geometry$transforms[[1]] (the first transform stored in the geometry, typically the ScannerAnat or voxel-to-world transform).

Breaking change: Earlier versions (before 0.2.6) of ravetools returned the raw surface$geometry$vertices without applying any transform, so downstream code often multiplied by surface$geometry$transforms[[1]] (or an equivalent) manually before working in world space. Such code will now double apply the transform and produce incorrect coordinates. If you previously applied a transform from surface$geometry$transforms by hand after calling a ravetools mesh function on an 'ieegio_surface', remove that manual step.

Surfaces with an empty or missing geometry$transforms list (for example, surfaces produced by ieegio's volume_to_surface, which stores an identity transform) are unaffected.

If geometry$transforms contains multiple transforms targeting different coordinate spaces, only the first one is used. Callers that need a specific target space should select and apply that transform themselves before calling ravetools mesh functions.

Examples


if(is_not_cran()) {

# Prepare mesh with no normals
data("left_hippocampus_mask")

# Grow 2mm on each direction to fill holes
volume <- grow_volume(left_hippocampus_mask, 2)

# Initial mesh
mesh <- vcg_isosurface(volume)

# Start: examples
rgl_view({
  rgl_call("mfrow3d", 2, 4)
  rgl_call("title3d", "Naive ISOSurface")
  rgl_call("shade3d", mesh, col = 2)

  rgl_call("next3d")
  rgl_call("title3d", "Implicit Smooth")
  rgl_call("shade3d", col = 2,
           x = vcg_smooth_implicit(mesh, degree = 2))

  rgl_call("next3d")
  rgl_call("title3d", "Explicit Smooth - taubin")
  rgl_call("shade3d", col = 2,
           x = vcg_smooth_explicit(mesh, "taubin"))

  rgl_call("next3d")
  rgl_call("title3d", "Explicit Smooth - laplace")
  rgl_call("shade3d", col = 2,
           x = vcg_smooth_explicit(mesh, "laplace"))

  rgl_call("next3d")
  rgl_call("title3d", "Explicit Smooth - angWeight")
  rgl_call("shade3d", col = 2,
           x = vcg_smooth_explicit(mesh, "angWeight"))

  rgl_call("next3d")
  rgl_call("title3d", "Explicit Smooth - HClaplace")
  rgl_call("shade3d", col = 2,
           x = vcg_smooth_explicit(mesh, "HClaplace"))

  rgl_call("next3d")
  rgl_call("title3d", "Explicit Smooth - fujiLaplace")
  rgl_call("shade3d", col = 2,
           x = vcg_smooth_explicit(mesh, "fujiLaplace"))

  rgl_call("next3d")
  rgl_call("title3d", "Explicit Smooth - surfPreserveLaplace")
  rgl_call("shade3d", col = 2,
           x = vcg_smooth_explicit(mesh, "surfPreserveLaplace"))
})

}
#> Package `rgl` is not installed. Please install `rgl` to use this function.
#> Error in loadNamespace(name): there is no package called ‘rgl’