Cross-modal coherent registration of whole mouse brains

Recent whole-brain mapping projects are collecting large-scale three-dimensional images using modalities such as serial two-photon tomography, fluorescence micro-optical sectioning tomography, light-sheet fluorescence microscopy, volumetric imaging with synchronous on-the-fly scan and readout or magnetic resonance imaging. Registration of these multi-dimensional whole-brain images onto a standard atlas is essential for characterizing neuron types and constructing brain wiring diagrams. However, cross-modal image registration is challenging due to intrinsic variations of brain anatomy and artifacts resulting from different sample preparation methods and imaging modalities. We introduce a cross-modal registration method, mBrainAligner, which uses coherent landmark mapping and deep neural networks to align whole mouse brain images to the standard Allen Common Coordinate Framework atlas. We build a brain atlas for the fluorescence micro-optical sectioning tomography modality to facilitate single-cell mapping, and used our method to generate a whole-brain map of three-dimensional single-neuron morphology and neuron cell types. mBrainAligner is a cross-modal registration platform for whole mouse brains imaged with different modalities. In addition, a fluorescence micro-optical sectioning tomography-based mouse brain atlas has been generated.

A critical technique required to use these data is registration (alignment) of brains, their compartments, and various neuronal structures including somas, dendrites and axon arbors. Data that are acquired across subjects, modalities, time or even species must be mapped onto a canonical coordinate space. Its importance is obvious in multiple ways. First, registration of datasets of different modalities generated in different projects provides a valuable resource for ever-expanding studies of neuron wiring and functions. Second, the registration to a common standard atlas enables digital delineation and quantification of brain anatomy and neuron morphology and other attributes such as transcriptomics. Third, registration allows quantitative comparison of features of neurons across samples under different conditions, thus enabling analysis of distinct aspects of neurobiology in a coordinated manner.
Established brain atlases such as the Allen Common Coordinate Framework atlas (CCFv3) (using STPT imaging) 29,30 and its derived LSFM-domain atlas 31 have helped neurobiologists to study gene expression patterns, brain anatomy and neural circuits. Currently, fMOST imaging is increasingly used for whole-brain single-cell analysis. High-resolution images in this modality are being generated (for example in the BICCN program) to produce large single-neuron morphology databases 32 . It is crucial to register such data to the Allen brain atlas to compare mapped neurons and populations of neurons generated using different experimental conditions and modalities. Despite a number of existing studies on registration of whole-brain images within one single imaging modality, including for smaller nervous systems from Caenorhabditis elegans 33 , fruitfly 34 , to zebrafish 35,36 , there has been little reported work on cross-modal registration for prevailing bioimaging methods, especially considering the increasing diversity of intensity, texture and anatomy of the structures imaged. In particular, nonuniform characteristics of brain structures makes it challenging to coregister datasets using a number of existing registration tools developed in the medical imaging field [37][38][39] or tools adapted from off-the-shelf registration toolkits [40][41][42] . Recently, deep neural network (DNN)-based methods were proposed to improve registration efficiency and accuracy further 43,44 . Yet, application of DNNs to cross-modal bioimaging registration remains unrealized due to the large scale of the problem and difficulties in designing effective similarity loss functions.
In this study, we develop a coherent landmark mapping (CLM) method, which integrates deep learning for reliable cross-modal whole-brain registration. With CLM as the core local registration module, we built a comprehensive registration pipeline, named mBrainAligner, to facilitate the whole-brain mapping of images, as well as their digitally reconstructed neuronal structures. Different from the 'detection-matching-screening' strategy that was often used in landmark-based registration methods, CLM uses a 'detection-mapping' strategy to systematically tackle challenges of cross-modal registration, such as problematic landmark detection across modalities, nontrivial identification and screening of unreliable matches, and inefficient utilization of shape priors of brain regions that are embedded in atlases. In addition, we also used a DNN to boost CLM further in terms of registration accuracy and robustness. We show that mBrainAligner can tackle the challenging registration case of fMOST data to the Allen brain atlas with improved accuracy compared to existing pipelines, and provides a foundation for current large-scale single-neuron analysis efforts 32 . We have also validated the generalization ability and the performance of mBrainAligner using other modalities including LSFM, VISoR and magnetic resonance imaging (MRI). We also report an fMOST-space mouse brain atlas and showcase the utility of our method for analyzing single cell types and registering partially imaged brain. results mBrianAligner algorithm and workflow. We designed mBrain-Aligner to align three-dimensional (3D) mouse brain images of different modalities to the CCFv3 (Fig. 1). Considering the need, importance and technical challenges of registering rapidly accumulating fMOST data to the CCFv3, we used fMOST image alignment as a specific example to describe our general cross-modal registration strategy. mBrainAligner also enables mapping digitally reconstructed dendritic, axonal and soma distributions to the target space (Fig. 1f), thus mBrainAligner faciliates the direct visualization, comparison and analysis of various neuronal structures at different resolutions, such as at the population, single-neuron or subcellular level.
Brain images of different modalities often vary substantially in their voxel intensity, image texture and brain anatomy (for example, due to uneven brain shrinkage) compared to the average template of the CCFv3 acquired using STPT (Fig. 1a). To address this challenge, we designed three modules in mBrainAligner: (1) an image preprocessing and robust point matching-based automatic global registration module (Fig. 1b); (2) a CLM-based automatic local registration module (Fig. 1c); and (3) an optional semiautomatic refinement module ( Fig. 1d and Extended Data Fig. 1). Using mBrainAligner, we also translated the CCFv3 to a 'normalized' fMOST space to produce a fMOST-domain atlas (Fig. 1g).
In an initial step, images are preprocessed to match in size and overall intensity profile. In the fMOST example, all images are first downsampled anisotropically (XYZ multipliers are 64 × 64 × 16) roughly to match the size of the target image (CCFv3 average template, 25 µm), followed by contrast enhancement and noise filtering to remove stripe artifacts that were generated during the knife cutting and imaging processes ( Fig. 1b and Supplementary Video 1). Next, we extract about 1,500 brain contour points in the subject and target brain, estimate an affine transform that maximizes the similarity of the brains' outer contour, and register the global orientation and scale of these brains ( Fig. 1b and Supplementary Video 2).
A key advance in mBrainAligner is the robust and accurate nonlinear local registration that warps a subject image to the target via a dense point-by-point correspondence between two images. As the cross-modal variations in intensity, texture and shape of globally registered brains often remain severe (Fig. 1b), we used a CLM algorithm to coherently deform the target landmark points obtained via two and a half dimensional (2.5D) corner detector to fit the subject brain ( Fig. 1c and Supplementary Video 3), so that each landmark and its neighbors in the target image find their best matches in a coordinated manner. As landmarks before and after deformation naturally form matching pairs, CLM intrinsically avoids error-prone landmark detection and identification of unreliable matches across modalities. To ensure a robust and coherent deformation of landmarks, the deformation process is driven by three complementary features (gradient, histogram of oriented gradient and segmentation probability) that encode the edge, shape context and semantic information, respectively, while being regularized by the shape priors of brain anatomy embedded in the CCFv3 (or the fMOST atlas). We used a DNN to estimate the segmentation probability of six main brain regions, including hypothalamus (HY), caudoputamen (CP), hippocampal formation (HPF), cerebral cortex (CTX), cerebellar cortex (CBX) and brain stem (BS) (Extended Data Fig. 2). As brain delineation or segmentation is a major application of atlas-based registration, we show that the registration task can also benefit from segmentation by incorporating the semantic information as a discriminative feature (Extended Data Fig. 2c,d). Using CLM, the joint efforts of a coherent mapping strategy, discriminative feature integration and effective shape prior utilization ensure the robustness and accuracy of mBrainAligner in challenging cross-modal scenarios. The computation time for CLM-based local registration was about 6 hours per brain with 2214 target landmarks on a 2.3 GHz dual-10-core Xeon workstation with 32 G RAM.
Registration accuracy. We first evaluated mBrainAligner's performance in registering the border of six brain regions of interest, that is, HY, HPF, CTX, CBX, BS (thalamus (TH), midbrain (MB) and hindbrain (HB)) and CP, which three neuroanatomists annotated independently (Fig. 2). The registration quality of mBrainAligner was visually evident for test images from different modalities including fMOST, LSFM, VISoR and MRI (Fig. 2a). For 31 fMOST whole mouse brain images, the region-wise median Dice score achieved by mBrainAligner's local alignment ranges from 0.82 to 0.90, while the average median Dice score for all brain regions of interest reaches 0.88 (Fig. 2b). We compared mBrainAligner with the mouse brain registration pipelines MIRACL 42 , ClearMap 41 and aMAP 40 . In our tests, mBrainAligner outperformed these methods in all six regions, and achieved 0.11 higher average median Dice score than the second-ranked method (for example, 0.77 for aMAP) (Fig. 2b). We also compared the performance of mBrainAligner and alternative methods in aligning LSFM, VISoR and MRI modality images to the CCFv3 (Extended Data Fig. 3), where we observed consistent superiority of mBrainAligner.
The anatomical location of somas is often used as a key feature to assign neuron types 45 . We considered 1,708 neurons that were reconstructed from 31 fMOST mouse brains and that had all somas manually annotated (Fig. 3a). These somas distribute mainly in the cortex and thalamus. Five brain regions, ventral posteromedial nucleus (VPM), CP, ventral posterolateral nucleus (VPL), dorsal part of the lateral geniculate complex (LGd) and medial geniculate complex (MG), contain more than 40 somas. We computationally identified the anatomical locations of these somas by mapping them to the CCFv3. Compared with the manually annotated results, mBrainAligner achieves an average soma mapping accuracy of 76% (Fig. 3b), which is at least 19% higher than other methods. A similar conclusion can be drawn by taking a closer look at the results of individual brain regions (Fig. 3c). For those somas whose locations mismatch between automatic registration and manual annotation, we calculated their shortest distance to the manually assigned region as distance error, and found mBrainAligner achieved an average distance error of 47.88 µm for all somas which is half the value achieved with the second-ranked method (Fig. 3d). Visualization of registered somas in the VPM (in thalamus) and primary somatosensory area nose (SSp-n) (in cortex) shows that mBrainAligner produces results more consistent with annotators' expectation ( Fig. 3e). We analyzed the distribution of wrongly localized somas ( Fig. 3f and Extended Data Fig. 4) and observed that most errors occurred in localizing VPL as VPM, which are two adjacent small regions without a clear discriminative boundary between them in fMOST imaging (Fig. 3g).
Different thalamic regions were hypothesized to have unique cortical projection patterns 46 . We examined this by computationally predicting the axon projections of 63 neurons from six different brains with somas in VPM or VPL. We found their axons terminated predominantly in layers 4 and 6 of the cortex and formed clear patterns (Fig. 3h). Compared with manual annotation by anatomists, mBrainAligner predicted 89.1% of the projection target regions when we used automatic local registration ( Fig. 3i and Extended Data Fig. 5). The accuracy was further improved to 96.29% after we applied semiautomatic registration (Extended Data Fig. 1e). For 162 neurons with apical dendrites that distributed in the L1 layer of the cortex (Fig. 3j), the prediction accuracy is more than 90%.
fMOST-domain atlas. The Allen CCFv3 is an increasingly used standard in mouse brain research. This atlas was developed using iteratively averaged STPT TissueCyte imaging. In a similar manner we used mBrainAligner to build a fMOST-domain atlas of mouse brain (Fig. 4a,b and Extended Data Fig. 6), including an average template and an annotation template (Fig. 4b). The fMOST atlas offers both biological and engineering advantages. First, we identified the morphology variation between two modalities in a quantitative way (Fig. 4c), which helps estimate the deformation caused by sample preparation (for example, different clearing methods) in various methods. Second, we established a reusable and invariant transformation between two atlases (Fig. 4b), thus brain annotation labels and reconstructed results provided in each modality can be effectively shared. Third, we reduced the technical difficulty to map images directly across these two modalities (Fig. 4d). The robust and accurate registration of mBrainAligner makes it possible to generate a high-quality average template using a small number of samples, that is 36 fMOST whole mouse brains in our case. This number is less than that used in building the CCFv3 (1675 STPT whole mouse brains) and the LSFM atlas (139 LSFM imaged brains). We built the fMOST annotation template by backward-mapping the annotation of the CCFv3 using the same warping parameters for image registration (Fig. 4b). In this way, we ensured consistent definition of brain anatomical regions in the two atlases, and facilitated the sharing of information across modalities. We also asked independent neuroanatomists to validate further the biological meaningfulness of the fMOST atlas, especially in the motor cortex, a region of recent focus in the BRAIN Initiative's consortium effort 47 .
The fMOST-domain atlas also affords improved registration accuracy. Using the same metrics described above, we re-evaluated the performance of mBrainAligner and alternative methods to register fMOST images to the fMOST atlas. Because the most notable differences of cross-modal anatomy and intensity are no longer present, all metrics, including the Dice score of brain delineation, soma localization accuracy, distance errors of mapped somas and axon projection prediction, show improvements especially for MIRACL, ClearMap and aMAP (Fig. 4d). mBrainAligner still outperforms the other methods in all metrics, and also outperforms its results achieved with the original CCFv3. With the fMOST atlas, mBrainAligner further reduces the average distance error of mapped somas from 47.88 µm to 29.88 µm, and improves the average median Dice score of brain regions to around 0.90. These results not only verified the quality and utility of the fMOST atlas, but also demonstrated the effectiveness of mBrainAligner in both intra and cross-modal registration.
Cell type analysis enabled by mBrainAligner. We used mBrain-Aligner in cell type analysis and also studied neuronal projection  pathways at single-neuron resolution. The topography of thalamic neurons was characterized previously at the neuron population level 48 . Here we used mBrainAligner to map neurons reconstructed via our BICCN consortium effort 32 to the CCFv3. We identified 1268 cortical axonal arbor clusters that projected to the major cortical areas (for example, primary somatosensory area barrel field/ mouth/nose (SSp-bfd/m/n) for VPM neurons) and studied their correlation with soma locations (Fig. 5a and Extended Data Fig. 7). Following the convention of the CCFv3, larger values of x, y and z dimensions represent more posterior, ventral and medial locations in the 3D brain atlas.
As an illustrative analysis, we created a two-dimensional (2D) top view of the Common Coordinate Framework (CCF) cortical structures by virtual flattening and quantified the axonal arbor locations within the 2D plane. We found that axonal arbors of VPM neurons exhibited a striped distribution, with the stripe approximately aligning with the x axis of the 2D map. For VPM neurons, the x coordinates (in the 2D cortex) of the axon centers are positively correlated with the x coordinates (in 3D CCF space) and anticorrelated with y and z coordinates of somata (Fig. 5b). The dominant pattern (P = 6.5×10 −79 , two-sided Fisher's exact test) of VPM topography is that neurons originating in the dorsal VPM project to the posterior part of the primary somatosensory cortex, namely the barrel field (SSp-bfd). This pathway is shorter than the pathway between ventral-VPM and SSp-m (Fig. 5c). The z coordinate of the axon is anticorrelated with the y and z coordinates of soma (Fig. 5b). We observed the same relationship between axon and soma locations of VPL neurons, except that their axon z axis is anticorrelated with soma x axis (P = 0.025, two-sided Fisher's exact test, Extended Data Fig. 7). For intergeniculate leaflet of the lateral geniculate complex (LGd) neurons, the x and z coordinates of the axon centers are anticorrelated with the x and y coordinates of somata (Extended Data Fig. 7). Among the four sensory relay nuclei, MG neurons show the lowest level of topography where only the x axes of soma and axon were correlated (P = 0.015, two-sided Fisher's exact test, Extended Data Fig. 7). For the nonrelay thalamic nuclei, especially high-order nuclei including ventral medial nucleus of the thalamus (VM) and nucleus of reuniens (RE), we integrated neurons with different nuclei origins and performed an overall topography test (Fig. 5a,d). The overall relationship of soma and axon locations is highly consistent with that of VPM and VPL, especially the anterior-posterior arrangement of axonal arbors. From anterior-ventral-medial nuclei, such as parataenial nucleus (PT) and rhomboid nucleus (RH), to posterior-dorsal-lateral nuclei, such as intergeniculate leaflet of the lateral geniculate complex (IGL), axonal arbors follow an anterior-posterior distribution pattern (Fig. 5e).

Discussion
Motivated in part by the ongoing large-scale collaboration of the BRAIN Initiative's BICCN program to characterize brain cell types across multiple modalities, we developed mBrainAligner, a registration pipeline to tackle challenges of cross-modal registration. A whole-brain map of full neuronal morphology in three dimensions generated in this way helped test a number of hypotheses about neuronal projection diversity at multiple organizational levels 32 .
With mBrainAligner, we generated a mouse brain atlas in the fMOST domain. As a large amount of fMOST images are being generated, and 'it is essential that brain atlases keep pace with large-scale data collection' 30 , an atlas in the fMOST domain is desirable. However, the atlas building process is traditionally computationally demanding and labor intensive. mBrainAligner enabled us to build a high-quality 'standard' average fMOST-based template from a small sample set without manual intervention. By precisely registering the fMOST average template to the CCFv3 with mBrainAligner and then back-projecting the annotation of the CCFv3 onto fMOST, a smooth transformation of annotations and identical parcellation of brain regions between two atlases are guaranteed. This atlas building pipeline is generic and can be applied to other modalities. Finally, for cross-modal brain registration, how to quantify the registration accuracy is often debated 49 . As soma locations, axon projection and neuronal morphology were used as defining features in the task of cell typing 50 , in a biologically meaningful way, we show that mBrainAligner produces improved results compared to alternative approaches. In addition to mBrainAligner's use cases for intra and cross-modal registration of whole mouse brain, this tool can also be used to register partially imaged brains. We have open-sourced mBrainAligner, as well as some useful tools including a 2.5D corner detector, stripe artifacts removal, and image or metadata warping tools, to promote interdisciplinary and reproducible researches.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41592-021-01334-w.

Methods
Animal care and use. All experiments related to the use of mice followed National Institutes of Health guidelines, and received approval from the Institutional Animal Care and Use Committee of the Allen Institute for Brain Science, the University of Science and Technology of China, the Icahn School of Medicine at Mount Sinai and the Center for Genomic Regulation, respectively. All mice used in this study were group-housed with a 12-hour light/dark cycle, controlled temperature (20-22 °C) and humidity (50-70%), and free access to food and water.
Data preparation. The average and annotation template of the CCFv3 were downloaded from Allen Institute web portal (http://atlas.brain-map.org/). The detailed sample preparation, imaging and atlas generation methods have been described previously 30 . As the foreground brain area in the CCFv3 average template is close to image boundaries, to eliminate computational ambiguity, we padded the CCFv3 average template (25 µm) by 20 pixels at both ends of the anterior-posterior axes and used the padded image as the target. Our fMOST dataset consisted of 57 brains of male and female transgenic mice (P56), with sparsely labeled cells in cortical, thalamic, claustral and striatal regions, and for cholinergic, noradrenergic and serotonergic neuronal types. The detailed sample preparation and imaging methods were described previously 32 . The auto-fluorescence channel was used for image registration. We used the TeraConverter module of Vaa3D 51 to perform image downsampling. The neurons and somas of 31 fMOST brains were reconstructed via the TeraVR module 20 of Vaa3D.
The MRI dataset was downloaded from https://www.nitrc.org/projects/ incfwhsmouse, which contains image volumes representing the canonical Waxholm Space adult C57BL/6J mouse brain. This dataset provides five images, including T1, T2* and T2-weighted magnetic resonance volumes, Nissl-stained optical histology and a volume of labels. All volumes are represented at 21.5 µm isotropic resolution. The labeled volume represents phase one delineations of the Waxholm Space mouse brain. The atlas contains 26 structures (+ inner ear) in accordance with the NeuroLex Brain Partonomy at: https://scicrunch.org/ scicrunch/interlex/dashboard. We chose the T2*-weighted MRI image dataset as test data, and removed brain tissues that did not match the CCFv3 based on the provided labels.
The VISoR dataset was generated from an 8-week-old C57BL/6J mouse. The brain sample was sectioned into 300 μm slices, cleared and then stained with 4,6-diamidino-2-phenylindole to label the nuclei of all cells. All slices were imaged with a custom VISoR microscope for high-throughput volumetric imaging at 1 × 1 × 2.5 μm 3 voxel resolution, and stitched together to reconstruct the whole mouse brain volume. The dataset was downsampled to 4 × 4 × 4 μm 3 resolution before registration.
The LSFM partial brain hippocampus dataset was generated from a 12-week-old Thy1-yellow fluorescent protein mouse (strain B6.Cg-Tg(Thy1-YFPH)2Jrs/J no. 003782; The Jackson Laboratories). The mice were killed and transcardially perfused with cold phosphate buffered saline (PBS), followed by cold CLARITY hydrogel solution. The brains were removed from the skull, immersed in 20 ml of hydrogel solution and stored at 4 °C shielded from light. After 2-3 days at 4 °C, hydrogel polymerization was induced by 3 h incubation at 37 °C. Contact with oxygen was prevented with a vegetable oil layer. Thereafter, the already polymerized surrounding hydrogel was removed. Entire brains were immersed in 50 ml of clearing solution and stored in a water bath at 45 °C. After 2 months of incubation in clearing solution for passive clearing, optically transparent brains were washed with 50 ml PBST (0.1% TritonX in 1× PBS) twice for 24 hours each, at 45 °C. After that, the samples were immersed in a refractive index matching (75 g diatrizoic acid, 70 g d-sorbitol and 23 g N-methyl-d-glucamine in 100 ml of water) medium for 2 days before imaging. After incubation in refractive index matching the samples were imaged in a custom-built single plane illumination microscopy. We used a detection lens with magnification 5×, a laser wavelength of 488 for fluorescence excitation and a 5 μm z step, obtaining a 1.29 × 1.29 × 5 μm 3 voxel size.
Removal of stripe artifacts. The stripe noise present in the raw fMOST images is mainly caused by fluorescent bleaching during the knife cutting and imaging process. We modeled this kind of stripe artifacts as multiplicative periodic noise and designed a log-space frequency notch filter to realize high-quality stripe artifacts removal (Supplementary Video 1). For each fMOST image, (1) one 2D coronal slice with the largest foreground brain area was extracted (Extended Data Fig. 8a) and converted to the frequency domain using fast Fourier transform (Extended Data Fig. 8b). (2) By examining the frequency spectrum, the bandwidth, orientation and cutoff frequency of the notch filter were identified and then used to construct a Gaussian notch filter (Extended Data Fig. 8c). (3) We imposed the log transform on the intensity of the fMOST image in the spatial domain to convert the multiplicative noise to additive one (Extended Data Fig. 8d). (4) The log-transformed fMOST image was filtered in the frequency domain in a slice-by-slice manner using the Gaussian notch filter generated in step 2 (Extended Data Fig. 8e). (5) Finally, we computed the inverse fast Fourier transform on the filtered spectrum and used inverse-log transform to bring the filtered image back from the log space (Extended Data Fig. 8f).

Robust point matching-based automatic global registration.
To minimize the impact of intensity variation between different modalities, we implemented the automatic global alignment by affine aligning brains' outer contour points through point-set matching. First, the foreground brain region was extracted via Otsu thresholding 53 , followed by morphology filtering to remove the small and thin artifacts. Then we uniformly sampled (7 × 7 × 7 grid) the contour of the foreground region to obtain a point set (~1,500 points per brain) that represents the shape of the brain's outer contour. To tackle the large direction variation between brains, the principal axes of target and subject point sets were extracted using principal component analysis and then rigidly aligned. Finally, following the deterministic annealing framework of robust point matching 54 , we iteratively refined an affine transform to maximize the overall shape of two point sets (Supplementary Video 2). We set the initial annealing temperature to 0.25, and gradually decreased the temperature with an annealing rate of 0.96 during iteration.

2.5D corner detection.
As CLM establishes the correspondence between two brains by deforming the landmarks defined in the target brain to fit the subject brain, the quantity and quality of target landmarks will have a big impact on the overall performance of registration. We designed a 2.5D corner detector automatically to extract points of high curvature (corners or junction of different brain regions) of a 3D image and applied it on the annotation template of the CCFv3 to generate anatomically meaningful target landmarks. We realized the 2.5D corner detector via a two-level detection and screening strategy.
In the first level, given a 3D image, we started by extracting all its coronal, sagittal and horizontal planes. For each plane, the Harris response 55 of each pixel was computed, thresholded (discarded if response is less than 0.1), and nonmaximum suppressed (with window size 7 × 7) to find the local maxima as candidate corners. We denoted these candidate corners as the first-level corners, and saved their Harris responses at coronal, sagittal and horizontal planes, respectively, for the second-level screening.
In the second level, we screened the first-level corners by only retaining ones that were detected as first-level corners in at least two orthogonal planes. We called the remained corners second-level ones and updated their Harris response by accumulating their responses at corresponding orthogonal directions. The final landmarks were obtained by applying 3D nonmaximum suppression (with windows size 31 × 31 × 31) on the second-level corners.
CLM-based automatic local registration. Coherent landmark mapping. CLM establishes a dense correspondence between two images by coherently deforming target landmarks to fit the subject brain (Supplementary Video 3). The landmarks before and after deforming naturally form matching pairs, and were used to calculate the displacement field for local image warping.
Suppose we have a subject image that consists of N voxels V = {v i , i=1, 2,…, N }, and a target landmark-set of M points L = {l j , j=1, 2,…, M}, where v i is the ith voxel and l j is the jth landmark in subject and target image, respectively. Let P = { pij, pij ∈ [0, 1] } be the classification probability of ith voxel to the jth landmark, we form the CLM as the following energy minimization problem: where f(l j , Θ) is a nonlinear spatial transform which maps the landmark l j to a new location with parameters Θ, si,j ∈ [0, 1] denotes the matching score between voxel v i and currently deformed landmark f(l j ,Θ), and r is the distance influence factor that controls the impact of voxel-to-landmark distance to their classification probability. Minimizing the first term encourages assigning larger classification probability to the voxel-to-landmark pair with higher matching score and smaller Euclidian distance. While minimizing the second term is equivalent to maximizing the entropy of P, which aims to enforce the convexity of the energy function and alleviate the risk of algorithm getting trapped in the local minimum during optimization. T acts as a balance factor between two terms. We adopted a dual update strategy to solve the joint classification and spatial transform optimization problem posed in energy function (1). In each iteration, we first fix the transformation Θ and calculate the voxel-to-landmark classification P with currently deformed landmarks, then update the transform and perform landmark deformation with classification probability fixed. We denote the above two steps as landmark-guided voxel classification and voxel classification-guided landmark deformation, respectively.
Landmark-guided voxel classification. Given spatial transform Θ, the deformed target landmarks can be obtained as f(L,Θ). For any i,j value, minimizing E(P,Θ) with respect to P will lead to: we chose T = 1, and normalize classification probability so that ∑ i pij = 1 for any i ∈ [0, N] once they are updated using equation (2). To speed up the iterative voxel classification, we only concerned a subdomain (21 × 21 × 21) around each landmark and chose r = 10.
Voxel classification-guided landmark deformation. With fixed voxel classification probability P, we realized the coherent landmark deforming in three steps. First, we preliminarily updated each landmark to the weighted center of mass of all voxels l pre j = ∑ N i=1 pijvi. Then we updated transform parameters Θ so that deformed landmarks f(L, Θ) can best approximate their preliminary updated ones { l pre j , j = 1, 2, …, M } under smoothness and brain shape prior constraints.
We chose smooth-thin-plate-spline (STPS) 56 to model the spatial transform f with parameters Θ, which can be decomposed into an affine transform and a smoothness parameter λ regularized nonlinear warping component. If λ = 0, an exact interpolation was performed. We chose λ = 0.001 to balance the goodness of fit between two point sets and the smoothness of deformation. The spatial transform parameters Θ were solved as a least-squares solution of STPS. Finally, we used LittleQuickWarp 57 (a fast and memory efficient implementation of STPS) with updated Θ to transform the initial target landmarks and realize the coherent landmark deforming. As target landmarks are obtained on the CCFv3 annotation template, we ensured that all landmarks are on the boundaries of the brain region. They naturally capture the shape prior of these regions and form a simplified point-set representation of the the CCFv3. The smooth (regularized by STPS) and coherent deformation of such a simplified atlas onto the subject brain implicitly utilizes the shape priors and spatial relationship of brain regions, while these priors are hard to incorporate in the commonly used B-spline based methods [37][38][39][40][41][42] . If provided, the modality-specific shape priors can also be accommodated under the CLM framework to enhance further the robustness and accuracy of registration. Given the fMOST atlas, as the transformation between the CCFv3 and fMOST atlas was known and fixed, we utilized these modality-specific shape priors by first mapping the initial landmarks of different brain regions to fMOST space, and then deforming them in a region-wise manner following the above-mentioned STPS deforming strategy. In our implementation, this region-wise constraint was applied before the brain-wise constraint in each iteration.
Matching score computation. We based the calculation of matching score between a voxel-to-landmark pair upon three complementary features: (1) gradient, (2) histogram of oriented gradients (HOG) 58 , and (3) segmentation probability. As these features encode the edge, shape context and semantic information, respectively, their integration is expected to provide CLM with more reliable matching metric to cope better with the challenges of cross-modal registration.
We generated the gradient feature by applying the Sobel filter (3 × 3 × 3) on Gaussian smoothed (3 × 3 × 3) volume. At a given location, its HOG feature was calculated by concatenating the normalized gradient histograms of eight blocks within the patch (21 × 21 × 21) that centered on that location. We divided each patch into 27 cells (7 × 7 × 7) and combined its neighboring (2 × 2 × 2) cells to create blocks. For each block, the gradient histogram (nine orientations) of each cell within it was computed, concatenated and normalized to form the block feature. This resulted in a 576-dimensional HOG feature for each location.
We chose to generate the segmentation probability feature using a semantic segmentation network. As a proof of principle, we adopted a slightly modified 3D U-Net 59 to generate the segmentation probability (with probability values range from 0 to 1) of each voxel to six main brain regions (HY, HPF, CTX, CBX, BS, CP) and background. Indeed, under the CLM framework, 3D U-Net can be readily substituted with other more sophisticated semantic segmentation networks to improve the registration performance further. The network architecture used in this study is shown in Extended Data Fig. 2a, which contains a classic encoder-decoder structure. We used four layers at the encoder stage and four layers at the decoder stage, and they have 32, 64, 128 and 256 channels, respectively. The seven-dimensional score maps generated at the final softmax layer were used as the segmentation probability feature. The network was trained from scratch on 57 brains with ground truth segmentation labels obtained via mBrainAligner's semiautomatic registration module. We used 23 brains as the training set, three brains as the validation set and the remainder for testing. The segmentation results of six brain regions are visualized in Extended Data Fig. 2b. Potential segmentation errors can be mitigated as the segmentation information is only one of the three discriminative features that guide the deformation of landmarks, for which the deformation is further constrained by the overall shape priors.
Finally, we computed the matching score of each landmark-to-voxel pair by accumulating the similarity score of their centered patches (21 × 21 × 21) on three feature maps. In our implementation, (1) inverse mean squared error, (2) mutual information and (3) normalized cross-correlation were used as similarity metrics.
Semiautomatic refinement. The semiautomatic registration module provides an iteractive interface for investigators to visualize, examine and fine-tune the registration results (Fig. 1d, Extended Data Fig. 1 and Supplementary Video 4). In this interface, the boundary of brain regions defined in the CCFv3 was vectorized as Bezier curves and visualized as a layer overlaid on the registered brain. The control nodes of Bezier curves were evenly sampled from brain region boundaries defined in the CCFv3 annotation template. The registration quality could be examined by checking how well the registered brain would fit the CCFv3 boundaries in coronal, sagittal and horizontal views at different scales. Different to the reported method that performed slice-by-slice adjustment in coronal view 60 , our method could fine-tune three-view panels synchronously by simply dragging the curves, thus distortions in Z direction could also be rectified. In addition, we also provided a magnetic lasso 61 mode to assist in the alignment of curves in a semiautomatic manner. Three levels of subregion granularity (coarse, medium, fine) were provided for each main brain region to accommodate different local granularity requirements (Extended Data Fig. 1b). We realized the final 3D image deforming by using STPS transformation that computed on the offsets of fine-tuned control nodes. We demonstrated that the semiautomatic registration module could register partially imaged brain (LSFM partially imaged brain of hippocampus in Extended Data Fig. 9), a known challenging registration problem for fully automatic methods.
3D fMOST atlas building of mouse brain. We built the fMOST average template from 36 mouse brains in an iterative way. Before entering the iteration, all fMOST brains were preprocessed and globally aligned to the CCFv3 following the methods described above. To encourage a symmetric atlas, the globally aligned brains were mirrored across the midline, which resulted in a total of 72 brains. We iterated the following two steps until convergence: (1) each brain was locally registered to the template (we chose a high-quality globally aligned fMOST brain as the template in the first iteration) by using the automatic local registration module of mBrainAligner and averaged together; (2) the average deformation field over all brains was computed, inverted, and used to perform shape normalization on the average image obtained in (1) to compensate the potential shape bias towards CCF. This shape-normalized average brain was then used as the anatomical template in step (1) for the next iteration. We stopped the iteration until the sum of the magnitude of the average deformation field was below the threshold of 1.0.
We generated the fMOST annotation template by first aligning the fMOST average template to the CCFv3 and then back-projecting the annotation template of the CCFv3 using the inverse of the deformation field generated. To ensure a highly accurate registration, we registered the fMOST average template to the CCFv3 by concatenating the automatic and semiautomatic registration module of mBrainAligner.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The full resolution fMOST image datasets (https://download.brainimagelibrary. org/biccn/zeng/luo/fMOST/) of all mouse brains used in this study, as well as the original and CCFv3 registered single neuron reconstructions (https:// doi.org/10.35077/g. 25), are available at BICCN's Brain Image Library (BIL) at Pittsburgh Supercomputing Center (www.brainimagelibrary.org). The single-neuron reconstructions, the CCFv3 registered version of these reconstructions, as well as 3D navigation movie-gallery of these data are also available at SEU-ALLEN Joint Center, Institute for Brain and Intelligence (https:// braintell.org/projects/fullmorpho/). The Allen CCFv3 is available at (http://atlas. brain-map.org/), and MRI mouse brian used in this study can be obtained through (https://www.nitrc.org/projects/incfwhsmouse). We also provided downsampled fMOST, LSFM, VISoR and MRI whole mouse brain data and LSFM partially imaged brain data along with code and testing scripts at Github page (https:// github.com/Vaa3D/vaa3d_tools/tree/master/hackathon/mBrainAligner). The full resolution LSFM and VISoR brain images are available on request. Source data are provided with this paper. Fig. 3 | Comparison of different methods in registering images of different modalities to CCFv3. Four green dashed boxes show the brain registration results of four different modalities (LSFM, VISoR, MRI, fMOST) generated by different methods (MIRACL, ClearMap, aMAP and mBrainAligner,). The first and second rows of each box show the coronal and sagittal view of registered brains respectively with zoom-in view of two subregions. Registered brains of different modalities with brain region boundaries (orange) defined in CCFv3 overlaid. In each box, arrows of the same color point to the same spatial coordinates.