-
Notifications
You must be signed in to change notification settings - Fork 1
/
split_lines.R
115 lines (104 loc) · 3.11 KB
/
split_lines.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
#' @include deprecated.R
# fm_split_lines ####
#' @title Split lines at triangle edges
#'
#' @description Compute intersections between line segments and triangle edges,
#' and filter out segment of length zero.
#'
#' @param mesh An `fm_mesh_2d` or `inla.mesh` object
#' @param segm An [fm_segm()] object with segments to be split
#' @param ... Unused.
#' @return An [fm_segm()] object with the same crs as the mesh,
#' with an added field `origin`, that for each new segment gives the
#' originator index into to original `segm` object for each new line segment.
#'
#' @author Finn Lindgren \email{finn.lindgren@@gmail.com}
#'
#' @export
#' @examples
#' mesh <- fm_mesh_2d(
#' boundary = fm_segm(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 1)), is.bnd = TRUE)
#' )
#' splitter <- fm_segm(rbind(c(0.8, 0.2), c(0.2, 0.8)))
#' segm_split <- fm_split_lines(mesh, splitter)
#'
#' plot(mesh)
#' lines(splitter)
#' points(segm_split$loc)
fm_split_lines <- function(mesh, ...) {
UseMethod("fm_split_lines")
}
#' @rdname fm_split_lines
#' @export
fm_split_lines.fm_mesh_2d <- function(mesh, segm, ...) {
segm <-
fm_transform(
fm_as_segm(segm),
crs = fm_crs(mesh),
crs0 = fm_crs(segm),
passthrough = TRUE
)
origin <- seq_len(NROW(segm$idx))
if (NROW(segm$loc) > 0) {
# Filter out segments not on the mesh
t1 <- fm_bary(mesh, loc = segm$loc, crs = fm_crs(segm))$t
keep <- !(is.na(t1[segm$idx[, 1]]) | is.na(t1[segm$idx[, 2]]))
# if (any(!keep)) { warning("points outside boundary! filtering...")}
segm <- fm_segm(
loc = segm$loc,
idx = segm$idx[keep, , drop = FALSE],
grp = segm$grp[keep],
is.bnd = all(segm$is.bnd),
crs = fm_crs(segm)
)
origin <- origin[keep]
}
if (NROW(segm$idx) == 0) {
segm$origin <- origin
return(segm)
}
# Split the segments into parts
if (NCOL(segm$loc) == 2) {
segm$loc <- cbind(segm$loc, rep(0, NROW(segm$loc)))
}
splt <- fmesher_split_lines(
mesh_loc = mesh$loc,
mesh_tv = mesh$graph$tv - 1L,
loc = segm$loc,
idx = segm$idx - 1L,
options = list()
)
indexoutput <- list("split.idx", "split.t", "split.origin")
for (name in intersect(names(splt), indexoutput)) {
splt[[name]] <- splt[[name]] + 1L
}
segm.split <- fm_segm(
loc = splt$split.loc,
idx = splt$split.idx,
grp = segm$grp[splt$split.origin],
is.bnd = all(segm$is.bnd),
crs = fm_crs(segm)
)
origin <- origin[splt$split.origin]
# plot(mesh)
# lines(segm)
# points(segm.split$loc[, 1:2], col="red", pch = 20)
# points(segm$loc[, 1:2], col="blue", pch = 20)
# Filter out zero length segments
keep <- rowSums((segm.split$loc[segm.split$idx[, 2], , drop = FALSE] -
segm.split$loc[segm.split$idx[, 1], , drop = FALSE])^2) > 0
segm.split <- fm_segm(
loc = segm.split$loc,
idx = segm.split$idx[keep, , drop = FALSE],
grp = segm.split$grp[keep],
is.bnd = all(segm.split$is.bnd),
crs = fm_crs(segm)
)
segm.split$origin <- origin[keep]
return(segm.split)
}
#' @rdname fm_split_lines
#' @export
fm_split_lines.inla.mesh <- function(mesh, ...) {
fm_split_lines(fm_as_mesh_2d(mesh), ...)
}