pairdistlpp.R
#
# pairdistlpp.R
#
# $Revision: 1.11 $ $Date: 2015/02/13 06:42:44 $
#
#
# pairdist.lpp
# Calculates the shortest-path distance between each pair of points
# in a point pattern on a linear network.
#
pairdist.lpp <- function(X, ..., method="C") {
stopifnot(inherits(X, "lpp"))
stopifnot(method %in% c("C", "interpreted"))
#
L <- as.linnet(X, sparse=FALSE)
Y <- as.ppp(X)
n <- npoints(Y)
#
# Lseg <- L$lines
Lvert <- L$vertices
from <- L$from
to <- L$to
dpath <- L$dpath
# nearest segment for each point
pro <- coords(X, local=TRUE, spatial=FALSE, temporal=FALSE)$seg
pairdistmat <- matrix(0,n,n)
if(method == "interpreted") {
# loop through all pairs of data points
for (i in 1:(n-1)) {
proi <- pro[i]
Xi <- Y[i]
nbi1 <- from[proi]
nbi2 <- to[proi]
vi1 <- Lvert[nbi1]
vi2 <- Lvert[nbi2]
dXi1 <- crossdist(Xi, vi1)
dXi2 <- crossdist(Xi, vi2)
for (j in (i+1):n) {
Xj <- Y[j]
proj <- pro[j]
if(proi == proj) {
# points i and j lie on the same segment
# use Euclidean distance
d <- crossdist(Xi, Xj)
} else {
# shortest path from i to j passes through ends of segments
nbj1 <- from[proj]
nbj2 <- to[proj]
vj1 <- Lvert[nbj1]
vj2 <- Lvert[nbj2]
# Calculate shortest of 4 possible paths from i to j
d1Xj <- crossdist(vj1,Xj)
d2Xj <- crossdist(vj2,Xj)
d11 <- dXi1 + dpath[nbi1,nbj1] + d1Xj
d12 <- dXi1 + dpath[nbi1,nbj2] + d2Xj
d21 <- dXi2 + dpath[nbi2,nbj1] + d1Xj
d22 <- dXi2 + dpath[nbi2,nbj2] + d2Xj
d <- min(d11,d12,d21,d22)
}
# store result
pairdistmat[i,j] <- pairdistmat[j,i] <- d
}
}
} else {
# C code
# convert indices to start at 0
from0 <- from - 1L
to0 <- to - 1L
segmap <- pro - 1L
zz <- .C("linpairdist",
np = as.integer(n),
xp = as.double(Y$x),
yp = as.double(Y$y),
nv = as.integer(Lvert$n),
xv = as.double(Lvert$x),
yv = as.double(Lvert$y),
ns = as.double(L$n),
from = as.integer(from0),
to = as.integer(to0),
dpath = as.double(dpath),
segmap = as.integer(segmap),
answer = as.double(pairdistmat))
pairdistmat <- matrix(zz$answer, n, n)
}
return(pairdistmat)
}