## MarmoNet: a pipeline for automated projection mapping of the common marmoset brain from whole-brain serial two-photon tomography

Brain Image Analysis Unit RIKEN Center for Brain Science, Japan http://bia.riken.jp

Henrik Skibbe, Akiya Watakabe, Ken Nakae, Carlos Enrique Gutierrez, Hiromichi Tsukada, Junichi Hata, Takashi Kawase, Rui Gong, Alexander Woodward, Kenji Doya, Hideyuki Okano, Tetsuo Yamamori, Shin Ishii

Abstract -Understanding the connectivity in the brain is an important prerequisite for understanding how the brain processes information. In the Brain/MINDS project, a connectivity study on marmoset brains uses two-photon microscopy fluorescence images of axonal projections to collect the neuron connectivity from defined brain regions at the mesoscopic scale. The processing of the images requires the detection and segmentation of the axonal tracer signal. The objective is to detect as much tracer signal as possible while not misclassifying other background structures as the signal. This can be challenging because of imaging noise, a cluttered image background, distortions or varying image contrast cause problems.

brain research projects [1], [4], [5]. Brain imaging data from primates such as macaque or marmoset monkeys can bridge the gap between the brain of simple vertebrate models, such as rodents, and the complexity of the human brain.

We are developing MarmoNet , a pipeline that processes and analyzes tracer image data of the common marmoset brain. The pipeline incorporates state-of-the-art machine learning techniques based on artificial convolutional neural networks (CNN) and image registration techniques to extract and map all relevant information in a robust manner. The pipeline processes new images in a fully automated way.

This report introduces the current state of the tracer signal analysis part of the pipeline.

## I. INTRODUCTION

A major challenge of contemporary neuroscience research projects is the analysis of brain image data to better understand the brain structure and function [1], [2], [3]. In particular, the primate brain is in focus of major national and international

- H. Skibbe is with the Brain Image Analysis Unit, RIKEN Center for Brain Science, Wako, Japan. (e-mail: henrik.skibbe@riken.jp).
- J. Hata and H. Okano are with the Laboratory for Marmoset Neural Architecture, RIKEN Center for Brain Science, Wako, Japan and the Department of Physiology, Keio University School of Medicine, Tokyo, Japan.
- C. E. Gutierrez, H. Tsukada and K. Doya are with the Neural Computation Unit, Okinawa Institute of Science and Technology, Okinawa, Japan.
- A. Watakabe and T. Yamamori are with the Laboratory for Molecular Analysis of Higher Brain Function, RIKEN Center for Brain Science, Wako, Japan.
- K. Nakae and S. Ishii are with the Department of Systems Science, Kyoto University, Kyoto, Japan.
- R. Gong and A.Woodward with the Connectome Analysis Unit, RIKEN Center for Brain Science, Wako, Japan.
- T. Kawase is with the Janelia Farm Research Campus in Ashburn, Virginia, USA.

This research is supported by the program for Brain Mapping by Integrated Neurotechnologies for Disease Studies (Brain/MINDS) from Japan Agency for Medical Research and Development, AMED

A structural map of the brain is an indispensable prerequisite for studies of the neural networks controlling higher brain functions. The Brain/MINDS project in Japan conducts multimodal studies to construct an integrated structural map of the marmoset brain [4], [6]. An ambitious goal of the project is to gain new insights into information processing and diseases of the primate brain. As part of the Brain/MINDS project, we conduct a systematic tracer mapping using viral-based anterograde tracers.

A large amount of tracer image data demands automated processing and analysis. To cope with this demand, we are developing MarmoNet, an image processing pipeline that automatically processes and analyses imaging data of the marmoset brain.

The goal in our anterograde tracer study is to track the fluorescent signal of the axons from their cell bodies to their point of termination. Figure 1 illustrates the concept of the anterograde tracer technique that we use. To trace the axon projections of neurons located in a specific region of the brain, a virus bearing a gene of a fluorescent protein is injected into the region's center; see Fig. 1(a). After viral infection, the neurons in that brain region become filled with the fluorescent protein. All fluorescent signal outside the infected area is considered to come from axons that had originated from cells located at the injection site; see Fig. 1(b). By connecting the injection site location with the axon tracer projection sites, we can infer brain connectivity. By systematic injections of many tracers, we can construct a map of neural connections across the brain regions of interest; see Fig. 1(c).

In recent years, several pipelines for the processing of neural anterograde tracer image data have been proposed [7], [8], [9]. Typically, hand-crafted heuristic features detect the tracer signal, and thresholds are being applied to segment the images into background and tracer signal. In [7] for example, edge enhancing filters have been applied to enhance the signal of axon structures. Then, morphological attributes have been

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-72efb483-0f61-47ef-9b92-fa5c1ddc9f0a.jpg)

(a) An adeno-associated virus bearing a gene of fluorescent protein has been injected to a defined region in the left frontal cortex to tag axonal projections. Following viral infection, the cells become filled with the fluorescent protein.

(b) All fluorescent signal outside the infected area is considered to come from axons that had originated from cells located at the injection site (red area in the image). By connecting the injection site location with the axon tracer projection sites, we can make assumptions about brain connectivity.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-19b69cd7-1f97-46a1-b275-08e0be8999cf.jpg)

(c) Injection locations in the prefrontal cortex. By systematic injections of many tracers, we can construct a map of neural connections across the brain regions of interest.We exemplarily show the pattern of 3D axonal projections corresponding to four different injection sites.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-3e53a943-689c-4493-9588-3726cae4994d.jpg)

Fig. 1. Anterograde tracer injection procedure for tagging axonal projections.

extracted and classified to further improve the segmentation results. In [8], an independent component analysis was used to separate the tracer signal from the background.

It is challenging to describe all kinds of patterns of axon tracer signal within the images. The axons in the images appear in various shapes and with varying contrasts and intensities. Some axons, particularly when sparsely sampled and running within the imaging plane, appear as elongated, bright structures. However, some other densely sampled axons build dense, blob-like patterns with a texture that can vary depending on whether they cross and bifurcate or they share a dominant direction as part of a thicker axon bundle. Hence it is oftentimes challenging to manually craft a feature extraction pipeline that can cope with all cases equally well so that no further manual interaction is necessary.

MarmoNet uses U-Net [10] to deal with these problems. U-Net is an artificial convolutional neural network (CNN) architecture that integrates both the feature extraction and the classification into one common trainable framework. U-Net has shown excellent performance in various segmentation tasks of neurite structures from biomedical image data [11], [12] and therefore was a nearby candidate for solving our tracer segmentation problem.

An advantage of integrating machine learning into the pipeline is that the feature extraction and classification can be constructed in a data-driven way. This is done by providing a set of representative examples to train an U-Net to solve the task. A single example is a pair of an input image and its corresponding segmentation result. The machine then learns to extract the features and learns the classification parameters

Fig. 2. MarmoNet processing pipeline. In this report, we focus on the tracer segmentation parts (b) and (c).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-9e6b0038-59e8-4ca7-b7af-c32c91867476.jpg)

that are best suited to generate the segmentation for an input image.

As the time of writing, the construction of our marmoset connectivity map is an ongoing project. We frequently presented our progress at conferences during the last years [13], [14], [15], [16], [17], [18], [19], [20], [21]. Figure 2 shows an outline of the entire tracer data processing and analysis pipeline, including the pre-processing and the post-processing that generates the output. The pipeline automatically detects the location of the injection site (connections origin), and the axonal projections (connection targets). The pipeline maps all image data into a common image reference space which allows the integration of connectivity data from different individuals. Figure 3 shows the result for 36 examples.

This report focuses on the tracer signal segmentation part including the training and application of the U-Net to process and analyse the data. The tracer segmentation can be divided into two major steps (red box around those steps in Fig. 2). The first step is the automatic determination of the location and boundary of the injection site (b), and the segmentation of the axon tracer signal outside the injection site (c).

## II. IMAGE ACQUISITION

Individual brain images have been acquired with one of our TissueCyte 1000 and TissueCyte 1100 (TissueVision, Cambridge, MA) serial two-photon tomographs (two-photon microscopes), which can generate optical section images from the block face surface of agarose-embedded biological material in coordination with a computer-driven microtome and stage movement [22]. We target the neurons with a mixture of clover and 1/4 amount of presynapse targeting mTFP1 . Figure 4 (a) illustrates the emitted wavelengths of the fluorescent signal as an overlay of a bar diagram representing the three filters of the tissue microscope. The bar diagram represents the wavelengths

Fig. 3. 3D image stacks of the detected injection site (first column, superior view) and the segmented axon tracer signal (2nd column superior view and 3rd column lateral view) for 36 injection sites. All images have been mapped to the Brain/MINDS reference space.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-15abacd9-3cbb-42c2-950b-6fb7f200c260.jpg)

captured by the three color channels: red ( C R ), green ( C G ) and blue ( C B ).

- (a) Applying an intensity correction that reduces inhomogeneous illumination effects and stitching the images.

Figure 4 (c) illustrates the key features of all three channels. The autofluorescent signal is similarly strong in channels C R and C G . The tracer signal, however, is much more stronger in C G than in channel C R . This makes channel C R an ideal candidate for background signal removal. In the injection site, we have the best contrast in channel C B . According to their characteristics, we will call C R the background channel (for background removal), C G the tracer channel (axonal tracer signal best visible), and C B the injection site channel (best image of cell bodies).

## III. MARMONET PIPELINE

The pipeline comprises three major components. The preprocessing part that performs the image stitching, the actual image analysis part, and a pipeline output part where connectivity calculations are being carried out. Figure 2 illustrates the workflow of the pipeline. The single steps are:

- (b) A CNN detects the location of the injection site.
- (c) A CNN segments the axon tracer signal outside the injection site from the background.
- (d) An image registration algorithm maps all image data to the reference image space of the Brain/MINDS marmoset brain atlas.
- (e) Connectivity calculations are being carried out.

All pipeline steps will be described in the remainder of this section, the tracer segmentation parts will be described in more detail.

## A. Pre-Processing

(a) Image stitching: The TissueCyte microscope acquires an image stack of the entire marmoset brain in a fully automated manner. Starting with the prefrontal cortex, the brain is imaged every 50 µm from its anterior part to its posterior regions (front to back).

Fig. 4. Emitted fluorescent tracer signal and its application for identifying the cell body signal, axon tracer signal and the auto-fluorescent background signal. The channels C R , C G and C B are representing the three imaging channels of the TissuCyte microscope.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-34894996-309d-478f-ae38-cf6ee81344a6.jpg)

The TissueCyte microscope scans the entire surface of the block-face in a zigzag pattern using a movable table while generating a new 2D image tile I j k every 750 µm , where j is the tile index, and k denotes one of the three color channels C R , C G or C B . The size of each image tile has been set to 720 × 720 pixels with a spatial in-plane resolution of 1 34 . × 1 34 . µm 2 .

Intensity Correction : The illumination of single image tiles is inhomogeneous due to distortions in the optical path. Typically, the intensity decreases in a nonlinear manner and is most bright in the image center, while getting darker towards the image boundaries. This visual effect is called vignetting . The vignetting effect differs from microscope to microscopes such as for our TC1000 and our TC1100. The different pathways of the various color channels cause different effects as well. The distortion is often modeled as a linear modulation

$$I _ { k } ^ { j } ( x, y ) = \overline { I _ { k } ^ { j } } ( x, y ) \hat { F } _ { k } ( x, y ) ; \quad \ \ ( 1 ) \ \stackrel { \lambda \ \dots } { [ 2 5 ] }.$$

see e.g. [23], [24]. The I j k are the observed image tiles, I j k the undistorted image, and ˆ F k the shading field , which is the image of the vignetting effect in channel k . We assume that the ˆ F k are constant throughout an entire marmoset brain image stack.

The technique to correct this distortion is called flat-field correction. The goal is to estimate the shading field ˆ F k so that we can recover the original image tile I j k by inverting (1).

Similar to [25], we use a data-driven heuristic to estimate the shading field ˆ F k by averaging over all image tiles of an entire marmoset brain image stack. Averaging over a sufficiently large number of image tiles I j k yields an estimate of the shading field F k (up to a constant factor). In contrast to [25], we exclude outliers from the average. Outliers are dark pixels with no information and bright pixels, most likely, can be attributed to tracer signal. We describe the details of our flatfield correction in section A in the appendix. Figure 5 and Fig. 21 of the appendix are show results of our flat-field correction.

It is worth noting that we only perform an intensity correction for tile images of the channels C R (background) and C G (tracer). Channel C B (Injection site) has a clear signal around the injection site which is an insufficiently small amount of imaging data for estimating the shading field.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-fc99540f-2da1-4a01-846e-3b4067ad50e1.jpg)

(a) before

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-8519c691-644b-4f1f-8785-31743b4ad038.jpg)

(b) after

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-0caafe2b-e506-45f8-ac67-21431fb767ec.jpg)

(d)

after

Fig. 5. The flatfield correction noticeably reduces the vignetting effect.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-473a2471-8aee-45bc-9e20-e3ab346a9251.jpg)

Image Stitching : The TissueCyte microscope provides the offset for each image tile in 3D world coordinates with micrometer resolution. The microscope generates optical section images from the block face surface before a built-in vibrating blade microtome mechanically cuts off the most upper tissue section. There is no spatial distortion within an image slice. Hence an entire image slice can be reconstructed by aligning and fusing all image tiles I j k according to their world coordinates. In our settings, there is a small overlap between adjacent image tiles (about 80 pixels). We crop 50 pixels from the image tile boundaries and fuse the remaining overlapping parts using a linear blending function (similar to [25]).

Full resolution image stacks : We reconstruct the entire brain image stack in full resolution. Each slice of the image stack that corresponds to a physical brain section has a spatial resolution of 1 34 . × 1 34 . µm 2 pixels. The spatial distance between two image slices is 50 µm . We denote the entire 3D image stacks of channel C k with a spatial resolution of 1 34 . × 1 34 . × 50 µm 3 by I high k . Each image stack is a function

$$\mathcal { I } _ { k } ^ { \text{high} } \colon \Omega \rightarrow R _ { \geq 0 },$$

where Ω defines the image domain. We address single voxels in I high k via I high k ( x, y, z ) , and entire 2D sections with I high k ( z ) . The index k defines the channel and z the section number that corresponds to a physical section in the TissueCyte microscopy data.

Low resolution image stacks : In addition, we scale the full resolution image stacks I high k down to lower resolution image stacks I low k with an isotropic voxel resolution of 50 3 µm 3 . The isotropic image stacks are used for image registration, injection site localization and connectivity calulations.

## B. Image Analysis

The extraction of brain connectivity information from the tracer data demands the segmentation of the tracer signal from

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-05f79b10-3815-465f-931f-dcd6046c89b2.jpg)

Fig. 6. Example of an estimated shading field.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-035c0dfa-bd50-4a90-9b61-8b9448575a06.jpg)

injection site detection in the 3D image stack

cell detection in high-resolution 2D image slices

Fig. 7. Injection site localization pipeline.

the image background. The pipeline detects the tracer signal in two steps. In the first step, it detects the signal of the cell bodies in the injection site. In the second step, it segments the axonal tracer signal in the remaining parts of the brain image. After tracer segmentation, the image stacks are mapped to a common population average brain reference image space to integrate brain connectivity data from different individuals.

(b) Injection site localization: The injection site location is determined in two steps. First, a heuristic approach roughly estimates the volume of the injection site region in the lowresolution image of the injection site channel C B . Further analysis is then limited to that volume to reduce computation costs. The second step identifies the location of the infected cell bodies in the injection site. A CNN determines the cell body locations in the high-resolution image section of the injection site. The cell body positions define a 3D point cloud that implicitly determines the exact injection site boundary. The procedure is illustrated in Fig. 7.

- 1) Rough Localization: This step is illustrated in Fig. 7 (a). The injection site appears significantly brighter than the remaining tissue regions in channel C B . We use this feature to roughly estimate the injection site location and boundary. In this step, all calculations are performed on the low-resolution 3D image stack I low C B . We first classify the image into background voxels and injection site candidate voxels using a constant intensity threshold. As threshold we determined 4500 in a heuristic manner. The width of the Gaussian is σ = 150 µm (3

voxels). Typically, most of the voxels in the injection site have been classified as tracer signal so that after convolution, the injection site appears as a homogeneous, bright structure in the image. We denote this image by ˆ S low . Figure 7 (a.ii) shows an example. We obtain the final injection site mask S low by thresholding ˆ S low a second time. The second threshold is determined dynamically. It has been set to t low = 0 5 max( ˆ . S low ) , half the maximum value of ˆ S low . The thresholded image is defined by

$$S ^ { \text{low} } ( x, y, z ) = \begin{cases} 1 & \text{if } \hat { S } ^ { \text{low} } ( x, y, z ) > t _ { \text{low} } \\ 0 & \text{else} \end{cases}. \quad ( 3 )$$

There can be more than one connected regions in S low . We used a connected component algorithm to remove all but the largest connected region in S low .

Figure 7 (a) shows the maximum intensity projection of a raw image I low C B of channel C B , the maximum intensity projection of ˆ S low , and a maximum intensity projection of the final mask S low .

- 2) Precise Localization: This step is illustrated in Fig. 7 (b). We use the rough localization of the injection site in the low-resolution image I low C B to determine the corresponding region of interest in the full resolution image stack I high C B . Then we determine the precise location of infected cells in the full resolution 2D image sections I high C B ( z ) . Sub image (i) in Fig. 7 (b) shows an example image.

The pipeline uses a CNN to detect the locations of the cell bodies in the 2D image section I high C B ( z ) . The network is applied to all 2D slices that intersect with the segmentation mask S .

The pipeline uses a U-Net [10] for cell detection, a widely used CNN architecture for image segmentation tasks. The network has been trained to mimic a function F S that takes an image slice I high C B ( z ) as input and maps it to a cell saliency image ˆ S high ( z ) : Ω → [0 , 1] . The extents of ˆ S high ( z ) are identical to I high C B ( z ) . The values in ˆ S high ( z ) correspond to the likeliness of cell body locations. As brighter a pixel in ˆ S high ( z ) , the more likely there is a cell at the corresponding position in the image slice I high C B ( z ) . Figure 7 (b.ii) shows an example. In a post-processing, we only keep the positions of local maxima in the image plane of ˆ S high ( z ) that exceed a given threshold t high &gt; 0 . The final set S high of detected cell positions is defined by the 3D point cloud

$$S ^ { \text{high} } \coloneqq & \{ \forall ( x, y, z ) | \\ & \hat { S } ^ { \text{high} } ( x, y, z ) \text{ is local maxima in } \hat { S } ^ { \text{high} } ( z ) \\ & \wedge \hat { S } ^ { \text{high} } ( x, y, z ) > t _ { \text{high} } \},$$

where the threshold t high is equal to 0.5. Figure 7 (b.iii) shows the detected cell positions in the corresponding image slice.

Details of the network architecture and the training are explained in the next section of this report.

Fig. 8. Segmented tracer signal from image slices of five different marmoset subjects acquired with different microscopes.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-0a82fef7-4a84-4f5a-85a2-c3bfe0761a43.jpg)

(c) Axon tracer signal detection: The detection of axons can be quite challenging. The axons in the images appear in various shapes and with varying contrasts and intensities. Some axons, particularly when sparsely sampled and running within the imaging plane, appear as elongated, bright structures. Some other densely sampled axons build dense, bloblike patterns with a texture that can vary depending on whether they cross and bifurcate or they share a dominant direction as part of a thicker axon bundle.

Further, in some images, there are various non-axonal structures, like the boundary of some blood vessels, that, locally, appear similar to weak axon tracer signal in the channels C R and C G , which makes them even harder to distinguish.

We experimented with intensity thresholding or edge detection approaches. While we could achieve good results for many individual image slices, it was impractical to find good parameters that work well for all images. The choice for parameters was always a trade-off between too many false positive detections or too many false negative detections. Section IV-2 shows some examples for the thresholding approach. We finally made use of a machine learning approach that can cope with the variation in the data.

Similar to the cell detection task in the previous section, we use a CNN based on the U-net architecture to perform the tracer signal segmentation. The network architecture is mostly identical to the network that we use for cell detection. Only the image input dimension is different. Figure 9 shows an outline of the axon tracer segmentation pipeline.

The trained network is a function F L that takes both an image I high C R ( z ) of the background channel C R and an image

Fig. 9. Axon tracer segmentation pipeline.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-7a3f17c4-920c-45af-9c6f-a8c030ae3f54.jpg)

I high C G ( z ) of the tracer channel C G as inputs (Fig. 9 (a)), and maps them to an axon tracer signal saliency image ˆ L high ( z ) : Ω → [0 , 1] ; see Fig. 9 (b). As brighter a value in ˆ L high ( z ) , the more likely there is an axon at the corresponding position in the image slices I high k ( z ) .

The intensity correlates strongly with the axon density. This is an important feature that we believe is correlated with the connection strength between the injection site and the projection targets of the axons. To reflect this important feature in the signal segmentation, we first compute the signal strength by subtracting the background channel from the tracer signal channel (removing the autofluorescence signal), where

$$T \{ t \} ( x, x, z ) = \begin{cases} 0 \text{ if } \mathcal { I } _ { C _ { G } } ^ { \text{high} } ( x, y, z ) < t \mathcal { I } _ { C _ { R } } ^ { \text{high} } ( x, y. z ) \\ \mathcal { I } _ { C _ { G } } ^ { \text{high} } ( x, y, z ) ) - t \mathcal { I } _ { C _ { R } } ^ { \text{high} } ( x, y, z ) \text{ otherwise} \end{cases}$$

is the tracer signal after background subtraction. The threshold t has been set to 1 1 . . We then threshold the saliency map ˆ ( L z ) , which is the output of the neural network, to obtain a mask that removes all remaining non-axonal structures. Finally, we weight the tracer signal with the raw tracer signal T z ( ) according to

$$L ( x, y, z ) \coloneqq ( \ddot { L } ( x, y, z ) > 0. 5 ) T \{ 1. 1 \} ( x, y, z ). \quad ( 5 )$$

The image L is the segmented tracer signal. The scaling factor 1 1 . has been determined experimentaly by comparing tracerless background image tiles of channel C R and C G . Figure 9 (c) shows an example slice. Other examples can be found in Fig. 8.

(d) Mapping to reference space: After pre-processing and signal analysis, the image data is mapped to the Brain/MINDS reference image space. The mapping to a common reference space is a required prerequisite for the comparison and integration of imaging data from different individuals. In the Brain/MINDS project, a refined version of the Brain/MINDS marmoset brain atlas [26] defines the common space that connects the imaging data of different modalities and data from different laboratories. The refined version of Brain/MINDS marmoset brain atlas defines 786 brain regions per hemisphere. The atlas is in alignment with a brain reference image space based on a T2 weighted MRI population average brain image. The spatial resolution of the T2 weighted MRI image is 100 3 µm 3 .

All images, including the raw images, but also the injection site and tracer labels are in the same image space as the images

Fig. 10. The TissueCyte population average serves as a proxy to map individual brain images to the Brain/MINDS reference image space.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-660f1598-0570-4d87-a5a8-152fcc9abd85.jpg)

of channel C R . In a first step, we determine the mapping between the background image C R and the T2 image of the reference brain using the ANTs registration toolkit [27]. ANTs is that provides the mapping from the individual image space to the reference image space and its inverse mapping. Once we have this mapping, we can apply this mapping to all remaining images to map them from the individual image space to the reference space and vice versa.

In a preliminary version of the pipeline, we directly mapped the 3D image stacks of the background channel C R to the reference brain image. There are two drawbacks with such approach. First, the spatial resolution of the target image is lower than the 50 3 µm 3 resolution of the isotropic TissuCyte image stacks. Further, the modalities differ. Both make an exact mapping challenging. While working fine for many, some images frequently made problems.

We decided to put a TissuCyte average brain image as a proxy in between the TissuCyte images of an individual marmoset and the Brain/MINDS marmoset brain atlas. Figure 10 illustrates the procedure. Due to the same contrast and same spatial resolution, the registration between two TissuCyte images works better than the registration between two different modalities. In this new setting, the more challenging registration between the Brain/MINDS marmoset brain atlas and the TissuCyte average image has to be determined only once. This has the advantage that even semi-automated or even manual approaches can be applied to improve the registration accuracy between the TissueCyte population average and the T2 MRI reference brain image.

We created a TissuCyte average image with a spatial resolution of 50 3 µm 3 based on images of channel C R from 34 different marmoset monkey brains. We have decided to make the TissuCyte population brain image axisymmetric. The main reason for that was that the number of samples was, at the time the Brain/MINDS program started, even smaller than the 34 images. Compared to the Allen mouse brain atlas template [7], [2] a rather small number. An axisymmetric population average image doubles the number of samples per voxel, which significantly improves the quality of the result. Further, it removes any bias in favor of the left or right hemisphere. Figure 11 shows an ortho view of the TissuCyte population average brain image.

2.45 cm

3.61 cm

Fig. 11. The axisymmetric population average image defines the TissueCyte reference image space.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-e0d79a9b-6004-4aa9-94be-130c7f3775a1.jpg)

## C. Post-Processing

(e) Brain connectivity calculations: After determining injection site location, tracer signal and identifying individual brain regions by mapping the data to the reference brain image space, we can determine the neural connectivity associated with the injection site.

Therefore, we first identify the connection source area by determining the brain gray matter region(s) that intersect with the injection site. Then, we identify the projection targets by determining gray matter regions that intersect with the axon tracer signal. The amount of signal that intersects with the brain regions is used to determine the connection strength between two brain regions.

Each brain provides complementary information about marmoset brain connectivity. A structural map of marmoset brain connectivity can be obtained by integrating the information of a large number of individual marmoset brains.

## IV. TRAINING THE U-NET IN MARMONET

For cell body detection and tracer signal segmentation, in both problems we seek a function F that maps entire images X : Ω in → R to output images Y : Ω out → [0 , 1] with identical domains Ω in = Ω out . Each point in Y ( x, y ) is a label, or a value that is associated with the image content at location X x,y ( ) . In case of the cell body detection the outputs Y are the cell body saliency maps ˆ S high ( z ) , in case of the axon tracer signal the functions ˆ L high ( z ) .

It is oftentimes challenging, or even impossible to explicitly define the function F . However, it is frequently possible to define a sufficiently large set of samples X i from the input domain, and their corresponding samples ˆ Y i from the output domain with ˆ Y i = F ( X i ) . Once such a set is given, it is often

possible to sufficiently approximate the function F by fitting a trainable function F{ } ξ to the data, where ξ are the free parameters that need to be determined. The set of pairs of samples ( X , Y i ˆ ) i is called a training set, because it is used to train F{ } ξ to approximate F .

In our case, the function F{ } ξ is an artificial neural network based on the U-Net network architecture. U-Net is an artificial neural network architecture popular for a broad variety of image segmentation tasks [11], [10], [28]. U-Net takes images X : Ω in → R as input and generates output images Y : Ω out → [0 , 1] .

The creation of the training set for the two different tasks, cell detection and axon signal segmentation, are discussed in more detail in the subsections IV-1 and IV-2 below. Further technical details regarding the U-Net architecture in this pipeline can be found in section B of the appendix.

(a) Network Training: Injection site localization: For cell body detection, we have trained an artificial neural network F{ } ξ that maps full resolution image slices I high C B ( z ) of channel C B to cell body saliency maps ˆ S high ( z ) .

We created a training set by manually labeling the cell body locations of 6068 cells in 44 2D image slices that intersect with the injection site. The images have been selected from 10 different marmoset brain images. The training images have been divided into smaller, overlapping image tiles that were sufficiently small to be fed into the neural network for training.

Figure 13 shows the manual labels of such a training image tile. The first column shows a training input image X . The middle column the ground truth labels ˆ Y . It is a binary image, where only a single pixel at the center of a cell body has been marked positive (=1). All the background has been set to zero.

It is difficult, error-prone, and often ambitious to manually find and label the exact center of a cell. Further, some bright axons may look similar to cell bodies, and we would like to put an extra penalty to incorrect predictions on those structures during the training. Therefore we automatically generated an additional, spatially dependent weight image W ˆ Y that scales the error during the network training based on the image content. The most right column in Fig. 13 shows the corresponding example.

The default weight in W ˆ Y was set to one. In order to tolerate a small displacement between the position of a manually placed label and the position predicted by the network, we masked out any prediction error in a small circular surrounding around a cell body label by setting its weight to zero. To ensure a clear boundary between labels and background, we set the weight of the circle boundaries to two. Further, there are much more negative examples (background pixels) than positive examples (cell body centers). We took this into account by assigning a high weight to the cell body labels (500). Finally, we used a Laplacian of Gaussian filter to detect edge and blob-like structures, thresholded the results and took the result to add an additional weight (=2) to structures, like axons, that may share some similarities with the cell bodies. We determined all parameters heuristically.

The samples in the training set may not sufficiently represent all variations in the shape and size of unseen images.

input image gamma 0.8

gamma 1.2

max warp, granularity 1

max warp, granularity 0.1

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-00293c43-6442-4410-b20c-aed74ffe45c1.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-59c7f31d-c091-4e67-bc60-c25f39da33bc.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-7369aebf-3220-42d8-b09c-65d2fc618a1f.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-1803296f-046d-4e17-83ec-8bb004a89bc3.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-95b7af64-5b74-4111-a536-4932aeb4261d.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-3e42db42-ddf1-4335-ac2c-5678d7beb2f5.jpg)

max warp, granularity 0.01

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-1fbd2b84-8bfa-4527-9202-e772d30cae09.jpg)

rotation

Fig. 12. Augmentation of training image data.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-989ae35b-b3eb-4284-9e85-3d8460d75cbe.jpg)

After training with an insufficient number of samples, the trained network may work well on the training set but then often performs poorly on unseen images. In this situation, data augmentation can help. Data augmentation increases variation in the training set by randomly altering the appearance of the existing images in the training set. In our training, the deformations are one or more combinations of rotations, gamma corrections, smooth, non-linear spatial deformations and global scaling. Figure 12 shows examples. An example of detected cells after applying the network to an unseen image is shown in Fig 15.

After training with the augmented data, the network was able to segment the axins in all of our 36 different marmosets brain images; see Figure 3 for the results.

Fig. 13. An exampled of an image for training a network for cell body detection

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-38b426f4-9ec8-46c0-8634-61d2b101a1de.jpg)

EVALUATION: We have created a test set by manually labeling 3524 cells in five unseen images. Further, we have implemented a hand-crafted cell detection approach and put its detection performance in relation to the neural network approach. Figure 16 shows the precision-recall graphs for both approaches when applied to the training set, and when applied to the test set. It clearly shows the superior performance of the neural network approach. Figure 14 shows a qualitative example.

For the neural network approach, we took the network predictions and detected the local maxima as described in section III-B2. Then we computed precision and recall by varying the threshold t high in (4).

Fig. 14. Cell body detection using a Hessian filter and the proposed artificial neural network approach.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-1e69fd0b-f006-475d-900d-dcc9e05d8f87.jpg)

Fig. 15. The most left figure shows a slice image of the tracer channel C G . As an overlay, the cell density is shown as a heat map at the injection site location. The images on the right show the detected cells.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-889e3e68-a2e5-40eb-ac07-34b0480e7a71.jpg)

As an alternative to the neural network, we implemented a blob detection filter capable to detect the roundish shape of the cell bodies within an image. We have implemented a multi-scale filter that analyzes the local curvature of the image based on the second order image derivatives of the Gaussian filtered image. We first applied an isotropic Gaussian filter to the image with standard deviation σ . Then we computed the Hessian matrix field for the entire image domain and computed the first and second eigenvalues of the matrices. We obtained two new images λ 1 and λ 2 of eigenvalues, where, without loss

(c) Performance on the traing dataset (6068 labeled cells in 44 images) and test dataset (3524 labeled cells in 5 images)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-96b5a119-67db-4f73-9c52-735a84f322c2.jpg)

Fig. 16. Cell detection performance of the U-Net in comparison with a handcrafted cell detectioon filter.

of generality, | λ 1 ( x, y, σ ) | &gt; λ | 2 ( x, y, σ ) | . With

$$\mathcal { H } ( \sigma ) = - \lambda _ { 1 } ( \sigma ) \frac { | \lambda _ { 2 } ( \sigma ) | } { | \lambda _ { 1 } ( \sigma ) | + \epsilon } ( \lambda _ { 2 } ( \sigma ) < 0 )$$

we define the cell detection filter response H ( σ ) : Ω → R , glyph[epsilon1] &gt; 0 is a small constant. The eigenvalues λ 1 ( x, y, σ ) and λ 2 ( x, y, σ ) of the Hessian matrix at location ( x, y ) are proportional to the two local main curvatures within the image plain.

In case of a bright, roundish plateau, like a cell body, both values are negative. As larger their magnitude, as higher the curvature, as more likely it is a cell. Furthermore, as more roundish the shape is, as more similar should be the two eigenvalues. We took this into account by weighting the term with -λ 1 ( σ ) (maximum curvature → cell contrast), and weight it with the fraction | λ 2 ( σ ) | | λ 1 ( σ ) + | glyph[epsilon1] (ratio between minimum and maximum curvature → cell shape score). The last term ( λ 2 ( σ ) &lt; 0) removes all saddle points such as edges of a cell bodies or axons.

We have evaluated the performance for all possible combinations of scales (Gaussian filter size) out of the interval σ ∈ { 2 , .., 5 } on the training set. In case of multiple scales, we took the maximum filter response over all scales. We achieved the best performance with a single scale of σ = 4 .

Figure 16(a) shows the precision-recall graph for the Hessian filter for σ = 4 and the neural network approach on the training set. Figure 16(b) shows the precision-recall graph on the test set. We have marked the threshold that performed best (precision equals recall) during training in the test graph indicating that both approaches can maintain their performance on the test set.

(b) Network Training: Axon tracer signal detection: For axon tracer signal detection, we have trained an artificial neural network F{ } ξ that takes two full resolution image slices I high C R ( z ) and I high C G ( z ) of channel C R and channel C G as input and maps them to the axon tracer saliency map ˆ L high ( z ) .

In contrast to the cell body detection task, it was impractical to manually annotate a sufficiently large amount of images for training the network for the axon tracer detection. The axons appear in various shapes and with varying contrasts, intensities and densities. Some axons appear as elongated, bright line-like structures with clear boundaries where a manual annotation may be feasible. However, many others build dense, blob-like or cluttered patterns with a texture that can vary depending on whether they cross and bifurcate or they share a dominant direction as part of a thicker axon bundle.

Instead of manual annotation, we used an approach based on intensity thresholding to obtain a large set of labeled brain image sections and afterwards manually selected the best candidates for training. We selected about 600 image sections from 20 different marmoset brain image stacks. Even within the 600 image slices, various non-axonal structures, like the boundary of some blood vessels, lead to falsepositive detections. After selecting the 600 best candidates, we manually annotated such structures in the images to provide some explicit negative examples to the network. Figures 17

DENSE

SPARSE

Fig. 17. The threshold based approach worked well for many images. We manually selected the best candidates for training the neural network.

and 18 show examples of the thresholded tracer signal, and manually annotated background structures, respectively.

From the 600 labeled full-brain image slices, we automatically generated a training set containing about 12000 smaller image tiles: We smoothed the labels (the thresholded tracer signal) with a large isotropic Gaussian function to obtain a tracer signal density map. We used the density map to generate two probabilistic density functions to randomly pick 20 smaller image tiles from each of the 600 full resolution image sections. The first function was designed to draw 10 tiles from axon dense regions while the second one was used to draw 10 tiles from axon sparse regions. About one-third of the full brain slices did not contain any tracer signal. When there was no tracer signal in a full brain slice then all 20 image tiles that were generated from that slice were only background signal (about 1/3 of the data). Further, all structures with manually labeled negative examples have been divided into smaller tiles and have been added to the training set. The size of an image tile was 810 × 810 pixels which then has been cropped, after augmentation, down to the network input size of 572 × 572 pixels.

The ratio of background to tracer signal in the training data was about eight to one. To balance the different classes, we

C1 (Background)

C2 (Tracer)

Tracer labels (RED)

Fig. 18. Bright structures, like blood vessels, frequently lead to false positive detections. We manually labeled many of them as negative examples for training the U-Net.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-dd81a422-0634-45d4-af22-d0584a8b8e6e.jpg)

weighted the error of incorrectly predicted tracer labels by eight. The weights for manually annotated negative examples have been set to 100.

THRESHOLDING: The axon tracer signal is most dominant in the green color channel C G ; see Fig. 4. While it is significantly stronger in C G than in C R , the autofluorescent background signal is mostly identical in both channel (up to a small constant factor).

Subtracting channel C R from C G increases the contrast between the tracer signal and the background; see Fig. 19 (b). With equation (5), we can apply a threshold to that high-contrast image to obtain a tracer signal segmentation mask. However, determining a good threshold that separates the background from the tracer signal is tricky. The varying brightness and contrast of the axon tracer signal within different brain regions make the determination of a good threshold value difficult. It is always a trade-off between a low falsenegative rate or a low false-positive rate. Automated methods like Otsu's method [29] did not perform well.

We obtain good results by using a double-threshold approach. A high threshold separates the bright tracer signal

Fig. 19. Morphological pipeline for creating axonal tracer training data

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-6e21cafc-1037-4339-bc87-06b12107c57a.jpg)

from the background. The high threshold has been chosen conservatively so that we were confident that no background structures remains in the image. A second, lower threshold has been determined that recovers as much tracer signal as possible while keeping the false positive detections low. We did a grid search over threshold parameters and have qualitatively compared the results by visual inspection. We determined 300 for the high threshold, and 100 for the low threshold. In Fig. 19 (c) we show examples.

The lower threshold retrieves the weak tracer signal thas has been classified as background by the higher threshold. However, it incorrectly classifies many other small structures as tracer signal as well. In the next step, we use morphological reconstruction to select those connected regions out of the lowthreshold image that have - at least - one pixel in common with the high threshold image. The idea is that the high threshold image selects the reliable regions in the low threshold image. The first image in Fig. 19 (d) shows an example. In a final step, we use morphological closing with a disk-shaped structure element of size three to smoothen the result.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-e22cc5a7-3d2c-4548-9d2d-86f4c5aa8c24.jpg)

recall: 0.94

recall: 1

Fig. 20. Axon tracer signal detection.

EVALUATION: We exemplary selected two unseen images with axon tracer signal. One image with densely labeled axons that share the same main direction as part of a thicker axon bundle in the corpus callosum. The other image from branching axons with varying density from within the brain gray matter. For both images, we subtracted channel C R from channel C G and manually tuned the image intensities for both images to increase the contrast between the axon signal and the background signal. Then, we applied a threshold to obtain a segmentation mask and manually corrected errors. We took the manually segmented image as ground truth and compared it with the two automated approaches, the thresholded output of the neural network, and the output of the double-threshold approach. Figure 19 shows the results. The last column illustrates the difficulty of finding a threshold parameter that works well for all images. While performing well for the image of the axon bundle, it is too sensitive for the other image leading to many false positive detections. The neural network approach has a clear advantage. Data augmentation provides an automated method to prepares the network for unseen images with a contrast or brightness that is different to the images in the training data.

## V. DISCUSSION &amp; CONCLUSION

We are developing a pipeline that processes, analyzes and maps tracer image data to a common marmoset brain image space. The pipeline incorporates state-of-the-art machine learning techniques based on artificial convolutional neural networks (CNN) and state-of-the-art image registration techniques to extract and map all relevant information in a robust manner. The pipeline processes images in a fully automated way.

In Fig. 3 we show results of the tracer segmentation and the corresponding cell densities in the injection site. All images have been mapped to the Brain/MINDS marmoset brain atlas reference image space.

In this report, we briefly introduced all steps of the tracer signal processing pipeline. The two major components of the image analysis part of the pipeline have been explained in more detail, which are: the injection site detection and the axon tracer signal segmentation.

In future versions of the image analysis pipeline, we plan to improve the axon tracer signal segmentation. Axon tracer signal segmentation is difficult due to the broad variety of shape and contrast of the axon signal pattern, and the density of labels. Completely manually creating a sufficiently large ground truth database for training a CNN remains impractical. However, some kind of interactive relevance feedback system might be appropriated, where an artificial neural network is capable to automatically detect ambiguous classification results and gets feedback from human experts to improve its performance. With such an approach we may be able to create a smaller and better performing training set with reasonable effort and time.

It is often impractical, error-prone, and/or strongly biased to analyze a large amount of data in a manual manner so that automated approaches are strongly demanded. Of course, the lack of correct ground truth data makes it difficult to perfectly master all those problems with automation. Nevertheless, automated processing is often the only option to gain new insights from large scale imagining data. We believe that the best practice to evaluate the data and the results is to share the results with a broader community of experts by making them publicly available. We plan to make our results publicly available in the near future.

## APPENDIX

## A. Flat-field correction

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-c9780f58-f086-421b-b43c-dcd49ae9f5d3.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-4f052f0b-0f54-4aa4-8364-31df135c4822.jpg)

(a) before

(b) after

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-4d284d55-c9a7-4e35-896e-ec8226a84668.jpg)

(c) before

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-97f448fd-faac-4210-a2c4-457f2b573ca2.jpg)

(d) after

Fig. 21. The flat-field correction noticeably reduces the vignetting effect.

The relation between the true image I and the observed image I is often modeled by

$$I = \overline { I } \cdot F + D ;$$

see for instance [23], [24]. The function I : Ω → R ≥ 0 is the image tile, D : Ω → R ≥ 0 the darkfield, and F : Ω → R ≥ 0 a multiplicative shading component which actually causes the vignette effect. With Ω we denote the image domain. The darkfield D represents the amount of signal given out by the detector when there is no incident light. A flat field can be corrected by inverting the model. The inverting is defined by

$$\overline { I } = \frac { I - D } { F }.$$

Removing the bias D is often necessary to calibrate CCD sensor based cameras and microscope systems. For our twophoton microscope, however, which scans the probe in a raster pattern pixel-by-pixel, it can be neglected. The flat-field correction then simplifies to

$$\bar { I } = \frac { I } { F }.$$

The vignetting effect slightly differs from microscope to microscopes such as for our TC1000 and our TC1100. The different pathways of the different color channels cause different effects as well.

Similar to [25], [30], [31], we use a data-driven heuristic to estimate the shading field F . When averaging over a sufficiently large number of image tiles I j we obtain an estimate of the shading field image F (up to a constant factor).

Input image

Output image

Fig. 22. The MarmoNet U-Net network architecture for cell body detection and tracer signal segmentation.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00876v1-figure-2cc033b2-582e-4498-955f-175df8accee8.jpg)

With (9), we have the relation

$$\frac { 1 } { M } \sum _ { i = 0 } ^ { M } I ^ { j } = F \frac { 1 } { M } \sum _ { i = 0 } ^ { M } \bar { I } ^ { j }. \quad \quad ( 1 0 ) \, \stackrel { \sqrt { \text{after} } } { \text{where} } \\ \text{be } c \, \text{c} \, \text{ }$$

We assume that all pixel intensities in the tiles I j are Poisson distributed in the same manner with mean value λ ∈ R &gt; 0 . The averaging over a large number of sample tiles I j leads to a flat mean intensity image with a constant value close to λ . In this case, equation (10) simplifies to

$$F \underbrace { \frac { 1 } { M } \sum _ { i = 0 } ^ { M } \overline { I } ^ { j } } _ { \approx \lambda } = F \cdot \lambda. \quad \quad \quad ( 1 1 ) \begin{array} { c } \text{How} \\ \text{basis} \\ \text{by tl} \\ \text{supe} \\ \text{outp} \end{array}$$

Let now

$$\hat { F } \coloneqq \frac { F \cdot \lambda | \Omega | } { \sum _ { \mathbf x \in \Omega } ( F ( \mathbf x ) \cdot \lambda ) } \quad \quad ( 1 2 ) \quad \text{func!}$$

be our normalized shading field approximation. Then we can approximate the flat-field corrected image I j with

$$\bar { I } ^ { j } \approx \frac { I ^ { j } } { \hat { F } }. \quad \quad \ \ ( 1 3 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

An entire marmoset brain image contains a large number of tiles. We maintain a low memory footprint by using a running average to approximate ˆ F . Further, we exclude outliers by discarding all dark pixels with an intensity lower than 2, and all bright pixels with an intensities greater than 2500 (axon tracer signal) during the average computation.

## B. U-Net

A U-Net consists of a feature extraction path and a voting path; see the left and right part of Fig. 22. With n we denote the number of channels of the input image. For the cell body detection problem, where only C B is fed into the network, it is n = 1 . For the tracer segmentation tasks, C R and C G are being concatenated so that the input dimension is n = 2 (two values per pixel).

The feature extraction path is a standard convolutional neural network. It consists of repeated applications of unpadded 3 × 3 convolutions followed by rectified linear units (light blue arrows) and 2 × 2 max-pooling (orange arrows). The free parameters are the convolution kernels. The objective of the network training is to determine the kernel weights that extract the image features from the input image that are best suited to solve the task.

Starting with 64 features after the first convolution layer, the number of features is doubled after each pooling operation to up to 512 features. Each pooling operation halves the extents of a feature image by taking the maximum over a four-pixel neighborhood. The pooling is an efficient way to enlarge the perceptive field of the neural network while keeping the computation costs on a reasonable level.

The voting consists of repeated applications of 2 × 2 up-convolutions (dark blue arrows) followed by applications of 3 × 3 convolutions and rectified linear units. Each upconvolution operation doubles the extent of a feature image. Starting with 512 features, we halve the number of features after each up-convolution down to 64.

While in the feature extraction path, the convolutions can be considered as an information extraction operation, the convolutions in the voting path can be considered as a drawing operation. During voting, each convolution kernel can be considered as a basis function that locally draws a specific pattern into the output image Y . The basis functions, determined by the kernel's free parameters, are tuned during the training. However, the weights that determine the contribution of each basis function in the image generation process are determined by the extracted features in the feature extraction path. The superposition of all those basis functions generates the final output.

The last activation function in the network is a sigmoid function that generates the label image Y by mapping all values to the interval [0 , 1] .

U-Net employs direct intermediate connections from feature layers to voting layers (gray connections in Fig. 22). The connections propagate spatial features directly from highresolution feature layers to deeper voting layers. Those features support a precise spatial localization of the output labels.

We further use drop out and batch normalization [32], [33] to improve the training performance; see Fig. 22, where we have labeled dropout layers with DrO and batch normalization layers with BaN .

The image patch input size of our network is 572 × 572 × n pixels. The output image is, due to the unpadded convolutions, cropped down to 484 × 484 pixels. Note that since all operations are convolutions, the image size can be changed during, or even after the training. For training, however, 572 × 572 × n was a good trade-off regarding GPU memory size and computation time.

When it comes to applications, an entire input image is rather large, and the amount of memory of a modern graphics card, which is used to perform the network processing, is limited. When applying the network to large image slices, a sliding window approach is used to move the receptive field of the network over the entire image to label its content one after another. This is done by dividing a large image into smaller image tiles, feed them into the network one after another, and stitch the network outputs together to form the entire output image.

For energy minimization during training, we used weighted logistic regression as a loss function.

A problem that appears quite often with biomedical image data is the small number of training samples. The number of free parameters in U-Net is large, and the samples in the training set may not sufficiently represent all variations in the shape and size of unseen images. As a result, often, the trained network works well on the training set but then performs poorly on unseen images. This is called overfitting. We increased the variation in the training set by altering the appearance of the existing images in the training set randomly. Figure 12 shows examples. The deformations are combinations of one or more transformations. We applied rotations, gamma corrections, smooth, non-linear spatial deformations and global scaling. The process is known as image augmentation.

## REFERENCES

- [1] H. Markram, 'The human brain project,' Scientific American , vol. 306, no. 6, pp. 50-55, 2012.
- [2] S. W. Oh, J. A. Harris, L. Ng, B. Winslow, N. Cain, S. Mihalas, Q. Wang, C. Lau, L. Kuan, A. M. Henry, et al. , 'A mesoscale connectome of the mouse brain,' Nature , vol. 508, no. 7495, pp. 207-214, 2014.
- [3] S. Grillner, N. Ip, C. Koch, W. Koroshetz, H. Okano, M. Polachek, M. Poo, and T. J. Sejnowski, 'Worldwide initiatives to advance brain research,' Nature neuroscience , vol. 19, no. 9, pp. 1118-1122, 2016.
- [4] H. Okano, E. Sasaki, T. Yamamori, A. Iriki, T. Shimogori, Y. Yamaguchi, K. Kasai, and A. Miyawaki, 'Brain/minds: A japanese national brain project for marmoset neuroscience,' Neuron , vol. 92, no. 3, pp. 582590, 2016.
- [5] M.-m. Poo, J.-l. Du, N. Y. Ip, Z.-Q. Xiong, B. Xu, and T. Tan, 'China brain project: basic neuroscience, brain diseases, and brain-inspired computing,' Neuron , vol. 92, no. 3, pp. 591-596, 2016.
- [6] H. Okano and P. Mitra, 'Brain-mapping projects using the common marmoset,' Neuroscience research , vol. 93, pp. 3-7, 2015.
- [7] L. Kuan, Y. Li, C. Lau, D. Feng, A. Bernard, S. M. Sunkin, H. Zeng, C. Dang, M. Hawrylycz, and L. Ng, 'Neuroinformatics of the allen mouse brain connectivity atlas,' Methods , vol. 73, pp. 4-17, 2015.
- [8] H. Abe, T. Tani, H. Mashiko, N. Kitamura, N. Miyakawa, K. Mimura, K. Sakai, W. Suzuki, T. Kurotani, H. Mizukami, et al. , '3d reconstruction of brain section images for creating axonal projection maps in marmosets,' Journal of neuroscience methods , vol. 286, pp. 102-113, 2017.
- [9] M. K. Lin, Y. S. Takahashi, B.-X. Huo, M. Hanada, J. Nagashima, J. Hata, A. S. Tolpygo, K. Ram, B. C. Lee, M. I. Miller, et al. , 'A high-throughput neurohistological pipeline for brain-wide mesoscale connectivity mapping of the common marmoset,' Elife , vol. 8, p. e40042, 2019.
- [10] O. Ronneberger, P. Fischer, and T. Brox, 'U-net: Convolutional networks for biomedical image segmentation,' in International Conference on Medical image computing and computer-assisted intervention , pp. 234241, Springer, 2015.
- [11] T. Falk, D. Mai, R. Bensch, Ö. Çiçek, A. Abdulkadir, Y. Marrakchi, A. Böhm, J. Deubner, Z. Jäckel, K. Seiwald, et al. , 'U-net: deep learning for cell counting, detection, and morphometry,' Nature methods , vol. 16, no. 1, p. 67, 2019.
- [12] Y. Kanaoka, H. Skibbe, Y. Hayashi, T. Uemura, and Y. Hattori, 'Determ: Software for automatic detection of neuronal dendritic branch terminals via an artificial neural network,' Genes to Cells , 2019.
- [13] A. Watakabe, J. Wang, M. Takaji, H. Mizukami, A. Woodward, T. Kawase, H. Skibbe, Y. Yamaguchi, S. Ishii, and T. Yamamori, 'Mapping connectivity of marmoset prefrontal cortex by serial twophoton tomography.' SFN, 2016. San Diego, USA.
- [14] A. Watakabe, J. Wang, M. Takaji, H. Mizukami, A. Woodward, T. Kawase, H. Skibbe, K. Nakae, Y. Yamaguchi, S. Ishii, and T. Yamamori, 'Columnar and diffuse connectivity of the marmoset pfc neurons.' SFN, 2017. Washington DC, USA.
- [15] A. Watakabe, J. Wang, M. Takaji, H. Mizukami, A. Woodward, T. Kawase, H. Skibbe, K. Naka, Y. Yamaguchi, S. Ishii, and T. Yamamori, 'Prefrontal projection mapping of the common marmoset.' SfN, 2018. San Diego, USA.
- [16] A. Watakabe, J. Wang, M. Takaji, H. Mizukami, A. Woodward, H. Skibbe, K. Nakae, Y. Yamaguchi, S. Ishii, and T. Yamamori, 'Prefrontal projection mapping of the marmoset.' Neuroscience, 2018. Kobe, Japan.
- [17] A. Watakabe, H. Skibbe, K. Nakae, J. Wang, M. Takaji, H. Mizukami, A. Woodward, R. Gong, Y. Yamaguchi, J. Hata, H. Okano, S. Ishii, and T. Yamamori, 'Mapping connectivity of the marmoset prefrontal cortex.' Neuroscience, 2019. Niigata, Japan.
- [18] H. Skibbe, K. Nakae, A. Watakabe, A. Woodward, C. E. Gutierrez, H. Tsukada, R. Gong, J. Hata, H. Okano, T. Yamamori, and S. Ishii, 'A pipeline for the automated integration of individual structural information from images of axonal tracer projections in marmoset brains into a structural map of full brain connectivity.' Neuroscience, 2018. Kobe, Japan.
- [19] H. Skibbe, N. Ken, W. Akiya, C. E. Gutierrez, A. Woodward, H. Tsukada, R. Gong, J. Hata, T. Yamamori, and S. Ishii, 'Marmu-net: Axon tracer signal detection using an artificial neural network.' INCF, 2018. Wako, Japan.
- [20] H. Skibbe, A. Watakabe, K. Nakae, C. E. Gutierrez, A. Woodward, H. Tsukada, R. Gong, J. Hata, H. Okano, T. Yamamori, and S. Ishii, 'A fully automated, ai-driven pipeline for the determination of the marmoset brain connectivity based on tracer data obtained with the tissuecyte microscope.' International Symposium of Brain/MINDS ISBM, 2019. Tokyo, Japan.
- [21] H. Skibbe, A. Watakabe, K. Nakae, C. E. Gutierrez, A. Woodward, H. Tsukada, R. Gong, J. Hata, K. Doya, H. Okano, T. Yamamori, and S. Ishii, 'Fully automated data processing for mapping connectivity of the marmoset prefrontal cortex.' Neuroscience, 2019. Niigata, Japan.
- [22] T. Ragan, L. R. Kadiri, K. U. Venkataraju, K. Bahlmann, J. Sutin, J. Taranda, I. Arganda-Carreras, Y. Kim, H. S. Seung, and P. Osten, 'Serial two-photon tomography for automated ex vivo mouse brain imaging,' Nature methods , vol. 9, no. 3, pp. 255-258, 2012.
- [23] P. Kask, K. Palo, C. Hinnah, and T. Pommerencke, 'Flat field correction for high-throughput imaging of fluorescent samples,' Journal of microscopy , vol. 263, no. 3, pp. 328-340, 2016.
- [24] B. Likar, J. a. Maintz, M. A. Viergever, and F. Pernus, 'Retrospective shading correction based on entropy minimization.,' Journal of Microscopy , vol. 197, no. Pt 3, pp. 285-295, 2000.
- [25] 'Allen mouse brain connectivity atlas, technical white paper: Informatics data processing,' Tech. Rep. v.9, Allen Institute for brain science, October 2017.
- [26] A. Woodward, T. Hashikawa, M. Maeda, T. Kaneko, K. Hikishima, A. Iriki, H. Okano, and Y. Yamaguchi, 'The brain/minds 3d digital marmoset brain atlas,' Scientific data , vol. 5, p. 180009, 2018.
- [27] B. B. Avants, N. J. Tustison, G. Song, P. A. Cook, A. Klein, and J. C. Gee, 'A reproducible evaluation of ants similarity metric performance in brain image registration,' Neuroimage , vol. 54, no. 3, pp. 2033-2044, 2011.
- [28] Ö. Çiçek, A. Abdulkadir, S. S. Lienkamp, T. Brox, and O. Ronneberger, '3d u-net: learning dense volumetric segmentation from sparse annotation,' in International conference on medical image computing and computer-assisted intervention , pp. 424-432, Springer, 2016.
- [29] N. Otsu, 'A threshold selection method from gray-level histograms,' IEEE transactions on systems, man, and cybernetics , vol. 9, no. 1, pp. 62-66, 1979.
- [30] L. Abdeladim, K. S. Matho, S. Clavreul, P. Mahou, J.-M. Sintes, X. Solinas, I. Arganda-Carreras, S. G. Turney, J. W. Lichtman, A. Chessel, et al. , 'Multicolor multiscale brain imaging with chromatic multiphoton serial microscopy,' Nature communications , vol. 10, no. 1, p. 1662, 2019.
- [31] K. Guesmi, L. Abdeladim, S. Tozer, P. Mahou, T. Kumamoto, K. Jurkus, P. Rigaud, K. Loulier, N. Dray, P. Georges, et al. , 'Dual-color deeptissue three-photon microscopy with a multiband infrared laser,' Light: Science &amp; Applications , vol. 7, no. 1, p. 12, 2018.

- [32] S. Ioffe and C. Szegedy, 'Batch normalization: Accelerating deep network training by reducing internal covariate shift,' arXiv preprint arXiv:1502.03167 , 2015.
- [33] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov, 'Dropout: a simple way to prevent neural networks from overfitting,' The journal of machine learning research , vol. 15, no. 1, pp. 19291958, 2014.
Henrik Skibbe Akiya Watakabe Ken Nakae Carlos Enrique Gutierrez Hiromichi Tsukada Junichi Hata Takashi Kawase Rui Gong Alexander Woodward Kenji Doya Hideyuki Okano Tetsuo Yamamori Shin Ishii
MarmoNet: a pipeline for automated projection mapping of the common marmoset brain from whole-brain serial two-photon tomography
Understanding the connectivity in the brain is an important prerequisite for
understanding how the brain processes information. In the Brain/MINDS project,
a connectivity study on marmoset brains uses two-photon microscopy fluorescence
images of axonal projections to collect the neuron connectivity from defined
brain regions at the mesoscopic scale. The processing of the images requires
the detection and segmentation of the axonal tracer signal. The objective is to
detect as much tracer signal as possible while not misclassifying other
background structures as the signal. This can be challenging because of imaging
noise, a cluttered image background, distortions or varying image contrast
cause problems.
  We are developing MarmoNet, a pipeline that processes and analyzes tracer
image data of the common marmoset brain. The pipeline incorporates
state-of-the-art machine learning techniques based on artificial convolutional
neural networks (CNN) and image registration techniques to extract and map all
relevant information in a robust manner. The pipeline processes new images in a
fully automated way.
  This report introduces the current state of the tracer signal analysis part
of the pipeline.