R/window.R

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

#
#	window.S
#
#	A class 'owin' to define the "observation window"
#
#	$Revision: 4.211 $	$Date: 2024/02/01 02:30:49 $
#
#
#	A window may be either
#
#		- rectangular:
#                       a rectangle in R^2
#                       (with sides parallel to the coordinate axes)
#
#		- polygonal:
#			delineated by 0, 1 or more non-self-intersecting
#                       polygons, possibly including polygonal holes.
#	
#		- digital mask:
#			defined by a binary image
#			whose pixel values are TRUE wherever the pixel
#                       is inside the window
#
#	Any window is an object of class 'owin', 
#       containing at least the following entries:	
#
#		$type:	a string ("rectangle", "polygonal" or "mask")
#
#		$xrange   
#		$yrange
#			vectors of length 2 giving the real dimensions 
#			of the enclosing box
#               $units
#                       name of the unit of length
#
#	The 'rectangle' type has only these entries.
#
#       The 'polygonal' type has an additional entry
#
#               $bdry
#                       a list of polygons.
#                       Each entry bdry[[i]] determines a closed polygon.
#
#                       bdry[[i]] has components $x and $y which are
#                       the cartesian coordinates of the vertices of
#                       the i-th boundary polygon (without repetition of
#                       the first vertex, i.e. same convention as in the
#                       plotting function polygon().)
#
#
#	The 'mask' type has entries
#
#		$m		logical matrix
#		$dim		its dimension array
#		$xstep,ystep	x and y dimensions of a pixel
#		$xcol	        vector of x values for each column
#               $yrow           vector of y values for each row
#	
#	(the row index corresponds to increasing y coordinate; 
#	 the column index "   "     "   "  "  "  x "   "    ".)
#
#
#-----------------------------------------------------------------------------
#

.Spatstat.Image.Warning <-
  c("Row index corresponds to increasing y coordinate; column to increasing x",
    "Transpose matrices to get the standard presentation in R",
    "Example: image(result$xcol,result$yrow,t(result$d))")

owin <- function(xrange=c(0,1), yrange=c(0,1),
                 ..., poly=NULL, mask=NULL, unitname=NULL, xy=NULL) {
  ## trap a common abuse of syntax: owin(window)
  if(nargs() == 1 && !missing(xrange) && is.owin(xrange))
    return(xrange)

  ## parse
  poly.given <- !is.null(poly)
  mask.given <- !is.null(mask)
  range.given <- !missing(xrange) || !missing(yrange)
  
  if(range.given) {
    if(missing(xrange) != missing(yrange))
      stop("If one of xrange, yrange is specified then both must be.")
    xrange <- unname(xrange)
    yrange <- unname(yrange)
  }
  
  if(poly.given && mask.given)
     stop("Ambiguous -- both polygonal boundary and digital mask supplied")

  if(!mask.given && !is.null(xy))
    warning("Argument xy ignored: it is only applicable when a mask is given")

  if(poly.given) {
    ## convert data frames to vanilla lists
    if(is.data.frame(poly))
      poly <- as.list(poly)
    else if(is.list(poly) && any(unlist(lapply(poly, is.data.frame))))
      poly <- lapply(poly, as.list)
  }

  if(!poly.given && !mask.given) {
    ## rectangular window
    owinInternalRect(xrange, yrange, ..., unitname=unitname)
  } else if(poly.given) {
    ## polygonal window
    if(range.given) {
      owinInternalPoly(xrange, yrange, ..., poly=poly, unitname=unitname)
    } else {
      owinInternalPoly(                ..., poly=poly, unitname=unitname)
    }
  } else {
    ## mask window
    if(range.given) {
      owinInternalMask(xrange, yrange, ..., mask=mask, xy=xy, unitname=unitname)
    } else {
      owinInternalMask(                ..., mask=mask, xy=xy, unitname=unitname)
    }
  }
}

owinInternalRect <- function(xrange=c(0,1), yrange=c(0,1),
                             ..., unitname=NULL, check = TRUE) {
  if(check) {
    if(!is.vector(xrange) || length(xrange) != 2 || xrange[2L] < xrange[1L])
      stop("xrange should be a vector of length 2 giving (xmin, xmax)")
    if(!is.vector(yrange) || length(yrange) != 2 || yrange[2L] < yrange[1L])
      stop("yrange should be a vector of length 2 giving (ymin, ymax)")
  }
  unitname <- as.unitname(unitname)
  w <- list(type="rectangle", xrange=xrange, yrange=yrange, units=unitname)
  class(w) <- "owin"
  return(w)
}


isXYdata <- function(x) { (is.matrix(x) || is.data.frame(x)) && ncol(x) == 2 }
asXYdata <- function(xy) { list(x=xy[,1], y=xy[,2]) }

owinInternalPoly <- function(xrange=c(0,1), yrange=c(0,1),
                             ..., poly=NULL, unitname=NULL, 
                             check     = TRUE,
                             calculate = check,
                             strict    = spatstat.options("checkpolygons"),
                             fix       = spatstat.options("fixpolygons")) {
  unitname <- as.unitname(unitname)
  if(length(poly) == 0) {
    ## empty polygon
    if(check) {
      if(!is.vector(xrange) || length(xrange) != 2 || xrange[2L] < xrange[1L])
        stop("xrange should be a vector of length 2 giving (xmin, xmax)")
      if(!is.vector(yrange) || length(yrange) != 2 || yrange[2L] < yrange[1L])
        stop("yrange should be a vector of length 2 giving (ymin, ymax)")
    }
    w <- list(type="polygonal", xrange=xrange, yrange=yrange,
              bdry=list(), units=unitname)
    class(w) <- "owin"
    return(w)
  }
  ## convert matrix or data frame to list(x,y)
  if(isXYdata(poly)) {
    poly <- asXYdata(poly)
  } else if(is.list(poly) && all(unlist(lapply(poly, isXYdata)))) {
    poly <- lapply(poly, asXYdata)
  }
  ## nonempty polygon  
  ## test whether it's a single polygon or multiple polygons
  if(verify.xypolygon(poly, fatal=FALSE))
    psingle <- TRUE
  else if(all(unlist(lapply(poly, verify.xypolygon, fatal=FALSE))))
    psingle <- FALSE
  else
    stop("poly must be either a list(x,y) or a list of list(x,y)")
  
  w.area <- NULL

  if(psingle) {
    ## single boundary polygon
    bdry <- unname(list(poly))
    if(check || calculate) {
      w.area <- Area.xypolygon(poly)
      if(w.area < 0)
        stop(paste("Area of polygon is negative -",
                   "maybe traversed in wrong direction?"))
    }
  } else {
    ## multiple boundary polygons
    bdry <- unname(poly)
    if(check || calculate) {
      w.area <- sapply(poly, Area.xypolygon)
      if(sum(w.area) < 0)
        stop(paste("Area of window is negative;\n",
                   "check that all polygons were traversed",
                   "in the right direction"))
    }
  }

  actual.xrange <- range(unlist(lapply(bdry, getElement, name="x")))
  if(missing(xrange))
    xrange <- actual.xrange
  else if(check) {
    if(!is.vector(xrange) || length(xrange) != 2 || xrange[2L] < xrange[1L])
      stop("xrange should be a vector of length 2 giving (xmin, xmax)")
    if(!all(xrange == range(c(xrange, actual.xrange))))
      stop("polygon's x coordinates outside xrange")
  }
    
  actual.yrange <- range(unlist(lapply(bdry, getElement, name="y")))
  if(missing(yrange))
    yrange <- actual.yrange
  else if(check) {
    if(!is.vector(yrange) || length(yrange) != 2 || yrange[2L] < yrange[1L])
      stop("yrange should be a vector of length 2 giving (ymin, ymax)")
    if(!all(yrange == range(c(yrange, actual.yrange))))
      stop("polygon's y coordinates outside yrange")
  }

  if(!is.null(w.area)) {
    ## tack on area and hole data
    holes <- (w.area < 0)
    for(i in seq_along(bdry)) 
      bdry[[i]] <- append(bdry[[i]], list(area=w.area[i], hole=holes[i]))
  }
    
  w <- list(type="polygonal",
            xrange=xrange, yrange=yrange, bdry=bdry, units=unitname)
  class(w) <- "owin"
      
  if(check && strict) { 
    ## strict checks on geometry (self-intersection etc)
    ok <- owinpolycheck(w)
    if(!ok) {
      errors <- attr(ok, "err")
      stop(paste("Polygon data contain", commasep(errors)))
    }
  }
  if(check && fix) {
    if(length(bdry) == 1 &&
       length(bx <- bdry[[1L]]$x) == 4 &&
       length(unique(bx)) == 2 &&
       length(unique(bdry[[1L]]$y)) == 2) {
      ## it's really a rectangle
      if(Area.xypolygon(bdry[[1L]]) < 0)
        w$bdry <- lapply(bdry, reverse.xypolygon)
    } else {
      ## repair polygon data by invoking polyclip
      ##        to intersect polygon with larger-than-bounding rectangle
      ##        (Streamlined version of intersect.owin)
      ww <- lapply(bdry, reverse.xypolygon)
      xrplus <- mean(xrange) + c(-1,1) * diff(xrange)
      yrplus <- mean(yrange) + c(-1,1) * diff(yrange)
      bignum <- (.Machine$integer.max^2)/2
      epsclip <- max(diff(xrange), diff(yrange))/bignum
      rr <- list(list(x=xrplus[c(1,2,2,1)], y=yrplus[c(2,2,1,1)]))
      bb <- polyclip::polyclip(ww, rr, "intersection",
                               fillA="nonzero", fillB="nonzero", eps=epsclip)
      ## ensure correct polarity
      totarea <- sum(unlist(lapply(bb, Area.xypolygon)))
      if(totarea < 0)
        bb <- lapply(bb, reverse.xypolygon)
      w$bdry <- bb
    }
  }
  return(w)
}

owinInternalMask <- function(xrange=c(0,1), yrange=c(0,1),
                             ..., mask=NULL, unitname=NULL, xy=NULL,
                             check     = TRUE) {

  unitname <- as.unitname(unitname)

  if(is.data.frame(mask) &&
     ncol(mask) %in% c(2,3) &&
     sum(sapply(mask, is.numeric)) == 2) {
    ## data frame with 2 columns of coordinates
    W <- as.owin(W=mask, xy=xy)
    unitname(W) <- unitname
    return(W)
  }
    
  if(!is.matrix(mask))
    stop(paste(sQuote("mask"), "must be a matrix"))
  if(!is.logical(mask))
    stop(paste("The entries of", sQuote("mask"), "must be logical"))
  if(anyNA(mask)) 
    mask[is.na(mask)] <- FALSE
      
  nc <- ncol(mask)
  nr <- nrow(mask)

  if(!is.null(xy)) {
    ## pixel coordinates given explicitly
    ## validate dimensions
    if(!is.list(xy) || !checkfields(xy, c("x","y")))
      stop("xy should be a list with entries x and y")
    xcol <- xy$x
    yrow <- xy$y
    if(length(xcol) != nc)
      stop(paste("length of xy$x =", length(xcol),
                 "!=", nc, "= number of columns of mask"))
    if(length(yrow) != nr)
      stop(paste("length of xy$y =", length(yrow),
                 "!=", nr, "= number of rows of mask"))
    ## x and y should be evenly spaced
    if(!evenly.spaced(xcol))
      stop("xy$x is not evenly spaced")
    if(!evenly.spaced(yrow))
      stop("xy$y is not evenly spaced")
    ## determine other parameters
    xstep <- diff(xcol)[1L]
    ystep <- diff(yrow)[1L]
    if(missing(xrange) && missing(yrange)) {
      xrange <- range(xcol) + c(-1,1) * xstep/2
      yrange <- range(yrow) + c(-1,1) * ystep/2
    }
  } else {
    ## determine pixel coordinates from xrange, yrange
    if(missing(xrange) && missing(yrange)) {
      ## take pixels to be 1 x 1 unit
      xrange <- c(0,nc)
      yrange <- c(0,nr)
    } else if(check) {
      if(!is.vector(xrange) || length(xrange) != 2 || xrange[2L] < xrange[1L])
        stop("xrange should be a vector of length 2 giving (xmin, xmax)")
      if(!is.vector(yrange) || length(yrange) != 2 || yrange[2L] < yrange[1L])
        stop("yrange should be a vector of length 2 giving (ymin, ymax)")
    }
    xstep <- diff(xrange)/nc
    ystep <- diff(yrange)/nr
    xcol  <- seq(from=xrange[1L]+xstep/2, to=xrange[2L]-xstep/2, length.out=nc)
    yrow  <- seq(from=yrange[1L]+ystep/2, to=yrange[2L]-ystep/2, length.out=nr)
  }

  out <- list(type     = "mask",
              xrange   = unname(xrange),
              yrange   = unname(yrange),
              dim      = c(nr, nc),
              xstep    = unname(xstep),
              ystep    = unname(ystep),
              warnings = .Spatstat.Image.Warning,
              xcol    = unname(xcol), 
              yrow    = unname(yrow),
              m       = mask,
              units   = unitname)
  class(out) <- "owin"
  return(out)
}

#
#-----------------------------------------------------------------------------
#

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

#
#-----------------------------------------------------------------------------
#

as.owin <- function(W, ..., fatal=TRUE) {
  UseMethod("as.owin")
}

as.owin.owin <- function(W, ..., fatal=TRUE) {
  if(verifyclass(W, "owin", fatal=fatal)) 
    return(owin(W$xrange, W$yrange, poly=W$bdry, mask=W$m, unitname=unitname(W), check=FALSE))
  else
    return(NULL)
}

as.owin.ppp <- function(W, ..., fatal=TRUE) {
  if(verifyclass(W, "ppp", fatal=fatal))
    return(W$window)
  else
    return(NULL)
}

as.owin.quad <- function(W, ..., fatal=TRUE) {
  if(verifyclass(W, "quad", fatal=fatal))
    return(W$data$window)
  else
    return(NULL)
}

as.owin.im <- function(W, ..., fatal=TRUE) {
  if(!verifyclass(W, "im", fatal=fatal))
    return(NULL)
  out <- list(type     = "mask",
              xrange   = W$xrange,
              yrange   = W$yrange,
              dim      = W$dim,
              xstep    = W$xstep,
              ystep    = W$ystep,
              warnings = .Spatstat.Image.Warning,
              xcol    = W$xcol,
              yrow    = W$yrow,
              m       = !is.na(W$v),
              units   = unitname(W))
  class(out) <- "owin"
  return(out)
}

as.owin.psp <- function(W, ..., fatal=TRUE) {
  if(!verifyclass(W, "psp", fatal=fatal))
    return(NULL)
  return(W$window)
}

as.owin.tess <- function(W, ..., fatal=TRUE) {
  if(!verifyclass(W, "tess", fatal=fatal))
    return(NULL)
  return(W$window)
}

as.owin.data.frame <- function(W, ..., step, fatal=TRUE) {
  if(!verifyclass(W, "data.frame", fatal=fatal))
    return(NULL)
  if(missing(step)) {
    xstep <- ystep <- NULL
  } else {
    step <- ensure2vector(step)
    xstep <- step[1L]
    ystep <- step[2L]
  }
  if(!(ncol(W) %in% c(2,3))) {
    whinge <- "need exactly 2 or 3 columns of data"
    if(fatal) stop(whinge)
    warning(whinge)
    return(NULL)
  }
  if(twocol <- (ncol(W) == 2)) {
    # assume data is a list of TRUE pixels
    W <- cbind(W, TRUE)
  } 
  mch <- matchNameOrPosition(c("x", "y", "z"), names(W))
  ix <- mch[1L]
  iy <- mch[2L]
  iz <- mch[3L]
  df <- data.frame(x=W[,ix], y=W[,iy], z=as.logical(W[,iz]))
  with(df, {
    xx <- sortunique(x)
    yy <- sortunique(y)
    jj <- match(x, xx)
    ii <- match(y, yy)
    ## make logical matrix (for incomplete x, y sequence)
    ok <- checkbigmatrix(length(xx), length(yy), fatal=fatal)
    if(!ok) return(NULL)
    mm <- matrix(FALSE, length(yy), length(xx))
    mm[cbind(ii,jj)] <- z
    ## ensure xx and yy are complete equally-spaced sequences
    fx <- fillseq(xx, step=xstep)
    fy <- fillseq(yy, step=ystep)
    xcol <- fx[[1L]]
    yrow <- fy[[1L]]
    ## trap very large matrices
    ok <- checkbigmatrix(length(xcol), length(yrow), fatal=fatal)
    if(!ok) return(NULL)
    ## mapping from xx to xcol, yy to yrow
    jjj <- fx[[2L]]
    iii <- fy[[2L]]
    ## make logical matrix for full sequence
    m <- matrix(FALSE, length(yrow), length(xcol))
    m[iii,jjj] <- mm
    ## make binary mask
    out <- owin(mask=m, xy=list(x=xcol, y=yrow))
    ## warn if area fraction is small: may be a misuse of as.owin
    if(twocol) {
      pcarea <- 100 * nrow(df)/prod(dim(m))
      if(pcarea < 1) 
       warning(paste("Window occupies only",
                     paste0(signif(pcarea, 2), "%"),
                     "of frame area. Did you mean owin(poly=df) ?"),
               call.=FALSE)
    }
    return(out)
  })
}

as.owin.default <- function(W, ..., fatal=TRUE) {
  ## Tries to interpret data as an object of class 'owin'
  ## W may be
  ##	a structure with entries xrange, yrange
  ##	a structure with entries xl, xu, yl, yu
  ##	a structure with entries xmin, xmax, ymin, ymax
  ##	a four-element vector (interpreted xmin, xmax, ymin, ymax)
  ##	an object with attribute "bbox"

  if(inherits(W, "box3")) {
    #' cannot be flattened
    if(fatal)
      stop("3D box cannot be converted to a 2D window")
    return(NULL)
  }
  
  if(checkfields(W, c("xrange", "yrange"))) {
    Z <- owinInternalRect(W$xrange, W$yrange)
    return(Z)
  } else if(checkfields(W, c("xmin", "xmax", "ymin", "ymax"))) {
    W <- as.list(W)
    Z <- owinInternalRect(c(W$xmin, W$xmax),c(W$ymin, W$ymax))
    return(Z)
  } else if(checkfields(W, c("xl", "xu", "yl", "yu"))) {
    W <- as.list(W)
    Z <- owinInternalRect(c(W$xl, W$xu),c(W$yl, W$yu))
    return(Z)
  } else if(checkfields(W, c("x", "y", "area"))
            && checkfields(W$area, c("xl", "xu", "yl", "yu"))) {
    V <- as.list(W$area)
    Z <- owinInternalRect(c(V$xl, V$xu),c(V$yl, V$yu))
    return(Z)
  } else if(is.vector(W) && is.numeric(W) && length(W) == 4) {
    Z <- owinInternalRect(W[1:2], W[3:4])
    return(Z)
  } else if(!is.null(Z <- attr(W, "bbox"))) {
    return(as.owin(Z, ..., fatal=fatal))
  } else if(inherits(W, c("SpatialPolygons", "SpatialPolygonsDataFrame"))) {
    gripe <- "The package 'maptools' is needed to convert this data type"
    if(fatal) stop(gripe, call.=FALSE) else warning(gripe, call.=FALSE)
    return(NULL)
  } else {
    #' no idea
    if(fatal) stop("Can't interpret W as a window", call.=FALSE)
    return(NULL)
  }
}		

#
#-----------------------------------------------------------------------------
#
#
Frame <- function(X) { UseMethod("Frame") }

"Frame<-" <- function(X, value) { UseMethod("Frame<-") }

Frame.default <- function(X) { as.rectangle(X) }

"Frame<-.default" <- function(X, value) {
  Frame(Window(X)) <- value
  return(X)
}

## .........................................................

as.rectangle <- function(w, ...) {
  if(inherits(w, "owin"))
    return(owinInternalRect(w$xrange, w$yrange, unitname=unitname(w), check=FALSE))
  if(inherits(w, "im"))
    return(owinInternalRect(w$xrange, w$yrange, unitname=unitname(w), check=FALSE))
  if(inherits(w, "ppp")) 
    return(owinInternalRect(w$window$xrange, w$window$yrange, unitname=unitname(w$window),
                            check=FALSE))
  if(inherits(w, "layered")) 
    return(do.call(boundingbox, unname(lapply(w, as.rectangle, ...))))
  w <- as.owin(w, ...)
  return(owinInternalRect(w$xrange, w$yrange, unitname=unitname(w), check=FALSE))
}

#
#-----------------------------------------------------------------------------
#
as.mask <- function(w, eps=NULL, dimyx=NULL, xy=NULL,
                    rule.eps=c("adjust.eps", "grow.frame", "shrink.frame")) {
  ## 	  eps:		   grid mesh (pixel) size
  ##	dimyx:		   dimensions of pixel raster
  ##       xy:             coordinates of pixel raster
  ## rule.eps:             rule for adjusting frame when 'eps' is given
  ###########################
  ##  First determine window
  ###########################
  w.absent <- missing(w) || is.null(w)
  if(w.absent) {
    ## w is not given
    ## Ensure there is some window information
    if(is.null(xy)) 
      stop("If w is missing, xy is required")
    uname <- unitname(xy) # could be null
  } else {
    ## w is given, but may require conversion
    ## Handle cases where w contains pixel data
    if(is.data.frame(w))
      return(owin(mask=w, xy=xy))
    if(is.matrix(w)) {
      w <- as.data.frame(w)
      colnames(w) <- c("x", "y")
      return(owin(mask=w, xy=xy))
    }
    ## Handle cases where w can be converted to a window
    w <- as.owin(w)
    ## Already a mask?
    if(is.mask(w) && is.null(eps) && is.null(dimyx) && is.null(xy)) {
      ## w contained raster data, and no further raster data is provided
      return(w)
    }
    uname <- unitname(w)
    xrange <- w$xrange
    yrange <- w$yrange
  }

  #####################################
  ##  Next determine pixel coordinates
  #####################################
  if(is.null(xy)) {
    ## Pixel coordinates to be computed from other dimensions
    ## First determine row & column dimensions
    if(!is.null(dimyx)) {
      ## Pixel dimensions are given
      dimyx <- ensure2vector(dimyx)
      nr <- dimyx[1L]
      nc <- dimyx[2L]
    } else {
      ## use pixel size 'eps'
      if(!is.null(eps)) {
        eps <- ensure2vector(eps)
        width <- diff(xrange)
        height <- diff(yrange)
        nc <- width/eps[1L]
        nr <- height/eps[2L]
        fc <- nc %% 1
        fr <- nr %% 1
        rule.eps <- match.arg(rule.eps)
        switch(rule.eps,
               adjust.eps = {
                 ## Frame size is fixed; pixel size will be adjusted
                 nc <- ceiling(nc)
                 nr <- ceiling(nr)
                 if(fc != 0) eps[1L] <- width/nc
                 if(fr != 0) eps[2L] <- height/nr
               },
               grow.frame = {
                 ## Pixel size is fixed; frame will be expanded
                 nc <- ceiling(nc)
                 nr <- ceiling(nr)
                 if(fc != 0) 
                   xrange <- mean(xrange) + c(-1,1) * nc * eps[1L]/2
                 if(fr != 0) 
                   yrange <- mean(yrange) + c(-1,1) * nr * eps[2L]/2
               },
               shrink.frame = {
                 ## Pixel size is fixed; frame will be contracted
                 nc <- floor(nc)
                 nr <- floor(nr)
                 if(nc <= 0 || nr <= 0) 
                   stop(paste("pixel size argument eps exceeds frame size;",
                              "unable to shrink frame"))
                 if(fc != 0) 
                   xrange <- mean(xrange) + c(-1,1) * nc * eps[1L]/2
                 if(fr != 0) 
                   yrange <- mean(yrange) + c(-1,1) * nr * eps[2L]/2
               })
      } else {
        ## use spatstat default dimensions
        np <- spatstat.options("npixel")
        if(length(np) == 1)
          nr <- nc <- np[1L]
        else {
          nr <- np[2L]  
          nc <- np[1L]
        }
      }
    }
    if((mpix <- (nr * nc)/1048576) >= 10) {
      whinge <- paste("Creating",
                      articlebeforenumber(mpix),
                      paste0(round(mpix, 1), "-megapixel"),
                      "window mask")
      message(whinge)
      warning(whinge, call.=FALSE)
    }
    ## Initialise mask with all entries TRUE
    nowidth  <- (diff(xrange) < .Machine$double.eps)
    noheight <- (diff(yrange) < .Machine$double.eps)
    if(nowidth || noheight) {
      if(nowidth && noheight) {
        nr <- nc <- 1
      } else if(noheight) {
        nr <- 1
        yrange <- yrange + c(-1/2,1/2) * diff(xrange)/(nc+1)
      } else if(nowidth) {
        ## ensure square pixels
        nc <- 1
        xrange <- xrange + c(-1/2,1/2) * diff(yrange)/(nr+1)
      }
    }
    rasta <- owin(xrange, yrange, mask=matrix(TRUE, nr, nc))
  } else {
    ## 
    ## Pixel coordinates given explicitly:
    ##    xy is an image, a mask, or a list(x,y)
    if(is.im(xy)) {
      rasta <- as.owin(xy)
      rasta$m[] <- TRUE
    } else if(is.owin(xy)) {
      if(xy$type != "mask")
        stop("argument xy does not contain raster coordinates.")
      rasta <- xy
      rasta$m[] <- TRUE
    } else {
      if(!checkfields(xy, c("x", "y")))
        stop(paste(sQuote("xy"),
                   "should be a list containing two vectors x and y"))
      x <- sortunique(xy$x)
      y <- sortunique(xy$y)
      ## derive other parameters
      nr <- length(y)
      nc <- length(x)
      ## check size
      if((mpix <- (nr * nc)/1048576) >= 10) {
        whinge <- paste("Creating",
                        articlebeforenumber(mpix),
                        paste0(round(mpix, 1), "-megapixel"),
                        "window mask")
        message(whinge)
        warning(whinge, call.=FALSE)
      }
      ## x and y pixel sizes
      dx <- diff(x)
      if(diff(range(dx)) > 0.01 * mean(dx))
        stop("x coordinates must be evenly spaced")
      xstep <- mean(dx)
      dy <- diff(y)
      if(diff(range(dy)) > 0.01 * mean(dy))
        stop("y coordinates must be evenly spaced")
      ystep <- mean(dy)
      xr <- range(x)
      yr <- range(y)
      xrange <-  xr + xstep * c(-1,1)/2
      yrange <-  yr + ystep * c(-1,1)/2
      ## initialise mask with all entries TRUE
      rasta <- list(type     = "mask",
                    xrange   = xrange,
                    yrange   = yrange,
                    dim      = c(nr, nc),
                    xstep    = xstep,
                    ystep    = ystep,
                    warnings = .Spatstat.Image.Warning,
                    xcol    = seq(from=xr[1L], to=xr[2L], length.out=nc),
                    yrow    = seq(from=yr[1L], to=yr[2L], length.out=nr),
                    m       = matrix(TRUE, nr, nc),
                    units   = uname)
      class(rasta) <- "owin"
    }
    if(w.absent) {
      ## No more window information
      out <- rasta
      if(!(identical(x, xy$x) && identical(y, xy$y))) {
        ## xy is an enumeration of the TRUE pixels
        out$m[] <- FALSE
        ij <- cbind(i=match(xy$y, y), j=match(xy$x, x))
        out$m[ij] <- TRUE
      }
      return(out)
    }
  }

  ####################################################
  ## Finally, mask pixel raster with existing window
  #################################################### 
  switch(w$type,
         rectangle = {
           out <- rasta
           wxrange <- w$xrange
           wyrange <- w$yrange
           if(!all(wxrange == rasta$xrange)
              || !all(wyrange == rasta$yrange)) {
             xcol <- rasta$xcol
             yrow <- rasta$yrow
             badrow <- which(yrow > wyrange[2L] | yrow < wyrange[1L])
             badcol <- which(xcol > wxrange[2L] | xcol < wxrange[1L])
             out$m[badrow , ] <- FALSE
             out$m[ , badcol] <- FALSE
           }
         },
         mask = {
           ## resample existing mask on new raster
           out <- rastersample(w, rasta)
         },
         polygonal = {
           ## use C code
           out <- owinpoly2mask(w, rasta, FALSE)
         })

  unitname(out) <- uname
  return(out)
}

as.matrix.owin <- function(x, ...) {
  m <- as.mask(x, ...)
  return(m$m)
}

#
#
#-----------------------------------------------------------------------------
#
as.polygonal <- function(W, repair=FALSE) {
  verifyclass(W, "owin")
  switch(W$type,
         rectangle = {
           xr <- W$xrange
           yr <- W$yrange
           return(owin(xr, yr, poly=list(x=xr[c(1,2,2,1)],y=yr[c(1,1,2,2)]),
                       unitname=unitname(W),
                       check=FALSE))
         },
         polygonal = {
           if(repair)
             W <- owin(poly=W$bdry, unitname=unitname(W))
           return(W)
         },
         mask = {
           # This could take a while
           M <- W$m
           nr <- nrow(M)
           notM <- !M
           xcol <- W$xcol
           yrow <- W$yrow
           ## determine resolution for polyclip operations
           eps <- max(W$xstep, W$ystep)/(2^31)
           eps <- max(eps, 4 * .Machine$double.eps)
           p <- list(x0 = xcol[1],
                     y0 = yrow[1],
                     eps = eps)
           ## pixels will be slightly expanded to avoid artefacts
           xbracket <- c(-1,1) * (W$xstep/2 + 4 * eps)
           ybracket <- c(-1,1) * (W$ystep/2 + 4 * eps)
           # identify runs of TRUE entries in each column
           start <- M & rbind(TRUE, notM[-nr, ])
           finish <- M & rbind(notM[-1, ], TRUE)
           #' build result
           out <- NULL
           for(j in 1:ncol(M)) {
             xj <- xcol[j]
             # identify start and end positions in column j
             starts <- which(start[,j])
             finishes <- which(finish[,j])
             ns <- length(starts)
             nf <- length(finishes)
             if(ns != nf)
               stop(paste("Internal error: length(starts)=", ns,
                          ", length(finishes)=", nf))
             if(ns > 0) {
               for(k in 1:ns) {
                 yfrom <- yrow[starts[k]]
                 yto   <- yrow[finishes[k]]
                 yk <- sort(c(yfrom,yto))
                 #' make rectangle boundary in reversed orientation
                 xrect <- xj + xbracket
                 yrect <- yk + ybracket
                 recto <- list(list(x = xrect[c(1,2,2,1)],
                                    y = yrect[c(2,2,1,1)]))
                 #' add to result
                 if(is.null(out)) {
                   out <- recto
                 } else {
                   out <- polyclip::polyclip(out, recto, "union",
                                             fillA="nonzero", fillB="nonzero",
                                             eps = p$eps,
                                             x0  = p$x0,
                                             y0  = p$y0)
                 }
               }
             }
           }
           if(is.null(out)) return(emptywindow(Frame(W)))
           totarea <- sum(sapply(out, Area.xypolygon))
           if(totarea < 0)
             out <- lapply(out, reverse.xypolygon)
           out <- owin(poly=out, check=FALSE, unitname=unitname(W))
           return(out)
         }
         )
}

#
# ----------------------------------------------------------------------

is.polygonal <- function(w) {
  return(inherits(w, "owin") && (w$type == "polygonal"))
}

is.rectangle <- function(w) {
  return(inherits(w, "owin") && (w$type == "rectangle"))
}

is.mask <- function(w) {
  return(inherits(w, "owin") && (w$type == "mask"))
}

validate.mask <- function(w, fatal=TRUE) {
  verifyclass(w, "owin", fatal=fatal)
  if(w$type == "mask")
    return(TRUE)
  if(fatal)
      stop(paste(short.deparse(substitute(w)), "is not a binary mask"))
  else {
      warning(paste(short.deparse(substitute(w)), "is not a binary mask"))
      return(FALSE)
  }
}

dim.owin <- function(x) { return(x$dim) } ## NULL unless it's a mask

## internal use only:

rasterx.mask <- function(w, drop=FALSE) {
  validate.mask(w)
  x <- w$xcol[col(w)]
  x <- if(drop) x[w$m, drop=TRUE] else array(x, dim=w$dim)
  return(x)
}

rastery.mask <- function(w, drop=FALSE) {
  validate.mask(w)
  y <- w$yrow[row(w)]
  y <- if(drop) y[w$m, drop=TRUE] else array(y, dim=w$dim)
  return(y)
}

rasterxy.mask <- function(w, drop=FALSE) {
  validate.mask(w)
  x <- w$xcol[col(w)]
  y <- w$yrow[row(w)]
  if(drop) {
    m <- w$m
    x <- x[m, drop=TRUE] 
    y <- y[m, drop=TRUE]
  }
  return(list(x=as.numeric(x),
              y=as.numeric(y)))
}


nearest.raster.point <- function(x,y,w, indices=TRUE) {
  stopifnot(is.mask(w) || is.im(w))
  nr <- w$dim[1L]
  nc <- w$dim[2L]
  if(length(x) == 0) {
    cc <- rr <- integer(0)
  } else {
    cc <- 1 + round((x - w$xcol[1L])/w$xstep)
    rr <- 1 + round((y - w$yrow[1L])/w$ystep)
    cc <- pmax.int(1,pmin.int(cc, nc))
    rr <- pmax.int(1,pmin.int(rr, nr))
  }
  if(indices) 
    return(list(row=rr, col=cc))
  else
    return(list(x=w$xcol[cc], y=w$yrow[rr]))
}

mask2df <- function(w) {
  stopifnot(is.owin(w) && w$type == "mask")
  xx <- raster.x(w)
  yy <- raster.y(w)
  ok <- w$m
  xx <- as.vector(xx[ok])
  yy <- as.vector(yy[ok])
  return(data.frame(x=xx, y=yy))
}

#------------------------------------------------------------------
		
complement.owin <- function(w, frame=as.rectangle(w)) {
  w <- as.owin(w)

  if(reframe <- !missing(frame)) {
    verifyclass(frame, "owin")
    w <- rebound.owin(w, frame)
    # if w was a rectangle, it's now polygonal
  }

  switch(w$type,
         mask = {
           w$m <- !(w$m)
         },
         rectangle = {
           # return empty window
           return(emptywindow(w))
         },
         polygonal = {
           bdry <- w$bdry
           if(length(bdry) == 0) {
             # w is empty
             return(frame)
           }
           # bounding box, in anticlockwise order
           box <- list(x=w$xrange[c(1,2,2,1)],
                       y=w$yrange[c(1,1,2,2)])
           boxarea <- Area.xypolygon(box)
                 
           # first check whether one of the current boundary polygons
           # is the bounding box itself (with + sign)
           if(reframe)
             is.box <- rep.int(FALSE, length(bdry))
           else {
             nvert <- lengths(lapply(bdry, getElement, name="x"))
             areas <- sapply(bdry, Area.xypolygon)
             boxarea.mineps <- boxarea * (0.99999)
             is.box <- (nvert == 4 & areas >= boxarea.mineps)
             if(sum(is.box) > 1)
               stop("Internal error: multiple copies of bounding box")
             if(all(is.box)) {
               return(emptywindow(box))
             }
           }
                 
           # if box is present (with + sign), remove it
           if(any(is.box))
             bdry <- bdry[!is.box]
                 
           # now reverse the direction of each polygon
           bdry <- lapply(bdry, reverse.xypolygon, adjust=TRUE)

           # if box was absent, add it
           if(!any(is.box))
             bdry <- c(bdry, list(box))   # sic

           # put back into w
           w$bdry <- bdry
         },
         stop("unrecognised window type", w$type)
         )
  return(w)
}

#-----------------------------------------------------------

inside.owin <- function(x, y, w) {
  # test whether (x,y) is inside window w
  # x, y may be vectors 

  if((missing(y) || is.null(y)) && all(c("x", "y") %in% names(x))) {
    y <- x$y
    x <- x$x
  }

  w <- as.owin(w)

  if(length(x)==0)
    return(logical(0))
  
  # test whether inside bounding rectangle
  xr <- w$xrange
  yr <- w$yrange
  eps <- sqrt(.Machine$double.eps)
  frameok <- (x >= xr[1L] - eps) & (x <= xr[2L] + eps) & 
             (y >= yr[1L] - eps) & (y <= yr[2L] + eps)
 
  if(!any(frameok))  # all points OUTSIDE window - no further work needed
    return(frameok)

  ok <- frameok
  switch(w$type,
         rectangle = {
           return(ok)
         },
         polygonal = {
           ## check scale
           framesize <- max(diff(xr), diff(yr))
           if(framesize > 1e6 || framesize < 1e-6) {
             ## rescale to avoid numerical overflow
             scalefac <- as.numeric(framesize)/100
             w <- as.polygonal(rescale(w, scalefac))
             x <- x/scalefac
             y <- y/scalefac
           }
           xy <- list(x=x,y=y)
           bdry <- w$bdry
           total <- numeric(length(x))
           on.bdry <- rep.int(FALSE, length(x))
           for(i in seq_along(bdry)) {
             score <- inside.xypolygon(xy, bdry[[i]], test01=FALSE)
             total <- total + score
             on.bdry <- on.bdry | attr(score, "on.boundary")
           }
           # any points identified as belonging to the boundary get score 1
           total[on.bdry] <- 1
           # check for sanity now..
           uhoh <- (total * (1-total) != 0)
           if(any(uhoh)) {
             nuh <- sum(uhoh)
             warning(paste("point-in-polygon test had difficulty with",
                           nuh,
                           ngettext(nuh, "point", "points"),
                           "(total score not 0 or 1)"),
                     call.=FALSE)
             total[uhoh] <- 0
           }
           return(ok & (total != 0))
         },
         mask = {
           # consider only those points which are inside the frame
           xf <- x[frameok]
           yf <- y[frameok]
           # map locations to raster (row,col) coordinates
           loc <- nearest.raster.point(xf,yf,w)
           # look up mask values
           okf <- (w$m)[cbind(loc$row, loc$col)]
           # insert into 'ok' vector
           ok[frameok] <- okf
           return(ok)
         },
         stop("unrecognised window type", w$type)
         )
}

#-------------------------------------------------------------------------
  
print.owin <- function(x, ..., prefix="window: ") {
  verifyclass(x, "owin")
  unitinfo <- summary(unitname(x))
  switch(x$type,
         rectangle={
           rectname <- paste0(prefix, "rectangle =")
         },
         polygonal={
           nonemp <- (length(x$bdry) != 0)
           splat(paste0(prefix,
                        if(nonemp) "polygonal boundary" else "empty"))
           rectname <- "enclosing rectangle:"
         },
         mask={
           splat(paste0(prefix, "binary image mask"))
           di <- x$dim
           splat(di[1L], "x", di[2L], "pixel array (ny, nx)")
           rectname <- "enclosing rectangle:"
         }
         )
  splat(rectname,
        prange(zapsmall(x$xrange)),
        "x",
        prange(zapsmall(x$yrange)),
        unitinfo$plural,
        unitinfo$explain)
  invisible(NULL)
}

summary.owin <- function(object, ...) {
  verifyclass(object, "owin")
  result <- list(xrange=object$xrange,
                 yrange=object$yrange,
                 type=object$type,
                 area=area(object),
                 units=unitname(object))
  result$areafraction <- with(result, area/(diff(xrange) * diff(yrange)))
  switch(object$type,
         rectangle={
         },
         polygonal={
           poly <- object$bdry
           result$npoly <- npoly <- length(poly)
           if(npoly == 0) {
             result$areas <- result$nvertices <- numeric(0)
           } else if(npoly == 1) {
             result$areas <- Area.xypolygon(poly[[1L]])
             result$nvertices <- length(poly[[1L]]$x)
           } else {
             result$areas <- unlist(lapply(poly, Area.xypolygon))
             result$nvertices <- lengths(lapply(poly, getElement, name="x"))
           }
           result$nhole <- sum(result$areas < 0)
         },
         mask={
           result$npixels <- object$dim
           result$xstep <- object$xstep
           result$ystep <- object$ystep
         }
         )
  class(result) <- "summary.owin"
  result
}

print.summary.owin <- function(x, ...) {
  verifyclass(x, "summary.owin")
  unitinfo <- summary(x$units)
  pluralunits <- unitinfo$plural
  singularunits <- unitinfo$singular
  switch(x$type,
         rectangle={
           rectname <- "Window: rectangle ="
         },
         polygonal={
           np <- x$npoly
           splat("Window: polygonal boundary")
           if(np == 0) {
             splat("window is empty")
           } else if(np == 1) {
             splat("single connected closed polygon with",
                   x$nvertices, "vertices")
           } else {
             nh <- x$nhole
             holy <- if(nh == 0) "(no holes)" else
                     if(nh == 1) "(1 hole)" else
                     paren(paste(nh, "holes"))
             splat(np, "separate polygons", holy)
             if(np > 0)
               print(data.frame(vertices=x$nvertices,
                                area=signif(x$areas, 6),
                                relative.area=signif(x$areas/x$area,3),
                                row.names=paste("polygon",
                                  1:np,
                                  ifelse(x$areas < 0, "(hole)", "")
                                  )))
           }
           rectname <- "enclosing rectangle:"
         },
         mask={
           splat("binary image mask")
           di <- x$npixels
           splat(di[1L], "x", di[2L], "pixel array (ny, nx)")
           splat("pixel size:",
                 signif(x$xstep,3), "by", signif(x$ystep,3),
                 pluralunits)
           rectname <- "enclosing rectangle:"
         }
         )
  splat(rectname,
        prange(zapsmall(x$xrange)),
        "x",
        prange(zapsmall(x$yrange)),
        pluralunits)
  if(x$xrange[1] != 0 || x$yrange[1] != 0) {
    width <- diff(x$xrange)
    height <- diff(x$yrange)
    blank <- paste(rep(" ", nchar(rectname)), collapse="") 
    splat(blank, paren(paste(signif(width, 4), "x", signif(height, 4),
                             pluralunits)))
  }
  Area <- signif(x$area, 6)
  splat("Window area =", Area, "square",
        if(Area == 1) singularunits else pluralunits)
  if(!is.null(ledge <- unitinfo$legend))
    splat(ledge)
  if(x$type != "rectangle")
    splat("Fraction of frame area:", signif(x$areafraction, 3))
  return(invisible(x))
}

as.data.frame.owin <- function(x, ..., drop=TRUE) {
  stopifnot(is.owin(x))
  switch(x$type,
         rectangle = { x <- as.polygonal(x) },
         polygonal = { },
         mask = {
           xy <- rasterxy.mask(x, drop=drop)
           if(!drop) xy <- append(xy, list(inside=as.vector(x$m)))
           return(as.data.frame(xy, ...))
         })
  b <- x$bdry
  ishole <- sapply(b, is.hole.xypolygon)
  sign <- (-1)^ishole
  b <- lapply(b, as.data.frame, ...)
  nb <- length(b)
  if(nb == 1)
    return(b[[1L]])
  dfs <- mapply(cbind, b, id=as.list(seq_len(nb)), sign=as.list(sign),
                SIMPLIFY=FALSE)
  df <- do.call(rbind, dfs)
  return(df)
}

discretise <- function(X, eps=NULL, dimyx=NULL, xy=NULL,
                       move.points=FALSE,
                       rule.eps=c("adjust.eps", "grow.frame", "shrink.frame")) {
  verifyclass(X,"ppp")
  W <- X$window
  ok <- inside.owin(X$x,X$y,W)
  if(!all(ok))
    stop("There are points of X outside the window of X")
  new.grid <- !is.null(eps) || !is.null(dimyx) || !is.null(xy)
  if(new.grid) rule.eps <- match.arg(rule.eps)
  new.mask <- new.grid || !is.mask(W)
  WM <- if(!new.mask) W else as.mask(W,
                             eps=eps,dimyx=dimyx,xy=xy,rule.eps=rule.eps) 
  if(move.points) {
    ## move points to pixel centres 
    if(new.mask) X$window <- WM
    indices <- as.data.frame(nearest.valid.pixel(X$x, X$y, WM, nsearch=2))
    xnew <- WM$xcol[indices$col]
    ynew <- WM$yrow[indices$row]
    ## insurance:
    if(any(notfound <- !complete.cases(indices))) {
      XP <- project2set(X[notfound], WM)
      xnew[notfound] <- XP$x
      ynew[notfound] <- XP$y
    }
    X$x <- xnew
    X$y <- ynew
  } else if(new.mask) {
    ## ensure points are inside new window
    nok <- !inside.owin(X$x,X$y,WM)
    if(any(nok)) {
      ifix <- nearest.raster.point(X$x[nok],X$y[nok], WM)
      ifix <- cbind(ifix$row,ifix$col)
      WM$m[ifix] <- TRUE
    }
    X$window <- WM
  }
  return(X)
}

pixelcentres <- function (X, W=NULL,...) {
  X <- as.mask(as.owin(X), ...)
  if(is.null(W)) W <- as.rectangle(X)
  Y <- as.ppp(raster.xy(X,drop=TRUE),W=W)
  return(Y)
}

owin2polypath <- function(w) {
  w <- as.polygonal(w)
  b <- w$bdry
  xvectors <- lapply(b, getElement, name="x")
  yvectors <- lapply(b, getElement, name="y")
  xx <- unlist(lapply(xvectors, append, values=NA, after=FALSE))[-1]
  yy <- unlist(lapply(yvectors, append, values=NA, after=FALSE))[-1]
  return(list(x=xx, y=yy))
}

## generics which extract and assign the window of some object

Window <- function(X, ...) { UseMethod("Window") }

"Window<-" <- function(X, ..., value) { UseMethod("Window<-") }