#
# tess.R
#
# support for tessellations
#
#   $Revision: 1.42 $ $Date: 2011/05/18 09:20:08 $
#
tess <- function(..., xgrid=NULL, ygrid=NULL, tiles=NULL, image=NULL,
                 window=NULL, keepempty=FALSE) {
  if(!is.null(window))
    win <- as.owin(window)
  else win <- NULL
  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(win)) win <- owin(range(xgrid), range(ygrid))
    ntiles <- (length(xgrid)-1) * (length(ygrid)-1)
    out <- list(type="rect", window=win, xgrid=xgrid, ygrid=ygrid, n=ntiles)
  } else if(istiled) {
    stopifnot(is.list(tiles))
    if(!all(unlist(lapply(tiles, is.owin))))
      stop("tiles must be a list of owin objects")
    if(!keepempty) {
      # remove empty tiles
      isempty <- unlist(lapply(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(win)) {
      for(i in 1:ntiles) {
        if(i == 1)
          win <- tiles[[1]]
        else
          win <- union.owin(win, tiles[[i]])
      }
    }
    ismask <- function(x) {x$type == "mask"}
    if(ismask(win) || any(unlist(lapply(tiles, ismask)))) {
      # convert to pixel image tessellation
      win <- as.mask(win)
      ima <- as.im(win)
      ima$v[] <- NA
      for(i in 1:ntiles)
        ima[tiles[[i]]] <- i
      ima <- ima[win, drop=FALSE]
      ima <- eval.im(factor(ima, levels=1:ntiles))
      levels(ima) <- lev
      out <- list(type="image",
                  window=win, image=ima, n=length(lev))
    } else {
      # tile list
      win <- rescue.rectangle(win)
      out <- list(type="tiled", window=win, tiles=tiles, n=length(tiles))
    }
  } else if(isimage) {
    # convert to factor valued image
    image <- as.im(image)
    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(win)) win <- as.owin(image)
    out <- list(type="image", window=win, image=image, n=length(levels(image)))
  } else stop("Internal error: unrecognised format")
  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))
             equispaced <- function(z) {
               dz <- diff(z)
               diff(range(dz))/mean(dz) < 0.01
             }
             if(equispaced(x$xgrid) && equispaced(x$ygrid)) 
               cat(paste("Tiles are equal rectangles, of dimension",
                         signif(mean(diff(x$xgrid)), 5),
                         "x",
                         signif(mean(diff(x$ygrid)), 5),
                         unitinfo$plural, " ", unitinfo$explain,
                         "\n"))
             else
               cat(paste("Tiles are unequal rectangles\n"))
           }
           cat(paste(length(x$xgrid)-1, "by", length(x$ygrid)-1,
                     "grid of tiles", "\n"))
         },
         tiled={
           if(full) {
             if(win$type == "polygonal")
               cat("Tiles are irregular polygons\n")
             else
               cat("Tiles are windows of general type\n")
           }
           cat(paste(length(x$tiles), "tiles (irregular windows)\n"))
         },
         image={
           nlev <- length(levels(x$image))
           if(full) {
             cat(paste("Tessellation is determined by",
                       "a factor-valued image",
                       "with", nlev, "levels\n"))
           } else cat(paste(nlev, "tiles (levels of a pixel image)\n"))
         })
  if(full) print(win)
  invisible(NULL)
}

plot.tess <- function(x, ..., main, add=FALSE, col=NULL) {
  xname <- deparse(substitute(x))
  if(missing(main))
    main <- xname
  switch(x$type,
         rect={
           win <- x$window
           if(!add)
             do.call.matched("plot.owin",
                             resolve.defaults(list(x=win, main=main),
                                              list(...)),
                             extrargs=c("sub", "lty", "lwd"))
           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=col),
                                            list(...),
                                            .StripNull=TRUE))
           do.call.matched("segments",
                           resolve.defaults(list(x0=win$xrange[1], y0=yg,
                                                 x1=win$xrange[2], y1=yg),
                                            list(col=col),
                                            list(...),
                                            .StripNull=TRUE))
         },
         tiled={
           if(!add)
             do.call.matched("plot.owin",
                             resolve.defaults(list(x=x$window, main=main),
                                              list(...)))
           til <- tiles(x)
           plotem <- function(z, ..., col=NULL) {
             if(is.null(col))
               plot(z, ..., add=TRUE)
             else if(z$type != "mask")
               plot(z, ..., border=col, add=TRUE)
             else plot(z, ..., col=col, add=TRUE)
           }
           lapply(til, plotem, ..., col=col)
         },
         image={
           plot(x$image, main=main, ..., add=add)
         })
  return(invisible(NULL))
}

"[<-.tess" <- function(x, ..., value) {
  switch(x$type,
         rect=,
         tiled={
           til <- tiles(x)
           til[...] <- 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, ...) {
  switch(x$type,
         rect=,
         tiled={
           til <- tiles(x)[...]
           return(tess(tiles=til))
         },
         image={
           img <- x$image
           oldlev <- levels(img)
           newlev <- unique(oldlev[...])
           img <- eval.im(factor(ifelse(img %in% newlev, img, NA)))
           levels(img) <- newlev
           return(tess(image=img))
         })
}

tiles <- function(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 <- owin(xg[c(i,i+1)], yg[c(j,j+1)])
               dout <- list(winij)
               names(dout) <- paste("Tile row ", ny-j+1, ", col ", i,
                                    sep="")
               out <- append(out, dout)
             }
         },
         tiled={
           out <- x$tiles
           if(is.null(names(out)))
             names(out) <- paste("Tile", seq_along(out))
         },
         image={
           out <- list()
           ima <- x$image
           lev <- levels(ima)
           for(i in seq_along(lev))
             out[[i]] <- solutionset(ima == lev[i])
           names(out) <- paste(lev)
         })
  out <- as.listof(out)
  return(out)
}

tilenames <- function(x) {
  stopifnot(is.tess(x))
  switch(x$type,
         rect={
           nx <- length(x$xgrid) - 1
           ny <- length(x$ygrid) - 1
           nam <- outer(rev(seq_len(ny)),
                        seq_len(nx),
                        function(j,i,ny) {
                          paste("Tile row ", ny-j+1, ", col ", i,
                                sep="")},
                        ny=ny)
           return(nam)
         },
         tiled={
           til <- x$tiles
           if(!is.null(names(til)))
             nam <- names(til)
           else 
             nam <- paste("Tile", seq_along(til))
         },
         image={
           ima <- x$image
           lev <- levels(ima)
           nam <- paste(lev)
         })
  return(nam)
}

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

         
as.im.tess <- function(X, W=NULL, ...,
                       eps=NULL, dimyx=NULL, xy=NULL,
                       na.replace=NULL) {
  # if W is present, it may have to be converted
  if(!is.null(W)) {
    stopifnot(is.owin(W))
    if(W$type != "mask")
      W <- as.mask(W, eps=eps, dimyx=dimyx, xy=xy)
  } 
  switch(X$type,
         image={
           out <- as.im(X$image, W=W, eps=eps, dimyx=dimyx, xy=xy,
                        na.replace=na.replace)
         },
         tiled={
           if(is.null(W))
             W <- as.mask(as.owin(X), eps=eps, dimyx=dimyx, xy=xy)
           til <- X$tiles
           ntil <- length(til)
           nama <- names(til)
           if(is.null(nama) || !all(nzchar(nama)))
             nama <- paste(seq_len(ntil))
           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 <- out$v
             } else {
               outv <- pmin(outv, tag$v, na.rm=TRUE)
             }
           }
           out <- im(factor(outv, levels=nama),
                     out$xcol, out$yrow)
           unitname(out) <- unitname(W)
         },
         rect={
           if(is.null(W))
             out <- as.im(as.rectangle(X), eps=eps, dimyx=dimyx, xy=xy)
           else
             out <- as.im(W)
           xg <- X$xgrid
           yg <- X$ygrid
           nrows <- length(yg) - 1
           ncols <- length(xg) - 1
           jx <- findInterval(out$xcol, xg, rightmost.closed=TRUE)
           iy <- findInterval(out$yrow, yg, rightmost.closed=TRUE)
           M <- as.matrix(out)
           Jcol <- jx[col(M)]
           Irow <- nrows - iy[row(M)] + 1
           Ktile <- Jcol + ncols * (Irow - 1)
           Ktile <- factor(Ktile, levels=seq_len(nrows * ncols))
           out <- im(Ktile, xcol=out$xcol, yrow=out$yrow,
                     unitname=unitname(W))
         }
         )
  return(out)
}

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"), fields)
  X <- unclass(X)[fields]
  class(X) <- c("tess", class(X))
  return(X)
}

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

as.tess.quadratcount <- function(X) {
  return(attr(X, "tess"))
}

as.tess.quadrattest <- function(X) {
  Y <- attr(X, "quadratcount")
  Z <- attr(Y, "tess")
  return(Z)
}

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

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

intersect.tess <- function(X, Y, ...) {
  X <- as.tess(X)
  if(is.owin(Y) && Y$type == "mask") {
    # special case
    # convert to pixel image 
    result <- as.im(Y)
    Xtiles <- tiles(X)
    for(i in seq_along(Xtiles)) {
      tilei <- Xtiles[[i]]
      result[tilei] <- i
    }
    result <- result[Y, drop=FALSE]
    return(tess(image=result, window=Y))
  }
  Y <- as.tess(Y)
  Xtiles <- tiles(X)
  Ytiles <- tiles(Y)
  Ztiles <- list()
  for(i in seq_along(Xtiles))
    for(j in seq_along(Ytiles)) {
      Tij <- intersect.owin(Xtiles[[i]], Ytiles[[j]], ..., fatal=FALSE)
      if(!is.null(Tij) && !is.empty(Tij))
        Ztiles <- append(Ztiles, list(Tij))
    }
  Xwin <- as.owin(X)
  Ywin <- as.owin(Y)
  Zwin <- intersect.owin(Xwin, Ywin)
  return(tess(tiles=Ztiles, window=Zwin))
}


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
             vdist <- function(x,w) {
               z <- as.ppp(vertices(x), W=w)
               min(bdist.points(z))
             }
             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
             edist <- function(x,b) {
               xd <- crossdist(as.psp(x), b, type="separation")
               min(xd)
             }
             d <- sapply(tt, edist, b=as.psp(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)
}