R/tess.R

In spatstat.geom: Geometrical Functionality of the 'spatstat' Family

#
# tess.R
#
# support for tessellations
#
#   $Revision: 1.129 $ $Date: 2026/01/17 04:23:46 $
#
tess <- function(..., xgrid=NULL, ygrid=NULL, tiles=NULL, image=NULL,
                 window=NULL, marks=NULL, keepempty=FALSE,
                 unitname=NULL, check=TRUE) {
  uname <- unitname
  if(!is.null(window)) {
    window <- as.owin(window)
    if(is.null(uname)) uname <- unitname(window) 
  }
  isrect <- !is.null(xgrid) && !is.null(ygrid)
  istiled <- !is.null(tiles)
  isimage <- !is.null(image)
  if(isrect + istiled + isimage != 1)
    stop("Must specify either (xgrid, ygrid) or tiles or img")
  if(isrect) {
    stopifnot(is.numeric(xgrid) && all(diff(xgrid) > 0))
    stopifnot(is.numeric(ygrid) && all(diff(ygrid) > 0))
    if(!is.null(window))
      warning("Argument 'window' ignored, because xgrid, grid are given")
    window <- owinInternalRect(range(xgrid), range(ygrid), unitname=uname)
    ntiles <- (length(xgrid)-1) * (length(ygrid)-1)
    out <- list(type="rect", window=window, xgrid=xgrid, ygrid=ygrid, n=ntiles)
  } else if(istiled) {
    stopifnot(is.list(tiles))
    if(check) {
      if(!all(sapply(tiles, is.owin)))
        stop("Tiles must be a list of owin objects")
      if(!is.null(uname)) {
        ## attach new unit name to each tile
        tiles <- lapply(tiles, "unitname<-", value=uname)
      } else {
        ## extract unit names from tiles, check agreement, use as unitname
        uu <- unique(lapply(tiles, unitname))
        uu <- uu[!sapply(uu, is.null)]
        nun <- length(uu)
        if(nun > 1)
          stop("Tiles have inconsistent names for the unit of length")
        if(nun == 1) {
          ## use this unit name
          uname <- uu[[1]]
          if(!is.null(window))
            unitname(window) <- uname
        }
      }
    }
    if(!keepempty && check) {
      # remove empty tiles
      isempty <- sapply(tiles, is.empty)
      if(all(isempty))
        stop("All tiles are empty")
      if(any(isempty))
        tiles <- tiles[!isempty]
    }
    ntiles <- length(tiles)
    nam <- names(tiles)
    lev <- if(!is.null(nam) && all(nzchar(nam))) nam else 1:ntiles
    if(is.null(window)) 
      window <- do.call(union.owin, unname(tiles))
    if(is.mask(window) || any(sapply(tiles, is.mask))) {
      # convert to pixel image tessellation
      Grid <- do.call(commonGrid, append(list(window), unname(tiles)))
      ima <- as.im(window, W=Grid)
      ima$v[] <- NA
      for(i in 1:ntiles)
        ima[tiles[[i]]] <- i
      ima <- ima[window, drop=FALSE]
      ima <- eval.im(factor(ima, levels=1:ntiles))
      levels(ima) <- lev
      out <- list(type="image",
                  window=window, image=ima, n=length(lev))
    } else {
      # tile list
      window <- rescue.rectangle(window)
      out <- list(type="tiled", window=window, tiles=tiles, n=length(tiles))
    }
  } else if(isimage) {
    # convert to factor valued image
    image <- as.im(image)
    if(!is.null(uname)) unitname(image) <- uname
    switch(image$type,
           logical={
             # convert to factor
             if(keepempty) 
               image <- eval.im(factor(image, levels=c(FALSE,TRUE)))
             else
               image <- eval.im(factor(image))
           },
           factor={
             # eradicate unused levels
             if(!keepempty) 
               image <- eval.im(factor(image))
           },
           {
             # convert to factor
             image <- eval.im(factor(image))
           })
               
    if(is.null(window)) window <- as.owin(image)
    out <- list(type="image", window=window, image=image, n=length(levels(image)))
  } else stop("Internal error: unrecognised format")
  ## add marks!
  if(!is.null(marks)) {
    mf <- markformat(marks)
    switch(mf,
           none = { marks <- NULL },
           list = { marks <- hyperframe(marks=marks, row.names=NULL) },
           vector = { marks <- data.frame(marks=marks, row.names=NULL) },
           dataframe = ,
           hyperframe = { row.names(marks) <- NULL }
           )
    if(nrow(marks) != out$n)
      stop(paste("wrong number of marks:",
                 nrow(marks), "should be", out$n),
           call.=FALSE)
    out$marks <- marks
  }
  class(out) <- c("tess", class(out))
  return(out)
}

is.tess <- function(x) { inherits(x, "tess") }

print.tess <- function(x, ..., brief=FALSE) {
  full <- !brief
  if(full) cat("Tessellation\n")
  win <- x$window
  switch(x$type,
         rect={
           if(full) {
             unitinfo <- summary(unitname(win))
             if(evenly.spaced(x$xgrid) && evenly.spaced(x$ygrid)) 
               splat("Tiles are equal rectangles, of dimension",
                     signif(mean(diff(x$xgrid)), 5),
                     "x",
                     signif(mean(diff(x$ygrid)), 5),
                     unitinfo$plural, " ", unitinfo$explain)
             else
               splat("Tiles are unequal rectangles")
           }
           splat(length(x$xgrid)-1, "by", length(x$ygrid)-1, "grid of tiles")
         },
         tiled={
           if(full) {
             if(win$type == "polygonal")
               splat("Tiles are irregular polygons")
             else
               splat("Tiles are windows of general type")
           }
           splat(length(x$tiles), "tiles (irregular windows)")
         },
         image={
           nlev <- length(levels(x$image))
           if(full) {
             splat("Tessellation is determined by a factor-valued image with",
                   nlev, "levels")
           } else splat(nlev, "tiles (levels of a pixel image)")
         })
  if(!is.null(marx <- x$marks)) {
    mf <- markformat(marx)
    switch(mf,
           none = { },
           vector = {
             splat("Tessellation has", paste0(typeof(marx), "-valued marks"))
           },
           list = {
             if(is.solist(marx)) {
               splat("Tessellation has a list of spatial objects as marks")
             } else {
               cls <- unique(sapply(marks, class))
               if(!is.character(cls)) {
                 splat("Tessellation has a list of marks")
               } else {
                 splat("Tessellation has a list of marks of class",
                       commasep(sQuote(cls)))
               }
             }
           },
           dataframe = {
             splat("Tessellation has a data frame of marks:")
             nc <- ncol(marx)
             cn <- colnames(marx)
             ty <- unname(sapply(marx, typeof))
             for(i in seq_len(nc)) {
               cat(paste0("\t$", cn[i], ":\t\t", ty[i], "\n"))
             }
           },
           hyperframe = {
             splat("Tessellation has a hyperframe of marks:")
             nc <- ncol(marx)
             cn <- colnames(marx)
             cls <- sQuote(unclass(marx)$vclass)
             for(i in seq_len(nc)) {
               cat(paste0("\t$", cn[i], ":\t\t", cls[i], "\n"))
             }
           })
  }
  if(full) print(win)
  invisible(NULL)
}

unitname.tess <- function(x) unitname(x$window)

"unitname<-.tess" <- function(x, value) {
  unitname(x$window) <- value
  switch(x$type,
         rect={},
         tiled={
           x$tiles <- lapply(x$tiles, "unitname<-", value)
         },
         image={
           unitname(x$image) <- value
         })
  return(x)
}

plot.tess <- local({

  plotpars <- c("sub", "lty", "lwd",
                "cex.main", "col.main", "font.main",
                "cex.sub", "col.sub", "font.sub", "border")

  plot.tess <- function(x, ..., main, add=FALSE, show.all=!add,
                        border=NULL,
                        do.plot=TRUE,
                        do.labels=!missing(labels), labels=tilenames(x),
                        labelargs=list(),
                        do.col=!missing(values), 
                        values=marks(x),
                        multiplot=TRUE,
                        col=NULL,
                        ribargs=list()) {
    if(missing(main) || is.null(main))
      main <- short.deparse(substitute(x))
    ntiles <- x$n
    if(!do.col) {
      #' Plot tiles, with adornment
      y <- NULL
      result <- NULL
      bbox <- NULL
      need.legend <- FALSE
    } else {
      #' Fill tiles with colours determined by 'values'
      if(markformat(values) == "hyperframe") 
        values <- as.data.frame(values) #' automatic warning
      #' Determine values associated with each tile
      switch(markformat(values),
             none = {
               #' no values assigned.
               #' default is tile name
               tn <- tilenames(x)
               values <- factor(tn, levels=tn)
             },
             vector = {
               #' vector of values.
               #' validate length of vector
               check.anyvector(values, ntiles, things="tiles", vname="values")
             },
             dataframe = {
               #' data frame or matrix of values.
               values <- as.data.frame(values)
               if(nrow(values) != ntiles)
                 stop(paste("Number of rows of values =", nrow(values),
                            "!=", ntiles, "= number of tiles"),
                      call.=FALSE)
               if(multiplot && ncol(values) > 1 && !add) {
                 #' Multiple Panel Plot
                 result <- multi.plot.tess(x, ...,
                                           main=main, show.all=show.all,
                                           border=border, do.plot=do.plot,
                                           do.labels=do.labels, labels=labels,
                                           labelargs=labelargs, do.col=do.col, 
                                           col=col, ribargs=ribargs)
                 return(invisible(result))
               }
               if(ncol(values) > 1)
                 warning("Using only the first column of values")
               values <- values[,1]
             },
             stop("Format of values is not understood")
             )
      #' Single Panel Plot
      #' Determine colour map and plan layout (including colour ribbon)
      #' using rules for pixel images
      y <- as.im(as.function(x, values=values))
      dont.complain.about(y)
      result <- do.call(plot.im,
                        resolve.defaults(
                          list(x=quote(y),
                               do.plot=FALSE,
                               show.all=show.all, add=add, main=main,
                               col=col, ribargs=ribargs),
                          list(...),
                          list(valuesAreColours=FALSE)
                          ))
      #' exit if not actually plotting
      if(!do.plot) return(invisible(result))
      #' extract info
      colmap <- result
      bbox <- attr(result, "bbox")
      bbox.legend <- attr(result, "bbox.legend")
      need.legend <- !is.null(bbox.legend)
    }
    #'      Start Plot 
    #' initialise plot region if it is determined
    if(do.plot && !is.null(bbox) && !add) {
      plot(bbox, main=" ", type="n")
      add <- TRUE
    }
    switch(x$type,
           rect={
             win <- x$window
             dont.complain.about(win)
             z <- do.call.matched(plot.owin,
                                  resolve.defaults(list(x=quote(win),
                                                        main=main,
                                                        add=add,
                                                        show.all=show.all,
                                                        do.plot=do.plot),
                                                   list(...)),
                                  extrargs=plotpars)
             if(is.null(result))
               result <- z
             if(do.plot) {
               #' actually plot
               if(do.col) {
                 #' fill tiles with colours
                 colours <- colmap(values)
                 til <- tiles(x)
                 for(i in seq_len(x$n))
                   plot(til[[i]], add=TRUE,
                        col=colours[i], border=border,
                        main="", ...)
               } else {
                 #' draw tile boundaries only
                 xg <- x$xgrid
                 yg <- x$ygrid
                 do.call.matched(segments,
                                 resolve.defaults(list(x0=xg, y0=win$yrange[1],
                                                       x1=xg, y1=win$yrange[2]),
                                                  list(col=border),
                                                  list(...),
                                                  .StripNull=TRUE))
                 do.call.matched(segments,
                                 resolve.defaults(list(x0=win$xrange[1], y0=yg,
                                                       x1=win$xrange[2], y1=yg),
                                                  list(col=border),
                                                  list(...),
                                                  .StripNull=TRUE))
               }
             }
           },
           tiled={
             xwin <- x$window
             dont.complain.about(xwin)
             z <- do.call.matched(plot.owin,
                                  resolve.defaults(list(x=quote(xwin),
                                                        main=main,
                                                        add=add,
                                                        show.all=show.all,
                                                        do.plot=do.plot),
                                                   list(...)),
                                  extrargs=plotpars)
             if(is.null(result)) result <- z
             if(do.plot) {
               #' plot each tile
               til <- tiles(x)
               if(!do.col) {
                 #' border only
                 lapply(til, plot.owin, ..., add=TRUE, border=border)
               } else {
                 #' fill with colour
                 colours <- colmap(values)
                 mapply(plot.owin, x=til, col=colours,
                        MoreArgs=list(add=TRUE, main="", border=border, ...))
               } 
             }
           },
           image={
             if(is.null(y)) y <- x$image
             dont.complain.about(y)
             result <-
               do.call(plot,
                       resolve.defaults(list(quote(y), add=add, main=main,
                                             show.all=show.all,
                                             do.plot=do.plot,
                                             col=col, ribargs=ribargs),
                                        list(...),
                                        list(valuesAreColours=FALSE)))
             need.legend <- FALSE
           })
    if(do.plot && do.labels) {
      labels <- paste(as.vector(labels))
      til <- tiles(x)
      incircles <- lapply(til, incircle)
      x0 <- sapply(incircles, getElement, name="x")
      y0 <- sapply(incircles, getElement, name="y")
      do.call.matched(text.default,
                      resolve.defaults(list(x=x0, y = y0),
                                       list(labels=labels),
                                       labelargs),
                      funargs=graphicsPars("text"))
    }
    if(do.plot && need.legend) {
      #' determine position of legend
      xlim <- bbox.legend$xrange
      ylim <- bbox.legend$yrange
      sidecode <- attr(colmap, "side.legend")
      vertical <- sidecode %in% c(2,4)
      do.call(plot.colourmap,
              resolve.defaults(list(x=quote(colmap),
                                    add=TRUE, main="",
                                    xlim=xlim, ylim=ylim,
                                    side=sidecode, vertical=vertical),
                               ribargs,
                               list(...)))
    }
    return(invisible(result))
  }

  multi.plot.tess <- function(x, ..., zlim=NULL, col=NULL, equal.ribbon=FALSE) {
    if(equal.ribbon && is.null(zlim) && !inherits(col, "colourmap"))
      zlim <- range(marks(x))
    if(!is.null(zlim)) {
      result <- plot(unstack(x), ..., zlim=zlim, col=col)
    } else {
      result <- plot(unstack(x), ..., col=col)
    }
    return(invisible(result))
  }
  
  plot.tess
})


"[<-.tess" <- function(x, i, ..., value) {
  switch(x$type,
         rect=,
         tiled={
           til <- tiles(x)
           til[i] <- value
           ok <- !unlist(lapply(til, is.null))
           x <- tess(tiles=til[ok])
         },
         image={
           stop("Cannot assign new values to subsets of a pixel image")
         })
  return(x)
}
  
"[.tess" <- function(x, i, ...) {
  trap.extra.arguments(..., .Context="in [.tess")
  if(missing(i)) return(x)
  if(is.owin(i))
    return(intersect.tess(x, i))
  switch(x$type,
         rect=,
         tiled={
           til <- tiles(x)[i]
           return(tess(tiles=til))
         },
         image={
           img <- x$image
           oldlev <- levels(img)
           newlev <- unique(oldlev[i])
           img <- eval.im(factor(img, levels=newlev))
           return(tess(image=img))
         })
}

tiles <- function(x) {
  stopifnot(is.tess(x))
  switch(x$type,
         rect={
           out <- list()
           xg <- x$xgrid
           yg <- x$ygrid
           nx <- length(xg) - 1
           ny <- length(yg) - 1
           for(j in rev(seq_len(ny))) {
             for(i in seq_len(nx)) {
               winij <- owinInternalRect(xg[c(i,i+1)], yg[c(j,j+1)])
               out <- append(out, list(winij))
             }
           }
         },
         tiled={
           out <- x$tiles
           },
         image={
           out <- list()
           ima <- x$image
           lev <- levels(ima)
           for(i in seq_along(lev))
             out[[i]] <- solutionset(ima == lev[i])
           })
  names(out) <- tilenames(x)
  out <- as.solist(out)
  return(out)
}

tiles.empty <- function(x) {
  stopifnot(is.tess(x))
  switch(x$type,
         rect = {
           nx <- length(x$xgrid) - 1
           ny <- length(x$ygrid) - 1
           ans <- rep(FALSE, nx * ny)
         },
         tiled = {
           ans <- sapply(x$tiles, is.empty)
         },
         image = {
           ans <- (table(x$image[]) == 0)
         })
  return(ans)
}
           
tilenames <- function(x) {
  UseMethod("tilenames")
}

tilenames.tess <- function(x) {
  switch(x$type,
         rect={
           if(!is.null(x$tilenames)) {
             out <- x$tilenames
           } else {
             nx <- length(x$xgrid) - 1
             ny <- length(x$ygrid) - 1
             ij <- expand.grid(1:nx, 1:ny)
             out <- paste0("Tile row ", ij[,2], ", col ", ij[,1])
           }
         },
         tiled={
           out <- names(x$tiles)
           if(sum(nzchar(out)) != x$n)
             out <- paste("Tile", seq_len(x$n))
         },
         image={
           out <- levels(x$image)
         }
         )
  return(as.character(out))
}

"tilenames<-" <- function(x, value) {
  UseMethod("tilenames<-")
}

"tilenames<-.tess" <- function(x, value) {
  if(!is.null(value)) {
    ## validate length
    value <- as.character(value)
    nv <- length(value)
    switch(x$type,
           rect = {
             nx <- length(x$xgrid) - 1
             ny <- length(x$ygrid) - 1
             n <- nx * ny
           },
           tiled = { n <- length(x$tiles) },
           image = { n <- length(levels(x$image)) })
    if(nv != n)
      stop("Replacement value has wrong length",
           paren(paste(nv, "instead of", n)))
  }
  switch(x$type,
         rect={
           x$tilenames <- value
         },
         tiled={
           names(x$tiles) <- value
         },
         image={
           levels(x$image) <- value %orifnull% (1:n)
         }
         )
  return(x)
}

marks.tess <- function(x, ...) {
  stopifnot(is.tess(x))
  return(x$marks)
}

"marks<-.tess" <- function(x, ..., value) {
  stopifnot(is.tess(x))
  if(!is.null(value)) {
    mf <- markformat(value)
    switch(mf,
           none = { value <- NULL },
           list = { value <- hyperframe(marks=value, row.names=NULL) },
           vector = { value <- data.frame(marks=value, row.names=NULL) },
           dataframe = ,
           hyperframe = { row.names(value) <- NULL }
           )
    ntil <- x$n
    if(nrow(value) != ntil)
      stop(paste("replacement value for marks has wrong length:",
                 nrow(value), "should be", ntil),
           call.=FALSE)
  }
  x$marks <- value
  return(x)
}

unmark.tess <- function(X) { marks(X) <- NULL; return(X) }

tile.areas <- function(x) {
  stopifnot(is.tess(x))
  switch(x$type,
         rect={
           xg <- x$xgrid
           yg <- x$ygrid
           a <- outer(rev(diff(yg)), diff(xg), "*")
           a <- as.vector(t(a))
           names(a) <- as.vector(t(tilenames(x)))
         },
         tiled={
           a <- as.numeric(sapply(x$tiles, area))
           names(a) <- tilenames(x)
         },
         image={
           z <- x$image
           a <- as.numeric(table(z$v)) * z$xstep * z$ystep
           names(a) <- tilenames(x)
         })
  return(a)
}

tile.centroids <- function(x) {
  stopifnot(is.tess(x))
  switch(x$type,
         rect={
           xg <- x$xgrid
           yg <- x$ygrid
           nx <- length(xg)
           ny <- length(yg)
           xmid <- (xg[-1] + xg[-nx])/2
           ymid <- (yg[-1] + yg[-ny])/2
           xy <- cbind(x=rep(xmid, ny-1),
                       y=rep(rev(ymid), each=nx-1))
         },
         tiled=,
         image = {
           til <- tiles(x)
           xy <- t(sapply(lapply(til, centroid.owin), unlist))
         })
  z <- as.ppp(xy, W=Frame(x))
  return(z)
}
           
as.im.tess <- function(X, W=NULL, ...,
                       eps=NULL, dimyx=NULL, xy=NULL,
                       rule.eps=c("adjust.eps", "grow.frame", "shrink.frame"),
                       rule.tile=c("sample", "majority"),
                       na.replace=NULL, values=NULL) {
  rule.eps <- match.arg(rule.eps)
  rule.tile <- match.arg(rule.tile)
  ## if W is present, it may have to be converted
  if(!is.null(W)) {
    stopifnot(is.owin(W))
    W <- as.mask(W, eps=eps, dimyx=dimyx, xy=xy, rule.eps=rule.eps)
  }
  ## map tile names/indices to other values?
  if(!is.null(values)) {
    if(!is.atomic(values))
      stop("Argument 'values' should contain numeric, logical or factor values",
           call.=FALSE)
    if(length(values) != nobjects(X))
      stop(paste("length(values) =", length(values), "!=",
                 nobjects(X), "= number of tiles"),
           call.=FALSE)
  }
  if(rule.tile == "majority") {
    ## Requires polygonal calculation (except *)
    switch(X$type,
           tiled = { },
           rect = {
             til <- tiles(X)
             WX <- Window(X)
             X <- tess(tiles=til, window=WX)
           },
           image = {
             ## (*) required only if there is a change of pixel grid
             if(!is.null(eps) || !is.null(dimyx) || !is.null(xy)) {
               til <- solapply(tiles(X), as.polygonal)
               WX <- Window(X)
               X <- tess(tiles=til, window=WX)
             }
           })
  }
  ## Do conversion
  switch(X$type,
         image={
           ## Tessellation is already given as a factor image
           ## Sample onto another grid 
           out <- as.im(X$image, W=W, eps=eps, dimyx=dimyx, xy=xy,
                        rule.eps=rule.eps)
           if(!is.null(values)) {
             ## map factor levels to 'values'
             out <- eval.im(values[as.integer(out)])
           }
         },
         tiled={
           ## Tessellation is a list of polygons
           if(is.null(W))
             W <- as.mask(as.owin(X), eps=eps, dimyx=dimyx, xy=xy,
                          rule.eps=rule.eps)
           til <- X$tiles
           ntil <- length(til)
           nama <- names(til)
           if(is.null(nama) || !all(nzchar(nama)))
             nama <- paste(seq_len(ntil))
           switch(rule.tile,
                  sample = {
                    ## sample at centre of each pixel
                    xy <- list(x=W$xcol, y=W$yrow)
                    for(i in seq_len(ntil)) {
                      indic <- as.mask(til[[i]], xy=xy)
                      tag <- as.im(indic, value=i)
                      if(i == 1) {
                        out <- tag
                        outv <- as.integer(out$v)
                      } else {
                        outv <- pmin.int(outv, tag$v, na.rm=TRUE)
                      }
                    }
                    if(is.null(values)) {
                      ## result is factor image
                      outv <- factor(outv, levels=seq_len(ntil), labels=nama)
                    } else {
                      ## map tiles to 'values'
                      outv <- values[outv]
                    }
                    out <- im(outv, out$xcol, out$yrow, unitname=unitname(W))
                  },
                  majority = {
                    ## find tile occuping largest fraction of pixel
                    A <- solapply(til, pixellate, W=W)
                    J <- im.apply(A, which.max)
                    if(is.null(values)) {
                      ## result is factor image
                      out <- eval.im(factor(J,
                                            levels=seq_len(ntil),
                                            labels=nama))
                    } else {
                      ## map tile indices to 'values'
                      out <- eval.im(values[J])
                    }
                  })
         },
         rect={
           ## discretise window or frame to required resolution
           M <- as.im(W %orifnull% as.rectangle(X),
                      eps=eps, dimyx=dimyx, xy=xy,
                      rule.eps=rule.eps)
           Mxcol <- M$xcol
           Myrow <- M$yrow
           ## extract information about rectangular tessellation
           xg <- X$xgrid
           yg <- X$ygrid
           nrows <- length(yg) - 1L
           ncols <- length(xg) - 1L
           ## find tile containing each pixel of M
           jx <- findInterval(Mxcol, xg, rightmost.closed=TRUE)
           iy <- findInterval(Myrow, yg, rightmost.closed=TRUE)
           MM <- as.matrix(M)
           Jcol <- jx[col(MM)]
           Irow <- nrows - iy[row(MM)] + 1L
           Ktile <- as.integer(Jcol + ncols * (Irow - 1L))
           ## 
           if(is.null(values)) {
             ## result is factor image
             outv <- factor(Ktile, levels=seq_len(nrows * ncols))
           } else {
             ## map tiles to 'values'
             outv <- values[Ktile]
           }
           out <- im(outv, xcol=Mxcol, yrow=Myrow,
                     unitname=unitname(W))
         })
  ## replace NA by other value
  out <- na.handle.im(out, na.replace=na.replace)
  return(out)
}

nobjects.tess <- function(x) {
  switch(x$type,
         image = length(levels(x$image)),
         rect = (length(x$xgrid)-1) * (length(x$ygrid)-1),
         tiled = length(x$tiles))
}
  
tileindex <- function(x, y, Z, close.gaps=TRUE, all.inside=FALSE) {
  stopifnot(is.tess(Z))
  if((missing(y) || is.null(y)) && all(c("x", "y") %in% names(x))) {
    y <- x$y
    x <- x$x
  }
  stopifnot(length(x) == length(y))
  switch(Z$type,
         rect={
           jx <- findInterval(x, Z$xgrid, rightmost.closed=TRUE)
           iy <- findInterval(y, Z$ygrid, rightmost.closed=TRUE)
           nrows <- length(Z$ygrid) - 1
           ncols <- length(Z$xgrid) - 1
           iy[iy < 1 | iy > nrows] <- NA
           jx[jx < 1 | jx > ncols] <- NA
           jcol <- jx
           irow <- nrows - iy + 1
           ktile <- jcol + ncols * (irow - 1)
           m <- factor(ktile, levels=seq_len(nrows*ncols))
           ij <- expand.grid(j=seq_len(ncols),i=seq_len(nrows))
           levels(m) <- paste("Tile row ", ij$i, ", col ", ij$j, sep="")
         },
         tiled={
           n <- length(x)
           todo <- seq_len(n)
           til <- Z$tiles
           nt <- length(til)
           m <- integer(n)
           for(i in 1:nt) {
             ti <- til[[i]]
             hit <- inside.owin(x[todo], y[todo], ti)
             if(any(hit)) {
               m[todo[hit]] <- i
               todo <- todo[!hit]
             }
             if(length(todo) == 0)
               break
           }
           m[m == 0] <- NA
           nama <- names(til)
           lev <- seq_len(nt)
           lab <- if(!is.null(nama) && all(nzchar(nama))) nama else paste("Tile", lev)
           m <- factor(m, levels=lev, labels=lab)
         },
         image={
           Zim   <- Z$image
           m <- lookup.im(Zim, x, y, naok=TRUE)
           if(anyNA(m)) {
             #' look up neighbouring pixels
             isna <- is.na(m)
             rc <- nearest.valid.pixel(x[isna], y[isna], Zim, nsearch=2)
             m[isna] <- Zim$v[cbind(rc$row, rc$col)]
           }
         }
         )
  if((close.gaps || all.inside) && anyNA(m)) {
    bad <- is.na(m)
    if(!all.inside) {
      ## Points lying outside the window itself remain classified as NA.
      ## Reclassify only those points which lie inside the window.
      bad <- inside.owin(x[bad], y[bad], Window(Z))
    }
    if(any(bad)) {
      ## assign bad points to nearest tile
      nbad <- sum(bad)
      xbad <- x[bad]
      ybad <- y[bad]
      til <- tiles(Z)
      for(i in seq_along(til)) {
        til.i <- as.polygonal(til[[i]])
        dtile.i <- distfun(til.i)(x,y)
        if(i == 1) {
          dtile.min <- dtile.i
          closest   <- rep(1L, nbad)
        } else {
          better <- (dtile.i < dtile.min)
          if(any(better)) {
            dtile.min[better] <- dtile.i[better]
            closest[better] <- i
          }
        }
      }
      m[bad] <- levels(m)[closest]
    }
  }
  return(m)
}
  
as.tess <- function(X) {
  UseMethod("as.tess")
}

as.tess.tess <- function(X) {
  fields <- 
    switch(X$type,
           rect={ c("xgrid", "ygrid") },
           tiled={ "tiles" },
           image={ "image" },
           stop(paste("Unrecognised tessellation type", sQuote(X$type))))
  fields <- c(c("type", "window", "n", "marks"), fields)
  X <- unclass(X)[fields]
  class(X) <- unique(c("tess", class(X)))
  return(X)
}

as.tess.im <- function(X) {
  return(tess(image = X))
}

as.tess.list <- function(X) {
  W <- lapply(X, as.owin)
  return(tess(tiles=W))
}

as.tess.owin <- function(X) {
  return(tess(tiles=list(X)))
}

domain.tess <- Window.tess <- function(X, ...) { as.owin(X) } 

intersect.tess <- function(X, Y, ..., keepempty=FALSE, keepmarks=FALSE, sep="x") {
  X <- as.tess(X)
  check.1.string(sep)
  keepmarks <- keepmarks && !is.null(marx <- marks(X))
  if(is.owin(Y)) {
    ## efficient code for intersection of a tessellation with a window
    if(Y$type == "mask" || X$type == "image") {
      ## pixel based computation
      if(Y$type == "mask") {
        ## Y is a binary pixel mask
        ## X is any tessellation
        ## Make pixel image by evaluating tile index of X at each pixel of Y
        XimY <- as.im(X, xy=Y)
      } else {
        #' Y is a polygonal or rectangular window
        #' X is a tessellation defined by a pixel image
        XimY <- as.im(X)
      }
      #' Restrict raster to Y
      XimY <- XimY[Y, drop=FALSE]
      ## result is a tessellation defined by a pixel image
      out <- tess(image=XimY, keepempty=keepempty)
      if(keepmarks) {
        if(keepempty) {
          marks(out) <- marx
        } else {
          #' identify non-empty tiles
          tab <- tabulate(XimY[], nbins=nobjects(X))
          marks(out) <- marksubset(marx, tab > 0)
        }
      }
      return(out)
    } else {
      ## polygonal computation
      Zwin <- intersect.owin(as.owin(X), Y)
      Ztiles <- lapply(tiles(X), intersect.owin, B=Y, ..., fatal=FALSE)
      if(!keepempty) {
        isempty <- sapply(Ztiles, is.empty)
        Ztiles <- Ztiles[!isempty]
      }
      out <- tess(tiles=Ztiles, window=Zwin, keepempty=keepempty)
      if(keepmarks) 
        marks(out) <- marksubset(marx, !isempty)
      return(out)
    }
  }
  
  ## general case: intersection of two tessellations
  Y <- as.tess(Y)
  Xtiles <- tiles(X)
  Ytiles <- tiles(Y)
  Ztiles <- list()
  namesX <- tilenames(X)
  namesY <- tilenames(Y)

  if(keepmarks) {
    ## initialise the mark variables to be inherited from parent tessellations
    Xmarks <- marks(X)
    Ymarks <- marks(Y)
    mfX <- markformat(Xmarks)
    mfY <- markformat(Ymarks)
    gotXmarks <- (mfX != "none")
    gotYmarks <- (mfY != "none")
    if(gotXmarks && gotYmarks) {
      ## marks from each input will be combined as separate columns
      switch(mfX,
             vector = { Xmarks <- data.frame(Xmarks=Xmarks) },
             list   = { Xmarks <- hyperframe(Xmarks=Xmarks) },
             hyperframe = ,
             dataframe = {
               colnames(Xmarks) <- paste0("X", colnames(Xmarks))
             })
      switch(mfY,
             vector = { Ymarks <- data.frame(Ymarks=Ymarks) },
             list   = { Ymarks <- hyperframe(Ymarks=Ymarks) },
             hyperframe = ,
             dataframe = {
               colnames(Ymarks) <- paste0("Y", colnames(Ymarks))
             })
      ## ensure hyperframe code is dispatched where required
      if(is.hyperframe(Xmarks) && !is.hyperframe(Ymarks)) 
        Ymarks <- as.hyperframe(Ymarks)
      if(!is.hyperframe(Xmarks) && is.hyperframe(Ymarks)) 
        Xmarks <- as.hyperframe(Xmarks)
    }
    ## initialise 
    if(gotXmarks || gotYmarks) {
      marx <- if(gotXmarks && gotYmarks) {
                cbind(marksubset(Xmarks, integer(0)),
                      marksubset(Ymarks, integer(0)))
              } else if(gotXmarks) {
                marksubset(Xmarks, integer(0))
              } else {
                marksubset(Ymarks, integer(0))
              }
    } else keepmarks <- FALSE
  }

  ## now compute intersection tiles
  Xtrivial <- (length(Xtiles) == 1)
  for(i in seq_along(Xtiles)) {
    Xi <- Xtiles[[i]]
    Ti <- lapply(Ytiles, intersect.owin, B=Xi, ..., fatal=FALSE)
    keep <- keepempty | !sapply(Ti, is.empty)
    if(any(keep)) {
      Ti <- Ti[keep]
      names(Ti) <- if(Xtrivial) namesY[keep] else
                   paste(namesX[i], namesY[keep], sep=sep)
      Ztiles <- append(Ztiles, Ti)
      if(keepmarks) {
        extra <- if(gotXmarks && gotYmarks) {
                   cbind(marksubset(Xmarks, i),
                         marksubset(Ymarks, keep),
                         row.names=NULL)
                 } else if(gotYmarks) {
                   marksubset(Ymarks, keep)
                 } else {
                   marksubset(Xmarks, rep(i, sum(keep)))
                 }
        marx <- rbind(marx, extra)
      }
    }
  }
  ## form tessellation object
  Xwin <- as.owin(X)
  Ywin <- as.owin(Y)
  Zwin <- intersect.owin(Xwin, Ywin)
  out <- tess(tiles=Ztiles, window=Zwin, keepempty=keepempty)
  if(keepmarks) 
    marks(out) <- marx
  return(out)
}

venn.tess <- function(..., window=NULL, labels=FALSE) {
  argh <- list(...)
  nargh <- length(argh)
  if(nargh == 0) stop("No arguments given")
  iswin <- sapply(argh, is.owin)
  istes <- sapply(argh, is.tess)
  if(!all(iswin | istes)) 
    stop("All arguments must be windows or tessellations", call.=FALSE)
  nama <- names(argh)
  if(sum(nzchar(nama)) < nargh)
    nama <- paste0("T", seq_len(nargh))
  W <- window %orifnull% do.call(union.owin, unname(lapply(argh, as.owin)))
  for(i in seq_len(nargh)) {
    A <- argh[[i]]
    if(is.owin(A)) {
      Z <- list(A, Out=setminus.owin(W, A))
      names(Z) <- paste0(c("", "not"), nama[i])
      A <- tess(tiles=Z, window=W)
      if(labels) marks(A) <- c(TRUE, FALSE)
    }
    Y <- if(i == 1) A else intersect.tess(Y, A, keepmarks=labels)
  }
  if(labels) colnames(marks(Y)) <- nama
  return(Y)
}

bdist.tiles <- local({

  vdist <- function(x,w) {
    z <- as.ppp(vertices(x), W=w, check=FALSE)
    min(bdist.points(z))
  }
  edist <- function(x,b) {
    xd <- crossdist(edges(x, check=FALSE), b, type="separation")
    min(xd)
  }

  bdist.tiles <-  function(X) {
    if(!is.tess(X))
      stop("X must be a tessellation")
    W <- as.owin(X)
    switch(X$type,
           rect=,
           tiled={
             tt <- tiles(X)
             if(is.convex(W)) {
               # distance is minimised at a tile vertex
               d <- sapply(tt, vdist, w=W)
             } else {
               # coerce everything to polygons
               W  <- as.polygonal(W)
               tt <- lapply(tt, as.polygonal)
               # compute min dist from tile edges to window edges
               d <- sapply(tt, edist, b=edges(W))
             }
           },
           image={
             Xim <- X$image
             # compute boundary distance for each pixel
             bd <- bdist.pixels(as.owin(Xim), style="image")
             bd <- bd[W, drop=FALSE]
             # split over tiles
             bX <- split(bd, X)
             # compute minimum distance over each level of factor
             d <- sapply(bX, function(z) { summary(z)$min })
           }
           )
    return(d)
  }
  bdist.tiles
})


## ......... geometrical transformations ..................

shift.tess <- function(X, ...) {
  Y <- X
  Y$window <- wY <- shift(X$window, ...)
  vec <- getlastshift(wY)
  switch(X$type,
         rect={
           Y$xgrid <- Y$xgrid + vec[1]
           Y$ygrid <- Y$ygrid + vec[2]
         },
         tiled={
           Y$tiles <- lapply(Y$tiles, shift, vec=vec)
         },
         image = {
           Y$image <- shift(Y$image, vec)
         })
  attr(Y, "lastshift") <- vec
  return(Y)
}

affine.tess <- function(X, mat=diag(c(1,1)), vec=c(0,0), ...) {
  Y <- X
  Y$window <- affine(X$window, mat=mat, vec=vec, ...)
  switch(Y$type,
         rect = {
           if(all(mat == diag(diag(mat)))) {
             ## result is rectangular
             Y$xgrid <- sort(mat[1,1] * X$xgrid + vec[1])
             Y$ygrid <- sort(mat[2,2] * X$ygrid + vec[2])
           } else {
             ## shear transformation; treat rectangles as general tiles
             Y <- tess(tiles=tiles(X), window=Y$window)
             Y$tiles <- lapply(Y$tiles, affine, mat=mat, vec=vec, ...)
           }
         },
         tiled={
           Y$tiles <- lapply(Y$tiles, affine, mat=mat, vec=vec, ...)
         },
         image = {
           Y$image <- affine(Y$image, mat=mat, vec=vec, ...)
         })
  return(Y)
}

reflect.tess <- function(X) {
  Y <- X
  Y$window <- reflect(Y$window)
  switch(X$type,
         rect = {
           Y$xgrid <- rev(- Y$xgrid)
           Y$ygrid <- rev(- Y$ygrid)
         },
         tiled = {
           Y$tiles <- lapply(Y$tiles, reflect)
         },
         image = {
           Y$image <- reflect(Y$image)
         })
  return(Y)
}

flipxy.tess <- function(X) {
  Y <- X
  Y$window <- flipxy(Y$window)
  switch(X$type,
         rect = {
           Y$xgrid <- X$ygrid
           Y$ygrid <- X$xgrid
         },
         tiled = {
           Y$tiles <- lapply(Y$tiles, flipxy)
         },
         image = {
           Y$image <- flipxy(Y$image)
         })
  return(Y)
}

scalardilate.tess <- function(X, f, ...) {
  Y <- X
  Y$window <- scalardilate(X$window, f, ...)
  switch(X$type,
         rect = {
           Y$xgrid <- f * Y$xgrid
           Y$ygrid <- f * Y$ygrid
         },
         tiled = {
           Y$tiles <- lapply(Y$tiles, scalardilate, f=f, ...)
         },
         image = {
           Y$image <- scalardilate(Y$image, f=f, ...)
         })
  return(Y)
}

rotate.tess <- function(X, angle=pi/2, ..., centre=NULL) {
  if(angle %% (2 * pi) == 0) return(X)
  if(!is.null(centre)) {
    X <- shift(X, origin=centre)
    negorigin <- getlastshift(X)
  } else negorigin <- NULL
  Y <- X
  Y$window <- rotate(X$window, angle=angle, ...)
  switch(X$type,
         rect = {
           if(angle %% (pi/2) == 0) {
             ## result is rectangular
             co <- round(cos(angle))
             si <- round(sin(angle))
             Y$xgrid <- sort((if(co == 0) 0 else (co * X$xgrid)) -
                             (if(si == 0) 0 else (si * X$ygrid)))
             Y$ygrid <- sort((if(si == 0) 0 else (si * X$xgrid)) +
                             (if(co == 0) 0 else (co * X$ygrid)))
           } else {
             ## general tessellation
             Y <- tess(tiles=lapply(tiles(X), rotate, angle=angle, ...),
                       window=Y$window)
           }
         },
         tiled = {
           Y$tiles <- lapply(X$tiles, rotate, angle=angle, ...)
         },
         image = {
           Y$image <- rotate(X$image, angle=angle, ...)
         })
  if(!is.null(negorigin))
    Y <- shift(Y, -negorigin)
  return(Y)
}
  
as.data.frame.tess <- function(x, ...) {
  switch(x$type,
         rect =,
         tiled = {
           y <- lapply(tiles(x), as.data.frame, ...)
           z <- mapply(assignDFcolumn,
                       x=y, value=tilenames(x),
                       MoreArgs=list(name="Tile", ...),
                       SIMPLIFY=FALSE)
           z <- do.call(rbind, z)
           row.names(z) <- NULL
         },
         image = {
           z <- as.data.frame(x$image, ...)
           if(!is.na(m <- match("value", colnames(z))))
             colnames(z)[m] <- "Tile"
         },
         {
           z <- NULL
           warning("Unrecognised type of tessellation")
         })
  return(z)
}

connected.tess <- function(X, ...) {
  Xim <- as.im(X, ...)
  X <- as.tess(Xim)
  tilesX <- tiles(X)
  namesX <- names(tilesX)
  shards <- lapply(tilesX, connected) # list of factor images
  shardnames <- lapply(shards, levels)
  nshards <- lengths(shardnames)
  broken <- (nshards > 1)
  #' unbroken tiles keep their original tile names
  shardnames[!broken] <- namesX[!broken]
  #' shards of broken tiles are named "tilename[i] shard j"
  shardnames[broken] <- mapply(paste,
                               namesX[broken], "shard", shardnames[broken],
                               SIMPLIFY=FALSE)
  #' rename them
  shards <- mapply("levels<-", shards, shardnames, SIMPLIFY=FALSE)
  #' separate them
  shards <- lapply(lapply(shards, as.tess), tiles)
  shards <- unlist(shards, recursive=FALSE, use.names=FALSE)
  names(shards) <- unlist(shardnames)
  #' form tessellation
  result <- tess(tiles=shards, window=as.owin(Xim))
  result
}

identify.tess <- function(x, ..., labels=tilenames(x),
                          n=nobjects(x), plot=TRUE,
                          paint=plot, paint.args=list()) {
  verifyclass(x, "tess")
  if (dev.cur() == 1 && interactive()) {
    eval(substitute(plot(X), list(X = substitute(x))))
  }
  if(!missing(labels)) {
    if(!is.null(dim(labels))) labels <- labels[,1]
    labels <- unname(unlist(labels))
    if(length(labels) != nobjects(x))
      stop(paste("The number of", sQuote("labels"),
                 "should equal the number of tiles"),
           call.=FALSE)
  }
  ## Find a representative point for each tile for plotting labels
  til <- tiles(x)
  incircles <- lapply(til, incircle)
  xc <- sapply(incircles, getElement, name="x")
  yc <- sapply(incircles, getElement, name="y")
  gpo <- graphicsPars("owin")
  gpt <- graphicsPars("text")
  ## start loop
  id <- integer(0)
  while(length(id) < n) {
    mouse <- spatstatLocator(1, type="n")
    ## check for interrupt exit
    if(length(mouse$x) == 0)
      break
    ## determine tile
    ident <- as.integer(tileindex(mouse$x, mouse$y, x))
    if(length(ident) == 0) {
      cat("Query location is too far away\n")
    } else if(ident %in% id) {
      cat(paste("Tile", ident, "already selected\n"))
    } else {
      ## add to list
      if(plot || paint) {
        ## Plot label
        mix <- xc[ident]
        miy <- yc[ident]
        li <- labels[ident]
        dont.complain.about(li, mix, miy)
        tili <- til[[ident]]
        if(paint) {
          if(is.rectangle(tili) || is.polygonal(tili)) {
            do.call.matched(plot.owin,
                            resolve.defaults(list(x=tili, add=TRUE),
                                             paint.args,
                                             list(...),
                                             list(col="#AAAAE6"),
                                             list(border=1, lwd=2)),
                            extrargs=gpo)
          } else {
            do.call.matched(plot.owin,
                            resolve.defaults(list(x=tili, add=TRUE),
                                             paint.args,
                                             list(...),
                                             list(col="#AAAAE6")),
                            extrargs=gpo)
            plot(edges(tili), add=TRUE, col=1, lwd=2)
          }
        }
        if(plot) 
          do.call.matched(graphics::text.default,
                          resolve.defaults(list(x=quote(mix), 
                                                y=quote(miy), 
                                                labels=quote(li)),
                                           list(...)),
                          extrargs=gpt)
      }
      id <- c(id, ident)
    }
  }
  marx <- marks(x)
  out <- data.frame(id=id, name=tilenames(x)[id])
  switch(markformat(marx),
         vector = {
           out <- cbind(out, data.frame(marks=marx[id]))
         },
         dataframe = {
           out <- cbind(out, marx[id, , drop=FALSE])
         },
         {})
  return(out)
}