neighbourhood.R

if (getRversion() >= "3.1.0") {
  utils::globalVariables(c("angles", "indices", "to", "x", "y", "rasterVal"))
}

##############################################################
#' Fast `adjacent` function, and Just In Time compiled version
#'
#' Faster function for determining the cells of the 4, 8 or bishop
#'  neighbours of the \code{cells}. This is a hybrid function that uses
#'  matrix for small numbers of loci (<1e4) and data.table for larger numbers of loci
#'
#' Between 4x (large number loci) to 200x (small number loci) speed gains over
#' \code{adjacent} in raster package. There is some extra speed gain if
#' \code{NumCol} and \code{NumCells} are passed rather than a raster.
#' Efficiency gains come from:
#'  1. use \code{data.table} internally
#'     - no need to remove NAs because wrapped or outside points are
#'       just removed directly with data.table
#'     - use data.table to sort and fast select (though not fastest possible)
#'  2. don't make intermediate objects; just put calculation into return statement
#'
#' The steps used in the algorithm are:
#' 1. Calculate indices of neighbouring cells
#' 2. Remove "to" cells that are
#'    - <1 or >numCells (i.e., they are above or below raster), using a single modulo
#'      calculation
#'    - where the modulo of "to" cells is equal to 1 if "from" cells are 0 (wrapped right
#'      to left)
#'    - or where the modulo of the "to" cells is equal to 0 if "from" cells are 1 (wrapped
#'      left to right)
#'
#' @param x \code{Raster*} object for which adjacency will be calculated.
#'
#' @param cells vector of cell numbers for which adjacent cells should be found.
#'              Cell numbers start with 1 in the upper-left corner and increase
#'              from left to right and from top to bottom.
#'
#' @param directions the number of directions in which cells should be connected:
#'                   4 (rook's case), 8 (queen's case), or \code{"bishop"} to connect
#'                   cells with one-cell diagonal moves.
#'                   Or a neighbourhood matrix (see Details).
#'
#' @param sort logical. Whether the outputs should be sorted or not, using cell ids
#'             of the \code{from} cells (and \code{to} cells, if \code{match.adjacent}
#'             is \code{TRUE}).
#'
#' @param pairs logical. If \code{TRUE}, a matrix of pairs of adjacent cells is returned.
#'              If \code{FALSE}, a vector of cells adjacent to cells is returned
#'
#' @param include logical. Should the focal cells be included in the result?
#'
#' @param target a vector of cells that can be spread to. This is the inverse of a mask.
#'
#' @param numCol numeric indicating number of columns in the raster.
#'               Using this with \code{numCell} is a bit faster execution time.
#'
#' @param numCell numeric indicating number of cells in the raster.
#'                Using this with \code{numCol} is a bit faster execution time.
#'
#' @param match.adjacent logical. Should the returned object be the same as
#'                       \code{raster::adjacent}.
#'                       Default \code{FALSE}, which is faster.
#'
#' @param cutoff.for.data.table numeric. If the number of cells is above this value,
#'                              the function uses data.table which is faster with
#'                              large numbers of cells. Default is 5000, which appears
#'                              to be the turning point where data.table becomes faster.
#'
#' @param torus Logical. Should the spread event wrap around to the other side of the raster?
#'                       Default is \code{FALSE}.
#'
#' @param id numeric If not \code{NULL} (default), then function will return \code{"id"} column.
#'
#' @param numNeighs A numeric scalar, indicating how many neighbours to return. Must be
#'                  less than or equal to \code{directions}; which neighbours are random
#'                  with equal probabilities.
#' @param returnDT A logical. If TRUE, then the function will return the result
#'                 as a \code{data.table}, if the internals used \code{data.table},
#'                 i.e., if number of cells is greater than \code{cutoff.for.data.table}.
#'                 User should be warned that this will therefore cause the output
#'                 format to change depending \code{cutoff.for.data.table}.
#'                 This will be faster for situations where \code{cutoff.for.data.table = TRUE}.
#'
#' @return Either a matrix (if more than 1 column, i.e., \code{pairs = TRUE},
#' and/or \code{id} is provided), a vector (if only one column), or a \code{data.table}
#' (if \code{cutoff.for.data.table} is less than \code{length(cells)} \emph{and}
#' \code{returnDT} is \code{TRUE}.
#' To get a consistent output, say a matrix, it would be wise to test the output
#' for its class.
#' The variable output is done to minimize coercion to maintain speed.
#' The columns will be one or more of \code{id}, \code{from}, \code{to}.
#'
#' @seealso \code{\link[raster]{adjacent}}
#'
#' @author Eliot McIntire
#' @export
#' @importFrom data.table := data.table key set setcolorder setkeyv
#' @importFrom raster ncell ncol nrow
#' @importFrom stats na.omit
#' @rdname adj
#'
#' @examples
#' library(raster)
#' a <- raster(extent(0, 1000, 0, 1000), res = 1)
#' sam <- sample(1:length(a), 1e4)
#' numCol <- ncol(a)
#' numCell <- ncell(a)
#' adj.new <- adj(numCol = numCol, numCell = numCell, cells = sam, directions = 8)
#' adj.new <- adj(numCol = numCol, numCell = numCell, cells = sam, directions = 8,
#'                include = TRUE)
#'
adj <- function(x = NULL, cells, directions = 8, sort = FALSE, pairs = TRUE,
                    include = FALSE, target = NULL, numCol = NULL, numCell = NULL,
                    match.adjacent = FALSE, cutoff.for.data.table = 2e3,
                    torus = FALSE, id = NULL, numNeighs = NULL, returnDT = FALSE) {
  to <- NULL
  J <- NULL # nolint
  cells <- as.integer(cells)

  if (is.null(numCol) | is.null(numCell)) {
    if (is.null(x)) stop("must provide either numCol & numCell or a x")
    numCol <- as.integer(ncol(x))
    numCell <- as.integer(ncell(x))
  }

  if (directions == "bishop")  {
    dirs <- 4
    needCorners <- TRUE
  } else {
    needCorners <- if (directions == 8) TRUE else FALSE
    dirs <- directions
  }

  numToCells <- dirs + include
  fromCells <- rep.int(cells, times = numToCells)

  if (is.numeric(directions)) {
    top <- cells - numCol
    lef <- cells - 1L
    rig <- cells + 1L
    bot <- cells + numCol
  }
  if (needCorners) {
    topl <- cells - numCol - 1L
    topr <- cells - numCol + 1L
    botl <- cells + numCol - 1L
    botr <- cells + numCol + 1L
  }

  toCells <- if (directions == 8) {

    if (match.adjacent)
      if (include)
        c(cells, topl, lef, botl, topr, rig, botr, top, bot)
      else
        c(topl, lef, botl, topr, rig, botr, top, bot)
    else
      if (include)
        c(topl, top, topr, lef, cells, rig, botl, bot, botr)
      else
        c(topl, top, topr, lef, rig, botl, bot, botr)
  } else if (directions == 4) {
    if (match.adjacent)
      if (include)
        c(cells, lef, rig, top, bot)
      else
        c(lef, rig, top, bot)
    else
      if (include)
        c(top, lef, cells, rig, bot)
      else
        c(top, lef, rig, bot)
  } else if (directions == "bishop") {
    if (match.adjacent)
      if (include)
        c(cells, topl, botl, topr, botr)
      else
        c(topl, botl, topr, botr)
    else
      if (include)
        c(topl, topr, cells, botl, botr)
      else
        c(topl, topr, botl, botr)
  } else {
    stop("directions must be 4 or 8 or \'bishop\'")
  }

  if (!is.null(numNeighs)) {
    lenCells <- length(cells)
    if (length(numNeighs) == 1) numNeighs <- rep(numNeighs, lenCells)
    ind <- unlist(sampleV(1:(directions + include), size = numNeighs))
    minusVal <- lenCells - rep.int(seq_along(cells), numNeighs)
    indFull2 <- ind * lenCells - minusVal

    toCells <- toCells[indFull2]
    fromCells <- fromCells[indFull2]
  }

  useMatrix <- (length(cells) < cutoff.for.data.table)
  if (useMatrix) {
    adj <- cbind(from = fromCells, to = toCells)
    if (!is.null(id)) adj <- cbind(adj, id = rep.int(id, times = numToCells))
  } else {
    adj <- data.table(from = fromCells, to = toCells)
    if (!is.null(id)) set(adj, , "id", rep.int(id, times = numToCells))
  }

  if (useMatrix) {
    # Remove all cells that are not target cells, if target is a vector of cells
    if (!is.null(target)) {
      adj <- adj[na.omit(adj[, "to"] %in% target), , drop = FALSE]
    }

    if (sort) {
      if (pairs) {
        if (match.adjacent) {
          adj <- adj[order(adj[, "from"], adj[, "to"]), , drop = FALSE]
        } else {
          adj <- adj[order(adj[, "from"]), , drop = FALSE]
        }
      } else {
        adj <- adj[order(adj[, "to"]), , drop = FALSE]
      }
    }

    # Remove the "from" column if pairs is FALSE
    # Good time savings if no intermediate object is created
    keepCols <- if (is.null(id)) "to" else c("to", "id")
    if (!torus) {
      if (pairs) {
        return(adj[
          !((((adj[, "to"] - 1) %% numCell + 1) != adj[, "to"]) | # top or bottom of raster
              ((adj[, "from"] %% numCol + adj[, "to"] %% numCol) == 1)) # | #right & left edge cells, with neighbours wrapped
          , , drop = FALSE])
      } else {
        adj <- adj[
          !((((adj[, "to"] - 1) %% numCell + 1) != adj[, "to"]) | # top or bottom of raster
              ((adj[, "from"] %% numCol + adj[, "to"] %% numCol) == 1)) # | #right & left edge cells, with neighbours wrapped
          , keepCols, drop = FALSE]
        if (match.adjacent) {
          adj <- unique(adj[, "to"])
        }
        return(adj)
      }
    } else {
      whLefRig <- (adj[, "from"] %% numCol + adj[, "to"] %% numCol) == 1
      adj[whLefRig, "to"] <- adj[whLefRig, "to"] +
        numCol * (adj[whLefRig, "from"] - adj[whLefRig, "to"])
      whBotTop <- ((adj[, "to"] - 1) %% numCell + 1) != adj[, "to"]
      adj[whBotTop, "to"] <- adj[whBotTop, "to"] +
        sign(adj[whBotTop, "from"] - adj[whBotTop, "to"]) * numCell
      if (pairs) {
        return(adj)
      } else {
        if (match.adjacent) {
          adj <- unique(adj[, "to", drop = TRUE])
        } else {
          adj <- adj[, keepCols, drop = FALSE]
        }
        return(adj)
      }
    }
  } else {
    ## use data.table
    # Remove all cells that are not target cells, if target is a vector of cells
    if (!is.null(target)) {
      set(adj, , "ord", seq_len(NROW(adj)))
      setkeyv(adj, "to")
      adj <- adj[J(target)]
      adj <- na.omit(adj)
      setkeyv(adj, "ord")
      set(adj, , "ord", NULL)
    }

    if (sort) {
      if (pairs) {
        if (match.adjacent) {
          setkeyv(adj, c("from", "to"))
        } else {
          setkeyv(adj, "from")
        }
      } else {
        setkeyv(adj, "to")
      }
    }

    # Remove the "from" column if pairs is FALSE
    if (!pairs) {
      from <- as.integer(adj$from)
      set(adj, , "from", NULL)
    }

    if (!torus) {
      if (!pairs) {
        adj <- adj[
          !((((to - 1) %% numCell + 1) != to) |  #top or bottom of raster
              ((from %% numCol + to %% numCol) == 1))# | #right & left edge cells, with neighbours wrapped
          ]
        if (match.adjacent) {
          if (returnDT)
            return(unique(adj[, list(to)]))
          else
            return(unique(adj$to))
        }
        if (returnDT)
          return(adj)
        else
          return(as.matrix(adj))
      } else {
        if (returnDT)
          return(adj[
            # top or bottom of raster | right & left edge cells, with neighbours wrapped
            !((((to - 1) %% numCell + 1) != to) | ((from %% numCol + to %% numCol) == 1))
          ])
        else
          return(as.matrix(adj[
            # top or bottom of raster | right & left edge cells, with neighbours wrapped
            !((((to - 1) %% numCell + 1) != to) | ((from %% numCol + to %% numCol) == 1))
          ]))
        return()
      }
    } else {
      if (!pairs) {
        whLefRig <- (from %% numCol + adj$to %% numCol) == 1
        toWhLefRig <- adj$to[whLefRig]
        set(adj, which(whLefRig), "to", toWhLefRig + numCol * (from[whLefRig] - toWhLefRig))
        whBotTop <- ((adj$to - 1) %% numCell + 1) != adj$to
        toWhBotTop <- adj$to[whBotTop]
        set(adj, which(whBotTop), "to", toWhBotTop +
              as.integer(sign(from[whBotTop] - toWhBotTop) * numCell))

        if (match.adjacent) {
          if (returnDT)
            return(unique(adj[, list(to)]))
          else
            return(unique(adj$to))
        }
      } else {
        whLefRig <- (adj$from %% numCol + adj$to %% numCol) == 1
        toWhLefRig <- adj$to[whLefRig]
        set(adj, which(whLefRig), "to", toWhLefRig + numCol * (adj$from[whLefRig] - toWhLefRig))
        whBotTop <- ((adj$to - 1) %% numCell + 1) != adj$to
        toWhBotTop <- adj$to[whBotTop]
        set(adj, which(whBotTop), "to", toWhBotTop +
              as.integer(sign(adj$from[whBotTop] - toWhBotTop) * numCell))
      }
      if (returnDT)
        return(adj)
      else
        return(as.matrix(adj))
    }
  }
}

##############################################################
#' Identify pixels in a circle or ring (donut) around an object.
#'
#' Identify the pixels and coordinates that are at a (set of) buffer distance(s)
#' of the objects passed into \code{coords}.
#' This is similar to \code{rgeos::gBuffer} but much faster and without
#' the geo referencing information.
#' In other words, it can be used for similar problems, but where speed is important.
#' This code is substantially adapted from \code{PlotRegionHighlighter::createCircle}.
#'
#' @param landscape    Raster on which the circles are built.
#'
#' @param coords Either a matrix with 2 (or 3) columns, x and y (and id), representing the
#'               coordinates (and an associated id, like cell index),
#'               or a \code{SpatialPoints*} object around which to make circles. Must be same
#'               coordinate system as the \code{landscape} argument. Default is missing,
#'               meaning it uses the default to \code{loci}
#'
#' @param loci   Numeric. An alternative to \code{coords}. These are the indices on
#'               \code{landscape} to initiate this function. See \code{coords}. Default is one
#'               point in centre of \code{landscape}..
#'
#' @param maxRadius  Numeric vector of length 1 or same length as coords
#'
#' @param minRadius  Numeric vector of length 1 or same length as \code{coords}. Default is
#'                   \code{maxRadius}, meaning return all cells that are touched
#'                   by the narrow ring at that exact radius. If smaller than \code{maxRadius},
#'                   then this will create a buffer or donut or ring.
#'
#' @param allowOverlap Logical. Should duplicates across id be removed or kept. Default TRUE.
#'
#' @param allowDuplicates Logical. Should duplicates within id be removed or kept. Default FALSE.
#'                        This is useful if the actual x, y coordinates are desired, rather
#'                        than the cell indices. This will increase the size of the returned
#'                        object.
#'
#' @param includeBehavior Character string. Currently accepts only "includePixels", the default,
#'                        and "excludePixels". See details.
#'
#' @param returnDistances Logical. If TRUE, then a column will be added to the returned
#'                        data.table that reports the distance from \code{coords} to every
#'                        point that was in the circle/donut surrounding \code{coords}. Default
#'                        FALSE, which is faster.
#'
#' @param angles Numeric. Optional vector of angles, in radians, to use. This will create
#'               "spokes" outward from coords. Default is NA, meaning, use internally
#'               derived angles that will "fill" the circle.
#'
#' @param returnAngles Logical. If TRUE, then a column will be added to the returned
#'                        data.table that reports the angle from \code{coords} to every
#'                        point that was in the circle/donut surrounding \code{coords}. Default
#'                        FALSE.
#'
#' @param closest Logical. When determining non-overlapping circles, should the function
#'                give preference to the closest \code{loci} or the first one (much faster).
#'                Default is FALSE, meaning the faster, though maybe not desired behaviour.
#'
#' @param simplify logical. If TRUE, then all duplicate pixels are removed. This means
#' that some x, y combinations will disappear.
#'
#' @inheritParams spread
#'
#' @details This function identifies all the pixels as defined by a donut
#' with inner radius minRadius and outer radius of maxRadius. The includeBehavior defines
#' whether the cells that intersect the radii but whose centres are not inside
#' the donut are included \code{includePixels} or not \code{excludePixels} in the returned
#' pixels identified. If this is \code{excludePixels}, and if a \code{minRadius} and
#' \code{maxRadius} are equal, this will return no pixels.
#'
#'
#' @return A \code{matrix} with 4 columns, \code{id}, \code{indices},
#' \code{x}, \code{y}. The \code{x} and \code{y} indicate the
#' exact coordinates of
#' the \code{indices} (i.e., cell number) of the \code{landscape}
#' associated with the ring or circle being identified by this function.
#'
#' @importFrom data.table data.table set setkeyv
#' @importFrom sp coordinates
#' @importFrom fpCompare %==%
#' @importFrom raster cellFromXY extract res xyFromCell ncell ncol
#' @export
#' @rdname cir
#' @seealso \code{\link{rings}} which uses \code{spread} internally.
#' \code{cir} tends to be faster when there are few starting points, \code{rings}
#' tends to be faster when there are many starting points. \code{cir} scales with
#' \code{maxRadius} ^ 2 and \code{coords}. Another difference
#' between the two functions is that \code{rings} takes the centre of the pixel
#' as the centre of a circle, whereas \code{cir} takes the exact coordinates.
#' See example. For the specific case of creating distance surfaces from specific
#' points, see \code{\link{distanceFromEachPoint}}, which is often faster.
#' For the more general GIS buffering, see \code{rgeos::gBuffer}.
#'
#'@example inst/examples/example_cir.R
#'
#' @export
#' @rdname cir
# setMethod(
#   "cir",
#   signature(landscape = "RasterLayer", coords = "matrix", loci = "missing"),
#   definition =
cir <- function(landscape, coords, loci,
                maxRadius = ncol(landscape) / 4, minRadius = maxRadius,
                allowOverlap = TRUE, allowDuplicates = FALSE,
                includeBehavior = "includePixels", returnDistances = FALSE,
                angles = NA_real_,
                returnAngles = FALSE, returnIndices = TRUE,
                closest = FALSE, simplify = TRUE)  {

  if (missing(coords)) {
    if (missing(loci)) {
      ncells <- ncell(landscape)
      middleCell <- if (identical(ncells / 2, floor(ncells / 2))) {
        ncells / 2 - ncol(landscape) / 2
      } else {
        round(ncells / 2)
      }
      coords <- xyFromCell(landscape, middleCell)
    } else if (is.numeric(loci)) {
      coords <- xyFromCell(landscape, loci)
      coords <- cbind(coords, id = loci)
    } else {
      stop("Need either a numeric loci or coords")
    }
  } else if (inherits(coords, "Spatial")) {
    coords <- coordinates(coords)
  } else if (!is.numeric(coords)) {
    stop("Need either a numeric loci or coords with matrix or SpatialPoints")
  }

  ### adapted from createCircle of the package PlotRegionHighlighter

  if (!all(c("x", "y") %in% colnames(coords))) {
    stop("coords must have columns named x and y")
  }
  suppliedAngles <- if (all(!is.na(angles))) TRUE else FALSE

  scaleRaster <- res(landscape)
  if (!isTRUE(all.equal(scaleRaster[1], scaleRaster[2]))) {
    stop("cir function only accepts rasters with identical resolution in x and y dimensions")
  }

  if (!any(includeBehavior == c("includePixels", "excludePixels"))) {
    stop("includeBehavior can only be \"includePixels\" or \"excludePixels\"")
  }

  scaleRaster <- scaleRaster[1]

  moreThanOne <- NROW(coords) > 1

  equalRadii <- TRUE
  if (suppliedAngles) {
    # if provided with angles, then problem is easier
    seqNumInd <- seq_len(NROW(coords))
    maxRadius <- c(seq(minRadius, maxRadius, by = max(0.68, 0.75 - maxRadius / 3e3)), maxRadius)
    numAngles <- length(angles)
    rads <- rep(rep(maxRadius, each = numAngles), NROW(coords))
    x <- kronecker(coords[, "x"], c(cos(angles) %o% maxRadius), "+")
    y <- kronecker(coords[, "y"], c(sin(angles) %o% maxRadius), "+")
    id <- rep(rep(seqNumInd, each = numAngles), each = length(maxRadius))
  } else {
    if (moreThanOne) {
      # create an index sequence for the number of individuals
      seqNumInd <- seq_len(NROW(coords))

      if (length(maxRadius) == 1) maxRadius <- rep(maxRadius, NROW(coords))
      if (length(minRadius) == 1) minRadius <- rep(minRadius, NROW(coords))
      equalRadii <- sum(maxRadius - maxRadius[1]) %==% 0

      # The goal of maxRadius and numAngles is to identify every cell within the circle
      #  The 0.68 and 0.75 were found by trial and error to minimize the number of
      #  pixels selected that are duplicates of each other.
      if (any(minRadius != maxRadius)) {
        if (any(minRadius > maxRadius)) stop("minRadius must be less than or equal to maxRadius")
        maxRadiusList <- lapply(seqNumInd, function(x) {
          ## 0.75 was the maximum that worked with 1e4 pixels, 1e2 maxRadius
          ## 0.66 was the maximum that worked with 4e6 pixels, 1.3e3 maxRadius
          a <- seq(minRadius[x], maxRadius[x], by = max(0.68, 0.75 - maxRadius[x] / 3e3))
          if (a[length(a)] != maxRadius[x]) a <- c(a, maxRadius[x])
          a
        })

        if (equalRadii) {
          maxRadius <- do.call(cbind, maxRadiusList)
        } else {
          lengths <- unlist(lapply(maxRadiusList, length))
          maxLen <- max(lengths)
          maxRadius <- do.call(cbind, lapply(seq_along(maxRadiusList), function(y) {
            c(maxRadiusList[[y]], rep(NA_real_, maxLen - lengths[y]))
          }))
        }
      }
    } else {
      seqNumInd <- 1
      if (any(minRadius != maxRadius)) {
        ## 0.66 was the maximum that worked with 4e6 pixels, 1.3e3 maxRadius
        a <- seq(minRadius, maxRadius, by = max(0.68, 0.75 - maxRadius / 3e3))
        if (a[length(a)] != maxRadius) a <- c(a, maxRadius)
        maxRadius <- a
      }
    }

    numAngles <- ceiling(maxRadius / scaleRaster * 2.6 * pi) + 1

    if (moreThanOne) {
      if (is.matrix(numAngles)) {
        nAngles <- apply(numAngles, 2, sum, na.rm = TRUE)
      } else {
        nAngles <- numAngles
      }
    } else {
      nAngles <- sum(numAngles)
    }

    # create individual IDs for the number of points that will be done for their circle
    if (!c("id") %in% colnames(coords)) {
      if (moreThanOne) {
        id <- rep.int(seqNumInd, times = nAngles)
      } else {
        id <- 1L
      }
    } else {
      id <- as.integer(rep(coords[, "id"], times = nAngles))
    }

    # create vector of radius for the number of points that will be done for
    # each individual circle
    if (equalRadii)
      rads <- rep.int(maxRadius, times = numAngles)
    else
      rads <- rep.int(na.omit(as.vector(maxRadius)), times = na.omit(as.vector(numAngles)))

    # extract the individuals' current coords
    xs <- rep.int(coords[, "x"], times = nAngles)
    ys <- rep.int(coords[, "y"], times = nAngles)

    angles <- if (all(is.na(angles))) {
      if (!is.null(dim(numAngles))) {
        if (equalRadii) {
          rep(unlist(lapply(numAngles[, 1], function(na) {
            seq_len(na) * (pi * 2 / na)
          })), ncol(numAngles))
        } else {
          unlist(lapply(na.omit(as.vector(numAngles)), function(na) {
            seq_len(na) * (pi * 2 / na)
          }))
        }
      } else {
        unlist(lapply(numAngles, function(na) seq.int(na) * (pi * 2 / na)))
      }
    } else {
      rep(angles, length(numAngles))
    }
    x <- cos(angles) * rads + xs
    y <- sin(angles) * rads + ys
  }

  indices <- as.integer(cellFromXY(landscape, cbind(x, y)))

  if (moreThanOne & allowOverlap & !closest) {
    matDT <- data.table(id, indices, rads, angles, x = x, y = y)
    setkeyv(matDT, c("id", "indices"))
    if (!equalRadii) {
      matDT[, maxRad := rep(apply(maxRadius, 2, max, na.rm = TRUE), nAngles)]
      matDT[, minRad := rep(apply(maxRadius, 2, min, na.rm = TRUE), nAngles)]
    }
    if (!allowDuplicates) {
      matDT <- unique(matDT, by = c("id", "indices"))
    }
    matDT <- na.omit(matDT)
    matDT <- as.matrix(matDT)
  } else {
    matDT <- cbind(id, rads, angles, x, y, indices)
    if (!closest & !allowDuplicates) {
      notDups <- !duplicatedInt(indices)
      matDT <- matDT[notDups, , drop = FALSE]
    }
    matDT <- na.omit(matDT)
  }
  rm(id, indices, rads, x, y)

  # only need to calculate distances for these two cases
  if (includeBehavior == "excludePixels" | returnDistances | closest) {
    if (equalRadii) {
      maxRad <- maxRadius[NROW(maxRadius)]
      minRad <- maxRadius[1]
    }

    # if distances are not required, then only need the inner circle and outer
    # circle distances. Don't waste resources on calculating all distances.
    if (returnDistances | closest) {
      matDT2 <- matDT
    } else {
      if (equalRadii) {
        # 0.71 is the sqrt of 1, so keep
        matDT2 <- matDT[matDT[, "rads"] >= (maxRad - 0.71) | matDT[, "rads"] <=
                      (minRad + 0.71), , drop = FALSE]
      } else {
        # 0.71 is the sqrt of 1, so keep
        matDT2 <- matDT[matDT[, "rads"] >= (matDT[, "maxRad"] - 0.71) | matDT[, "rads"] <=
                      (matDT[, "minRad"] + 0.71), , drop = FALSE]
      }
    } #  only pixels that are in inner or outer ring of pixels

    if (suppliedAngles) {
      a <- cbind(id = matDT2[, "id"], rads = matDT2[, "rads"], angles = matDT2[, "angles"],
                 x = matDT2[, "x"], y = matDT2[, "y"], to = matDT2[, "indices"])

    } else {
      xyC <- xyFromCell(landscape, matDT2[, "indices"]);
      a <- cbind(id = matDT2[, "id"], rads = matDT2[, "rads"], angles = matDT2[, "angles"],
                 x = xyC[, "x"], y = xyC[, "y"], to = matDT2[, "indices"])
    }
    if (!equalRadii)
      a <- cbind(a, maxRad = matDT2[, "maxRad"], minRad = matDT2[, "minRad"])

    b <- cbind(coords, id = 1:NROW(coords))

    colnames(b)[1:2] <- c("x", "y")
    d <- distanceFromEachPoint(b, a)

    if (closest) {
      d <- d[order(d[, "rads"]), , drop = FALSE]
      dups <- duplicated(d[, "to", drop = FALSE])
      d <- d[!dups, , drop = FALSE]

    }

    if (includeBehavior == "excludePixels")
      if (equalRadii)
        d <- d[d[, "dists"] %<=% maxRad & d[, "dists"] %>=% minRad, , drop = FALSE]
    else
      d <- d[d[, "dists"] %<=% d[, "maxRad"] & d[, "dists"] %>=% d[, "minRad"], , drop = FALSE]

    colnames(d)[which(colnames(d) == "to")] <- "indices"
    if (!returnDistances)
      d <- d[, -which(colnames(d) == "dists"), drop = FALSE]

    if (!returnAngles) {
      d <- d[, -which(colnames(d) == "angles"), drop = FALSE]
      matDT <- matDT[, -which(colnames(matDT) == "angles"), drop = FALSE]
    } else {
      ## convert 'd' and 'matDT' to geographic
      d[, "angles"] <- (pi / 2 - d[, "angles"]) %% (2 * pi)
      matDT[, "angles"] <- pi / 2 -  matDT[, "angles", drop = FALSE] %% (2 * pi)
    }

    if (returnDistances) {
      wh <- na.omit(match("rads", colnames(d)))
      if (length(wh) > 0) matDT <- d[, -wh, drop = FALSE]
    } else if (closest) {
      wh <- na.omit(match(c("rads", "dists"), colnames(d)))
      if (length(wh) > 0) matDT <- d[, -wh, drop = FALSE]
    } else {
      if (equalRadii)
        matDTinterior <- matDT[matDT[, "rads"] < (maxRad - 0.71) &
                                 matDT[, "rads"] > (minRad + 0.71), , drop = FALSE]
      else
        matDTinterior <- matDT[matDT[, "rads"] < (matDT[, "maxRad"] - 0.71) &
                                 matDT[, "rads"] > (matDT[, "minRad"] + 0.71), , drop = FALSE]

      matDT <- rbind(d[, colnames(matDTinterior), drop = FALSE], matDTinterior)
      matDT <- matDT[, -which(colnames(matDT) == "rads"), drop = FALSE]
    }
  } else {
    if (!returnAngles) {
      matDT <- matDT[, -which(colnames(matDT) == "angles"), drop = FALSE]
    }
    matDT <- matDT[, -which(colnames(matDT) == "rads"), drop = FALSE]
  }
  if (!returnIndices) {
    ras <- raster(landscape)
    ras[] <- 0
    if (!allowOverlap) {
      if (!returnDistances) {
        ras[matDT[, "indices"]] <- matDT[, "id"]
      } else {
        ras[matDT[, "indices"]] <- matDT[, "dists"]
      }
    } else {
      matDT <- data.table(matDT, key = "indices")
      if (!returnDistances) {
        matDT <- matDT[, sum(id), by = indices]
      } else {
        matDT <- matDT[, sum(1 / dists), by = indices]
      }
      ras[matDT$indices] <- matDT$V1
    }
    return(ras)
  }
  return(matDT)
}

################################################################################
#' Wrap coordinates or pixels in a torus-like fashion
#'
#' Generally useful for model development purposes.
#'
#' If \code{withHeading} used, then \code{X} must be a \code{SpatialPointsDataFrame}
#' that contains two columns, \code{x1} and \code{y1}, with the immediately
#' previous agent locations.
#'
#' @param X A \code{SpatialPoints*} object, or matrix of coordinates.
#'
#' @param bounds Either a \code{Raster*}, \code{Extent}, or \code{bbox} object
#'               defining bounds to wrap around.
#'
#' @param withHeading logical. If \code{TRUE}, the previous points must be wrapped
#'                    also so that the subsequent heading calculation will work.
#'                    Default \code{FALSE}. See details.
#'
#' @return Object of the same class as \code{X}, but with coordinates updated to
#'         reflect the wrapping.
#'
#' @author Eliot McIntire
#' @export
#' @rdname wrap
#'
#' @examples
#' library(raster)
#' library(quickPlot)
#'
#' xrange <- yrange <- c(-50, 50)
#' hab <- raster(extent(c(xrange, yrange)))
#' hab[] <- 0
#'
#' # initialize agents
#' N <- 10
#'
#' # previous points
#' x1 <- rep(0, N)
#' y1 <- rep(0, N)
#' # initial points
#' starts <- cbind(x = stats::runif(N, xrange[1], xrange[2]),
#'                 y = stats::runif(N, yrange[1], yrange[2]))
#'
#' # create the agent object
#' agent <- SpatialPointsDataFrame(coords = starts, data = data.frame(x1, y1))
#'
#'
#' ln <- rlnorm(N, 1, 0.02) # log normal step length
#' sd <- 30 # could be specified globally in params
#'
#' if (interactive()) {
#'   clearPlot()
#'   Plot(hab, zero.color = "white", axes = "L")
#' }
#' for (i in 1:10) {
#'   agent <- crw(agent = agent, extent = extent(hab), stepLength = ln,
#'                stddev = sd, lonlat = FALSE, torus = TRUE)
#'   if (interactive()) Plot(agent, addTo = "hab", axes = TRUE)
#' }
setGeneric("wrap", function(X, bounds, withHeading) {
  standardGeneric("wrap")
})

#' @export
#' @rdname wrap
setMethod(
  "wrap",
  signature(X = "matrix", bounds = "Extent", withHeading = "missing"),
  definition = function(X, bounds) {
    if (identical(colnames(X), c("x", "y"))) {
      return(cbind(
        x = (X[, "x"] - bounds@xmin) %% (bounds@xmax - bounds@xmin) + bounds@xmin,
        y = (X[, "y"] - bounds@ymin) %% (bounds@ymax - bounds@ymin) + bounds@ymin
      ))
    } else {
      stop("When X is a matrix, it must have 2 columns, x and y,",
           "as from say, coordinates(SpatialPointsObj)")
    }
})

#' @export
#' @rdname wrap
setMethod(
  "wrap",
  signature(X = "SpatialPoints", bounds = "ANY", withHeading = "missing"),
  definition = function(X, bounds) {
    X@coords <- wrap(X@coords, bounds = bounds)
    return(X)
})

#' @export
#' @rdname wrap
setMethod(
  "wrap",
  signature(X = "matrix", bounds = "Raster", withHeading = "missing"),
  definition = function(X, bounds) {
    X <- wrap(X, bounds = extent(bounds))
    return(X)
})

#' @export
#' @rdname wrap
setMethod(
  "wrap",
  signature(X = "matrix", bounds = "Raster", withHeading = "missing"),
  definition = function(X, bounds) {
    X <- wrap(X, bounds = extent(bounds))
    return(X)
})

#' @export
#' @rdname wrap
setMethod(
  "wrap",
  signature(X = "matrix", bounds = "matrix", withHeading = "missing"),
  definition = function(X, bounds) {
    if (identical(colnames(bounds), c("min", "max")) &
         (identical(rownames(bounds), c("s1", "s2")))) {
      X <- wrap(X, bounds = extent(bounds))
      return(X)
    } else {
      stop("Must use either a bbox, Raster*, or Extent for 'bounds'")
    }
})

#' @export
#' @rdname wrap
setMethod(
  "wrap",
  signature(X = "SpatialPointsDataFrame", bounds = "Extent", withHeading = "logical"),
  definition = function(X, bounds, withHeading) {
    if (withHeading) {
      # This requires that previous points be "moved" as if they are
      #  off the bounds, so that the heading is correct
      X@data[coordinates(X)[, "x"] < bounds@xmin, "x1"] <-
        (X@data[coordinates(X)[, "x"] < bounds@xmin, "x1"] - bounds@xmin) %%
        (bounds@xmax - bounds@xmin) + bounds@xmax
      X@data[coordinates(X)[, "x"] > bounds@xmax, "x1"] <-
        (X@data[coordinates(X)[, "x"] > bounds@xmax, "x1"] - bounds@xmax) %%
        (bounds@xmin - bounds@xmax) + bounds@xmin
      X@data[coordinates(X)[, "y"] < bounds@ymin, "y1"] <-
        (X@data[coordinates(X)[, "y"] < bounds@ymin, "y1"] - bounds@ymin) %%
        (bounds@ymax - bounds@ymin) + bounds@ymax
      X@data[coordinates(X)[, "y"] > bounds@ymax, "y1"] <-
        (X@data[coordinates(X)[, "y"] > bounds@ymax, "y1"] - bounds@ymax) %%
        (bounds@ymin - bounds@ymax) + bounds@ymin
    }
    return(wrap(X, bounds = bounds))
})

#' @export
#' @rdname wrap
setMethod(
  "wrap",
  signature(X = "SpatialPointsDataFrame", bounds = "Raster", withHeading = "logical"),
  definition = function(X, bounds, withHeading) {
      X <- wrap(X, bounds = extent(bounds), withHeading = withHeading)
      return(X)
})

#' @export
#' @rdname wrap
setMethod(
  "wrap",
  signature(X = "SpatialPointsDataFrame", bounds = "matrix", withHeading = "logical"),
  definition = function(X, bounds, withHeading) {
    if (identical(colnames(bounds), c("min", "max")) &
        identical(rownames(bounds), c("s1", "s2"))) {
      X <- wrap(X, bounds = extent(bounds), withHeading = withHeading)
      return(X)
    } else {
      stop("Must use either a bbox, Raster*, or Extent for 'bounds'")
    }
})

################################################################################
#' Identify outward radiating spokes from initial points
#'
#' This is a generalized version of a notion of a viewshed.
#' The main difference is that there can be many "viewpoints".
#'
#' @inheritParams cir
#'
#' @param stopRule A function. If the spokes are to stop. This can be a function
#'                 of \code{landscape}, \code{fromCell}, \code{toCell}, \code{x}
#'                 (distance from coords cell), or any other named argument passed
#'                 into the \code{...} of this function. See examples.
#'
#' @param nAngles Numeric, length one. Alternative to angles. If provided, the function
#'                will create a sequence of angles from \code{0} to \code{2*pi},
#'                with a length \code{nAngles}, and not including \code{2*pi}.
#'                Will not be used if \code{angles} is provided, and will show
#'                warning of both are given.
#'
#' @param ... Objects to be used by \code{stopRule()}. See examples.
#'
#' @return A matrix containing columns id (representing the row numbers of \code{coords}),
#' angles (from \code{coords} to each point along the spokes), x and y coordinates
#' of each point along the spokes, the corresponding indices on the \code{landscape}
#' Raster, dists (the distances between each \code{coords} and each point along the
#' spokes), and stop, indicating if it was a point that caused a spoke to stop
#' going outwards due to \code{stopRule}.
#'
#' @author Eliot McIntire
#' @export
#' @importFrom fpCompare %<<%
#' @rdname spokes
#'
#' @example inst/examples/example_spokes.R
#'
setGeneric(
  "spokes",
  function(landscape, coords, loci, maxRadius = ncol(landscape) / 4,
           minRadius = maxRadius, allowOverlap = TRUE, stopRule = NULL,
           includeBehavior = "includePixels", returnDistances = FALSE,
           angles = NA_real_, nAngles = NA_real_, returnAngles = FALSE,
           returnIndices = TRUE, ...) {
  standardGeneric("spokes")
})

#' @export
#' @rdname spokes
setMethod(
  "spokes",
  signature(landscape = "RasterLayer", coords = "SpatialPoints", loci = "missing"),
  definition = function(landscape, coords, loci, maxRadius, minRadius = maxRadius,
                        allowOverlap, stopRule, includeBehavior, returnDistances,
                        angles, nAngles, returnAngles, returnIndices, ...) {
  if (!missing(nAngles)) {
    if (missing(angles)) {
    angles <- seq(0, pi * 2, length.out = 17)
    angles <- angles[-length(angles)]
    } else {
      warning("Both angles and nAngles are provided. Using angles only.")
    }
  }

  aCir <- cir(landscape, coords = coords, minRadius = minRadius, maxRadius = maxRadius,
              returnAngles = TRUE, returnDistances = TRUE,
              allowOverlap = allowOverlap, allowDuplicates = TRUE,
              angles = angles, returnIndices = returnIndices)

  if (!is.null(stopRule)) {
    forms <- names(formals(stopRule))
    fromC <- "fromCell" %in% forms
    if (fromC) fromCell <- cellFromXY(landscape, coordinates(coords))
    toC <- "toCell" %in% forms
    if (toC) toCell <- cellFromXY(landscape, to[, c("x", "y")])
    land <- "landscape" %in% forms
    listArgs <- if (land) list(landscape = landscape[aCir[, "indices"]]) else NULL
    if (length(list(...)) > 0) listArgs <- append(listArgs, list(...))
    xDist <- "x" %in% forms

    a <- cbind(aCir, stop = do.call(stopRule, args = listArgs))
    a <- cbind(a, stopDist = a[, "stop"] * a[, "dists"])
    a[a[, "stop"] %==% 0, "stopDist"] <- maxRadius

    sortedUniqAngles <- sort(unique(a[, "angles"]))
    dxx <- lapply(sort(unique(a[, "id"])), function(id) {
      aID <- a[a[, "id"] == id, , drop = FALSE]
      b <- tapply(aID[, "stopDist"], aID[, "angles"], min, na.rm = TRUE)
      d1 <- lapply(sortedUniqAngles, function(x) {
        a1 <- aID[aID[, "angles"] %==% x, , drop = FALSE]
        if (includeBehavior == "excludePixels")
          a1[a1[, "dists"] %<<% b[as.numeric(names(b)) %==% x], , drop = FALSE]
        else
          a1[a1[, "dists"] %<=% b[as.numeric(names(b)) %==% x], , drop = FALSE]
      })
      do.call(rbind, d1)
    })
    d2xx <- do.call(rbind, dxx)
    whDrop <- match(c("stopDist"), colnames(d2xx))
    d2xx[, -whDrop, drop = FALSE]
  }
})

#' This is a very fast version of cir with allowOverlap = TRUE, allowDuplicates = FALSE,
#' returnIndices = TRUE, returnDistances = TRUE, and includeBehavior = "excludePixels".
#' It is used inside spread2, when asymmetry is active. The basic algorithm is to run cir
#' just once, then add to the xy coordinates of every locus
#' @name cirSpecialQuick
#' @inheritParams cir
.cirSpecialQuick <- function(landscape,
                            loci,
                            maxRadius,
                            minRadius) {

  bb <- xyFromCell(landscape, loci)
  middleCell <- if (identical(ncell(landscape) / 2, floor(ncell(landscape) / 2))) {
    ncell(landscape) / 2 - ncol(landscape) / 2
  } else {
    round(ncell(landscape) / 2)
  }
  xy <- xyFromCell(landscape, middleCell)

  # just run one, central locus with cir
  pureCircle2 <- cir(landscape,
                     #loci = attributes(dt)$spreadState$clusterDT$initialPixels,
                     allowOverlap = TRUE, allowDuplicates = FALSE,
                     maxRadius = maxRadius,
                     minRadius = minRadius,
                     returnIndices = TRUE,
                     returnDistances = TRUE,
                     includeBehavior = "excludePixels")
  pureCircle2 <- pureCircle2[order(pureCircle2[, "indices"]), ]
  cc <- cbind(pureCircle2[, "x"] - xy[,"x"], pureCircle2[, "y"] - xy[, "y"])
  dd <- cbind(x = rep(bb[, "x"], each = NROW(pureCircle2)),
              y = rep(bb[, "y"], each = NROW(pureCircle2))) +
    matrix(rep(t(cc), NROW(bb)), ncol = 2, byrow = TRUE)
  lociAll <- rep(loci, each = NROW(pureCircle2))
  distsAll <- rep(pureCircle2[, "dists"], nrow(bb))
  dd <- cbind(id = lociAll, dd, indices = cellFromXY(landscape, dd[, c("x", "y")]),
              dists = distsAll)

  dd[!as.logical(dd[, "x"] > xmax(landscape) | dd[, "x"] < xmin(landscape) |
                       dd[, "y"] > ymax(landscape) | dd[, "y"] < ymin(landscape)), ]

}