Revision 56d752cf6d1a8e2775e1acc46482c206cfd5b269 authored by Adrian Baddeley on 14 March 2012, 07:09:51 UTC, committed by cran-robot on 14 March 2012, 07:09:51 UTC
1 parent a4d492a
deldir.R
#
# deldir.R
#
# Interface to deldir package
#
# $Revision: 1.10 $ $Date: 2011/07/31 10:08:13 $
#
.spst.triEnv <- new.env()
assign("use.trigraf", TRUE, envir=.spst.triEnv)
assign("use.trigrafS", TRUE, envir=.spst.triEnv)
assign("debug.delaunay", FALSE, envir=.spst.triEnv)
dirichlet <- function(X) {
w <- X$window
dd <- deldir(X$x, X$y, rw=c(w$xrange,w$yrange))
pp <- lapply(tile.list(dd), function(z) { owin(poly=z[c("x","y")]) })
dir <- tess(tiles=pp, window=as.rectangle(w))
if(w$type != "rectangle")
dir <- intersect.tess(dir, w)
return(dir)
}
delaunay <- function(X) {
stopifnot(is.ppp(X))
X <- unique(X)
nX <- npoints(X)
if(nX < 3) return(NULL)
w <- X$window
dd <- deldir(X$x, X$y, rw=c(w$xrange, w$yrange))
a <- dd$delsgs[,5]
b <- dd$delsgs[,6]
use.trigraf <- get("use.trigraf", envir=.spst.triEnv)
use.trigrafS <- get("use.trigrafS", envir=.spst.triEnv)
debug.delaunay <- get("debug.delaunay", envir=.spst.triEnv)
if(use.trigrafS) {
# first ensure a[] < b[]
swap <- (a > b)
if(any(swap)) {
oldb <- b
b[swap] <- a[swap]
a[swap] <- oldb[swap]
}
# next ensure a is sorted
o <- order(a, b)
a <- a[o]
b <- b[o]
#
nv <- nX
ne <- length(a)
ntmax <- ne
z <- .C("trigrafS",
nv = as.integer(nv),
ne = as.integer(ne),
ie = as.integer(a - 1),
je = as.integer(b - 1),
ntmax = as.integer(ntmax),
nt = as.integer(integer(1)),
it = as.integer(integer(ne)),
jt = as.integer(integer(ne)),
kt = as.integer(integer(ne)),
status = as.integer(integer(1)),
PACKAGE="spatstat")
if(z$status != 0)
stop("Internal error: overflow in trigrafS")
tlist <- with(z, cbind(it, jt, kt)[1:nt, ]) + 1
} else if(use.trigraf) {
nv <- nX
ne <- length(a)
ntmax <- ne
z <- .C("trigraf",
nv = as.integer(nv),
ne = as.integer(ne),
ie = as.integer(a - 1),
je = as.integer(b - 1),
ntmax = as.integer(ntmax),
nt = as.integer(integer(1)),
it = as.integer(integer(ntmax)),
jt = as.integer(integer(ntmax)),
kt = as.integer(integer(ntmax)),
status = as.integer(integer(1)),
PACKAGE="spatstat")
if(z$status != 0)
stop("Internal error: overflow in trigraf")
tlist <- with(z, cbind(it, jt, kt)[1:nt, ]) + 1
} else {
tlist <- matrix(integer(0), 0, 3)
for(i in seq_len(nX)) {
# find all Delaunay neighbours of i
jj <- c(b[a==i], a[b==i])
jj <- sort(unique(jj))
# select those with a higher index than i
jj <- jj[jj > i]
# find pairs of neighbours which are Delaunay neighbours
# (thus, triangles where the first numbered vertex is i)
if(length(jj) > 0)
for(j in jj) {
kk <- c(b[a == j], a[b == j])
kk <- kk[(kk %in% jj) & (kk > j)]
if(length(kk) > 0)
for(k in kk)
# add (i,j,k) to list of triangles (i < j < k)
tlist <- rbind(tlist, c(i, j, k))
}
}
}
# At this point, `tlist' contains all triangles formed by the Delaunay edges,
# with vertices given in ascending order i < j < k in the 3 columns of tlist.
# Some of these triangles may not belong to the Delaunay triangulation.
# They will be weeded out later.
# Assemble coordinates of triangles
x <- X$x
y <- X$y
xtri <- matrix(x[tlist], nrow(tlist), 3)
ytri <- matrix(y[tlist], nrow(tlist), 3)
# ensure triangle vertices are in anticlockwise order
ztri <- ytri - min(y)
dx <- cbind(xtri[,2]-xtri[,1], xtri[,3]-xtri[,2], xtri[,1]-xtri[,3])
zm <- cbind(ztri[,1]+ztri[,2], ztri[,2]+ztri[,3], ztri[,3]+ztri[,1])
negareas <- apply(dx * zm, 1, sum)
clockwise <- (negareas > 0)
#
if(any(clockwise)) {
xc <- xtri[clockwise, , drop=FALSE]
yc <- ytri[clockwise, , drop=FALSE]
tc <- tlist[clockwise, , drop=FALSE]
xtri[clockwise,] <- xc[,c(1,3,2)]
ytri[clockwise,] <- yc[,c(1,3,2)]
tlist[clockwise,] <- tc[, c(1,3,2)]
}
# At this point, triangle vertices are listed in anticlockwise order.
# The same directed edge (i, j) cannot appear twice.
# To weed out invalid triangles, check for such duplication
triedges <- rbind(tlist[, c(1,2)],
tlist[, c(2,3)],
tlist[, c(3,1)])
if(any(bad <- duplicated(triedges))) {
badedges <- unique(triedges[bad, , drop=FALSE])
ntri <- nrow(tlist)
triid <- rep(seq_len(ntri), 3)
illegal <- rep(FALSE, ntri)
for(j in seq_len(nrow(badedges))) {
from <- badedges[j, 1]
to <- badedges[j, 2]
if(debug.delaunay)
cat(paste("Suspect edge from vertex", from, "to vertex", to, "\n"))
# find all triangles sharing this edge in this orientation
sustri <- triid[(triedges[,1] == from) & (triedges[,2] == to)]
if(debug.delaunay)
cat(paste("\tInvestigating triangles", commasep(sustri), "\n"))
# list all vertices associated with the suspect triangles
susvert <- sort(unique(as.vector(tlist[sustri, ])))
if(debug.delaunay)
cat(paste("\tInvestigating vertices", commasep(susvert), "\n"))
xsusvert <- x[susvert]
ysusvert <- y[susvert]
# take each triangle in turn and check whether it contains a data point
for(k in sustri) {
if(!illegal[k] &&
any(inside.triangle(xsusvert, ysusvert, xtri[k,], ytri[k,]))) {
if(debug.delaunay)
cat(paste("Triangle", k, "is illegal\n"))
illegal[k] <- TRUE
}
}
}
if(!any(illegal)) {
if(debug.delaunay)
cat("No illegal triangles found\n")
} else {
if(debug.delaunay)
cat(paste("Removing", sum(illegal), "triangles\n"))
tlist <- tlist[!illegal, , drop=FALSE]
xtri <- xtri[!illegal, , drop=FALSE]
ytri <- ytri[!illegal, , drop=FALSE]
}
}
# make tile list
tiles <- list()
for(m in seq_len(nrow(tlist))) {
p <- list(x=xtri[m,], y=ytri[m,])
tiles[[m]] <- owin(poly=p, check=FALSE)
}
wc <- convexhull.xy(x, y)
del <- tess(tiles=tiles, window=wc)
if(w$type != "rectangle")
del <- intersect.tess(del, w)
return(del)
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...