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
139
140
141
142
143
144
145
146
147
148
149
150
|
#
# edgeTrans.R
#
# $Revision: 1.16 $ $Date: 2019/01/18 02:26:41 $
#
# Translation edge correction weights
#
# edge.Trans(X) compute translation correction weights
# for each pair of points from point pattern X
#
# edge.Trans(X, Y, W) compute translation correction weights
# for all pairs of points X[i] and Y[j]
# (i.e. one point from X and one from Y)
# in window W
#
# edge.Trans(X, Y, W, paired=TRUE)
# compute translation correction weights
# for each corresponding pair X[i], Y[i].
#
# To estimate the K-function see the idiom in "Kest.R"
#
#######################################################################
edge.Trans <- function(X, Y=X, W=Window(X), exact=FALSE, paired=FALSE,
...,
trim=spatstat.options("maxedgewt"),
dx=NULL, dy=NULL,
give.rmax=FALSE,
gW = NULL) {
given.dxdy <- !is.null(dx) && !is.null(dy)
if(!given.dxdy) {
## dx, dy will be computed from X, Y
X <- as.ppp(X, W)
W <- X$window
Y <- if(!missing(Y)) as.ppp(Y, W) else X
nX <- X$n
nY <- Y$n
if(paired) {
if(nX != nY)
stop("X and Y should have equal length when paired=TRUE")
dx <- Y$x - X$x
dy <- Y$y - X$y
} else {
dx <- outer(X$x, Y$x, "-")
dy <- outer(X$y, Y$y, "-")
}
} else {
## dx, dy given
if(paired) {
## dx, dy are vectors
check.nvector(dx, vname="dx")
check.nvector(dy, vname="dy")
stopifnot(length(dx) == length(dy))
} else {
## dx, dy are matrices
check.nmatrix(dx)
check.nmatrix(dy)
stopifnot(all(dim(dx) == dim(dy)))
nX <- nrow(dx)
nY <- ncol(dx)
}
stopifnot(is.owin(W))
}
## For irregular polygons, exact evaluation is very slow;
## so use pixel approximation, unless exact=TRUE
if(W$type == "polygonal" && !exact)
W <- as.mask(W)
## compute
if(!paired) {
dx <- as.vector(dx)
dy <- as.vector(dy)
}
switch(W$type,
rectangle={
## Fast code for this case
wide <- diff(W$xrange)
high <- diff(W$yrange)
weight <- wide * high / ((wide - abs(dx)) * (high - abs(dy)))
},
polygonal={
## This code is SLOW
n <- length(dx)
weight <- numeric(n)
if(n > 0) {
for(i in seq_len(n)) {
Wshift <- shift(W, c(dx[i], dy[i]))
weight[i] <- overlap.owin(W, Wshift)
}
weight <- area(W)/weight
}
},
mask={
## compute set covariance of window
if(is.null(gW)) gW <- setcov(W)
## evaluate set covariance at these vectors
gvalues <- lookup.im(gW, dx, dy, naok=TRUE, strict=FALSE)
weight <- area(W)/gvalues
}
)
## clip high values
if(length(weight) > 0)
weight <- pmin.int(weight, trim)
if(!paired)
weight <- matrix(weight, nrow=nX, ncol=nY)
if(give.rmax)
attr(weight, "rmax") <- rmax.Trans(W, gW)
return(weight)
}
## maximum radius for translation correction
## = radius of largest circle centred at 0 contained in W + ^W
rmax.Trans <- function(W, g=setcov(W)) {
## calculate maximum permissible 'r' value
## for validity of translation correction
W <- as.owin(W)
if(is.rectangle(W))
return(shortside(W))
## find support of set covariance
if(is.null(g)) g <- setcov(W)
eps <- 2 * max(1, max(g)) * .Machine$double.eps
gsupport <- solutionset(g > eps)
gboundary <- bdry.mask(gsupport)
xy <- rasterxy.mask(gboundary, drop=TRUE)
rmax <- with(xy, sqrt(min(x^2 + y^2)))
return(rmax)
}
## maximum radius for rigid motion correction
## = radius of smallest circle centred at 0 containing W + ^W
rmax.Rigid <- function(X, g=setcov(as.owin(X))) {
stopifnot(is.ppp(X) || is.owin(X))
if(is.ppp(X))
return(max(pairdist(X[chull(X)])))
W <- X
if(is.rectangle(W)) return(diameter(W))
if(is.null(g)) g <- setcov(W)
eps <- 2 * max(1, max(g)) * .Machine$double.eps
gsupport <- solutionset(g > eps)
gboundary <- bdry.mask(gsupport)
xy <- rasterxy.mask(gboundary, drop=TRUE)
rmax <- with(xy, sqrt(max(x^2 + y^2)))
return(rmax)
}
|
