Skip to content

Commit

Permalink
Diffusion smoothing
Browse files Browse the repository at this point in the history
  • Loading branch information
baddstats committed Oct 6, 2024
1 parent 32dc551 commit fe2ba51
Show file tree
Hide file tree
Showing 14 changed files with 744 additions and 5 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: spatstat.explore
Version: 3.3-2.001
Date: 2024-09-06
Version: 3.3-2.002
Date: 2024-10-06
Title: Exploratory Data Analysis for the 'spatstat' Family
Authors@R: c(person("Adrian", "Baddeley",
role = c("aut", "cre", "cph"),
Expand Down
12 changes: 12 additions & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,8 @@ export("bind.ratfv")
export("bits.envelope")
export("bits.test")
export("blur")
export("blurHeat")
export("blurHeat.im")
export("boyce")
export("bw.abram.ppp")
export("bw.CvL")
Expand All @@ -78,6 +80,7 @@ export("bw.pcf")
export("bw.ppl")
export("bw.pplHeat")
export("bw.relrisk")
export("bw.relriskHeatppp")
export("bw.relrisk.ppp")
export("bw.scott")
export("bw.scott.iso")
Expand Down Expand Up @@ -424,6 +427,8 @@ export("rectcontact")
export("RelevantDeviation")
export("reload.or.compute")
export("relrisk")
export("relriskHeat")
export("relriskHeat.ppp")
export("relrisk.ppp")
export("rename.fv")
export("resolve.2D.kernel")
Expand Down Expand Up @@ -482,6 +487,9 @@ export("smoothcrossEngine")
export("Smoothfun")
export("Smoothfun.ppp")
export("Smooth.fv")
export("SmoothHeat")
export("SmoothHeat.im")
export("SmoothHeat.ppp")
export("Smooth.im")
export("smoothpointsEngine")
export("Smooth.ppp")
Expand Down Expand Up @@ -568,6 +576,7 @@ S3method("as.data.frame", "fv")
S3method("as.tess", "quadrattest")
S3method("auc", "ppp")
S3method("berman.test", "ppp")
S3method("blurHeat", "im")
S3method("bw.abram", "ppp")
S3method("bw.relrisk", "ppp")
S3method("cbind", "fv")
Expand Down Expand Up @@ -671,6 +680,7 @@ S3method("quadrat.test", "quadratcount")
S3method("quadrat.test", "splitppp")
S3method("range", "ssf")
S3method("[", "rat")
S3method("relriskHeat", "ppp")
S3method("relrisk", "ppp")
S3method("resolve.lambdacross", "ppp")
S3method("resolve.lambda", "ppp")
Expand All @@ -691,6 +701,8 @@ S3method("shift", "quadrattest")
S3method("simulate", "rhohat")
S3method("Smoothfun", "ppp")
S3method("Smooth", "fv")
S3method("SmoothHeat", "im")
S3method("SmoothHeat", "ppp")
S3method("Smooth", "im")
S3method("Smooth", "ppp")
S3method("Smooth", "solist")
Expand Down
20 changes: 19 additions & 1 deletion NEWS
Original file line number Diff line number Diff line change
@@ -1,10 +1,28 @@

CHANGES IN spatstat.explore VERSION 3.3-2.001
CHANGES IN spatstat.explore VERSION 3.3-2.002

OVERVIEW

o relative risk estimation using diffusion.

o smoothing using diffusion.

o Tweaks to bandwidth selection.

NEW FUNCTIONS

o relriskHeat, relriskHeat.ppp
Relative risk estimation using diffusion.

o blurHeat, blurHeat.im
Image smoothing using diffusion.

o SmoothHeat, SmoothHeat.ppp
Smoothing numerical values observed at points, using diffusion.

o bw.relriskHeatppp
Bandwidth selection for relriskHeat.ppp

SIGNIFICANT USER-VISIBLE CHANGES

o bw.ppl
Expand Down
78 changes: 78 additions & 0 deletions R/SmoothHeat.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,78 @@
#'
#' SmoothHeat.R
#'
#' Nadaraya-Watson style smooth regression using diffusion
#'
#' Copyright (C) 2018-2024 Adrian Baddeley, Tilman Davies and Suman Rakshit
#'
#' $Revision: 1.3 $ $Date: 2024/10/06 01:26:29 $

SmoothHeat <- function(X, ...) {
UseMethod("SmoothHeat")
}

SmoothHeat.im <- function(X, sigma, ...) {
blurHeat(X, sigma, ...)
}

SmoothHeat.ppp <- function(X, sigma, ..., weights=NULL) {
stopifnot(is.ppp(X))
stopifnot(is.marked(X))
marx <- marks(X)
if(!is.vector(marx)) stop("Marks of X should be a numeric vector")
marx <- as.numeric(marx)
if(is.null(weights)) {
numwt <- marx
denwt <- NULL
} else {
check.nvector(weights, npoints(X), oneok=TRUE)
if(length(weights) == 1) weights <- rep(weights, npoints(X))
numwt <- marx * weights
denwt <- weights
}
Y <- unmark(X)
numer <- densityHeat(Y, sigma, weights=numwt, ...)
denom <- densityHeat(Y, sigma, weights=denwt, ...)
return(numer/denom)
}

bw.SmoothHeatppp <- function(X, ..., weights=NULL,
srange=NULL, ns=16, sigma=NULL,
leaveoneout=TRUE, verbose=TRUE) {
stopifnot(is.ppp(X))
stopifnot(is.marked(X))
marx <- marks(X)
if(!is.vector(marx)) stop("Marks of X should be a numeric vector")
marx <- as.numeric(marx)
if(is.null(weights)) {
numwt <- marx
denwt <- NULL
} else {
check.nvector(weights, npoints(X), oneok=TRUE)
if(length(weights) == 1) weights <- rep(weights, npoints(X))
numwt <- marx * weights
denwt <- weights
}
#' compute weighted and unweighted intensity estimates
U <- unmark(X)
aNumer <- HeatEstimates.ppp(U, ..., weights=numwt,
srange=srange, ns=ns, sigma=sigma,
leaveoneout=leaveoneout, verbose=verbose)
aDenom <- HeatEstimates.ppp(U, ..., weights=denwt,
srange=srange, ns=ns, sigma=sigma,
leaveoneout=leaveoneout, verbose=verbose)
h <- aDenom$h
hname <- aDenom$hname
#' compute smoother
zhat <- aNumer$lambda/aDenom$lambda
#' compute least squares cross-validation criterion
zobs <- matrix(marx, nrow(zhat), ncol(zhat), byrow=TRUE)
CV <- rowSums((zhat - zobs)^2)
iopt <- which.min(CV)
result <- bw.optim(CV, h, iopt,
criterion="Least squares cross-validation",
hname=hname,
unitname=unitname(X))
return(result)
}

137 changes: 137 additions & 0 deletions R/blurHeat.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,137 @@
#'
#' blurHeat.R
#'
#' Image blurring by diffusion
#'
#' Copyright (C) 2018-2024 Adrian Baddeley, Tilman Davies and Suman Rakshit
#'
#' Licence: GNU Public Licence >= 2
#'
#' $Revision: 1.3 $ $Date: 2024/10/06 02:28:55 $

blurHeat <- function(X, ...) {
UseMethod("blurHeat")
}

blurHeat.im <- function(X, sigma, ..., connect=8,
symmetric=FALSE, k=1, show=FALSE) {
Y <- as.im(X)
check.1.integer(k)
stopifnot(k >= 1)
if(!(connect %in% c(4,8)))
stop("connectivity must be 4 or 8")
if(is.im(sigma)) {
# ensure Y and sigma are on the same grid
A <- harmonise(Y=Y, sigma=sigma)
Y <- A$Y
sigma <- A$sigma
} else if(is.function(sigma)) {
sigma <- as.im(sigma, as.owin(Y))
} else check.1.real(sigma)
#' initial state
v <- as.matrix(Y)
u <- as.vector(v)
#' symmetric random walk?
if(symmetric) {
asprat <- with(Y, ystep/xstep)
if(abs(asprat-1) > 0.01)
warning(paste("Symmetric random walk on a non-square grid",
paren(paste("aspect ratio", asprat))),
call.=FALSE)
}
#' determine appropriate jump probabilities & time step
pmax <- 1/(connect+1) # maximum permitted jump probability
xstep <- Y$xstep
ystep <- Y$ystep
minstep <- min(xstep, ystep)
if(symmetric) {
#' all permissible transitions have the same probability 'pjump'.
#' Determine Nstep, and dt=sigma^2/Nstep, such that
#' Nstep >= 16 and M * pjump * minstep^2 = dt
M <- if(connect == 4) 2 else 6
Nstep <- max(16, ceiling(max(sigma)^2/(M * pmax * minstep^2)))
sn <- (sigma^2)/Nstep
px <- py <- pxy <- sn/(M * minstep^2)
} else {
#' px is the probability of jumping 1 step to the right
#' py is the probability of jumping 1 step up
#' if connect=4, horizontal and vertical jumps are exclusive.
#' if connect=8, horizontal and vertical increments are independent
#' Determine Nstep, and dt = sigma^2/Nstep, such that
#' Nstep >= 16 and 2 * pmax * minstep^2 = dt
Nstep <- max(16, ceiling(max(sigma)^2/(2 * pmax * minstep^2)))
sn <- (sigma^2)/Nstep
px <- sn/(2 * xstep^2)
py <- sn/(2 * ystep^2)
if(max(px) > pmax) stop("Internal error: px exceeds pmax")
if(max(py) > pmax) stop("Internal error: py exceeds pmax")
if(connect == 8) pxy <- px * py
}
#' construct adjacency matrices
dimv <- dim(v)
my <- gridadjacencymatrix(dimv, across=FALSE, down=TRUE, diagonal=FALSE)
mx <- gridadjacencymatrix(dimv, across=TRUE, down=FALSE, diagonal=FALSE)
if(connect == 8)
mxy <- gridadjacencymatrix(dimv, across=FALSE, down=FALSE, diagonal=TRUE)
#' restrict to window
if(anyNA(u)) {
ok <- !is.na(u)
u <- u[ok]
mx <- mx[ok,ok,drop=FALSE]
my <- my[ok,ok,drop=FALSE]
if(connect == 8)
mxy <- mxy[ok,ok,drop=FALSE]
if(is.im(sigma)) {
px <- px[ok]
py <- py[ok]
if(connect == 8)
pxy <- pxy[ok]
}
} else ok <- TRUE
#' construct iteration matrix
if(connect == 4) {
A <- px * mx + py * my
} else {
A <- px * (1 - 2 * py) * mx + py * (1 - 2 * px) * my + pxy * mxy
}
#' debug
stopifnot(min(rowSums(A)) >= 0)
stopifnot(max(rowSums(A)) <= 1)
#'
diag(A) <- 1 - rowSums(A)
#' k-step transition probabilities
if(k > 1) {
Ak <- A
for(j in 2:k) Ak <- Ak %*% A
} else Ak <- A
k <- as.integer(k)
Nstep <- as.integer(Nstep)
Nblock <- Nstep/k
Nrump <- Nstep - Nblock * k
#' run
U <- u
Z <- Y
if(!show) {
for(istep in 1:Nblock) U <- U %*% Ak
} else {
opa <- par(ask=FALSE)
each <- max(1, round(Nblock/60))
for(istep in 1:Nblock) {
U <- U %*% Ak
if(istep %% each == 0) {
Z[] <- as.vector(U)
f <- sqrt(istep/Nstep)
main <- if(is.im(sigma)) paste(signif(f, 3), "* sigma") else
paste("sigma =", signif(f * sigma, 3))
plot(Z, main=main)
Sys.sleep(0.4)
}
}
par(opa)
}
if(Nrump > 0) for(istep in 1:Nrump) U <- U %*% A
#' pack up
Z[] <- as.vector(U)
return(Z)
}

4 changes: 4 additions & 0 deletions R/densityHeat.ppp.R
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,10 @@
#'
#' Diffusion estimator of density/intensity
#'
#' Copyright (C) 2018-2024 Adrian Baddeley, Tilman Davies and Suman Rakshit
#'
#' Licence: GNU Public Licence >= 2
#'

densityHeat <- function(x, sigma, ...) {
UseMethod("densityHeat")
Expand Down
Loading

0 comments on commit fe2ba51

Please sign in to comment.