R/window.R

Defines functions Window owin2polypath discretise as.data.frame.owin print.summary.owin summary.owin print.owin inside.owin complement.owin mask2df nearest.raster.point rasterxy.mask rastery.mask rasterx.mask dim.owin validate.mask is.mask is.rectangle is.polygonal as.polygonal as.matrix.owin as.mask as.rectangle Frame.default Frame as.owin.default as.owin.data.frame as.owin.tess as.owin.psp as.owin.im as.owin.quad as.owin.ppp as.owin.owin as.owin is.owin owinInternalMask owinInternalPoly asXYdata isXYdata owinInternalRect owin

Documented in as.data.frame.owin as.mask as.matrix.owin as.owin as.owin.data.frame as.owin.default as.owin.im as.owin.owin as.owin.ppp as.owin.psp as.owin.quad as.owin.tess as.polygonal as.rectangle complement.owin dim.owin discretise Frame Frame.default inside.owin is.mask is.owin is.polygonal is.rectangle mask2df nearest.raster.point owin owin2polypath owinInternalMask owinInternalPoly owinInternalRect print.owin print.summary.owin rasterx.mask rasterxy.mask rastery.mask summary.owin validate.mask Window

#
#	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<-") }

Try the spatstat.geom package in your browser

Any scripts or data that you put into this service are public.

spatstat.geom documentation built on Sept. 18, 2024, 9:08 a.m.