#
# linnet.R
#
# Linear networks
#
# $Revision: 1.89 $ $Date: 2024/05/01 06:26:34 $
#
# An object of class 'linnet' defines a linear network.
# It includes the following components
#
# vertices (ppp) vertices of network
#
# m (matrix) adjacency matrix
#
# lines (psp) edges of network
#
# dpath (matrix) matrix of shortest path distances
# between each pair of vertices
#
# from, to (vectors) map from edges to vertices.
# The endpoints of the i-th segment lines[i]
# are vertices[from[i]] and vertices[to[i]]
#
#
# FUNCTIONS PROVIDED:
# linnet creates an object of class "linnet" from data
# print.linnet print an object of class "linnet"
# plot.linnet plot an object of class "linnet"
#
# Make an object of class "linnet" from the minimal data
<- (vertices, m, edges, sparse=, warn=) {
((m) && (edges))
("specify either m or edges")
(!(m) && !(edges))
("do not specify both m and edges")
# validate inputs
(is.ppp(vertices))
nv <- npoints(vertices)
(nv <= 1) {
m <- (, nv, nv)
from <- to <- (0)
} (!(m)) {
# check logical matrix or logical sparse matrix
(!(m) && !(m, ("lgCMatrix", "lgTMatrix")))
("m should be a matrix or sparse matrix")
((m) && (m))
((m) != vertices$n)
("dimensions of matrix m do not match number of vertices")
(((m))) {
("diagonal entries of the matrix m should not be TRUE; ignored")
(m) <-
}
sparse <- !(m)
## determine 'from' and 'to' vectors
ij <- (m, arr.ind=)
ij <- ij ij,1L] < ij,2L], , =]
from <- ij,1L]
to <- ij,2L]
} {
## check (from, to) pairs
((edges) && (edges) == 2)
(((edges %% 1) != 0))
("Entries of edges list should be integers")
((self <- (edges,1L] == edges,2L]))) {
("edge list should not join a vertex to itself; ignored")
edges <- edges!self, , =]
}
np <- npoints(vertices)
((edges > np))
("index out-of-bounds in edges list")
from <- edges,1L]
to <- edges,2L]
## avoid duplication in either sense
up <- (from < to)
ee <- ((up, from , to), (up, to, from))
((ee)) {
("Duplicated segments were ignored", call.=)
ok <- !(ee)
from <- from[ok]
to <- to[ok]
}
## convert to adjacency matrix
(!sparse) {
m <- (, np, np)
m[edges] <-
}
m <- sparseMatrix(i=from, j=to, x=, dims=(np, np))
m <- m | (m)
}
# create line segments
xx <- vertices$x
yy <- vertices$y
<- psp(xx[from], yy[from], xx[to], yy[to], =vertices$,
=)
# tolerance
toler <- ()
## pack up
out <- (vertices=vertices, m=m, =, from=from, to=to,
sparse=sparse, =vertices$,
toler=toler)
(out) <- ("linnet", (out))
## finish ?
(sparse)
(out)
# compute matrix of distances between adjacent vertices
n <- (m)
d <- (, n, n)
(d) <- 0
d[m] <- pairdist(vertices)[m]
## now compute shortest-path distances between each pair of vertices
out$dpath <- dpath <- dist2dpath(d)
(warn && ((dpath)))
("Network is not connected", call.=)
# pre-compute bounding radius
out$boundingradius <- boundingradius(out)
(out)
}
<- (x, of=("segments", "vertices"), ) {
of <- (of)
m <- (of,
vertices = marks(x$vertices, ),
= marks(x$, ))
(m)
}
"marks<-.linnet" <- (x, of=("segments", "vertices"), , value) {
of <- (of)
(of,
vertices = {
marks(x$vertices, ) <- value
},
= {
marks(x$, ) <- value
})
(x)
}
<- (x, ) {
nv <- x$vertices$n
nl <- x$$n
splat("Linear network with",
nv, (nv, "vertex", "vertices"),
"and",
nl, (nl, "line", "lines"))
(!(br <- x$boundingradius) && (br))
splat("[Network is not connected]")
(!(x$vertices$marks)) splat("Vertices are marked")
(!(x$$marks)) splat("Line segments are marked")
(as.owin(x), prefix="Enclosing window: ")
(())
}
<- (object, ) {
deg <- (object)
sparse <- object$sparse %orifnull% (object$dpath)
result <- (nvert = object$vertices$n,
nline = object$$n,
nedge = (deg)/2,
vertmarks = !(object$vertices$marks),
segmarks = !(object$$marks),
unitinfo = (unitname(object)),
totlength = (lengths_psp(object$)),
maxdegree = (deg),
ncomponents = ((connected(object, what="labels"))),
win = as.owin(object),
sparse = sparse)
(!sparse) {
result$diam <- diameter(object)
result$boundrad <- boundingradius(object)
}
result$toler <- object$toler
(result) <- ("summary.linnet", (result))
result
}
<- (x, ) {
dig <- ('digits')
(x, {
splat("Linear network with",
nvert, (nvert, "vertex", "vertices"),
"and",
nline, (nline, "line", "lines"))
(vertmarks) splat("Vertices are marked")
(segmarks) splat("Line segments are marked")
splat("Total length", (totlength, dig),
unitinfo$plural, unitinfo$explain)
splat("Maximum vertex degree:", maxdegree)
(sparse) splat("[Sparse matrix representation]")
splat("[Non-sparse matrix representation]")
(ncomponents > 1) {
splat("Network is disconnected: ", ncomponents, "connected components")
} {
splat("Network is connected")
(!sparse) {
splat("Diameter:", (diam, dig), unitinfo$plural)
splat("Bounding radius:", (boundrad, dig), unitinfo$plural)
}
}
(!(x$toler))
splat("Numerical tolerance:", (x$toler, dig), unitinfo$plural)
(win, prefix="Enclosing window: ")
})
(())
}
<- (x, , main=, add=,
do.plot=,
show.vertices=, show.window=,
args.vertices=(), args.segments=()) {
((main))
main <- short.deparse((x))
((x, "linnet"))
argh <- ()
<- as.psp(x)
(show.vertices) vert <- vertices(x)
#' plan layout and save symbolmaps
RS <- (,
resolve.defaults((x=(),
do.plot=,
show.window=show.window,
main=main,
multiplot=),
args.segments,
(),
(use.marks=, =)))
B <- as.owin(RS)
(show.vertices) {
RV <- (,
resolve.defaults((x=(vert),
do.plot=),
args.vertices,
(),
(use.marks=, =)))
BV <- as.owin(RV)
B <- boundingbox(B, BV)
} {
RV <-
}
## initialise plot
(!add) {
args.main <- argh[names(argh) %in% ("cex.main", "adj.main", "col")]
(, ((x=(B), ="n", main=main),
args.main))
}
## plot segments and (optionally) vertices
(,
resolve.defaults((x=(),
add=,
show.window=show.window,
show.all=,
main="",
multiplot=),
args.segments,
(),
(use.marks=, =)))
(show.vertices) {
(,
resolve.defaults((x=(vert),
add=,
show.window=,
show.all=,
main=""),
args.vertices,
(),
(use.marks=, =)))
}
result <- (=RS, vertices=RV)
(result, "bbox") <- B
((result))
}
<- (x, , fatal=) {
verifyclass(x, "linnet", fatal=fatal)
(x$)
}
<- (w) {
verifyclass(w, "linnet")
(w$vertices)
}
<- (x, ) {
verifyclass(x, "linnet")
(x$vertices$n)
}
<- (x) {
(x$$n)
}
<- (X, ) {
(X$)
}
"Window<-.linnet" <- (X, , =, value) {
() {
X <- X[value]
} {
X$ <- value
X$$ <- value
X$vertices$ <- value
}
(X)
}
<- (W, ) {
(Window(W))
}
<- (X, ) {
("as.linnet")
}
<- (X, , sparse, maxsize=30000) {
((sparse)) (X)
((X$sparse)) X$sparse <- (X$dpath)
(sparse && !(X$sparse)) {
# delete distance matrix
X$dpath <-
# convert adjacency matrix to sparse matrix
X$m <- (X$m, "sparseMatrix")
X$sparse <-
} (!sparse && X$sparse) {
# convert adjacency matrix to path-distance matrix
nv <- nvertices(X)
(nv > maxsize) {
(("Unable to create a matrix of size", nv, "x", nv,
paren(("max permitted size", maxsize, "x", maxsize))),
call.=)
}
X$m <- m <- (X$m)
edges <- (m, arr.ind=)
from <- edges,1L]
to <- edges,2L]
# compute distances to one-step neighbours
n <- (m)
d <- (, n, n)
(d) <- 0
coo <- coords(vertices(X))
d[edges] <- (((coo[from, 1:2] - coo[to, 1:2])^2))
# compute shortest path distance matrix
X$dpath <- dist2dpath(d)
# compute bounding radius
X$boundingradius <- boundingradius(X)
X$sparse <-
} (!sparse) {
# possibly update internals
X$boundingradius <- boundingradius(X)
}
# possibly update internals
X$circumradius <-
X$toler <- (X)
(X)
}
<- (X, , eps, sparse=) {
X <- selfcut.psp(X)
camefrom <- (X, "camefrom")
V <- (endpoints.psp(X))
((eps) || (eps)) {
eps <- ($double.eps) * diameter(Frame(X))
} {
check.1.real(eps)
(eps >= 0)
}
(eps > 0 && minnndist(V) <= eps) {
gV <- marks(connected(V, eps))
xx <- ((V$x, gV, ))
yy <- ((V$y, gV, ))
V <- ppp(xx, yy, =Window(X))
}
first <- endpoints.psp(X, "first")
second <- endpoints.psp(X, "second")
from <- nncross(first, V, what="which")
to <- nncross(second, V, what="which")
((reverse <- (from > to))) {
newfrom <- (reverse, to, from)
newto <- (reverse, from, to)
from <- newfrom
to <- newto
}
fromto <- (from, to)
nontrivial <- (from != to) & !(fromto)
join <- fromto[nontrivial, , =]
result <- (V, edges=join, sparse=sparse)
(is.marked(X)) marks(result$) <- marks(X[nontrivial])
(result, "camefrom") <- camefrom[nontrivial]
(result)
}
<- (x) {
unitname(x$)
}
"unitname<-.linnet" <- (x, value) {
w <- x$
v <- x$vertices
l <- x$
unitname(w) <- unitname(v) <- unitname(l) <- value
x$ <- w
x$vertices <- v
x$ <- l
(x)
}
<- (x) {
((x, "linnet"))
dpath <- x$dpath
((dpath)) () ((0, dpath))
}
<- (x) {
(lengths_psp(x$))
}
<- (x) {
verifyclass(x, "linnet")
((x$m))
}
<- (x, ) {
("boundingradius.linnet")
(x, )
}
<- (x, ) {
((x, "linnet"))
cr <- x$boundingradius %orifnull% x$circumradius
(!(cr))
(cr)
dpath <- x$dpath
((dpath)) ()
(((dpath))) ()
((dpath) <= 1)
((0,dpath))
from <- x$from
to <- x$to
<- x$
nseg <- $n
leng <- lengths_psp()
(spatstat.options("Clinearradius")) {
fromC <- from - 1L
toC <- to - 1L
nv <- npoints(vertices(x))
huge <- (leng)
z <- (SL_linearradius,
= (nseg),
from = (fromC),
to = (toC),
= (leng),
nv = (nv),
dpath = (dpath),
huge = (huge),
result = ((1)),
PACKAGE="spatstat.linnet")
(z$result)
}
sA <- sB <- (, nseg, nseg)
(i 1:nseg) {
# endpoints of segment i
A <- from[i]
B <- to[i]
AB <- leng[i]
sA[i,i] <- sB[i,i] <- AB/2
(j in (1:nseg)-i]) {
# endpoints of segment j
<- from[j]
<- to[j]
CD <- leng[j]
AC <- dpath[A,]
AD <- dpath[A,]
BC <- dpath[B,]
BD <- dpath[B,]
# max dist from A to any point in segment j
sA[i,j] <- (AD > AC + CD) AC + CD
(AC > AD + CD) AD + CD else
(AC + AD + CD)/2
# max dist from B to any point in segment j
sB[i,j] <- (BD > BC + CD) BC + CD
(BC > BD + CD) BD + CD else
(BC + BD + CD)/2
}
}
# max dist from each A to any point in another segment
mA <- (sA, 1, )
# max dist from each B to any point in another segment
mB <- (sB, 1, )
# min of these
(mA, mB)
}
####################################################
# affine transformations
####################################################
<- (X, f, ) {
trap.extra.arguments(, .Context="In scalardilate(X,f)")
check.1.real(f, "In scalardilate(X,f)")
((f) && f > 0)
Y <- X
Y$vertices <- scalardilate(X$vertices, f=f)
Y$ <- scalardilate(X$, f=f)
Y$ <- scalardilate(X$, f=f)
(!(X$dpath)) {
Y$dpath <- f * X$dpath
Y$boundingradius <- f * (X$boundingradius %orifnull% X$circumradius)
Y$circumradius <-
}
(!(X$toler))
X$toler <- (f * X$toler)
(Y)
}
<- (X, mat=((1,1)), vec=(0,0), ) {
verifyclass(X, "linnet")
((((mat)$values)) == 1) {
# transformation is an isometry
scal <- (((mat)))
Y <- X
Y$vertices <- affine(X$vertices, mat=mat, vec=vec, )
Y$ <- affine(X$, mat=mat, vec=vec, )
Y$ <- affine(X$, mat=mat, vec=vec, )
(!(X$dpath)) {
Y$dpath <- scal * X$dpath
Y$boundingradius <- scal * (X$boundingradius %orifnull% X$circumradius)
X$circumradius <-
}
(!(Y$toler))
Y$toler <- (scal * Y$toler)
} {
# general case
vertices <- affine(X$vertices, mat=mat, vec=vec, )
Y <- (vertices, edges=(X$from, X$to))
}
(Y)
}
<- (X, vec=(0,0), , origin=) {
verifyclass(X, "linnet")
Y <- X
(!(origin)) {
(!(vec))
("argument vec ignored; argument origin has precedence")
locn <- interpretAsOrigin(origin, Window(X))
vec <- -locn
}
Y$ <- W <- shift(X$, vec=vec, )
v <- getlastshift(W)
Y$vertices <- shift(X$vertices, vec=v, )
Y$ <- shift(X$, vec=v, )
# tack on shift vector
(Y, "lastshift") <- v
(Y)
}
<- (X, angle=/2, , centre=) {
verifyclass(X, "linnet")
(!(centre)) {
X <- shift(X, origin=centre)
negorigin <- getlastshift(X)
} negorigin <-
Y <- X
Y$vertices <- rotate(X$vertices, angle=angle, )
Y$ <- rotate(X$, angle=angle, )
Y$ <- rotate(X$, angle=angle, )
(!(negorigin))
Y <- shift(Y, -negorigin)
(Y)
}
<- (X, s, unitname) {
((unitname)) unitname <-
((s) || (s)) s <- 1/unitname(X)$multiplier
Y <- scalardilate(X, f=1/s)
unitname(Y) <- rescale(unitname(X), s, unitname)
(Y)
}
"[.linnet" <- (x, i, , snip=) {
(!is.owin(i))
("In [.linnet: the index i should be a window", call.=)
x <- (x)
w <- i
wp <- as.polygonal(w)
(is.mask(w)) {
## protect against pixellation artefacts
pixel <- owinInternalRect(w$xstep * (-1,1)/2, w$ystep * (-1,1)/2)
wp <- MinkowskiSum(wp, pixel)
wp <- intersect.owin(wp, Frame(w))
}
## Find vertices that lie inside window
vertinside <- inside.owin(x$vertices, w=wp)
from <- x$from
to <- x$to
(snip) {
(!(marks(vertices(x)))) {
(("Marks attached to the network vertices were removed,",
"because no mark information is available",
"for the new vertices created by [.linnet when snip=TRUE"),
call.=)
marks(x$vertices) <-
}
## For efficiency, first restrict network to relevant segments.
## Find segments EITHER OF whose endpoints lie in 'w' ...
okedge <- vertinside[from] | vertinside[to]
## ... or which cross the boundary of 'w'
xlines <- as.psp(x)
wlines <- edges(wp)
((miss <- !okedge)) {
hits <- (test.crossing.psp(xlines[miss], wlines), 1, )
okedge[miss] <- hits
}
## extract relevant subset of network graph
x <- (x, retainedges=okedge)
## Now add vertices at crossing points with boundary of 'w'
b <- (crossing.psp(xlines, wlines))
novel <- (nncross(b, x$vertices, what="dist") > 0)
x <- (x, b[novel])
boundarypoints <- (x, "id")
## update data
from <- x$from
to <- x$to
vertinside <- inside.owin(x$vertices, w=wp)
vertinside[boundarypoints] <-
}
## find segments whose endpoints BOTH lie in 'w'
edgeinside <- vertinside[from] & vertinside[to]
## .. and which are not trivial
umap <- uniquemap(x$vertices)
nontrivial <- (umap[from] != umap[to])
## extract relevant subset of network
xnew <- (x, retainedges=edgeinside & nontrivial)
## adjust window efficiently
Window(xnew, =) <- w
(xnew)
}
<- (x, ) {
pixellate(as.psp(x), )
}
<- (X, , what=("labels", "components")) {
verifyclass(X, "linnet")
what <- (what)
nv <- npoints(vertices(X))
lab0 <- cocoEngine(nv, X$from - 1L, X$to - 1L, "connected.linnet")
lab <- lab0 + 1L
lab <- (((lab)))
(what == "labels")
(lab)
nets <- ()
subsets <- ((nv), lab)
(i (subsets))
nets[i]] <- (X, retainvertices=subsets[i]])
(nets)
}
<- (X, ) {
(!(dpath <- X$dpath))
(((dpath)))
lab <- connected(X, what="labels")
npieces <- ((lab))
(npieces == 1)
}
<- (X, Y) {
X <- (X)
(!(Y, ("linnet", "infline", "psp")))
("L should be an object of class psp, linnet or infline", call.=)
## convert infinite lines to segments
((Y, "linnet")) Y <- as.psp(Y)
((Y, "infline")) {
Y <- clip.infline(Y, Frame(X))
id <- marks(Y)
lev <- (id)
} {
id <- lev <- (nsegments(Y))
}
## extract segments of network
S <- as.psp(X)
## find crossing points
SY <- crossing.psp(S, Y, fatal=, details=)
((SY) || npoints(SY) == 0)
((L=X))
SY <- (SY)
Z <- ((SY),
(x=x, y=y, seg=iA, tp=tA, L=X,
marks=(id[as.integer(jB)], =lev)))
(Z)
}
<- (x, ) {
density.psp(as.psp(x), )
}
<- (L) {
## identify terminal vertices and return them as an lpp object
verifyclass(L, "linnet")
ind <- ((L) == 1)
nV <- (ind)
(nV == 0) {
## return empty pattern
((L=L))
}
V <- vertices(L)[ind]
## construct network coordinates
seg <- (, nV)
tp <- (, nV)
## look for vertices mentioned as the start of a segment
left <- (ind, L$from)
found <- !(left)
seg[found] <- left[found]
tp[found] <- 0
## look for vertices mentioned as the end of a segment
right <- (ind, L$to)
found <- !(right)
seg[found] <- right[found]
tp[found] <- 1
## pack up
<- (coords(V), (seg=seg, tp=tp))
B <- (, L)
(B)
}
