Reconstruction of micron resolution mouse brain surface from large-scale imaging dataset using resampling-based variational model

Brain surface profile is essential for brain studies, including registration, segmentation of brain structure and drawing neuronal circuits. Recent advances in high-throughput imaging techniques enable imaging whole mouse brain at micron spatial resolution and provide a basis for more fine quantitative studies in neuroscience. However, reconstructing micron resolution brain surface from newly produced neuronal dataset still faces challenges. Most current methods apply global analysis, which are neither applicable to a large imaging dataset nor to a brain surface with an inhomogeneous signal intensity. Here, we proposed a resampling-based variational model for this purpose. In this model, the movement directions of the initial boundary elements are fixed, the final positions of the initial boundary elements that form the brain surface are determined by the local signal intensity. These features assure an effective reconstruction of the brain surface from a new brain dataset. Compared with conventional typical methods, such as level set based method and active contour method, our method significantly increases the recall and precision rates above 97% and is approximately hundreds-fold faster. We demonstrated a fast reconstruction at micron level of the whole brain surface from a large dataset of hundreds of GB in size within 6 hours.

Many methods have been proposed for digital reconstruction of the brain surface [21][22][23][24][25][26][27][28][29][30] , such as an improved geometric active contour model 21 for brain surface extraction, which solved the boundary leakage problem and is less sensitive to intensity inhomogeneity to some extent; the deformable surface model 22 , which adapts the popular human brain extraction tool (BET) through the incorporation of information on the brain geometry and MR image characteristics of the rat brain; the deformable model 24 , which evolves to fit the brain's surface by the application of a set of locally adaptive model forces; the locally linear representation-based classification method 25 for brain extraction; and the brain extraction method 30 using a hybrid level set based active contour neighborhood model. All of those methods are able to accurately reconstruct the brain surface, each for a specific purpose. However, these methods cannot be applied on a large-scale fluorescence imaging dataset of a mouse brain. On the one hand, these reconstruction methods are global methods. They are suitable for small data, such as MRI datasets (hundreds of MB size in 0.1 mm resolution); however, the computational efficiency of these methods makes applications in large brain datasets impractical. Although downsampling can reduce the data volume and make global methods feasible, the surface is not fine enough due to downsampling, which may reduce the resolution simultaneously. On the other hand, most of these methods are based either on the active contour model [21][22][23][24] or on the level set method 30 ; in principle, all of them require homogeneous signal intensity and are therefore unsuitable for the new imaging datasets. The reason is as follows: the active contour model 31,32 reconstructs the shape of an object by minimizing the variational model that drives the initial surface elements move to the real surface. Because of a lack of a strong constraint in the movement direction of the initial surface elements, the path of the initial surface elements is seriously disrupted by the intensity inhomogeneity of the signals near the surface, resulting in a low-precision reconstruction. The level set based method 33 is essential for using one or several optimal threshold values to distinguish between the foreground and background. Therefore, for the strong-intensity inhomogeneity of signals near the brain surface, one or several thresholds may meet difficulties in identifying the surface signals, and reconstructing the brain surface using the level set method may not be the optimal selection.
This obstacle propels us to establish a local method that does not lower the resolution and is not affected by the intensity inhomogeneity of the signals. Here, we proposed a resampling-based variational model (RSVM) for reconstructing the brain surface from a new imaging brain dataset. The resampling strategy makes our method local and does not lower the resolution; furthermore, fixes in the evolution direction can overcome the intensity inhomogeneity of the signals. In this model, the signals in the rays that are perpendicular to the initial brain boundary line form the resampling dataset. This model has two features: the initial boundary element must move along its corresponding ray; the final position of the initial boundary element is determined by the signals in its corresponding and adjacent rays, rather than one or several fixed thresholds. These two features make RSVM suitable for dealing with the surface signals with strong-intensity inhomogeneity. The results show that RSVM can effectively reconstruct the brain surface with recall and precision rates that are above 97% (manual reconstruction is regarded as the gold standard), which are obviously superior to the level set based method 33 and the active contour methods 31,32 . The proposed model only analyzes the resampling dataset rather than the whole dataset and is suitable for addressing a large dataset. From the whole brain dataset (hundreds of GB in size), RSVM can reconstruct the brain surface within 6 hours, approximately 300-fold faster than the level set based method.

Results
The availability of RSVM. RSVM can effectively reconstruct the surface of a mouse brain. To verify this statement, two brain datasets produced by the 2p-fMOST imaging system 11 were chosen. These two datasets contain 3478 × 4500 × 5000 and 4180 × 6000 × 4500 volume pixels, and the size of each pixel is 2 μ m × 2 μ m × 2 μ m. Fig. 1a presents the horizontal plane, coronal plane and sagittal plane of the reconstructed surface of a mouse brain (from left to right) respectively, which are derived from RSVM to analyze one dataset. To display the reconstructions in detail, we extracted two cross-sections from the 3D mouse brain dataset, labeled by the dark green line in Fig. 1a. For each cross-section, the brain boundary lines reconstructed using RSVM (red) and manual (green) are provided, as shown in Fig. 1b,c. From a comparison of the reconstructed results derived from RSVM and manually, we conclude that RSVM can provide a high-precision reconstruction of the mouse brain surface. This conclusion has been further verified quantitatively, as shown in Fig. 1d. In a quantitative evaluation, the recall and precision rates are calculated. The recall rate is defined as the ratio of the area of the overlap region between the two regions enclosed by the boundary lines derived manually and by RSVM (briefly, called the overlap area) and the area of the region enclosed by the boundary line derived manually. The precision rate is defined as the ratio of the overlap area and the area of that derived from RSVM. In Fig. 1d, 8 cross-sections were selected from the dataset in Fig. 1a for a quantitative evaluation, and the recall and precision rates are all more than 97%. Fig. 1e-h that are from the analysis of the other brain dataset reveal the same information as Fig. 1a-d and further verify the effectiveness of RSVM in the reconstruction of a mouse brain surface.
Resistance of signal interference during reconstruction. RSVM is resistant to signal interference from the mouse brain boundary in reconstructing the mouse brain surface. To illustrate this point, we use the reconstruction results from the brain dataset with 3478 × 4500 × 5000 volume pixels. Using RSVM, the whole surface of the mouse brain was reconstructed and shown in Fig. 2a. From this reconstruction, three cross-sections labeled by the red lines in Fig. 2a were selected, and their reconstructed boundary lines (red curves) are shown in Fig. 2b-d. In Fig. 2b-d, the signals labeled by the yellow squares, which are useless and may be generated by some improper operation in sampling preparation, are very close to the boundary of a mouse brain and thus severely interfere with the brain surface reconstruction. However, RSVM can almost obtain perfect brain boundary lines of the cross-sections. This feature of RSVM can be revealed by the following two facts. On the one hand, similar to other variational models that have been used for reconstructing the boundary of an object, RSVM also requires a smoothness of  the reconstructed boundary. The smoothness force reduces the influence of signal interference on the movement of the boundary elements. On the other hand, the design and direction of evolution in RSVM for each boundary element can be fixed, and the final position of the boundary element located at a predetermined region can be easily made. The predetermined region can also exclude the case in which a boundary element moves along a region of undesired signal.
The robustness of RSVM to the sampling point density. RSVM has a strong robustness to the sampling point density. Here, the sampling point density is measured by the number of discretized points on the initial boundary. Noting that RSVM produces an optimal boundary by minimizing the resampling-based variational model, and the initial boundary is necessary. We used the dataset shown in Fig. 1e for testing the robustness of RSVM to its direction of evolution in the sampling point density. In Fig. 3a,c, the reconstructed brain surfaces correspond to the initial boundaries, which are discretized into 360, 540 and 720 points. For each type of initial boundary, the corresponding reconstruction is displayed three ways, by the horizontal, coronal and sagittal plane (from left to right). Obviously, the difference between these reconstructions is indiscernible. We also quantified this difference. Compared with the reconstructions shown in Fig. 3a and Fig. 3b, we regarded the reconstruction in Fig. 3b as the gold standard and computed the recall and precision rates of the reconstructions in Fig. 3a. The result of the comparison is shown in Fig. 3d. Both of these reconstructed brain surfaces (Fig. 3a,b) contain 4500 cross-section boundary lines that were equally divided into 9 parts. In each part, every cross-sections boundary line from Fig. 3a and its gold standard was compared, and the recall and precision rates were generated. We use the average recall and precision rates for measuring the difference of this part reconstruction and the gold standard. Similarly, Fig. 3e reveals a comparison of the reconstructions in Fig. 3b,c under the assumption of the reconstruction in Fig. 3c, regarded as the gold standard, and Fig. 3f reveals a comparison of the reconstructions in Fig. 3c and Fig. 3a (gold standard). From Fig. 3d,f, we conclude that RSVM can obtain similar brain surfaces even if the sampling point density is set to a wide range.
Comparison with the spherical-coordinated variational model. Compared with our previous variation model based on the sphere coordinate system (VM_SCS) 34 , RSVM can obtain the brain surface with a higher precision. Considering the huge time consumption of VM_SCS in the reconstructing of the whole brain surface, we selected a part of the whole brain dataset, which contains 4180 × 6000 × 100 volume pixels. For these 100 cross-sections, their boundary lines were derived from RSVM and VM_SCS, as shown in Fig. 4a,b. We conclude that RSVM has obvious advantages in driving the reconstructed boundary elements to reach the concave positions of real boundary lines (indicated by the arrows in Fig. 4a,b). This conclusion is clearly illustrated by a comparison of the reconstructions of a single cross-section (Fig. 4c,d).
A comparison with other methods. Compared with the classical level set based method (LSM) 33 and the active contour method (Snake) 32 , RSVM has a stronger ability to address boundary signals with a variable intensity. We used the dataset shown in Fig. 1e to verify this result. For each cross-section, reconstructing its boundary line will require approximately 700 seconds using LSM and approximately 2300 seconds using Snake; this speed is impractical for the reconstruction of a whole brain surface. Therefore, the size of the cross-section was decreased by merging 2 × 2 × 2 pixels into a single pixel for LSM and by merging 8 × 8 × 2 pixels into a single pixel for Snake; the rule used to merge multiple pixels into a single pixel (downsampling) is assigning the average pixel value of multiple pixels to the single pixel. From the dataset and two merging datasets, the brain surfaces derived from RSVM, LSM and Snake are shown in Fig. 5a,c, respectively. To display the comparison in detail, we extracted two cross-sections from the 3D mouse brain dataset, labeled by the red and dark green lines in Fig. 5a,c. The reconstructions derived from RSVM, LSM and Snake are shown in Fig. 5a1,a2, Fig. 5b1,b2 and Fig. 5c1,c2, respectively. We can see that RSVM produces the best brain boundary lines, which are close to their corresponding real boundaries. In a quantitative evaluation, the recall and precision rates are calculated, and all the frequencies of those methods are given on each extracted cross-sections.
Additionally, RSVM analyzes the resampling dataset rather than the whole dataset, and the time efficiency of RSVM is mostly determined by the size of the resampling dataset. For the classical level set based method (LSM) and active contour method (Snake), the whole dataset needs to be analyzed. Therefore, RSVM is faster than LSM and Snake in analyzing a large dataset. To illustrate this premise, 10 cross-sections from the brain datasets that all have a size of 4180 × 6000 pixels were selected. From the selected dataset, three groups of datasets were also generated by merging 2 × 2, 4 × 4, and 8 × 8 pixels into a single pixel. For each group of the dataset, the average time consumption of the reconstructions derived from RSVM, LSM and Snake were presented in Table 1. In analyzing the large size of the dataset (4180 × 6000 pixels), RSVM is 300-fold and 960-fold faster than LSM and Snake, respectively; for a relatively large size of the dataset (522 × 750 pixels), RSVM still has a 6-fold and 38-fold speed enhancement compared with LSM and Snake, respectively.

Discussion
In this paper, we have proposed the RSVM method to reconstruct the brain surface from a new imaging brain dataset. The RSVM method is applicable to a large dataset and a brain surface with inhomogeneous signal intensity. The assumption of the RSVM method is that the mouse brain surface is smooth to some extent; this assumption is reflected in two aspects: first, there are no large changes between the adjacent layers; second, the boundary of each layer is smooth, i.e., a sharp change of the boundary is impossible when the position varies slowly. This is a normal and universal assumption for the mouse brain surface because the mouse brain surface is continuous and should not have large changes when the positions vary slowly. Therefore, the assumption is satisfied for other mouse brains. The limitation of the RSVM method is that we did not consider the situation of a shape's separation and integration; therefore, RSVM cannot reconstruct a shape that is composed of multiple connected domains.
Similar to the spherical-coordinate-based variational model (VM_SCS) 34 , RSVM drives the boundary elements to move toward the real boundary under the combined action of the gradient force and the smoothness force. The convergent positions of the boundary elements correspond to the real boundary. It is clear that without an interference signal, the gradient force is dominant; on the contrary, when the interference signal is present, the smoothness force is dominant. Based on this premise, RSVM can reconstruct the brain surface under signal interference.
Different from the spherical coordinate based variational model (VM_SCS) 34 , RSVM can reconstruct the surface of nonconvex objects. In VM_SCS, the rays that the initial boundary elements move along are all started at the original point of the spherical coordinated system. This design constrains the diversity of the evolution of the initial boundary elements, causing a large loss in reconstructing a nonconvex shape. RSVM also has a requirement that the rays orthogonal to the initial boundary should be used to fix the evolution directions of the initial boundary elements. This requirement indicates RSVM can reconstruct any type of shape providing that the difference between the initial boundary and the real boundary is small. In fact, when RSVM obtains the boundary line from the cross-section, the boundary line is regarded as the initial boundary in the next reconstruction.
The level set based method (LSM) 33 designs a variational model in the reconstruction and by minimizing the variational model, the optimal threshold can be obtained and used to distinguish between the foreground and background. Although the designed variational model requires a smoothness of the reconstructed boundary, a single threshold experiences difficulties in identifying the boundary with a variable signal-intensity. RSVM minimizes the variational model for obtaining the final positions of the initial boundary elements in the rays, and these positions form the reconstructed boundary line. The position of an initial boundary element is determined by the signals of the corresponding ray and its several adjacent positions in other rays, i.e., the threshold for identifying the boundary is variable. Therefore, RSVM is powerful for reconstructing a boundary line whose signal intensity is significantly changed.
Active contour methods (Snake) 31,32 require an initial curve, and then, under the combined effect of internal force and external force, the initial curve moves toward the boundary of the shape. Generally speaking, the internal force is a representation of the smoothness of the shape boundary and is obtained by computing the average values. The external force of classical active contour methods is expressed as the gradient or the gradient vector flow (GVF) of the image. In active contour models, boundary elements move along the gradient or the gradient vector flow. However, when the signal near the brain surface is significantly changed, the gradient or the gradient vector flow is affected by the interference signal; therefore, the reconstruction is dependent on the initial curve, and the circumcircle of the boundary is used as the initial curve of Snake. While in the RSVM method, the boundary element moves along its ray direction, this design allows the RSVM method to overcome the situation of a significantly variable signal near the brain surface.
To address the high-throughput image stacks, computational efficiency is an essential factor. RSVM only analyzes the resampling dataset rather than the whole dataset; the resampling strategy makes our method local and does not lower the resolution and is suitable for addressing a large dataset. From the whole brain dataset (hundreds of GB in size), RSVM can reconstruct the mouse brain surface within 6 hours in an Inter(R)Xeon(R)CPU 3.10 GHz computing platform.

Methods
The resampling-based variational model for the reconstruction of a brain surface. The position of a voxel in brain image stacks is labeled by the coordinate (x, y, z) in an orthogonal coordinate system. The x-y plane of the imaging dataset is briefly called the coronal plane, and the boundaries of a series of coronal planes associated with the z direction form the brain surface. The boundary of a coronal plane can be discretized and expressed by a group of nose-to-tail points. To acquire the boundary of a coronal plane, the initial boundary needs to be pre-provided. In the boundary reconstruction of the first coronal plane (z = 1), the initial boundary is given manually. The reconstructed boundary of the coronal plane (z = 1), regarded as the initial boundary, is used for reconstructing the boundary of the coronal plane (z = 2). By repeating this procedure, a series of reconstructed boundaries that form the brain surface are produced, as shown in Fig. 6a. Considering the continuous property of the brain surface, the boundary differences between the two adjacent coronal planes are subtle, and this setting of the initial boundary is reasonable.
After ascertaining the initial boundary, we designed a resampling-based variational model that can drive the initial boundary to evolve to the real boundary. For the given initial boundary element denoted by s, we assume that its movement directions should be perpendicular to the initial boundary, i.e. the exterior or interior normal direction. According to this assumption, the position of s can be denoted by r 0 + r(s)*n. Here, r 0 is the initial position of the element s, n is set to the exterior normal direction, and r(s) is an unknown scale. When r(s) changes in a given interval, the set of the corresponding positions is called the ray determined by s. The notes of s and n are also further illustrated in Fig. 6b. The real boundary can be modeled as the function r(s) when the initial boundary is set. To resolve r(s), we assumed that the boundary of the coronal plane is smooth to some extent, i.e., a sharp change of r is impossible when s varies slowly; and that the intensity of the differential signal on the boundary is stronger than that in the inside and outside of the coronal plane. Based on these two assumptions, we design and minimize the following variational model for resolving r(s), written as in ext Here, E in (r(s)) and E ext (r(s)) are functions with respect to r(s). E in is the internal force that represents the smoothness of the boundary, i.e., the distance between the adjacent boundary elements cannot be too large, which allows our method to overcome the situation of a significantly changed signal near the brain surface; E ext is the external force, which represents a valid gradient of resampling rays and makes the initial boundary elements move toward the real boundary. Under the combined effect of the internal force (E in ) and external force (E ext ), we can reconstruct the boundary of a coronal plane. The small integration value of E in (r(s)) indicates that the function r(s) is slowly changed when s varies, corresponding to   Table 1. A contrast of the computation time among RSVM, LSM and Snake. The 1st column is the pixel binning mode. The 2nd column is the size (pixels) of the image after pixel binning. The 3rd column is the computation time to calculate the image with the size in the 2nd column using RSVM. The 4th column is the computation time to calculate the image with the size in the 2nd column using LSM. The 5th column is the computation time to calculate the image with the size in the 2nd column using Snake. The units of time in the 3rd to 5th columns are second (s).
Scientific RepoRts | 5:12782 | DOi: 10.1038/srep12782 the smoothness of the boundary. Maximizing the integration value of E ext (r(s)) will drive the element s to move to the boundary position where the intensity of the differential signal is the strongest. Therefore, in the parameters setting, a large value of α and a small value β are much more likely to stress the smoothness of r(s) than to stress the detailed reconstructions. Conversely, the more detailed reconstructions can be obtained. From the above analysis, E in (r(s)) and E ext (r(s)) are designed as Here, g(p) denotes the differential signal of the position r 0 + p*n, in which r 0 is the initial position of s, and p is perpendicular to the initial boundary and determined by r 0 . σ is the parameter used to measure the difference between p(s) and r(s). Obviously, finding the peak position of the function g(p) is equivalent to maximizing E ext (r(s)), and using the optimal r* from maximizing E ext (r(s)), the final position of s can be determined by r 0 + r*n and regarded as the reconstructed boundary position.
The numerical solution of the designed variational model. To resolve r(s) from the variational model (1), we introduce the procedure used in the active contour method 31,32 or the solving spherical-coordinated variational model 34 . It is well known that if r(s) is the optimal solution that minimizes the variational model (1), r(s) must satisfy the Euler equation In fact, when r(s, t) stabilizes, i.e., the right hand side of the equation (5) is near zero, Equation (5) approaches equation (4). Equation (5) can be given as  (3) is substituted into equation (5). To obtain the stable r(s, t), we discretize the equation (6), and the corresponding differential equation is given as signal intensity in the enclosed region is generally more than that outside of the enclosed region, indicating that L1 is far less than L2. This setting is based on following facts: minimalizing the variational model (1) requires that the signal value of g(p m, s ) is more than zero and as large as possible when p m, s is on the real boundary, and the signal of the position inside the boundary is more than that outside the boundary, indicating that the value of the differential signal on the boundary is less than zero. Therefore, we take the converse of the differential signal. In addition, from the above analysis, we conclude that the nonnegative values of the differential signal are useless for the minimalizing variational model (1), and thus we set these nonnegative values to zeroes. Note that the signals of g(p m, s ) are shown in Fig. 6d.
Step 3. Brain surface reconstruction. For all the initial boundary elements, the initial corresponding scale r(s i ) are all set to zeros, and by repeatedly calculating the equation (8), the convergent values of these scales can be ascertained. Using these convergent values, we can calculate the final positions of these boundary elements using the expression in section Step 1. These final positions form the reconstructed boundary of the given coronal plane associated with the z-direction, and these reconstructed boundaries, arrayed along with the z-direction, form a reconstructed brain surface, as shown in Fig. 6a,e. In equation (8), there are some parameters that need to be predetermined. We set the parameter β to With σ being 0.5 and β 0 being 0.2, and set Δ t and α to 1 and 0.8, respectively. RSVM reconstructed the mouse brain surface of two datasets in Fig. 6a,e by giving two initial boundary lines.

Sample preparation.
All the experiments were performed in accordance with the guidelines of the Experimental Animal Ethics Committee at Huazhong University of Science and Technology. All the experimental protocols were approved by the Institutional Animal Ethics Committee of Huazhong University of Science and Technology. Thy1-EGFP M line transgenic mice were used in this study. The sample-preparation procedures were as follows: 1) deeply anesthetized the mice; 2) intracardial perfusion; 3) extracted the entire brain and post-fixed it; 4) rinsed the brain; 5) dehydrated the brain; 6) infiltrated the entire brain; and 7) embedded the brain in gelatin capsules filled with pre-polymerized GMA and polymerized it. A detailed description of the parameters, such as the solutions and its concentrations and the timing and frequency in every sample preparation procedure is provided elsewhere 8,10 .
Imaging system. The imaging system, called two-photon fluorescence micro-optical sectioning tomography (2p-fMOST) 11 , can be used to section and image the entire embedded mouse brains and to obtain the large-scale imaging datasets. All the image tiles of every coronal plane were recorded in a folder automatically, and then through image pre-processing procedures 35 , the raw image tiles were montaged into the coronal plane sequence images. This step removed the inhomogeneous illumination pattern and redundant portions at the same time. The original size of the volume pixel is 0.5 μ m × 0.5 μ m × 2 μ m. The total amount of uncompressed volume of each brain was approximately 1 TB and contained approximately 5000 coronal sections. To balance the resolution of the different directions, the size of the cross-section was decreased by merging 4 × 4 pixels into a single pixel, and the original brain dataset was converted to the testing dataset with the size of the volume pixel of 2 μ m × 2 μ m × 2 μ m.