-
Notifications
You must be signed in to change notification settings - Fork 0
/
Kest_directional.R
142 lines (130 loc) · 4.03 KB
/
Kest_directional.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
#' Directed K function
#'
#' We use translation edge correction if we have an axis-oriented bounding box. Otherwise minus-correction.
#'
#' @param x pp, spatstat's ppp-object, or a coordinate matrix, or a list with $x~coordinates $bbox~bounding box
#' @param u unit direction(s), if many then one direction per row.
#' @param epsilon Central angle for the directed cone (total angle is 2 epsilon)
#' @param r radius vector at which to evaluate K
#' @param ... Passed on to \code{\link{Kest_anin}}
#' @param cylindrical If TRUE, compute the cylindrical version using \code{\link{Kest_anin_cylinder}}
#'
#' @details
#'
#' Compute the sector/cone/cylindrical K function. This version uses the more general anisotropic-inhomonogeneous \link{Kest_anin} and \code{\link{Kest_anin_cylinder}}, by setting the intensity = constant.
#'
#' @return
#' Returns a dataframe.
#'
#' @examples
#'
#' x <- matrix(runif(200), ncol=2)
#' k <- Kest_directional(x)
#' plot(k, rmax = 0.1)
#'
#' @useDynLib Kdirectional
#' @export
Kest_directional <- function(x, u, epsilon, r, ..., cylindrical = FALSE) {
x <- check_pp(x)
bbox <- x$bbox
n <- nrow(x$x)
lambda1 <- n/bbox_volume(bbox)
lambda <- rep(lambda1, n)
Kest_f <- if(cylindrical) Kest_anin_cylinder else Kest_anin
Kest <- Kest_f(x, u, epsilon, r, lambda, ...)
if(cylindrical) attr(Kest, "fun_name") <- "Cylindrical K"
else attr(Kest, "fun_name") <- "Conical K"
Kest
}
#' Directed K function, old and slow
#'
#' We use translation edge correction.
#'
#' @param x pp, list with $x~coordinates $bbox~bounding box
#' @param u unit vector of direction
#' @param epsilon Central angle for the directed cone (total angle is 2 epsilon)
#' @param r radius vector at which to evaluate K
#' @param pregraph A super graph for determining neighbourhoods
#'
#' @details
#'
#'
#' Compute the sector K function.
#'
#' Try \link{Kest_anin} for a bit better version that does also second order inhomogeneous patterns.
#'
#' @return
#' Returns a dataframe.
#'
#' @useDynLib Kdirectional
#' @export
Kest_directional_old <- function(x, u, epsilon, r, pregraph) {
x <- check_pp(x)
bbox <- x$bbox
if(is.bbquad(bbox)) stop("bbquad window not yet supported.")
dim <- ncol(bbox)
sidelengths <- apply(bbox, 2, diff)
if(missing(u)) stop("u needed")
if(missing(epsilon)) stop("epsilon needed")
if(abs(epsilon)>pi/2) stop("epsilon should be in range [0, pi/2]")
if(missing(r)) {
b <- min(sidelengths)*0.3
r <- seq(0, b, length=50)
}
# make sure unit vector
u <- rbind(u)
u <- t(apply(u, 1, function(ui) ui/c(sqrt(t(ui)%*%ui) )))
# start:
xc <- as.matrix(x$x)
from <- 1:nrow(xc)
to <- 1:nrow(xc)
Nr <- length(r)
N <- nrow(xc)
# lambda:
lambda <- N/prod(apply(bbox, 2, diff))
# the largest graph:
edges <- directed_geom(xc, u, epsilon, r[Nr], from, to, TRUE, pregraph)
#
# compute weights for translation correction
W <- translation_weights(x, asMatrix=TRUE)
# the K value
Kv <- function(edges) {
has_neighs <- which(sapply(edges, length)>0)
S <- if(length(has_neighs))
sum(sapply(has_neighs,
function(i) sum( 1/W[i, edges[[i]] ] )
)
) else 0
S/lambda^2
}
#
res <- c(Kv(edges), r[Nr])
for(i in (Nr-1):1) {
edges <- c_cutgeom(xc, edges, r[i])
res <- rbind(c(Kv(edges), r[i]), res)
}
#
bd <- pi * if(dim==3) 4/3 else 1
#
list(trans=res[,1], r=r, theo=bd*r^dim)
}
#' Compute all directed K functions along axis
#'
#' @param x pp, list with $x~coordinates $bbox~bounding box
#' @param ... parameters except u for Kest_directional
#' @export
Kest_along_axis <- function(x, ...) {
if(is.null(x$bbox)) stop("x should be list(x=coordinates-matrix, bbox=bounding-box)")
bbox <- x$bbox
dim <- ncol(x$x)
Kvals <- NULL
vecs <- rbind(c(1,0,0), c(0,1,0), c(0,0,1))
vecs <- vecs[1:dim, 1:dim]
# compute Kdir in different directions along axis
for(i in 1:dim){
kest <- Kest_directional(x, u=vecs[i,], ...)
Kvals <- rbind(Kvals, kest$trans)
}
r <- kest$r
list(r=r, vals=Kvals, dim=dim, theo=kest$theo)
}