-
Notifications
You must be signed in to change notification settings - Fork 4
/
flywire-coord.R
165 lines (151 loc) · 6.14 KB
/
flywire-coord.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
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
#' Map points between FlyWire v1 and FAFB14 space (xyz nm)
#'
#' @description \code{flywire2fafb} maps points FlyWire->FAFB
#'
#' @details Note that you can also access FlyWire<->FAFB bridging registrations
#' via the \code{\link{xform_brain}} series of functions. This will allow you
#' to transform most kinds of 3D data objects, whereas the \code{flywire2fafb}
#' function is restricted to plain 3D coordinates. See examples.
#'
#' Mapping single points is unlikely to be useful, but you may wish to adjust
#' the \code{chunksize} argument to send more points at once at the risk of
#' possible server timeouts. The default value is quite conservative.
#'
#' When \code{swap=TRUE} displacements will be applied in the opposite
#' direction to what is intended. This can be used to provide a coarse inverse
#' mapping if you feed in FAFB points. This is wrong but was useful before the
#' inverse mapping was available as it can get you closer to the right place
#' in FlyWire than just assuming that FAFB14 and FlyWire are in the same
#' space. This works because deformations are mostly fairly smooth at the
#' scale of FAFB-FlyWire displacements. Operationally we find that residual
#' displacements are typically of the order 100 nm using this procedure. Since
#' 2020-08-08 we have a real inverse available, so this is now only of
#' historical interest.
#'
#' @param xyz A Nx3 matrix of points
#' @param method Whether to map many points at once (default) or just one
#' @param chunksize The number of points to send to the server when mapping many
#' points at once.
#' @param swap When \code{TRUE} applies the deformation field in the opposite
#' direction e.g. to give a coarse mapping of points FAFB->FlyWire. This is
#' wrong but may be useful.
#' @param ... Additional arguments passed to \code{httr::GET}/\code{POST}
#' operation
#'
#' @return an Nx3 matrix of points
#' @export
#'
#' @examples
#' # identified location in FAFB14
#' p.fafb.nm <- cbind(477042, 284535, 90680)
#' p.fafb.raw <- p.fafb.nm/c(4,4,40)
#' # corresponding location in FlyWire
#' p.flywire.raw <- cbind(118865, 71338, 2267)
#' p.flywire.nm <- p.flywire.raw * c(4,4,40)
#'
#' # check displacement
#' flywire2fafb(p.flywire.nm)-p.fafb.nm
#'
#' # check what happens when you apply the inverse
#' fafb2flywire(p.fafb.nm)-p.flywire.nm
#'
#' data("AV4b1", package='catmaid')
#' set.seed(42)
#' before=xyzmatrix(AV4b1)[sample(nvertices(AV4b1), size=2000), ]
#' after=fafb2flywire(before)
#' d=sqrt(rowSums((before-after)^2))
#' hist(d, br=20, main="FAFB14 - FlyWire Displacements", xlab="d /nm")
#'
#' \dontrun{
#' AV4b1.flywire <- xform_brain(AV4b1, reference="FlyWire", sample="FAFB14")
#' plot3d(neuronlist(AV4b1.flywire, AV4b1))
#' }
flywire2fafb <- function(xyz, method=c("mapmany", "map1"), chunksize=40e3,
swap=FALSE, ...) {
method=match.arg(method)
if(swap)
warn_hourly("Please use fafb2flywire for more accurate FAFB->FlyWire transforms. See ?flywire2fafb")
baseurl <- "https://services.itanna.io/app/transform-service/dataset/flywire_v1"
mapwrapper(xyz, baseurl=baseurl, method=method, chunksize=chunksize, swap=swap, ...)
}
#' @rdname flywire2fafb
#' @export
#' @description \code{fafb2flywire} maps points FAFB->FlyWire
fafb2flywire <- function(xyz, method=c("mapmany", "map1"), chunksize=40e3, swap=FALSE, ...) {
method=match.arg(method)
if(swap)
warn_hourly("Please use flywire2fafb for more accurate FlyWire->FAFB transforms. See ?flywire2fafb")
baseurl <- "https://services.itanna.io/app/transform-service/dataset/flywire_v1_inverse"
mapwrapper(xyz, baseurl=baseurl, method=method, chunksize=chunksize, swap=swap, scale=4, ...)
}
# Private function to make bridging registration available to xform and friends
register_fafb_flywire <- function() {
flywire2fafb.reg <- nat::reglist(function(xyz, ...) flywire2fafb(xyz, ...))
fafb2flywire.reg <- nat::reglist(function(xyz, ...) fafb2flywire(xyz, ...))
nat.templatebrains::add_reglist(flywire2fafb.reg, sample = 'FlyWire',
reference='FAFB14')
nat.templatebrains::add_reglist(fafb2flywire.reg, reference = 'FlyWire',
sample='FAFB14')
}
#' Handle raw and nm calibrated flywire coordinates
#'
#' @description \code{flywire_voxdims} returns the image voxel dimensions which
#' are normally used to scale between \bold{raw} and \bold{nm} coordinates.
#'
#' @param url Optional neuroglancer URL containing voxel size. Defaults to
#' \code{getOption("fafbseg.sampleurl")} as set by
#' \code{\link{choose_segmentation}}.
#'
#' @return For \code{flywire_voxdims} A 3-vector
#' @export
#'
#' @examples
#' flywire_voxdims()
#' # ensure that we use default production flywire scene
#' with_segmentation('flywire', flywire_voxdims())
flywire_voxdims <- memoise::memoise(function(url=getOption("fafbseg.sampleurl")) {
sc=ngl_blank_scene(url)
voxdims(sc)
})
#' @param x 3D coordinates in any form compatible with \code{\link{xyzmatrix}}
#'
#' @return for \code{flywire_raw2nm} and \code{flywire_nm2raw} an Nx3 matrix of
#' coordinates
#' @param vd The voxel dimensions in nm. Expert use only. Normally found
#' automatically.
#' @export
#' @rdname flywire_voxdims
#' @details relies on nat >= 1.10.4
#' @examples
#' flywire_raw2nm(c(159144, 22192, 3560))
#' flywire_raw2nm('159144 22192 3560')
#' \dontrun{
#' flywire_nm2raw(clipr::read_clip())
#' }
flywire_nm2raw <- function(x, vd=flywire_voxdims()) {
xyz=xyzmatrix(x)
xyz[,1]=xyz[,1]/vd[1]
xyz[,2]=xyz[,2]/vd[2]
xyz[,3]=xyz[,3]/vd[3]
xyz
}
#' @export
#' @rdname flywire_voxdims
flywire_raw2nm <- function(x, vd=flywire_voxdims()) {
xyz=xyzmatrix(x)
xyz[,1]=xyz[,1]*vd[1]
xyz[,2]=xyz[,2]*vd[2]
xyz[,3]=xyz[,3]*vd[3]
xyz
}
# FIXME try to make a generic version of this for different stacks?
is_rawcoord <- function(xyz) {
# dput(boundingbox(elmr::FAFB14)/c(4,4,40))
rawbb=makeboundingbox(c(0, 253951, 0, 155647, 0, 7062))
pointsinside(xyz, rawbb)
}
.spine_baseurl <- "https://services.itanna.io"
spine_ok <- memoise::memoise(~memoise::timeout(10*60), f=function() {
status=try(httr::status_code(httr::HEAD(.spine_baseurl, httr::timeout(2))), silent = T)
identical(status, 200L)
})