Decimate 3D surface mesh evenly and preserving bilaterally symmetry

Question

Here is a sample of 3D facial surface mesh. As can been seen from Fig 1 below, the landmarks are bilaterially symmetric. I wish to reduce the number of landmarks.

Here face is the coordinates of the vertices, triang is the triangulation file, landPairs is a two-column file containing vertices pairing information. landPairs is not used for plotting later but is provided in case it is needed. All data are available from here:

Here is the code to plot the original vertices before decimation:

library(rgl)
library(Rvcg)

# Customized function to convert vb and it information to 3D mesh
lm2mesh <- function(vb, it) {
    vb <- t(vb)
    vb <- rbind(vb, 1)
    rownames(vb) <- c("xpts", "ypts", "zpts", "")

    it_mat <- t(as.matrix(it))
    rownames(it_mat) <- NULL

    vertices <- c(vb)
    indices <- c(it_mat)

    tmesh3d(vertices = vertices, indices = indices, homogeneous = TRUE, 
            material = NULL, normals = NULL, texcoords = NULL)
}
# Load `face` and `triang`    
face <- as.matrix(read.csv("<PATH>\\SampleFace.csv", header=F))
triang <- as.matrix(read.csv("<PATH>\\triangulation.csv", header=F))

facemesh <- lm2mesh(face,triang)

# Plot the undecimated mesh
shade3d(facemesh, col="steelblue", specular = "#202020", alpha = 0.7)
plot3d(face, type = "s", col = "red", xlab = "x", ylab = "y", zlab = "z", 
       size = 0.2, aspect = FALSE, alpha = 0.8, add=T)

Here is Fig 1 below (vertices are evenly spaced and are absolutely bilaterally symmetric):

# Plot the decimated mesh
open3d()
facemeshdecim <- vcgQEdecim(facemesh,percent=0.1)
shade3d(facemeshdecim, col="steelblue", specular = "#202020", alpha = 0.7)
plot3d(t(facemeshdecim$vb[-4, ]), type = "s", col = "red", xlab = "x", ylab = "y", zlab = "z", 
       size = 0.4, aspect = FALSE, alpha = 0.8, add=T)

Here is Fig 2 below (vertices are no evenly spaced and are no longer symmetric):

It can be seen that in the decimated face, the vertices are NOT as evenly spaced as before decimation and the originally symmetric vertices became NO LONGER symmetric. My question is if there is a way to reduce the number of vertices while ensuring that the reduced vertices are as evenly spaced as possible AND preserving the bilateral symmetry of the vertices?

Answer 1

Here's a way. Start with your code, then add this:

# Get the positive part of the face
posface <- clipMesh3d(facemesh, fn="y")

# Decimate it, keeping the boundary
posdeci <- vcgQEdecim(posface, percent=0.1, bound = TRUE)

# Duplicate it in a reflection
negdeci <- posdeci
negdeci$vb[2,] <- -negdeci$vb[2,]

# Join them together
fulldeci <- merge(posdeci, negdeci)

# Plot it
open3d()
shade3d(fulldeci, col="steelblue", specular = "#202020", alpha = 0.7)
plot3d(t(fulldeci$vb[-4, ]), type = "s", col = "red", xlab = "x", ylab = "y", zlab = "z", 
   size = 0.4, aspect = FALSE, alpha = 0.8, add=T)

This has too many points along the mid-line, but otherwise does what you want.

Edited to add:

Getting the points more uniform is a little tricky. If you don't use bound = TRUE in the call to vcgQEdecim(), it will leave a gap down the middle of the face. To fill that, you need to add quads joining the two sides of the edge, but figuring out which vertices form the edge needs a new function:

getBorder <- function(mesh) {
  border <- which(vcgBorder(mesh)$bordervb)
  inorder <- NULL
  repeat{
    i <- 1
    inorder <- c(inorder, border[i])
    repeat{
      found <- FALSE
      tris <- which(apply(mesh$it, 2, function(col) border[i] %in% col))
      for (j in tris) {
        tri <- mesh$it[,j]
        i0 <- which(tri == border[i])
        i1 <- i0 %% 3 + 1
        # keep tri[i1] if the edge from tri[i0] to tri[i1] is external
        tris1 <- which(apply(mesh$it[,tris,drop=FALSE], 2, function(col) all(tri[c(i0, i1)] %in% col)))
        if (length(tris1) == 1) {
          if (tri[i1] %in% inorder)
            break
          inorder <- c(inorder, tri[i1])
          i <- which(border == tri[i1])
          found <- TRUE
          break
        }
      }
      if (!found) break
    }
    border <- setdiff(border, inorder)
    if (!length(border)) break
    inorder <- c(inorder, NA)
  }
  inorder
}

Using that function, the following code does a reasonable job:

# Try joining halves using quads

posdeci2 <- vcgQEdecim(posface,percent=0.1, bound = FALSE)
negdeci2 <- posdeci2
negdeci2$vb[2,] <- -negdeci2$vb[2,]

# This one has the gap
fulldeci2 <- merge(posdeci2, negdeci2)

# Fill in the gap with quads
# Keep the ones in the middle, but not the outside edge
border <- getBorder(posdeci2)
border <- border[posdeci2$vb[2, border] < 0.005]

borderverts <- posdeci2$vb[, border]
negverts <- negdeci2$vb[, border]

# The quads have both sets of vertices
quadverts <- cbind(borderverts, negverts)
n <- ncol(borderverts)

# We'll assume n > 1
indices <- rbind(1:(n-1), 2:n, n + 2:n, n + 1:(n-1))
quads <- mesh3d(vertices = quadverts, quads = indices)
fulldeci3 <- merge(fulldeci2, quads)

# plot it
open3d()
shade3d(fulldeci3, col="steelblue", specular = "#202020", alpha = 0.7)
plot3d(t(fulldeci3$vb[-4, ]), type = "s", col = "red", xlab = "x", ylab = "y", zlab = "z", 
       size = 0.4, aspect = FALSE, alpha = 0.8, add=T)