-
Notifications
You must be signed in to change notification settings - Fork 4
/
read_merge_info.R
143 lines (138 loc) · 4.82 KB
/
read_merge_info.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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#' Read CSV files containing raw agglomeration (merge) information
#'
#' @details Peter Li provided 10 CSV files with merge information in Aug 2018
#' against which this function has been tested. This merge information has
#' been processed to identify which groups segments should be merged together
#' in the \code{\link[fafbsegdata]{fafbsegdata}} package, which can be
#' accessed using the \code{\link{find_merged_segments}} function.
#'
#' @param x Input directory
#' @param pattern Optional pattern to identify CSV files
#' @param voxdims Optional voxel dimensions to convert raw pixel coordinates
#' into physical (nm) coordinates
#' @param ... Additional arguments passed to \code{readr::\link{read_csv}}
#'
#' @return a \code{data.frame} with columns \itemize{
#'
#' \item \code{id1,id2} Segment ids to be merged
#'
#' \item \code{x,y,z} Location (in nm) of merge point
#'
#' }
#' @export
#' @seealso \code{\link{find_merged_segments}},
#' \code{\link[fafbsegdata]{fafbsegdata}} package
#'
#' @examples
#' \dontrun{
#' csvdir="~/projects/fafbseg/fafb14_v00c_split3xfill2x_merges/ffnreseg16nm_ms1000_md0.02_c0.6_iou0.7/"
#' mergeinfo <- read_mergeinfo(csvdir)
#'
#' # Show 3D position of merge locations for a large merge group containing a
#' # tangential cell (widefield visual interneuron)
#' tang_segs <- find_merged_segments(22139217567)
#' tang_mi <- subset(mergeinfo, id1%in% tang_segs | id2 %in%tang_segs)
#' points3d(xyzmatrix(tang_mi))
#'
#' # plot brain surface for context
#' library(elmr)
#' plot3d(FAFB)
#' }
#' @importFrom readr read_csv col_double col_integer cols
read_mergeinfo <- function(x, pattern = "10\\.csv$", voxdims = c(16, 16, 40), ...) {
read_merge_csv <- function(x, voxdims) {
res = readr::read_csv(
x,
col_names = c("id1", 'id2', 'x', 'y', 'z'),
col_types = readr::cols(
`id1` = readr::col_double(),
`id2` = readr::col_double(),
`x` = readr::col_integer(),
`y` = readr::col_integer(),
`z` = readr::col_integer()
)
)
res$x = res$x * voxdims[1]
res$y = res$y * voxdims[2]
res$z = res$z * voxdims[3]
res
}
ff = dir(x, pattern = pattern, full.names = TRUE)
stopifnot(length(ff)>0)
mergel = lapply(ff, read_merge_csv, voxdims = voxdims)
dplyr::bind_rows(mergel)
}
#' Find all merged segment ids for a given starting google segment id
#'
#' This is a naive implementation that only depends only on the \code{mergeinfo}
#' object created by \code{\link{read_mergeinfo}}. It is retained only to verify
#' correctness. It is slow but does not require any pre-processing of the table
#' of merge information.
#'
#' @param x a segment id or any other input that can be interpreted by
#' \code{\link{ngl_segments}}
#' @param mergeinfo The merge information data.frame created by
#' \code{\link{read_mergeinfo}}
#' @return vector of segment ids
#' @seealso \code{\link{find_merged_segments}}
find_merged_segments_slow <- function(x, mergeinfo) {
x=ngl_segments(x)
while(TRUE) {
cat('.')
w1=which(mergeinfo$id1%in%x)
w2=which(mergeinfo$id2%in%x)
w=union(w1,w2)
if(length(w) < 1) break
segs=union(mergeinfo$id1[w], mergeinfo$id2[w])
newsegs=setdiff(segs, x)
if(length(newsegs)==0) break
x=union(x, newsegs)
}
x
}
#' Graph representation of the agglomeration (merge) between raw segments
#'
#' @description \code{make_merge_graph} generates a graph representation of all
#' agglomeration merges.
#'
#' @param x The \code{mergeinfo} \code{data.frame} generated by
#' \code{\link{read_mergeinfo}}
#' @param g Either a graph generated by \code{make_merge_graph} or the
#' \code{mergeinfo} \code{data.frame} generated by
#' \code{\link{read_mergeinfo}}.
#' @param n For testing purposes, the number of rows of merge information to
#' process - the default processes all rows.
#'
#' @return \code{make_merge_graph} returns an \code{igraph} object,
#' \code{merge_graph_components} returns a list with elements describing the
#' mapping between segment ids and merge groups.
#' @export
#'
#' @examples
#' \dontrun{
#' g=make_merge_graph(mergeinfo)
#' gc=merge_graph_components(g)
#' str(gc)
#' }
make_merge_graph <- function(x, n=NA) {
el=x[,1:2]
if(!is.na(n))
el=el[seq_len(n),,drop=FALSE]
el=data.matrix(el)
vertexnames=sort(unique(c(el[,1], el[,2])))
rawel=match(t(el), vertexnames)
if(any(is.na(rawel))) stop("Bad vertices!")
g=igraph::graph(rawel, n=length(vertexnames), directed=FALSE)
igraph::V(g)$id=vertexnames
g
}
#' @rdname make_merge_graph
#' @description \code{merge_graph_components} uses the full merge graph to find
#' each connected subgraph corresponding to isolated merge groups.
merge_graph_components <- function(g, n=NA) {
if(!igraph::is.igraph(g)) g=make_merge_graph(g, n=n)
gc=igraph::components(g)
gc$id=igraph::vertex_attr(g, "id")
gc$membership=as.integer(gc$membership)
gc
}