Large-scale 3D histology reconstruction of the human

brain with MRI as volumetric reference generally requires

reassembling the tissue blocks into the MRI space, prior to

any further reconstruction of the histology of the individual

blocks. This is a challenging registration problem, particularly in the frequent case that blockface photographs of

paraffin embedded tissue are used as intermediate modality, as their contrast between white and gray matter is rather

low. Here we propose a registration framework to address

this problem, relying on two key components. First, blocks

are simultaneously aligned to the MRI while exploiting the

spatial constraints that they impose on each other, by means

of a customized soft shape constraint (similarly to a jigsaw

puzzle). And second, we adopt a hierarchical optimization

strategy that capitalizes on our prior knowledge on the slicing

and blocking procedure. Our framework is validated quantitatively on synthetic data, and qualitatively on the histology

of a whole human hemisphere

Crampsie, S

Holton, JL

Iglesias, JE

Jaunmuktane, Z

Mancini, M

Thomas, DL

English

Crossref

Hierarchical Joint Registration of Tissue Blocks With Soft Shape Constraints For Large-Scale Histology of The Human Brain

UCL Discovery

HIERARCHICAL JOINT REGISTRATION OF TISSUE BLOCKS WITH SOFT SHAPECONSTRAINTS FOR LARGE-SCALE HISTOLOGY OF THE HUMAN BRAINM. Mancini 1, S. Crampsie 2, D.L. Thomas 3,4, Z. Jaunmuktane 2, J.L. Holton 2, J.E. Iglesias 11 Center for Medical Image Computing, Dept. Med. Phys. & Biomed. Eng., University College London2 Queen Square Brain Bank, Institute of Neurology, University College London3 Neuroradiological Academic Unit, Institute of Neurology, University College London4 Leonard Wolfson Experimental Neurology Centre, Institute of Neurology, University College LondonABSTRACTLarge-scale 3D histology reconstruction of the humanbrain with MRI as volumetric reference generally requiresreassembling the tissue blocks into the MRI space, prior toany further reconstruction of the histology of the individualblocks. This is a challenging registration problem, partic-ularly in the frequent case that blockface photographs ofparaffin embedded tissue are used as intermediate modal-ity, as their contrast between white and gray matter is ratherlow. Here we propose a registration framework to addressthis problem, relying on two key components. First, blocksare simultaneously aligned to the MRI while exploiting thespatial constraints that they impose on each other, by meansof a customized soft shape constraint (similarly to a jigsawpuzzle). And second, we adopt a hierarchical optimizationstrategy that capitalizes on our prior knowledge on the slicingand blocking procedure. Our framework is validated quanti-tatively on synthetic data, and qualitatively on the histologyof a whole human hemisphere.Index Terms— Registration, ex vivo MR, histology, brain1. INTRODUCTIONHistology and magnetic resonance imaging (MRI) are twopowerful tools to study biological systems, in particular thebrain. However, they target very different scales: on the onehand, histology allows microstructural examination of smallsamples with high levels of detail; on the other hand, MRIoffers a global overview of the entire brain. These comple-mentary approaches hold great potential in building whole-brain atlases with histology-level resolution, the goal of sev-eral large-scale histology initiatives [1, 2].Building such atlases presents several challenges. First,the fact that histological analysis yields a set of spatially dis-torted 2D sections, due to cutting, processing, sectioning, andmounting of tissue. Recovering the 3D structure of histolog-ical sections, which often involves registration to a referencevolume (typically MRI), is a difficult, largely studied problemin the literature [3]. A related challenge is the large differ-ences in contrast and resolution between MRI and histology,which, along with histological artifacts (e.g., tears, folding)make registration between the two modalities difficult. An-other obstacle in large-scale histology studies of the humanbrain is the size of the human brain, which requires tailoredwhole-brain microtomes that are only available in very fewselected sites around the world.An alternative to whole-brain microtomes is cutting thebrain into blocks, which enables tissue processing with stan-dard equipment and techniques – and hence applicabilityat nearly every research site. The cutting procedure is con-strained by both the necessity to preserve structures of interest(e.g., hippocampus [4], subthalamic nuclei [5]), and the prac-tical limits imposed by the standard size of cassettes andglass slides (i.e., a maximum of 74× 52 mm), for subsequenttissue processing, sectioning and staining. A resulting issueis hence the need to solve a “jigsaw puzzle” to reassemble theblocks, before further reconstruction of individual blocks.The registration between blocks and MRI is commonlycarried out by means of an intermediate modality, whichdoes not suffer from nonlinear geometric distortions, mostfrequently blockface photographs. The literature on MRI-blockface registration is rather sparse. A frequent approachis to align each volume of blockface photographs to the MRIindependently. However, this registration is not easy (de-spite its linear nature, and even with reasonable initialization)because of its whole-to-part nature, particularly when thecontrast of the blockface images is poor – which is the casefor the ubiquitous paraffin embedding.In this paper, we present a hierarchical approach to theblockface-MRI registration problem, which jointly registersall blocks, while ensuring that they do not overlap or leavegaps in between (i.e., the pieces in the jigsaw puzzle fit). Thisis achieved by a soft constraint imposed on the shapes of theblocks, and which is integrated in the loss function of theregistration. The optimization is carried out in a hierarchi-cal fashion, which takes advantage of our prior knowledge onthe geometry of the brain cut and blocking procedures.2. METHODS2.1. General registration frameworkLet V (x) be a volumetric MRI scan, defined in a spatial do-main x ∈ ΩV . Likewise, let {Bn(x)}n=1,...,N beN volumesof blockface photographs to be registered to the MRI scan, de-fined in their own domains {Ωn}. In addition, we assume theavailability of binary masks MV (x) and {Mn(x)}n=1,...,Nfor the MRI scan and blockface volumes, respectively, whereMV ,Mn ∈ {0, 1}. These masks are assumed to be obtainedwith an automated segmentation algorithm, and hence to benoisy. Our goal is then to find a set of linear spatial trans-forms, parameterized by their 4 × 4 transformation matrices{Θn}n=1,...,N , that best align the blocks. This is achieved byminimizing the following general loss function:L = 1NN∑n=1D [V (x), Bn(T (x,Θn))] +γ1NN∑n=1R(Θn)+γ2C [MV (x),M1(T (x,Θ1)), . . . ,MN (T (x,ΘN ))] , (1)where D is a dissimilarity metric between the MRI and eachof the registered blocks; R is a regularizer that encodes priorknowledge on the transforms; C is a soft constraint term thatpenalizes implausible spatial arrangement of the blockfacevolumes based on their masks; and γ1, γ2 are relative weights.2.2. Modeling of spatial constraints with CThe term C encodes our prior knowledge on the spatial rela-tionships between the different blocks. If the blockface vol-umes were perfectly cut (i.e, without any deformation) andsegmented, then we would expect them to fit exactly, andcould use hard constraints to prevent them from overlappingor being disjoint. In practice, this is never the case, so weresort to a soft, continuous regularizer, which simultaneouslypenalizes overlap and gaps, while encouraging the externalboundary of the blocks to fit that of the MRI:C =∑x∈ΩV∣∣∣∣∣MV (x)−∑nMn[T (x,Θn)]∣∣∣∣∣ . (2)When a gap between blocks exists at a given voxel, the sumof masks is zero, and a penalty of one is incurred. When twomasks overlap at a voxel, the sum is two, and a penalty ofone is incurred as well. Higher penalty is paid when mul-tiple blocks overlap at the same voxel. Also, by consider-ing all voxels inside the MRI mask in the sum, penalty isalso incurred when the union of the blockface masks does notmatch MV (x) well, effectively fitting the outer boundaries ofthe blocks and the MRI mask. We note that we also tried aquadratic function in Eq. 2, but more robust results were ob-tained with the L1 norm in pilot experiments.2.3. HierarchyRather than optimizing all {Θn} directly in Equation 1, analternative is force groups of transforms to share the sameparameters in the earlier iterations of the algorithm, to reflectour prior knowledge on the cutting procedure. For instance,one could group the blocks from the same slice in the firstlevels, and then allow each block to have its own transform infiner levels.This hierarchy can be reflected by replacing Θn byΘhj(n;h) in Eqns. 1 and 2, where j(n;h) groups the blocksthat share parameters Θhj at hierarchy level h (where h =1, . . . ,H). We use Jh to refer to the number of differenttransforms at hierarchy level h, such that Jh < Jh+1, and(typically) JH = N . We also note that the number of param-eters of Θhj can be different at every level, e.g., could increasefrom 3 (translation transform) to 6 (rigid), 7 (similarity) or 12(affine), as we progress through the hierarchy.2.4. Model instantiation and implementation detailsWe complete the framework with choices for the dissimilar-ity metric D and the regularizer R. For the dissimilarity, weuse the negated cross correlation between the edge maps (gra-dient magnitudes) of the two images. Even though mutualinformation might be a more natural alternative, our choiceproved to be more robust and perform better in preliminaryexperiments. For the regularizer, we first computed the logJacobian linear transformation matrices, and quadratically pe-nalized deviations from an a priori value η that describes theexpected shrinkage due to tissue processing (which can be de-termined empirically):R(Θn) = [log |Θn| − η]2. (3)We can then combine Eqns. 1-3 to obtain the final loss func-tion at hierarchy level h:L =− 1NN∑n=1NCC[V E(x), BEn(T (x,Θhj(n;h)))]+γ1JhJh∑j=1(log |Θhj(n;h)| − η)2+ γ2∑x∈ΩV∣∣∣∣∣MV (x)−∑nMn[T (x,Θhj(n;h))]∣∣∣∣∣ , (4)where V E and {BEn} denote the edge maps of the MRI vol-ume and blocks, respectively, and NCC is the normalizedcross correlation. In order to minimize Eq. 4, we use the L-BFGS algorithm [6]. We note that we use linear interpolationfor Mn[T (x,Θhj(n;h))], otherwise derivatives would be zeroalmost everywhere.3. EXPERIMENTS AND RESULTS3.1. Data3.1.1. Real dataWe used a left human hemisphere donated for research tothe UCL Queen Square Brain Bank, through a program withapproval from Research Ethics Committes. MRI data wereacquired ex vivo with a T2 CISS sequence (TR = 6.37 ms,TE = 2.88ms, BW = 449 Hz/Px, ten averages), with Flu-orinert immersion to avoid susceptibility artifacts. For thehistological processing, the hemisphere was first separatedinto cerebrum, cerebellum and brainstem. Each of them wassubsequently cut into 10mm thick slices, following the coro-nal, sagittal and axial directions, respectively. The cerebrumslices were further cut in blocks able to fit large cassettes(74×52 mm), while preserving the main subcortical struc-tures and taking pictures before and after this cut (brain cutphotographs). Each block was then processed and embeddedin paraffin. The blocks were sectioned at 25μm thickness,while taking blockface photographs of every section. Theblockface photographs were stacked in volumes, correctedfor perspective [3], and co-registered with the SIFT algorithm[7]. The rearranging of the blockface volumes was initializedwith registrations derived from the brain cut photographs.3.1.2. Synthetic dataSynthetic data were used to evaluate the registration methodquantitatively, as obtaining perfectly aligned MRI and block-face volumes is extremely challenging. Here we used 50randomly selected multimodal (T1/T2) scans from the IXIdataset (brain-development.org/ixi-dataset/),letting the T1 play the role of reference (V ), and deformingthe T2 to play the role of blockface volumes ({Bn}).The synthetic deformation mimicked the brain cut proce-dure as follows: 1. Segmentation of the T1 with FreeSurfer [8]to extract the left hemisphere and subdivide it into cerebrum,cerebellum and brainstem; the union of these three masks,corrupted with noise at the boundary, provides the mask MV ,and the background of the T1 volume is corrupted by Gaus-sian noise. 2. Subdivision of the cerebrum, cerebellum andbrainstem of the T2 into 10 mm thick slices, with coronal,sagittal, and axial orientation, respectively. 3. Further sub-division of cerebrum slices into blocks, using a samplingprocedure to attempt to minimize both the number of blocks(constrained to fit in a 74×52 mm cassette) and the number ofcrossed structure boundaries. 4. Corrupting blocks and cor-responding masks the same way as the T1. 5. Corrupting thepositioning of the blocks with random similarity transforms,first in groups corresponding to slices, and then each blockindividually – attempting to reproduce the type of initializa-tion that one typically achieves with brain cut photographs.Sample T1 and corrupted T2 volumes can be seen in Fig. 1a,b.3.2. Experimental setupFor the quantitive evaluation, we attempted to recover theoriginal location of the T2 blocks by minimizing Eq. 4 withfour competing approaches. The baseline approach consistedof independent registration of the blocks, i.e., with γ2 = 0 andno hierarchy. The second approach was the full model withγ2 > 0, and a hierarchy with four levels: 1. Blocks grouped inslices, allowed to translate and rotate only on the plane paral-lel to the slice; 2. Blocks may individually translate and rotateon the plane parallel to the block, and a global scaling factor isshared among blocks; 3. Same as 2, but each block has its ownscaling factor; and 4. Full similarity transform for each block.Finally, we attempted to separate the contributions of the hi-erarchy and the soft spatial constraints by evaluating a modelwith γ2 > 0 by only the finest hierarchy level (“Overlap”),and a model with full hierarchy but γ2 = 0 (“Hierarchy”).Models were evaluated with statistics on the registration error(in mm) across all voxels in the masks of the T2 blocks.Throughout the experiments, we used the following pa-rameters: γ1 = 1, γ2 = 1 (or 0 if for the baseline and “Hierar-chy”), and η = log(3/4) – which was empirically determinedto be the approximate shrinking factor due to tissue process-ing in our setup. The full model with γ1 = γ2 = 1 was usedfor the real data. In all cases, we operated at 1 mm resolu-tion, but the resulting linear transforms can be applied to theoriginal data at any resolution.3.3. Results on synthetic dataTable 1 shows the mean, standard deviation, and robust max-imum (95th percentile) of the registration error across the 50volumes, for the four competing methods. Using the hierar-chy alone provides modest reductions (∼5%) on the mean andmaximum errors. When spatial constraints are used with nohierarchy, larger improvement is achieved with respect to thebaseline, both in terms of mean and maximum errors. Whenthe hierarchy is added on top of spatial constraints, the meanerror barely decreases, but the robustness of the method fairlyimproves, as reflected by the reductions in the standard de-viation of the error and, especially, its robust maximum (0.4mm decrease). Compared with the baseline method, the fullmodel achieves a 14% reduction in mean error and 26% de-crease in maximum error, thanks to the robustness providedby the spatial constraints. Sample reconstructions producedby the different methods are shown in Figs. 1 and 2.3.4. Results on real dataFigure 3 shows the results obtained on the real data. Evenwith the already mentioned differences between MRI andblockface photographs, the final alignment between the twomodalities shows good qualitative details, e.g. the alignmentof the ventricle, the smooth separation between cerebellumand cerebrum.Method Mean Std. dev. Robust Max.Baseline 2.15 1.11 4.48Hierarchy 2.08 1.08 4.25Overlap 1.86 0.93 3.72Full model 1.85 0.82 3.34Table 1. Distance error statistics for each method, in mm.The robust maximum is the 95th percentile of the error.(a)																																																																				(b)				(c)																																																																					(d)				Fig. 1. Sample sagittal, coronal and axial slices of: (a) Ref-erence T1 volume; (b) Distorted T2 blocks; (c) T2 blocks re-constructed with baseline method; and (3) Reconstructed withproposed (full) model. Notice how our model is able, e.g., torecover the largely displaced block above the occipital hornof the lateral ventricle.INITIALIZATION																													BASELINE																																PROPOSEDFig. 2. 3D renderings of the block masks {Mn} correspond-ing to the synthetic blockface volumes from Fig. 1. Note theimprovement in, among others, the block mentioned in thecaption of the aforementioned figure, displayed in white here.4. CONCLUSIONA method for jointly registering histological volumes to an exvivo MRI scan by means of blockface photographs has beenproposed in this paper. Future work will be focused on a mul-tiresolution version of the algorithm, and on evaluation onreal data using pairs of corresponding, manually placed land-marks, or using manual segmentations.This method will be used as part of a pipeline for full his-tological reconstruction for whole brain hemispheres, eventu-ally leading to the construction of a high-resolution compu-tational atlas of the human brain based on 3D reconstructedhistology.a)b)Fig. 3. Sample sagittal, coronal and axial slices of: (a) Refer-ence T2 volume; (b) Reconstructed blockface volume. Noticehow most structures, in particular subcortical ones, present asmooth profile across blocks.Acknowledgement: Research primarily supported by the Eu-ropean Research Council (ERC) through the Starting GrantNo 677697 (project “BUNGEE-TOOLS”), awarded to JEI.5. REFERENCES[1] K. Amunts, C. Lepage, L. Borgeat, et al., “BigBrain: Anultrahigh-resolution 3D human brain model,” Science,vol. 340, no. 6139, pp. 1472–1475, 2013.[2] S. Ding, J.J. Royall, S.M. Sunkin, et al., “Comprehen-sive cellular-resolution atlas of the adult human brain,” J.Comp. Neurol., vol. 524, no. 16, pp. 3127–3481, 2016.[3] J. Pichat, J.E. Iglesias, T. Yousry, et al., “A survey ofmethods for 3D histology reconstruction,” Medical ImageAnalysis, vol. 46, pp. 73–105, 2018.[4] J.E. Iglesias, J.C. Augustinack, K. Nguyen, et al., “Acomputational atlas of the hippocampal formation usingex vivo, ultra-high resolution MRI: Application to adap-tive segmentation of in vivo MRI,” NeuroImage, vol. 115,pp. 117 – 137, 2015.[5] L.A. Massey, M.A. Miranda, L. Zrinzo, et al., “High reso-lution MR anatomy of the subthalamic nucleus: Imagingat 9.4T with histological validation,” NeuroImage, vol.59, no. 3, pp. 2035 – 2044, 2012.[6] D.C. Liu and J. Nocedal, “On the limited memory BFGSmethod for large scale optimization,” Mathematical pro-gramming, vol. 45, no. 1-3, pp. 503–528, 1989.[7] D.G. Lowe, “Object recognition from local scale invari-ant features,” in Proc. ICCV, 1999, vol. 2, pp. 1150–1157.[8] B. Fischl, D. Salat, E. Busa, et al., “Whole brain segmen-tation: automated labeling of neuroanatomical structuresin the human brain,” Neuron, vol. 33, pp. 341–355, 2002.

https://discovery.ucl.ac.uk/10083232/1/ISBI19_0690_FI.pdf

Hierarchical Joint Registration of Tissue Blocks With Soft Shape Constraints For Large-Scale Histology of The Human Brain

Abstract

Similar works

Full text

Available Versions

Crossref

UCL Discovery