Intro

One of the core features of elmr is to enable data in the FAFB EM space to be transformed into one of a number of standard light level template brains. This transformation depends on manually selected landmarks. We can inspect these landmarks and to try to

Setup

library(elmr)
library(knitr)
opts_chunk$set(fig.width=5, fig.height=5)
# set up for 3d plots based on rgl package
rgl::setupKnitr()
# frontal view
view3d(userMatrix=rgl::rotationMatrix(angle = pi, 1,0,0), zoom=0.7)

Landmarks

The standard landmarks used by elmr were specified by Davi Bock on the FAFB and JFRC2013 brains. We can plot the landmarks in the context of the JFRC2013 brain as follows.

plot3d(JFRC2013)
# note that the landmarks are in raw voxel coordinates and must be scaled
xyz=scale(elm.landmarks[,c("X","Y","Z")], scale = 1/voxdims(JFRC2013), center = FALSE)
xyz=as.data.frame(xyz)
spheres3d(xyz, col = ifelse(elm.landmarks$Use, "green", 'red'), radius = 4)

Landmamarks based transform

The standard elm.landmarks set is used within to the elmr package to register a transformation that can map FAFB->JFRC2013 space in conjunction with the xform_brain functions provided by the nat.templatebrains packages.

Let’s first look at the FAFB landmarks

# note that these positions are in nm units and do not need scaling
xyz.fafb=elm.landmarks[,c("X1","Y1","Z1")]
kable(head(xyz.fafb))
X1 Y1 Z1
591360.1 287783.6 97510.11
525666.5 172470.4 80994.73
595391.6 263523.1 84156.77
501080.2 253465.3 98233.34
528543.2 338099.9 19349.01
670999.9 179097.9 67561.69

Now we can try to look at the consistency of the transform between each landmark pair as follows:

  1. Drop out each landmark pair in turn
  2. Calculate position of that landmark in target space using remaining landmark pairs
  3. Measure difference between calculated position and manually defined position
# record transformed position with leave one out
xyzt.loo=xyz
for(i in 1:nrow(xyz)) {
  thisreg <- reglist(tpsreg(xyz.fafb[-i, ], xyz[-i, ]))
  xyzt.loo[i,]=xform(xyz.fafb[i, , drop=FALSE], thisreg)
}
deltas=sqrt(rowSums((xyzt.loo-xyz)^2))
hist(deltas)

Let’s looks at the deltas (distance between leave one out vs manually defined landmark position) for each point.

kable(cbind(elm.landmarks, delta=deltas))
Label Use X Y Z X1 Y1 Z1 delta
Pt-0 TRUE 702.1487 285.40275 209.73076 591360.1 287783.6 97510.11 4.495873
Pt-1 TRUE 571.4001 38.85996 287.05954 525666.5 172470.4 80994.73 4.439362
Pt-2 TRUE 715.8113 213.29936 217.39349 595391.6 263523.1 84156.77 2.110814
Pt-3 TRUE 513.0022 198.00197 217.79409 501080.2 253465.3 98233.34 6.576990
Pt-5 TRUE 602.0898 271.07851 17.30991 528543.2 338099.9 19349.01 11.787314
Pt-6 TRUE 867.0125 31.91925 276.22344 670999.9 179097.9 67561.69 4.625576
Pt-7 TRUE 935.2109 234.22952 351.51807 702703.9 251846.4 127865.89 7.413386
Pt-8 TRUE 935.2109 232.62607 405.50079 702633.2 250322.0 162468.60 11.278113
Pt-9 TRUE 935.2109 143.90200 414.05251 702633.2 201876.1 160251.58 11.335541
Pt-10 TRUE 509.7629 244.91917 362.20772 504182.4 244953.7 148700.67 10.118559
Pt-11 TRUE 509.2284 281.12828 403.48102 504050.2 252491.4 181628.86 5.534175
Pt-12 TRUE 509.2284 165.87835 401.67082 504050.2 204223.2 168536.93 8.871823
Pt-14 TRUE 1106.8503 372.72023 222.01044 784687.6 334367.1 98303.01 12.835334
Pt-17 TRUE 718.4719 242.59406 355.69466 604164.4 244944.2 144100.18 3.665553
Pt-18 TRUE 776.9947 213.33267 378.56799 639572.0 227939.5 149633.26 6.290179
Pt-19 TRUE 839.6388 287.10433 382.68932 660850.9 259437.0 158289.08 5.201948
Pt-20 TRUE 601.2203 288.54680 382.68932 552292.5 256070.9 166824.68 4.807353
Pt-21 TRUE 653.3550 221.57531 380.62865 573691.6 227458.6 153600.51 2.656390
Pt-23 TRUE 646.4883 588.60129 222.43882 575568.7 414728.5 162492.58 14.515807
Pt-24 TRUE 769.5130 605.94266 204.25736 628557.3 429391.3 150296.75 2.592245
Pt-25 TRUE 1150.4870 201.70284 365.00699 806460.8 229605.7 130248.64 13.362692
Pt-26 TRUE 1045.9179 168.05378 215.08408 758726.8 253891.8 60372.56 2.592961
Pt-27 TRUE 389.1564 177.03002 206.10784 443130.5 243562.0 92329.64 8.591221
Pt-39 TRUE 294.9652 198.55742 355.90246 397489.3 219837.0 154312.24 7.061507
Pt-40 TRUE 241.6744 211.20539 350.83541 374566.7 240128.5 154291.60 7.200079
Pt-41 TRUE 370.7507 378.18580 217.87953 427042.8 327135.6 129327.00 2.181738
Pt-42 TRUE 356.8177 360.25507 239.47496 417693.5 316017.5 134633.36 4.847158
Pt-43 TRUE 311.3398 452.29866 394.44954 442448.7 338329.6 223884.56 24.673536
Pt-44 TRUE 435.6423 511.76526 303.65791 476555.7 338329.6 192565.15 19.237017
Pt-45 TRUE 325.9132 597.71786 158.23298 413428.0 417820.1 149932.91 21.649187
Pt-46 TRUE 197.7791 449.11545 319.63260 379721.1 365104.0 194728.93 15.990273
Pt-47 TRUE 120.6732 274.52403 323.29304 306977.7 296201.8 158477.26 14.886355
Pt-48 TRUE 162.5141 275.24054 425.51338 358354.3 244585.2 203131.63 16.122111
Pt-49 TRUE 413.3632 161.16906 316.82330 455995.0 219330.2 129508.15 4.643610
Pt-50 TRUE 542.6337 173.56902 150.26896 501073.6 258902.8 68029.46 9.397889
Pt-51 TRUE 732.8606 262.06885 208.87088 605899.8 285350.0 89498.26 3.135272
Pt-52 TRUE 763.4152 262.35759 225.77664 622311.7 278256.8 91584.51 5.725675
Pt-53 TRUE 673.4191 235.90844 112.99012 579791.1 293613.3 48661.11 8.988793
Pt-54 TRUE 522.3014 213.16795 94.27100 501220.4 290557.8 48661.11 8.007639
Pt-55 TRUE 918.1357 217.61283 104.17984 698302.0 314347.2 24435.13 12.130931
Pt-56 TRUE 1021.5445 379.96095 136.36467 742073.3 362250.9 64369.84 4.104513
Pt-57 TRUE 1212.4784 422.50773 128.28765 817108.7 368563.9 87998.78 18.177230
Pt-58 TRUE 1249.0229 401.20869 365.68265 829769.9 313224.9 176053.38 4.574636
Pt-59 TRUE 1174.4594 631.83024 248.45789 818414.3 433119.5 173412.53 22.373505
Pt-60 TRUE 1077.1313 573.66563 299.84596 747375.4 404862.4 173412.53 17.742575
Pt-61 TRUE 1112.7070 339.48705 401.72907 756090.2 286552.3 173412.53 19.347517
Pt-62 TRUE 1306.1523 263.30353 411.07188 859611.5 246939.6 173412.53 7.594850
Pt-63 TRUE 1328.8859 469.20296 354.47855 859687.9 333999.7 192126.24 12.448687
Pt-64 TRUE 549.6954 583.20048 162.32429 509921.5 403489.6 147127.08 12.925469
Pt-65 TRUE 765.1116 188.70754 178.33252 616575.5 271178.7 58176.36 7.625301
Pt-66 TRUE 785.3950 425.64347 15.25080 632860.0 399076.2 37528.28 6.358484
Pt-67 TRUE 491.9622 447.69119 309.37694 489103.1 329705.4 176917.26 13.309858
Pt-68 TRUE 636.6000 582.23100 304.45079 548824.8 380895.4 183068.00 18.791822
Pt-69 TRUE 751.5368 622.07216 324.54432 622633.7 397958.8 198940.87 10.800934
Pt-70 TRUE 248.6964 386.77664 133.64117 375957.0 360354.0 115925.73 14.752904
Pt-71 TRUE 371.3233 241.13746 246.95856 425172.9 260481.9 115925.73 6.101949
Pt-72 TRUE 125.3698 356.45723 424.67776 349795.9 298764.2 229860.66 13.397052
Pt-73 TRUE 366.0449 417.33611 390.53431 447017.2 298565.8 211606.86 8.953431
Pt-74 TRUE 297.1745 286.05935 401.86582 401391.2 252122.0 196178.45 12.326181
Pt-75 TRUE 120.6515 261.44368 422.62512 307822.1 252122.0 202694.69 19.538183
Pt-76 TRUE 371.7625 362.30048 386.60631 445986.7 279612.4 193734.86 5.209407
Pt-77 TRUE 863.6581 20.51800 305.28403 666626.3 172167.2 73062.24 6.422350
Pt-78 TRUE 877.3450 69.27787 262.22208 675880.3 195927.0 63101.97 4.434200
Pt-79 TRUE 887.6924 116.68930 198.29597 683285.2 236316.9 46525.29 2.511654
Pt-80 TRUE 456.7993 163.75185 410.38070 486884.6 200114.5 167707.33 4.602266
Pt-81 TRUE 457.2808 92.26043 414.81480 486884.6 170297.8 145557.76 9.982881
Pt-82 TRUE 543.1608 46.12096 309.74346 514997.5 166251.2 96573.13 5.349113
Pt-83 TRUE 610.5969 134.29354 334.08200 548304.2 204859.4 120887.29 3.590523
Pt-84 TRUE 687.3255 150.56257 339.67920 590275.8 207454.9 120887.29 2.896198
Pt-85 TRUE 637.3402 37.40585 246.20518 558273.3 175584.6 62700.93 6.317059
Pt-86 TRUE 601.4145 48.82580 235.55270 530820.9 183745.3 62633.07 6.675230

And then visualise them in 3D, colouring each point by the magnitude of the delta.

clear3d()
view3d(userMatrix=rgl::rotationMatrix(angle = pi, 1,0,0), zoom=0.7)
jet.colors <-
  colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
                     "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))

plot3d(JFRC2013)
spheres3d(xyz, col = jet.colors(10)[cut(deltas, 10)], radius = 4)
# nb these are numbered from 0 like the elm landmarks
texts3d(xyz, texts = 1:nrow(xyz)-1, adj = c(1,1))

Based on this, I would definitely suggest checking the following landmarks

59/32 27 30 52 18

Compare the leave one out landmark position with manually defined position

clear3d()
view3d(userMatrix=rgl::rotationMatrix(angle = pi, 1,0,0), zoom=0.7)
for(i in 1:nrow(xyz)){
  segments3d(rbind(xyz[i, , drop=FALSE], xyzt.loo[i, , drop=FALSE]), col=jet.colors(10)[cut(deltas, 10)][i], lwd=3)
}
points3d(xyz)
texts3d(xyz, texts = 1:nrow(xyz)-1, adj = c(1,1))
plot3d(JFRC2013)

Now of course this is not just a measure of consistency in each landmark point, but it will also depend strongly on the extent to which the neighbouring points can define the required registration; in areas requiring substantial non-rigid deformation this may be a problem; this will be further emphasised when the neighbouring points are distant. We can check the latter point like so:

library(nabor)
nearest3=rowMeans(nabor::knn(xyz, k=4)$nn.dists[,-1])
plot(deltas~nearest3)
deltabynearest3=lm(deltas~nearest3)
abline(deltabynearest3)

summary(deltabynearest3)
## 
## Call:
## lm(formula = deltas ~ nearest3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.9280  -2.5874  -0.5281   2.4646  16.0718 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.26996    1.69207   1.342    0.184    
## nearest3     0.16138    0.03635   4.440 3.33e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.993 on 69 degrees of freedom
## Multiple R-squared:  0.2222, Adjusted R-squared:  0.211 
## F-statistic: 19.71 on 1 and 69 DF,  p-value: 3.331e-05

So there is a correlation between neighbour distance and our calculated delta, but with an R^2 of about 0.21 it’s not super strong. Furthermore it could also be explained by Davi placing more landmarks in areas where he is more certain.

Perhaps the ideal thing would be to find a set of landmarks on one side of the brain and then find exactly the same landmarks on the other side of the brain. For the JFRC2013 brain we can exactly mirror points from one side to the other using an image-based registration distributed with the nat.flybrains package.

We can find the mirror image position of all the points that we used like so:

#  nb this is conditioned on being able to find CMTK
if(isTRUE(nzchar(cmtk.bindir()))){
  xyzm=mirror_brain(xyz, brain = JFRC2013)
  
  clear3d()
  spheres3d(xyz, col='red', rad=2)
  spheres3d(xyzm, col='green', rad=2)
  plot3d(JFRC2013)
  view3d(userMatrix=rgl::rotationMatrix(angle = pi, 1,0,0), zoom=0.7)

}

One could then write out a table of landmarks including these mirrored positions like so:

xyzm.pixels=scale(xyzm, scale = voxdims(JFRC2013), center = FALSE)
colnames(xyzm.pixels)=c("XM","YM","ZM")
write.table(cbind(elm.landmarks, xyzm.pixels), 'elm.landmarks-with-mirror.tsv',
            col.names = T, sep='\t', row.names = F)

If one were then to take the mirror image landmarks in JFRC2013 space and manually pick the corresponding locations in EM space, one should be able to compare two different paths to the same point (one mirroring in FAFB, the other mirroring in JFRC2013) and see whether there is a difference in the calculated position. Differences would presumably reflect the uncertainty in picking positions in the EM volume (and a small contribution from the mirroring uncertainty in JFRC2013, which nevertheless expect to be <1 µm based on analysis in Manton et al 2014).