forked from artxz/FAFB-photoreceptor-connectivity
-
Notifications
You must be signed in to change notification settings - Fork 0
/
layer_DRA_anterior.R
276 lines (231 loc) · 11.6 KB
/
layer_DRA_anterior.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
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
# DRA use Mi1 to define a coord trans -----------------------------------------
# - 2 col of DRA R7/8
DRAR7_me <- nlapply(DRAR7, function(x) subset(x, pointsinside(x, ME_msh))) #me portion
DRAR8_me <- nlapply(DRAR8, function(x) subset(x, pointsinside(x, ME_msh)))
# This is the difference, use the 3rd (anterior) R78 pair instead of the 1st (polar)
DRA_ref_com <- rbind(xyzmatrix(DRAR7_me[[3]]), xyzmatrix(DRAR8_me[[3]])) %>% colMeans() #change back to polar col
ii <- sweep(xyz_M5_avg, 2, DRA_ref_com)^2 %>% rowSums() %>% sqrt() %>% order() %>% head(.,7)
DRA_Mi1 <- nlapply(Mi1[ii], function(x) subset(x, pointsinside(x,ME_msh,rval='distance') > 0))
# pca
node_xyz <- xyzmatrix(DRA_Mi1)
DRA_me_pca <- prcomp(node_xyz)
if (DRA_me_pca$rotation[,1] %*% c(-0.84, 0.20, -0.49) < 0) {
DRA_me_pca$rotation <- - DRA_me_pca$rotation
}
if (t(cross3D(DRA_me_pca$rotation[,1],DRA_me_pca$rotation[,2])) %*% DRA_me_pca$rotation[,3] < 0 ) {
DRA_me_pca$rotation[,3] <- - DRA_me_pca$rotation[,3]
}
DRAR7_xform <- xEucl_neu(DRAR7, DRA_me_pca$rotation, DRA_me_pca$center)
DRAR8_xform <- xEucl_neu(DRAR8, DRA_me_pca$rotation, DRA_me_pca$center)
# DRA Mi1
anno_Mi1_DRA <- catmaid_query_by_annotation('DRA_column') # DRA col, light blue, or darker gray
anno_Mi1_DRA_excl <- catmaid_query_by_annotation('non_PR_column') # DRA col excluded
ind_Mi1_DRA <- which(anno_Mi1$skid %in% anno_Mi1_DRA$skid)
ind_Mi1_DRA_hcol <- c(133, 268, 42)
ind_Mi1_DRA <- ind_Mi1_DRA[!(ind_Mi1_DRA %in% ind_Mi1_DRA_hcol)] #exclude home column
ind_Mi1_DRA_excl <- which(anno_Mi1$skid %in% anno_Mi1_DRA_excl$skid)
# - meshes
# - transfomr medulla mesh
DRA_ME_msh_xform <- ME_msh
DRA_ME_msh_xform$vb[1:3,] <- sweep(t(ME_msh$vb[1:3,]), 2, DRA_me_pca$center) %*% DRA_me_pca$rotation %>% t()
DRA_AOTU_msh_xform <- AOTU_msh
DRA_AOTU_msh_xform$vb[1:3,] <- sweep(t(AOTU_msh$vb[1:3,]), 2, DRA_me_pca$center) %*% DRA_me_pca$rotation %>% t()
DRA_ME_L_msh_xform <- ME_L_msh
DRA_ME_L_msh_xform$vb[1:3,] <- sweep(t(ME_L_msh$vb[1:3,]), 2, DRA_me_pca$center) %*% DRA_me_pca$rotation %>% t()
DRA_PLP_R_msh_xform <- PLP_R_msh
DRA_PLP_R_msh_xform$vb[1:3,] <- sweep(t(PLP_R_msh$vb[1:3,]), 2, DRA_me_pca$center) %*% DRA_me_pca$rotation %>% t()
# - axes
DRA_axis_ori <- c(-110000, -50000, -40000)
DRA_axis_lat <- -DRA_me_pca$rotation[1,] * 15000
DRA_axis_dor <- DRA_me_pca$rotation[2,] * -15000
DRA_axis_post <- DRA_me_pca$rotation[3,] * 15000
# - scale bar
DRA_xAng <- -30 #rotation angle around x-axis, align with long axis of ME
# DRA_xAng <- 20 # align col
DRA_xAngRot <- (DRA_xAng - 90)/180*pi
DRA_xRot <- matrix(c(1,0,0,
0, cos(DRA_xAngRot), sin(DRA_xAngRot),
0, -sin(DRA_xAngRot), cos(DRA_xAngRot)), ncol = 3, byrow = T)
DRA_xAngRot2 <- DRA_xAng/180*pi
DRA_xRot2 <- matrix(c(1,0,0,
0, cos(DRA_xAngRot2), sin(DRA_xAngRot2),
0, -sin(DRA_xAngRot2), cos(DRA_xAngRot2)), ncol = 3, byrow = T)
# top view
DRA_xAng_top <- 30
DRA_xAngRot_top <- (DRA_xAng_top)/180*pi
DRA_xRot_top <- matrix(c(1,0,0,
0, cos(DRA_xAngRot_top), sin(DRA_xAngRot_top),
0, -sin(DRA_xAngRot_top), cos(DRA_xAngRot_top)), ncol = 3, byrow = T)
DRA_yAng_top <- 90
DRA_yAngRot_top <- DRA_yAng_top/180*pi
DRA_yRot_top <- matrix(c(1,0,0,
0, cos(DRA_yAngRot_top), sin(DRA_yAngRot_top),
0, -sin(DRA_yAngRot_top), cos(DRA_yAngRot_top)), ncol = 3, byrow = T)
DRA_zAng_top <- 80
DRA_zAngRot_top <- (DRA_zAng_top)/180*pi
DRA_zRot_top <- matrix(c(1,0,0,
0, cos(DRA_zAngRot_top), sin(DRA_zAngRot_top),
0, -sin(DRA_zAngRot_top), cos(DRA_zAngRot_top)), ncol = 3, byrow = T)
DRA_sbar <- matrix(c(30000,-5000,40000, 30000,5000,40000), ncol=3,byrow=T) %*%
t(matrix(c(1,0,0,
0,0,1,
0,-1,0), ncol=3, byrow=T))
# -- side, anterior view
DRA_sbar_side <- matrix(c(30000,-5000,40000, 30000,5000,40000), ncol=3,byrow=T) %*% DRA_xRot_top
# #rotation angle around y-axis to level layer 5
# DRA_yAng <- 0
# DRA_yAngRot <- -DRA_yAng/180*pi
# DRA_yRot <- matrix(c(cos(DRA_yAngRot), 0, sin(DRA_yAngRot),
# 0,1,0,
# -sin(DRA_yAngRot), 0, cos(DRA_yAngRot)), ncol = 3, byrow = T)
#
# # DRA_sbar_rot90 <- DRA_sbar %*% t(matrix(c(1,0,0,
# # 0,0,1,
# # 0,-1,0), ncol=3, byrow=T)) %*% t(DRA_xRot) %*% t(DRA_yRot)
#
# DRA_sbar_rot90 <- matrix(c(30000,5000,10000, 30000,5000,20000), ncol=3,byrow=T)
# col for top view --------------------------------------------------------
aMe12_ind_Mi1 <- list()
Mi1_ind_aMe12 <- list()
for (ii_n in 1:3) {
tar <- aMe12[[ii_n]]
ii_col = match(tar$tags$`columnar branch`, tar$d$PointNo)
xyz_col = xyzmatrix(tar$d[ii_col,])
aMe12_ind_Mi1[[ii_n]] <- apply(xyz_col, 1, function(x) order(rowSums(sweep(Mi1_M5_xyz,2,x)^2))[1] )
Mi1_ind_aMe12[[ii_n]] <- apply(Mi1_M5_xyz, 1, function(x) order(rowSums(sweep(xyz_col,2,x)^2))[1] )
}
i1 <- sqrt(rowSums(sweep(xyz_M5_avg, 2, xyz_M5_avg[ind_Mi1_DRA_hcol[1],], '-')^2)) <60000
i2 <- sqrt(rowSums(sweep(xyz_M5_avg, 2, xyz_M5_avg[ind_Mi1_DRA_hcol[2],], '-')^2)) <30000
i3 <- sqrt(rowSums(sweep(xyz_M5_avg, 2, xyz_M5_avg[ind_Mi1_DRA_hcol[3],], '-')^2)) <30000
ind_Mi1_top <- which(i1|i2|i3)
ind_Mi1_top_grey <- ind_Mi1_top[!(ind_Mi1_top %in% c(ind_Mi1_DRA_excl, ind_Mi1_DRA))]
ind_Mi1_top_red <- ind_Mi1_top[ind_Mi1_top %in% c(ind_Mi1_DRA)]
ind_Mi1_top_y <- ind_Mi1_top_grey[!(ind_Mi1_top_grey %in% unlist(aMe12_ind_Mi1))]
ind_Mi1_top_p <- ind_Mi1_top_grey[(ind_Mi1_top_grey %in% unlist(aMe12_ind_Mi1))]
# meshes ------------------------------------------------------------------
# - Mi1 M5 nodes
xyz_M5_avg_xform_DRA <- sweep(xyz_M5_avg, 2, DRA_me_pca$center) %*% DRA_me_pca$rotation
xyz_M5_avg_xform_DRA_yz <- xyz_M5_avg_xform_DRA
xyz_M5_avg_xform_DRA_yz[,1] <- 0
Nnb <- 16
ind_nb <- sweep(xyz_M5_avg_xform_DRA, 2,
xyz_M5_avg_xform_DRA[ind_Mi1_DRA_hcol[3],])^2 %>%
rowSums() %>%
order() %>%
.[1:Nnb]
nopen3d()
points3d(xyz_M5_avg_xform_DRA_yz, col='grey', size = 20)
# points3d(xyz_M5_avg_xform_DRA, col='grey', size = 10)
# points3d(ref_com4,col='red',size=10)
points3d(xyz_M5_avg_xform_DRA_yz[ind_nb,], size=25, col='blue')
plot3d(DRAR7_xform, col='green', lwd=2)
plot3d(DRAR8_xform, col='gold2', lwd=2)
# -- make local mesh
xyz <- t(DRA_ME_msh_xform$vb[1:3,])
dd <- sweep(xyz[,2:3], 2, colMeans(xyz_M5_avg_xform_DRA_yz[ind_nb,])[-1]) %>% .^2 %>% rowSums() %>% sqrt()
xyz_msh <- xyz[dd < 1.2*max(dist(xyz_M5_avg_xform_DRA_yz[ind_nb,])), ]
DRA_ME_msh_local_3 <- ashape3d(xyz_msh, alpha = 60000) %>% as.mesh3d()
# - transfomr LO mesh
DRA_LO_msh_xform <- LO_msh
DRA_LO_msh_xform$vb[1:3,] <- sweep(t(LO_msh$vb[1:3,]), 2, DRA_me_pca$center) %*% DRA_me_pca$rotation %>% t()
# ME layer side DRAR7, 3 - anterior ---------------------------------------
# --layers from 7-column
# bp7c <- c(0, 5.5,14.5,20,23,28,36,42,48,56,61) #layer 1 to 10
bp7c <- c(0, 5.5, 14.5, 20, 23, 28, 33, 39, 45, 56, 61) #layer 1 to 10, modify 6-7-8-9 by Mi1, C2, Tm5 and Tm20
bp7c <- max(bp7c) - bp7c %>% rev() #layer 10 to 1
bp7c_prob <- bp7c/max(bp7c)
# -- ME boundary
yz_msh <- t(DRA_ME_msh_local_3$vb)
# yz_dorsal <- data.frame(y = seq(80000, 35000, length.out = 20),
# z = seq(-40000, -0000, length.out = 20))
#
# yz_posterior <- data.frame(y = seq(30000, 70000, length.out = 16),
# z = seq(-40000, -8000, length.out = 16))
yz_posterior <- data.frame(y = seq(-40000, 0000, length.out = 20),
z = seq(20000, -20000, length.out = 20))
yz_dorsal <- data.frame(y = seq(-25000, 25000, length.out = 16),
z = seq(-25000, 25000, length.out = 16))
# choose yz
nopen3d()
points3d(cbind(0, yz_dorsal))
points3d(yz_posterior, col='red')
# points3d(yz_msh[,1:3])
shade3d(DRA_ME_msh_local_3, alpha=0.1, col='gray')
shade3d(DRA_ME_msh_xform, alpha=0.1, col='gold')
axes3d(c('x','y','z')); title3d('','','x','y','z')
rgl.viewpoint(fov=0)
# points3d(xyz_bot)
# points3d(top_bd_dorsal_3)
# points3d(bot_bd_dorsal_3)
# points on mesh
xyz <- t(DRA_ME_msh_xform$vb[1:3,])
dd <- sweep(xyz[,2:3], 2, colMeans(xyz_M5_avg_xform_DRA_yz[ind_nb,])[-1]) %>% .^2 %>% rowSums() %>% sqrt()
xyz_msh <- xyz[dd < max(dist(xyz_M5_avg_xform_DRA_yz[ind_nb,])), ]
# top
xyz_top <- xyz_msh[xyz_msh[,1] > 0, ]
x <- xyz_top[,1]; y <- xyz_top[,2]; z <- xyz_top[,3]
fitlm <- lm(x ~ poly(y, z, degree = 2, raw = T))
valfit <- predict(fitlm, yz_dorsal) #generate values from the fit
top_bd_dorsal_3 <- cbind(valfit, yz_dorsal)
valfit <- predict(fitlm, yz_posterior) #generate values from the fit
top_bd_posterior_3 <- cbind(valfit, yz_posterior)
# bottom
# xyz_bot <- xyz_msh[xyz_msh[,1] < -24000, ]
xyz_bot <- xyz_msh[xyz_msh[,1] < -24000 & xyz_msh[,2] < 68000, ]
x <- xyz_bot[,1]; y <- xyz_bot[,2]; z <- xyz_bot[,3]
fitlm <- lm(x ~ poly(y, z, degree = 2, raw = T))
valfit <- predict(fitlm, yz_dorsal) #generate values from the fit
bot_bd_dorsal_3 <- cbind(valfit, yz_dorsal)
valfit <- predict(fitlm, yz_posterior) #generate values from the fit
bot_bd_posterior_3 <- cbind(valfit, yz_posterior)
x_qua <- apply(cbind(top_bd_dorsal_3$valfit, bot_bd_dorsal_3$valfit), MARGIN = 1,
function(x) quantile(x, probs = bp7c_prob) )
DRA_layers_ME_dorsal_3 <- matrix(ncol = 3, nrow = 0)
for (j in 1:length(bp7c_prob)) {
cc <- cbind(x_qua[j,], yz_dorsal)
c2 <- cbind(cc[-nrow(cc),], cc[-1,]) %>%
t() %>%
matrix(., ncol = 3, byrow = T)
DRA_layers_ME_dorsal_3 <- rbind(DRA_layers_ME_dorsal_3, c2)
}
# manual translation
DRA_layers_ME_dorsal_3[,1] <- DRA_layers_ME_dorsal_3[,1] - 4000
x_qua <- apply(cbind(top_bd_posterior_3$valfit, bot_bd_posterior_3$valfit), MARGIN = 1,
function(x) quantile(x, probs = bp7c_prob) )
DRA_layers_ME_posterior_3 <- matrix(ncol = 3, nrow = 0)
for (j in 1:length(bp7c_prob)) {
cc <- cbind(x_qua[j,], yz_posterior)
c2 <- cbind(cc[-nrow(cc),], cc[-1,]) %>%
t() %>%
matrix(., ncol = 3, byrow = T)
DRA_layers_ME_posterior_3 <- rbind(DRA_layers_ME_posterior_3, c2)
}
bot2 <- DRA_layers_ME_posterior_3[(nrow(yz_posterior)*2-3):(nrow(yz_posterior)*2-2),]
top2 <- DRA_layers_ME_posterior_3[(nrow(DRA_layers_ME_posterior_3)-1):(nrow(DRA_layers_ME_posterior_3)), ]
DRA_layers_ME_posterior_3 <- DRA_layers_ME_posterior_3[-c(
(nrow(yz_posterior)*2-3):(nrow(yz_posterior)*2-2),
(nrow(DRA_layers_ME_posterior_3)-1):(nrow(DRA_layers_ME_posterior_3))), ]
cc <- cbind(x_qua[2:10,nrow(yz_posterior)], yz_posterior[nrow(yz_posterior),])
c2 <- cbind(cc[-nrow(cc),], cc[-1,]) %>%
t() %>%
matrix(., ncol = 3, byrow = T)
DRA_layers_ME_posterior_3 <- rbind(DRA_layers_ME_posterior_3, c2)
# add round corner
pt_corner <- rbind(bot2,
as.numeric(c(x_qua[2,nrow(yz_posterior)], yz_posterior[nrow(yz_posterior),])))
ptArc <- round_corner(pt_corner)
cc <- rbind(pt_corner[3,], ptArc, pt_corner[1,])
c2 <- cbind(cc[-nrow(cc),], cc[-1,]) %>%
t() %>%
matrix(., ncol = 3, byrow = T)
DRA_layers_ME_posterior_3 <- rbind(DRA_layers_ME_posterior_3, c2)
pt_corner <- rbind(top2,
as.numeric(c(x_qua[10,nrow(yz_posterior)], yz_posterior[nrow(yz_posterior),])))
ptArc <- round_corner(pt_corner)
cc <- rbind(pt_corner[3,], ptArc, pt_corner[1,])
c2 <- cbind(cc[-nrow(cc),], cc[-1,]) %>%
t() %>%
matrix(., ncol = 3, byrow = T)
DRA_layers_ME_posterior_3 <- rbind(DRA_layers_ME_posterior_3, c2)
# manual adjust
DRA_layers_ME_posterior_3[,1] <- DRA_layers_ME_posterior_3[,1] - 5500