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
|
#'
#' hotbox.R
#'
#' Heat kernel for a one-dimensional rod
#' and two-dimensional rectangle
#'
#' Code from Greg McSwiggan and Adrian Baddeley
hotbox <- function(Xsource, Xquery, sigma,
..., W=NULL, squared=FALSE, nmax=20) {
#' heat kernel in a rectangle
check.1.real(sigma)
if(is.null(W)) {
if(is.ppp(Xsource)) W <- Window(Xsource) else
if(is.sob(Xquery)) W <- Window(Xquery) else
stop("No window information is present")
} else {
stopifnot(is.owin(W))
if(!is.sob(Xsource)) Xsource <- as.ppp(Xsource, W)
if(!is.sob(Xquery)) Xquery <- as.ppp(Xquery, W)
}
if(!is.rectangle(W))
stop("The window must be a rectangle")
slen <- sidelengths(W)
Xsource <- shift(Xsource, origin="bottomleft")
Xquery <- shift(Xquery, origin="bottomleft")
nsource <- npoints(Xsource)
if(is.ppp(Xquery)) {
nquery <- npoints(Xquery)
answer <- numeric(nquery)
for(i in seq_len(nsource)) {
cx <- hotrod(slen[1], Xsource$x[i], Xquery$x, sigma,
ends="insulated", nmax=nmax)
cy <- hotrod(slen[2], Xsource$y[i], Xquery$y, sigma,
ends="insulated", nmax=nmax)
contrib <- cx * cy
if(squared) contrib <- contrib^2
answer <- answer + contrib
}
} else if(is.im(Xquery) || is.owin(Xquery)) {
Xquery <- as.im(Xquery, ...)
if(anyNA(Xquery)) stop("Image must be a full rectangle")
ansmat <- matrix(0, nrow(Xquery), ncol(Xquery))
xx <- Xquery$xcol
yy <- Xquery$yrow
for(i in seq_len(nsource)) {
cx <- hotrod(slen[1], Xsource$x[i], xx, sigma,
ends="insulated", nmax=nmax)
cy <- hotrod(slen[2], Xsource$y[i], yy, sigma,
ends="insulated", nmax=nmax)
contrib <- outer(cy, cx, "*")
if(squared) contrib <- contrib^2
ansmat <- ansmat + contrib
}
answer <- Xquery
answer[] <- ansmat
} else stop("Unrecognised format for Xquery")
return(answer)
}
|
