R/geometry.R

In davidbolin/excursions: Excursion Sets and Contour Credibility Regions for Random Fields

## geometry.R
##
##   Copyright (C) 2014-2018 Finn Lindgren, David Bolin
##
##   This program is free software: you can redistribute it and/or modify
##   it under the terms of the GNU General Public License as published by
##   the Free Software Foundation, either version 3 of the License, or
##   (at your option) any later version.
##
##   This program is distributed in the hope that it will be useful,
##   but WITHOUT ANY WARRANTY; without even the implied warranty of
##   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##   GNU General Public License for more details.
##
##   You should have received a copy of the GNU General Public License
##   along with this program.  If not, see <http://www.gnu.org/licenses/>.



## Trace a length 2 contour segment through pixels
##
contour_segment_pixels <- function(c.x, c.y,
                                   edge.x, edge.y,
                                   pairs = FALSE) {
  nx <- length(edge.x)
  ny <- length(edge.y)
  d.x <- diff(c.x)
  d.y <- diff(c.y)
  if (abs(d.x) < abs(d.y)) {
    ## Swap x and y
    idx <- contour_segment_pixels(c.y, c.x, edge.y, edge.x, TRUE)
    idx <- list(x = idx$y, y = idx$x)
  } else if (d.y < 0) {
    ## Flip ordering of segment
    idx <- contour_segment_pixels(c.x[2:1], c.y[2:1], edge.x, edge.y, TRUE)
    idx <- list(x = rev(idx$x), y = rev(idx$y))
  } else {
    idx <- list(x = numeric(0), y = numeric(0))
    ## |d.x| >= d.y >= 0
    dir.x <- sign(d.x) ## Can be 0 only if dir.y is 0
    dir.y <- sign(d.y) ## Must now be >= 0
    if (dir.y == 0) {
      start.y <- min(which((c.y[1] >= edge.y[-ny]) &
        (c.y[1] <= edge.y[-1])))
      if (dir.x == 0) {
        start.x <- min(which((c.x[1] >= edge.x[-nx]) &
          (c.x[1] <= edge.x[-1])))
        end.x <- start.x
      } else if (dir.x > 0) {
        start.x <- min(which((c.x[1] >= edge.x[-nx]) &
          (c.x[1] < edge.x[-1])))
        end.x <- min(which((c.x[2] > edge.x[-nx]) &
          (c.x[2] <= edge.x[-1])))
      } else { ## dir.x < 0
        start.x <- min(which((c.x[1] > edge.x[-nx]) &
          (c.x[1] <= edge.x[-1])))
        end.x <- min(which((c.x[2] >= edge.x[-nx]) &
          (c.x[2] < edge.x[-1])))
      }
      idx$x <- start.x:end.x
      idx$y <- rep(start.y, length(idx$x))
    } else { ## dir.y > 0
      start.y <- min(which((c.y[1] >= edge.y[-ny]) &
        (c.y[1] < edge.y[-1])))
      end.y <- min(which((c.y[2] > edge.y[-ny]) &
        (c.y[2] <= edge.y[-1])))
      if (dir.x > 0) {
        start.x <- min(which((c.x[1] >= edge.x[-nx]) &
          (c.x[1] < edge.x[-1])))
        end.x <- min(which((c.x[2] > edge.x[-nx]) &
          (c.x[2] <= edge.x[-1])))
      } else { ## dir.x < 0
        start.x <- min(which((c.x[1] > edge.x[-nx]) &
          (c.x[1] <= edge.x[-1])))
        end.x <- min(which((c.x[2] >= edge.x[-nx]) &
          (c.x[2] < edge.x[-1])))
      }
      if (start.y == end.y) {
        idx$x <- start.x:end.x
        idx$y <- rep(start.y, length(idx$x))
      } else {
        ## Need to step through each y-pixel-row.
        curr.y <- start.y
        curr.x <- start.x
        while (curr.y < end.y) {
          next.y <- curr.y + 1
          ## Find intersection with edge.y[next.y]
          ## (x-c.x[1])*d.y/d.x + c.y[1] = edge.y[next.y]
          intersect.x <- c.x[1] + d.x / d.y * (edge.y[next.y] - c.y[1])
          if (dir.x > 0) {
            next.x <- min(which((intersect.x >= edge.x[-nx]) &
              (intersect.x < edge.x[-1])))
          } else { ## dir.x < 0
            next.x <- min(which((intersect.x > edge.x[-nx]) &
              (intersect.x <= edge.x[-1])))
          }
          idx$x <- c(idx$x, curr.x:next.x)
          idx$y <- c(idx$y, rep(curr.y, length(curr.x:next.x)))
          curr.y <- next.y
          curr.x <- next.x
        }
        idx$x <- c(idx$x, curr.x:end.x)
        idx$y <- c(idx$y, rep(curr.y, length(curr.x:end.x)))
      }
    }
  }
  if (pairs) {
    return(idx)
  } else {
    return((idx$y - 1) * (nx - 1) + idx$x)
  }
}


contour_pixels <- function(contourlines, pixelgrid, do.plot = 0) {
  if (inherits(pixelgrid, "pgrid")) {
    mid.x <- pixelgrid$upx
    mid.y <- pixelgrid$upy
    edge.x <- pixelgrid$ubx
    edge.y <- pixelgrid$uby
  } else if (inherits(pixelgrid, c("fm_lattice_2d", "inla.mesh.lattice"))) {
    mid.x <- pixelgrid$x
    mid.y <- pixelgrid$y
    step.x <- mid.x[2] - mid.x[1]
    step.y <- mid.y[2] - mid.y[1]
    edge.x <- c(mid.x[1] - step.x / 2, mid.x + step.x / 2)
    edge.y <- c(mid.y[1] - step.y / 2, mid.y + step.y / 2)
  } else {
    stop("Unsupported grid specification class.")
  }
  nx <- length(mid.x)
  ny <- length(mid.y)
  which.pixel <- vector("list", length(contourlines))
  for (level in seq_along(contourlines)) {
    c.x <- contourlines[[level]]$x
    c.y <- contourlines[[level]]$y
    for (segment in seq_len(length(c.x) - 1)) {
      tmp <- contour_segment_pixels(
        c.x[segment + (0:1)],
        c.y[segment + (0:1)],
        edge.x,
        edge.y
      )
      which.pixel[[level]] <- c(which.pixel[[level]], tmp)
      if (do.plot > 1) {
        setvec <- rep(0, nx * ny)
        setvec[sort(unique(unlist(which.pixel)))] <- 1
        image(mid.x, mid.y, matrix(setvec, nx, ny), col = c(0, 3))
        setvec <- rep(0, nx * ny)
        setvec[tmp] <- 1
        image(mid.x, mid.y, matrix(setvec, nx, ny),
          col = c(0, 4),
          add = TRUE
        )
        # plot.contourLines(contourlines, add=TRUE)
        lines(c.x[segment + (0:1)], c.y[segment + (0:1)], col = 2)
        if (do.plot > 2) {
          readline("next")
        }
      }
    }
  }
  if (do.plot > 0) {
    setvec <- rep(0, nx * ny)
    setvec[sort(unique(unlist(which.pixel)))] <- 1
    image(mid.x, mid.y, matrix(setvec, nx, ny), col = c(0, 3))
    # plot.contourLines(contourlines, add=TRUE)
    lines(c.x[segment + (0:1)], c.y[segment + (0:1)], col = 2)
  }

  sort(unique(unlist(which.pixel)))
}


## Connect segments into sequences looping around regions,
## in counterclockwise order.
## Initial segments have inner on their left hand side.
## grp.ccw keeps its orientation (multiple groups allowed)
## grp.cw is traversed in reverse orientation (multiple groups allowed)
##
## while (any segments remain) do
##   while (forward connected segments are found) do
##     follow segments forward
##   if (sequence not closed loop)
##     while (backward connected segments are found) do
##       follow segments backward
##
## sequences : list
## sequences[[k]] : node index vector for a single sequence
## seg : list
## seg[[k]] : segment index vector for a single sequence
## grp : list
## grp[[k]] : group index vector for a single sequence
connect.segments <- function(segment.set,
                             segment.grp = rep(0L, nrow(segment.set)),
                             grp.ccw = unique(segment.grp),
                             grp.cw = integer(0),
                             ccw = TRUE,
                             ambiguous.warning = FALSE) {
  ## Remove unneeded segments
  segment.idx <- seq_len(nrow(segment.set))
  segment.idx <- c(
    segment.idx[segment.grp %in% grp.ccw],
    segment.idx[segment.grp %in% grp.cw]
  )
  segment.set <- segment.set[segment.idx, , drop = FALSE]
  segment.grp <- as.integer(segment.grp[segment.idx])
  ## Unneeded segment removal done
  ## Reverse the direction of segments in grp.cw
  segment.set[segment.grp %in% grp.cw, ] <-
    segment.set[segment.grp %in% grp.cw, 2:1, drop = FALSE]
  ## Segment reversion done
  nE <- nrow(segment.set)
  if (nE == 0) {
    return(list(sequences = list(), seg = list(), grp = list()))
  }
  ## Remap nodes into 1...nV
  segments <- sort(unique(as.vector(segment.set)))
  nV <- length(segments)
  segments.reo <- Matrix::sparseMatrix(
    i = segments,
    j = rep(1, nV),
    x = seq_len(nV),
    dims = c(max(segments), 1)
  )
  segment.set <- matrix(
    segments.reo[as.vector(segment.set)],
    nrow(segment.set), ncol(segment.set)
  )
  ## Node remapping done

  segment.unused <- rep(TRUE, nE)
  segment.unused.idx <- which(segment.unused)
  segment.VV <- Matrix::sparseMatrix(
    i = segment.set[, 1],
    j = segment.set[, 2],
    x = seq_len(nE),
    dims = c(nV, nV)
  )
  loops.seg <- list()
  loops <- list()
  grp <- list()
  while (length(segment.unused.idx) > 0) {
    n <- 0
    loop.seg <- integer(0)
    ## Forward loop
    ##        message("Forwards")
    while (length(segment.unused.idx) > 0) {
      if (n == 0) {
        si <- segment.unused.idx[1]
      } else {
        si <- as.vector(as.matrix(segment.VV[segment.set[si, 2], ]))
        si <- si[si %in% segment.unused.idx]
        if (length(si) == 0) {
          ## End of sequence
          break
        } else {
          if ((length(si) > 1) && ambiguous.warning) {
            warning("Ambiguous segment sequence.")
          }
          si <- si[1]
        }
      }
      segment.unused.idx <-
        segment.unused.idx[segment.unused.idx != si]
      segment.unused[si] <- FALSE
      loop.seg <- c(loop.seg, si)
      n <- n + 1

      ##            print(loop.seg)
      ##            loop <- c(segment.set[loop.seg,1],
      ##                      segment.set[loop.seg[n],2])
      ##            print(loop)
    }
    if ((segment.set[loop.seg[n], 2] != segment.set[loop.seg[1], 1]) &&
      (length(segment.unused.idx) > 0)) {
      ## No closed sequence found
      ## Backward loop
      ##            message("Backwards")
      si <- loop.seg[1]
      while (length(segment.unused.idx) > 0) {
        si <- as.vector(as.matrix(segment.VV[, segment.set[si, 1]]))
        si <- si[si %in% segment.unused.idx]
        if (length(si) == 0) {
          ## End of sequence
          break
        } else if (length(si) > 1) {
          warning("Ambiguous segment sequence.")
          si <- si[1]
        } else {
          si <- si[1]
        }
        segment.unused.idx <-
          segment.unused.idx[segment.unused.idx != si]
        segment.unused[si] <- FALSE
        loop.seg <- c(si, loop.seg)
        n <- n + 1

        ##                print(loop.seg)
        ##                loop <- c(segment.set[loop.seg,1],
        ##                          segment.set[loop.seg[n],2])
        ##                print(loop)
      }
    }

    loop <- c(
      segment.set[loop.seg, 1],
      segment.set[loop.seg[n], 2]
    )
    loop.grp <- segment.grp[loop.seg]

    loops.seg <- c(loops.seg, list(loop.seg))
    loops <- c(loops, list(loop))
    grp <- c(grp, list(loop.grp))
  }

  ## Remap nodes and segments back to original indices
  for (k in seq_along(loops)) {
    loops[[k]] <- segments[loops[[k]]]
    loops.seg[[k]] <- segment.idx[loops.seg[[k]]]
  }
  ## Node and segment index remapping done

  if (!ccw) {
    for (k in seq_along(loops)) {
      loops[[k]] <- rev(loops[[k]])
      loops.seg[[k]] <- rev(loops.seg[[k]])
      grp[[k]] <- rev(grp[[k]])
    }
  }

  return(list(sequences = loops, seg = loops.seg, grp = grp))
}

## Compute simple outline of 1/0 set on a grid, eliminating spikes.
outline.on.grid <- function(z, grid) {
  if (missing(grid) || is.null(grid)) {
    if (!is.matrix(z)) {
      stop("grid not supplied, and z is not a matrix")
    }
    ni <- nrow(z)
    nj <- ncol(z)
    z <- (z != FALSE)
    grid <- list(x = seq(0, 1, length = ni), y = seq(0, 1, length = nj))
    grid$loc <- cbind(rep(grid$x, times = nj), rep(grid$y, each = ni))
  } else {
    ni <- grid$dims[1]
    nj <- grid$dims[2]
    z <- matrix(z != FALSE, ni, nj)
  }

  ij2k <- function(i, j) {
    return((j - 1) * ni + i)
  }

  seg <- matrix(integer(), 0, 2)
  bnd.seg <- matrix(integer(), 0, 2)

  ## Extract horizontal segment locations:
  zz.p <- z[-ni, -c(1, 2), drop = FALSE] + z[-1, -c(1, 2), drop = FALSE]
  zz.n <- z[-ni, -c(1, nj), drop = FALSE] + z[-1, -c(1, nj), drop = FALSE]
  zz.m <- z[-ni, -c(nj - 1, nj), drop = FALSE] + z[-1, -c(nj - 1, nj), drop = FALSE]
  zz <- (zz.n == 2) * (zz.p + zz.m > 0) * ((zz.p == 0) * 1 + (zz.m == 0) * 2)
  ## zz=0 : No segment, zz=1 : set is below, zz=2 : set is above
  ijv <- private.sparse.gettriplet(zz == 1)
  idx <- which(ijv$x > 0)
  seg <- rbind(
    seg,
    cbind(
      ij2k(ijv$i[idx] + 1, ijv$j[idx] + 1),
      ij2k(ijv$i[idx], ijv$j[idx] + 1)
    )
  )
  ijv <- private.sparse.gettriplet(zz == 2)
  idx <- which(ijv$x > 0)
  seg <- rbind(
    seg,
    cbind(
      ij2k(ijv$i[idx], ijv$j[idx] + 1),
      ij2k(ijv$i[idx] + 1, ijv$j[idx] + 1)
    )
  )

  ## Extract vertical segment locations:
  zz.p <- z[-c(1, 2), -nj, drop = FALSE] + z[-c(1, 2), -1, drop = FALSE]
  zz.n <- z[-c(1, ni), -nj, drop = FALSE] + z[-c(1, ni), -1, drop = FALSE]
  zz.m <- z[-c(ni - 1, ni), -nj, drop = FALSE] + z[-c(ni - 1, ni), -1, drop = FALSE]
  zz <- (zz.n == 2) * (zz.p + zz.m > 0) * ((zz.p == 0) * 1 + (zz.m == 0) * 2)
  ## zz=0 : No segment, zz=1 : set is on left, zz=2 : set is on right
  ijv <- private.sparse.gettriplet(zz == 1)
  idx <- which(ijv$x > 0)
  seg <- rbind(
    seg,
    cbind(
      ij2k(ijv$i[idx] + 1, ijv$j[idx]),
      ij2k(ijv$i[idx] + 1, ijv$j[idx] + 1)
    )
  )
  ijv <- private.sparse.gettriplet(zz == 2)
  idx <- which(ijv$x > 0)
  seg <- rbind(
    seg,
    cbind(
      ij2k(ijv$i[idx] + 1, ijv$j[idx] + 1),
      ij2k(ijv$i[idx] + 1, ijv$j[idx])
    )
  )

  ## Extract diagonal segment locations;
  ## Quadruples with three 1, one 0
  zz <- (z[-ni, -nj, drop = FALSE] + z[-1, -nj, drop = FALSE] +
    z[-ni, -1, drop = FALSE] + z[-1, -1, drop = FALSE])
  ## Which element was 0?
  zz <- (zz == 3) *
    (15 - (z[-ni, -nj, drop = FALSE] * 1 + z[-1, -nj, drop = FALSE] * 2 +
      z[-ni, -1, drop = FALSE] * 4 + z[-1, -1, drop = FALSE] * 8))
  ## zz=0 : No diagonal
  ## zz=1 : (0,0), zz=2 : (1,0), zz=4 : (0,1), zz=8 : (1,1)
  ijv <- private.sparse.gettriplet(zz == 1)
  idx <- which(ijv$x > 0)
  seg <- rbind(
    seg,
    cbind(
      ij2k(ijv$i[idx], ijv$j[idx] + 1),
      ij2k(ijv$i[idx] + 1, ijv$j[idx])
    )
  )
  ijv <- private.sparse.gettriplet(zz == 2)
  idx <- which(ijv$x > 0)
  seg <- rbind(
    seg,
    cbind(
      ij2k(ijv$i[idx], ijv$j[idx]),
      ij2k(ijv$i[idx] + 1, ijv$j[idx] + 1)
    )
  )
  ijv <- private.sparse.gettriplet(zz == 4)
  idx <- which(ijv$x > 0)
  seg <- rbind(
    seg,
    cbind(
      ij2k(ijv$i[idx] + 1, ijv$j[idx] + 1),
      ij2k(ijv$i[idx], ijv$j[idx])
    )
  )
  ijv <- private.sparse.gettriplet(zz == 8)
  idx <- which(ijv$x > 0)
  seg <- rbind(
    seg,
    cbind(
      ij2k(ijv$i[idx] + 1, ijv$j[idx]),
      ij2k(ijv$i[idx], ijv$j[idx] + 1)
    )
  )

  ## Extract horizontal boundary segment locations:
  zz.pm <- z[-ni, c(2, nj - 1), drop = FALSE] + z[-1, c(2, nj - 1), drop = FALSE]
  zz.n <- z[-ni, c(1, nj), drop = FALSE] + z[-1, c(1, nj), drop = FALSE]
  zz <- (zz.n == 2) * (zz.pm > 0) * matrix(rep(c(2, 1), each = ni - 1), ni - 1, 2)
  ## zz=0 : No segment, zz=1 : set is below, zz=2 : set is above
  ijv <- private.sparse.gettriplet(zz == 1)
  idx <- which(ijv$x > 0)
  bnd.seg <- rbind(
    bnd.seg,
    cbind(
      ij2k(ijv$i[idx] + 1, nj),
      ij2k(ijv$i[idx], nj)
    )
  )
  ijv <- private.sparse.gettriplet(zz == 2)
  idx <- which(ijv$x > 0)
  bnd.seg <- rbind(
    bnd.seg,
    cbind(
      ij2k(ijv$i[idx], 1),
      ij2k(ijv$i[idx] + 1, 1)
    )
  )

  ## Extract vertical boundary segment locations:
  zz.pm <- z[c(2, ni - 1), -nj, drop = FALSE] + z[c(2, ni - 1), -1, drop = FALSE]
  zz.n <- z[c(1, ni), -nj, drop = FALSE] + z[c(1, ni), -1, drop = FALSE]
  zz <- (zz.n == 2) * (zz.pm > 0) * matrix(rep(c(2, 1), times = nj - 1), 2, nj - 1)
  ## zz=0 : No segment, zz=1 : set is on left, zz=2 : set is on right
  ijv <- private.sparse.gettriplet(zz == 1)
  idx <- which(ijv$x > 0)
  bnd.seg <- rbind(
    bnd.seg,
    cbind(
      ij2k(ni, ijv$j[idx]),
      ij2k(ni, ijv$j[idx] + 1)
    )
  )
  ijv <- private.sparse.gettriplet(zz == 2)
  idx <- which(ijv$x > 0)
  bnd.seg <- rbind(
    bnd.seg,
    cbind(
      ij2k(1, ijv$j[idx] + 1),
      ij2k(1, ijv$j[idx])
    )
  )

  segment.grp <- rep(c(1L, 0L), c(nrow(seg), nrow(bnd.seg)))
  segment.set <- rbind(seg, bnd.seg)

  list(loc = grid$loc, idx = segment.set, grp = segment.grp)
}

## Compute simple outline of 1/0 set on a mesh, eliminating spikes.
## Returns a vector of triangle indices for the indicator set.
outlinetri.on.mesh <- function(z, mesh, complement = FALSE) {
  t.count <- rowSums(matrix((z >= 0.5)[mesh$graph$tv], nrow(mesh$graph$tv), 3))

  if (complement) {
    t.keep <- which(t.count < 3)
  } else {
    t.keep <- which(t.count == 3)
  }

  t.keep
}

## Compute simple outline of 1/0 set on a mesh, eliminating spikes.
outline.on.mesh <- function(z, mesh, complement = FALSE) {
  t.keep <- outlinetri.on.mesh(z, mesh, complement)

  if (length(t.keep) > 0) {
    graph <-
      generate_trigraph_properties(
        list(tv = mesh$graph$tv[t.keep, , drop = FALSE]),
        Nv = nrow(mesh$loc)
      )
    idx <- graph$ev[is.na(graph$ee), , drop = FALSE]
  } else {
    idx <- matrix(0L, 0, 2)
  }

  if (complement) {
    grp <- rep(1L, nrow(idx))
  } else {
    grp <- rep(3L, nrow(idx))
  }

  list(loc = mesh$loc, idx = idx, grp = grp)
}


#' Extract a part of a grid
#'
#' Extracts a part of a grid.
#'
#' @param z A matrix with values indicating which nodes that should be
#' present in the submesh.
#' @param grid A list with locations and dimensions of the grid.
#'
#' @return An `fm_mesh_2d` object.
#' @export
#' @author Finn Lindgren \email{finn.lindgren@@gmail.com}
#' @examples
#' \dontrun{
#' if (require("fmesher")) {
#'   nxy <- 40
#'   x <- seq(from = 0, to = 4, length.out = nxy)
#'   lattice <- fm_lattice_2d(x = x, y = x)
#'   mesh <- fm_rcdt_2d_inla(lattice = lattice, extend = FALSE, refine = FALSE)
#'
#'   # extract a part of the mesh inside a circle
#'   xy.in <- rowSums((mesh$loc[, 1:2] - 2)^2) < 1
#'   submesh <- submesh.grid(
#'     matrix(xy.in, nxy, nxy),
#'     list(loc = mesh$loc, dim = c(nxy, nxy))
#'   )
#'   plot(mesh$loc[, 1:2])
#'   lines(2 + cos(seq(0, 2 * pi, length.out = 100)), 2 + sin(seq(0, 2 * pi, length.out = 100)))
#'   plot(submesh, add = TRUE)
#'   points(mesh$loc[xy.in, 1:2], col = "2")
#' }
#' }
#'
submesh.grid <- function(z, grid = NULL) {
  outline <- outline.on.grid(z, grid)
  fmesher::fm_rcdt_2d_inla(
    loc = grid$loc,
    boundary = as.fm_segm.outline(outline),
    refine = FALSE
  )
}


#' Extract a part of a mesh
#'
#' Extracts a part of a mesh
#'
#' @param z A matrix with values indicating which nodes that should be
#' present in the submesh.
#' @param mesh An `fm_mesh_2d` object.
#'
#' @return An `fm_mesh_2d` object.
#' @export
#' @author Finn Lindgren \email{finn.lindgren@@gmail.com}
#' @examples
#' \dontrun{
#' if (require(fmesher)) {
#'   nxy <- 30
#'   x <- seq(from = 0, to = 4, length.out = nxy)
#'   lattice <- fm_lattice_2d(x = x, y = x)
#'   mesh <- fm_mesh_2d_inla(lattice = lattice, extend = FALSE, refine = FALSE)
#'
#'   # extract a part of the mesh inside a circle
#'   xy.in <- rowSums((mesh$loc[, 1:2] - 2)^2) < 1
#'   submesh <- excursions:::submesh.mesh(matrix(xy.in, nxy, nxy), mesh)
#'   plot(mesh$loc[, 1:2])
#'   lines(2 + cos(seq(0, 2 * pi, length.out = 100)), 2 + sin(seq(0, 2 * pi, length.out = 100)))
#'   plot(submesh, add = TRUE)
#'   points(mesh$loc[xy.in, 1:2], col = "2")
#' }
#' }
#'
submesh.mesh <- function(z, mesh) {
  submesh.mesh.tri(outlinetri.on.mesh(z, mesh), mesh)
}
submesh.mesh.tri <- function(tri, mesh) {
  tv <- mesh$graph$tv[tri, , drop = FALSE]
  v <- sort(unique(as.vector(tv)))
  idx <- rep(as.integer(NA), nrow(mesh$loc))
  idx[v] <- seq_len(length(v))
  tv <- matrix(idx[tv], nrow(tv), 3)
  loc <- mesh$loc[v, , drop = FALSE]

  mesh <- fmesher::fm_rcdt_2d_inla(loc = loc, tv = tv, refine = FALSE)

  idx <- rep(as.integer(NA), length(idx))
  idx[v] <- mesh$idx$loc
  mesh$idx$loc <- idx

  mesh
}


as.sp.outline <- function(outline,
                          grp.ccw = unique(outline$grp),
                          grp.cw = integer(0),
                          ccw = FALSE,
                          ambiguous.warning = FALSE,
                          ID = "outline",
                          closed = TRUE,
                          ...) {
  seg <- connect.segments(outline$idx, outline$grp,
    grp.ccw = grp.ccw,
    grp.cw = grp.cw,
    ccw = ccw,
    ambiguous.warning = ambiguous.warning
  )

  coords <- list()
  for (k in seq_along(seg$sequences)) {
    coords <- c(coords, list(outline$loc[seg$sequences[[k]],
      1:2,
      drop = FALSE
    ]))
  }

  if (closed) {
    as.Polygons.raw(coords, ID = ID)
  } else {
    as.Lines.raw(coords, ID = ID)
  }
}


as.fm_segm.outline <- function(outline,
                               grp.ccw = unique(outline$grp),
                               grp.cw = integer(0),
                               grp,
                               ...) {
  ik.ccw <- outline$grp %in% grp.ccw
  ik.cw <- outline$grp %in% grp.cw
  if (missing(grp)) {
    grp <- c(outline$grp[ik.ccw], outline$grp[ik.cw])
  }
  fmesher::fm_segm(
    loc = outline$loc,
    idx = rbind(
      outline$idx[ik.ccw, , drop = FALSE],
      outline$idx[ik.cw, 2:1, drop = FALSE]
    ),
    grp = grp
  )
}



as.Polygons.raw <- function(sequences, ID = " ") {
  if (!requireNamespace("sp", quietly = TRUE)) {
    stop("The 'sp' package is needed.")
  }
  polys <- lapply(
    sequences,
    function(x) {
      if (is.list(x)) {
        p <- sp::Polygon(cbind(x$x, x$y))
      } else {
        p <- sp::Polygon(x)
      }
      if (p@area > 0) {
        p
      } else {
        warning("Skipping zero area polygon, in as.Polygons.raw.")
        NULL
      }
    }
  )
  if (length(polys) == 0) {
    sp <- NULL
  } else {
    ok <- unlist(lapply(polys, function(x) is(x, "Polygon")))
    sp <- sp::Polygons(polys[ok], ID = ID)
  }
  sp
}

as.Lines.raw <- function(cl, ID = " ") {
  if (!requireNamespace("sp", quietly = TRUE)) {
    stop("The 'sp' package is needed.")
  }
  polys <- lapply(
    cl,
    function(x) {
      if (is.list(x)) {
        p <- sp::Line(cbind(x$x, x$y))
      } else {
        p <- sp::Line(x)
      }
      p
    }
  )
  if (length(polys) == 0) {
    sp <- NULL
  } else {
    sp <- sp::Lines(polys, ID = ID)
  }
  sp
}

#' Calculate contour curves on a triangulation
#'
#' Calculates contour curves and/or regions between them,
#' for functions defined on a triangulation
#'
#' @export
#' @param x An object generated by a call to `fm_mesh_2d` or
#'   `fm_rcdt_2d`, a triangle-vertex index matrix, or a list
#'   of triangulation information, `list(loc, graph=list(tv))`.
#' @param z A vector containing the values to be contoured
#'   (`NA`s are allowed).
#' @param nlevels Number of contour levels desired, if and only if
#'   `levels` is not supplied.
#' @param levels Numeric vector of levels at which to calculate contour lines.
#' @param ... Additional arguments passed to the other methods.
#' @return For `tricontour`, a list with some of the fields that
#'   `fm_segm` objects have:
#'   \item{loc}{A coordinate matrix}
#'   \item{idx}{Contour segment indices, as a 2-column matrix, each row
#'     indexing a single segment}
#'   \item{grp}{A vector of group labels.  Each segment has a label, in
#'     `1,...,nlevels*2+1`, where even labels indicate interior
#'     on-level contour segments, and odd labels indicate boundary segments
#'     between levels.}
#'   For `tricontourmap`, a list:
#'   \item{contour}{A list of `sp` or `fm_segm` objects
#'     defining countour curves (level sets)}
#'   \item{map}{A list of `sp` or `fm_segm` objects
#'     enclosing regions between level sets}
#'
#' @author Finn Lindgren \email{finn.lindgren@@gmail.com}
#'
#' @examples
#' \dontrun{
#' if (require("fmesher") &&
#'   require("sp")) {
#'   ## Generate mesh and SPDE model
#'   n.lattice <- 20 # increase for more interesting, but slower, examples
#'   x <- seq(from = 0, to = 10, length.out = n.lattice)
#'   lattice <- fm_lattice_2d(x = x, y = x)
#'   mesh <- fm_rcdt_2d_inla(lattice = lattice, extend = FALSE, refine = FALSE)
#'
#'   ## Generate an artificial sample
#'   sigma2.e <- 0.1
#'   n.obs <- 1000
#'   obs.loc <- cbind(
#'     runif(n.obs) * diff(range(x)) + min(x),
#'     runif(n.obs) * diff(range(x)) + min(x)
#'   )
#'   Q <- fm_matern_precision(mesh, alpha = 2, rho = 3, sigma = 1)
#'   x <- fm_sample(n = 1, Q = Q)
#'   A <- fm_basis(mesh, loc = obs.loc)
#'   Y <- as.vector(A %*% x + rnorm(n.obs) * sqrt(sigma2.e))
#'
#'   ## Calculate posterior
#'   Q.post <- (Q + (t(A) %*% A) / sigma2.e)
#'   mu.post <- as.vector(solve(Q.post, (t(A) %*% Y) / sigma2.e))
#'
#'   ## Calculate continuous contours
#'   tric <- tricontour(mesh,
#'     z = mu.post,
#'     levels = as.vector(quantile(x, c(0.25, 0.75)))
#'   )
#'
#'   ## Discrete domain contours
#'   map <- contourmap(
#'     n.levels = 2, mu = mu.post, Q = Q.post,
#'     alpha = 0.1, compute = list(F = FALSE), max.threads = 1
#'   )
#'
#'   ## Calculate continuous contour map
#'   setsc <- tricontourmap(mesh,
#'     z = mu.post,
#'     levels = as.vector(quantile(x, c(0.25, 0.75)))
#'   )
#'
#'   ## Plot the results
#'   reo <- mesh$idx$lattice
#'   idx.setsc <- setdiff(names(setsc$map), "-1")
#'   cols2 <- contourmap.colors(map,
#'     col = heat.colors(100, 0.5, rev = TRUE),
#'     credible.col = grey(0.5, 0)
#'   )
#'   names(cols2) <- as.character(-1:2)
#'
#'   par(mfrow = c(1, 2))
#'   image(matrix(mu.post[reo], n.lattice, n.lattice),
#'     main = "mean", axes = FALSE, asp = 1
#'   )
#'   plot(setsc$map[idx.setsc], col = cols2[idx.setsc])
#'   par(mfrow = c(1, 1))
#' }
#' }
tricontour <-
  function(x, z, nlevels = 10,
           levels = pretty(range(z, na.rm = TRUE), nlevels),
           ...) {
    UseMethod("tricontour")
  }

#' @rdname tricontour
#' @export
tricontour.fm_mesh_2d <-
  function(x, z, nlevels = 10,
           levels = pretty(range(z, na.rm = TRUE), nlevels),
           ...) {
    tricontour.list(x$graph,
      z = z,
      nlevels = nlevels, levels = levels,
      loc = x$loc, ...
    )
  }

#' @rdname tricontour
#' @export
#' @param loc coordinate matrix, to be supplied when `x` is given as a
#'   triangle-vertex index matrix only.
tricontour.matrix <-
  function(x, z, nlevels = 10,
           levels = pretty(range(z, na.rm = TRUE), nlevels),
           loc, ...) {
    tricontour.list(list(tv = x),
      z = z,
      nlevels = nlevels, levels = levels,
      loc = loc, ...
    )
  }



## Generate triangulation graph properties
## Nt,Ne,Nv,ev,et,eti,ee,te,tt,tti
generate_trigraph_properties <- function(x, Nv = NULL) {
  stopifnot(is.list(x))
  stopifnot("tv" %in% names(x))

  x$Nt <- nrow(x$tv)
  x$Ne <- 3 * x$Nt ## number of unidirectional edges
  if (is.null(Nv)) {
    x$Nv <- max(as.vector(x$tv))
  } else {
    x$Nv <- Nv
  }

  x$ev <- cbind(
    as.vector(x$tv[, c(2, 3, 1)]),
    as.vector(x$tv[, c(3, 1, 2)])
  )
  x$et <- rep(seq_len(x$Nt), times = 3)
  x$eti <- rep(1:3, each = x$Nt) ## Opposing vertex within-triangle-indices
  x$te <- matrix(seq_len(x$Ne), x$Nt, 3)
  ev <- Matrix::sparseMatrix(
    i = rep(seq_len(x$Ne), times = 2),
    j = as.vector(x$ev),
    x = rep(1, x$Ne * 2),
    dims = c(x$Ne, x$Nv)
  )
  ev.tr <- private.sparse.gettriplet(ev %*% t(ev))
  ee <- cbind(
    ev.tr$i[(ev.tr$x == 2) & (ev.tr$i != ev.tr$j)],
    ev.tr$j[(ev.tr$x == 2) & (ev.tr$i != ev.tr$j)]
  )
  x$ee <- rep(NA, x$Ne)
  x$ee[ee[, 1]] <- ee[, 2]

  if (is.null(x$tt) ||
    (nrow(x$tt) != x$Nt)) { ## Workaround for bug in fmesher < 2014-09-12
    x$tt <- matrix(x$et[x$ee], x$Nt, 3)
  }
  if (is.null(x$tti)) {
    x$tti <- matrix(x$eti[x$ee], x$Nt, 3)
  }

  x
}



display.dim.list <- function(x) {
  lapply(
    as.list(sort(names(x))),
    function(xx) {
      sz <- dim(x[[xx]])
      type <- mode(x[[xx]])
      cl <- class(x[[xx]])[[1]]
      if (is.null(sz)) {
        sz <- length(x[[xx]])
      }
      message(paste(xx, " = ", paste(sz, collapse = " x "),
        " (", type, ", ", cl, ")",
        sep = ""
      ))
    }
  )
  invisible()
}

#' @rdname tricontour
#' @export
#' @param type `"+"` or `"-"`, indicating positive or negative
#'  association. For `+`, the generated contours enclose regions
#'  where \eqn{u_1 \leq z < u_2}, for `-` the regions fulfil \eqn{u_1
#'    < z \leq u_2}.
#' @param tol tolerance for determining if the value at a vertex lies on a level.
tricontour.list <- function(x, z, nlevels = 10,
                            levels = pretty(range(z, na.rm = TRUE), nlevels),
                            loc, type = c("+", "-"), tol = 1e-7, ...) {
  ## Returns val=list(loc, idx, grp), where
  ##   grp = 1,...,nlevels*2+1, level groups are even, 2,4,...
  ## Suitable for
  ##   fm_segm(val$loc, val$idx[val$grp==k], val$idx[val$grp==k])
  ##     (supports R2 and S2)
  ## and, for odd k=1,3,...,nlevels*2-1,nlevels*2+1,
  ##   seg <- as.fm_segm.outline(val, grp.ccw=c(k-1,k), grp.cw=c(k+1))
  ##   sp <- as.sp.outline(val, grp.ccw=c(k-1,k), grp.cw=c(k+1), ccw=FALSE)

  type <- match.arg(type)
  nlevels <- length(levels)

  x <- generate_trigraph_properties(x)
  Nv <- x$Nv

  ## Find vertices on levels
  ## For each edge on a level, store edge if
  ##   it is a boundary edge, or
  ##   opposing vertices are +/-, and either
  ##     0/- (type="+", u1 <= z < u2) or
  ##     +/0 (type="-", u1 < z <= u2)
  ## For each edge crossing at least one level,
  ##   calculate splitting vertices
  ##   if boundary edge, store new split edges
  ## For each triangle, find non-level edge crossings, and
  ##   store new vertex-edge crossing edges
  ##   store new edge-edge crossing edges
  idx <- matrix(0, 0, 2)
  grp <- integer(0)

  ## Find vertices on levels
  vcross.lev <- integer(length(z))
  for (lev in seq_along(levels)) {
    signv <- (z > levels[lev] + tol) - (z < levels[lev] - tol)
    vcross.lev[signv == 0] <- lev
  }
  ## Find level crossing span for each edge (includes flat edges in levels)
  ecross.grp.lower <- rep(1L, x$Ne)
  ecross.grp.upper <- rep(2L * length(levels) + 1L, x$Ne)
  for (lev in seq_along(levels)) {
    signv <- (z > levels[lev] + tol) - (z < levels[lev] - tol)
    lev.grp <- 2L * lev
    i <- pmin(signv[x$ev[, 1]], signv[x$ev[, 2]])
    ecross.grp.lower[i == 0] <- lev.grp
    ecross.grp.lower[i > 0] <- lev.grp + 1L
  }
  for (lev in rev(seq_along(levels))) {
    signv <- (z > levels[lev] + tol) - (z < levels[lev] - tol)
    lev.grp <- 2L * lev
    i <- pmax(signv[x$ev[, 1]], signv[x$ev[, 2]])
    ecross.grp.upper[i == 0] <- lev.grp
    ecross.grp.upper[i < 0] <- lev.grp - 1L
  }

  ## For each edge on a level, store edge if ...
  ##   it is a boundary edge, or
  ##   opposing vertices are +/-, and either
  ##     0/- (type="+", u1 <= z < u2) or
  ##     +/0 (type="-", u1 < z <= u2)
  cross1 <- vcross.lev[x$ev[, 1]] ## left neighbour
  cross2 <- vcross.lev[x$ev[, 2]] ## right neighbour
  vv.edges <- which((cross1 > 0) & (cross2 > 0) & (cross1 == cross2))
  for (edge in vv.edges) {
    lev <- cross1[edge]
    v1 <- x$tv[x$et[edge], x$eti[edge]]
    sign1 <- (z[v1] > levels[lev] + tol) - (z[v1] < levels[lev] - tol)
    neighb.t <- x$tt[x$et[edge] + (x$eti[edge] - 1) * x$Nt]
    if (is.na(neighb.t)) {
      v2 <- NA
      sign2 <- NA
    } else {
      v2 <- x$tv[neighb.t, x$tti[x$et[edge] + (x$eti[edge] - 1) * x$Nt]]
      sign2 <- (z[v2] > levels[lev] + tol) - (z[v2] < levels[lev] - tol)
    }
    if (is.na(neighb.t)) {
      ## Make sure the edge gets the right label
      if (sign1 != 0) {
        idx <- rbind(idx, x$ev[edge, ])
        ## Associate with the neighbour
        grp <- c(grp, lev * 2L + sign1)
      } else {
        idx <- rbind(idx, x$ev[edge, ])
        grp <- c(grp, lev * 2L + (type == "+") - (type == "-"))
      }
    } else if (((sign1 > 0) && (sign2 < 0)) ||
      ((type == "+") && ((sign1 == 0) && (sign2 < 0))) ||
      ((type == "-") && ((sign1 > 0) && (sign2 == 0)))) {
      idx <- rbind(idx, x$ev[edge, ])
      grp <- c(grp, lev * 2L)
    }
  }

  ## For each boundary edge entirely between levels
  ##   store edge
  e.lower <- ecross.grp.lower + ((ecross.grp.lower - 1) %% 2)
  e.upper <- ecross.grp.upper - ((ecross.grp.upper - 1) %% 2)
  e.on.bnd <- which(is.na(x$tti[x$et + (x$eti - 1) * x$Nt]))
  e.noncrossing <- e.on.bnd[e.lower[e.on.bnd] == e.upper[e.on.bnd]]
  idx <- rbind(idx, x$ev[e.noncrossing, , drop = FALSE])
  grp <- c(grp, as.integer(e.lower[e.noncrossing]))

  ## For each edge crossing at least one level,
  ##   calculate splitting vertices
  ##   if boundary edge, store new split edges
  e.crossing <- which(e.lower < e.upper)
  e.bndcrossing <- e.on.bnd[e.lower[e.on.bnd] < e.upper[e.on.bnd]]
  Nnewv <- (sum(e.upper[e.crossing] - e.lower[e.crossing]) / 2 +
    sum(e.upper[e.bndcrossing] - e.lower[e.bndcrossing]) / 2) / 2
  loc.new <- matrix(NA, Nnewv, ncol(loc))
  loc.last <- 0L
  e.newv <- sparseMatrix(
    i = integer(0), j = integer(0), x = double(0),
    dims = c(x$Ne, length(levels))
  )
  for (edge in e.crossing) {
    is.boundary.edge <- is.na(x$tt[x$et[edge], x$eti[edge]])
    if (is.boundary.edge) {
      edge.reverse <- NA
    } else {
      ## Interior edges appear twice; handle each only once.
      if (x$ev[edge, 1] > x$ev[edge, 2]) {
        next
      }
      edge.reverse <- x$te[
        x$tt[x$et[edge], x$eti[edge]],
        x$tti[x$et[edge], x$eti[edge]]
      ]
    }
    ## lev = (1-beta_lev) * z1 + beta_lev * z2
    ##     = z1 + beta_lev * (z2-z1)
    ## beta_lev = (lev-z1)/(z2-z1)
    e.levels <- ((e.lower[edge] + 1) / 2):((e.upper[edge] - 1) / 2)
    loc.new.idx <- loc.last + seq_along(e.levels)
    beta <- ((levels[e.levels] - z[x$ev[edge, 1]]) /
      (z[x$ev[edge, 2]] - z[x$ev[edge, 1]]))
    ## print(str(loc.new))
    ## print(e.levels)
    ## print(c(loc.last, loc.new.idx))
    ##    message("???")
    ##    message(paste(paste(dim(loc.new), collapse=", "),
    ##                  loc.new.idx, collapse="; "))
    ## Temporary workaround for boundary part error:
    if (max(loc.new.idx) > nrow(loc.new)) {
      ##      browser()
      loc.new <- rbind(
        loc.new,
        matrix(
          NA,
          max(loc.new.idx) - nrow(loc.new),
          ncol(loc.new)
        )
      )
    }

    ##    xxx <-
    ##      (as.matrix(1-beta) %*% loc[x$ev[edge,1],,drop=FALSE]+
    ##       as.matrix(beta) %*% loc[x$ev[edge,2],,drop=FALSE])

    loc.new[loc.new.idx, ] <-
      (as.matrix(1 - beta) %*% loc[x$ev[edge, 1], , drop = FALSE] +
        as.matrix(beta) %*% loc[x$ev[edge, 2], , drop = FALSE])

    e.newv[edge, e.levels] <- Nv + loc.new.idx
    if (!is.boundary.edge) {
      e.newv[edge.reverse, e.levels] <- Nv + loc.new.idx
    } else { ## edge is a boundary edge, handle now
      ev <- x$ev[edge, ]
      the.levels <- which(e.newv[edge, ] > 0)
      the.loc.idx <- e.newv[edge, the.levels]
      the.levels <- c(min(the.levels) - 1L, the.levels) * 2L + 1L
      if (z[ev[1]] > z[ev[2]]) {
        the.loc.idx <- rev(the.loc.idx)
        the.levels <- rev(the.levels)
      }

      idx <- rbind(idx, cbind(
        c(ev[1], the.loc.idx),
        c(the.loc.idx, ev[2])
      ))
      grp <- c(grp, the.levels)
    }
    loc.last <- loc.last + length(e.levels)
  }
  loc <- rbind(loc, loc.new)
  Nv <- nrow(loc)

  ## For each triangle, find non-level edge crossings, and
  ##   store new vertex-edge crossing edges
  ##   store new edge-edge crossing edges
  tris <- unique(x$et[e.crossing])
  for (tri in tris) {
    ## connect vertex-edge
    v.lev <- vcross.lev[x$tv[tri, ]]
    for (vi in which(v.lev > 0)) {
      opposite.edge <- x$te[tri, vi]
      opposite.v <- e.newv[opposite.edge, v.lev[vi]]
      if (opposite.v > 0) {
        ## v2 on the right, v3 on the left
        v123 <- x$tv[tri, ((vi + (0:2) - 1) %% 3) + 1]
        if (z[v123[3]] > z[v123[1]]) {
          idx <- rbind(idx, cbind(v123[1], opposite.v))
        } else {
          idx <- rbind(idx, cbind(opposite.v, v123[1]))
        }
        grp <- c(grp, v.lev[vi] * 2L)
      }
    }
    ## connect edge-edge
    for (ei in 1:3) {
      edge <- x$te[tri, ei]
      next.edge <- x$te[tri, (ei %% 3L) + 1L]
      e.lev <- which(e.newv[edge, ] > 0)
      e.lev <- e.lev[e.newv[next.edge, e.lev] > 0]
      if (length(e.lev) > 0) {
        ## v1 on the left, v2 on the right
        v12 <- x$ev[edge, ]
        if (z[v12[1]] > z[v12[2]]) {
          idx <- rbind(idx, cbind(
            e.newv[edge, e.lev],
            e.newv[next.edge, e.lev]
          ))
        } else {
          idx <- rbind(idx, cbind(
            e.newv[next.edge, e.lev],
            e.newv[edge, e.lev]
          ))
        }
        grp <- c(grp, e.lev * 2L)
      }
    }
  }

  ## Filter out unused nodes
  reo <- sort(unique(as.vector(idx)))
  loc <- loc[reo, , drop = FALSE]
  ireo <- integer(Nv)
  ireo[reo] <- seq_along(reo)
  idx <- matrix(ireo[idx], nrow(idx), 2)

  list(loc = loc, idx = idx, grp = grp)
}






#' @rdname tricontour
#' @export
tricontourmap <- function(x, z, nlevels = 10,
                          levels = pretty(range(z, na.rm = TRUE), nlevels),
                          ...) {
  UseMethod("tricontourmap")
}

#' @rdname tricontour
#' @export
tricontourmap.fm_mesh_2d <-
  function(x, z, nlevels = 10,
           levels = pretty(range(z, na.rm = TRUE), nlevels),
           ...) {
    tricontourmap.list(x$graph,
      z = z,
      nlevels = nlevels, levels = levels,
      loc = x$loc, ...
    )
  }

#' @rdname tricontour
#' @export
tricontourmap.matrix <-
  function(x, z, nlevels = 10,
           levels = pretty(range(z, na.rm = TRUE), nlevels),
           loc, ...) {
    tricontourmap.list(list(tv = x),
      z = z,
      nlevels = nlevels, levels = levels,
      loc = loc, ...
    )
  }



#' @rdname tricontour
#' @param output The format of the generated output.  Implemented options
#'   are `"sp"` (default) and `"fm"` (and deprecated
#'   `"inla.mesh.segment"`, converted to `"fm"`).
#' @importFrom fmesher fm_segm
#' @export
tricontourmap.list <-
  function(x, z, nlevels = 10,
           levels = pretty(range(z, na.rm = TRUE), nlevels),
           loc, type = c("+", "-"), tol = 1e-7,
           output = c("sp", "fm", "inla.mesh.segment"), ...) {
    type <- match.arg(type)
    output <- match.arg(output)
    if (output == "inla.mesh.segment") {
      output <- "fm"
    }
    nlevels <- length(levels)

    if (output == "sp") {
      if (!requireNamespace("sp", quietly = TRUE)) {
        stop("The 'sp' package is needed.")
      }
    }

    tric <- tricontour(
      x = x, z = z, nlevels = nlevels, levels = levels,
      loc = loc, type = type, tol = tol, ...
    )

    out <- list(map = list(), contour = list())
    for (k in seq_len(nlevels + 1L) * 2L - 1L) {
      ID <- as.character((k - 1) / 2)
      if (output == "sp") {
        spobj <- tryCatch(
          as.sp.outline(tric,
            grp.ccw = c(k - 1, k),
            grp.cw = c(k + 1),
            ccw = FALSE,
            closed = TRUE,
            ID = ID
          ),
          error = function(e) NULL
        )
        if (!is.null(spobj)) {
          out$map[[ID]] <- spobj
        }
      } else {
        out$map[[ID]] <-
          as.fm_segm.outline(tric,
            grp.ccw = c(k - 1, k),
            grp.cw = c(k + 1),
            grp = (k - 1) / 2
          )
      }
    }
    if (output == "sp") {
      out$map <- sp::SpatialPolygons(out$map)
    } else {
      out$map <- do.call(fm_segm, out$map)
    }
    for (k in seq_len(nlevels)) {
      if (output == "sp") {
        ID <- as.character(k)
        spobj <- tryCatch(
          as.sp.outline(tric,
            grp.ccw = k * 2L,
            ccw = FALSE,
            closed = FALSE,
            ID = ID
          ),
          error = function(e) NULL
        )
        if (!is.null(spobj)) {
          out$contour[[ID]] <- spobj
        }
      } else {
        out$contour[[as.character(k)]] <-
          as.fm_segm.outline(tric,
            grp.ccw = k * 2L,
            grp = k
          )
      }
    }
    if (output == "sp") {
      if (length(out$contour) > 0) {
        out$contour <- sp::SpatialLines(out$contour)
      } else {
        out$contour <- NULL
      }
    } else {
      out$contour <- do.call(fm_segm, out$contour)
    }

    out
  }


tricontour_step <- function(x, z, levels, loc, ...) {
  stopifnot(length(levels) == 1)

  x <- list(loc = loc, graph = x)
  outline.lower <- outline.on.mesh(z >= levels, x, TRUE)
  outline.upper <- outline.on.mesh(z >= levels, x, FALSE)

  list(
    loc = loc,
    idx = rbind(outline.lower$idx, outline.upper$idx),
    grp = c(outline.lower$grp, outline.upper$grp)
  )
}





## excursions --> E,F
## contourmap --> E,P0123,F
## simconf
## continterp(excurobj, grid or mesh, outputgrid(opt), alpha, method)
## gaussint


## Input: One of
##   fm_mesh_2d
##   fm_lattice_2d
##   list(loc, dims, ...)
##   list(x, y, ...)
## The last 3 are all treated as topological lattices, and code in
## submesh.grid() assumes that the lattice boxes are convex.
## Output:
##   list(loc, dims, geometry, manifold)
get.geometry <- function(geometry) {
  geometrytype <- ""
  manifoldtype <- ""
  if (inherits(geometry, c("fm_mesh_2d", "inla.mesh"))) {
    loc <- geometry$loc
    dims <- nrow(loc)
    geometrytype <- "mesh"
    manifoldtype <- geometry$manifold
  } else if (inherits(geometry, c("fm_lattice_2d", "inla.mesh.lattice")) ||
    is.list(geometry)) {
    if (("loc" %in% names(geometry)) && ("dims" %in% names(geometry))) {
      loc <- geometry$loc
      dims <- geometry$dims
      geometrytype <- "lattice"
      manifoldtype <- "R2"
    } else if ("x" %in% names(geometry)) {
      geometrytype <- "lattice"
      if ("y" %in% names(geometry)) {
        loc <- as.matrix(expand.grid(geometry$x, geometry$y))
        dims <- c(length(geometry$x), length(geometry$y))
        manifoldtype <- "R2"
      } else {
        loc <- geometry$x
        dims <- length(loc)
        manifoldtype <- "R1"
      }
    }
  }
  geometrytype <- match.arg(geometrytype, c("mesh", "lattice"))
  manifoldtype <- match.arg(manifoldtype, c("M", "R1", "R2", "S2"))
  list(loc = loc, dims = dims, geometry = geometrytype, manifold = manifoldtype)
}


## Input:
##   list(loc, graph=list(tv, ...)) or an fm_mesh_2d
## Output:
##   list(loc, graph=list(tv), A, idx=list(loc))
## TODO: switch to using fm_subdivide()
subdivide.mesh <- function(mesh) {
  graph <- generate_trigraph_properties(mesh$graph, nrow(mesh$loc))

  v1 <- seq_len(graph$Nv)
  edges.boundary <- which(is.na(graph$ee))
  edges.interior <- which(!is.na(graph$ee))
  edges.interior.main <- edges.interior[graph$ev[edges.interior, 1] <
    graph$ev[edges.interior, 2]]
  edges.interior.secondary <- graph$ee[edges.interior.main]
  Neb <- length(edges.boundary)
  Nei <- length(edges.interior.main)
  Nv2 <- graph$Nv + Neb + Nei
  v2.boundary <- graph$Nv + seq_len(Neb)
  v2.interior <- graph$Nv + Neb + seq_len(Nei)
  edge.split.v <- integer(length(graph$ee))
  edge.split.v[edges.boundary] <- v2.boundary
  edge.split.v[edges.interior.main] <- v2.interior
  edge.split.v[edges.interior.secondary] <- v2.interior

  ## Mapping matrix from mesh nodes to sub-mesh nodes
  idx <- seq_len(graph$Nv) ## Same indices as in input mesh
  ridx <- seq_len(graph$Nv)
  ##  ridx[mesh$idx$loc[!is.na(mesh$idx$loc)]] <-
  ##    which(!is.na(mesh$idx$loc))
  A <- sparseMatrix(
    i = (c(
      v1,
      rep(v2.boundary, times = 2),
      rep(v2.interior, times = 2)
    )),
    j = (ridx[c(
      v1,
      as.vector(graph$ev[edges.boundary, ]),
      as.vector(graph$ev[edges.interior.main, ])
    )]),
    x = (c(
      rep(1, graph$Nv),
      rep(1 / 2, 2 * Neb),
      rep(1 / 2, 2 * Nei)
    )),
    dims = c(Nv2, length(ridx))
  )
  loc <- as.matrix(A[, ridx, drop = FALSE] %*% mesh$loc)
  tv <- rbind(
    cbind(
      edge.split.v[graph$te[, 1]],
      edge.split.v[graph$te[, 2]],
      edge.split.v[graph$te[, 3]]
    ),
    cbind(
      graph$tv[, 1],
      edge.split.v[graph$te[, 3]],
      edge.split.v[graph$te[, 2]]
    ),
    cbind(
      graph$tv[, 2],
      edge.split.v[graph$te[, 1]],
      edge.split.v[graph$te[, 3]]
    ),
    cbind(
      graph$tv[, 3],
      edge.split.v[graph$te[, 2]],
      edge.split.v[graph$te[, 1]]
    )
  )

  newmesh <- fmesher::fm_rcdt_2d_inla(loc = loc, tv = tv)
  ## Handle possible node reordering in fm_rcdt_2d_inla()
  newmesh$idx$loc <- newmesh$idx$loc[idx]
  ## Add mapping matrix
  newmesh$A <- A

  newmesh
}







## New geometry implementation ####

# Copy 'G' and interpolate 'F' within coherent single level regions.
F.interpolation <- function(F.geometry, F, G, type, method, subdivisions = 1) {
  # Construct interpolation mesh
  F.geometry.A <- list()
  for (subdivision in seq_len(subdivisions)) {
    F.geometry <- subdivide.mesh(F.geometry)
    F.geometry.A <- c(F.geometry.A, list(F.geometry$A))
  }

  if (method == "log") {
    F.zero <- -Inf
    F <- log(F)
    F[is.infinite(F) & F < 0] <- F.zero
  } else if (method == "linear") {
    F.zero <- 0
  } else {
    ## 'step'
    F.zero <- 0
  }

  ## For ordinary excursions, the input set/group information is not used.
  if (type == ">") {
    G <- rep(1, length(G))
  } else if (type == "<") {
    G <- rep(0, length(G))
  }
  ## Copy 'G' and interpolate 'F' within coherent single level regions.
  G.interp <- G
  F.interp <- F
  for (subdivision in seq_len(subdivisions)) {
    G.input <- G.interp
    F.input <- F.interp
    G.interp <- rep(-1, nrow(F.geometry.A[[subdivision]]))
    F.interp <- rep(F.zero, nrow(F.geometry.A[[subdivision]]))
    for (k in unique(G[G >= 0])) {
      ok.in <- (G.input == k)
      ok.out <- (rowSums(F.geometry.A[[subdivision]][, ok.in, drop = FALSE]) >
        1 - 1e-12)
      G.interp[ok.out] <- k
    }
    ok.in <- (G.input >= 0)
    ok.out <- (G.interp >= 0)
    if (method == "step") {
      for (vtx in which(ok.out)) {
        F.interp[vtx] <- min(F.input[F.geometry.A[[subdivision]][vtx, ] > 0])
      }
    } else {
      F.interp[ok.out] <-
        as.vector(F.geometry.A[[subdivision]][ok.out, ok.in, drop = FALSE] %*%
          F.input[ok.in])
    }
  }

  if (method == "log") {
    F <- exp(F.interp)
  } else if (method == "linear") {
    F <- F.interp
  } else {
    ## 'step'
    F <- F.interp
  }
  F[G.interp == -1] <- 0
  G <- G.interp

  F.geometry <-
    fmesher::fm_rcdt_2d_inla(
      loc = F.geometry$loc,
      tv = F.geometry$graph$tv
    )
  ## Handle possible node reordering in fm_rcdt_2d_inla()
  F[F.geometry$idx$loc] <- F
  G[F.geometry$idx$loc] <- G

  list(F = F, G = G, F.geometry = F.geometry)
}








## In-filled points at transitions should have G[i] == -1
## to get a conservative approximation.
## To get only an "over/under set", use a constant non-negative integer G
##   and let calc.complement=FALSE
probabilitymap <-
  function(mesh, F, level, G,
           calc.complement = TRUE,
           tol = 1e-7,
           output = c("sp", "fm", "inla.mesh.segment"),
           method, ...) {
    output <- match.arg(output)
    if (output == "inla.mesh.segment") {
      output <- "fm"
    }

    if (output == "sp") {
      if (!requireNamespace("sp", quietly = TRUE)) {
        stop("The 'sp' package is needed.")
      }
    }

    if (calc.complement &&
      (output == "fm")) {
      stop("Output format '", output, "' not supported for 'calc.complement = TRUE'.")
    }

    spout <- list()
    inlaout <- list()

    ## Find individual avoidance/between-level/under/over sets.
    for (k in sort(unique(G[G >= 0]))) {
      active.triangles <-
        which((rowSums(matrix(
          G[mesh$graph$tv] == k,
          nrow(mesh$graph$tv), 3
        )) == 3) &
          (rowSums(matrix(
            is.finite(F[mesh$graph$tv]),
            nrow(mesh$graph$tv), 3
          )) == 3))

      if (length(active.triangles) > 0) {
        submesh <- submesh.mesh.tri(active.triangles, mesh)
        active.nodes.idx <- which(!is.na(submesh$idx$loc))
        subF <- rep(NA, length(active.nodes.idx))
        subF[submesh$idx$loc[active.nodes.idx]] <- F[active.nodes.idx]

        if (nrow(submesh$loc) > 0) {
          ##      submesh$n <- nrow(submesh$loc)
          ##      submesh$graph <-
          ##        generate_trigraph_properties(
          ##          list(tv=submesh$graph$tv),
          ##          Nv=nrow(submesh$loc))

          #      if (FALSE) { ## Debugging plots
          #        op <- par(mfrow=c(2,1))
          #        on.exit(par(op))

          #        class(mesh) <- "inla.mesh"
          #        mesh$n <- nrow(mesh$loc)
          #        mesh$manifold <- "R2"
          #        proj <- inla.mesh.projector(mesh)
          #        image(proj$x, proj$y, inla.mesh.project(proj, field=exp(F)),
          #              zlim=range(exp(F)))
          #        if (length(spout) > 0) {
          #          plot(sp::SpatialPolygons(spout), add=TRUE, col="blue")
          #        }
          #        proj <- inla.mesh.projector(submesh)
          #        image.plot(proj$x, proj$y, inla.mesh.project(proj, field=exp(subF)),
          #                   xlim=range(mesh$loc[,1]), ylim=range(mesh$loc[,1]),
          #                   zlim=range(exp(F)))
          #        plot(submesh, add=TRUE)
          #      }

          subF[is.na(subF)] <- 0

          if (method == "step") {
            tric <- tricontour_step(
              x = submesh$graph, z = subF, levels = level,
              loc = submesh$loc
            )
          } else {
            tric <- tricontour(
              x = submesh$graph, z = subF, levels = level,
              loc = submesh$loc, type = "+", tol = tol, ...
            )
          }
          ID <- as.character(k)

          if (output == "sp" || calc.complement) {
            spobj <- tryCatch(
              as.sp.outline(tric,
                grp.ccw = c(2, 3),
                grp.cw = c(),
                ccw = FALSE,
                closed = TRUE,
                ID = ID
              ),
              error = function(e) NULL
            )
            if (!is.null(spobj)) {
              if (spobj@area == 0) {
                warning("Skipping zero area polygon in probabilitymap.")
              } else {
                spout[[ID]] <- spobj
              }
            }
          }
          if (output == "fm") {
            inlaout[[ID]] <-
              as.fm_segm.outline(tric,
                grp.ccw = c(2, 3),
                grp.cw = c(),
                grp = k
              )
          }
        }
      }
    }

    if (calc.complement) {
      ## Find contour set
      if (!requireNamespace("sf", quietly = TRUE)) {
        stop("Package 'sf' required for set complement calculations.")
      }

      ID <- "-1"
      outline <- fmesher::fm_segm(mesh, boundary = TRUE)
      sp.domain <- as.sp.outline(outline,
        grp.ccw = unique(outline$grp),
        grp.cw = integer(0),
        ID = ID,
        closed = TRUE
      )
      sp.domain <- sp::SpatialPolygons(list(sp.domain))

      if (length(spout) == 0) {
        ## Complement is the entire domain
        spout[[ID]] <- sp.domain@polygons[[1]]
      } else {
        spout.joined <- sp::SpatialPolygons(spout)
        spout.union <- sf::st_union(sf::st_as_sfc(spout.joined))
        sp.diff <- sf::st_difference(sf::st_as_sfc(sp.domain), spout.union)
        sp.diff <- sf::as_Spatial(sp.diff)
        if (!is.null(sp.diff)) {
          spout[[ID]] <- sp.diff@polygons[[1]]
        }
      }
      if (!is.null(spout[[ID]])) {
        spout[[ID]]@ID <- ID

        if (output == "fm") {
          inlaout[[ID]] <- fmesher::fm_as_segm(spout[[ID]])
        }
      }
    }

    if ((output == "sp") && (length(spout) > 0)) {
      out <- sp::SpatialPolygons(spout)
    } else if ((output == "fm") && (length(inlaout) > 0)) {
      out <- do.call(fm_segm, do.call(c, inlaout))
    } else {
      out <- NULL
    }

    out
  }




gaussquad <- function(mesh, method = c("direct", "make.A")) {
  method <- match.arg(method)
  if (mesh$manifold == "M" && method == "make.A") {
    stop("Guassian quadrature method 'make.A' not supported for 'M' manifolds.")
  }

  ## Construct cubic Gauss-quadrature points and weights
  nT <- nrow(mesh$graph$tv)
  I.w <- fmesher::fm_fem(mesh, order = 0)$ta
  if (method == "make.A") {
    ## Old method, kept for sanity checking.
    I.loc <-
      rbind(
        (mesh$loc[mesh$graph$tv[, 1], , drop = FALSE] +
          mesh$loc[mesh$graph$tv[, 2], , drop = FALSE] +
          mesh$loc[mesh$graph$tv[, 3], , drop = FALSE]) / 3,
        (mesh$loc[mesh$graph$tv[, 1], , drop = FALSE] * 3 +
          mesh$loc[mesh$graph$tv[, 2], , drop = FALSE] +
          mesh$loc[mesh$graph$tv[, 3], , drop = FALSE]) / 5,
        (mesh$loc[mesh$graph$tv[, 1], , drop = FALSE] +
          mesh$loc[mesh$graph$tv[, 2], , drop = FALSE] * 3 +
          mesh$loc[mesh$graph$tv[, 3], , drop = FALSE]) / 5,
        (mesh$loc[mesh$graph$tv[, 1], , drop = FALSE] +
          mesh$loc[mesh$graph$tv[, 2], , drop = FALSE] +
          mesh$loc[mesh$graph$tv[, 3], , drop = FALSE] * 3) / 5
      )
    if (mesh$manifold == "S2") {
      I.loc <- I.loc / sqrt(rowSums(I.loc^2))
    }
    A <- fmesher::fm_basis(mesh, I.loc)
  } else {
    A <- sparseMatrix(
      i = (rep(seq_len(nT), times = 3 * 4) +
        rep(c(0, 1, 2, 3) * nT, each = nT * 3)),
      j = rep(
        c(
          mesh$graph$tv[, 1],
          mesh$graph$tv[, 2],
          mesh$graph$tv[, 3]
        ),
        times = 4
      ),
      x = rep(
        c(
          c(1, 1, 1) / 3,
          c(3, 1, 1) / 5,
          c(1, 3, 1) / 5,
          c(1, 1, 3) / 5
        ),
        each = nT
      ),
      dims = c(nT * 4, mesh$n)
    )
  }

  I.w <- (rep(I.w, times = 4) *
    rep(c(-27, 25, 25, 25) / 48, each = nT))

  list(A = A, w = I.w)
}


calc.continuous.P0 <- function(F, G, F.geometry, method) {
  tri <-
    which((G[F.geometry$graph$tv[, 1]] >= 0) &
      (G[F.geometry$graph$tv[, 1]] == G[F.geometry$graph$tv[, 2]]) &
      (G[F.geometry$graph$tv[, 1]] == G[F.geometry$graph$tv[, 3]]) &
      (rowSums(matrix(
        is.finite(F[F.geometry$graph$tv]),
        nrow(F.geometry$graph$tv), 3
      )) == 3))
  submesh <- submesh.mesh.tri(tri, F.geometry)
  active.nodes.idx <- which(!is.na(submesh$idx$loc))
  subF <- rep(NA, length(active.nodes.idx))
  subF[submesh$idx$loc[active.nodes.idx]] <- F[active.nodes.idx]

  tot.area <- sum(fmesher::fm_fem(F.geometry, order = 0)$ta)

  if (method == "linear") {
    I.w <- fmesher::fm_fem(submesh, order = 0)$va
    P0 <- sum(I.w * subF) / tot.area
  } else if (method == "log") {
    II <- gaussquad(submesh)
    tmp <- exp(as.vector(II$A %*% log(subF)))
    P0 <- sum(II$w * tmp) / tot.area
  } else {
    ## "step"
    I.w <- fmesher::fm_fem(submesh, order = 0)$ta
    tmp <- matrix(
      subF[submesh$graph$tv],
      nrow(submesh$graph$tv),
      ncol(submesh$graph$tv)
    )
    tmp <- apply(tmp, 1, min)
    P0 <- sum(I.w * tmp) / tot.area
  }
  P0
}


#' Calculate continuous domain excursion and credible contour sets
#'
#' Calculates continuous domain excursion and credible contour sets
#'
#' @param ex An `excurobj` object generated by a call to [excursions()]
#' or [contourmap()].
#' @param geometry Specification of the lattice or triangulation geometry of the input.
#' One of `list(x, y)`, `list(loc, dims)`, `fm_lattice_2d`,
#' `inla.mesh.lattice`, `fm_mesh_2d`, or
#' `inla.mesh`, where `x` and `y` are vectors, `loc` is
#' a two-column matrix of coordinates, and `dims` is the lattice size vector.
#' The first three versions are all treated topologically as lattices, and the
#' lattice boxes are assumed convex.
#' @param alpha The target error probability.  A warning is given if it is detected
#' that the information `ex` isn't sufficient for the given `alpha`.
#' Defaults to the value used when calculating `ex`.
#' @param method The spatial probability interpolation transformation method to use.
#' One of `log`, `linear`, or `step`.  For `log`, the probabilities
#' are interpolated linearly in the transformed scale. For `step`, a conservative
#' step function is used.
#' @param output Specifies what type of object should be generated. `sp` gives a
#' `SpatialPolygons` object, and `fm` or `inla` gives a `fm_segm` object.
#' @param subdivisions The number of mesh triangle subdivisions to perform for the
#' interpolation of the excursions or contour function. 0 is no subdivision.
#' The setting has a small effect on the evaluation of `P0` for the `log`
#' method (higher values giving higher accuracy) but the main effect is on the visual
#' appearance of the interpolation. Default=1.
#' @param calc.credible Logical, if TRUE (default), calculate credible contour region
#' objects in addition to avoidance sets.
#'
#' @return A list with elements
#' \item{M}{`SpatialPolygons` or `fm_segm` object. The subsets
#' are tagged, so that credible regions are tagged `"-1"`, and regions between
#' levels are tagged `as.character(0:nlevels)`.}
#' \item{F}{Interpolated F function.}
#' \item{G}{Contour and inter-level set indices for the interpolation.}
#' \item{F.geometry}{Mesh geometry for the interpolation.}
#' \item{P0}{P0 measure based on interpolated F function (only for `contourmap`
#' input).}
#' @export
#' @author Finn Lindgren \email{finn.lindgren@gmail.com}
#' @references Bolin, D. and Lindgren, F. (2017) *Quantifying the uncertainty of contour maps*, Journal of Computational and Graphical Statistics, vol 26, no 3, pp 513-524.
#'
#' Bolin, D. and Lindgren, F. (2018), *Calculating Probabilistic Excursion Sets and Related Quantities Using excursions*, Journal of Statistical Software, vol 86, no 1, pp 1-20.
#'
#' @examples
#' \dontrun{
#' if (require("fmesher") &&
#'   require("sp")) {
#'   # Generate mesh and SPDE model
#'   n.lattice <- 10 # Increase for more interesting, but slower, examples
#'   x <- seq(from = 0, to = 10, length.out = n.lattice)
#'   lattice <- fm_lattice_2d(x = x, y = x)
#'   mesh <- fm_rcdt_2d_inla(lattice = lattice, extend = FALSE, refine = FALSE)
#'
#'   # Generate an artificial sample
#'   sigma2.e <- 0.1
#'   n.obs <- 100
#'   obs.loc <- cbind(
#'     runif(n.obs) * diff(range(x)) + min(x),
#'     runif(n.obs) * diff(range(x)) + min(x)
#'   )
#'   Q <- fm_matern_precision(mesh, alpha = 2, rho = 3, sigma = 1)
#'   x <- fm_sample(n = 1, Q = Q)
#'   A <- fm_basis(mesh, loc = obs.loc)
#'   Y <- as.vector(A %*% x + rnorm(n.obs) * sqrt(sigma2.e))
#'
#'   ## Calculate posterior
#'   Q.post <- (Q + (t(A) %*% A) / sigma2.e)
#'   mu.post <- as.vector(solve(Q.post, (t(A) %*% Y) / sigma2.e))
#'   vars.post <- excursions.variances(chol(Q.post), max.threads = 1)
#'
#'   ## Calculate contour map with two levels
#'   map <- contourmap(
#'     n.levels = 2, mu = mu.post, Q = Q.post,
#'     alpha = 0.1, F.limit = 0.1, max.threads = 1
#'   )
#'
#'   ## Calculate the continuous representation
#'   sets <- continuous(map, mesh, alpha = 0.1)
#'
#'   ## Plot the results
#'   reo <- mesh$idx$lattice
#'   cols <- contourmap.colors(map,
#'     col = heat.colors(100, 1, rev = TRUE),
#'     credible.col = grey(0.5, 1)
#'   )
#'   names(cols) <- as.character(-1:2)
#'
#'   par(mfrow = c(2, 2))
#'   image(matrix(mu.post[reo], n.lattice, n.lattice),
#'     main = "mean", axes = FALSE, asp = 1
#'   )
#'   image(matrix(sqrt(vars.post[reo]), n.lattice, n.lattice),
#'     main = "sd", axes = FALSE, asp = 1
#'   )
#'   image(matrix(map$M[reo], n.lattice, n.lattice),
#'     col = cols, axes = FALSE, asp = 1
#'   )
#'   idx.M <- setdiff(names(sets$M), "-1")
#'   plot(sets$M[idx.M], col = cols[idx.M])
#' }
#' }
#'
#' @importFrom lifecycle deprecated
continuous <- function(ex,
                       geometry,
                       alpha,
                       method = c("log", "linear", "step"),
                       output = c("sp", "fm", "inla"),
                       subdivisions = 1,
                       calc.credible = TRUE) {
  stopifnot(inherits(ex, "excurobj"))
  method <- match.arg(method)
  output <- match.arg(output)
  if (output == "inla") {
    lifecycle::deprecate_warn(
      "2.5.8.9001",
      'continuous(output = "inla")',
      'continuous(output = "mf")',
      "The 'inla' output format is deprecated and may be removed in the future. Use 'fm' instead.")
    output <- "fm"
  }

  if (!(ex$meta$calculation %in% c(
    "excursions",
    "contourmap"
  ))) {
    stop(paste("Unsupported calculation '",
      ex$meta$calculation, "'.",
      sep = ""
    ))
  }

  if (missing(alpha)) {
    alpha <- ex$meta$alpha
  }
  if (alpha > ex$meta$F.limit) {
    warning(paste("Insufficient data: alpha = ", alpha,
      " > F.limit = ", ex$meta$F.limit,
      sep = ""
    ))
  }

  info <- get.geometry(geometry)
  if (!fmesher::fm_manifold(info, c("M", "R2", "S2"))) {
    stop(paste("Unsupported manifold type '", fmesher::fm_manifold(info), "'.", sep = ""))
  }
  if ((output == "sp") && !fmesher::fm_manifold(info, c("R2"))) {
    stop(paste("Unsupported manifold type '", fmesher::fm_manifold(info), "' for 'sp' output.", sep = ""))
  }
  if (calc.credible &&
    (output == "fm")) {
    stop("Output format 'fm' not supported for 'calc.credible = TRUE'.")
  }

  if (length(ex$F) != prod(info$dims)) {
    stop(paste("The number of computed F-values (", length(ex$F), ") must match \n",
      "the number of elements of the continuous geometry definition (",
      prod(info$dims), ").",
      sep = ""
    ))
  }

  if (ex$meta$type == "=") {
    type <- "!="
    F.ex <- 1 - ex$F
  } else {
    type <- ex$meta$type
    F.ex <- ex$F
  }
  F.ex[is.na(F.ex)] <- 0

  if (is.null(ex$meta$ind)) {
    active.nodes <- rep(TRUE, length(ex$F))
  } else {
    active.nodes <- logical(length(ex$F))
    active.nodes[ex$meta$ind] <- TRUE
  }
  if (info$geometry == "mesh") {
    mesh <- submesh.mesh(active.nodes, geometry)
  } else if (info$geometry == "lattice") {
    mesh <- submesh.grid(active.nodes, geometry)
  }
  mesh$graph <-
    generate_trigraph_properties(mesh$graph, Nv = nrow(mesh$loc))

  active.nodes <- !is.na(mesh$idx$loc)
  F.ex[mesh$idx$loc[active.nodes]] <- F.ex[active.nodes]
  G.ex <- rep(-1, nrow(mesh$loc))
  G.ex[mesh$idx$loc[active.nodes]] <- ex$G[active.nodes]

  ## Construct interpolation
  interpolated <-
    F.interpolation(mesh, F.ex, G.ex, ex$meta$type, method,
      subdivisions = subdivisions
    )

  if (ex$meta$type == "=") {
    interpolated$F <- 1 - interpolated$F
    ##    F.ex <- 1-F.ex   ????????????????????
  }

  F.geometry <- mesh
  if (method == "log") {
    F.ex <- log(F.ex)
    level <- log(1 - alpha)
    F.ex[!is.finite(F.ex)] <- NA
  } else if (method == "linear") {
    level <- 1 - alpha
  } else {
    ## 'step'
    level <- 1 - alpha
  }

  ## For ordinary excursions, the input set/group information is not used.
  if (type == ">") {
    G.ex <- rep(1, length(G.ex))
  } else if (type == "<") {
    G.ex <- rep(0, length(G.ex))
  }

  M <- probabilitymap(F.geometry,
    F = F.ex,
    level = level,
    G = G.ex,
    calc.complement = calc.credible,
    method = method,
    output = output
  )

  if (method == "log") {
    F.ex <- exp(F.ex)
  }

  F.ex[is.na(F.ex)] <- 0
  interpolated$F[is.na(interpolated$F)] <- 0
  out <- list(
    F = interpolated$F, G = interpolated$G,
    M = M, F.geometry = interpolated$F.geometry
  )

  if (!is.null(ex$P0)) {
    out$P0 <- calc.continuous.P0(
      interpolated$F, interpolated$G,
      interpolated$F.geometry, method
    )
  }

  out
}