On variational solutions for whole brain 1 serial-section histology using the 2 computational anatomy random orbit model 3

10 This paper presents a variational framework for dense diffeomorphic atlas11 mapping onto high-throughput histology stacks at the 20 μm meso-scale. The 12 observed sections are modelled as Gaussian random fields conditioned on a 13 sequence of unknown section by section rigid motions and unknown diffeo14 morphic transformation of a three-dimensional atlas. To regularize over the 15 high-dimensionality of our parameter space (which is a product space of the 16 rigid motion dimensions and the diffeomorphism dimensions), the histology 17 stacks are modelled as arising from a first order Sobolev space smoothness 18 prior. We show that the joint maximum a-posteriori, penalized-likelihood 19 estimator of our high dimensional parameter space emerges as a joint op20 timization interleaving rigid motion estimation for histology restacking and 21 large deformation diffeomorphic metric mapping to atlas coordinates. We 22 show that joint optimization in this parameter space solves the classical cur23 vature non-identifiability of the histology stacking problem. The algorithms 24 are demonstrated on a collection of whole-brain histological image stacks 25 from the Mouse Brain Architecture Project. 26 not peer-reviewed) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was . http://dx.doi.org/10.1101/271692 doi: bioRxiv preprint first posted online Feb. 25, 2018;

shows the components of the model for the histology stacking prob-139 lem. We define the mouse brain to be sectioned as a dense three-dimensional 140 (3D) object I(x, y, z), (x, y, z) ∈ R 3 , modelled to be a smooth deformation 141 of a known, given template I 0 so that I = I 0 • ϕ −1 for some invertible dif-142 feomorphic transformation ϕ. The Allen Institute's mouse brain atlas [26] 143 (CCF 2017) is taken as the template.
144 Figure 1: The histological sectioning model; the template I 0 , the mouse brain in the orbit I ∈ I and observed histological sections J i , i = 1, . . . , n. The Sobolev image intensity prior and the shape prior are depicted in the top row. The model shows the template and mouse brain as elements of the same orbit I 0 , I ∈ I, so there exists diffeomorphism I = I 0 • ϕ −1 , ϕ ∈ Diff.
Distinct from volumetric imaging such as MRI which delivers a dense 3D metric of the brain, the histology procedure (bottom row, Figure 1) consisting of sectioning, staining, and imaging generates a jitter process which randomly 147 translates and rotates the stack sections. Denote the rigid motions acting on 148 the 2D sectioning planes R i : R 2 → R 2 , 149 R i (x, y) = (cos θ i x + sin θ i y + t x i , − sin θ i x + cos θ i y + t y i ) , (x, y) ∈ R 2 , (1) with θ i the rotation angle and (t x i , t y i ) ∈ R 2 the translation vector in section 150 i. The histology stack J i (x, y), (x, y) ∈ R 2 , i = 1, . . . , n, is a sequence of 2D 151 jittered image sections under smooth deformation of the atlas in noise: Modeling the photographic noise as Gaussian and conditioning on the Here α i is a weighting factor dependent on the noise of each section such that with v t the Eulerian velocity taking values in R 3 , identity the identity map-174 ping. The top row of Figure 1a shows that each ϕ has an inverse and that the 175 random orbit model assumes any individual brain I ∈ I can be generated 176 from the exemplar under the action of the diffeomorphism, so that for some This square-metric is used as a quadratic potential for the smoothness prior 186 between images I, I ∈ I [28, 29] minimizes the action 187 ρ 2 (I, I ) = min  for the z-partial derivatives, The gradient is forced to 0 at the boundaries of the image. Model the random sectioning with section-independent jitter as a product density π(R) = i π(θ i , t x i , t y i ), the priors centered at identity. Generating MAP estimates of the rigid motions generates the MAP estimator of the histology restacking problem denoted as Since the diffeomorphisms are infinite dimensional, the maximization of the 213 log-likelihood function with respect to a function with the deformation penalty 214 is termed the "penalized-likelihood estimator". Conditioned on the known 215 atlas, the augmented random variables to be estimated are (R 1 , . . . , R n , ϕ) ∈ 216 (R 3n × Diff).

218
Given histology stack J i (x, y), (x, y) ∈ R 2 , i = 1, . . . and reconstructed stack 219 I R (·, z i ) = J i • R i (·), i = 1, . . . , n modelled as conditionally Gaussian random 220 fields conditioned on jitter and smooth dormation of the template. The joint 221 MAP, Penalized-Likelihood estimators arg max R,v log π(R, v|J) given by The MAP, Penalized-Likelihood estimators satisfy with · 2 2 denoting the norm per z-axis section: We call this the atlas-informed model. The first two prior terms of 225 (9) control the smoothness of template deformation and the realigned target 226 image stack, with the third keeping the rigid motions close to the identity.

227
The last term is the "log-likelihood" conditioned on the other variables.

228
The optimization for the R * rigid-motions is not decoupled across sec-   When there is no atlas available this is equivalent to setting α i small and becomes a MAP rigid motion restacking of the sections: We term this the atlas-free model. The gradient of the rigid motions with   2. Update LDDMM for diffeomorphic transformation of atlas coordinates: 3. Deform atlas I 0 • ϕ new−1 and generate new histology image stack:

Return to
Step 1 until convergence criterion met.

253
The form of the gradients for the rigid motions is given in Appendix D.

254
The LDDMM update solutions are given by Beg [25].

256
The algorithm described above is applied to Nissl histological stacks using 257 the Allen Institute's mouse brain atlas as a template. The Allen Mouse Brain

258
Atlas is a micron-scale atlas that includes annotated Nissl-stained images at this is a masking procedure on the cost function that allows matching between 289 a whole atlas brain and some target which is a partial, or subset of a whole 290 brain.  The model was applied to binary image phantoms in order to examine the 302 "curvature" problem in which a 3D curved object cannot be accurately re-303 constructed after being sectioned. This is illustrated in Figure 3. We pro- 327 Figure 4: Transformed grids illustrating the difference in the mapping deformation from the atlas-free methods from (A) to histology stack target (B) versus the atlas-informed algorithm which produces (C).

328
A similar experiment was performed using the Allen mouse brain atlas as 329 the 3D phantom. A target histology stack was generated by sectioning the 330 Allen atlas in simulation and applying random rigid transforms to its coronal 331 sections. The atlas images were sampled at 40 µm isotropic voxels. This is 332 depicted in Figure 5a. A simulated atlas was generated by applying a given 333 random diffeomorphism to the Allen atlas. This random diffeomorphism is 334 depicted in Figure 5c. The histology stacks were then reconstructed and   Figure 5: Atlas phantom simulation to validate recovery of sectioning parameters and diffeomorphic shape difference. a) The ground truth target I is sectioned to generate the observed target J i . b) Transformed grids illustrating the brain phantom atlas (top) shown mapped onto the histological stack using the atlas-free algorithm (bottom left) and the atlas-informed algorithm (bottom right). c) The ground truth diffeomorphism to be recovered.   mapping is much smoother and has many fewer wrinkles than the atlas-free 394 mapping. This is seen clearly in Figure 9.   The Green's kernel is translation invariant and takes the form

553
We can write the gradient with respect to the components of R (translation vector t and rotation matrix r parametrized by rotation angle θ and section number z), where ∇ X is the 2D in-plane gradient, σ JJ is the weighting factor on the image smoothness prior. Rotations and translations are penalized by a regularization prior centered at identity ( θ σ 2 regR and t(z) σ 2 regT , respectively), where σ regR and σ regT are weighting factors on the rotation and translation priors.