closepairs.R
#
# closepairs.R
#
# $Revision: 1.19 $ $Date: 2012/03/18 08:53:45 $
#
# simply extract the r-close pairs from a dataset
#
# Less memory-hungry for large patterns
#
closepairs <- function(X, rmax, ordered=TRUE) {
verifyclass(X, "ppp")
stopifnot(is.numeric(rmax) && length(rmax) == 1 && rmax >= 0)
npts <- npoints(X)
if(npts == 0)
return(list(i=integer(0),
j=integer(0),
xi=numeric(0),
yi=numeric(0),
xj=numeric(0),
yj=numeric(0),
dx=numeric(0),
dy=numeric(0),
d=numeric(0)))
# sort points by increasing x coordinate
oo <- fave.order(X$x)
Xsort <- X[oo]
# First make an OVERESTIMATE of the number of pairs
nsize <- ceiling(4 * pi * (npts^2) * (rmax^2)/area.owin(X$window))
nsize <- max(1024, nsize)
# Now extract pairs
if(spatstat.options("closepairs.newcode")) {
# ------------------- use new faster code ---------------------
x <- Xsort$x
y <- Xsort$y
r <- rmax
ng <- nsize
storage.mode(x) <- "double"
storage.mode(y) <- "double"
storage.mode(r) <- "double"
storage.mode(ng) <- "integer"
z <- .Call("Vclosepairs",
xx=x, yy=y, rr=r, nguess=ng,
PACKAGE="spatstat")
if(length(z) != 9)
stop("Internal error: incorrect format returned from Vclosepairs")
i <- z[[1]] # NB no increment required
j <- z[[2]]
xi <- z[[3]]
yi <- z[[4]]
xj <- z[[5]]
yj <- z[[6]]
dx <- z[[7]]
dy <- z[[8]]
d <- z[[9]]
} else {
# ------------------- use older code --------------------------
DUP <- spatstat.options("dupC")
z <-
.C("closepairs",
nxy=as.integer(npts),
x=as.double(Xsort$x),
y=as.double(Xsort$y),
r=as.double(rmax),
noutmax=as.integer(nsize),
nout=as.integer(integer(1)),
iout=as.integer(integer(nsize)),
jout=as.integer(integer(nsize)),
xiout=as.double(numeric(nsize)),
yiout=as.double(numeric(nsize)),
xjout=as.double(numeric(nsize)),
yjout=as.double(numeric(nsize)),
dxout=as.double(numeric(nsize)),
dyout=as.double(numeric(nsize)),
dout=as.double(numeric(nsize)),
status=as.integer(integer(1)),
DUP=DUP,
PACKAGE="spatstat")
if(z$status != 0) {
# Guess was insufficient
# Obtain an OVERCOUNT of the number of pairs
# (to work around gcc bug #323)
rmaxplus <- 1.25 * rmax
nsize <- .C("paircount",
nxy=as.integer(npts),
x=as.double(Xsort$x),
y=as.double(Xsort$y),
rmaxi=as.double(rmaxplus),
count=as.integer(integer(1)),
DUP=DUP,
PACKAGE="spatstat")$count
if(nsize <= 0)
return(list(i=integer(0),
j=integer(0),
xi=numeric(0),
yi=numeric(0),
xj=numeric(0),
yj=numeric(0),
dx=numeric(0),
dy=numeric(0),
d=numeric(0)))
# add a bit more for safety
nsize <- ceiling(1.1 * nsize) + 2 * npts
# now extract points
z <-
.C("closepairs",
nxy=as.integer(npts),
x=as.double(Xsort$x),
y=as.double(Xsort$y),
r=as.double(rmax),
noutmax=as.integer(nsize),
nout=as.integer(integer(1)),
iout=as.integer(integer(nsize)),
jout=as.integer(integer(nsize)),
xiout=as.double(numeric(nsize)),
yiout=as.double(numeric(nsize)),
xjout=as.double(numeric(nsize)),
yjout=as.double(numeric(nsize)),
dxout=as.double(numeric(nsize)),
dyout=as.double(numeric(nsize)),
dout=as.double(numeric(nsize)),
status=as.integer(integer(1)),
DUP=DUP,
PACKAGE="spatstat")
if(z$status != 0)
stop(paste("Internal error: C routine complains that insufficient space was allocated:", nsize))
}
# trim vectors to the length indicated
npairs <- z$nout
if(npairs <= 0)
return(list(i=integer(0),
j=integer(0),
xi=numeric(0),
yi=numeric(0),
xj=numeric(0),
yj=numeric(0),
dx=numeric(0),
dy=numeric(0),
d=numeric(0)))
actual <- seq_len(npairs)
i <- z$iout[actual] # sic
j <- z$jout[actual]
xi <- z$xiout[actual]
yi <- z$yiout[actual]
xj <- z$xjout[actual]
yj <- z$yjout[actual]
dx <- z$dxout[actual]
dy <- z$dyout[actual]
d <- z$dout[actual]
# ------------------- end code switch ------------------------
}
# convert i,j indices to original sequence
i <- oo[i]
j <- oo[j]
# are (i, j) and (j, i) equivalent?
if(!ordered) {
ok <- (i < j)
i <- i[ok]
j <- j[ok]
xi <- xi[ok]
yi <- yi[ok]
xj <- xj[ok]
yj <- yj[ok]
dx <- dx[ok]
dy <- dy[ok]
d <- d[ok]
}
# done
return(list(i=i,
j=j,
xi=xi,
yi=yi,
xj=xj,
yj=yj,
dx=dx,
dy=dy,
d=d))
}
#######################
crosspairs <- function(X, Y, rmax) {
verifyclass(X, "ppp")
verifyclass(Y, "ppp")
stopifnot(is.numeric(rmax) && length(rmax) == 1 && rmax >= 0)
# order patterns by increasing x coordinate
ooX <- fave.order(X$x)
Xsort <- X[ooX]
ooY <- fave.order(Y$x)
Ysort <- Y[ooY]
if(spatstat.options("crosspairs.newcode")) {
# ------------------- use new faster code ---------------------
# First (over)estimate the number of pairs
nX <- npoints(X)
nY <- npoints(Y)
nsize <- ceiling(2 * pi * (rmax^2) * nX * nY/area.owin(Y$window))
nsize <- max(1024, nsize)
# .Call
Xx <- Xsort$x
Xy <- Xsort$y
Yx <- Ysort$x
Yy <- Ysort$y
r <- rmax
ng <- nsize
storage.mode(Xx) <- storage.mode(Xy) <- "double"
storage.mode(Yx) <- storage.mode(Yy) <- "double"
storage.mode(r) <- "double"
storage.mode(ng) <- "integer"
z <- .Call("Vcrosspairs",
xx1=Xx, yy1=Xy,
xx2=Yx, yy2=Yy,
rr=r, nguess=ng,
PACKAGE="spatstat")
if(length(z) != 9)
stop("Internal error: incorrect format returned from Vclosepairs")
i <- z[[1]] # NB no increment required
j <- z[[2]]
xi <- z[[3]]
yi <- z[[4]]
xj <- z[[5]]
yj <- z[[6]]
dx <- z[[7]]
dy <- z[[8]]
d <- z[[9]]
} else {
# Older code
# obtain upper estimate of number of pairs
# (to work around gcc bug 323)
DUP <- spatstat.options("dupC")
rmaxplus <- 1.25 * rmax
nsize <- .C("crosscount",
nn1=as.integer(X$n),
x1=as.double(Xsort$x),
y1=as.double(Xsort$y),
nn2=as.integer(Ysort$n),
x2=as.double(Ysort$x),
y2=as.double(Ysort$y),
rmaxi=as.double(rmaxplus),
count=as.integer(integer(1)),
DUP=DUP,
PACKAGE="spatstat")$count
if(nsize <= 0)
return(list(i=integer(0),
j=integer(0),
xi=numeric(0),
yi=numeric(0),
xj=numeric(0),
yj=numeric(0),
dx=numeric(0),
dy=numeric(0),
d=numeric(0)))
# allow slightly more space to work around gcc bug #323
nsize <- ceiling(1.1 * nsize) + X$n + Y$n
# now extract pairs
z <-
.C("crosspairs",
nn1=as.integer(X$n),
x1=as.double(Xsort$x),
y1=as.double(Xsort$y),
nn2=as.integer(Y$n),
x2=as.double(Ysort$x),
y2=as.double(Ysort$y),
r=as.double(rmax),
noutmax=as.integer(nsize),
nout=as.integer(integer(1)),
iout=as.integer(integer(nsize)),
jout=as.integer(integer(nsize)),
xiout=as.double(numeric(nsize)),
yiout=as.double(numeric(nsize)),
xjout=as.double(numeric(nsize)),
yjout=as.double(numeric(nsize)),
dxout=as.double(numeric(nsize)),
dyout=as.double(numeric(nsize)),
dout=as.double(numeric(nsize)),
status=as.integer(integer(1)),
DUP=DUP,
PACKAGE="spatstat")
if(z$status != 0)
stop(paste("Internal error: C routine complains that insufficient space was allocated:", nsize))
# trim vectors to the length indicated
npairs <- z$nout
if(npairs <= 0)
return(list(i=integer(0),
j=integer(0),
xi=numeric(0),
yi=numeric(0),
xj=numeric(0),
yj=numeric(0),
dx=numeric(0),
dy=numeric(0),
d=numeric(0)))
actual <- seq_len(npairs)
i <- z$iout[actual] # sic
j <- z$jout[actual]
xi <- z$xiout[actual]
yi <- z$yiout[actual]
xj <- z$xjout[actual]
yj <- z$yjout[actual]
dx <- z$dxout[actual]
dy <- z$dyout[actual]
d <- z$dout[actual]
}
# convert i,j indices to original sequences
i <- ooX[i]
j <- ooY[j]
# done
return(list(i=i,
j=j,
xi=xi,
yi=yi,
xj=xj,
yj=yj,
dx=dx,
dy=dy,
d=d))
}