R/mesh.R
In inbo/INLA: Full Bayesian Analysis of Latent Gaussian Models using Integrated Nested Laplace Approximations
Defines functions
inla.mesh.segment
inla.mesh.segment.default
inla.mesh.segment.inla.mesh.segment
lines.inla.mesh.segment
`inla.generate.colors`
`plot.inla.trimesh`
plot.inla.mesh
inla.mesh.map.lim
inla.mesh.map
inla.mesh.lattice
extract.groups
extract.groups.inla.mesh.segment
inla.mesh.parse.segm.input
inla.mesh.filter.locations
inla.mesh
inla.mesh.create
inla.mesh.extract.segments
inla.mesh.boundary
inla.mesh.interior
inla.mesh.2d
inla.mesh.create.helper
inla.delaunay
inla.mesh.query
summary.inla.mesh
print.summary.inla.mesh
inla.mesh.project
inla.mesh.project.inla.mesh
inla.mesh.project.inla.mesh.1d
inla.mesh.project.inla.mesh.projector
inla.mesh.projector
inla.mesh.projector.inla.mesh
inla.mesh.projector.inla.mesh.1d
inla.internal.make.spline.mesh
inla.mesh.basis
inla.parse.queries
`inla.fmesher.smorg`
inla.mesh.1d
inla.mesh.1d.bary
inla.mesh.1d.A
inla.mesh.1d.fem
inla.diameter
inla.diameter.default
inla.diameter.inla.mesh.1d
inla.diameter.inla.mesh
inla.diameter.inla.mesh.segment
inla.diameter.inla.mesh.lattice
inla.mesh.fem
inla.mesh.deriv
inla.simplify.curve
inla.contour.segment
inla.nonconvex.hull.basic
inla.nonconvex.hull
Documented in
extract.groups
extract.groups.inla.mesh.segment
inla.contour.segment
inla.delaunay
inla.diameter
inla.diameter.default
inla.diameter.inla.mesh
inla.diameter.inla.mesh.1d
inla.diameter.inla.mesh.lattice
inla.diameter.inla.mesh.segment
inla.mesh
inla.mesh.1d
inla.mesh.1d.A
inla.mesh.1d.bary
inla.mesh.1d.fem
inla.mesh.2d
inla.mesh.basis
inla.mesh.boundary
inla.mesh.create
inla.mesh.create.helper
inla.mesh.deriv
inla.mesh.fem
inla.mesh.interior
inla.mesh.lattice
inla.mesh.map
inla.mesh.map.lim
inla.mesh.project
inla.mesh.project.inla.mesh
inla.mesh.project.inla.mesh.1d
inla.mesh.project.inla.mesh.projector
inla.mesh.projector
inla.mesh.projector.inla.mesh
inla.mesh.projector.inla.mesh.1d
inla.mesh.query
inla.mesh.segment
inla.mesh.segment.default
inla.mesh.segment.inla.mesh.segment
inla.nonconvex.hull
inla.nonconvex.hull.basic
inla.simplify.curve
lines.inla.mesh.segment
plot.inla.mesh
print.summary.inla.mesh
summary.inla.mesh
## Export: inla.contour.segment inla.delaunay inla.fmesher.smorg
## Export: inla.generate.colors inla.mesh inla.mesh.1d inla.mesh.1d.A
## Export: inla.mesh.1d.bary inla.mesh.1d.fem inla.mesh.2d inla.mesh.basis
## Export: inla.mesh.boundary inla.mesh.create inla.mesh.create.helper
## Export: inla.mesh.deriv inla.mesh.fem
## Export: inla.mesh.interior inla.mesh.lattice inla.mesh.map
## Export: inla.mesh.map.lim inla.mesh.query
## Export: inla.nonconvex.hull inla.nonconvex.hull.basic
## Export: inla.simplify.curve plot.inla.trimesh plot!inla.trimesh
## Internal: inla.mesh.filter.locations
## Internal: inla.mesh.parse.segm.input inla.mesh.extract.segments
##
## S3methods; also export some methods explicitly
## Export: extract.groups inla.mesh.project inla.mesh.projector
## Export: extract.groups!inla.mesh.segment
## Export: inla.mesh.segment
## Export: inla.mesh.segment.default
## Export: inla.mesh.segment.inla.mesh.segment
## Export: inla.mesh.segment!default
## Export: inla.mesh.segment!inla.mesh.segment
## Export: inla.mesh.project.inla.mesh inla.mesh.project.inla.mesh.1d
## Export: inla.mesh.project.inla.mesh.projector
## Export: inla.mesh.projector.inla.mesh inla.mesh.projector.inla.mesh.1d
## Export: inla.mesh.project!inla.mesh inla.mesh.project!inla.mesh.1d
## Export: inla.mesh.project!inla.mesh.projector
## Export: inla.mesh.projector!inla.mesh inla.mesh.projector!inla.mesh.1d
## Export: lines.inla.mesh.segment plot.inla.mesh
## Export: print.summary.inla.mesh summary.inla.mesh
## Export: lines!inla.mesh.segment plot!inla.mesh
## Export: print!summary.inla.mesh summary!inla.mesh
## Export: inla.diameter
## Export: inla.diameter!default
## Export: inla.diameter!inla.mesh.1d
## Export: inla.diameter!inla.mesh
## Export: inla.diameter!inla.mesh.segment
## Export: inla.diameter!inla.mesh.lattice
inla.mesh.segment <- function(...) {
UseMethod("inla.mesh.segment")
}
inla.mesh.segment.default <-
function(loc = NULL, idx = NULL, grp = NULL, is.bnd = TRUE,
crs=NULL, ...)
{
if ((missing(loc) || is.null(loc)) &&
(missing(idx) || is.null(idx))) {
stop("At most one of 'loc' and 'idx' may be missing or null.")
}
if (!missing(loc) && !is.null(loc)) {
## Handle loc given as SpatialPoints or SpatialPointsDataFrame object
if (inherits(loc, "SpatialPoints") ||
inherits(loc, "SpatialPointsDataFrame")) {
loc <- inla.spTransform(coordinates(loc),
CRS(proj4string(loc)),
crs,
passthrough=TRUE)
}
if (!is.matrix(loc)) {
loc = as.matrix(loc)
}
if (!is.double(loc)) {
storage.mode(loc) = "double"
}
if (missing(idx) || is.null(idx))
idx = (inla.ifelse(is.bnd,
c(1:nrow(loc),1),
c(1:nrow(loc))))
} else {
loc = NULL
}
if (!missing(idx) && !is.null(idx)) {
if (!is.vector(idx) && !is.matrix(idx))
stop("'idx' must be a vector or a matrix")
if (is.vector(idx))
idx = as.matrix(idx, nrow=length(idx), ncol=1)
if (ncol(idx) == 1) {
if (nrow(idx) < 2) {
if (nrow(idx) == 1) {
warning("Segment specification must have at least 2, or 0, indices.")
}
idx <- matrix(0L, 0, 2)
} else {
idx = matrix(c(idx[-nrow(idx)],idx[-1]), nrow(idx)-1, 2)
}
}
storage.mode(idx) <- "integer"
if (!is.null(loc) &&
(nrow(idx) > 0) &&
(max(idx, na.rm=TRUE) > nrow(loc))) {
warning("Segment indices (max=", max(idx, na.rm=TRUE),
") exceed specified location list length (",
nrow(loc), ").")
}
}
if (!missing(grp) && !is.null(grp)) {
if (!is.vector(grp) && !is.matrix(grp))
stop("'grp' must be a vector or a matrix")
grp = matrix(grp, min(length(grp), nrow(idx)), 1)
if (nrow(grp)<nrow(idx))
grp = (matrix(c(as.vector(grp),
rep(grp[nrow(grp)], nrow(idx)-length(grp))),
nrow(idx), 1))
storage.mode(grp) <- "integer"
} else
grp = NULL
## Filter away NAs in loc and idx
if (!is.null(loc)) {
idx[is.na(idx)] = 0L ## Avoid R annoyances with logical+NA indexing
while (sum(is.na(loc))>0) {
i = min(which(rowSums(is.na(loc))>0))
loc = loc[-i,,drop=FALSE]
idx[idx==i] = 0L
idx[idx>i] = idx[idx>i]-1L
}
idx[idx==0L] = NA
}
while (sum(is.na(idx))>0) {
i = min(which(rowSums(is.na(idx))>0))
idx = idx[-i,,drop=FALSE]
if (!is.null(grp))
grp = grp[-i,,drop=FALSE]
}
if (!is.null(loc)) {
## Identify unused locations and remap indices accordingly.
idx.new = rep(0L, nrow(loc))
idx.new[as.vector(idx)] = 1L
loc = loc[idx.new==1L,, drop=FALSE]
idx.new[idx.new==1L] = seq_len(sum(idx.new))
idx = (matrix(idx.new[as.vector(idx)],
nrow=nrow(idx),
ncol=ncol(idx)))
}
ret = list(loc=loc, idx=idx, grp=grp, is.bnd=is.bnd, crs=crs)
class(ret) <- "inla.mesh.segment"
return(invisible(ret))
}
inla.mesh.segment.inla.mesh.segment <- function(..., grp.default=0) {
segm <- list(...)
if (!all(unlist(lapply(segm,
function(x) inherits(x,"inla.mesh.segment"))))) {
stop("All objects must be of class 'inla.mesh.segment'.")
}
Nloc <- unlist(lapply(segm, function(x) nrow(x$loc)))
cumNloc <- c(0, cumsum(Nloc))
Nidx <- unlist(lapply(segm, function(x) nrow(x$idx)))
loc <- do.call(rbind, lapply(segm, function(x) x$loc))
idx <- do.call(rbind, lapply(seq_along(segm),
function(x) segm[[x]]$idx+cumNloc[x]))
grp <- unlist(lapply(seq_along(segm),
function(x) {
if (is.null(segm[[x]]$grp)) {
rep(grp.default, Nidx[x])
} else {
segm[[x]]$grp
}
}))
is.bnd <- unlist(lapply(segm, function(x) x$is.bnd))
if (!all(is.bnd) || (any(!is.bnd) && !all(!is.bnd))) {
warning("Inconsistent 'is.bnd' attributes. Setting 'is.bnd=FALSE'.")
is.bnd <- FALSE
} else {
is.bnd <- all(is.bnd)
}
crs <- lapply(segm, function(x) x$crs)
if (!is.null(crs)) {
crs <- crs[unlist(lapply(crs,
function(x) !is.null(x)))]
if (length(crs) > 0) {
if (!all(unlist(lapply(crs,
function(x) identical(crs[[1]], x))))) {
lapply(crs, function(x) show(x))
stop("Inconsistent 'crs' attributes.")
} else {
crs <- crs[[1]]
}
} else {
crs <- NULL
}
}
invisible(inla.mesh.segment(loc=loc, idx=idx, grp=grp, is.bnd=is.bnd,
crs=crs))
}
lines.inla.mesh.segment <- function(x, loc=NULL, col=NULL,
colors=c("black", "blue", "red", "green"),
add=TRUE, xlim=NULL, ylim=NULL,
rgl=FALSE, ...)
{
segm = x
if (!is.null(segm$loc))
loc = segm$loc
stopifnot(!is.null(loc), ncol(loc)>=2)
if (ncol(loc) < 3) {
loc = cbind(loc, 0.0)
}
color = col
if (rgl) {
stopifnot(inla.require("rgl"))
if (!add) {
dev = rgl::open3d()
rgl::view3d(0, 0, fov=0)
} else {
dev = NULL
}
} else if (!add) {
idx = unique(as.vector(segm$idx))
if (is.null(xlim))
xlim=range(loc[idx,1])
if (is.null(ylim))
ylim=range(loc[idx,2])
plot.new()
plot.window(xlim=xlim, ylim=ylim, ...)
}
grps = inla.ifelse(is.null(segm$grp), rep(0L,nrow(segm$idx)), segm$grp)
for (grp in unique(grps)) {
idx = which(grps==grp)
if (is.null(col)) {
color=colors[1+(grp%%length(colors))]
}
if (rgl) {
rgl::segments3d(loc[as.vector(t(segm$idx[idx,, drop=FALSE])),,drop=FALSE],
color=color,
...)
} else {
lines(loc[t(cbind(segm$idx[idx,, drop=FALSE], NA)), 1],
loc[t(cbind(segm$idx[idx,, drop=FALSE], NA)), 2],
col=color,
...)
}
}
}
`inla.generate.colors` <- function(color,
color.axis = NULL,
color.n=512,
color.palette = cm.colors,
color.truncate=FALSE,
alpha=NULL)
{
if (is.character(color)) {
colors = color
} else if (is.vector(color) || (is.matrix(color) && (ncol(color)==1))) {
if (is.null(color.axis))
color.axis = c(min(color, na.rm=TRUE), max(color, na.rm=TRUE))
if (color.truncate) {
not.ok = ((color<color.axis[1]) |
(color>color.axis[2]))
} else {
not.ok = rep(FALSE, length(color))
}
cs = (pmax(color.axis[1],
pmin(color.axis[2], color, na.rm=TRUE), na.rm=TRUE))
cs = (cs-color.axis[1])/(color.axis[2]-color.axis[1])
not.ok = not.ok | is.na(cs)
cs[not.ok] = 0.5
if (is.null(alpha)) {
alpha = as.numeric(!not.ok)
} else {
alpha[not.ok] = 0
}
ics = (as.numeric(cut(cs, seq(0, 1, length.out=color.n+1),
include.lowest=TRUE)))
colors = color.palette(color.n)[ics]
## Todo: handle alpha, combining "input alpha" with "not.ok-alpha"
} else if (is.matrix(color) && (ncol(color)==3)) {
if (is.null(color.axis))
color.axis = c(min(color, na.rm=TRUE), max(color, na.rm=TRUE))
if (color.truncate) {
not.ok = ((color[, 1]<color.axis[1]) |
(color[, 2]<color.axis[1]) |
(color[, 3]<color.axis[1]) |
(color[, 1]>color.axis[2]) |
(color[, 2]>color.axis[2]) |
(color[, 3]>color.axis[2]))
} else {
not.ok = rep(FALSE, nrow(color))
}
cs = matrix(
pmax(color.axis[1],
pmin(color.axis[2], color, na.rm=TRUE), na.rm=TRUE), dim(color))
cs = (cs-color.axis[1])/(color.axis[2]-color.axis[1])
not.ok = not.ok | is.na(cs[, 1]) | is.na(cs[, 2]) | is.na(cs[, 3])
cs[not.ok,] = c(0.5, 0.5, 0.5)
if (is.null(alpha)) {
alpha = as.numeric(!not.ok)
} else {
alpha[not.ok] = 0
}
colors = rgb(cs[, 1], cs[, 2], cs[, 3])
} else {
stop("color specification must be character, matrix, or vector.")
}
return (list(colors=colors, alpha=alpha))
}
## Export as plot!inla.trimesh as well as plot.inla.trimesh,
## even though it's not actually an S3 method.
`plot.inla.trimesh` <- function(x, S, color = NULL, color.axis = NULL,
color.n=512, color.palette = cm.colors,
color.truncate=FALSE, alpha=NULL,
lwd = 1, specular = "black",
draw.vertices=TRUE, draw.edges=TRUE,
edge.color=rgb(0.3, 0.3, 0.3), ...)
{
TV = x
stopifnot(inla.require("rgl"))
## Make indices 1 based. Deprecated and is deactivated.
if (min(TV) == 0) {
stop("Zero-based indices in TV are not supported.")
}
colors = (inla.generate.colors(color, color.axis, color.n,
color.palette, color.truncate, alpha))
tTV = t(TV);
tETV = t(TV[, c(1, 2, 3, 1, NA)]);
Tx = S[tTV, 1]
Ty = S[tTV, 2]
Tz = S[tTV, 3]
if (length(colors$colors) == 1) {
## One color
Tcol = colors$colors
Talpha = colors$alpha
} else if (length(colors$colors) == nrow(S)) {
## One color per vertex
Tcol = colors$colors[tTV]
Talpha = colors$alpha[tTV]
} else {
## One color per triangle
stopifnot(length(colors$colors) == nrow(TV))
Tcol = colors$colors[t(matrix(rep(1:nrow(TV), 3), dim(TV)))]
Talpha = colors$alpha[t(matrix(rep(1:nrow(TV), 3), dim(TV)))]
}
Ex = S[tETV, 1]
Ey = S[tETV, 2]
Ez = S[tETV, 3]
Ecol = edge.color
if (draw.vertices) {
rgl::points3d(S, color="black", ...)
}
if (draw.edges) {
rgl::lines3d(Ex, Ey, Ez, color=Ecol, lwd=lwd, ...)
}
rgl::triangles3d(Tx, Ty, Tz, color=Tcol, specular=specular, alpha=Talpha, ...)
return (invisible())
}
## library(geometry)
## S = cbind(x=rnorm(30), y=rnorm(30), z=0)
## TV = delaunayn(S[, 1:2]) # NOTE: inconsistent triangle orders, only for test.
## trimesh(TV, S)
##
## colors = rgb(runif(30), runif(30), runif(30))
## rgl.viewpoint(0, 0, fov=20)
## plot.inla.trimesh(TV, S, colors)
## Ecol = col2rgb(color)/256
## Ecol = Ecol*0.5+(1-0.5)*0 # Rescale towards black
## Ecol = 1-Ecol # Invert
## Ecol = Ecol[, c(2, 3, 1)] # Permute
## Ecol = rgb(Ecol[1,], Ecol[2,], Ecol[3,], maxColorValue = 1)
## Ecol = Ecol[tETV]
plot.inla.mesh <- function(x,
col="white",
t.sub=1:nrow(mesh$graph$tv),
add=FALSE,
lwd=1,
xlim = range(mesh$loc[,1]),
ylim = range(mesh$loc[,2]),
main = NULL,
rgl = FALSE,
size = 2,
draw.vertices=FALSE,
vertex.color="black",
draw.edges=TRUE,
edge.color=rgb(0.3, 0.3, 0.3),
draw.segments=draw.edges,
...)
{
inla.require.inherits(x, "inla.mesh", "'mesh'")
mesh = x
if (rgl) {
stopifnot(inla.require("rgl"))
if (!add) {
dev = rgl::open3d()
rgl::view3d(0, 0, fov=0)
} else {
dev = NULL
}
tv = mesh$graph$tv[t.sub,,drop=FALSE]
if (draw.vertices) {
idx = intersect(unique(as.vector(tv)), mesh$idx$loc)
rgl::points3d(mesh$loc[idx,,drop=FALSE],
size=2*size, lwd=lwd, color = "blue", ...)
}
if (draw.segments) {
if (!is.null(mesh$segm$bnd))
lines(mesh$segm$bnd, mesh$loc, lwd=lwd+1,
rgl=TRUE, add=TRUE, ...)
if (!is.null(mesh$segm$int))
lines(mesh$segm$int, mesh$loc, lwd=lwd+1,
rgl=TRUE, add=TRUE, ...)
}
plot.inla.trimesh(tv, mesh$loc, color = col,
size=size, lwd=lwd,
draw.vertices=draw.vertices,
vertex.color=vertex.color,
draw.edges=draw.edges,
edge.color=edge.color,
...)
return(invisible(dev))
} else {
idx = cbind(mesh$graph$tv[t.sub,c(1:3,1), drop=FALSE], NA)
x = mesh$loc[t(idx), 1]
y = mesh$loc[t(idx), 2]
if (!add) {
plot.new()
plot.window(xlim=xlim, ylim=ylim, ...)
}
if (draw.edges) {
lines(x, y, type="l", col=edge.color, lwd=lwd)
}
tv = mesh$graph$tv[t.sub,,drop=FALSE]
if (draw.vertices) {
idx = unique(as.vector(tv))
points(mesh$loc[idx,,drop=FALSE],
pch=20, col=vertex.color, ...)
idx = intersect(idx, mesh$idx$loc)
points(mesh$loc[idx,,drop=FALSE],
pch=20, col="blue", ...)
}
if (draw.segments) {
if (!is.null(mesh$segm$bnd))
lines(mesh$segm$bnd, mesh$loc, lwd=lwd+1, ...)
if (!is.null(mesh$segm$int))
lines(mesh$segm$int, mesh$loc, lwd=lwd+1, ...)
}
if (!add && missing(main)) {
if (mesh$meta$is.refined) {
title("Constrained refined Delaunay triangulation", ...)
} else {
title("Constrained Delaunay triangulation", ...)
}
}
return(invisible())
}
}
inla.mesh.map.lim <-
function(loc=NULL,
projection=
c("default", "longlat", "longsinlat", "mollweide"))
{
projection = match.arg(projection)
if (identical(projection, "default")) {
if (is.null(loc)) {
lim = list(xlim=c(0, 1), ylim=c(0, 1))
} else {
lim =
list(xlim=range(loc[,1], na.rm=TRUE),
ylim=range(loc[,2], na.rm=TRUE))
}
} else if (identical(projection, "longlat")) {
lim = list(xlim=c(-180,180), ylim=c(-90,90))
} else if (identical(projection, "longsinlat")) {
lim = list(xlim=c(-180,180), ylim=c(-1,1))
} else if (identical(projection, "mollweide")) {
lim = list(xlim=c(-2,2), ylim=c(-1,1))
} else {
stop(paste("Unknown projection '", projection, "'.", sep=""))
}
return(lim)
}
inla.mesh.map <-
function(loc,
projection=
c("default", "longlat", "longsinlat", "mollweide"),
inverse=TRUE)
{
projection = match.arg(projection)
if (identical(projection, "default")) {
return(loc)
} else if (identical(projection, "longlat")) {
if (inverse) {
proj =
cbind(cos(loc[,1]*pi/180)*cos(loc[,2]*pi/180),
sin(loc[,1]*pi/180)*cos(loc[,2]*pi/180),
sin(loc[,2]*pi/180))
} else {
proj =
cbind(atan2(loc[,2], loc[,1])*180/pi,
asin(pmax(-1, pmin(+1, loc[,3])))*180/pi)
}
} else if (identical(projection, "longsinlat")) {
if (inverse) {
coslat = sqrt(pmax(0, 1-loc[,2]^2))
proj =
cbind(cos(loc[,1]*pi/180)*coslat,
sin(loc[,1]*pi/180)*coslat,
loc[,2]
)
} else {
proj =
cbind(atan2(loc[,2], loc[,1])*180/pi,
loc[,3])
}
} else if (identical(projection, "mollweide")) {
if (inverse) {
ok = ((loc[,1]^2+4*loc[,2]^2) <= 4)
cos.theta = sqrt(pmax(0, 1-loc[ok,2]^2))
theta = atan2(loc[ok,2], cos.theta)
sin.lat = (2*theta+sin(2*theta))/pi
cos.lat = sqrt(pmax(0, 1-sin.lat^2))
lon = loc[ok,1]*pi/2/(cos.theta+(cos.theta==0))
lon[cos.theta==0] = pi/2*sign(theta[cos.theta==0])
proj = matrix(NA, nrow(loc), 3)
proj[ok,] = cbind(cos(lon)*cos.lat, sin(lon)*cos.lat, sin.lat)
} else {
lon = atan2(loc[,2],loc[,1])
z = pmin(1, pmax(-1, loc[,3]))
sin.theta = z
cos.theta = sqrt(pmax(0, 1-sin.theta^2))
## NR-solver for sin.theta.
## Typically finishes after at most 7 iterations.
## When cos.theta=0, sin.theta is already correct, +/- 1.
nook = (cos.theta > 0)
for (k in 1:20) {
if (any(nook)) {
delta =
(atan2(sin.theta[nook], cos.theta[nook]) +
sin.theta[nook]*cos.theta[nook] - pi/2*z[nook])/
(2*cos.theta[nook])
sin.theta[nook] = sin.theta[nook] - delta
cos.theta[nook] = sqrt(1-sin.theta[nook]^2)
nook[nook] = (abs(delta) > 1e-14)
}
}
proj = cbind(2*lon/pi*cos.theta, sin.theta)
}
} else {
stop(paste("Unknown projection '", projection, "'.", sep=""))
}
return(proj)
}
inla.mesh.lattice <- function(x=seq(0, 1, length.out=2),
y=seq(0, 1, length.out=2),
z=NULL,
dims =
if (is.matrix(x)) {
dim(x)
} else {
c(length(x), length(y))
},
units = NULL,
crs = NULL)
{
if (is.null(crs)) {
units = match.arg(units, c("default", "longlat", "longsinlat", "mollweide"))
lim = inla.mesh.map.lim(projection=units)
xlim = lim$xlim
ylim = lim$ylim
} else { ## !is.null(crs)
if (!is.null(units)) {
stop("Only one of 'units' and 'crs' can be non-null.")
}
bounds <- inla.crs.bounds(crs)
xlim = bounds$xlim
ylim = bounds$ylim
}
if (missing(x) && !missing(dims)) {
x = seq(xlim[1], xlim[2], length.out=dims[1])
}
if (missing(y) && !missing(dims)) {
y = seq(ylim[1], ylim[2], length.out=dims[2])
}
dims = as.integer(dims)
if (is.matrix(x)) {
if (!identical(dims, dim(x)) ||
!identical(dims, dim(y)) ||
(is.matrix(z) && !identical(dims, dim(z))))
stop("The size of matrices 'x', 'y', and 'z' must match 'dims'.")
loc = cbind(as.vector(x), as.vector(y), as.vector(z))
x = NULL
y = NULL
} else {
if (!identical(dims[1], length(x)) ||
!identical(dims[2], length(y)))
stop(paste("The lengths of vectors 'x' and 'y' (",
length(x),",",length(y),
") must match 'dims' (",dims[1],",",dims[2],").",
sep=""))
loc = (cbind(rep(x, times = dims[2]),
rep(y, each = dims[1])))
}
if (!is.double(loc))
storage.mode(loc) = "double"
if (is.null(crs)) {
loc = inla.mesh.map(loc=loc, projection=units, inverse=TRUE)
}
## Construct lattice boundary
segm.idx = (c(1:(dims[1]-1),
dims[1]*(1:(dims[2]-1)),
dims[1]*dims[2]-(0:(dims[1]-2)),
dims[1]*((dims[2]-1):1)+1))
segm.grp = (c(rep(1L, dims[1]-1),
rep(2L, dims[2]-1),
rep(3L, dims[1]-1),
rep(4L, dims[2]-1)))
segm = (inla.mesh.segment(loc=loc[segm.idx,, drop=FALSE],
grp=segm.grp,
is.bnd=TRUE,
crs=crs))
lattice = list(dims=dims, x=x, y=y, loc=loc, segm=segm, crs=crs)
class(lattice) = "inla.mesh.lattice"
return(lattice)
}
extract.groups <- function(...)
{
UseMethod("extract.groups")
}
extract.groups.inla.mesh.segment <- function(segm,
groups,
groups.new=groups,
...)
{
inla.require.inherits(segm, "inla.mesh.segment", "'segm'")
if (length(groups.new)==1L) {
groups.new = rep(groups.new, length(groups))
}
if (length(groups.new)!=length(groups)) {
stop("Length of 'groups.new' (", length(groups.new),
") does not match length of 'groups' (",length(groups),")")
}
idx = c()
segm.grp = c()
for (k in 1:length(groups)) {
extract.idx = which(segm$grp==groups[k])
idx = c(idx, extract.idx)
segm.grp = c(segm.grp, rep(groups.new[k], length(extract.idx)))
}
segm.idx = segm$idx[idx,, drop=FALSE]
return(inla.mesh.segment(loc=segm$loc,
idx=segm.idx,
grp=segm.grp,
segm$is.bnd,
crs=segm$crs))
}
inla.mesh.parse.segm.input <- function(boundary=NULL,
interior=NULL,
segm.offset=0L,
loc.offset=0L,
crs=NULL)
{
###########################################
homogenise.segm.input <- function(x, is.bnd, crs=NULL)
{
if (is.matrix(x) || is.vector(x)) { ## Coordinates or indices
x = (inla.ifelse(is.matrix(x),x,
as.matrix(x,nrow=length(x),ncol=1)))
if (is.integer(x)) { ## Indices
ret = inla.mesh.segment(NULL, x, NULL, is.bnd, crs=crs)
} else if (is.numeric(x)) { ## Coordinates
ret = inla.mesh.segment(x, NULL, NULL, is.bnd, crs=crs)
} else {
stop("Segment info matrix must be numeric or integer.")
}
} else if (inherits(x, "inla.mesh.segment")) {
## Override x$is.bnd:
ret = inla.mesh.segment(x$loc, x$idx, x$grp, is.bnd, x$crs)
} else if (inherits(x, c("SpatialPoints", "SpatialPointsDataFrame",
"Line", "Lines",
"SpatialLines", "SpatialLinesDataFrame",
"Polygon", "Polygons",
"SpatialPolygons",
"SpatialPolygonsDataFrame"))) {
x <- as.inla.mesh.segment(x)
## Override x$is.bnd:
ret <- inla.mesh.segment(x$loc, x$idx, x$grp, is.bnd, x$crs)
} else if (!is.null(x)) {
inla.require.inherits(NULL,
c("matrix", "inla.mesh.segment"),
"Segment info")
} else {
ret = NULL
}
ret
}
##################################################
homogenise.segm.grp <- function(input) {
grp.idx = 0L
for (k in 1:length(input)) if (!is.null(input[[k]])) {
inla.require.inherits(input[[k]],
"inla.mesh.segment",
"Segment info list members ")
if (is.null(input[[k]]$grp)) {
grp.idx = grp.idx+1L
input[[k]] = (inla.mesh.segment(input[[k]]$loc,
input[[k]]$idx,
grp.idx,
input[[k]]$is.bnd,
input[[k]]$crs))
} else {
grp.idx = max(grp.idx, input[[k]]$grp, na.rm=TRUE)
}
}
return(input)
}
##################################################
join.segm.input <- function(segm, segm.offset=0L, loc.offset=0L,
crs=NULL)
{
loc = NULL
bnd = list(loc=NULL, idx = matrix(,0,2), grp = matrix(,0,1))
int = list(loc=NULL, idx = matrix(,0,2), grp = matrix(,0,1))
storage.mode(bnd$idx) <- "integer"
storage.mode(bnd$grp) <- "integer"
storage.mode(int$idx) <- "integer"
storage.mode(int$grp) <- "integer"
for (k in seq_along(segm)) {
if (!is.null(segm[[k]])) {
inla.require.inherits(segm[[k]],
"inla.mesh.segment",
"Segment info list members ")
if (!is.null(segm[[k]]$loc)) {
local.loc <- inla.spTransform(segm[[k]], crs,
passthrough=TRUE)$loc
extra.loc.n = nrow(local.loc)
idx.offset = segm.offset
segm.offset = segm.offset + extra.loc.n
loc = (inla.ifelse(is.null(loc),
local.loc,
rbind(loc, local.loc)))
} else {
idx.offset = loc.offset
}
if (segm[[k]]$is.bnd) {
bnd$idx = rbind(bnd$idx, segm[[k]]$idx+idx.offset)
bnd$grp = rbind(bnd$grp, segm[[k]]$grp)
} else {
int$idx = rbind(int$idx, segm[[k]]$idx+idx.offset)
int$grp = rbind(int$grp, segm[[k]]$grp)
}
}
}
if (nrow(bnd$idx)==0)
bnd = NULL
if (nrow(int$idx)==0)
int = NULL
return(list(loc = loc,
bnd = (inla.ifelse(is.null(bnd),
NULL,
inla.mesh.segment(bnd$loc,
bnd$idx,
bnd$grp,
TRUE,
crs))),
int = (inla.ifelse(is.null(int),
NULL,
inla.mesh.segment(int$loc,
int$idx,
int$grp,
FALSE,
crs)))))
}
###########################
segm = (c(lapply(inla.ifelse(inherits(boundary, "list"),
boundary, list(boundary)),
function(x){homogenise.segm.input(x, TRUE, crs=crs)}),
lapply(inla.ifelse(inherits(interior, "list"),
interior, list(interior)),
function(x){homogenise.segm.input(x, FALSE, crs=crs)})))
segm = homogenise.segm.grp(segm)
join.segm.input(segm=segm, segm.offset=segm.offset, loc.offset=loc.offset,
crs=crs)
}
################################
##
## Old code. Filtering is now done in fmesher itself.
## Retained for now so that we can check if the results are the same.
##
inla.mesh.filter.locations <- function(loc, cutoff)
{
## Map locations to nodes, avoiding near-duplicates.
loc.n = nrow(loc)
loc.dim = ncol(loc)
loc.is.na = (rowSums(is.na(loc))>0)
if (sum(loc.is.na)>0)
stop("NAs in locations not yet supported.")
node.coord = matrix(nrow=loc.n, ncol=loc.dim)
map.loc.to.node = rep(0L, nrow=loc.n)
excluded = c()
loc.i = 1L
node.i.max = 1L
node.coord[node.i.max,] = loc[loc.i,]
map.loc.to.node[[loc.i]] = node.i.max
for (loc.i in 2:loc.n) {
loc.to.node.dist =
sqrt(rowSums((as.matrix(rep(1, node.i.max)) %*%
loc[loc.i,, drop=FALSE] -
node.coord[1:node.i.max,, drop=FALSE])^2))
if (min(loc.to.node.dist) > cutoff) {
node.i.max = node.i.max+1L
node.coord[node.i.max,] = loc[loc.i,, drop=FALSE]
map.loc.to.node[[loc.i]] = node.i.max
} else {
excluded = c(excluded, loc.i)
}
}
## Remove excess nodes.
node.coord = node.coord[1:node.i.max,]
## Identify nearest nodes for excluded locations.
for (loc.i in excluded) {
loc.to.node.dist =
sqrt(rowSums((as.matrix(rep(1, node.i.max)) %*%
loc[loc.i,, drop=FALSE] -
node.coord)^2))
node.i = which.min(loc.to.node.dist);
map.loc.to.node[[loc.i]] = node.i
}
return(list(loc = node.coord, node.idx = map.loc.to.node))
}
############################
inla.mesh <- function(...)
{
args = list(...)
if (length(args)>0) {
if (inherits(args[[1]],"inla.mesh")) {
warning("'inla.mesh(mesh, ...)' is deprecated. Use 'inla.mesh.query(mesh, ...)' instead.")
return(inla.mesh.query(...))
}
}
warning("'inla.mesh(...)' is deprecated. Use 'inla.mesh.create(...)' instead.")
return(inla.mesh.create(...))
}
inla.mesh.create <- function(loc=NULL, tv=NULL,
boundary=NULL, interior=NULL,
extend = (missing(tv) || is.null(tv)),
refine=FALSE,
lattice=NULL,
globe=NULL,
cutoff = 1e-12,
plot.delay = NULL,
data.dir,
keep = (!missing(data.dir) && !is.null(data.dir)),
timings = FALSE,
quality.spec=NULL,
crs=NULL)
{
if (!timings) {
system.time <- function(expr) {
expr
structure(c(0,0,0,0,0), class = "proc_time")
}
}
time.pre = system.time({ ## Pre-processing timing start
if (!is.null(crs) &&
identical(inla.as.list.CRS(crs)$proj, "geocent")) {
## Build all geocentric meshes on a sphere, and transform afterwards,
## to allow general geoids.
crs.target <- crs
crs <- inla.CRS("sphere")
}
if (!(missing(loc) || is.null(loc))) {
## Handle loc given as SpatialPoints or SpatialPointsDataFrame object
if (inherits(loc, "SpatialPoints") ||
inherits(loc, "SpatialPointsDataFrame")) {
loc <- inla.spTransform(coordinates(loc),
CRS(proj4string(loc)),
crs,
passthrough=TRUE)
}
if (!is.matrix(loc)) {
loc = as.matrix(loc)
}
if (!is.double(loc)) {
storage.mode(loc) = "double"
}
}
loc.n = max(0L,nrow(loc))
if (is.logical(extend) && extend) extend = list()
if (is.logical(refine) && refine) refine = list()
if (missing(lattice) || is.null(lattice)) {
lattice = list(loc=NULL, segm=NULL)
lattice.n = 0L
} else {
inla.require.inherits(lattice, "inla.mesh.lattice", "'lattice'")
if (!is.null(tv)) {
warning("Both 'lattice' and 'tv' specified. Ignoring 'tv'.")
tv = NULL
}
if (!inherits(extend, "list")) {
boundary = (c(inla.ifelse(inherits(boundary, "list"),
boundary, list(boundary)),
list(lattice$segm)))
}
if (!is.null(lattice[["crs"]])) {
lattice$loc <- inla.spTransform(lattice$loc,
lattice$crs,
crs,
passthrough=TRUE)
}
if (!is.matrix(lattice$loc)) {
lattice$loc = as.matrix(lattice$loc)
}
if (!is.double(lattice$loc)) {
storage.mode(lattice$loc) = "double"
}
}
lattice.n = max(0L,nrow(lattice$loc))
segm = (inla.mesh.parse.segm.input(boundary,
interior,
loc.n,
0L,
crs=crs))
segm.n = max(0,nrow(segm$loc))
## Run parse again now that we know where the indices should point:
segm = (inla.mesh.parse.segm.input(boundary,
interior,
0L,
segm.n+lattice.n,
crs=crs))
loc0 = rbind(segm$loc, lattice$loc, loc)
if ((!is.null(loc0)) && (nrow(loc0)>0))
idx0 = 1:nrow(loc0)
else
idx0 = c()
if (is.null(quality.spec)) {
quality <- NULL
} else {
quality <- rep(NA, segm.n+lattice.n+loc.n)
## Order must be same as for loc0 above.
if (!is.null(quality.spec$segm)) {
quality[seq_len(segm.n)] <- quality.spec$segm
}
if (!is.null(quality.spec$lattice)) {
quality[segm.n+seq_len(lattice.n)] <- quality.spec$lattice
}
if (!is.null(quality.spec$loc)) {
quality[segm.n+lattice.n+seq_len(loc.n)] <- quality.spec$loc
}
## NA:s will be replaced with max.edge settings below.
}
## Where to put the files?
if (keep) {
if (missing(data.dir)) {
data.dir = inla.fmesher.make.dir("inla.mesh.data")
}
prefix = paste(data.dir, "/mesh.", sep="")
keep.dir = TRUE
} else {
if (missing(data.dir)) {
prefix = inla.tempfile(pattern="fmesher", tmpdir=tempdir())
prefix = paste(prefix, ".", sep="")
keep.dir = TRUE
} else {
data.dir = inla.fmesher.make.dir(data.dir)
prefix = paste(data.dir, "/mesh.", sep="")
## We should not try to delete the session tempdir...
keep.dir = identical(tempdir(), dirname(prefix))
}
}
prefix = inla.fmesher.make.prefix(NULL, prefix)
if (is.null(loc0)) {
all.args = ""
} else {
loc.file = fmesher.write(inla.affirm.double(loc0), prefix, "input.s")
all.args = "--input=input.s"
if (!is.null(tv)) {
fmesher.write(inla.affirm.integer(tv)-1L, prefix, "input.tv")
all.args = paste(all.args, ",input.tv", sep="")
}
}
if (!missing(cutoff)) {
all.args = paste(all.args, " --cutoff=", cutoff, sep="")
}
if (!is.null(segm$bnd)) {
fmesher.write(inla.affirm.integer(segm$bnd$idx)-1L, prefix, "input.segm.bnd.idx")
fmesher.write(inla.affirm.integer(segm$bnd$grp), prefix, "input.segm.bnd.grp")
all.args = paste(all.args," --boundary=input.segm.bnd.idx")
all.args = paste(all.args," --boundarygrp=input.segm.bnd.grp")
}
if (!is.null(segm$int)) {
fmesher.write(inla.affirm.integer(segm$int$idx)-1L, prefix, "input.segm.int.idx")
fmesher.write(inla.affirm.integer(segm$int$grp), prefix, "input.segm.int.grp")
all.args = paste(all.args," --interior=input.segm.int.idx")
all.args = paste(all.args," --interiorgrp=input.segm.int.grp")
}
if (!missing(globe) && !is.null(globe)) {
all.args = paste(all.args, " --globe=", globe, sep="")
}
if (inherits(extend,"list")) {
cet = c(0,0)
cet[1] = inla.ifelse(is.null(extend$n), 16, extend$n)
cet[2] = inla.ifelse(is.null(extend$offset), -0.1, extend$offset)
all.args = (paste(all.args," --cet=",
cet[1],",", cet[2], sep=""))
}
if (inherits(refine,"list")) {
rcdt = c(0,0,0)
if (!missing(globe) && !is.null(globe)) {
max.edge.default <- pi
} else {
## Multiply by to ensure cover S2 domains
max.edge.default <- inla.diameter(loc0) * 2
}
if ((inherits(extend,"list")) && (!is.null(extend$offset))) {
max.edge.default = (max.edge.default +
max(0,2*extend$offset))
max.edge.default = (max.edge.default *
(1+max(0,-2*extend$offset)))
}
rcdt[1] = inla.ifelse(is.null(refine$min.angle), 21, refine$min.angle)
rcdt[2] = (inla.ifelse(is.null(refine$max.edge) ||
is.na(refine$max.edge),
max.edge.default,
refine$max.edge))
rcdt[3] = (inla.ifelse(is.null(refine$max.edge) ||
is.na(refine$max.edge),
max.edge.default,
refine$max.edge))
rcdt[2] = (inla.ifelse(is.null(refine$max.edge.extra) ||
is.na(refine$max.edge.extra),
rcdt[2], refine$max.edge.extra))
rcdt[3] = (inla.ifelse(is.null(refine$max.edge.data) ||
is.na(refine$max.edge.date),
rcdt[3], refine$max.edge.data))
all.args = (paste(all.args," --rcdt=",
rcdt[1],",", rcdt[2],",", rcdt[3], sep=""))
if (!is.null(refine[["max.n.strict"]]) &&
!is.na(refine$max.n.strict)) {
rcdt_max_n0 <- refine$max.n.strict
} else {
rcdt_max_n0 <- -1
}
if (!is.null(refine[["max.n"]]) &&
!is.na(refine$max.n)) {
rcdt_max_n1 <- refine$max.n
} else {
rcdt_max_n1 <- -1
}
all.args = (paste(all.args,
" --max_n0=", rcdt_max_n0,
" --max_n1=", rcdt_max_n1,
sep=""))
is.refined = TRUE
if (!is.null(quality)) {
quality[is.na(quality)] <- rcdt[3]
quality.file <-
fmesher.write(inla.affirm.double(as.matrix(quality)),
prefix, "input.quality")
all.args <- (paste(all.args,
" --input=input.quality",
" --quality=input.quality",
sep=""))
}
} else {
is.refined = FALSE
}
if (!is.null(plot.delay)) {
all.args = paste(all.args," --x11=", plot.delay, sep="")
}
all.args = paste(all.args, inla.getOption("fmesher.arg"))
}) ## Pre-processing timing end
## Call fmesher:
time.fmesher = system.time({
echoc = (inla.fmesher.call(all.args=all.args,
prefix=prefix))
})
time.post = system.time({ ## Post-processing timing start
## Read the mesh:
manifold = 1L+fmesher.read(prefix, "manifold")
manifold = list("M", "R2", "S2")[[manifold]]
loc = fmesher.read(prefix, "s")
graph = (list(tv = 1L+fmesher.read(prefix, "tv"),
vt = 1L+fmesher.read(prefix, "vt"),
tt = 1L+fmesher.read(prefix, "tt"),
tti = 1L+fmesher.read(prefix, "tti"),
vv = fmesher.read(prefix, "vv")))
graph$tv[graph$tv==0L] = NA
graph$vt[graph$vt==0L] = NA
graph$tt[graph$tt==0L] = NA
graph$tti[graph$tti==0L] = NA
## Read the vertex input/output mapping:
idx.all = 1L+fmesher.read(prefix, "idx")
idx = (list(loc = (inla.ifelse(loc.n>0,
idx.all[idx0[segm.n+lattice.n+(1:loc.n)]],
NULL)),
lattice = (inla.ifelse(lattice.n>0,
idx.all[idx0[segm.n+(1:lattice.n)]],
NULL)),
segm = (inla.ifelse(segm.n>0,
idx.all[idx0[(1:segm.n)]],
NULL))))
if (!is.null(idx$loc)) idx$loc[idx$loc == 0L] = NA
if (!is.null(idx$lattice)) idx$lattice[idx$lattice == 0L] = NA
if (!is.null(idx$segm)) idx$segm[idx$segm == 0L] = NA
## Read constraint segment information:
segm.bnd = (inla.mesh.segment(NULL,
1L+fmesher.read(prefix, "segm.bnd.idx"),
fmesher.read(prefix, "segm.bnd.grp"),
TRUE))
segm.int = (inla.mesh.segment(NULL,
1L+fmesher.read(prefix, "segm.int.idx"),
fmesher.read(prefix, "segm.int.grp"),
FALSE))
## Remap indices to remove unused vertices
used = !is.na(graph$vt)
if (!all(used)) {
used = which(used)
idx.map = rep(NA, nrow(loc))
idx.map[used] = seq_len(length(used))
loc = loc[used,,drop=FALSE]
graph$tv = matrix(idx.map[as.vector(graph$tv)], nrow(graph$tv), 3)
graph$vt = graph$vt[used,,drop=FALSE]
## graph$tt ## No change needed
## graph$tti ## No change needed
graph$vv = graph$vv[used, used, drop=FALSE]
if (!is.null(idx$loc)) idx$loc = idx.map[idx$loc]
if (!is.null(idx$lattice)) idx$lattice = idx.map[idx$lattice]
if (!is.null(idx$segm)) idx$segm = idx.map[idx$segm]
segm.bnd$idx = matrix(idx.map[segm.bnd$idx], nrow(segm.bnd$idx), 2)
segm.int$idx = matrix(idx.map[segm.int$idx], nrow(segm.int$idx), 2)
segm.bnd$idx[segm.bnd$idx == 0L] = NA
segm.int$idx[segm.int$idx == 0L] = NA
}
if (!keep)
unlink(paste(prefix, "*", sep=""), recursive=FALSE)
if (!keep.dir) {
unlink(dirname(prefix), recursive=TRUE)
}
if (!is.null(crs) &&
identical(inla.as.list.CRS(crs)$proj, "geocent")) {
if (!inla.identical.CRS(crs, crs.target)) {
## Target is a non-spherical geoid
loc <- inla.spTransform(loc, crs, crs.target)
crs <- crs.target
}
}
}) ## Post-processing timing end
time.object = system.time({ ## Object construction timing start
mesh = (list(meta = (list(call=match.call(),
fmesher.args = all.args,
time = (rbind(pre = time.pre,
fmesher = time.fmesher,
post = time.post)),
prefix = prefix,
is.refined = is.refined)),
manifold = manifold,
n = nrow(loc),
loc = loc,
graph = graph,
segm = list(bnd=segm.bnd, int=segm.int),
idx = idx,
crs = crs))
class(mesh) <- "inla.mesh"
}) ## Object construction timing end
## Entire function timing
time.total = time.pre+time.fmesher+time.post+time.object
mesh$meta$time = (rbind(mesh$meta$time,
object=time.object,
total=time.total))
return(invisible(mesh))
}
inla.mesh.extract.segments <- function(mesh.loc,
mesh.idx,
mesh.grp,
grp=NULL,
is.bnd,
crs=NULL)
{
segments = list()
if (nrow(mesh.idx)>0) {
if (is.null(grp)) {
grp = unique(sort(mesh.grp))
}
for (g in grp) {
extract = (mesh.grp==g)
segments =
c(segments,
list(inla.mesh.segment(mesh.loc,
idx=mesh.idx[extract,,drop=FALSE],
grp=mesh.grp[extract,drop=FALSE],
is.bnd=is.bnd,
crs=crs)) )
}
}
if (length(segments)>0)
return(segments)
else
return(NULL)
}
inla.mesh.boundary <- function(mesh, grp=NULL)
{
inla.require.inherits(mesh, "inla.mesh", "'mesh'")
return(inla.mesh.extract.segments(mesh$loc,
mesh$segm$bnd$idx,
mesh$segm$bnd$grp,
grp,
TRUE,
mesh$crs))
}
inla.mesh.interior <- function(mesh, grp=NULL)
{
inla.require.inherits(mesh, "inla.mesh", "'mesh'")
return(inla.mesh.extract.segments(mesh$loc,
mesh$segm$int$idx,
mesh$segm$int$grp,
grp,
FALSE,
mesh$crs))
}
## Generate nice triangulation, with an inner domain strictly enclosed
## by an outer domain
## The inner and outer domains can have different quality parameters
## At least one of loc, loc.domain, boundary[[1]], boundary[[2]], interior
## must be non-NULL
## For more complicated multi-step meshings, study the code and write your own.
inla.mesh.2d <-
function(loc=NULL, ## Points to include in final triangulation
loc.domain=NULL, ## Points that determine the automatic domain
offset=NULL, ## Size of automatic extensions
n=NULL, ## Sides of automatic extension polygons
boundary=NULL, ## User-specified domains (list of length 2)
interior=NULL, ## User-specified constraints for the inner domain
max.edge=NULL,
min.angle=NULL, ## Angle constraint for the entire domain
cutoff=1e-12, ## Only add input points further apart than this
max.n.strict=NULL,
max.n=NULL,
plot.delay=NULL,
crs=NULL) ## Coordinate Reference System
## plot.delay: Do plotting.
## NULL --> No plotting
## <0 --> Intermediate meshes displayed at the end
## >0 --> Dynamical fmesher plotting
{
###########################
unify.one.segm <- function(segm, crs=NULL) {
if (inherits(segm, "inla.mesh.segment")) {
segm <- inla.spTransform(segm, crs, passthrough=TRUE)
} else if (inherits(segm, "matrix")) {
segm <- inla.mesh.segment(loc=segm, crs=crs)
} else {
segm <- inla.spTransform(as.inla.mesh.segment(segm), crs,
passthrough=TRUE)
}
if (ncol(segm$loc)==2) {
segm$loc = cbind(segm$loc, 0.0)
}
segm
}
unify.segm.coords <- function(segm, crs=NULL) {
if (!is.null(segm)) {
if (inherits(segm, "list")) {
for (j in seq_along(segm)) {
segm[[j]] <- unify.one.segm(segm[[j]], crs=crs)
}
} else {
segm <- unify.one.segm(segm, crs=crs)
}
}
segm
}
###########################
if ((missing(max.edge) || is.null(max.edge)) &&
(missing(max.n.strict) || is.null(max.n.strict)) &&
(missing(max.n) || is.null(max.n))) {
stop("At least one of max.edge, max.n.strict, and max.n must be specified")
}
if (!is.null(crs)) {
issphere <- inla.identical.CRS(crs, inla.CRS("sphere"))
isgeocentric <- identical(inla.as.list.CRS(crs)[["proj"]], "geocent")
if (isgeocentric) {
crs.target <- crs
crs <- inla.CRS("sphere")
}
}
## Handle loc given as SpatialPoints or SpatialPointsDataFrame object
if (!(missing(loc) || is.null(loc)) &&
(inherits(loc, "SpatialPoints") ||
inherits(loc, "SpatialPointsDataFrame"))) {
loc = inla.spTransform(coordinates(loc),
CRS(proj4string(loc)),
crs,
passthrough=TRUE)
}
if (!(missing(loc.domain) || is.null(loc.domain)) &&
(inherits(loc.domain, "SpatialPoints") ||
inherits(loc.domain, "SpatialPointsDataFrame"))) {
loc.domain = inla.spTransform(coordinates(loc.domain),
CRS(proj4string(loc.domain)),
crs,
passthrough=TRUE)
}
if (missing(loc) || is.null(loc)) {
loc = matrix(c(0.0), 0, 3)
} else if (!is.matrix(loc)) {
loc = as.matrix(loc)
}
if (missing(loc.domain) || is.null(loc.domain)) {
loc.domain = loc
} else if (!is.matrix(loc.domain)) {
loc.domain = as.matrix(loc.domain)
}
if (missing(boundary)) {
boundary = list(NULL)
} else {
if (!inherits(boundary, "list"))
boundary = list(boundary)
}
if (missing(interior))
interior = NULL
if (missing(offset) || is.null(offset)) {
if (length(boundary) < 2)
offset = -0.05
else
offset = c(-0.05, -0.15)
}
if (missing(n) || is.null(n))
n = c(8)
if (missing(max.edge) || is.null(max.edge))
max.edge = c(NA)
if (missing(min.angle) || is.null(min.angle))
min.angle = c(21)
if (missing(max.n.strict) || is.null(max.n.strict))
max.n.strict = c(NA)
if (missing(max.n) || is.null(max.n))
max.n = c(NA)
if (missing(cutoff) || is.null(cutoff))
cutoff = 1e-12
if (missing(plot.delay) || is.null(plot.delay))
plot.delay = NULL
num.layers =
max(c(length(boundary), length(offset), length(n),
length(min.angle), length(max.edge),
length(max.n.strict), length(max.n)))
if (num.layers > 2) {
warning(paste("num.layers=", num.layers, " > 2 detected. ",
"Excess information ignored.", sep=""))
num.layers = 2
}
if (length(boundary) < num.layers)
boundary = c(boundary, list(NULL))
if (length(min.angle) < num.layers)
min.angle = c(min.angle, min.angle)
if (length(max.n.strict) < num.layers)
max.n0 = c(max.n.strict, max.n.strict)
if (length(max.n) < num.layers)
max.n = c(max.n, max.n)
if (length(max.edge) < num.layers)
max.edge = c(max.edge, max.edge)
if (length(offset) < num.layers)
offset = c(offset, -0.15)
if (length(n) < num.layers)
n = c(n, 16)
if (length(n) < num.layers)
n = c(n, 16)
## Unify the dimensionality of the point input.
if (!is.null(loc) && (ncol(loc)==2))
loc = cbind(loc, 0.0)
if (!is.null(loc.domain) && (ncol(loc.domain)==2))
loc.domain = cbind(loc.domain, 0.0)
## Unify the dimensionality of the boundary&interior segments input
## and optionally transform coordinates.
for (k in seq_len(num.layers)) {
if (!is.null(boundary[[k]])) {
boundary[[k]] <- unify.segm.coords(boundary[[k]], crs=crs)
}
}
interior <- unify.segm.coords(interior, crs=crs)
## Triangulate to get inner domain boundary
## Constraints included only to get proper domain extent
## First, attach the loc points to the domain definition set
if (!is.null(loc) && !is.null(loc.domain)) {
loc.domain = rbind(loc.domain, loc)
}
mesh1 =
inla.mesh.create(loc=loc.domain,
boundary=boundary[[1]],
interior=interior,
cutoff=cutoff,
extend=list(n=n[1], offset=offset[1]),
refine=FALSE,
plot.delay=plot.delay,
crs=crs)
## Save the resulting boundary
boundary1 = inla.mesh.boundary(mesh1)
interior1 = inla.mesh.interior(mesh1)
if (!is.null(plot.delay) && (plot.delay<0)) {
inla.dev.new()
plot(mesh1)
}
## Triangulate inner domain
mesh2 =
inla.mesh.create(loc=loc,
boundary=boundary1,
interior=interior1,
cutoff=cutoff,
extend=FALSE, ## Should have no effect
refine=
list(min.angle=min.angle[1],
max.edge=max.edge[1],
max.edge.extra=max.edge[1],
max.n.strict=max.n.strict[1],
max.n=max.n[1]),
plot.delay=plot.delay,
crs=crs)
boundary2 = inla.mesh.boundary(mesh2)
interior2 = inla.mesh.interior(mesh2)
if (!is.null(plot.delay) && (plot.delay<0)) {
inla.dev.new()
plot(mesh2)
}
if (num.layers == 1) {
if (!is.null(crs) && isgeocentric && !issphere) {
mesh2$loc <- inla.spTransform(mesh2$loc, mesh2$crs, crs.target)
mesh2$crs <- crs.target
}
return(invisible(mesh2))
}
## Triangulate inner+outer domain
mesh3 =
inla.mesh.create(loc=rbind(loc, mesh2$loc),
boundary=boundary[[2]],
interior=c(boundary2, interior2),
cutoff=cutoff,
extend=list(n=n[2], offset=offset[2]),
refine=
list(min.angle=min.angle[2],
max.edge=max.edge[2],
max.edge.extra=max.edge[2],
max.n.strict=mesh2$n + max.n.strict[2],
max.n=mesh2$n + max.n[2]),
plot.delay=plot.delay,
crs=crs)
## Hide generated points, to match regular inla.mesh.create output
mesh3$idx$loc = mesh3$idx$loc[seq_len(nrow(loc))]
## Obtain the corresponding segm indices.
segm.loc = matrix(0.0, 0, 3)
for (k in seq_along(boundary)) {
if (!is.null(boundary[[k]])) {
segm.loc = rbind(segm.loc, boundary[[k]]$loc)
}
}
for (k in seq_along(interior)) {
if (!is.null(interior[[k]])) {
segm.loc = rbind(segm.loc, interior[[k]]$loc)
}
}
if (nrow(segm.loc) > 0) {
proj = inla.mesh.project(mesh3, loc=segm.loc)
mesh3$idx$segm = rep(NA, nrow(segm.loc))
if (any(proj$ok)) {
t.idx <- proj$t[proj$ok]
tv.idx <- max.col(proj$bary[proj$ok,,drop=FALSE],
ties.method="first")
mesh3$idx$segm[proj$ok] <-
mesh3$graph$tv[t.idx + nrow(mesh3$graph$tv)*(tv.idx-1)]
}
} else {
mesh3$idx$segm = NULL
}
if (!is.null(crs) && isgeocentric && !issphere) {
mesh3$loc <- inla.spTransform(mesh3$loc, mesh3$crs, crs.target)
mesh3$crs <- crs.target
}
if (!is.null(plot.delay) && (plot.delay<0)) {
inla.dev.new()
plot(mesh3)
}
return(invisible(mesh3))
}
## Support for legacy code:
inla.mesh.create.helper <- function(points=NULL, points.domain=NULL, ...)
{
return(invisible(inla.mesh.2d(loc=points, loc.domain=points.domain, ...)))
}
## Example:
##if (FALSE) {
## mesh =
## inla.mesh.create.helper(loc=matrix(runif(20),10,2)*200,
## n=c(8,16),
## offset=c(10,140),
## max.edge=c(25,1000),
## min.angle=26,
## cutoff=0,
## plot.delay=-1
## )
##}
##
inla.delaunay <- function(loc, ...)
{
## Handle loc given as SpatialPoints or SpatialPointsDataFrame object
if (!(missing(loc) || is.null(loc)) &&
(inherits(loc, "SpatialPoints") ||
inherits(loc, "SpatialPointsDataFrame"))) {
crs = CRS(proj4string(loc))
loc = coordinates(loc)
} else {
crs = NULL
}
hull = chull(loc[,1],loc[,2])
bnd = inla.mesh.segment(loc=loc[hull[length(hull):1],],is.bnd=TRUE)
mesh =
inla.mesh.create(loc=loc,
boundary=bnd,
extend=list(n=3),
refine=FALSE,
crs=crs,
...)
return(invisible(mesh))
}
inla.mesh.query <- function(mesh, ...)
{
inla.require.inherits(mesh, "inla.mesh", "'mesh'")
not.known <- function(mesh, queryname)
{
stop(paste("Query '", queryname,
"' unknown.", sep=""))
}
not.implemented <- function(mesh, queryname)
{
stop(paste("Query '", queryname,
"' not implemented for inla.mesh.", sep=""))
## stop(paste("Query '", queryname,
## "' not implemented for inla.mesh for mesh type '",
## model$type, "'.", sep=""))
}
result = list()
queries = inla.parse.queries(...)
if (length(queries)==0L)
return(result)
for (query.idx in 1:length(queries)) {
query = names(queries)[query.idx]
param = queries[[query.idx]]
answer = NULL
query = (match.arg(query, c("tt.neighbours",
"vt.neighbours"
)))
if (identical(query, "tt.neighbours")) {
if (is.null(param))
param = list(c(1))
if (length(param)<2)
param = c(as.list(param), list(c(1)))
nT = nrow(mesh$graph$tv)
i = rep(1:nT,3)
j = as.vector(mesh$graph$tt)
i = i[!is.na(j)]
j = j[!is.na(j)]
tt = sparseMatrix(i=i, j=j, x=rep(1,length(i)), dims=c(nT,nT))
answer = (sparseMatrix(i=param[[1]],
j=rep(1,length(param[[1]])),
x=1,
dims=c(nT,1)))
tt0 = (answer == 0.5)*1
tt1 = answer
order = 0
while (order<min(param[[2]])) {
order = order+1
tt0 = tt1
tt1 = answer
answer = ((((tt %*% answer) > 0) - tt1 - tt0) > 0)
}
while (order<max(param[[2]])) {
order = order+1
answer = ((((tt %*% answer) > 0) - tt1 - tt0) > 0)
}
## not.implemented(mesh,query)
} else if (identical(query, "vt.neighbours")) {
if (is.null(param))
param = list(1)
if (length(param)<2)
param = c(as.list(param), list(c(1)))
nV = nrow(mesh$loc)
nT = nrow(mesh$graph$tv)
i = rep(1:nT,3)
j = as.vector(mesh$graph$tv)
# i = i[!is.na(j)]
# j = j[!is.na(j)]
tv = sparseMatrix(i=i, j=j, x=rep(1,length(i)), dims=c(nT,nV))
vv = (sparseMatrix(i=param[[1]],
j=rep(1,length(param[[1]])),
x=1,
dims=c(nV,1)))
vt = (tv %*% vv ) > 0
vv0 = (vv == 0.5)*1
vv1 = vv
vt0 = (tv %*% vv0 ) > 0
vt1 = vt
order = 0
while (order<min(param[[2]])) {
order = order+1
vv0 = vv1
vv1 = vv
vv = ((((mesh$graph$vv %*% vv) > 0) - vv1 - vv0) > 0)
vt0 = vt1
vt1 = vt
vt = (((tv %*% vv) > 0) - vt1 - vt0 ) > 0
}
while (order<max(param[[2]])) {
order = order+1
vv = ((((mesh$graph$vv %*% vv) > 0) - vv1 - vv0) > 0)
vt = (((((tv %*% vv) > 0) - vt1 - vt0 ) > 0 ) + vt) > 0
}
answer = vt
## not.implemented(mesh,query)
} else if (!identical(query, "")) {
not.known(mesh,query)
}
## Expand the result list:
result[query.idx] = list(NULL)
names(result)[query.idx] = query
## Set the answer:
if (!is.null(answer))
result[[query.idx]] = answer
}
return(result)
}
summary.inla.mesh <- function(object, verbose=FALSE, ...)
{
x=object
## provides a summary for a mesh object
inla.require.inherits(x, "inla.mesh", "'x'")
ret = list(verbose=verbose)
if (verbose) {
ret = (c(ret, list(call=x$meta$call,
fmesher.args=x$meta$fmesher.args,
prefix=x$meta$prefix,
time = x$meta$time,
is.refined = x$meta$is.refined)))
}
ret = (c(ret, list(manifold=x$manifold,
nV=x$n,
nT=nrow(x$graph$tv),
xlim=range(x$loc[,1]),
ylim=range(x$loc[,2]),
zlim=range(x$loc[,3]))))
if (is.na(inla.CRSargs(x$crs))) {
ret <- c(ret, list(crs="N/A"))
} else {
ret <- c(ret, list(crs=inla.CRSargs(x$crs)))
}
my.segm <- function(x) {
if (is.null(x))
return(list(n=0, grps=NULL))
n = max(0, nrow(x$idx))
if (max(0, length(unique(x$grp)))>0) {
grps = unique(x$grp)
} else {
grps = NULL
}
return(list(n=n, grps=grps))
}
if(!is.null(x$segm)) {
ret = (c(ret, list(segm.bnd=my.segm(x$segm$bnd),
segm.int=my.segm(x$segm$int))))
} else {
ret = (c(ret, list(segm.bnd=my.segm(NULL),
segm.int=my.segm(NULL))))
}
class(ret) <- "summary.inla.mesh"
return (ret)
}
print.summary.inla.mesh <- function(x, ...)
{
my.print.proc_time <- function (x, ...)
{
if (!is.matrix(x)) {
y = matrix(x,1,5)
} else {
y = x
}
for (k in 1:nrow(x)) {
if (!is.na(y[k,4L]))
y[k,1L] <- y[k,1L] + y[k,4L]
if (!is.na(y[k,5L]))
y[k,2L] <- y[k,2L] + y[k,5L]
}
y <- y[,1L:3L]
colnames(y) <- c(gettext("user"), gettext("system"), gettext("elapsed"))
print(y, ...)
invisible(x)
}
inla.require.inherits(x, "summary.inla.mesh", "'x'")
if (x$verbose) {
cat("\nCall:\n")
print(x$call)
cat("\nfmesher:\t", x$fmesher.args, "\n", sep = "")
cat("prefix:\t\t", x$prefix, "\n", sep = "")
cat("\nTimings:\n")
my.print.proc_time(x$time)
}
cat("\nManifold:\t", x$manifold, "\n", sep="")
cat("CRS:\t", x$crs, "\n", sep="")
if (x$verbose) {
cat("Refined:\t", x$is.refined, "\n", sep="")
}
cat("Vertices:\t", as.character(x$nV), "\n", sep="")
cat("Triangles:\t", as.character(x$nT), "\n", sep="")
my.print.segm <- function(x) {
cat(as.character(x$n))
if (!is.null(x$grps)) {
n = length(x$grps)
cat(" (", n, " group", inla.ifelse(n==1, "", "s"), sep="")
if (n <= 10) {
cat(":", x$grps, sep=" ")
} else {
cat(":", x$grps[1:10], "...", sep=" ")
}
cat(")")
}
cat("\n", sep="")
return(invisible())
}
cat("Boundary segm.:\t")
my.print.segm(x$segm.bnd)
cat("Interior segm.:\t")
my.print.segm(x$segm.int)
cat("xlim:\t", x$xlim[1], " ", x$xlim[2], "\n", sep="")
cat("ylim:\t", x$ylim[1], " ", x$ylim[2], "\n", sep="")
cat("zlim:\t", x$zlim[1], " ", x$zlim[2], "\n", sep="")
cat("\n")
invisible(x)
}
inla.mesh.project <- function(...)
{
UseMethod("inla.mesh.project")
}
inla.mesh.project.inla.mesh <- function(mesh, loc=NULL, field=NULL,
crs=NULL, ...)
{
inla.require.inherits(mesh, "inla.mesh", "'mesh'")
if (!is.null(mesh$crs) &&
identical(inla.as.list.CRS(mesh$crs)[["proj"]], "geocent")) {
crs.sphere <- inla.CRS("sphere")
if (!inla.identical.CRS(mesh$crs, crs.sphere)) {
## Convert the mesh to a perfect sphere.
mesh$loc <- inla.spTransform(mesh$loc, mesh$crs, crs.sphere)
mesh$manifold <- "S2"
mesh$crs <- crs.sphere
}
}
## Handle loc given as SpatialPoints or SpatialPointsDataFrame object
## or an explicit crs
if (!(missing(loc) || is.null(loc))) {
if (inherits(loc, "SpatialPoints") ||
inherits(loc, "SpatialPointsDataFrame")) {
if (is.null(mesh$crs)) {
loc <- coordinates(loc)
} else {
loc = inla.spTransform(coordinates(loc),
CRS(proj4string(loc)),
mesh$crs,
passthrough=FALSE)
}
} else if (!is.null(crs)) {
if (!is.null(mesh$crs)) {
loc = inla.spTransform(loc,
crs,
mesh$crs,
passthrough=FALSE)
}
}
}
if (!missing(field) && !is.null(field)) {
proj <- inla.mesh.projector(mesh, loc=loc, ...)
return(inla.mesh.project(proj, field=field))
}
jj =
which(rowSums(matrix(is.na(as.vector(loc)),
nrow=nrow(loc),
ncol=ncol(loc))) == 0)
smorg = (inla.fmesher.smorg(mesh$loc,
mesh$graph$tv,
points2mesh=loc[jj,,drop=FALSE]))
ti = matrix(0L, nrow(loc), 1)
b = matrix(0, nrow(loc), 3)
ti[jj,1L] = smorg$p2m.t
b[jj,] = smorg$p2m.b
ok = (ti[,1L] > 0L)
ti[ti[,1L] == 0L,1L] = NA
ii = which(ok)
A = (sparseMatrix(dims=c(nrow(loc),mesh$n),
i = rep(ii, 3),
j = as.vector(mesh$graph$tv[ti[ii,1L],]),
x = as.vector(b[ii,]) ))
list(t=ti, bary=b, A=A, ok=ok)
}
inla.mesh.project.inla.mesh.1d <- function(mesh, loc, field=NULL, ...)
{
inla.require.inherits(mesh, "inla.mesh.1d", "'mesh'")
if (!missing(field) && !is.null(field)) {
proj = inla.mesh.projector(mesh, loc, ...)
return(inla.mesh.project(proj, field))
}
A = inla.mesh.1d.A(mesh, loc=loc)
return (list(A=A,
ok=(loc >= mesh$interval[1]) & (loc <= mesh$interval[2])))
}
inla.mesh.project.inla.mesh.projector <-
function(projector, field, ...)
{
inla.require.inherits(projector, "inla.mesh.projector", "'projector'")
if (is.data.frame(field)) {
field = as.matrix(field)
}
if (is.null(dim(field))) {
if (is.null(projector$lattice)) {
data = as.vector(projector$proj$A %*% as.vector(field))
data[!projector$proj$ok] = NA
return(data)
} else {
data = as.vector(projector$proj$A %*% as.vector(field))
data[!projector$proj$ok] = NA
return(matrix(data,
projector$lattice$dims[1],
projector$lattice$dims[2]))
}
} else {
data = projector$proj$A %*% field
data[!projector$proj$ok,] = NA
return(data)
}
}
inla.mesh.projector <- function(...)
{
UseMethod("inla.mesh.projector")
}
inla.mesh.projector.inla.mesh <-
function(mesh,
loc=NULL,
lattice=NULL,
xlim=NULL,
ylim=NULL,
dims=c(100,100),
projection=NULL,
crs=NULL,
...)
{
inla.require.inherits(mesh, "inla.mesh", "'mesh'")
if (missing(loc) || is.null(loc)) {
if (missing(lattice) || is.null(lattice)) {
if (identical(mesh$manifold, "R2") &&
(is.null(mesh$crs) || is.null(crs))) {
units = "default"
lim <- list(xlim=range(mesh$loc[,1]),
ylim=range(mesh$loc[,2]))
} else if (identical(mesh$manifold, "S2") &&
(is.null(mesh$crs) || is.null(crs))) {
projection =
match.arg(projection, c("longlat", "longsinlat",
"mollweide"))
units = projection
lim = inla.mesh.map.lim(loc=mesh$loc, projection=projection)
} else {
lim <- inla.crs.bounds(crs)
if (identical(mesh$manifold, "R2")) {
lim0 <- list(xlim=range(mesh$loc[,1]),
ylim=range(mesh$loc[,2]))
lim$xlim[1] <- max(lim$xlim[1], lim0$xlim[1])
lim$xlim[2] <- min(lim$xlim[2], lim0$xlim[2])
lim$ylim[1] <- max(lim$ylim[1], lim0$ylim[1])
lim$ylim[2] <- min(lim$ylim[2], lim0$ylim[2])
}
}
if (missing(xlim) && is.null(xlim)) {
xlim = lim$xlim
}
if (missing(ylim) && is.null(ylim)) {
ylim = lim$ylim
}
x = seq(xlim[1], xlim[2], length.out=dims[1])
y = seq(ylim[1], ylim[2], length.out=dims[2])
if (is.null(mesh$crs) || is.null(crs)) {
lattice <- inla.mesh.lattice(x=x, y=y, units = units)
} else {
lattice <- inla.mesh.lattice(x=x, y=y, crs=crs)
}
} else {
dims = lattice$dims
x = lattice$x
y = lattice$y
}
if (is.null(mesh$crs) || is.null(lattice$crs)) {
proj = inla.mesh.project(mesh, lattice$loc)
} else {
proj <- inla.mesh.project(mesh,
loc=lattice$loc,
crs=lattice$crs)
}
projector = list(x=x, y=y, lattice=lattice, loc=NULL, proj=proj, crs=crs)
class(projector) = "inla.mesh.projector"
} else {
proj = inla.mesh.project(mesh, loc=loc, crs=crs)
projector = list(x=NULL, y=NULL, lattice=NULL, loc=loc, proj=proj, crs=crs)
class(projector) = "inla.mesh.projector"
}
return (projector)
}
inla.mesh.projector.inla.mesh.1d <-
function(mesh,
loc=NULL,
xlim=mesh$interval,
dims=100,
...)
{
inla.require.inherits(mesh, "inla.mesh.1d", "'mesh'")
if (missing(loc) || is.null(loc)) {
loc = seq(xlim[1], xlim[2], length.out=dims[1])
}
proj = inla.mesh.project(mesh, loc)
projector = list(x=loc, lattice=NULL, loc=loc, proj=proj)
class(projector) = "inla.mesh.projector"
return (projector)
}
inla.internal.make.spline.mesh <-
function(interval, m, degree, boundary, free.clamped)
{
boundary =
match.arg(boundary,
c("neumann", "dirichlet", "free", "cyclic"))
if (degree <= 1) {
n = (switch(boundary,
neumann = m,
dirichlet = m+2,
free = m,
cyclic = m+1))
if (n<2) {
n = 2
degree = 0
boundary = "c"
}
} else {
stopifnot(degree==2)
n = (switch(boundary,
neumann = m+1,
dirichlet = m+1,
free = m-1,
cyclic = m))
if (boundary=="free") {
if (m <= 1) {
n = 2
degree = 0
boundary = "c"
} else if (m == 2) {
n = 2
degree = 1
}
} else if (boundary=="cyclic") {
if (m <= 1) {
n = 2
degree = 0
}
}
}
return(inla.mesh.1d(seq(interval[1], interval[2], length=n),
degree=degree,
boundary=boundary,
free.clamped=free.clamped))
}
inla.mesh.basis <- function(mesh,
type="b.spline",
n=3,
degree=2,
knot.placement="uniform.area",
rot.inv=TRUE,
boundary="free",
free.clamped=TRUE,
...)
{
inla.require.inherits(mesh, c("inla.mesh", "inla.mesh.1d"), "'mesh'")
type = match.arg(type, c("b.spline", "sph.harm"))
knot.placement = (match.arg(knot.placement,
c("uniform.area",
"uniform.latitude")))
if (identical(type, "b.spline")) {
if (identical(mesh$manifold, "R1") || identical(mesh$manifold, "S1")) {
mesh1 =
inla.internal.make.spline.mesh(mesh$interval, n, degree,
boundary, free.clamped)
basis = inla.mesh.1d.A(mesh1, mesh$loc)
} else if (identical(mesh$manifold, "R2")) {
if (length(n) == 2) {
if (length(degree)==1) {
degree = rep(degree, 2)
}
if (length(boundary)==1) {
boundary = rep(boundary, 2)
}
mesh1x =
inla.internal.make.spline.mesh(range(mesh$loc[,1]),
n[1], degree[1],
boundary[1], free.clamped)
mesh1y =
inla.internal.make.spline.mesh(range(mesh$loc[,2]),
n[2], degree[2],
boundary[2], free.clamped)
basis =
inla.row.kron(inla.mesh.1d.A(mesh1y, mesh$loc[,2]),
inla.mesh.1d.A(mesh1x, mesh$loc[,1]))
} else {
warning("Old style call for an R2 mesh detected.\n Please supply a 2-element n-vector instead.")
long = ((mesh$loc[,1]-min(mesh$loc[,1]))/
diff(range(mesh$loc[,1])))*90*pi/180
sinlat = (((mesh$loc[,2]-min(mesh$loc[,2]))/
diff(range(mesh$loc[,2])))*2-1)
coslat = sapply(sinlat, function(x) sqrt(max(0.0,1.0-x^2)))
loc = (matrix(c(cos(long)*coslat,
sin(long)*coslat,
sinlat), mesh$n, 3))
knots = 0
degree = max(0L, min(n-1L, degree))
basis = (inla.fmesher.smorg(loc,
mesh$graph$tv,
bspline = c(n, degree, knots),
fem=-1)$bspline)
}
} else if (identical(mesh$manifold, "S2")) {
loc = mesh$loc
knots = identical(knot.placement, "uniform.latitude")
degree = max(0L, min(n-1L, degree))
basis = (inla.fmesher.smorg(loc,
mesh$graph$tv,
bspline = c(n, degree, knots),
fem=-1)$bspline)
if (!rot.inv) {
warning("Currently only 'rot.inv=TRUE' is supported for B-splines.")
}
} else {
stop("Only know how to make B-splines on R2 and S2.")
}
} else if (identical(type, "sph.harm")) {
if (!identical(mesh$manifold, "S2")) {
stop("Only know how to make spherical harmonics on S2.")
}
if (rot.inv) {
basis = (inla.fmesher.smorg(mesh$loc,
mesh$graph$tv,
sph0=n)$sph0)
} else {
basis = (inla.fmesher.smorg(mesh$loc,
mesh$graph$tv,
sph=n)$sph)
}
}
return(basis)
}
inla.parse.queries <-function(...)
{
queries = list(...)
if (length(queries)==0)
return(queries)
## Make sure that we have a list of names, empty or not:
if (is.null(names(queries))) {
q.names = rep("", length(queries))
} else {
q.names = names(queries)
}
## All nameless entries must be strings with query names. Replace
## empty names with those names, and set those entries to NULL.
for (query.idx in 1:length(queries)) {
if (q.names[[query.idx]]=="") {
if (is.character(queries[[query.idx]])) {
names(queries)[query.idx] = queries[[query.idx]]
queries[query.idx] = list(NULL)
} else {
queries[query.idx] = list(NULL)
warning(paste("Unnamed query ignored. Check query #",
query.idx, sep=""))
}
}
}
return(queries)
}
`inla.fmesher.smorg` <- function(loc, tv,
fem=NULL,
aniso=NULL,
gradients=FALSE,
sph0=NULL,
sph=NULL,
bspline=NULL,
points2mesh=NULL,
splitlines=NULL,
output=NULL,
keep=FALSE)
{
prefix = inla.fmesher.make.prefix(NULL, NULL)
n = nrow(loc)
s.dim = ncol(loc)
if (s.dim==1)
stop("1-D models not implemented yet.")
stopifnot(s.dim>=2, s.dim<=3)
if (missing(output) || is.null(output)) {
output.given = FALSE
output = NULL
} else {
output.given = TRUE
}
output.fem = list("c0", "g1", "g2")
output.aniso = list("g1aniso", "g2aniso")
output.gradients = list("dx", "dy", "dz")
output.sph0 = list("sph0")
output.sph = list("sph")
output.bspline = list("bspline")
output.p2m = list("p2m.t", "p2m.b")
output.splitlines <- list("split.loc", "split.idx", "split.origin",
"split.t", "split.b1", "split.b2")
## Outputs that need +1L index adjustment
indexoutput <- list("split.idx", "split.t", "split.origin")
fmesher.write(inla.affirm.double(loc), prefix, "s")
fmesher.write(inla.affirm.integer(tv)-1L, prefix, "tv")
all.args = "--smorg --input=s,tv"
## additional arguments
if (!is.null(fem)) {
all.args = paste(all.args," --fem=", max(2,fem), sep="")
if (!output.given) output = c(output, output.fem)
}
if (!is.null(aniso)) {
fmesher.write(inla.affirm.double(aniso[[1]]), prefix, "aniso.gamma")
fmesher.write(inla.affirm.double(aniso[[2]]), prefix, "aniso.vec")
all.args = paste(all.args," --aniso=aniso.gamma,aniso.vec", sep="")
if (!output.given) output = c(output, output.aniso)
}
if (!is.null(gradients) && gradients) {
all.args = paste(all.args," --grad", sep="")
if (!output.given) output = c(output, output.gradients)
}
if (!is.null(sph0)) {
all.args = paste(all.args," --sph0=", sph0, sep="")
if (!output.given) output = c(output, output.sph0)
}
if (!is.null(sph)) {
all.args = paste(all.args," --sph=", sph, sep="")
if (!output.given) output = c(output, output.sph)
}
if (!is.null(bspline)) {
all.args = (paste(all.args, " --bspline=",
bspline[1], ",", bspline[2], ",", bspline[3],
sep=""))
if (!output.given) output = c(output, output.bspline)
}
if (!is.null(points2mesh)) {
fmesher.write(inla.affirm.double(points2mesh), prefix, "p2m")
all.args = paste(all.args," --points2mesh=p2m", sep="")
if (!output.given) output = c(output, output.p2m)
}
if (!is.null(splitlines)) {
fmesher.write(inla.affirm.double(splitlines$loc), prefix, "splitloc")
fmesher.write(inla.affirm.integer(splitlines$idx)-1L, prefix, "splitidx")
all.args = paste(all.args," --splitlines=splitloc,splitidx", sep="")
if (!output.given) output = c(output, output.splitlines)
}
all.args = paste(all.args, inla.getOption("fmesher.arg"))
echoc = inla.fmesher.call(all.args=all.args, prefix=prefix)
result = list()
for (name in output) {
if (identical(name, "p2m.t"))
if (!file.exists(paste(prefix, name, sep="")))
result[[name]] = fmesher.read(prefix, "points2mesh.t")+1L
else
result[[name]] = fmesher.read(prefix, name) + 1L
else if (identical(name, "p2m.b"))
if (!file.exists(paste(prefix, name, sep="")))
result[[name]] = fmesher.read(prefix, "points2mesh.b")
else
result[[name]] = fmesher.read(prefix, name)
else if (name %in% indexoutput)
result[[name]] = fmesher.read(prefix, name) + 1L
else
result[[name]] = fmesher.read(prefix, name)
}
if (!keep)
unlink(paste(prefix, "*", sep=""), recursive=FALSE)
return (result)
}
## Deprecated: cyclic
inla.mesh.1d <- function(loc,
interval=range(loc),
boundary=NULL,
degree=1,
free.clamped=FALSE,
...)
{
## Note: do not change the order of these options without also
## changing 'basis.reduction' below.
boundary.options = c("neumann", "dirichlet", "free", "cyclic")
if ("cyclic" %in% names(list(...))) {
warning(paste("Option 'cyclic' is deprecated.",
" Use 'boundary=\"cyclic\"' instead.", sep=""))
if (cyclic) {
if (!missing(boundary)) {
boundary = match.arg.vector(boundary, boundary.options, length=2)
warning(paste("'cyclic=TRUE' overrides 'boundary=c(\"",
paste(boundary, collapse="\",\""), "\")'.", sep=""))
}
boundary = "cyclic"
}
}
boundary = match.arg.vector(boundary, boundary.options, length=2)
cyclic = !is.na(pmatch(boundary[1], "cyclic"))
if (cyclic && is.na(pmatch(boundary[2], "cyclic"))) {
stop("Inconsistent boundary specification 'boundary=c(",
paste(boundary, collapse=","), ")'.", sep="")
}
loc.orig = loc
if (cyclic)
loc =
(sort(unique(c(0, loc-interval[1]) %% diff(interval))) +
interval[1])
else
loc =
(sort(unique(c(interval,
pmax(interval[1], pmin(interval[2], loc)) ))))
n = length(loc)
if (loc[1]<interval[1])
stop("All 'loc' must be >= interval[1].")
if (loc[n]>interval[2])
stop("All 'loc' must be <= interval[2].")
if ( (degree<0) || (degree>2)) {
stop(paste("'degree' must be 0, 1, or 2. 'degree=",
degree,
"' is not supported.", sep=""))
}
if (length(free.clamped) == 1L) {
free.clamped = rep(free.clamped, 2)
}
## Number of basis functions
if (degree==0) {
basis.reduction = c(0, 1, 0, 1/2) ## neu, dir, free, cyclic
} else if (degree==1) {
basis.reduction = c(0, 1, 0, 1/2) ## neu, dir, free, cyclic
} else {
basis.reduction = c(1, 1, 0, 1)
}
m = (n+cyclic+(degree==2)*1
-basis.reduction[pmatch(boundary[1], boundary.options)]
-basis.reduction[pmatch(boundary[2], boundary.options)])
## if (m < 1+max(1,degree)) {
if (m < 1L) {
stop("Degree ", degree,
" meshes must have at least ", 1L,
" basis functions, not 'm=", m, "'.", sep="")
}
## Compute representative basis midpoints.
if ((degree==0) || (degree==1)) {
mid = loc
if (boundary[1] == "dirichlet") {
mid = mid[-1]
}
if (boundary[2] == "dirichlet") {
mid = mid[-(m+1)]
}
} else { ## degree==2
if (cyclic) {
mid = (loc + c(loc[-1], interval[2]))/2
} else {
mid = c(loc[1], (loc[-n]+loc[-1])/2, loc[n])
mid =
switch(boundary[1],
neumann = mid[-1],
dirichlet = mid[-1],
free = mid)
mid =
switch(boundary[2],
neumann = mid[-(m+1)],
dirichlet = mid[-(m+1)],
free = mid)
}
}
mesh =
list(n=n,
m=m,
loc=loc,
mid=mid,
interval=interval,
boundary=boundary,
cyclic=cyclic,
manifold = inla.ifelse(cyclic, "S1", "R1"),
degree=degree,
free.clamped=free.clamped,
idx = list(loc=NULL))
class(mesh) = "inla.mesh.1d"
if (degree<2) {
mesh$idx$loc =
inla.mesh.1d.bary(mesh, loc.orig, method="nearest")$index[,1]
} else {
if (length(mid) >= 2) {
mesh$idx$loc =
inla.mesh.1d.bary(inla.mesh.1d(mid, degree=0),
loc.orig,
method="nearest")$index[,1]
} else {
mesh$idx$loc = rep(1, length(loc.orig))
}
}
return(invisible(mesh))
}
inla.mesh.1d.bary <- function(mesh, loc, method=c("linear", "nearest"))
{
inla.require.inherits(mesh, "inla.mesh.1d", "'mesh'")
method = match.arg(method)
if (method=="linear") {
if (mesh$cyclic) {
mloc = c(mesh$loc-mesh$loc[1], diff(mesh$interval))
loc = (loc-mesh$loc[1]) %% diff(mesh$interval)
} else {
mloc = c(mesh$loc-mesh$loc[1], diff(mesh$interval))
loc = pmax(0, pmin(diff(mesh$interval), loc-mesh$loc[1]))
}
} else {
if (mesh$cyclic) {
mloc =
c(mesh$loc[mesh$n]-diff(mesh$interval),
mesh$loc,
diff(mesh$interval))
mloc = (mloc[-(mesh$n+2)] + mloc[-1])/2
loc = (loc-mloc[1]) %% diff(mesh$interval)
mloc = mloc-mloc[1]
} else {
mloc =
c(0,
(mesh$loc[1:(mesh$n-1L)] +
mesh$loc[2:mesh$n])/2-mesh$loc[1],
diff(mesh$interval))
loc = pmax(0, pmin(diff(mesh$interval), loc-mesh$loc[1]))
}
}
## Binary split method:
do.the.split <- function(knots, loc) {
n = length(knots)
if (n <= 2L) {
return(rep(1L, length(loc)))
}
split = 1L + (n-1L) %/% 2L ## Split point
upper = (loc >= knots[split])
idx = rep(0, length(loc))
idx[!upper] = do.the.split(knots[1:split], loc[!upper])
idx[upper] = split - 1L + do.the.split(knots[split:n], loc[upper])
return(idx)
}
idx = do.the.split(mloc, loc)
if (method=="nearest") {
u = rep(0, length(loc))
if (mesh$cyclic) {
found = which(idx==(mesh$n+1L))
idx[found] = 1L
}
} else { ## (method=="linear") {
u = pmax(0, pmin(1, (loc-mloc[idx])/(mloc[idx+1L]-mloc[idx]) ))
if (!mesh$cyclic) {
found = which(idx==mesh$n)
idx[found] = mesh$n-1L
u[found] = 1
}
}
if (mesh$cyclic) {
index = matrix(c(idx, (idx %% mesh$n)+1L), length(idx), 2)
bary = matrix(c(1-u, u), length(idx), 2)
} else {
index = matrix(c(idx, idx+1L), length(idx), 2)
bary = matrix(c(1-u, u), length(idx), 2)
}
return(list(index=index, bary=bary))
}
inla.mesh.1d.A <- function(mesh, loc,
weights=NULL,
derivatives=NULL,
method=c("linear", "nearest", "quadratic")
)
{
inla.require.inherits(mesh, "inla.mesh.1d", "'mesh'")
if (missing(method)) {
## Compute basis based on mesh$degree and mesh$boundary
if (mesh$degree==0) {
info = (inla.mesh.1d.A(mesh, loc, method="nearest",
weights=weights,
derivatives=
inla.ifelse(is.null(derivatives),
FALSE, TRUE)))
if (mesh$boundary[1] == "dirichlet") {
info$A = info$A[,-1,drop=FALSE]
if (!is.null(derivatives) && derivatives)
info$dA = info$dA[,-1,drop=FALSE]
}
if (mesh$boundary[2] == "dirichlet") {
info$A = info$A[,-(mesh$m+1),drop=FALSE]
if (!is.null(derivatives) && derivatives)
info$dA = info$dA[,-(mesh$m+1),drop=FALSE]
}
if (!is.null(derivatives) && derivatives)
return(info)
else
return(info$A)
} else if (mesh$degree==1) {
info = (inla.mesh.1d.A(mesh, loc, method="linear",
weights=weights,
derivatives=
inla.ifelse(is.null(derivatives),
FALSE, TRUE)))
if (mesh$boundary[1] == "dirichlet") {
info$A = info$A[,-1,drop=FALSE]
if (!is.null(derivatives) && derivatives)
info$dA = info$dA[,-1,drop=FALSE]
}
if (mesh$boundary[2] == "dirichlet") {
info$A = info$A[,-(mesh$m+1),drop=FALSE]
if (!is.null(derivatives) && derivatives)
info$dA = info$dA[,-(mesh$m+1),drop=FALSE]
}
if (!is.null(derivatives) && derivatives)
return(info)
else
return(info$A)
} else if (mesh$degree==2) {
info =
inla.mesh.1d.A(mesh, loc,
method="quadratic",
weights=weights,
derivatives=
inla.ifelse(
is.null(derivatives),
FALSE, TRUE))
adjust.matrix <- function(mesh, A) {
A = inla.as.dgTMatrix(A)
i = A@i+1L
j = A@j+1L
x = A@x
if (mesh$cyclic) {
j[j==1L] = mesh$m+1L
j[j==mesh$m+2L] = 2L
j = j-1L
} else {
if (mesh$boundary[1] == "neumann") {
j[j==1L] = 2L
j = j-1L
} else if (mesh$boundary[1] == "dirichlet") {
x[j==1L] = -x[j==1L]
j[j==1L] = 2L
j = j-1L
} else if ((mesh$boundary[1] == "free") &&
mesh$free.clamped[1]) {
j1 = which(j==1L)
i = c(i, i[j1])
j = c(j, j[j1]+1L)
x = c(x, -x[j1])
x[j1] = 2*x[j1]
}
if (mesh$boundary[2] == "neumann") {
j[j==(mesh$m+1L)] = mesh$m
} else if (mesh$boundary[2] == "dirichlet") {
x[j==(mesh$m+1L)] = -x[j==(mesh$m+1L)]
j[j==(mesh$m+1L)] = mesh$m
} else if ((mesh$boundary[2] == "free") &&
mesh$free.clamped[2]) {
j1 = which(j==(mesh$m))
i = c(i, i[j1])
j = c(j, j[j1]-1L)
x = c(x, -x[j1])
x[j1] = 2*x[j1]
}
}
return(sparseMatrix(i=i, j=j, x=x,
dims=c(length(loc), mesh$m)))
}
if (!is.null(derivatives) && derivatives) {
return(list(A=adjust.matrix(mesh, info$A),
dA=adjust.matrix(mesh, info$dA),
d2A=adjust.matrix(mesh, info$d2A)))
} else {
return(adjust.matrix(mesh, info$A))
}
} else {
stop(paste("'degree' must be 0, 1, or 2. 'degree=",
mesh$degree,
"' is not supported.", sep=""))
}
} else {
method = match.arg(method)
if (missing(weights) || is.null(weights)) {
weights = rep(1, length(loc))
}
if (!is.na(pmatch(method, c("linear", "nearest")))) {
idx = inla.mesh.1d.bary(mesh, loc, method)
dA = NULL
if (method=="linear") {
## Compute the n 1st order B-splines.
i = rep(1:length(loc), times=2)
weights.i = weights[i]
A = (sparseMatrix(i=i,
j=as.vector(idx$index),
x=weights.i*as.vector(idx$bary),
dims=c(length(loc), mesh$n)))
if (!is.null(derivatives) && derivatives) {
if (mesh$cyclic) {
d = (c(mesh$loc[-1], mesh$interval[2]) -
mesh$loc)
} else {
d = (mesh$loc[-1] - mesh$loc[-mesh$n])
}
dA = (sparseMatrix(i=i,
j=as.vector(idx$index),
x=(weights.i*c(-1/d[idx$index[,1]],
1/d[idx$index[,1]])),
dims=c(length(loc), mesh$n)))
}
} else {
## Nearest neighbours.
A = (sparseMatrix(i=1:length(loc),
j=idx$index[,1],
x=weights*idx$bary[,1],
dims=c(length(loc), mesh$n)))
}
if (missing(derivatives) || is.null(derivatives)) {
return(A)
} else {
return(list(A=A, dA=dA))
}
} else { ## Quadratic
## Compute the n+1 2nd order B-splines.
if (mesh$cyclic) {
Boundary.knots =
(c(mesh$loc, mesh$interval[2])[c(mesh$n,2)] +
diff(mesh$interval)*c(-1,1))
knots = c(mesh$loc, mesh$interval[2])
} else {
Boundary.knots =
c(mesh$loc[1] - (mesh$loc[2]-mesh$loc[1]),
mesh$loc[mesh$n] + (mesh$loc[mesh$n]-mesh$loc[mesh$n-1]))
knots = mesh$loc
}
if (FALSE) { ## Only for debugging.
## Using bs():
## Note: Intermediate step constructs dense matrix.
bsobj =
bs(x=loc,
knots=knots,
degree=2,
Boundary.knots=Boundary.knots)
bsobj =
Matrix(as.vector(bsobj[,1:(mesh$n+mesh$cyclic+1)]),
nrow(bsobj), mesh$n+mesh$cyclic+1)
}
## Direct calculation:
knots = c(Boundary.knots[1], knots)
if (mesh$cyclic) {
idx =
inla.mesh.1d.bary(inla.mesh.1d(knots, boundary="free"),
(loc-mesh$interval[1]) %%
diff(mesh$interval) + mesh$interval[1],
method="linear")
} else {
idx =
inla.mesh.1d.bary(inla.mesh.1d(knots, boundary="free"),
loc, method="linear")
}
knots = c(knots, Boundary.knots[2])
idx$index = idx$index[,1]-1L ## Indices into mesh intervals.
d = knots[2:length(knots)]-knots[1:(length(knots)-1)]
d2 = knots[3:length(knots)]-knots[1:(length(knots)-2)]
## Left intervals for each basis function:
i.l = 1:length(idx$index)
j.l = idx$index+2L
x.l = (idx$bary[,2]*d[idx$index+1]/d2[idx$index+1] * idx$bary[,2])
## Right intervals for each basis function:
i.r = 1:length(idx$index)
j.r = idx$index
x.r = (idx$bary[,1]*d[idx$index+1]/d2[idx$index] * idx$bary[,1])
## Middle intervals for each basis function:
i.m = 1:length(idx$index)
j.m = idx$index+1L
x.m = (1-(idx$bary[,1]*d[idx$index+1]/d2[idx$index]* idx$bary[,1] +
idx$bary[,2]*d[idx$index+1]/d2[idx$index+1] * idx$bary[,2]
))
i = c(i.l, i.r, i.m)
j = c(j.l, j.r, j.m)
weights.i = weights[i]
A = (sparseMatrix(i=i, j=j,
x=weights.i*c(x.l, x.r, x.m),
dims=c(length(idx$index), mesh$n+mesh$cyclic+1L)
))
dA = NULL
d2A = NULL
if (!is.null(derivatives) && derivatives) {
## dA:
## Left, right, middle intervals for each basis function:
x.l = (2 / d2[idx$index+1] * idx$bary[,2])
x.r = (-2 / d2[idx$index] * idx$bary[,1])
x.m = (-(-2/ d2[idx$index] * idx$bary[,1] +
2/ d2[idx$index+1] * idx$bary[,2]
))
dA = (sparseMatrix(i=i, j=j,
x=weights.i*c(x.l, x.r, x.m),
dims=(c(length(idx$index),
mesh$n+mesh$cyclic+1L))
))
## d2A:
## Left, right, middle intervals for each basis function:
x.l = (2 / d[idx$index+1] / d2[idx$index+1])
x.r = (2 / d[idx$index+1] / d2[idx$index])
x.m = (-(2/d[idx$index+1] / d2[idx$index] +
2/d[idx$index+1] / d2[idx$index+1] ))
d2A = (sparseMatrix(i=i, j=j,
x=weights.i*c(x.l, x.r, x.m),
dims=(c(length(idx$index),
mesh$n+mesh$cyclic+1L))
))
}
if (missing(derivatives) || is.null(derivatives)) {
return(A)
} else {
return(list(A=A, dA=dA, d2A=d2A))
}
}
}
}
inla.mesh.1d.fem <- function(mesh)
{
inla.require.inherits(mesh, "inla.mesh.1d", "'mesh'")
## Use the same matrices for degree 0 as for degree 1
if ((mesh$degree==0) || (mesh$degree==1)) {
if (mesh$cyclic) {
loc =
c(mesh$loc[mesh$n]-diff(mesh$interval),
mesh$loc,
mesh$loc[1]+diff(mesh$interval))
c0 = (loc[3:length(loc)]-loc[1:(length(loc)-2)])/2
c1.l = (loc[2:(length(loc)-1)]-loc[1:(length(loc)-2)])/6
c1.r = (loc[3:length(loc)]-loc[2:(length(loc)-1)])/6
c1.0 = (c1.l+c1.r)*2
g1.l = -1/(loc[2:(length(loc)-1)]-loc[1:(length(loc)-2)])
g1.r = -1/(loc[3:length(loc)]-loc[2:(length(loc)-1)])
g1.0 = -g1.l-g1.r
i.l = 1:mesh$n
i.r = 1:mesh$n
i.0 = 1:mesh$n
if (mesh$n > 1) {
j.l = c(mesh$n, 1:(mesh$n-1))
j.r = c(2:mesh$n, 1)
j.0 = 1:mesh$n
} else {
j.l = 1L
j.r = 1L
j.0 = 1L
}
} else {
c0 =
c((mesh$loc[2]-mesh$loc[1])/2,
(mesh$loc[mesh$n]-mesh$loc[mesh$n-1])/2
)
if (mesh$n>2) {
c0 =
c(c0[1],
(mesh$loc[3:mesh$n]-mesh$loc[1:(mesh$n-2)])/2,
c0[2]
)
}
c1.l = (mesh$loc[2:mesh$n]-mesh$loc[1:(mesh$n-1)])/6
c1.r = c1.l
c1.0 = (c(0, c1.l)+c(c1.r,0))*2
g1.l = -1/(mesh$loc[2:mesh$n]-mesh$loc[1:(mesh$n-1)])
g1.r = g1.l
g1.0 = -c(0, g1.l)-c(g1.r,0)
i.l = 2:mesh$n
i.r = 1:(mesh$n-1)
i.0 = 1:mesh$n
j.l = 1:(mesh$n-1)
j.r = 2:mesh$n
j.0 = 1:mesh$n
if (mesh$boundary[1] == "dirichlet") {
g1.0 = g1.0[-1]; g1.l = g1.l[-1]; g1.r = g1.r[-1]
c1.0 = c1.0[-1]; c1.l = c1.l[-1]; c1.r = c1.r[-1]
c0 = c0[-1]
i.l = i.l[-1]-1; i.r = i.r[-1]-1; i.0 = i.0[-1]-1
j.l = j.l[-1]-1; j.r = j.r[-1]-1; j.0 = j.0[-1]-1
} else if (mesh$boundary[1] == "free") {
g1.0[1] = 0
g1.r[1] = 0
}
if (mesh$boundary[2] == "dirichlet") {
m = mesh$m
g1.0 = g1.0[-(m+1)]; g1.l = g1.l[-m]; g1.r = g1.r[-m]
c1.0 = c1.0[-(m+1)]; c1.l = c1.l[-m]; c1.r = c1.r[-m]
c0 = c0[-(m+1)]
i.l = i.l[-m]; i.r = i.r[-m]; i.0 = i.0[-(m+1)]
j.l = j.l[-m]; j.r = j.r[-m]; j.0 = j.0[-(m+1)]
} else if (mesh$boundary[2] == "free") {
g1.0[mesh$m] = 0
g1.l[mesh$m-1] = 0
}
}
g1 =
sparseMatrix(i=c(i.l, i.r, i.0),
j=c(j.l, j.r, j.0),
x=c(g1.l, g1.r, g1.0),
dims=c(mesh$m, mesh$m))
c1 =
sparseMatrix(i=c(i.l, i.r, i.0),
j=c(j.l, j.r, j.0),
x=c(c1.l, c1.r, c1.0),
dims=c(mesh$m, mesh$m))
g2 = t(g1) %*% Diagonal(mesh$m, 1/c0) %*% g1
c0 = Diagonal(mesh$m, c0)
} else if (mesh$degree==2) {
if (mesh$cyclic) {
knots1 = mesh$loc
knots2 = c(mesh$loc[-1], mesh$interval[2])
} else {
knots1 = mesh$loc[-mesh$n]
knots2 = mesh$loc[-1]
}
knots.m = (knots1+knots2)/2
knots.d = (knots2-knots1)/2
## 3-point Gaussian quadrature
info =
inla.mesh.1d.A(mesh,
loc=(c(knots.m,
knots.m - knots.d*sqrt(3/5),
knots.m + knots.d*sqrt(3/5))),
weights =
c(knots.d*8/9, knots.d*5/9, knots.d*5/9)^0.5,
derivatives = TRUE
)
c1 = t(info$A) %*% info$A
g1 = t(info$dA) %*% info$dA
g2 = t(info$d2A) %*% info$d2A
g01 = t(info$A) %*% info$dA
g02 = t(info$A) %*% info$d2A
g12 = t(info$dA) %*% info$d2A
c0 = Diagonal(nrow(c1), rowSums(c1))
return(list(c0=c0, c1=c1, g1=g1, g2=g2, g01=g01, g02=g02, g12=g12))
} else {
stop(paste("Mesh basis degree=", mesh$degree,
" is not supported.", sep=""))
}
return(list(c0=c0, c1=c1, g1=g1, g2=g2))
}
inla.diameter <- function(x, ...) {
UseMethod("inla.diameter")
}
## Calculate upper bound for the diameter of a point set,
## by encapsulating in a circular domain.
inla.diameter.default <- function(x, manifold="", ...) {
if (nrow(x) <= 1) {
0
} else {
if (identical(manifold, "S2")) {
distance <- function(u,v) {
2*asin(pmin(1,
((u[1]-v[,1])^2 + (u[2]-v[,2])^2 + (u[3]-v[,3])^2)^0.5 / 2))
}
center <- colMeans(x)
tmp <- sqrt(sum(center^2))
if (tmp < 1e-6) {
pi
} else {
center <- center/tmp
min(pi, 2 * max(distance(center, x)))
}
} else {
distance <- function(u,v) {
d <- 0
for (k in seq_len(ncol(v))) {
d <- d + (u[k]-v[,k])^2
}
d^0.5
}
center <- rep(0, ncol(x))
for (k in seq_len(ncol(x))) {
center[k] <- mean(range(x[,k]))
}
2 * max(distance(center, x))
}
}
}
inla.diameter.inla.mesh.1d <- function(x, ...) {
diff(x$interval)
}
inla.diameter.inla.mesh <- function(x, ...) {
inla.diameter(x$loc, manifold=x$manifold, ...)
}
inla.diameter.inla.mesh.segment <- function(x, ...) {
inla.diameter(x$loc, ...)
}
inla.diameter.inla.mesh.lattice <- function(x, ...) {
inla.diameter(x$loc, ...)
}
inla.mesh.fem <- function(mesh, order=2)
{
inla.require.inherits(mesh, c("inla.mesh", "inla.mesh.1d"), "'mesh'")
if (inherits(mesh, "inla.mesh.1d")) {
if (order > 2) {
stop("order>2 not supported for 1d meshes.")
## TODO: Compute higher order matrices based on order 2 output.
}
return(inla.mesh.1d.fem(mesh))
} else {
## output name list:
if (order > 0) {
output <- paste("g", seq_len(order), sep="")
} else {
output <- NULL
}
output = c("c0", "c1", output,
"va", "ta")
return(inla.fmesher.smorg(mesh$loc, mesh$graph$tv,
fem=order, output=output))
}
}
inla.mesh.deriv <- function(mesh, loc)
{
inla.require.inherits(mesh, c("inla.mesh", "inla.mesh.1d"), "'mesh'")
if (inherits(mesh, "inla.mesh.1d")) {
stop("1d derivatives are not implemented.")
}
n.mesh = mesh$n
info = inla.mesh.project(mesh, loc=loc)
ii = which(info$ok)
n.ok = sum(info$ok)
tv = mesh$graph$tv[info$t[ii,1L],]
e1 = mesh$loc[tv[,3],]-mesh$loc[tv[,2],]
e2 = mesh$loc[tv[,1],]-mesh$loc[tv[,3],]
e3 = mesh$loc[tv[,2],]-mesh$loc[tv[,1],]
n1 = e2-e1*matrix(rowSums(e1*e2)/rowSums(e1*e1),n.ok,3)
n2 = e3-e2*matrix(rowSums(e2*e3)/rowSums(e2*e2),n.ok,3)
n3 = e1-e3*matrix(rowSums(e3*e1)/rowSums(e3*e3),n.ok,3)
g1 = n1/matrix(rowSums(n1*n1),n.ok,3)
g2 = n2/matrix(rowSums(n2*n2),n.ok,3)
g3 = n3/matrix(rowSums(n3*n3),n.ok,3)
x = cbind(g1[,1],g2[,1],g3[,1])
y = cbind(g1[,2],g2[,2],g3[,2])
z = cbind(g1[,3],g2[,3],g3[,3])
dx = (sparseMatrix(dims=c(nrow(loc),n.mesh),
i = rep(ii, 3),
j = as.vector(tv),
x = as.vector(x) ))
dy = (sparseMatrix(dims=c(nrow(loc),n.mesh),
i = rep(ii, 3),
j = as.vector(tv),
x = as.vector(y) ))
dz = (sparseMatrix(dims=c(nrow(loc),n.mesh),
i = rep(ii, 3),
j = as.vector(tv),
x = as.vector(z) ))
return (list(A=info$A, dx=dx, dy=dy, dz=dz))
}
inla.simplify.curve <- function(loc, idx, eps) {
## Variation of Ramer-Douglas-Peucker
## Uses width epsilon ellipse instead of rectangle,
## motivated by prediction ellipse for Brownian bridge
n = length(idx)
if ((n==2) || (eps==0)) {
return(idx)
}
segm = loc[idx[n],]-loc[idx[1],]
segm.len = sum(segm^2)^0.5
if (segm.len <= 1e-12) {
## End point same as start; closed curve. Split.
len2 = ((loc[idx[2:(n-1)], 1] - loc[idx[1], 1])^2 +
(loc[idx[2:(n-1)], 2] - loc[idx[1], 2])^2 )
split = which.max(len2)+1L
} else {
segm.mid = (loc[idx[n],]+loc[idx[1],])/2
segm = segm/segm.len
segm.perp = c(-segm[2], segm[1])
vec = (cbind(loc[idx[2:(n-1)], 1] - segm.mid[1],
loc[idx[2:(n-1)], 2] - segm.mid[2]))
## Always split if any point is outside the circle
epsi = min(c(eps, segm.len/2))
dist1 = abs(vec[,1]*segm[1] + vec[,2]*segm[2])/(segm.len/2)*epsi
dist2 = abs(vec[,1]*segm.perp[1] + vec[,2]*segm.perp[2])
dist = (dist1^2 + dist2^2)^0.5
## Find the furthest point, in the ellipse metric, and
## check if it inside the radius (radius=segm.len/2)
split = which.max(dist)+1L
if (dist[split-1L] < epsi) {
## Flat segment, eliminate.
return(idx[c(1, n)])
}
## Split at the furthest point.
split = which.max(dist)+1L
}
## Do the split recursively:
return(c(inla.simplify.curve(loc, idx[1L:split], eps),
inla.simplify.curve(loc, idx[split:n], eps)[-1L]
))
}
inla.contour.segment <-
function(x = seq(0, 1, length.out = nrow(z)),
y = seq(0, 1, length.out = ncol(z)),
z, nlevels = 10,
levels = pretty(range(z, na.rm=TRUE), nlevels),
groups = seq_len(length(levels)),
positive = TRUE,
eps = NULL,
crs = NULL)
{
## Input checking from contourLines:
if (missing(z)) {
if (!missing(x)) {
if (is.list(x)) {
z <- x$z
y <- x$y
x <- x$x
}
else {
z <- x
x <- seq.int(0, 1, length.out = nrow(z))
}
}
else stop("no 'z' matrix specified")
}
else if (is.list(x)) {
y <- x$y
x <- x$x
}
if (is.null(eps)) {
eps = min(c(min(diff(x)), min(diff(y))))/8
}
## End of input checking.
## Get contour pieces
curves = contourLines(x,y,z,levels=levels)
## Make a mesh for easy gradient interpolation:
latt = inla.mesh.lattice(x,y)
mesh =
inla.mesh.create(lattice=latt,
boundary=latt$segm,
extend=(list(n=3,
offset=(max(diff(range(x)),
diff(range(y)))*0.1)
))
)
## Map function values to mesh indexing:
zz = rep(0, prod(dim(z)))
zz[mesh$idx$lattice] = as.vector(z)
## Mapping from level to group value:
level2grp = function(level) {
if (length(groups)==1)
return(groups)
for (k in seq_along(groups)) {
if (levels[k] == level) {
return(groups[k])
}
}
return(0)
}
## Join all contour pieces into a single mesh.segment
## Different levels can later be identified via the grp indices.
loc = matrix(0,0,2)
idx = matrix(0,0,2)
grp = c()
for (k in seq_len(length(curves))) {
curve.loc = cbind(curves[[k]]$x, curves[[k]]$y)
curve.n = nrow(curve.loc)
## Extract the rotated gradients along the curve
curve.mid = (curve.loc[1:(curve.n-1),]+curve.loc[2:curve.n,])/2
A = inla.mesh.deriv(mesh, loc=curve.mid)
## Gradients rotated 90 degrees CW, i.e. to the direction
## of CCW curves around positive excursions:
grid.diff = cbind(A$dy %*% zz, -A$dx %*% zz)
## Determine the CCW/CW orientation
curve.diff = diff(curve.loc)
ccw = (sum(curve.diff * grid.diff) >= 0) ## True if in CCW direction
if ((ccw && positive) | (!ccw && !positive)) {
curve.idx = 1:curve.n
} else {
curve.idx = curve.n:1
}
## Filter short line segments:
curve.idx =
inla.simplify.curve(curve.loc,
curve.idx,
eps = eps)
## Reorder, making sure any unused points are removed:
curve.loc = curve.loc[curve.idx,,drop=FALSE]
curve.n = nrow(curve.loc)
curve.idx = cbind(1L:(curve.n-1L), 2L:curve.n)
## Check if the curve is closed, and adjust if it is:
if (max(abs(curve.loc[1,] - curve.loc[curve.n,])) < 1e-12) {
curve.loc = curve.loc[-curve.n,,drop=FALSE]
curve.n = nrow(curve.loc)
curve.idx = cbind(1L:curve.n, c(2L:curve.n, 1L))
}
## Add the curve:
offset = nrow(loc)
loc = rbind(loc, curve.loc)
idx = rbind(idx, curve.idx+offset)
grp = c(grp, rep(level2grp(curves[[k]]$level), curve.n-1L))
}
inla.mesh.segment(loc=loc, idx=idx, grp=grp, is.bnd=FALSE, crs=crs)
}
## Based on an idea from Elias Teixeira Krainski
## Requires splancs::nndistF
inla.nonconvex.hull.basic <-
function(points, convex=-0.15, resolution=40, eps=NULL, crs=NULL)
{
if (!(missing(points) || is.null(points)) &&
(inherits(points, "SpatialPoints") ||
inherits(points, "SpatialPointsDataFrame"))) {
points <- inla.spTransform(coordinates(points),
CRS(proj4string(points)),
crs,
passthrough=TRUE)
}
if (length(convex)==1)
convex = rep(convex,2)
if (length(resolution)==1)
resolution = rep(resolution,2)
lim = rbind(range(points[,1]), range(points[,2]))
ex = convex
if (convex[1]<0) {ex[1] = -convex[1]*diff(lim[1,])}
if (convex[2]<0) {ex[2] = -convex[2]*diff(lim[2,])}
domain = c(diff(lim[1,]), diff(lim[2,])) + 2*ex
dif = domain/(resolution-1)
if (any(dif > min(convex))) {
req.res = ceiling(domain/convex+1)
warning(paste("Resolution (",
paste(resolution,collapse=","),
") too small for convex (",
paste(convex,collapse=","),
").\n",
"Resolution >=(",
paste(req.res,collapse=","),
") required for more accurate results.",
sep=""))
}
ax = list(seq(lim[1,1] - ex[1], lim[1,2] + ex[1], length=resolution[1]),
seq(lim[2,1] - ex[2], lim[2,2] + ex[2], length=resolution[2]))
xy = as.matrix(expand.grid(ax[[1]], ax[[2]]))
tr = diag(c(1/ex[1],1/ex[2]))
stopifnot(inla.require("splancs"))
z = (matrix(splancs::nndistF(points%*%tr, xy%*%tr),
resolution[1], resolution[2]))
segm = inla.contour.segment(ax[[1]], ax[[2]], z,
levels=c(1), positive=FALSE, eps=eps, crs=crs)
return(segm)
}
## Morphological dilation by "convex",
## followed by closing by "concave", with
## minimum concave curvature radius "concave".
## If the dilated set has no gaps of width between
## 2*convex*(sqrt(1+2*concave/convex) - 1)
## and
## 2*concave,
## then the minimum convex curvature radius is "convex".
## Default is concave=convex
## Special case concave=0 delegates to inla.nonconvex.hull.basic()
##
## The implementation is based on the identity
## dilation(a) & closing(b) = dilation(a+b) & erosion(b)
## where all operations are with respect to disks with the specified radii.
inla.nonconvex.hull <-
function(points, convex=-0.15, concave=convex, resolution=40, eps=NULL,
crs=NULL)
{
if (!(missing(points) || is.null(points)) &&
(inherits(points, "SpatialPoints") ||
inherits(points, "SpatialPointsDataFrame"))) {
points <- inla.spTransform(coordinates(points),
CRS(proj4string(points)),
crs,
passthrough=TRUE)
}
if (length(resolution)==1)
resolution = rep(resolution,2)
lim = rbind(range(points[,1]), range(points[,2]))
approx.diam = max(diff(lim[1,]),diff(lim[2,]))
if (convex<0) {convex = -convex*approx.diam}
if (concave<0) {concave = -concave*approx.diam}
if (concave==0) {
return(inla.nonconvex.hull.basic(points, convex, resolution, eps,
crs=crs))
}
ex = convex+concave
domain = c(diff(lim[1,]), diff(lim[2,])) + 2*ex
dif = domain/(resolution-1)
if (max(dif) > min(convex,concave)) {
req.res = ceiling(domain/min(convex,concave)+1)
warning(paste("Resolution (",
paste(resolution,collapse=","),
") too small for convex/concave radius (",
convex,",",concave,
").\n",
"Resolution >=(",
paste(req.res,collapse=","),
") required for more accurate results.",
sep=""))
}
ax =
list(
seq(lim[1,1] - ex, lim[1,2] + ex, length=resolution[1]),
seq(lim[2,1] - ex, lim[2,2] + ex, length=resolution[2])
)
xy = as.matrix(expand.grid(ax[[1]], ax[[2]]))
stopifnot(inla.require("splancs"))
z = (matrix(splancs::nndistF(points, xy),
resolution[1],resolution[2]))
segm.dilation =
inla.contour.segment(ax[[1]], ax[[2]], z,
levels=c(convex+concave),
positive=TRUE,
eps=0) ## Don't simplify curve at this stage
mesh.dilation =
inla.mesh.create(loc=xy,
boundary=segm.dilation,
extend=(list(n=3,
offset=(max(diff(ax[[1]]),
diff(ax[[2]]))*0.1)
))
)
## This filtering is not necessary; the inla.mesh.create() should
## have removed all unused points. 2013-02-24 /FL
## points.dilation =
## mesh.dilation$loc[unique(as.vector(mesh.dilation$graph$tv)),]
z = (matrix(splancs::nndistF(mesh.dilation$loc, xy),
resolution[1],resolution[2]))
segm.closing =
inla.contour.segment(ax[[1]],ax[[2]],z,
levels=c(concave),
positive=TRUE,
eps=eps)
segm.closing$crs <- crs
segm.closing
}
inbo/INLA documentation built on Dec. 6, 2019, 9:51 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions inla.mesh.segment inla.mesh.segment.default inla.mesh.segment.inla.mesh.segment lines.inla.mesh.segment `inla.generate.colors` `plot.inla.trimesh` plot.inla.mesh inla.mesh.map.lim inla.mesh.map inla.mesh.lattice extract.groups extract.groups.inla.mesh.segment inla.mesh.parse.segm.input inla.mesh.filter.locations inla.mesh inla.mesh.create inla.mesh.extract.segments inla.mesh.boundary inla.mesh.interior inla.mesh.2d inla.mesh.create.helper inla.delaunay inla.mesh.query summary.inla.mesh print.summary.inla.mesh inla.mesh.project inla.mesh.project.inla.mesh inla.mesh.project.inla.mesh.1d inla.mesh.project.inla.mesh.projector inla.mesh.projector inla.mesh.projector.inla.mesh inla.mesh.projector.inla.mesh.1d inla.internal.make.spline.mesh inla.mesh.basis inla.parse.queries `inla.fmesher.smorg` inla.mesh.1d inla.mesh.1d.bary inla.mesh.1d.A inla.mesh.1d.fem inla.diameter inla.diameter.default inla.diameter.inla.mesh.1d inla.diameter.inla.mesh inla.diameter.inla.mesh.segment inla.diameter.inla.mesh.lattice inla.mesh.fem inla.mesh.deriv inla.simplify.curve inla.contour.segment inla.nonconvex.hull.basic inla.nonconvex.hull
Documented in extract.groups extract.groups.inla.mesh.segment inla.contour.segment inla.delaunay inla.diameter inla.diameter.default inla.diameter.inla.mesh inla.diameter.inla.mesh.1d inla.diameter.inla.mesh.lattice inla.diameter.inla.mesh.segment inla.mesh inla.mesh.1d inla.mesh.1d.A inla.mesh.1d.bary inla.mesh.1d.fem inla.mesh.2d inla.mesh.basis inla.mesh.boundary inla.mesh.create inla.mesh.create.helper inla.mesh.deriv inla.mesh.fem inla.mesh.interior inla.mesh.lattice inla.mesh.map inla.mesh.map.lim inla.mesh.project inla.mesh.project.inla.mesh inla.mesh.project.inla.mesh.1d inla.mesh.project.inla.mesh.projector inla.mesh.projector inla.mesh.projector.inla.mesh inla.mesh.projector.inla.mesh.1d inla.mesh.query inla.mesh.segment inla.mesh.segment.default inla.mesh.segment.inla.mesh.segment inla.nonconvex.hull inla.nonconvex.hull.basic inla.simplify.curve lines.inla.mesh.segment plot.inla.mesh print.summary.inla.mesh summary.inla.mesh
## Export: inla.contour.segment inla.delaunay inla.fmesher.smorg
## Export: inla.generate.colors inla.mesh inla.mesh.1d inla.mesh.1d.A
## Export: inla.mesh.1d.bary inla.mesh.1d.fem inla.mesh.2d inla.mesh.basis
## Export: inla.mesh.boundary inla.mesh.create inla.mesh.create.helper
## Export: inla.mesh.deriv inla.mesh.fem
## Export: inla.mesh.interior inla.mesh.lattice inla.mesh.map
## Export: inla.mesh.map.lim inla.mesh.query
## Export: inla.nonconvex.hull inla.nonconvex.hull.basic
## Export: inla.simplify.curve plot.inla.trimesh plot!inla.trimesh
## Internal: inla.mesh.filter.locations
## Internal: inla.mesh.parse.segm.input inla.mesh.extract.segments
##
## S3methods; also export some methods explicitly
## Export: extract.groups inla.mesh.project inla.mesh.projector
## Export: extract.groups!inla.mesh.segment
## Export: inla.mesh.segment
## Export: inla.mesh.segment.default
## Export: inla.mesh.segment.inla.mesh.segment
## Export: inla.mesh.segment!default
## Export: inla.mesh.segment!inla.mesh.segment
## Export: inla.mesh.project.inla.mesh inla.mesh.project.inla.mesh.1d
## Export: inla.mesh.project.inla.mesh.projector
## Export: inla.mesh.projector.inla.mesh inla.mesh.projector.inla.mesh.1d
## Export: inla.mesh.project!inla.mesh inla.mesh.project!inla.mesh.1d
## Export: inla.mesh.project!inla.mesh.projector
## Export: inla.mesh.projector!inla.mesh inla.mesh.projector!inla.mesh.1d
## Export: lines.inla.mesh.segment plot.inla.mesh
## Export: print.summary.inla.mesh summary.inla.mesh
## Export: lines!inla.mesh.segment plot!inla.mesh
## Export: print!summary.inla.mesh summary!inla.mesh
## Export: inla.diameter
## Export: inla.diameter!default
## Export: inla.diameter!inla.mesh.1d
## Export: inla.diameter!inla.mesh
## Export: inla.diameter!inla.mesh.segment
## Export: inla.diameter!inla.mesh.lattice
inla.mesh.segment <- function(...) {
UseMethod("inla.mesh.segment")
}
inla.mesh.segment.default <-
function(loc = NULL, idx = NULL, grp = NULL, is.bnd = TRUE,
crs=NULL, ...)
{
if ((missing(loc) || is.null(loc)) &&
(missing(idx) || is.null(idx))) {
stop("At most one of 'loc' and 'idx' may be missing or null.")
}
if (!missing(loc) && !is.null(loc)) {
## Handle loc given as SpatialPoints or SpatialPointsDataFrame object
if (inherits(loc, "SpatialPoints") ||
inherits(loc, "SpatialPointsDataFrame")) {
loc <- inla.spTransform(coordinates(loc),
CRS(proj4string(loc)),
crs,
passthrough=TRUE)
}
if (!is.matrix(loc)) {
loc = as.matrix(loc)
}
if (!is.double(loc)) {
storage.mode(loc) = "double"
}
if (missing(idx) || is.null(idx))
idx = (inla.ifelse(is.bnd,
c(1:nrow(loc),1),
c(1:nrow(loc))))
} else {
loc = NULL
}
if (!missing(idx) && !is.null(idx)) {
if (!is.vector(idx) && !is.matrix(idx))
stop("'idx' must be a vector or a matrix")
if (is.vector(idx))
idx = as.matrix(idx, nrow=length(idx), ncol=1)
if (ncol(idx) == 1) {
if (nrow(idx) < 2) {
if (nrow(idx) == 1) {
warning("Segment specification must have at least 2, or 0, indices.")
}
idx <- matrix(0L, 0, 2)
} else {
idx = matrix(c(idx[-nrow(idx)],idx[-1]), nrow(idx)-1, 2)
}
}
storage.mode(idx) <- "integer"
if (!is.null(loc) &&
(nrow(idx) > 0) &&
(max(idx, na.rm=TRUE) > nrow(loc))) {
warning("Segment indices (max=", max(idx, na.rm=TRUE),
") exceed specified location list length (",
nrow(loc), ").")
}
}
if (!missing(grp) && !is.null(grp)) {
if (!is.vector(grp) && !is.matrix(grp))
stop("'grp' must be a vector or a matrix")
grp = matrix(grp, min(length(grp), nrow(idx)), 1)
if (nrow(grp)<nrow(idx))
grp = (matrix(c(as.vector(grp),
rep(grp[nrow(grp)], nrow(idx)-length(grp))),
nrow(idx), 1))
storage.mode(grp) <- "integer"
} else
grp = NULL
## Filter away NAs in loc and idx
if (!is.null(loc)) {
idx[is.na(idx)] = 0L ## Avoid R annoyances with logical+NA indexing
while (sum(is.na(loc))>0) {
i = min(which(rowSums(is.na(loc))>0))
loc = loc[-i,,drop=FALSE]
idx[idx==i] = 0L
idx[idx>i] = idx[idx>i]-1L
}
idx[idx==0L] = NA
}
while (sum(is.na(idx))>0) {
i = min(which(rowSums(is.na(idx))>0))
idx = idx[-i,,drop=FALSE]
if (!is.null(grp))
grp = grp[-i,,drop=FALSE]
}
if (!is.null(loc)) {
## Identify unused locations and remap indices accordingly.
idx.new = rep(0L, nrow(loc))
idx.new[as.vector(idx)] = 1L
loc = loc[idx.new==1L,, drop=FALSE]
idx.new[idx.new==1L] = seq_len(sum(idx.new))
idx = (matrix(idx.new[as.vector(idx)],
nrow=nrow(idx),
ncol=ncol(idx)))
}
ret = list(loc=loc, idx=idx, grp=grp, is.bnd=is.bnd, crs=crs)
class(ret) <- "inla.mesh.segment"
return(invisible(ret))
}
inla.mesh.segment.inla.mesh.segment <- function(..., grp.default=0) {
segm <- list(...)
if (!all(unlist(lapply(segm,
function(x) inherits(x,"inla.mesh.segment"))))) {
stop("All objects must be of class 'inla.mesh.segment'.")
}
Nloc <- unlist(lapply(segm, function(x) nrow(x$loc)))
cumNloc <- c(0, cumsum(Nloc))
Nidx <- unlist(lapply(segm, function(x) nrow(x$idx)))
loc <- do.call(rbind, lapply(segm, function(x) x$loc))
idx <- do.call(rbind, lapply(seq_along(segm),
function(x) segm[[x]]$idx+cumNloc[x]))
grp <- unlist(lapply(seq_along(segm),
function(x) {
if (is.null(segm[[x]]$grp)) {
rep(grp.default, Nidx[x])
} else {
segm[[x]]$grp
}
}))
is.bnd <- unlist(lapply(segm, function(x) x$is.bnd))
if (!all(is.bnd) || (any(!is.bnd) && !all(!is.bnd))) {
warning("Inconsistent 'is.bnd' attributes. Setting 'is.bnd=FALSE'.")
is.bnd <- FALSE
} else {
is.bnd <- all(is.bnd)
}
crs <- lapply(segm, function(x) x$crs)
if (!is.null(crs)) {
crs <- crs[unlist(lapply(crs,
function(x) !is.null(x)))]
if (length(crs) > 0) {
if (!all(unlist(lapply(crs,
function(x) identical(crs[[1]], x))))) {
lapply(crs, function(x) show(x))
stop("Inconsistent 'crs' attributes.")
} else {
crs <- crs[[1]]
}
} else {
crs <- NULL
}
}
invisible(inla.mesh.segment(loc=loc, idx=idx, grp=grp, is.bnd=is.bnd,
crs=crs))
}
lines.inla.mesh.segment <- function(x, loc=NULL, col=NULL,
colors=c("black", "blue", "red", "green"),
add=TRUE, xlim=NULL, ylim=NULL,
rgl=FALSE, ...)
{
segm = x
if (!is.null(segm$loc))
loc = segm$loc
stopifnot(!is.null(loc), ncol(loc)>=2)
if (ncol(loc) < 3) {
loc = cbind(loc, 0.0)
}
color = col
if (rgl) {
stopifnot(inla.require("rgl"))
if (!add) {
dev = rgl::open3d()
rgl::view3d(0, 0, fov=0)
} else {
dev = NULL
}
} else if (!add) {
idx = unique(as.vector(segm$idx))
if (is.null(xlim))
xlim=range(loc[idx,1])
if (is.null(ylim))
ylim=range(loc[idx,2])
plot.new()
plot.window(xlim=xlim, ylim=ylim, ...)
}
grps = inla.ifelse(is.null(segm$grp), rep(0L,nrow(segm$idx)), segm$grp)
for (grp in unique(grps)) {
idx = which(grps==grp)
if (is.null(col)) {
color=colors[1+(grp%%length(colors))]
}
if (rgl) {
rgl::segments3d(loc[as.vector(t(segm$idx[idx,, drop=FALSE])),,drop=FALSE],
color=color,
...)
} else {
lines(loc[t(cbind(segm$idx[idx,, drop=FALSE], NA)), 1],
loc[t(cbind(segm$idx[idx,, drop=FALSE], NA)), 2],
col=color,
...)
}
}
}
`inla.generate.colors` <- function(color,
color.axis = NULL,
color.n=512,
color.palette = cm.colors,
color.truncate=FALSE,
alpha=NULL)
{
if (is.character(color)) {
colors = color
} else if (is.vector(color) || (is.matrix(color) && (ncol(color)==1))) {
if (is.null(color.axis))
color.axis = c(min(color, na.rm=TRUE), max(color, na.rm=TRUE))
if (color.truncate) {
not.ok = ((color<color.axis[1]) |
(color>color.axis[2]))
} else {
not.ok = rep(FALSE, length(color))
}
cs = (pmax(color.axis[1],
pmin(color.axis[2], color, na.rm=TRUE), na.rm=TRUE))
cs = (cs-color.axis[1])/(color.axis[2]-color.axis[1])
not.ok = not.ok | is.na(cs)
cs[not.ok] = 0.5
if (is.null(alpha)) {
alpha = as.numeric(!not.ok)
} else {
alpha[not.ok] = 0
}
ics = (as.numeric(cut(cs, seq(0, 1, length.out=color.n+1),
include.lowest=TRUE)))
colors = color.palette(color.n)[ics]
## Todo: handle alpha, combining "input alpha" with "not.ok-alpha"
} else if (is.matrix(color) && (ncol(color)==3)) {
if (is.null(color.axis))
color.axis = c(min(color, na.rm=TRUE), max(color, na.rm=TRUE))
if (color.truncate) {
not.ok = ((color[, 1]<color.axis[1]) |
(color[, 2]<color.axis[1]) |
(color[, 3]<color.axis[1]) |
(color[, 1]>color.axis[2]) |
(color[, 2]>color.axis[2]) |
(color[, 3]>color.axis[2]))
} else {
not.ok = rep(FALSE, nrow(color))
}
cs = matrix(
pmax(color.axis[1],
pmin(color.axis[2], color, na.rm=TRUE), na.rm=TRUE), dim(color))
cs = (cs-color.axis[1])/(color.axis[2]-color.axis[1])
not.ok = not.ok | is.na(cs[, 1]) | is.na(cs[, 2]) | is.na(cs[, 3])
cs[not.ok,] = c(0.5, 0.5, 0.5)
if (is.null(alpha)) {
alpha = as.numeric(!not.ok)
} else {
alpha[not.ok] = 0
}
colors = rgb(cs[, 1], cs[, 2], cs[, 3])
} else {
stop("color specification must be character, matrix, or vector.")
}
return (list(colors=colors, alpha=alpha))
}
## Export as plot!inla.trimesh as well as plot.inla.trimesh,
## even though it's not actually an S3 method.
`plot.inla.trimesh` <- function(x, S, color = NULL, color.axis = NULL,
color.n=512, color.palette = cm.colors,
color.truncate=FALSE, alpha=NULL,
lwd = 1, specular = "black",
draw.vertices=TRUE, draw.edges=TRUE,
edge.color=rgb(0.3, 0.3, 0.3), ...)
{
TV = x
stopifnot(inla.require("rgl"))
## Make indices 1 based. Deprecated and is deactivated.
if (min(TV) == 0) {
stop("Zero-based indices in TV are not supported.")
}
colors = (inla.generate.colors(color, color.axis, color.n,
color.palette, color.truncate, alpha))
tTV = t(TV);
tETV = t(TV[, c(1, 2, 3, 1, NA)]);
Tx = S[tTV, 1]
Ty = S[tTV, 2]
Tz = S[tTV, 3]
if (length(colors$colors) == 1) {
## One color
Tcol = colors$colors
Talpha = colors$alpha
} else if (length(colors$colors) == nrow(S)) {
## One color per vertex
Tcol = colors$colors[tTV]
Talpha = colors$alpha[tTV]
} else {
## One color per triangle
stopifnot(length(colors$colors) == nrow(TV))
Tcol = colors$colors[t(matrix(rep(1:nrow(TV), 3), dim(TV)))]
Talpha = colors$alpha[t(matrix(rep(1:nrow(TV), 3), dim(TV)))]
}
Ex = S[tETV, 1]
Ey = S[tETV, 2]
Ez = S[tETV, 3]
Ecol = edge.color
if (draw.vertices) {
rgl::points3d(S, color="black", ...)
}
if (draw.edges) {
rgl::lines3d(Ex, Ey, Ez, color=Ecol, lwd=lwd, ...)
}
rgl::triangles3d(Tx, Ty, Tz, color=Tcol, specular=specular, alpha=Talpha, ...)
return (invisible())
}
## library(geometry)
## S = cbind(x=rnorm(30), y=rnorm(30), z=0)
## TV = delaunayn(S[, 1:2]) # NOTE: inconsistent triangle orders, only for test.
## trimesh(TV, S)
##
## colors = rgb(runif(30), runif(30), runif(30))
## rgl.viewpoint(0, 0, fov=20)
## plot.inla.trimesh(TV, S, colors)
## Ecol = col2rgb(color)/256
## Ecol = Ecol*0.5+(1-0.5)*0 # Rescale towards black
## Ecol = 1-Ecol # Invert
## Ecol = Ecol[, c(2, 3, 1)] # Permute
## Ecol = rgb(Ecol[1,], Ecol[2,], Ecol[3,], maxColorValue = 1)
## Ecol = Ecol[tETV]
plot.inla.mesh <- function(x,
col="white",
t.sub=1:nrow(mesh$graph$tv),
add=FALSE,
lwd=1,
xlim = range(mesh$loc[,1]),
ylim = range(mesh$loc[,2]),
main = NULL,
rgl = FALSE,
size = 2,
draw.vertices=FALSE,
vertex.color="black",
draw.edges=TRUE,
edge.color=rgb(0.3, 0.3, 0.3),
draw.segments=draw.edges,
...)
{
inla.require.inherits(x, "inla.mesh", "'mesh'")
mesh = x
if (rgl) {
stopifnot(inla.require("rgl"))
if (!add) {
dev = rgl::open3d()
rgl::view3d(0, 0, fov=0)
} else {
dev = NULL
}
tv = mesh$graph$tv[t.sub,,drop=FALSE]
if (draw.vertices) {
idx = intersect(unique(as.vector(tv)), mesh$idx$loc)
rgl::points3d(mesh$loc[idx,,drop=FALSE],
size=2*size, lwd=lwd, color = "blue", ...)
}
if (draw.segments) {
if (!is.null(mesh$segm$bnd))
lines(mesh$segm$bnd, mesh$loc, lwd=lwd+1,
rgl=TRUE, add=TRUE, ...)
if (!is.null(mesh$segm$int))
lines(mesh$segm$int, mesh$loc, lwd=lwd+1,
rgl=TRUE, add=TRUE, ...)
}
plot.inla.trimesh(tv, mesh$loc, color = col,
size=size, lwd=lwd,
draw.vertices=draw.vertices,
vertex.color=vertex.color,
draw.edges=draw.edges,
edge.color=edge.color,
...)
return(invisible(dev))
} else {
idx = cbind(mesh$graph$tv[t.sub,c(1:3,1), drop=FALSE], NA)
x = mesh$loc[t(idx), 1]
y = mesh$loc[t(idx), 2]
if (!add) {
plot.new()
plot.window(xlim=xlim, ylim=ylim, ...)
}
if (draw.edges) {
lines(x, y, type="l", col=edge.color, lwd=lwd)
}
tv = mesh$graph$tv[t.sub,,drop=FALSE]
if (draw.vertices) {
idx = unique(as.vector(tv))
points(mesh$loc[idx,,drop=FALSE],
pch=20, col=vertex.color, ...)
idx = intersect(idx, mesh$idx$loc)
points(mesh$loc[idx,,drop=FALSE],
pch=20, col="blue", ...)
}
if (draw.segments) {
if (!is.null(mesh$segm$bnd))
lines(mesh$segm$bnd, mesh$loc, lwd=lwd+1, ...)
if (!is.null(mesh$segm$int))
lines(mesh$segm$int, mesh$loc, lwd=lwd+1, ...)
}
if (!add && missing(main)) {
if (mesh$meta$is.refined) {
title("Constrained refined Delaunay triangulation", ...)
} else {
title("Constrained Delaunay triangulation", ...)
}
}
return(invisible())
}
}
inla.mesh.map.lim <-
function(loc=NULL,
projection=
c("default", "longlat", "longsinlat", "mollweide"))
{
projection = match.arg(projection)
if (identical(projection, "default")) {
if (is.null(loc)) {
lim = list(xlim=c(0, 1), ylim=c(0, 1))
} else {
lim =
list(xlim=range(loc[,1], na.rm=TRUE),
ylim=range(loc[,2], na.rm=TRUE))
}
} else if (identical(projection, "longlat")) {
lim = list(xlim=c(-180,180), ylim=c(-90,90))
} else if (identical(projection, "longsinlat")) {
lim = list(xlim=c(-180,180), ylim=c(-1,1))
} else if (identical(projection, "mollweide")) {
lim = list(xlim=c(-2,2), ylim=c(-1,1))
} else {
stop(paste("Unknown projection '", projection, "'.", sep=""))
}
return(lim)
}
inla.mesh.map <-
function(loc,
projection=
c("default", "longlat", "longsinlat", "mollweide"),
inverse=TRUE)
{
projection = match.arg(projection)
if (identical(projection, "default")) {
return(loc)
} else if (identical(projection, "longlat")) {
if (inverse) {
proj =
cbind(cos(loc[,1]*pi/180)*cos(loc[,2]*pi/180),
sin(loc[,1]*pi/180)*cos(loc[,2]*pi/180),
sin(loc[,2]*pi/180))
} else {
proj =
cbind(atan2(loc[,2], loc[,1])*180/pi,
asin(pmax(-1, pmin(+1, loc[,3])))*180/pi)
}
} else if (identical(projection, "longsinlat")) {
if (inverse) {
coslat = sqrt(pmax(0, 1-loc[,2]^2))
proj =
cbind(cos(loc[,1]*pi/180)*coslat,
sin(loc[,1]*pi/180)*coslat,
loc[,2]
)
} else {
proj =
cbind(atan2(loc[,2], loc[,1])*180/pi,
loc[,3])
}
} else if (identical(projection, "mollweide")) {
if (inverse) {
ok = ((loc[,1]^2+4*loc[,2]^2) <= 4)
cos.theta = sqrt(pmax(0, 1-loc[ok,2]^2))
theta = atan2(loc[ok,2], cos.theta)
sin.lat = (2*theta+sin(2*theta))/pi
cos.lat = sqrt(pmax(0, 1-sin.lat^2))
lon = loc[ok,1]*pi/2/(cos.theta+(cos.theta==0))
lon[cos.theta==0] = pi/2*sign(theta[cos.theta==0])
proj = matrix(NA, nrow(loc), 3)
proj[ok,] = cbind(cos(lon)*cos.lat, sin(lon)*cos.lat, sin.lat)
} else {
lon = atan2(loc[,2],loc[,1])
z = pmin(1, pmax(-1, loc[,3]))
sin.theta = z
cos.theta = sqrt(pmax(0, 1-sin.theta^2))
## NR-solver for sin.theta.
## Typically finishes after at most 7 iterations.
## When cos.theta=0, sin.theta is already correct, +/- 1.
nook = (cos.theta > 0)
for (k in 1:20) {
if (any(nook)) {
delta =
(atan2(sin.theta[nook], cos.theta[nook]) +
sin.theta[nook]*cos.theta[nook] - pi/2*z[nook])/
(2*cos.theta[nook])
sin.theta[nook] = sin.theta[nook] - delta
cos.theta[nook] = sqrt(1-sin.theta[nook]^2)
nook[nook] = (abs(delta) > 1e-14)
}
}
proj = cbind(2*lon/pi*cos.theta, sin.theta)
}
} else {
stop(paste("Unknown projection '", projection, "'.", sep=""))
}
return(proj)
}
inla.mesh.lattice <- function(x=seq(0, 1, length.out=2),
y=seq(0, 1, length.out=2),
z=NULL,
dims =
if (is.matrix(x)) {
dim(x)
} else {
c(length(x), length(y))
},
units = NULL,
crs = NULL)
{
if (is.null(crs)) {
units = match.arg(units, c("default", "longlat", "longsinlat", "mollweide"))
lim = inla.mesh.map.lim(projection=units)
xlim = lim$xlim
ylim = lim$ylim
} else { ## !is.null(crs)
if (!is.null(units)) {
stop("Only one of 'units' and 'crs' can be non-null.")
}
bounds <- inla.crs.bounds(crs)
xlim = bounds$xlim
ylim = bounds$ylim
}
if (missing(x) && !missing(dims)) {
x = seq(xlim[1], xlim[2], length.out=dims[1])
}
if (missing(y) && !missing(dims)) {
y = seq(ylim[1], ylim[2], length.out=dims[2])
}
dims = as.integer(dims)
if (is.matrix(x)) {
if (!identical(dims, dim(x)) ||
!identical(dims, dim(y)) ||
(is.matrix(z) && !identical(dims, dim(z))))
stop("The size of matrices 'x', 'y', and 'z' must match 'dims'.")
loc = cbind(as.vector(x), as.vector(y), as.vector(z))
x = NULL
y = NULL
} else {
if (!identical(dims[1], length(x)) ||
!identical(dims[2], length(y)))
stop(paste("The lengths of vectors 'x' and 'y' (",
length(x),",",length(y),
") must match 'dims' (",dims[1],",",dims[2],").",
sep=""))
loc = (cbind(rep(x, times = dims[2]),
rep(y, each = dims[1])))
}
if (!is.double(loc))
storage.mode(loc) = "double"
if (is.null(crs)) {
loc = inla.mesh.map(loc=loc, projection=units, inverse=TRUE)
}
## Construct lattice boundary
segm.idx = (c(1:(dims[1]-1),
dims[1]*(1:(dims[2]-1)),
dims[1]*dims[2]-(0:(dims[1]-2)),
dims[1]*((dims[2]-1):1)+1))
segm.grp = (c(rep(1L, dims[1]-1),
rep(2L, dims[2]-1),
rep(3L, dims[1]-1),
rep(4L, dims[2]-1)))
segm = (inla.mesh.segment(loc=loc[segm.idx,, drop=FALSE],
grp=segm.grp,
is.bnd=TRUE,
crs=crs))
lattice = list(dims=dims, x=x, y=y, loc=loc, segm=segm, crs=crs)
class(lattice) = "inla.mesh.lattice"
return(lattice)
}
extract.groups <- function(...)
{
UseMethod("extract.groups")
}
extract.groups.inla.mesh.segment <- function(segm,
groups,
groups.new=groups,
...)
{
inla.require.inherits(segm, "inla.mesh.segment", "'segm'")
if (length(groups.new)==1L) {
groups.new = rep(groups.new, length(groups))
}
if (length(groups.new)!=length(groups)) {
stop("Length of 'groups.new' (", length(groups.new),
") does not match length of 'groups' (",length(groups),")")
}
idx = c()
segm.grp = c()
for (k in 1:length(groups)) {
extract.idx = which(segm$grp==groups[k])
idx = c(idx, extract.idx)
segm.grp = c(segm.grp, rep(groups.new[k], length(extract.idx)))
}
segm.idx = segm$idx[idx,, drop=FALSE]
return(inla.mesh.segment(loc=segm$loc,
idx=segm.idx,
grp=segm.grp,
segm$is.bnd,
crs=segm$crs))
}
inla.mesh.parse.segm.input <- function(boundary=NULL,
interior=NULL,
segm.offset=0L,
loc.offset=0L,
crs=NULL)
{
###########################################
homogenise.segm.input <- function(x, is.bnd, crs=NULL)
{
if (is.matrix(x) || is.vector(x)) { ## Coordinates or indices
x = (inla.ifelse(is.matrix(x),x,
as.matrix(x,nrow=length(x),ncol=1)))
if (is.integer(x)) { ## Indices
ret = inla.mesh.segment(NULL, x, NULL, is.bnd, crs=crs)
} else if (is.numeric(x)) { ## Coordinates
ret = inla.mesh.segment(x, NULL, NULL, is.bnd, crs=crs)
} else {
stop("Segment info matrix must be numeric or integer.")
}
} else if (inherits(x, "inla.mesh.segment")) {
## Override x$is.bnd:
ret = inla.mesh.segment(x$loc, x$idx, x$grp, is.bnd, x$crs)
} else if (inherits(x, c("SpatialPoints", "SpatialPointsDataFrame",
"Line", "Lines",
"SpatialLines", "SpatialLinesDataFrame",
"Polygon", "Polygons",
"SpatialPolygons",
"SpatialPolygonsDataFrame"))) {
x <- as.inla.mesh.segment(x)
## Override x$is.bnd:
ret <- inla.mesh.segment(x$loc, x$idx, x$grp, is.bnd, x$crs)
} else if (!is.null(x)) {
inla.require.inherits(NULL,
c("matrix", "inla.mesh.segment"),
"Segment info")
} else {
ret = NULL
}
ret
}
##################################################
homogenise.segm.grp <- function(input) {
grp.idx = 0L
for (k in 1:length(input)) if (!is.null(input[[k]])) {
inla.require.inherits(input[[k]],
"inla.mesh.segment",
"Segment info list members ")
if (is.null(input[[k]]$grp)) {
grp.idx = grp.idx+1L
input[[k]] = (inla.mesh.segment(input[[k]]$loc,
input[[k]]$idx,
grp.idx,
input[[k]]$is.bnd,
input[[k]]$crs))
} else {
grp.idx = max(grp.idx, input[[k]]$grp, na.rm=TRUE)
}
}
return(input)
}
##################################################
join.segm.input <- function(segm, segm.offset=0L, loc.offset=0L,
crs=NULL)
{
loc = NULL
bnd = list(loc=NULL, idx = matrix(,0,2), grp = matrix(,0,1))
int = list(loc=NULL, idx = matrix(,0,2), grp = matrix(,0,1))
storage.mode(bnd$idx) <- "integer"
storage.mode(bnd$grp) <- "integer"
storage.mode(int$idx) <- "integer"
storage.mode(int$grp) <- "integer"
for (k in seq_along(segm)) {
if (!is.null(segm[[k]])) {
inla.require.inherits(segm[[k]],
"inla.mesh.segment",
"Segment info list members ")
if (!is.null(segm[[k]]$loc)) {
local.loc <- inla.spTransform(segm[[k]], crs,
passthrough=TRUE)$loc
extra.loc.n = nrow(local.loc)
idx.offset = segm.offset
segm.offset = segm.offset + extra.loc.n
loc = (inla.ifelse(is.null(loc),
local.loc,
rbind(loc, local.loc)))
} else {
idx.offset = loc.offset
}
if (segm[[k]]$is.bnd) {
bnd$idx = rbind(bnd$idx, segm[[k]]$idx+idx.offset)
bnd$grp = rbind(bnd$grp, segm[[k]]$grp)
} else {
int$idx = rbind(int$idx, segm[[k]]$idx+idx.offset)
int$grp = rbind(int$grp, segm[[k]]$grp)
}
}
}
if (nrow(bnd$idx)==0)
bnd = NULL
if (nrow(int$idx)==0)
int = NULL
return(list(loc = loc,
bnd = (inla.ifelse(is.null(bnd),
NULL,
inla.mesh.segment(bnd$loc,
bnd$idx,
bnd$grp,
TRUE,
crs))),
int = (inla.ifelse(is.null(int),
NULL,
inla.mesh.segment(int$loc,
int$idx,
int$grp,
FALSE,
crs)))))
}
###########################
segm = (c(lapply(inla.ifelse(inherits(boundary, "list"),
boundary, list(boundary)),
function(x){homogenise.segm.input(x, TRUE, crs=crs)}),
lapply(inla.ifelse(inherits(interior, "list"),
interior, list(interior)),
function(x){homogenise.segm.input(x, FALSE, crs=crs)})))
segm = homogenise.segm.grp(segm)
join.segm.input(segm=segm, segm.offset=segm.offset, loc.offset=loc.offset,
crs=crs)
}
################################
##
## Old code. Filtering is now done in fmesher itself.
## Retained for now so that we can check if the results are the same.
##
inla.mesh.filter.locations <- function(loc, cutoff)
{
## Map locations to nodes, avoiding near-duplicates.
loc.n = nrow(loc)
loc.dim = ncol(loc)
loc.is.na = (rowSums(is.na(loc))>0)
if (sum(loc.is.na)>0)
stop("NAs in locations not yet supported.")
node.coord = matrix(nrow=loc.n, ncol=loc.dim)
map.loc.to.node = rep(0L, nrow=loc.n)
excluded = c()
loc.i = 1L
node.i.max = 1L
node.coord[node.i.max,] = loc[loc.i,]
map.loc.to.node[[loc.i]] = node.i.max
for (loc.i in 2:loc.n) {
loc.to.node.dist =
sqrt(rowSums((as.matrix(rep(1, node.i.max)) %*%
loc[loc.i,, drop=FALSE] -
node.coord[1:node.i.max,, drop=FALSE])^2))
if (min(loc.to.node.dist) > cutoff) {
node.i.max = node.i.max+1L
node.coord[node.i.max,] = loc[loc.i,, drop=FALSE]
map.loc.to.node[[loc.i]] = node.i.max
} else {
excluded = c(excluded, loc.i)
}
}
## Remove excess nodes.
node.coord = node.coord[1:node.i.max,]
## Identify nearest nodes for excluded locations.
for (loc.i in excluded) {
loc.to.node.dist =
sqrt(rowSums((as.matrix(rep(1, node.i.max)) %*%
loc[loc.i,, drop=FALSE] -
node.coord)^2))
node.i = which.min(loc.to.node.dist);
map.loc.to.node[[loc.i]] = node.i
}
return(list(loc = node.coord, node.idx = map.loc.to.node))
}
############################
inla.mesh <- function(...)
{
args = list(...)
if (length(args)>0) {
if (inherits(args[[1]],"inla.mesh")) {
warning("'inla.mesh(mesh, ...)' is deprecated. Use 'inla.mesh.query(mesh, ...)' instead.")
return(inla.mesh.query(...))
}
}
warning("'inla.mesh(...)' is deprecated. Use 'inla.mesh.create(...)' instead.")
return(inla.mesh.create(...))
}
inla.mesh.create <- function(loc=NULL, tv=NULL,
boundary=NULL, interior=NULL,
extend = (missing(tv) || is.null(tv)),
refine=FALSE,
lattice=NULL,
globe=NULL,
cutoff = 1e-12,
plot.delay = NULL,
data.dir,
keep = (!missing(data.dir) && !is.null(data.dir)),
timings = FALSE,
quality.spec=NULL,
crs=NULL)
{
if (!timings) {
system.time <- function(expr) {
expr
structure(c(0,0,0,0,0), class = "proc_time")
}
}
time.pre = system.time({ ## Pre-processing timing start
if (!is.null(crs) &&
identical(inla.as.list.CRS(crs)$proj, "geocent")) {
## Build all geocentric meshes on a sphere, and transform afterwards,
## to allow general geoids.
crs.target <- crs
crs <- inla.CRS("sphere")
}
if (!(missing(loc) || is.null(loc))) {
## Handle loc given as SpatialPoints or SpatialPointsDataFrame object
if (inherits(loc, "SpatialPoints") ||
inherits(loc, "SpatialPointsDataFrame")) {
loc <- inla.spTransform(coordinates(loc),
CRS(proj4string(loc)),
crs,
passthrough=TRUE)
}
if (!is.matrix(loc)) {
loc = as.matrix(loc)
}
if (!is.double(loc)) {
storage.mode(loc) = "double"
}
}
loc.n = max(0L,nrow(loc))
if (is.logical(extend) && extend) extend = list()
if (is.logical(refine) && refine) refine = list()
if (missing(lattice) || is.null(lattice)) {
lattice = list(loc=NULL, segm=NULL)
lattice.n = 0L
} else {
inla.require.inherits(lattice, "inla.mesh.lattice", "'lattice'")
if (!is.null(tv)) {
warning("Both 'lattice' and 'tv' specified. Ignoring 'tv'.")
tv = NULL
}
if (!inherits(extend, "list")) {
boundary = (c(inla.ifelse(inherits(boundary, "list"),
boundary, list(boundary)),
list(lattice$segm)))
}
if (!is.null(lattice[["crs"]])) {
lattice$loc <- inla.spTransform(lattice$loc,
lattice$crs,
crs,
passthrough=TRUE)
}
if (!is.matrix(lattice$loc)) {
lattice$loc = as.matrix(lattice$loc)
}
if (!is.double(lattice$loc)) {
storage.mode(lattice$loc) = "double"
}
}
lattice.n = max(0L,nrow(lattice$loc))
segm = (inla.mesh.parse.segm.input(boundary,
interior,
loc.n,
0L,
crs=crs))
segm.n = max(0,nrow(segm$loc))
## Run parse again now that we know where the indices should point:
segm = (inla.mesh.parse.segm.input(boundary,
interior,
0L,
segm.n+lattice.n,
crs=crs))
loc0 = rbind(segm$loc, lattice$loc, loc)
if ((!is.null(loc0)) && (nrow(loc0)>0))
idx0 = 1:nrow(loc0)
else
idx0 = c()
if (is.null(quality.spec)) {
quality <- NULL
} else {
quality <- rep(NA, segm.n+lattice.n+loc.n)
## Order must be same as for loc0 above.
if (!is.null(quality.spec$segm)) {
quality[seq_len(segm.n)] <- quality.spec$segm
}
if (!is.null(quality.spec$lattice)) {
quality[segm.n+seq_len(lattice.n)] <- quality.spec$lattice
}
if (!is.null(quality.spec$loc)) {
quality[segm.n+lattice.n+seq_len(loc.n)] <- quality.spec$loc
}
## NA:s will be replaced with max.edge settings below.
}
## Where to put the files?
if (keep) {
if (missing(data.dir)) {
data.dir = inla.fmesher.make.dir("inla.mesh.data")
}
prefix = paste(data.dir, "/mesh.", sep="")
keep.dir = TRUE
} else {
if (missing(data.dir)) {
prefix = inla.tempfile(pattern="fmesher", tmpdir=tempdir())
prefix = paste(prefix, ".", sep="")
keep.dir = TRUE
} else {
data.dir = inla.fmesher.make.dir(data.dir)
prefix = paste(data.dir, "/mesh.", sep="")
## We should not try to delete the session tempdir...
keep.dir = identical(tempdir(), dirname(prefix))
}
}
prefix = inla.fmesher.make.prefix(NULL, prefix)
if (is.null(loc0)) {
all.args = ""
} else {
loc.file = fmesher.write(inla.affirm.double(loc0), prefix, "input.s")
all.args = "--input=input.s"
if (!is.null(tv)) {
fmesher.write(inla.affirm.integer(tv)-1L, prefix, "input.tv")
all.args = paste(all.args, ",input.tv", sep="")
}
}
if (!missing(cutoff)) {
all.args = paste(all.args, " --cutoff=", cutoff, sep="")
}
if (!is.null(segm$bnd)) {
fmesher.write(inla.affirm.integer(segm$bnd$idx)-1L, prefix, "input.segm.bnd.idx")
fmesher.write(inla.affirm.integer(segm$bnd$grp), prefix, "input.segm.bnd.grp")
all.args = paste(all.args," --boundary=input.segm.bnd.idx")
all.args = paste(all.args," --boundarygrp=input.segm.bnd.grp")
}
if (!is.null(segm$int)) {
fmesher.write(inla.affirm.integer(segm$int$idx)-1L, prefix, "input.segm.int.idx")
fmesher.write(inla.affirm.integer(segm$int$grp), prefix, "input.segm.int.grp")
all.args = paste(all.args," --interior=input.segm.int.idx")
all.args = paste(all.args," --interiorgrp=input.segm.int.grp")
}
if (!missing(globe) && !is.null(globe)) {
all.args = paste(all.args, " --globe=", globe, sep="")
}
if (inherits(extend,"list")) {
cet = c(0,0)
cet[1] = inla.ifelse(is.null(extend$n), 16, extend$n)
cet[2] = inla.ifelse(is.null(extend$offset), -0.1, extend$offset)
all.args = (paste(all.args," --cet=",
cet[1],",", cet[2], sep=""))
}
if (inherits(refine,"list")) {
rcdt = c(0,0,0)
if (!missing(globe) && !is.null(globe)) {
max.edge.default <- pi
} else {
## Multiply by to ensure cover S2 domains
max.edge.default <- inla.diameter(loc0) * 2
}
if ((inherits(extend,"list")) && (!is.null(extend$offset))) {
max.edge.default = (max.edge.default +
max(0,2*extend$offset))
max.edge.default = (max.edge.default *
(1+max(0,-2*extend$offset)))
}
rcdt[1] = inla.ifelse(is.null(refine$min.angle), 21, refine$min.angle)
rcdt[2] = (inla.ifelse(is.null(refine$max.edge) ||
is.na(refine$max.edge),
max.edge.default,
refine$max.edge))
rcdt[3] = (inla.ifelse(is.null(refine$max.edge) ||
is.na(refine$max.edge),
max.edge.default,
refine$max.edge))
rcdt[2] = (inla.ifelse(is.null(refine$max.edge.extra) ||
is.na(refine$max.edge.extra),
rcdt[2], refine$max.edge.extra))
rcdt[3] = (inla.ifelse(is.null(refine$max.edge.data) ||
is.na(refine$max.edge.date),
rcdt[3], refine$max.edge.data))
all.args = (paste(all.args," --rcdt=",
rcdt[1],",", rcdt[2],",", rcdt[3], sep=""))
if (!is.null(refine[["max.n.strict"]]) &&
!is.na(refine$max.n.strict)) {
rcdt_max_n0 <- refine$max.n.strict
} else {
rcdt_max_n0 <- -1
}
if (!is.null(refine[["max.n"]]) &&
!is.na(refine$max.n)) {
rcdt_max_n1 <- refine$max.n
} else {
rcdt_max_n1 <- -1
}
all.args = (paste(all.args,
" --max_n0=", rcdt_max_n0,
" --max_n1=", rcdt_max_n1,
sep=""))
is.refined = TRUE
if (!is.null(quality)) {
quality[is.na(quality)] <- rcdt[3]
quality.file <-
fmesher.write(inla.affirm.double(as.matrix(quality)),
prefix, "input.quality")
all.args <- (paste(all.args,
" --input=input.quality",
" --quality=input.quality",
sep=""))
}
} else {
is.refined = FALSE
}
if (!is.null(plot.delay)) {
all.args = paste(all.args," --x11=", plot.delay, sep="")
}
all.args = paste(all.args, inla.getOption("fmesher.arg"))
}) ## Pre-processing timing end
## Call fmesher:
time.fmesher = system.time({
echoc = (inla.fmesher.call(all.args=all.args,
prefix=prefix))
})
time.post = system.time({ ## Post-processing timing start
## Read the mesh:
manifold = 1L+fmesher.read(prefix, "manifold")
manifold = list("M", "R2", "S2")[[manifold]]
loc = fmesher.read(prefix, "s")
graph = (list(tv = 1L+fmesher.read(prefix, "tv"),
vt = 1L+fmesher.read(prefix, "vt"),
tt = 1L+fmesher.read(prefix, "tt"),
tti = 1L+fmesher.read(prefix, "tti"),
vv = fmesher.read(prefix, "vv")))
graph$tv[graph$tv==0L] = NA
graph$vt[graph$vt==0L] = NA
graph$tt[graph$tt==0L] = NA
graph$tti[graph$tti==0L] = NA
## Read the vertex input/output mapping:
idx.all = 1L+fmesher.read(prefix, "idx")
idx = (list(loc = (inla.ifelse(loc.n>0,
idx.all[idx0[segm.n+lattice.n+(1:loc.n)]],
NULL)),
lattice = (inla.ifelse(lattice.n>0,
idx.all[idx0[segm.n+(1:lattice.n)]],
NULL)),
segm = (inla.ifelse(segm.n>0,
idx.all[idx0[(1:segm.n)]],
NULL))))
if (!is.null(idx$loc)) idx$loc[idx$loc == 0L] = NA
if (!is.null(idx$lattice)) idx$lattice[idx$lattice == 0L] = NA
if (!is.null(idx$segm)) idx$segm[idx$segm == 0L] = NA
## Read constraint segment information:
segm.bnd = (inla.mesh.segment(NULL,
1L+fmesher.read(prefix, "segm.bnd.idx"),
fmesher.read(prefix, "segm.bnd.grp"),
TRUE))
segm.int = (inla.mesh.segment(NULL,
1L+fmesher.read(prefix, "segm.int.idx"),
fmesher.read(prefix, "segm.int.grp"),
FALSE))
## Remap indices to remove unused vertices
used = !is.na(graph$vt)
if (!all(used)) {
used = which(used)
idx.map = rep(NA, nrow(loc))
idx.map[used] = seq_len(length(used))
loc = loc[used,,drop=FALSE]
graph$tv = matrix(idx.map[as.vector(graph$tv)], nrow(graph$tv), 3)
graph$vt = graph$vt[used,,drop=FALSE]
## graph$tt ## No change needed
## graph$tti ## No change needed
graph$vv = graph$vv[used, used, drop=FALSE]
if (!is.null(idx$loc)) idx$loc = idx.map[idx$loc]
if (!is.null(idx$lattice)) idx$lattice = idx.map[idx$lattice]
if (!is.null(idx$segm)) idx$segm = idx.map[idx$segm]
segm.bnd$idx = matrix(idx.map[segm.bnd$idx], nrow(segm.bnd$idx), 2)
segm.int$idx = matrix(idx.map[segm.int$idx], nrow(segm.int$idx), 2)
segm.bnd$idx[segm.bnd$idx == 0L] = NA
segm.int$idx[segm.int$idx == 0L] = NA
}
if (!keep)
unlink(paste(prefix, "*", sep=""), recursive=FALSE)
if (!keep.dir) {
unlink(dirname(prefix), recursive=TRUE)
}
if (!is.null(crs) &&
identical(inla.as.list.CRS(crs)$proj, "geocent")) {
if (!inla.identical.CRS(crs, crs.target)) {
## Target is a non-spherical geoid
loc <- inla.spTransform(loc, crs, crs.target)
crs <- crs.target
}
}
}) ## Post-processing timing end
time.object = system.time({ ## Object construction timing start
mesh = (list(meta = (list(call=match.call(),
fmesher.args = all.args,
time = (rbind(pre = time.pre,
fmesher = time.fmesher,
post = time.post)),
prefix = prefix,
is.refined = is.refined)),
manifold = manifold,
n = nrow(loc),
loc = loc,
graph = graph,
segm = list(bnd=segm.bnd, int=segm.int),
idx = idx,
crs = crs))
class(mesh) <- "inla.mesh"
}) ## Object construction timing end
## Entire function timing
time.total = time.pre+time.fmesher+time.post+time.object
mesh$meta$time = (rbind(mesh$meta$time,
object=time.object,
total=time.total))
return(invisible(mesh))
}
inla.mesh.extract.segments <- function(mesh.loc,
mesh.idx,
mesh.grp,
grp=NULL,
is.bnd,
crs=NULL)
{
segments = list()
if (nrow(mesh.idx)>0) {
if (is.null(grp)) {
grp = unique(sort(mesh.grp))
}
for (g in grp) {
extract = (mesh.grp==g)
segments =
c(segments,
list(inla.mesh.segment(mesh.loc,
idx=mesh.idx[extract,,drop=FALSE],
grp=mesh.grp[extract,drop=FALSE],
is.bnd=is.bnd,
crs=crs)) )
}
}
if (length(segments)>0)
return(segments)
else
return(NULL)
}
inla.mesh.boundary <- function(mesh, grp=NULL)
{
inla.require.inherits(mesh, "inla.mesh", "'mesh'")
return(inla.mesh.extract.segments(mesh$loc,
mesh$segm$bnd$idx,
mesh$segm$bnd$grp,
grp,
TRUE,
mesh$crs))
}
inla.mesh.interior <- function(mesh, grp=NULL)
{
inla.require.inherits(mesh, "inla.mesh", "'mesh'")
return(inla.mesh.extract.segments(mesh$loc,
mesh$segm$int$idx,
mesh$segm$int$grp,
grp,
FALSE,
mesh$crs))
}
## Generate nice triangulation, with an inner domain strictly enclosed
## by an outer domain
## The inner and outer domains can have different quality parameters
## At least one of loc, loc.domain, boundary[[1]], boundary[[2]], interior
## must be non-NULL
## For more complicated multi-step meshings, study the code and write your own.
inla.mesh.2d <-
function(loc=NULL, ## Points to include in final triangulation
loc.domain=NULL, ## Points that determine the automatic domain
offset=NULL, ## Size of automatic extensions
n=NULL, ## Sides of automatic extension polygons
boundary=NULL, ## User-specified domains (list of length 2)
interior=NULL, ## User-specified constraints for the inner domain
max.edge=NULL,
min.angle=NULL, ## Angle constraint for the entire domain
cutoff=1e-12, ## Only add input points further apart than this
max.n.strict=NULL,
max.n=NULL,
plot.delay=NULL,
crs=NULL) ## Coordinate Reference System
## plot.delay: Do plotting.
## NULL --> No plotting
## <0 --> Intermediate meshes displayed at the end
## >0 --> Dynamical fmesher plotting
{
###########################
unify.one.segm <- function(segm, crs=NULL) {
if (inherits(segm, "inla.mesh.segment")) {
segm <- inla.spTransform(segm, crs, passthrough=TRUE)
} else if (inherits(segm, "matrix")) {
segm <- inla.mesh.segment(loc=segm, crs=crs)
} else {
segm <- inla.spTransform(as.inla.mesh.segment(segm), crs,
passthrough=TRUE)
}
if (ncol(segm$loc)==2) {
segm$loc = cbind(segm$loc, 0.0)
}
segm
}
unify.segm.coords <- function(segm, crs=NULL) {
if (!is.null(segm)) {
if (inherits(segm, "list")) {
for (j in seq_along(segm)) {
segm[[j]] <- unify.one.segm(segm[[j]], crs=crs)
}
} else {
segm <- unify.one.segm(segm, crs=crs)
}
}
segm
}
###########################
if ((missing(max.edge) || is.null(max.edge)) &&
(missing(max.n.strict) || is.null(max.n.strict)) &&
(missing(max.n) || is.null(max.n))) {
stop("At least one of max.edge, max.n.strict, and max.n must be specified")
}
if (!is.null(crs)) {
issphere <- inla.identical.CRS(crs, inla.CRS("sphere"))
isgeocentric <- identical(inla.as.list.CRS(crs)[["proj"]], "geocent")
if (isgeocentric) {
crs.target <- crs
crs <- inla.CRS("sphere")
}
}
## Handle loc given as SpatialPoints or SpatialPointsDataFrame object
if (!(missing(loc) || is.null(loc)) &&
(inherits(loc, "SpatialPoints") ||
inherits(loc, "SpatialPointsDataFrame"))) {
loc = inla.spTransform(coordinates(loc),
CRS(proj4string(loc)),
crs,
passthrough=TRUE)
}
if (!(missing(loc.domain) || is.null(loc.domain)) &&
(inherits(loc.domain, "SpatialPoints") ||
inherits(loc.domain, "SpatialPointsDataFrame"))) {
loc.domain = inla.spTransform(coordinates(loc.domain),
CRS(proj4string(loc.domain)),
crs,
passthrough=TRUE)
}
if (missing(loc) || is.null(loc)) {
loc = matrix(c(0.0), 0, 3)
} else if (!is.matrix(loc)) {
loc = as.matrix(loc)
}
if (missing(loc.domain) || is.null(loc.domain)) {
loc.domain = loc
} else if (!is.matrix(loc.domain)) {
loc.domain = as.matrix(loc.domain)
}
if (missing(boundary)) {
boundary = list(NULL)
} else {
if (!inherits(boundary, "list"))
boundary = list(boundary)
}
if (missing(interior))
interior = NULL
if (missing(offset) || is.null(offset)) {
if (length(boundary) < 2)
offset = -0.05
else
offset = c(-0.05, -0.15)
}
if (missing(n) || is.null(n))
n = c(8)
if (missing(max.edge) || is.null(max.edge))
max.edge = c(NA)
if (missing(min.angle) || is.null(min.angle))
min.angle = c(21)
if (missing(max.n.strict) || is.null(max.n.strict))
max.n.strict = c(NA)
if (missing(max.n) || is.null(max.n))
max.n = c(NA)
if (missing(cutoff) || is.null(cutoff))
cutoff = 1e-12
if (missing(plot.delay) || is.null(plot.delay))
plot.delay = NULL
num.layers =
max(c(length(boundary), length(offset), length(n),
length(min.angle), length(max.edge),
length(max.n.strict), length(max.n)))
if (num.layers > 2) {
warning(paste("num.layers=", num.layers, " > 2 detected. ",
"Excess information ignored.", sep=""))
num.layers = 2
}
if (length(boundary) < num.layers)
boundary = c(boundary, list(NULL))
if (length(min.angle) < num.layers)
min.angle = c(min.angle, min.angle)
if (length(max.n.strict) < num.layers)
max.n0 = c(max.n.strict, max.n.strict)
if (length(max.n) < num.layers)
max.n = c(max.n, max.n)
if (length(max.edge) < num.layers)
max.edge = c(max.edge, max.edge)
if (length(offset) < num.layers)
offset = c(offset, -0.15)
if (length(n) < num.layers)
n = c(n, 16)
if (length(n) < num.layers)
n = c(n, 16)
## Unify the dimensionality of the point input.
if (!is.null(loc) && (ncol(loc)==2))
loc = cbind(loc, 0.0)
if (!is.null(loc.domain) && (ncol(loc.domain)==2))
loc.domain = cbind(loc.domain, 0.0)
## Unify the dimensionality of the boundary&interior segments input
## and optionally transform coordinates.
for (k in seq_len(num.layers)) {
if (!is.null(boundary[[k]])) {
boundary[[k]] <- unify.segm.coords(boundary[[k]], crs=crs)
}
}
interior <- unify.segm.coords(interior, crs=crs)
## Triangulate to get inner domain boundary
## Constraints included only to get proper domain extent
## First, attach the loc points to the domain definition set
if (!is.null(loc) && !is.null(loc.domain)) {
loc.domain = rbind(loc.domain, loc)
}
mesh1 =
inla.mesh.create(loc=loc.domain,
boundary=boundary[[1]],
interior=interior,
cutoff=cutoff,
extend=list(n=n[1], offset=offset[1]),
refine=FALSE,
plot.delay=plot.delay,
crs=crs)
## Save the resulting boundary
boundary1 = inla.mesh.boundary(mesh1)
interior1 = inla.mesh.interior(mesh1)
if (!is.null(plot.delay) && (plot.delay<0)) {
inla.dev.new()
plot(mesh1)
}
## Triangulate inner domain
mesh2 =
inla.mesh.create(loc=loc,
boundary=boundary1,
interior=interior1,
cutoff=cutoff,
extend=FALSE, ## Should have no effect
refine=
list(min.angle=min.angle[1],
max.edge=max.edge[1],
max.edge.extra=max.edge[1],
max.n.strict=max.n.strict[1],
max.n=max.n[1]),
plot.delay=plot.delay,
crs=crs)
boundary2 = inla.mesh.boundary(mesh2)
interior2 = inla.mesh.interior(mesh2)
if (!is.null(plot.delay) && (plot.delay<0)) {
inla.dev.new()
plot(mesh2)
}
if (num.layers == 1) {
if (!is.null(crs) && isgeocentric && !issphere) {
mesh2$loc <- inla.spTransform(mesh2$loc, mesh2$crs, crs.target)
mesh2$crs <- crs.target
}
return(invisible(mesh2))
}
## Triangulate inner+outer domain
mesh3 =
inla.mesh.create(loc=rbind(loc, mesh2$loc),
boundary=boundary[[2]],
interior=c(boundary2, interior2),
cutoff=cutoff,
extend=list(n=n[2], offset=offset[2]),
refine=
list(min.angle=min.angle[2],
max.edge=max.edge[2],
max.edge.extra=max.edge[2],
max.n.strict=mesh2$n + max.n.strict[2],
max.n=mesh2$n + max.n[2]),
plot.delay=plot.delay,
crs=crs)
## Hide generated points, to match regular inla.mesh.create output
mesh3$idx$loc = mesh3$idx$loc[seq_len(nrow(loc))]
## Obtain the corresponding segm indices.
segm.loc = matrix(0.0, 0, 3)
for (k in seq_along(boundary)) {
if (!is.null(boundary[[k]])) {
segm.loc = rbind(segm.loc, boundary[[k]]$loc)
}
}
for (k in seq_along(interior)) {
if (!is.null(interior[[k]])) {
segm.loc = rbind(segm.loc, interior[[k]]$loc)
}
}
if (nrow(segm.loc) > 0) {
proj = inla.mesh.project(mesh3, loc=segm.loc)
mesh3$idx$segm = rep(NA, nrow(segm.loc))
if (any(proj$ok)) {
t.idx <- proj$t[proj$ok]
tv.idx <- max.col(proj$bary[proj$ok,,drop=FALSE],
ties.method="first")
mesh3$idx$segm[proj$ok] <-
mesh3$graph$tv[t.idx + nrow(mesh3$graph$tv)*(tv.idx-1)]
}
} else {
mesh3$idx$segm = NULL
}
if (!is.null(crs) && isgeocentric && !issphere) {
mesh3$loc <- inla.spTransform(mesh3$loc, mesh3$crs, crs.target)
mesh3$crs <- crs.target
}
if (!is.null(plot.delay) && (plot.delay<0)) {
inla.dev.new()
plot(mesh3)
}
return(invisible(mesh3))
}
## Support for legacy code:
inla.mesh.create.helper <- function(points=NULL, points.domain=NULL, ...)
{
return(invisible(inla.mesh.2d(loc=points, loc.domain=points.domain, ...)))
}
## Example:
##if (FALSE) {
## mesh =
## inla.mesh.create.helper(loc=matrix(runif(20),10,2)*200,
## n=c(8,16),
## offset=c(10,140),
## max.edge=c(25,1000),
## min.angle=26,
## cutoff=0,
## plot.delay=-1
## )
##}
##
inla.delaunay <- function(loc, ...)
{
## Handle loc given as SpatialPoints or SpatialPointsDataFrame object
if (!(missing(loc) || is.null(loc)) &&
(inherits(loc, "SpatialPoints") ||
inherits(loc, "SpatialPointsDataFrame"))) {
crs = CRS(proj4string(loc))
loc = coordinates(loc)
} else {
crs = NULL
}
hull = chull(loc[,1],loc[,2])
bnd = inla.mesh.segment(loc=loc[hull[length(hull):1],],is.bnd=TRUE)
mesh =
inla.mesh.create(loc=loc,
boundary=bnd,
extend=list(n=3),
refine=FALSE,
crs=crs,
...)
return(invisible(mesh))
}
inla.mesh.query <- function(mesh, ...)
{
inla.require.inherits(mesh, "inla.mesh", "'mesh'")
not.known <- function(mesh, queryname)
{
stop(paste("Query '", queryname,
"' unknown.", sep=""))
}
not.implemented <- function(mesh, queryname)
{
stop(paste("Query '", queryname,
"' not implemented for inla.mesh.", sep=""))
## stop(paste("Query '", queryname,
## "' not implemented for inla.mesh for mesh type '",
## model$type, "'.", sep=""))
}
result = list()
queries = inla.parse.queries(...)
if (length(queries)==0L)
return(result)
for (query.idx in 1:length(queries)) {
query = names(queries)[query.idx]
param = queries[[query.idx]]
answer = NULL
query = (match.arg(query, c("tt.neighbours",
"vt.neighbours"
)))
if (identical(query, "tt.neighbours")) {
if (is.null(param))
param = list(c(1))
if (length(param)<2)
param = c(as.list(param), list(c(1)))
nT = nrow(mesh$graph$tv)
i = rep(1:nT,3)
j = as.vector(mesh$graph$tt)
i = i[!is.na(j)]
j = j[!is.na(j)]
tt = sparseMatrix(i=i, j=j, x=rep(1,length(i)), dims=c(nT,nT))
answer = (sparseMatrix(i=param[[1]],
j=rep(1,length(param[[1]])),
x=1,
dims=c(nT,1)))
tt0 = (answer == 0.5)*1
tt1 = answer
order = 0
while (order<min(param[[2]])) {
order = order+1
tt0 = tt1
tt1 = answer
answer = ((((tt %*% answer) > 0) - tt1 - tt0) > 0)
}
while (order<max(param[[2]])) {
order = order+1
answer = ((((tt %*% answer) > 0) - tt1 - tt0) > 0)
}
## not.implemented(mesh,query)
} else if (identical(query, "vt.neighbours")) {
if (is.null(param))
param = list(1)
if (length(param)<2)
param = c(as.list(param), list(c(1)))
nV = nrow(mesh$loc)
nT = nrow(mesh$graph$tv)
i = rep(1:nT,3)
j = as.vector(mesh$graph$tv)
# i = i[!is.na(j)]
# j = j[!is.na(j)]
tv = sparseMatrix(i=i, j=j, x=rep(1,length(i)), dims=c(nT,nV))
vv = (sparseMatrix(i=param[[1]],
j=rep(1,length(param[[1]])),
x=1,
dims=c(nV,1)))
vt = (tv %*% vv ) > 0
vv0 = (vv == 0.5)*1
vv1 = vv
vt0 = (tv %*% vv0 ) > 0
vt1 = vt
order = 0
while (order<min(param[[2]])) {
order = order+1
vv0 = vv1
vv1 = vv
vv = ((((mesh$graph$vv %*% vv) > 0) - vv1 - vv0) > 0)
vt0 = vt1
vt1 = vt
vt = (((tv %*% vv) > 0) - vt1 - vt0 ) > 0
}
while (order<max(param[[2]])) {
order = order+1
vv = ((((mesh$graph$vv %*% vv) > 0) - vv1 - vv0) > 0)
vt = (((((tv %*% vv) > 0) - vt1 - vt0 ) > 0 ) + vt) > 0
}
answer = vt
## not.implemented(mesh,query)
} else if (!identical(query, "")) {
not.known(mesh,query)
}
## Expand the result list:
result[query.idx] = list(NULL)
names(result)[query.idx] = query
## Set the answer:
if (!is.null(answer))
result[[query.idx]] = answer
}
return(result)
}
summary.inla.mesh <- function(object, verbose=FALSE, ...)
{
x=object
## provides a summary for a mesh object
inla.require.inherits(x, "inla.mesh", "'x'")
ret = list(verbose=verbose)
if (verbose) {
ret = (c(ret, list(call=x$meta$call,
fmesher.args=x$meta$fmesher.args,
prefix=x$meta$prefix,
time = x$meta$time,
is.refined = x$meta$is.refined)))
}
ret = (c(ret, list(manifold=x$manifold,
nV=x$n,
nT=nrow(x$graph$tv),
xlim=range(x$loc[,1]),
ylim=range(x$loc[,2]),
zlim=range(x$loc[,3]))))
if (is.na(inla.CRSargs(x$crs))) {
ret <- c(ret, list(crs="N/A"))
} else {
ret <- c(ret, list(crs=inla.CRSargs(x$crs)))
}
my.segm <- function(x) {
if (is.null(x))
return(list(n=0, grps=NULL))
n = max(0, nrow(x$idx))
if (max(0, length(unique(x$grp)))>0) {
grps = unique(x$grp)
} else {
grps = NULL
}
return(list(n=n, grps=grps))
}
if(!is.null(x$segm)) {
ret = (c(ret, list(segm.bnd=my.segm(x$segm$bnd),
segm.int=my.segm(x$segm$int))))
} else {
ret = (c(ret, list(segm.bnd=my.segm(NULL),
segm.int=my.segm(NULL))))
}
class(ret) <- "summary.inla.mesh"
return (ret)
}
print.summary.inla.mesh <- function(x, ...)
{
my.print.proc_time <- function (x, ...)
{
if (!is.matrix(x)) {
y = matrix(x,1,5)
} else {
y = x
}
for (k in 1:nrow(x)) {
if (!is.na(y[k,4L]))
y[k,1L] <- y[k,1L] + y[k,4L]
if (!is.na(y[k,5L]))
y[k,2L] <- y[k,2L] + y[k,5L]
}
y <- y[,1L:3L]
colnames(y) <- c(gettext("user"), gettext("system"), gettext("elapsed"))
print(y, ...)
invisible(x)
}
inla.require.inherits(x, "summary.inla.mesh", "'x'")
if (x$verbose) {
cat("\nCall:\n")
print(x$call)
cat("\nfmesher:\t", x$fmesher.args, "\n", sep = "")
cat("prefix:\t\t", x$prefix, "\n", sep = "")
cat("\nTimings:\n")
my.print.proc_time(x$time)
}
cat("\nManifold:\t", x$manifold, "\n", sep="")
cat("CRS:\t", x$crs, "\n", sep="")
if (x$verbose) {
cat("Refined:\t", x$is.refined, "\n", sep="")
}
cat("Vertices:\t", as.character(x$nV), "\n", sep="")
cat("Triangles:\t", as.character(x$nT), "\n", sep="")
my.print.segm <- function(x) {
cat(as.character(x$n))
if (!is.null(x$grps)) {
n = length(x$grps)
cat(" (", n, " group", inla.ifelse(n==1, "", "s"), sep="")
if (n <= 10) {
cat(":", x$grps, sep=" ")
} else {
cat(":", x$grps[1:10], "...", sep=" ")
}
cat(")")
}
cat("\n", sep="")
return(invisible())
}
cat("Boundary segm.:\t")
my.print.segm(x$segm.bnd)
cat("Interior segm.:\t")
my.print.segm(x$segm.int)
cat("xlim:\t", x$xlim[1], " ", x$xlim[2], "\n", sep="")
cat("ylim:\t", x$ylim[1], " ", x$ylim[2], "\n", sep="")
cat("zlim:\t", x$zlim[1], " ", x$zlim[2], "\n", sep="")
cat("\n")
invisible(x)
}
inla.mesh.project <- function(...)
{
UseMethod("inla.mesh.project")
}
inla.mesh.project.inla.mesh <- function(mesh, loc=NULL, field=NULL,
crs=NULL, ...)
{
inla.require.inherits(mesh, "inla.mesh", "'mesh'")
if (!is.null(mesh$crs) &&
identical(inla.as.list.CRS(mesh$crs)[["proj"]], "geocent")) {
crs.sphere <- inla.CRS("sphere")
if (!inla.identical.CRS(mesh$crs, crs.sphere)) {
## Convert the mesh to a perfect sphere.
mesh$loc <- inla.spTransform(mesh$loc, mesh$crs, crs.sphere)
mesh$manifold <- "S2"
mesh$crs <- crs.sphere
}
}
## Handle loc given as SpatialPoints or SpatialPointsDataFrame object
## or an explicit crs
if (!(missing(loc) || is.null(loc))) {
if (inherits(loc, "SpatialPoints") ||
inherits(loc, "SpatialPointsDataFrame")) {
if (is.null(mesh$crs)) {
loc <- coordinates(loc)
} else {
loc = inla.spTransform(coordinates(loc),
CRS(proj4string(loc)),
mesh$crs,
passthrough=FALSE)
}
} else if (!is.null(crs)) {
if (!is.null(mesh$crs)) {
loc = inla.spTransform(loc,
crs,
mesh$crs,
passthrough=FALSE)
}
}
}
if (!missing(field) && !is.null(field)) {
proj <- inla.mesh.projector(mesh, loc=loc, ...)
return(inla.mesh.project(proj, field=field))
}
jj =
which(rowSums(matrix(is.na(as.vector(loc)),
nrow=nrow(loc),
ncol=ncol(loc))) == 0)
smorg = (inla.fmesher.smorg(mesh$loc,
mesh$graph$tv,
points2mesh=loc[jj,,drop=FALSE]))
ti = matrix(0L, nrow(loc), 1)
b = matrix(0, nrow(loc), 3)
ti[jj,1L] = smorg$p2m.t
b[jj,] = smorg$p2m.b
ok = (ti[,1L] > 0L)
ti[ti[,1L] == 0L,1L] = NA
ii = which(ok)
A = (sparseMatrix(dims=c(nrow(loc),mesh$n),
i = rep(ii, 3),
j = as.vector(mesh$graph$tv[ti[ii,1L],]),
x = as.vector(b[ii,]) ))
list(t=ti, bary=b, A=A, ok=ok)
}
inla.mesh.project.inla.mesh.1d <- function(mesh, loc, field=NULL, ...)
{
inla.require.inherits(mesh, "inla.mesh.1d", "'mesh'")
if (!missing(field) && !is.null(field)) {
proj = inla.mesh.projector(mesh, loc, ...)
return(inla.mesh.project(proj, field))
}
A = inla.mesh.1d.A(mesh, loc=loc)
return (list(A=A,
ok=(loc >= mesh$interval[1]) & (loc <= mesh$interval[2])))
}
inla.mesh.project.inla.mesh.projector <-
function(projector, field, ...)
{
inla.require.inherits(projector, "inla.mesh.projector", "'projector'")
if (is.data.frame(field)) {
field = as.matrix(field)
}
if (is.null(dim(field))) {
if (is.null(projector$lattice)) {
data = as.vector(projector$proj$A %*% as.vector(field))
data[!projector$proj$ok] = NA
return(data)
} else {
data = as.vector(projector$proj$A %*% as.vector(field))
data[!projector$proj$ok] = NA
return(matrix(data,
projector$lattice$dims[1],
projector$lattice$dims[2]))
}
} else {
data = projector$proj$A %*% field
data[!projector$proj$ok,] = NA
return(data)
}
}
inla.mesh.projector <- function(...)
{
UseMethod("inla.mesh.projector")
}
inla.mesh.projector.inla.mesh <-
function(mesh,
loc=NULL,
lattice=NULL,
xlim=NULL,
ylim=NULL,
dims=c(100,100),
projection=NULL,
crs=NULL,
...)
{
inla.require.inherits(mesh, "inla.mesh", "'mesh'")
if (missing(loc) || is.null(loc)) {
if (missing(lattice) || is.null(lattice)) {
if (identical(mesh$manifold, "R2") &&
(is.null(mesh$crs) || is.null(crs))) {
units = "default"
lim <- list(xlim=range(mesh$loc[,1]),
ylim=range(mesh$loc[,2]))
} else if (identical(mesh$manifold, "S2") &&
(is.null(mesh$crs) || is.null(crs))) {
projection =
match.arg(projection, c("longlat", "longsinlat",
"mollweide"))
units = projection
lim = inla.mesh.map.lim(loc=mesh$loc, projection=projection)
} else {
lim <- inla.crs.bounds(crs)
if (identical(mesh$manifold, "R2")) {
lim0 <- list(xlim=range(mesh$loc[,1]),
ylim=range(mesh$loc[,2]))
lim$xlim[1] <- max(lim$xlim[1], lim0$xlim[1])
lim$xlim[2] <- min(lim$xlim[2], lim0$xlim[2])
lim$ylim[1] <- max(lim$ylim[1], lim0$ylim[1])
lim$ylim[2] <- min(lim$ylim[2], lim0$ylim[2])
}
}
if (missing(xlim) && is.null(xlim)) {
xlim = lim$xlim
}
if (missing(ylim) && is.null(ylim)) {
ylim = lim$ylim
}
x = seq(xlim[1], xlim[2], length.out=dims[1])
y = seq(ylim[1], ylim[2], length.out=dims[2])
if (is.null(mesh$crs) || is.null(crs)) {
lattice <- inla.mesh.lattice(x=x, y=y, units = units)
} else {
lattice <- inla.mesh.lattice(x=x, y=y, crs=crs)
}
} else {
dims = lattice$dims
x = lattice$x
y = lattice$y
}
if (is.null(mesh$crs) || is.null(lattice$crs)) {
proj = inla.mesh.project(mesh, lattice$loc)
} else {
proj <- inla.mesh.project(mesh,
loc=lattice$loc,
crs=lattice$crs)
}
projector = list(x=x, y=y, lattice=lattice, loc=NULL, proj=proj, crs=crs)
class(projector) = "inla.mesh.projector"
} else {
proj = inla.mesh.project(mesh, loc=loc, crs=crs)
projector = list(x=NULL, y=NULL, lattice=NULL, loc=loc, proj=proj, crs=crs)
class(projector) = "inla.mesh.projector"
}
return (projector)
}
inla.mesh.projector.inla.mesh.1d <-
function(mesh,
loc=NULL,
xlim=mesh$interval,
dims=100,
...)
{
inla.require.inherits(mesh, "inla.mesh.1d", "'mesh'")
if (missing(loc) || is.null(loc)) {
loc = seq(xlim[1], xlim[2], length.out=dims[1])
}
proj = inla.mesh.project(mesh, loc)
projector = list(x=loc, lattice=NULL, loc=loc, proj=proj)
class(projector) = "inla.mesh.projector"
return (projector)
}
inla.internal.make.spline.mesh <-
function(interval, m, degree, boundary, free.clamped)
{
boundary =
match.arg(boundary,
c("neumann", "dirichlet", "free", "cyclic"))
if (degree <= 1) {
n = (switch(boundary,
neumann = m,
dirichlet = m+2,
free = m,
cyclic = m+1))
if (n<2) {
n = 2
degree = 0
boundary = "c"
}
} else {
stopifnot(degree==2)
n = (switch(boundary,
neumann = m+1,
dirichlet = m+1,
free = m-1,
cyclic = m))
if (boundary=="free") {
if (m <= 1) {
n = 2
degree = 0
boundary = "c"
} else if (m == 2) {
n = 2
degree = 1
}
} else if (boundary=="cyclic") {
if (m <= 1) {
n = 2
degree = 0
}
}
}
return(inla.mesh.1d(seq(interval[1], interval[2], length=n),
degree=degree,
boundary=boundary,
free.clamped=free.clamped))
}
inla.mesh.basis <- function(mesh,
type="b.spline",
n=3,
degree=2,
knot.placement="uniform.area",
rot.inv=TRUE,
boundary="free",
free.clamped=TRUE,
...)
{
inla.require.inherits(mesh, c("inla.mesh", "inla.mesh.1d"), "'mesh'")
type = match.arg(type, c("b.spline", "sph.harm"))
knot.placement = (match.arg(knot.placement,
c("uniform.area",
"uniform.latitude")))
if (identical(type, "b.spline")) {
if (identical(mesh$manifold, "R1") || identical(mesh$manifold, "S1")) {
mesh1 =
inla.internal.make.spline.mesh(mesh$interval, n, degree,
boundary, free.clamped)
basis = inla.mesh.1d.A(mesh1, mesh$loc)
} else if (identical(mesh$manifold, "R2")) {
if (length(n) == 2) {
if (length(degree)==1) {
degree = rep(degree, 2)
}
if (length(boundary)==1) {
boundary = rep(boundary, 2)
}
mesh1x =
inla.internal.make.spline.mesh(range(mesh$loc[,1]),
n[1], degree[1],
boundary[1], free.clamped)
mesh1y =
inla.internal.make.spline.mesh(range(mesh$loc[,2]),
n[2], degree[2],
boundary[2], free.clamped)
basis =
inla.row.kron(inla.mesh.1d.A(mesh1y, mesh$loc[,2]),
inla.mesh.1d.A(mesh1x, mesh$loc[,1]))
} else {
warning("Old style call for an R2 mesh detected.\n Please supply a 2-element n-vector instead.")
long = ((mesh$loc[,1]-min(mesh$loc[,1]))/
diff(range(mesh$loc[,1])))*90*pi/180
sinlat = (((mesh$loc[,2]-min(mesh$loc[,2]))/
diff(range(mesh$loc[,2])))*2-1)
coslat = sapply(sinlat, function(x) sqrt(max(0.0,1.0-x^2)))
loc = (matrix(c(cos(long)*coslat,
sin(long)*coslat,
sinlat), mesh$n, 3))
knots = 0
degree = max(0L, min(n-1L, degree))
basis = (inla.fmesher.smorg(loc,
mesh$graph$tv,
bspline = c(n, degree, knots),
fem=-1)$bspline)
}
} else if (identical(mesh$manifold, "S2")) {
loc = mesh$loc
knots = identical(knot.placement, "uniform.latitude")
degree = max(0L, min(n-1L, degree))
basis = (inla.fmesher.smorg(loc,
mesh$graph$tv,
bspline = c(n, degree, knots),
fem=-1)$bspline)
if (!rot.inv) {
warning("Currently only 'rot.inv=TRUE' is supported for B-splines.")
}
} else {
stop("Only know how to make B-splines on R2 and S2.")
}
} else if (identical(type, "sph.harm")) {
if (!identical(mesh$manifold, "S2")) {
stop("Only know how to make spherical harmonics on S2.")
}
if (rot.inv) {
basis = (inla.fmesher.smorg(mesh$loc,
mesh$graph$tv,
sph0=n)$sph0)
} else {
basis = (inla.fmesher.smorg(mesh$loc,
mesh$graph$tv,
sph=n)$sph)
}
}
return(basis)
}
inla.parse.queries <-function(...)
{
queries = list(...)
if (length(queries)==0)
return(queries)
## Make sure that we have a list of names, empty or not:
if (is.null(names(queries))) {
q.names = rep("", length(queries))
} else {
q.names = names(queries)
}
## All nameless entries must be strings with query names. Replace
## empty names with those names, and set those entries to NULL.
for (query.idx in 1:length(queries)) {
if (q.names[[query.idx]]=="") {
if (is.character(queries[[query.idx]])) {
names(queries)[query.idx] = queries[[query.idx]]
queries[query.idx] = list(NULL)
} else {
queries[query.idx] = list(NULL)
warning(paste("Unnamed query ignored. Check query #",
query.idx, sep=""))
}
}
}
return(queries)
}
`inla.fmesher.smorg` <- function(loc, tv,
fem=NULL,
aniso=NULL,
gradients=FALSE,
sph0=NULL,
sph=NULL,
bspline=NULL,
points2mesh=NULL,
splitlines=NULL,
output=NULL,
keep=FALSE)
{
prefix = inla.fmesher.make.prefix(NULL, NULL)
n = nrow(loc)
s.dim = ncol(loc)
if (s.dim==1)
stop("1-D models not implemented yet.")
stopifnot(s.dim>=2, s.dim<=3)
if (missing(output) || is.null(output)) {
output.given = FALSE
output = NULL
} else {
output.given = TRUE
}
output.fem = list("c0", "g1", "g2")
output.aniso = list("g1aniso", "g2aniso")
output.gradients = list("dx", "dy", "dz")
output.sph0 = list("sph0")
output.sph = list("sph")
output.bspline = list("bspline")
output.p2m = list("p2m.t", "p2m.b")
output.splitlines <- list("split.loc", "split.idx", "split.origin",
"split.t", "split.b1", "split.b2")
## Outputs that need +1L index adjustment
indexoutput <- list("split.idx", "split.t", "split.origin")
fmesher.write(inla.affirm.double(loc), prefix, "s")
fmesher.write(inla.affirm.integer(tv)-1L, prefix, "tv")
all.args = "--smorg --input=s,tv"
## additional arguments
if (!is.null(fem)) {
all.args = paste(all.args," --fem=", max(2,fem), sep="")
if (!output.given) output = c(output, output.fem)
}
if (!is.null(aniso)) {
fmesher.write(inla.affirm.double(aniso[[1]]), prefix, "aniso.gamma")
fmesher.write(inla.affirm.double(aniso[[2]]), prefix, "aniso.vec")
all.args = paste(all.args," --aniso=aniso.gamma,aniso.vec", sep="")
if (!output.given) output = c(output, output.aniso)
}
if (!is.null(gradients) && gradients) {
all.args = paste(all.args," --grad", sep="")
if (!output.given) output = c(output, output.gradients)
}
if (!is.null(sph0)) {
all.args = paste(all.args," --sph0=", sph0, sep="")
if (!output.given) output = c(output, output.sph0)
}
if (!is.null(sph)) {
all.args = paste(all.args," --sph=", sph, sep="")
if (!output.given) output = c(output, output.sph)
}
if (!is.null(bspline)) {
all.args = (paste(all.args, " --bspline=",
bspline[1], ",", bspline[2], ",", bspline[3],
sep=""))
if (!output.given) output = c(output, output.bspline)
}
if (!is.null(points2mesh)) {
fmesher.write(inla.affirm.double(points2mesh), prefix, "p2m")
all.args = paste(all.args," --points2mesh=p2m", sep="")
if (!output.given) output = c(output, output.p2m)
}
if (!is.null(splitlines)) {
fmesher.write(inla.affirm.double(splitlines$loc), prefix, "splitloc")
fmesher.write(inla.affirm.integer(splitlines$idx)-1L, prefix, "splitidx")
all.args = paste(all.args," --splitlines=splitloc,splitidx", sep="")
if (!output.given) output = c(output, output.splitlines)
}
all.args = paste(all.args, inla.getOption("fmesher.arg"))
echoc = inla.fmesher.call(all.args=all.args, prefix=prefix)
result = list()
for (name in output) {
if (identical(name, "p2m.t"))
if (!file.exists(paste(prefix, name, sep="")))
result[[name]] = fmesher.read(prefix, "points2mesh.t")+1L
else
result[[name]] = fmesher.read(prefix, name) + 1L
else if (identical(name, "p2m.b"))
if (!file.exists(paste(prefix, name, sep="")))
result[[name]] = fmesher.read(prefix, "points2mesh.b")
else
result[[name]] = fmesher.read(prefix, name)
else if (name %in% indexoutput)
result[[name]] = fmesher.read(prefix, name) + 1L
else
result[[name]] = fmesher.read(prefix, name)
}
if (!keep)
unlink(paste(prefix, "*", sep=""), recursive=FALSE)
return (result)
}
## Deprecated: cyclic
inla.mesh.1d <- function(loc,
interval=range(loc),
boundary=NULL,
degree=1,
free.clamped=FALSE,
...)
{
## Note: do not change the order of these options without also
## changing 'basis.reduction' below.
boundary.options = c("neumann", "dirichlet", "free", "cyclic")
if ("cyclic" %in% names(list(...))) {
warning(paste("Option 'cyclic' is deprecated.",
" Use 'boundary=\"cyclic\"' instead.", sep=""))
if (cyclic) {
if (!missing(boundary)) {
boundary = match.arg.vector(boundary, boundary.options, length=2)
warning(paste("'cyclic=TRUE' overrides 'boundary=c(\"",
paste(boundary, collapse="\",\""), "\")'.", sep=""))
}
boundary = "cyclic"
}
}
boundary = match.arg.vector(boundary, boundary.options, length=2)
cyclic = !is.na(pmatch(boundary[1], "cyclic"))
if (cyclic && is.na(pmatch(boundary[2], "cyclic"))) {
stop("Inconsistent boundary specification 'boundary=c(",
paste(boundary, collapse=","), ")'.", sep="")
}
loc.orig = loc
if (cyclic)
loc =
(sort(unique(c(0, loc-interval[1]) %% diff(interval))) +
interval[1])
else
loc =
(sort(unique(c(interval,
pmax(interval[1], pmin(interval[2], loc)) ))))
n = length(loc)
if (loc[1]<interval[1])
stop("All 'loc' must be >= interval[1].")
if (loc[n]>interval[2])
stop("All 'loc' must be <= interval[2].")
if ( (degree<0) || (degree>2)) {
stop(paste("'degree' must be 0, 1, or 2. 'degree=",
degree,
"' is not supported.", sep=""))
}
if (length(free.clamped) == 1L) {
free.clamped = rep(free.clamped, 2)
}
## Number of basis functions
if (degree==0) {
basis.reduction = c(0, 1, 0, 1/2) ## neu, dir, free, cyclic
} else if (degree==1) {
basis.reduction = c(0, 1, 0, 1/2) ## neu, dir, free, cyclic
} else {
basis.reduction = c(1, 1, 0, 1)
}
m = (n+cyclic+(degree==2)*1
-basis.reduction[pmatch(boundary[1], boundary.options)]
-basis.reduction[pmatch(boundary[2], boundary.options)])
## if (m < 1+max(1,degree)) {
if (m < 1L) {
stop("Degree ", degree,
" meshes must have at least ", 1L,
" basis functions, not 'm=", m, "'.", sep="")
}
## Compute representative basis midpoints.
if ((degree==0) || (degree==1)) {
mid = loc
if (boundary[1] == "dirichlet") {
mid = mid[-1]
}
if (boundary[2] == "dirichlet") {
mid = mid[-(m+1)]
}
} else { ## degree==2
if (cyclic) {
mid = (loc + c(loc[-1], interval[2]))/2
} else {
mid = c(loc[1], (loc[-n]+loc[-1])/2, loc[n])
mid =
switch(boundary[1],
neumann = mid[-1],
dirichlet = mid[-1],
free = mid)
mid =
switch(boundary[2],
neumann = mid[-(m+1)],
dirichlet = mid[-(m+1)],
free = mid)
}
}
mesh =
list(n=n,
m=m,
loc=loc,
mid=mid,
interval=interval,
boundary=boundary,
cyclic=cyclic,
manifold = inla.ifelse(cyclic, "S1", "R1"),
degree=degree,
free.clamped=free.clamped,
idx = list(loc=NULL))
class(mesh) = "inla.mesh.1d"
if (degree<2) {
mesh$idx$loc =
inla.mesh.1d.bary(mesh, loc.orig, method="nearest")$index[,1]
} else {
if (length(mid) >= 2) {
mesh$idx$loc =
inla.mesh.1d.bary(inla.mesh.1d(mid, degree=0),
loc.orig,
method="nearest")$index[,1]
} else {
mesh$idx$loc = rep(1, length(loc.orig))
}
}
return(invisible(mesh))
}
inla.mesh.1d.bary <- function(mesh, loc, method=c("linear", "nearest"))
{
inla.require.inherits(mesh, "inla.mesh.1d", "'mesh'")
method = match.arg(method)
if (method=="linear") {
if (mesh$cyclic) {
mloc = c(mesh$loc-mesh$loc[1], diff(mesh$interval))
loc = (loc-mesh$loc[1]) %% diff(mesh$interval)
} else {
mloc = c(mesh$loc-mesh$loc[1], diff(mesh$interval))
loc = pmax(0, pmin(diff(mesh$interval), loc-mesh$loc[1]))
}
} else {
if (mesh$cyclic) {
mloc =
c(mesh$loc[mesh$n]-diff(mesh$interval),
mesh$loc,
diff(mesh$interval))
mloc = (mloc[-(mesh$n+2)] + mloc[-1])/2
loc = (loc-mloc[1]) %% diff(mesh$interval)
mloc = mloc-mloc[1]
} else {
mloc =
c(0,
(mesh$loc[1:(mesh$n-1L)] +
mesh$loc[2:mesh$n])/2-mesh$loc[1],
diff(mesh$interval))
loc = pmax(0, pmin(diff(mesh$interval), loc-mesh$loc[1]))
}
}
## Binary split method:
do.the.split <- function(knots, loc) {
n = length(knots)
if (n <= 2L) {
return(rep(1L, length(loc)))
}
split = 1L + (n-1L) %/% 2L ## Split point
upper = (loc >= knots[split])
idx = rep(0, length(loc))
idx[!upper] = do.the.split(knots[1:split], loc[!upper])
idx[upper] = split - 1L + do.the.split(knots[split:n], loc[upper])
return(idx)
}
idx = do.the.split(mloc, loc)
if (method=="nearest") {
u = rep(0, length(loc))
if (mesh$cyclic) {
found = which(idx==(mesh$n+1L))
idx[found] = 1L
}
} else { ## (method=="linear") {
u = pmax(0, pmin(1, (loc-mloc[idx])/(mloc[idx+1L]-mloc[idx]) ))
if (!mesh$cyclic) {
found = which(idx==mesh$n)
idx[found] = mesh$n-1L
u[found] = 1
}
}
if (mesh$cyclic) {
index = matrix(c(idx, (idx %% mesh$n)+1L), length(idx), 2)
bary = matrix(c(1-u, u), length(idx), 2)
} else {
index = matrix(c(idx, idx+1L), length(idx), 2)
bary = matrix(c(1-u, u), length(idx), 2)
}
return(list(index=index, bary=bary))
}
inla.mesh.1d.A <- function(mesh, loc,
weights=NULL,
derivatives=NULL,
method=c("linear", "nearest", "quadratic")
)
{
inla.require.inherits(mesh, "inla.mesh.1d", "'mesh'")
if (missing(method)) {
## Compute basis based on mesh$degree and mesh$boundary
if (mesh$degree==0) {
info = (inla.mesh.1d.A(mesh, loc, method="nearest",
weights=weights,
derivatives=
inla.ifelse(is.null(derivatives),
FALSE, TRUE)))
if (mesh$boundary[1] == "dirichlet") {
info$A = info$A[,-1,drop=FALSE]
if (!is.null(derivatives) && derivatives)
info$dA = info$dA[,-1,drop=FALSE]
}
if (mesh$boundary[2] == "dirichlet") {
info$A = info$A[,-(mesh$m+1),drop=FALSE]
if (!is.null(derivatives) && derivatives)
info$dA = info$dA[,-(mesh$m+1),drop=FALSE]
}
if (!is.null(derivatives) && derivatives)
return(info)
else
return(info$A)
} else if (mesh$degree==1) {
info = (inla.mesh.1d.A(mesh, loc, method="linear",
weights=weights,
derivatives=
inla.ifelse(is.null(derivatives),
FALSE, TRUE)))
if (mesh$boundary[1] == "dirichlet") {
info$A = info$A[,-1,drop=FALSE]
if (!is.null(derivatives) && derivatives)
info$dA = info$dA[,-1,drop=FALSE]
}
if (mesh$boundary[2] == "dirichlet") {
info$A = info$A[,-(mesh$m+1),drop=FALSE]
if (!is.null(derivatives) && derivatives)
info$dA = info$dA[,-(mesh$m+1),drop=FALSE]
}
if (!is.null(derivatives) && derivatives)
return(info)
else
return(info$A)
} else if (mesh$degree==2) {
info =
inla.mesh.1d.A(mesh, loc,
method="quadratic",
weights=weights,
derivatives=
inla.ifelse(
is.null(derivatives),
FALSE, TRUE))
adjust.matrix <- function(mesh, A) {
A = inla.as.dgTMatrix(A)
i = A@i+1L
j = A@j+1L
x = A@x
if (mesh$cyclic) {
j[j==1L] = mesh$m+1L
j[j==mesh$m+2L] = 2L
j = j-1L
} else {
if (mesh$boundary[1] == "neumann") {
j[j==1L] = 2L
j = j-1L
} else if (mesh$boundary[1] == "dirichlet") {
x[j==1L] = -x[j==1L]
j[j==1L] = 2L
j = j-1L
} else if ((mesh$boundary[1] == "free") &&
mesh$free.clamped[1]) {
j1 = which(j==1L)
i = c(i, i[j1])
j = c(j, j[j1]+1L)
x = c(x, -x[j1])
x[j1] = 2*x[j1]
}
if (mesh$boundary[2] == "neumann") {
j[j==(mesh$m+1L)] = mesh$m
} else if (mesh$boundary[2] == "dirichlet") {
x[j==(mesh$m+1L)] = -x[j==(mesh$m+1L)]
j[j==(mesh$m+1L)] = mesh$m
} else if ((mesh$boundary[2] == "free") &&
mesh$free.clamped[2]) {
j1 = which(j==(mesh$m))
i = c(i, i[j1])
j = c(j, j[j1]-1L)
x = c(x, -x[j1])
x[j1] = 2*x[j1]
}
}
return(sparseMatrix(i=i, j=j, x=x,
dims=c(length(loc), mesh$m)))
}
if (!is.null(derivatives) && derivatives) {
return(list(A=adjust.matrix(mesh, info$A),
dA=adjust.matrix(mesh, info$dA),
d2A=adjust.matrix(mesh, info$d2A)))
} else {
return(adjust.matrix(mesh, info$A))
}
} else {
stop(paste("'degree' must be 0, 1, or 2. 'degree=",
mesh$degree,
"' is not supported.", sep=""))
}
} else {
method = match.arg(method)
if (missing(weights) || is.null(weights)) {
weights = rep(1, length(loc))
}
if (!is.na(pmatch(method, c("linear", "nearest")))) {
idx = inla.mesh.1d.bary(mesh, loc, method)
dA = NULL
if (method=="linear") {
## Compute the n 1st order B-splines.
i = rep(1:length(loc), times=2)
weights.i = weights[i]
A = (sparseMatrix(i=i,
j=as.vector(idx$index),
x=weights.i*as.vector(idx$bary),
dims=c(length(loc), mesh$n)))
if (!is.null(derivatives) && derivatives) {
if (mesh$cyclic) {
d = (c(mesh$loc[-1], mesh$interval[2]) -
mesh$loc)
} else {
d = (mesh$loc[-1] - mesh$loc[-mesh$n])
}
dA = (sparseMatrix(i=i,
j=as.vector(idx$index),
x=(weights.i*c(-1/d[idx$index[,1]],
1/d[idx$index[,1]])),
dims=c(length(loc), mesh$n)))
}
} else {
## Nearest neighbours.
A = (sparseMatrix(i=1:length(loc),
j=idx$index[,1],
x=weights*idx$bary[,1],
dims=c(length(loc), mesh$n)))
}
if (missing(derivatives) || is.null(derivatives)) {
return(A)
} else {
return(list(A=A, dA=dA))
}
} else { ## Quadratic
## Compute the n+1 2nd order B-splines.
if (mesh$cyclic) {
Boundary.knots =
(c(mesh$loc, mesh$interval[2])[c(mesh$n,2)] +
diff(mesh$interval)*c(-1,1))
knots = c(mesh$loc, mesh$interval[2])
} else {
Boundary.knots =
c(mesh$loc[1] - (mesh$loc[2]-mesh$loc[1]),
mesh$loc[mesh$n] + (mesh$loc[mesh$n]-mesh$loc[mesh$n-1]))
knots = mesh$loc
}
if (FALSE) { ## Only for debugging.
## Using bs():
## Note: Intermediate step constructs dense matrix.
bsobj =
bs(x=loc,
knots=knots,
degree=2,
Boundary.knots=Boundary.knots)
bsobj =
Matrix(as.vector(bsobj[,1:(mesh$n+mesh$cyclic+1)]),
nrow(bsobj), mesh$n+mesh$cyclic+1)
}
## Direct calculation:
knots = c(Boundary.knots[1], knots)
if (mesh$cyclic) {
idx =
inla.mesh.1d.bary(inla.mesh.1d(knots, boundary="free"),
(loc-mesh$interval[1]) %%
diff(mesh$interval) + mesh$interval[1],
method="linear")
} else {
idx =
inla.mesh.1d.bary(inla.mesh.1d(knots, boundary="free"),
loc, method="linear")
}
knots = c(knots, Boundary.knots[2])
idx$index = idx$index[,1]-1L ## Indices into mesh intervals.
d = knots[2:length(knots)]-knots[1:(length(knots)-1)]
d2 = knots[3:length(knots)]-knots[1:(length(knots)-2)]
## Left intervals for each basis function:
i.l = 1:length(idx$index)
j.l = idx$index+2L
x.l = (idx$bary[,2]*d[idx$index+1]/d2[idx$index+1] * idx$bary[,2])
## Right intervals for each basis function:
i.r = 1:length(idx$index)
j.r = idx$index
x.r = (idx$bary[,1]*d[idx$index+1]/d2[idx$index] * idx$bary[,1])
## Middle intervals for each basis function:
i.m = 1:length(idx$index)
j.m = idx$index+1L
x.m = (1-(idx$bary[,1]*d[idx$index+1]/d2[idx$index]* idx$bary[,1] +
idx$bary[,2]*d[idx$index+1]/d2[idx$index+1] * idx$bary[,2]
))
i = c(i.l, i.r, i.m)
j = c(j.l, j.r, j.m)
weights.i = weights[i]
A = (sparseMatrix(i=i, j=j,
x=weights.i*c(x.l, x.r, x.m),
dims=c(length(idx$index), mesh$n+mesh$cyclic+1L)
))
dA = NULL
d2A = NULL
if (!is.null(derivatives) && derivatives) {
## dA:
## Left, right, middle intervals for each basis function:
x.l = (2 / d2[idx$index+1] * idx$bary[,2])
x.r = (-2 / d2[idx$index] * idx$bary[,1])
x.m = (-(-2/ d2[idx$index] * idx$bary[,1] +
2/ d2[idx$index+1] * idx$bary[,2]
))
dA = (sparseMatrix(i=i, j=j,
x=weights.i*c(x.l, x.r, x.m),
dims=(c(length(idx$index),
mesh$n+mesh$cyclic+1L))
))
## d2A:
## Left, right, middle intervals for each basis function:
x.l = (2 / d[idx$index+1] / d2[idx$index+1])
x.r = (2 / d[idx$index+1] / d2[idx$index])
x.m = (-(2/d[idx$index+1] / d2[idx$index] +
2/d[idx$index+1] / d2[idx$index+1] ))
d2A = (sparseMatrix(i=i, j=j,
x=weights.i*c(x.l, x.r, x.m),
dims=(c(length(idx$index),
mesh$n+mesh$cyclic+1L))
))
}
if (missing(derivatives) || is.null(derivatives)) {
return(A)
} else {
return(list(A=A, dA=dA, d2A=d2A))
}
}
}
}
inla.mesh.1d.fem <- function(mesh)
{
inla.require.inherits(mesh, "inla.mesh.1d", "'mesh'")
## Use the same matrices for degree 0 as for degree 1
if ((mesh$degree==0) || (mesh$degree==1)) {
if (mesh$cyclic) {
loc =
c(mesh$loc[mesh$n]-diff(mesh$interval),
mesh$loc,
mesh$loc[1]+diff(mesh$interval))
c0 = (loc[3:length(loc)]-loc[1:(length(loc)-2)])/2
c1.l = (loc[2:(length(loc)-1)]-loc[1:(length(loc)-2)])/6
c1.r = (loc[3:length(loc)]-loc[2:(length(loc)-1)])/6
c1.0 = (c1.l+c1.r)*2
g1.l = -1/(loc[2:(length(loc)-1)]-loc[1:(length(loc)-2)])
g1.r = -1/(loc[3:length(loc)]-loc[2:(length(loc)-1)])
g1.0 = -g1.l-g1.r
i.l = 1:mesh$n
i.r = 1:mesh$n
i.0 = 1:mesh$n
if (mesh$n > 1) {
j.l = c(mesh$n, 1:(mesh$n-1))
j.r = c(2:mesh$n, 1)
j.0 = 1:mesh$n
} else {
j.l = 1L
j.r = 1L
j.0 = 1L
}
} else {
c0 =
c((mesh$loc[2]-mesh$loc[1])/2,
(mesh$loc[mesh$n]-mesh$loc[mesh$n-1])/2
)
if (mesh$n>2) {
c0 =
c(c0[1],
(mesh$loc[3:mesh$n]-mesh$loc[1:(mesh$n-2)])/2,
c0[2]
)
}
c1.l = (mesh$loc[2:mesh$n]-mesh$loc[1:(mesh$n-1)])/6
c1.r = c1.l
c1.0 = (c(0, c1.l)+c(c1.r,0))*2
g1.l = -1/(mesh$loc[2:mesh$n]-mesh$loc[1:(mesh$n-1)])
g1.r = g1.l
g1.0 = -c(0, g1.l)-c(g1.r,0)
i.l = 2:mesh$n
i.r = 1:(mesh$n-1)
i.0 = 1:mesh$n
j.l = 1:(mesh$n-1)
j.r = 2:mesh$n
j.0 = 1:mesh$n
if (mesh$boundary[1] == "dirichlet") {
g1.0 = g1.0[-1]; g1.l = g1.l[-1]; g1.r = g1.r[-1]
c1.0 = c1.0[-1]; c1.l = c1.l[-1]; c1.r = c1.r[-1]
c0 = c0[-1]
i.l = i.l[-1]-1; i.r = i.r[-1]-1; i.0 = i.0[-1]-1
j.l = j.l[-1]-1; j.r = j.r[-1]-1; j.0 = j.0[-1]-1
} else if (mesh$boundary[1] == "free") {
g1.0[1] = 0
g1.r[1] = 0
}
if (mesh$boundary[2] == "dirichlet") {
m = mesh$m
g1.0 = g1.0[-(m+1)]; g1.l = g1.l[-m]; g1.r = g1.r[-m]
c1.0 = c1.0[-(m+1)]; c1.l = c1.l[-m]; c1.r = c1.r[-m]
c0 = c0[-(m+1)]
i.l = i.l[-m]; i.r = i.r[-m]; i.0 = i.0[-(m+1)]
j.l = j.l[-m]; j.r = j.r[-m]; j.0 = j.0[-(m+1)]
} else if (mesh$boundary[2] == "free") {
g1.0[mesh$m] = 0
g1.l[mesh$m-1] = 0
}
}
g1 =
sparseMatrix(i=c(i.l, i.r, i.0),
j=c(j.l, j.r, j.0),
x=c(g1.l, g1.r, g1.0),
dims=c(mesh$m, mesh$m))
c1 =
sparseMatrix(i=c(i.l, i.r, i.0),
j=c(j.l, j.r, j.0),
x=c(c1.l, c1.r, c1.0),
dims=c(mesh$m, mesh$m))
g2 = t(g1) %*% Diagonal(mesh$m, 1/c0) %*% g1
c0 = Diagonal(mesh$m, c0)
} else if (mesh$degree==2) {
if (mesh$cyclic) {
knots1 = mesh$loc
knots2 = c(mesh$loc[-1], mesh$interval[2])
} else {
knots1 = mesh$loc[-mesh$n]
knots2 = mesh$loc[-1]
}
knots.m = (knots1+knots2)/2
knots.d = (knots2-knots1)/2
## 3-point Gaussian quadrature
info =
inla.mesh.1d.A(mesh,
loc=(c(knots.m,
knots.m - knots.d*sqrt(3/5),
knots.m + knots.d*sqrt(3/5))),
weights =
c(knots.d*8/9, knots.d*5/9, knots.d*5/9)^0.5,
derivatives = TRUE
)
c1 = t(info$A) %*% info$A
g1 = t(info$dA) %*% info$dA
g2 = t(info$d2A) %*% info$d2A
g01 = t(info$A) %*% info$dA
g02 = t(info$A) %*% info$d2A
g12 = t(info$dA) %*% info$d2A
c0 = Diagonal(nrow(c1), rowSums(c1))
return(list(c0=c0, c1=c1, g1=g1, g2=g2, g01=g01, g02=g02, g12=g12))
} else {
stop(paste("Mesh basis degree=", mesh$degree,
" is not supported.", sep=""))
}
return(list(c0=c0, c1=c1, g1=g1, g2=g2))
}
inla.diameter <- function(x, ...) {
UseMethod("inla.diameter")
}
## Calculate upper bound for the diameter of a point set,
## by encapsulating in a circular domain.
inla.diameter.default <- function(x, manifold="", ...) {
if (nrow(x) <= 1) {
0
} else {
if (identical(manifold, "S2")) {
distance <- function(u,v) {
2*asin(pmin(1,
((u[1]-v[,1])^2 + (u[2]-v[,2])^2 + (u[3]-v[,3])^2)^0.5 / 2))
}
center <- colMeans(x)
tmp <- sqrt(sum(center^2))
if (tmp < 1e-6) {
pi
} else {
center <- center/tmp
min(pi, 2 * max(distance(center, x)))
}
} else {
distance <- function(u,v) {
d <- 0
for (k in seq_len(ncol(v))) {
d <- d + (u[k]-v[,k])^2
}
d^0.5
}
center <- rep(0, ncol(x))
for (k in seq_len(ncol(x))) {
center[k] <- mean(range(x[,k]))
}
2 * max(distance(center, x))
}
}
}
inla.diameter.inla.mesh.1d <- function(x, ...) {
diff(x$interval)
}
inla.diameter.inla.mesh <- function(x, ...) {
inla.diameter(x$loc, manifold=x$manifold, ...)
}
inla.diameter.inla.mesh.segment <- function(x, ...) {
inla.diameter(x$loc, ...)
}
inla.diameter.inla.mesh.lattice <- function(x, ...) {
inla.diameter(x$loc, ...)
}
inla.mesh.fem <- function(mesh, order=2)
{
inla.require.inherits(mesh, c("inla.mesh", "inla.mesh.1d"), "'mesh'")
if (inherits(mesh, "inla.mesh.1d")) {
if (order > 2) {
stop("order>2 not supported for 1d meshes.")
## TODO: Compute higher order matrices based on order 2 output.
}
return(inla.mesh.1d.fem(mesh))
} else {
## output name list:
if (order > 0) {
output <- paste("g", seq_len(order), sep="")
} else {
output <- NULL
}
output = c("c0", "c1", output,
"va", "ta")
return(inla.fmesher.smorg(mesh$loc, mesh$graph$tv,
fem=order, output=output))
}
}
inla.mesh.deriv <- function(mesh, loc)
{
inla.require.inherits(mesh, c("inla.mesh", "inla.mesh.1d"), "'mesh'")
if (inherits(mesh, "inla.mesh.1d")) {
stop("1d derivatives are not implemented.")
}
n.mesh = mesh$n
info = inla.mesh.project(mesh, loc=loc)
ii = which(info$ok)
n.ok = sum(info$ok)
tv = mesh$graph$tv[info$t[ii,1L],]
e1 = mesh$loc[tv[,3],]-mesh$loc[tv[,2],]
e2 = mesh$loc[tv[,1],]-mesh$loc[tv[,3],]
e3 = mesh$loc[tv[,2],]-mesh$loc[tv[,1],]
n1 = e2-e1*matrix(rowSums(e1*e2)/rowSums(e1*e1),n.ok,3)
n2 = e3-e2*matrix(rowSums(e2*e3)/rowSums(e2*e2),n.ok,3)
n3 = e1-e3*matrix(rowSums(e3*e1)/rowSums(e3*e3),n.ok,3)
g1 = n1/matrix(rowSums(n1*n1),n.ok,3)
g2 = n2/matrix(rowSums(n2*n2),n.ok,3)
g3 = n3/matrix(rowSums(n3*n3),n.ok,3)
x = cbind(g1[,1],g2[,1],g3[,1])
y = cbind(g1[,2],g2[,2],g3[,2])
z = cbind(g1[,3],g2[,3],g3[,3])
dx = (sparseMatrix(dims=c(nrow(loc),n.mesh),
i = rep(ii, 3),
j = as.vector(tv),
x = as.vector(x) ))
dy = (sparseMatrix(dims=c(nrow(loc),n.mesh),
i = rep(ii, 3),
j = as.vector(tv),
x = as.vector(y) ))
dz = (sparseMatrix(dims=c(nrow(loc),n.mesh),
i = rep(ii, 3),
j = as.vector(tv),
x = as.vector(z) ))
return (list(A=info$A, dx=dx, dy=dy, dz=dz))
}
inla.simplify.curve <- function(loc, idx, eps) {
## Variation of Ramer-Douglas-Peucker
## Uses width epsilon ellipse instead of rectangle,
## motivated by prediction ellipse for Brownian bridge
n = length(idx)
if ((n==2) || (eps==0)) {
return(idx)
}
segm = loc[idx[n],]-loc[idx[1],]
segm.len = sum(segm^2)^0.5
if (segm.len <= 1e-12) {
## End point same as start; closed curve. Split.
len2 = ((loc[idx[2:(n-1)], 1] - loc[idx[1], 1])^2 +
(loc[idx[2:(n-1)], 2] - loc[idx[1], 2])^2 )
split = which.max(len2)+1L
} else {
segm.mid = (loc[idx[n],]+loc[idx[1],])/2
segm = segm/segm.len
segm.perp = c(-segm[2], segm[1])
vec = (cbind(loc[idx[2:(n-1)], 1] - segm.mid[1],
loc[idx[2:(n-1)], 2] - segm.mid[2]))
## Always split if any point is outside the circle
epsi = min(c(eps, segm.len/2))
dist1 = abs(vec[,1]*segm[1] + vec[,2]*segm[2])/(segm.len/2)*epsi
dist2 = abs(vec[,1]*segm.perp[1] + vec[,2]*segm.perp[2])
dist = (dist1^2 + dist2^2)^0.5
## Find the furthest point, in the ellipse metric, and
## check if it inside the radius (radius=segm.len/2)
split = which.max(dist)+1L
if (dist[split-1L] < epsi) {
## Flat segment, eliminate.
return(idx[c(1, n)])
}
## Split at the furthest point.
split = which.max(dist)+1L
}
## Do the split recursively:
return(c(inla.simplify.curve(loc, idx[1L:split], eps),
inla.simplify.curve(loc, idx[split:n], eps)[-1L]
))
}
inla.contour.segment <-
function(x = seq(0, 1, length.out = nrow(z)),
y = seq(0, 1, length.out = ncol(z)),
z, nlevels = 10,
levels = pretty(range(z, na.rm=TRUE), nlevels),
groups = seq_len(length(levels)),
positive = TRUE,
eps = NULL,
crs = NULL)
{
## Input checking from contourLines:
if (missing(z)) {
if (!missing(x)) {
if (is.list(x)) {
z <- x$z
y <- x$y
x <- x$x
}
else {
z <- x
x <- seq.int(0, 1, length.out = nrow(z))
}
}
else stop("no 'z' matrix specified")
}
else if (is.list(x)) {
y <- x$y
x <- x$x
}
if (is.null(eps)) {
eps = min(c(min(diff(x)), min(diff(y))))/8
}
## End of input checking.
## Get contour pieces
curves = contourLines(x,y,z,levels=levels)
## Make a mesh for easy gradient interpolation:
latt = inla.mesh.lattice(x,y)
mesh =
inla.mesh.create(lattice=latt,
boundary=latt$segm,
extend=(list(n=3,
offset=(max(diff(range(x)),
diff(range(y)))*0.1)
))
)
## Map function values to mesh indexing:
zz = rep(0, prod(dim(z)))
zz[mesh$idx$lattice] = as.vector(z)
## Mapping from level to group value:
level2grp = function(level) {
if (length(groups)==1)
return(groups)
for (k in seq_along(groups)) {
if (levels[k] == level) {
return(groups[k])
}
}
return(0)
}
## Join all contour pieces into a single mesh.segment
## Different levels can later be identified via the grp indices.
loc = matrix(0,0,2)
idx = matrix(0,0,2)
grp = c()
for (k in seq_len(length(curves))) {
curve.loc = cbind(curves[[k]]$x, curves[[k]]$y)
curve.n = nrow(curve.loc)
## Extract the rotated gradients along the curve
curve.mid = (curve.loc[1:(curve.n-1),]+curve.loc[2:curve.n,])/2
A = inla.mesh.deriv(mesh, loc=curve.mid)
## Gradients rotated 90 degrees CW, i.e. to the direction
## of CCW curves around positive excursions:
grid.diff = cbind(A$dy %*% zz, -A$dx %*% zz)
## Determine the CCW/CW orientation
curve.diff = diff(curve.loc)
ccw = (sum(curve.diff * grid.diff) >= 0) ## True if in CCW direction
if ((ccw && positive) | (!ccw && !positive)) {
curve.idx = 1:curve.n
} else {
curve.idx = curve.n:1
}
## Filter short line segments:
curve.idx =
inla.simplify.curve(curve.loc,
curve.idx,
eps = eps)
## Reorder, making sure any unused points are removed:
curve.loc = curve.loc[curve.idx,,drop=FALSE]
curve.n = nrow(curve.loc)
curve.idx = cbind(1L:(curve.n-1L), 2L:curve.n)
## Check if the curve is closed, and adjust if it is:
if (max(abs(curve.loc[1,] - curve.loc[curve.n,])) < 1e-12) {
curve.loc = curve.loc[-curve.n,,drop=FALSE]
curve.n = nrow(curve.loc)
curve.idx = cbind(1L:curve.n, c(2L:curve.n, 1L))
}
## Add the curve:
offset = nrow(loc)
loc = rbind(loc, curve.loc)
idx = rbind(idx, curve.idx+offset)
grp = c(grp, rep(level2grp(curves[[k]]$level), curve.n-1L))
}
inla.mesh.segment(loc=loc, idx=idx, grp=grp, is.bnd=FALSE, crs=crs)
}
## Based on an idea from Elias Teixeira Krainski
## Requires splancs::nndistF
inla.nonconvex.hull.basic <-
function(points, convex=-0.15, resolution=40, eps=NULL, crs=NULL)
{
if (!(missing(points) || is.null(points)) &&
(inherits(points, "SpatialPoints") ||
inherits(points, "SpatialPointsDataFrame"))) {
points <- inla.spTransform(coordinates(points),
CRS(proj4string(points)),
crs,
passthrough=TRUE)
}
if (length(convex)==1)
convex = rep(convex,2)
if (length(resolution)==1)
resolution = rep(resolution,2)
lim = rbind(range(points[,1]), range(points[,2]))
ex = convex
if (convex[1]<0) {ex[1] = -convex[1]*diff(lim[1,])}
if (convex[2]<0) {ex[2] = -convex[2]*diff(lim[2,])}
domain = c(diff(lim[1,]), diff(lim[2,])) + 2*ex
dif = domain/(resolution-1)
if (any(dif > min(convex))) {
req.res = ceiling(domain/convex+1)
warning(paste("Resolution (",
paste(resolution,collapse=","),
") too small for convex (",
paste(convex,collapse=","),
").\n",
"Resolution >=(",
paste(req.res,collapse=","),
") required for more accurate results.",
sep=""))
}
ax = list(seq(lim[1,1] - ex[1], lim[1,2] + ex[1], length=resolution[1]),
seq(lim[2,1] - ex[2], lim[2,2] + ex[2], length=resolution[2]))
xy = as.matrix(expand.grid(ax[[1]], ax[[2]]))
tr = diag(c(1/ex[1],1/ex[2]))
stopifnot(inla.require("splancs"))
z = (matrix(splancs::nndistF(points%*%tr, xy%*%tr),
resolution[1], resolution[2]))
segm = inla.contour.segment(ax[[1]], ax[[2]], z,
levels=c(1), positive=FALSE, eps=eps, crs=crs)
return(segm)
}
## Morphological dilation by "convex",
## followed by closing by "concave", with
## minimum concave curvature radius "concave".
## If the dilated set has no gaps of width between
## 2*convex*(sqrt(1+2*concave/convex) - 1)
## and
## 2*concave,
## then the minimum convex curvature radius is "convex".
## Default is concave=convex
## Special case concave=0 delegates to inla.nonconvex.hull.basic()
##
## The implementation is based on the identity
## dilation(a) & closing(b) = dilation(a+b) & erosion(b)
## where all operations are with respect to disks with the specified radii.
inla.nonconvex.hull <-
function(points, convex=-0.15, concave=convex, resolution=40, eps=NULL,
crs=NULL)
{
if (!(missing(points) || is.null(points)) &&
(inherits(points, "SpatialPoints") ||
inherits(points, "SpatialPointsDataFrame"))) {
points <- inla.spTransform(coordinates(points),
CRS(proj4string(points)),
crs,
passthrough=TRUE)
}
if (length(resolution)==1)
resolution = rep(resolution,2)
lim = rbind(range(points[,1]), range(points[,2]))
approx.diam = max(diff(lim[1,]),diff(lim[2,]))
if (convex<0) {convex = -convex*approx.diam}
if (concave<0) {concave = -concave*approx.diam}
if (concave==0) {
return(inla.nonconvex.hull.basic(points, convex, resolution, eps,
crs=crs))
}
ex = convex+concave
domain = c(diff(lim[1,]), diff(lim[2,])) + 2*ex
dif = domain/(resolution-1)
if (max(dif) > min(convex,concave)) {
req.res = ceiling(domain/min(convex,concave)+1)
warning(paste("Resolution (",
paste(resolution,collapse=","),
") too small for convex/concave radius (",
convex,",",concave,
").\n",
"Resolution >=(",
paste(req.res,collapse=","),
") required for more accurate results.",
sep=""))
}
ax =
list(
seq(lim[1,1] - ex, lim[1,2] + ex, length=resolution[1]),
seq(lim[2,1] - ex, lim[2,2] + ex, length=resolution[2])
)
xy = as.matrix(expand.grid(ax[[1]], ax[[2]]))
stopifnot(inla.require("splancs"))
z = (matrix(splancs::nndistF(points, xy),
resolution[1],resolution[2]))
segm.dilation =
inla.contour.segment(ax[[1]], ax[[2]], z,
levels=c(convex+concave),
positive=TRUE,
eps=0) ## Don't simplify curve at this stage
mesh.dilation =
inla.mesh.create(loc=xy,
boundary=segm.dilation,
extend=(list(n=3,
offset=(max(diff(ax[[1]]),
diff(ax[[2]]))*0.1)
))
)
## This filtering is not necessary; the inla.mesh.create() should
## have removed all unused points. 2013-02-24 /FL
## points.dilation =
## mesh.dilation$loc[unique(as.vector(mesh.dilation$graph$tv)),]
z = (matrix(splancs::nndistF(mesh.dilation$loc, xy),
resolution[1],resolution[2]))
segm.closing =
inla.contour.segment(ax[[1]],ax[[2]],z,
levels=c(concave),
positive=TRUE,
eps=eps)
segm.closing$crs <- crs
segm.closing
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.