quantileTPD(x, thresholdPlot=0.99)

quantileTPD <- function(x, thresholdPlot = 0.99) {
  TPDList <- x$TPDc$TPDc
  quantileTPD <- list()
  for (comm in 1:length(TPDList)) {
    percentile <- rep(NA, length(TPDList[[comm]]))
    TPDList[[comm]] <- cbind(index = 1:length(TPDList[[comm]]),
                             prob = TPDList[[comm]], percentile)
    orderTPD <- order(TPDList[[comm]][, "prob"], decreasing = TRUE)
    TPDList[[comm]] <- TPDList[[comm]][orderTPD, ]
    TPDList[[comm]][, "percentile"] <- cumsum(TPDList[[comm]][, "prob"])
    TPDList[[comm]] <- TPDList[[comm]][order(TPDList[[comm]][, "index"]), ]
    quantileTPD[[comm]] <- TPDList[[comm]]
  }
  names(quantileTPD) <- names(TPDList)

  trait1Edges <- unique(x$data$evaluation_grid[, 1])
  trait2Edges <- unique(x$data$evaluation_grid[, 2])
  imageMat <- array(
    NA,
    c(length(trait1Edges), length(trait2Edges), length(quantileTPD)),
    dimnames = list(trait1Edges, trait2Edges, names(TPDList))
  )
  for (comm in 1:length(TPDList)) {
    percentileSpace <- x$data$evaluation_grid
    percentileSpace$percentile <- quantileTPD[[comm]][, "percentile"]
    for (i in 1:length(trait2Edges)) {
      colAux <- subset(percentileSpace, percentileSpace[, 2] == trait2Edges[i])
      imageMat[, i, comm] <- 1 - colAux$percentile
    }
    imageMat[imageMat < (1 - thresholdPlot)] <- 0
  }
  return(imageMat)
}