R/mesh_repair.R

In espadon: Easy Study of Patient DICOM Data in Oncology

#' Repair of a mesh
#' @description The \code{mesh.repair} function repairs holes in a \code{mesh} 
#' class object.
#' @param mesh "mesh" class object.
#' @param verbose Boolean, by default set to \code{FALSE}. Allows you to inhibit 
#' comments.
#' @return Returns a mesh, repaired by removing degenerated triangles and filling holes.
#' @examples
#' # loading of toy-patient objects (decrease dxyz for better result)
#' step <- 4
#' patient <- toy.load.patient(modality = c("ct", "rtstruct"), 
#'                             roi.name = "gizzard", dxyz = rep(step, 3))
#' CT <- patient$ct[[1]]
#' S <- patient$rtstruct[[1]]
#' 
#' # creation of the gizzard mesh
#' bin <- bin.from.roi(CT, struct = S, roi.name = "gizzard", verbose = FALSE)
#' mesh.gizzard <- mesh.from.bin(bin, alias = "gizzard", verbose = FALSE)
#' 
#' repair.mesh.gizzard <- mesh.repair(mesh.gizzard, verbose = FALSE)
#' str(mesh.gizzard)
#' str(repair.mesh.gizzard)

#' @importFrom Rvcg vcgBorder vcgIsolated vcgSmooth vcgVFadj
#' @importFrom methods is
#' @export
mesh.repair <- function(mesh, verbose = TRUE) {
  
  if (is.null(mesh)) return(NULL)
  if (!is(mesh, "mesh")) stop("mesh should be a mesh class object.")

  mesh_ <- mesh
  mesh <- mesh_$mesh
  ##
  ## Removal of degenerate triangles
  ##
  
  if (verbose) cat("mesh.repair INFO : removing degenerated triangles\n")
  repeat {
    ## connection lists
    f.list <- vcgVFadj(mesh) ## liste des triangles connectés à chaque vertex
    v.list <- lapply(1:length(f.list), function(v) { ## liste des vertex connectés à chaque vertex
      sel <- unique(as.vector(mesh$it[,f.list[[v]]]))
      return(list(v = v, n.v = sel[sel != v]))
    })
    v <- sort(which(sapply(v.list, function(v) length(v$n.v)) == 2), decreasing = TRUE)[1]
    if (is.na(v)) break
    
    if (verbose) cat("mesh.repair INFO : cleaning", v, "\n")
    mesh$vb <- mesh$vb[, -v]
    t <- which(mesh$it == v, arr.ind = TRUE)[,2]
    mesh$it <- mesh$it[, -t]
    sel <- which(mesh$it > v)
    mesh$it[sel] <- mesh$it[sel] - 1
    mesh$normals <- mesh$remvert <- NULL
  }
  mesh <- vcgIsolated(mesh, silent = !verbose)
  
  if (verbose) cat("mesh.repair INFO : removing 3 borders triangles\n")
  borders <- vcgBorder(mesh)$bordervb
  t.b <- apply(mesh$it, 2, function(t) !all(borders[t]))
  if (verbose) cat("mesh.repair INFO : number of tr :", ncol(mesh$it) - sum(t.b) ,"\n")
  mesh$it <- mesh$it[, t.b]
  mesh$normals <- mesh$remvert <- NULL
  mesh <- vcgIsolated(mesh, silent = !verbose)

  if (verbose) cat("mesh.repair INFO : removing >2 borders border vertex\n")
  
  repeat {
    f.list <- vcgVFadj(mesh) ## list of triangles connected to each vertex
    v.list <- lapply(1:length(f.list), function(v) { ## list of vertices connected to each vertex
      sel <- unique(as.vector(mesh$it[,f.list[[v]]]))
      return(list(v = v, n.v = sel[sel != v]))
    })
    
    borders <- vcgBorder(mesh)$bordervb
    b.count <- sapply(which(borders), function(v) {
      sum(borders[v.list[[v]]$n.v])
    })
    
    v <- sort(which(borders)[b.count > 2], decreasing = TRUE)[1]
    if (is.na(v)) break
    
    if (verbose) cat("mesh.repair INFO : cleaning", v, "\n")
    mesh$vb <- mesh$vb[, -v]
    t <- which(mesh$it == v, arr.ind = TRUE)[,2]
    mesh$it <- mesh$it[, -t]
    sel <- which(mesh$it > v)
    mesh$it[sel] <- mesh$it[sel] - 1
    mesh$normals <- mesh$remvert <- NULL
  }
  mesh <- vcgIsolated(mesh, silent = !verbose)
  
  
  ##
  ## repairing holes
  ##
  
  if (verbose) cat("mesh.repair INFO : mesh holes filling\n")
  
  ## listes de connexions
  f.list <- vcgVFadj(mesh) 
  v.list <- lapply(1:length(f.list), function(v) { 
    sel <- unique(as.vector(mesh$it[,f.list[[v]]]))
    return(list(v = v, n.v = sel[sel != v]))
  })
  # table(sapply(f.list, function(v) length(v)))
  # table(sapply(v.list, function(v) length(v$n.v)))
  
  repeat {
    borders <- vcgBorder(mesh)$bordervb
    if (verbose) cat("mesh.repair INFO : remaining borders :", sum(borders), "\n")
    b.list <- which(borders)[1]
    if (is.na(b.list)) break
    repeat {
      if (rev(b.list)[1] %in% v.list[[b.list[1]]]$n.v & length(b.list) > 2) break
      T1 <- which(mesh$it == b.list[1], arr.ind = TRUE)
      for (Tr in T1[, 2]) {
        id1 <- which(mesh$it[, Tr] == b.list[1])
        id2 <- which(mesh$it[, Tr] != b.list[1] & borders[mesh$it[, Tr]])
        if (length(id2) == 1) {
          if (id1 == 1 & id2 == 3) break
          if (id1 == 2 & id2 == 1) break
          if (id1 == 3 & id2 == 2) break
        } else if (length(id2) == 2) {
          stop()
        }
      }
      b.list <- c(mesh$it[id2, Tr], b.list)
    }
    G <- apply(mesh$vb[, b.list], 1, mean)
    mesh$vb <- cbind(mesh$vb, G)
    dimnames(mesh$vb) <- NULL
    tr <- matrix(c(rep(ncol(mesh$vb), length(b.list)), b.list[-1], b.list[1], b.list), nrow = 3, byrow = TRUE)
    mesh$it <- cbind(mesh$it, tr)
    mesh$normals <- mesh$remvert <- NULL
    
    f.list <- vcgVFadj(mesh)
    v.list <- lapply(1:length(f.list), function(v) {
      sel <- unique(as.vector(mesh$it[,f.list[[v]]]))
      return(list(v = v, n.v = sel[sel != v]))
    })
  }
  
  mesh_$mesh <- mesh
  mesh_$nb.faces <- dim(mesh$it)[2]
  class(mesh_)
  return(mesh_)
}