R/dist_mesh.R

In espadon: Easy Study of Patient DICOM Data in Oncology

#' @importFrom graphics boxplot
#' @importFrom rgl open3d points3d wire3d
#' @importFrom grDevices colorRampPalette
#' @importFrom Rvcg vcgClost vcgRaySearch
dist.mesh <- function(mesh.test, mesh.ref, plot=TRUE, 
                      col.in = "blue", col.out = "green", pt.size = 10, in.front = FALSE) {
  
  pt <- t(mesh.test$mesh$vb)
  # it <- t(mesh.test$mesh$it)
  min.max.test <- get.extreme.pt(mesh.test)
  max.min.ref <- get.extreme.pt(mesh.ref)[,2:1]

  tol <- max(sqrt(apply((min.max.test - max.min.ref)^2, 2,sum)))

  m <- vcgClost(mesh.test$mesh, mesh.ref$mesh,tol = tol)
  u <- t(m$vb)[,1:3] - pt[,1:3]
  dum <- u * u
  dum <- sqrt(dum[,1] + dum[,2] + dum[,3])
  fd <- dum != 0
  m$distance <- sign(m$quality) * dum
  s <- vcgArea(mesh.test$mesh, perface = TRUE)$pertriangle
  it.per.vertex <- vcgVFadj(mesh.test$mesh)
  m$surface <- sapply(it.per.vertex, function(v) sum(s[v]) / 3)

  if (in.front) {
    in.front.v <- rep(TRUE, length(m$distance))
    diff  <- rep(0,length(m$distance))
    success <- in.front.v

    m2 <- list()
    m2$vb <- m$vb[,fd]
    m2$normals <- t(as.matrix(cbind(-u[fd,]/dum[fd],1)))
    class(m2) <- "mesh3d"
    d <- vcgRaySearch(m2,mesh.test$mesh)

    diff[fd] <- abs(m$distance[fd]) - d$distance
    success[fd] <- d$quality == 1
    f.nok <- (diff <= -tol |  !success | (diff >= tol & diff <= 1e-2))
    in.front.v[!f.nok] <- abs(diff[!f.nok]) < tol
    if (any(f.nok))
      in.front.v[f.nok] <- .meshinfront( pt_x = as.numeric(pt[,1]),
                                         pt_y = as.numeric(pt[,2]),
                                         pt_z = as.numeric(pt[,3]),
                                         p2_x = as.numeric(t(m$vb)[f.nok,1]),
                                         p2_y = as.numeric(t(m$vb)[f.nok,2]),
                                         p2_z = as.numeric(t(m$vb)[f.nok,3]),
                                         u_x = -as.numeric(u[f.nok,1]),
                                         u_y = -as.numeric(u[f.nok,2]),
                                         u_z = -as.numeric(u[f.nok,3]),
                                         n_A = as.numeric(mesh.test$mesh$it[1,]) - 1,
                                         n_B = as.numeric(mesh.test$mesh$it[2,]) - 1,
                                         n_C = as.numeric(mesh.test$mesh$it[3,]) - 1)
    m$in.front <- in.front.v
  }
  m$quality <- NULL

  if (plot) {
    palet <- colorRampPalette(c(col.in, "grey90",col.out))(100)
    dummy <- factor(sign(m$distance), levels = c(-1, 1), labels = c("in", "out"))
    dev.new(width = 7, height = 5, noRStudioGD = TRUE)
    boxplot(d~in.out, data.frame(d = abs(m$distance), in.out = dummy), col = c(col.in, col.out),
             xlab = paste(mesh.test$object.alias, "vs", mesh.ref$object.alias ), ylab = "d (mm)")
    ## display of measurements
    d.max <- max(abs(m$distance))
    open3d(windowRect = c(100, 100, 600, 600))
    wire3d(mesh.test$mesh, color = "grey90", lit = FALSE)
    if (in.front) {
      f <- (m$in.front == 1)
      points3d(as.matrix(pt[f,1:3]), col = palet[1 + 99 * (m$distance[f] / d.max / 2 + 0.5)], size = pt.size)
    }

    display.palette(palet, breaks = seq(-d.max, d.max, length.out = 101))
    dev.new(width = 7, height = 5, noRStudioGD = TRUE)
    H <- hist(m$distance, breaks = seq(-d.max, d.max, length.out = 100))
  }
  return(m)
}