Nothing
#' @title Picka's Reduced Window Estimator of Coverage Probability
#' @author Kassel Liam Hingee
#' @import spatstat
#' @export cppicka
#'
#' @description
#' This function provides estimates of coverage probability from subsets of the observation window,
#' which are a key component of balanced estimators of covariance, centred covariance, pair-correlation and gliding box lacunarity.
#' @param xi An observation of a RACS of interest as a full binary map (as an \code{im} object) or as the foreground set (as an \code{owin} object).
#' In the latter case the observation window, \code{obswin}, must be supplied.
#' @param obswin If \code{xi} is an \code{owin} object then \code{obswin} is an
#' \code{owin} object that specifies the observation window.
#' @param setcov_boundarythresh To avoid instabilities caused by dividing by very small quantities, if the set covariance of the observation window
#' is smaller than \code{setcov_boundarythresh}, then the returned pixel value is NA.
#' @return An \code{im} object. Pixel values correspond to estimates of the coverage probability
#' from the subregion of the observation window, \eqn{W}, that is the intersection of \eqn{W} and \eqn{W} shifted by vector \eqn{v}, where \eqn{v} is the pixel location.
#'
#' @details
#' The plug-in moment covariance estimator (\code{\link{plugincvc}}) uses less of the observation window than the usual coverage probability estimators.
#' Picka (1997, 2000) created new 'balanced' estimators of centred covariance and pair-correlation
#' that accounted for this difference.
#' A key component of Picka's estimators is an estimate of the coverage probability from the subregion of the binary map that is
#' the intersection between \eqn{W} and \eqn{W} shifted by vector \eqn{v}, where \eqn{W} is the observation window (p.~687, Picka, 2000).
#' If we treat \eqn{X} and \eqn{W} as indicator functions representing the foreground and observation window respectively,
#' this coverage probability estimator used by Picka is
#' \deqn{ \frac{\int X(u) W(u) W(u - v) du} {\int W(u) W(u - v) du}. }{integral(X(u) W(u) W(u - v) du) / integral(W(u) W(u - v) du).}
#'
#' \code{cppicka} produces these estimates for an array of vectors \eqn{v} using fast Fourier transforms.
#' @references
#' Picka, J.D. (1997) \emph{Variance-Reducing Modifications for Estimators of Dependence in Random Sets}. Ph.D.: Illinois, USA: The University of Chicago.
#'
#' Picka, J.D. (2000) Variance reducing modifications for estimators of standardized moments of random sets. \emph{Advances in Applied Probability}, 32, 682-700.
#' @examples
#' xi <- heather$coarse
#' obswindow <- Frame(heather$coarse)
#' cp <- coverageprob(xi, obswindow)
#' cpp1 <- cppicka(xi, obswindow)
#' @keywords spatial nonparametric
cppicka <- function(xi, obswin = NULL,
setcov_boundarythresh = NULL) {
if (is.im(xi)){
if(!is.null(obswin)){xi[setminus.owin(Frame(xi), obswin)] <- 0}
#check that xi is only 1s, 0s and NAs
isbinarymap(xi, requiretrue = TRUE)
if(!is.null(obswin)){obswin <- as.im(obswin, W = Frame(obswin), xy = xi) * xi}
else {obswin <- xi}
#saving all NAs as 0 and everything else as 1 in obswin
obswin[!is.na(as.matrix(obswin))] <- 1
obswin[is.na(as.matrix(obswin))] <- 0
#saving all NAs as 0 and everything else as 1 in xi
xi[is.na(as.matrix(xi))] <- 0 #turn all NA's in xi to 0s
}
if (is.owin(xi)){
stopifnot(is.owin(obswin))
if (!is.subset.owin(boundingbox(xi), boundingbox(obswin))){
warning("Bounding box of xi is not inside bounding box of obswin. This can be caused by xi being a pixel mask, and obswin being a polygon that is not quite aligned with the pixels. The Frame(xi) will be enlarged")
}
#the following makes sure convolutions generated when the input frame of xi is much smaller than obswin
Frame(xi) <- boundingbox(obswin, Frame(xi)) #that obswin and xi are needed on right is when the binary grid in xi doesn't quite match with obswin
xi <- as.im(xi, W = obswin, na.replace = 0)
obswin <- as.im(obswin, na.replace = 0, xy = xi)
}
#numerator
top <- convolve.im(xi, obswin, reflectY = TRUE) #u in (obswin + v) iff (u - v) in obswin. Thus Y needs to have argument reflected
unitname(top) <- unitname(xi)
#denominator
bot <- convolve.im(obswin, obswin, reflectY = TRUE)
#now xi and obswin should be of the same style regardless of input
if (is.null(setcov_boundarythresh)){
setcov_boundarythresh <- 0.1 * sum(obswin)*obswin$xstep*obswin$ystep
} else if (setcov_boundarythresh < bot$xstep * bot$ystep * 1E-8){
warning("setcov_boundarythresh is smaller than A*1E-8 where A is the size of a pixel.
This might be smaller than the precision of the set covariance computations.
Consider setting setcov_boundarythresh higher.")
}
bot[bot < setcov_boundarythresh] <- NA #to remove small denominators
unitname(bot) <- unitname(obswin)
return(eval.im( top / bot))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.