Revision 8718b3c4133e9eecaa99e6e95a040fcf22b80c30 authored by Adrian Baddeley on 01 February 2007, 16:09:56 UTC, committed by cran-robot on 01 February 2007, 16:09:56 UTC
1 parent 79b88fd
distances.R
#
# distances.R
#
# $Revision: 1.17 $ $Date: 2006/10/09 03:15:32 $
#
#
# Interpoint distances
#
#
pairdist <- function(X, ...) {
UseMethod("pairdist")
}
pairdist.ppp <- function(X, ..., periodic=FALSE, method="C") {
verifyclass(X, "ppp")
if(!periodic)
return(pairdist.default(X$x, X$y, method=method))
# periodic case
W <- X$window
if(W$type != "rectangle")
stop(paste("periodic edge correction can't be applied",
"in a non-rectangular window"))
wide <- diff(W$xrange)
high <- diff(W$yrange)
return(pairdist.default(X$x, X$y, period=c(wide,high), method=method))
}
pairdist.default <-
function(X, Y=NULL, ..., period=NULL, method="C")
{
xy <- xy.coords(X,Y)[c("x","y")]
x <- xy$x
y <- xy$y
n <- length(x)
if(length(y) != n)
stop("lengths of x and y do not match")
# special cases
if(n == 0)
return(matrix(numeric(0), nrow=0, ncol=0))
else if(n == 1)
return(matrix(0,nrow=1,ncol=1))
if((periodic<- !is.null(period))) {
stopifnot(is.numeric(period))
stopifnot(length(period) == 2 || length(period) == 1)
stopifnot(all(period > 0))
if(length(period) == 1) period <- rep(period, 2)
wide <- period[1]
high <- period[2]
}
switch(method,
interpreted={
xx <- matrix(rep(x, n), nrow = n)
yy <- matrix(rep(y, n), nrow = n)
if(!periodic) {
d <- sqrt((xx - t(xx))^2 + (yy - t(yy))^2)
} else {
dx <- xx - t(xx)
dy <- yy - t(yy)
dx2 <- pmin(dx^2, (dx + wide)^2, (dx - wide)^2)
dy2 <- pmin(dy^2, (dy + high)^2, (dy - high)^2)
d <- sqrt(dx2 + dy2)
}
},
C={
d <- numeric( n * n)
if(!periodic) {
z<- .C("pairdist", n = as.integer(n),
x= as.double(x), y= as.double(y), d= as.double(d),
PACKAGE="spatstat")
} else {
z<- .C("pairPdist", n = as.integer(n),
x= as.double(x), y= as.double(y),
xwidth=as.double(wide), yheight=as.double(high),
d= as.double(d),
PACKAGE="spatstat")
}
d <- matrix(z$d, nrow=n, ncol=n)
},
stop(paste("Unrecognised method", sQuote(method)))
)
invisible(d)
}
nndist <- function(X, ..., method="C") {
UseMethod("nndist")
}
nndist.ppp <- function(X, ..., method="C") {
verifyclass(X, "ppp")
return(nndist.default(X$x, X$y, method=method))
}
nndist.default <-
function(X, Y=NULL, ..., method="C")
{
# computes the vector of nearest-neighbour distances
# for the pattern of points (x[i],y[i])
#
xy <- xy.coords(X,Y)[c("x","y")]
x <- xy$x
y <- xy$y
# validate
n <- length(x)
if(length(y) != n)
stop("lengths of x and y do not match")
# special cases
if(n == 0)
return(numeric(0))
else if(n == 1)
return(Inf)
switch(method,
interpreted={
# matrix of squared distances between all pairs of points
sq <- function(a, b) { (a-b)^2 }
squd <- outer(x, x, sq) + outer(y, y, sq)
# reset diagonal to a large value so it is excluded from minimum
diag(squd) <- Inf
# nearest neighbour distances
nnd <- sqrt(apply(squd,1,min))
},
C={
n <- length(x)
nnd<-numeric(n)
o <- order(y)
big <- sqrt(.Machine$double.xmax)
z<- .C("nndistsort",
n= as.integer(n),
x= as.double(x[o]), y= as.double(y[o]), nnd= as.double(nnd),
as.double(big),
PACKAGE="spatstat")
nnd[o] <- z$nnd
},
stop(paste("Unrecognised method", sQuote(method)))
)
invisible(nnd)
}
nnwhich <- function(X, ..., method="C") {
UseMethod("nnwhich")
}
nnwhich.ppp <- function(X, ..., method="C") {
verifyclass(X, "ppp")
return(nnwhich.default(X$x, X$y, method=method))
}
nnwhich.default <-
function(X, Y=NULL, ..., method="C")
{
# identifies nearest neighbour of each point in
# the pattern of points (x[i],y[i])
#
xy <- xy.coords(X,Y)[c("x","y")]
x <- xy$x
y <- xy$y
# validate
n <- length(x)
if(length(y) != n)
stop("lengths of x and y do not match")
# special cases
if(n == 0)
return(numeric(0))
else if(n == 1)
return(NA)
switch(method,
interpreted={
# matrix of squared distances between all pairs of points
sq <- function(a, b) { (a-b)^2 }
squd <- outer(x, x, sq) + outer(y, y, sq)
# reset diagonal to a large value so it is excluded from minimum
diag(squd) <- Inf
# nearest neighbour distances
nnwhich <- sqrt(apply(squd,1,which.min))
},
C={
n <- length(x)
nnd<-numeric(n)
nnwhich<-integer(n)
o <- order(y)
big <- sqrt(.Machine$double.xmax)
z<- .C("nnwhichsort",
n= as.integer(n),
x= as.double(x[o]), y= as.double(y[o]),
nnd= as.double(nnd),
nnwhich=as.integer(nnwhich),
as.double(big),
PACKAGE="spatstat")
# convert from C to R indexing
witch <- z$nnwhich + 1
if(any(witch <= 0))
stop("Internal error: non-positive index returned from C code")
if(any(witch > n))
stop("Internal error: index returned from C code exceeds n")
nnwhich[o] <- o[witch]
},
stop(paste("Unrecognised method", sQuote(method)))
)
return(nnwhich)
}
crossdist <- function(X, Y, ...) {
UseMethod("crossdist")
}
crossdist.ppp <- function(X, Y, ..., periodic=FALSE, method="C") {
verifyclass(X, "ppp")
Y <- as.ppp(Y)
if(!periodic)
return(crossdist.default(X$x, X$y, Y$x, Y$y, method=method))
# periodic case
WX <- X$window
WY <- Y$window
if(WX$type != "rectangle" || WY$type != "rectangle")
stop(paste("periodic edge correction can't be applied",
"in non-rectangular windows"))
if(!is.subset.owin(WX,WY) || !is.subset.owin(WY, WX))
stop(paste("periodic edge correction is not implemented",
"for the case when X and Y lie in different rectangles"))
wide <- diff(WX$xrange)
high <- diff(WX$yrange)
return(crossdist.default(X$x, X$y, Y$x, Y$y,
period=c(wide,high), method=method))
}
crossdist.default <-
function(X, Y, x2, y2, ..., period=NULL, method="C")
{
x1 <- X
y1 <- Y
# returns matrix[i,j] = distance from (x1[i],y1[i]) to (x2[j],y2[j])
if(length(x1) != length(y1))
stop("lengths of x and y do not match")
if(length(x2) != length(y2))
stop("lengths of x2 and y2 do not match")
n1 <- length(x1)
n2 <- length(x2)
if(n1 == 0 || n2 == 0)
return(matrix(numeric(0), nrow=n1, ncol=n2))
if((periodic<- !is.null(period))) {
stopifnot(is.numeric(period))
stopifnot(length(period) == 2 || length(period) == 1)
stopifnot(all(period > 0))
if(length(period) == 1) period <- rep(period, 2)
wide <- period[1]
high <- period[2]
}
switch(method,
interpreted = {
X1 <- matrix(rep(x1, n2), ncol = n2)
Y1 <- matrix(rep(y1, n2), ncol = n2)
X2 <- matrix(rep(x2, n1), ncol = n1)
Y2 <- matrix(rep(y2, n1), ncol = n1)
if(!periodic)
d <- sqrt((X1 - t(X2))^2 + (Y1 - t(Y2))^2)
else {
dx <- X1 - t(X2)
dy <- Y1 - t(Y2)
dx2 <- pmin(dx^2, (dx + wide)^2, (dx - wide)^2)
dy2 <- pmin(dy^2, (dy + high)^2, (dy - high)^2)
d <- sqrt(dx2 + dy2)
}
return(d)
},
C = {
if(!periodic) {
z<- .C("crossdist",
nfrom = as.integer(n1),
xfrom = as.double(x1),
yfrom = as.double(y1),
nto = as.integer(n2),
xto = as.double(x2),
yto = as.double(y2),
d = as.double(matrix(0, nrow=n1, ncol=n2)),
PACKAGE="spatstat")
} else {
z<- .C("crossPdist",
nfrom = as.integer(n1),
xfrom = as.double(x1),
yfrom = as.double(y1),
nto = as.integer(n2),
xto = as.double(x2),
yto = as.double(y2),
xwidth = as.double(wide),
yheight = as.double(high),
d = as.double(matrix(0, nrow=n1, ncol=n2)),
PACKAGE="spatstat")
}
return(matrix(z$d, nrow=n1, ncol=n2))
},
stop(paste("Unrecognised method", method))
)
}
Computing file changes ...