Skip to main content

METHODS article

Front. Hum. Neurosci., 05 May 2022
Sec. Brain Imaging and Stimulation
This article is part of the Research Topic Brain mapping: From Digital Microscopic to Gross Macroscopic Brain Imaging View all 4 articles

Neuronal Morphological Model-Driven Image Registration for Serial Electron Microscopy Sections

  • 1School of Future Technology, University of Chinese Academy of Sciences, Beijing, China
  • 2Institute of Automation, Chinese Academy of Sciences, Beijing, China
  • 3School of Artificial Intelligence, University of Chinese Academy of Sciences, Beijing, China
  • 4Center for Excellence in Brain Science and Intelligence Technology, Chinese Academy of Sciences, Shanghai, China
  • 5National Laboratory of Pattern Recognition, Institute of Automation, Chinese Academy of Sciences, Beijing, China

Registration of a series of the two-dimensional electron microscope (EM) images of the brain tissue into volumetric form is an important technique that can be used for neuronal circuit reconstruction. However, complex appearance changes of neuronal morphology in adjacent sections bring difficulty in finding correct correspondences, making serial section neural image registration challenging. To solve this problem, we consider whether there are such stable "markers" in the neural images to alleviate registration difficulty. In this paper, we employ the spherical deformation model to simulate the local neuron structure and analyze the relationship between registration accuracy and neuronal structure shapes in two adjacent sections. The relevant analysis proves that regular circular structures in the section images are instrumental in seeking robust corresponding relationships. Then, we design a new serial section image registration framework driven by this neuronal morphological model, fully utilizing the characteristics of the anatomical structure of nerve tissue and obtaining more reasonable corresponding relationships. Specifically, we leverage a deep membrane segmentation network and neural morphological physical selection model to select the stable rounded regions in neural images. Then, we combine feature extraction and global optimization of correspondence position to obtain the deformation field of multiple images. Experiments on real and synthetic serial EM section neural image datasets have demonstrated that our proposed method could achieve more reasonable and reliable registration results, outperforming the state-of-the-art approaches in qualitative and quantitative analysis.

1. Introduction

Connectomics is an important research field in exploring and understanding the nanoscale neuronal structures of the brain, which is of great significance for research on neuronal function, neurogenic diseases, and brain-inspired intelligence (Paul Alivisatos et al., 2012; Luo, 2017; Sun and Du, 2018). Serial sections electron microscopy (ssEM) imaging is one of the important techniques to obtain three-dimensional (3D) neural tissue image data at nanometer resolution (Wanner et al., 2015). For serial sections imaging (Kasthuri et al., 2007), the sections of brain tissues are cut into 30–50 nm-thick by an ultramicrotome and collected onto conductive supporters, such as electron opaque support tape (e.g., ATUM-SEM Hayworth et al., 2014) or metal support grids (e.g., ssTEM Harris et al., 2006). The sections can be imaged in parallel to save data acquisition time, which facilitates studies of large-scale neural circuit reconstruction. However, such brain imaging technique inevitably brings the loss of 3D integrity. Therefore, the image registration technique is crucial for getting these serial two-dimensional (2D) EM images into 3D volumetric form for subsequent neural circuit reconstruction and analysis (Briggman and Bock, 2012).

Images registration of neural tissue sections is not an easy task. This is because the neuronal structures change dynamically along the sections. In addition, sections cut sequentially at specific intervals are unique and contiguous. This means that the contents of the adjacent images are not identical but similar, and the similarity depends on neuronal structure variation and section thickness. Most existing image registration approaches consist of three steps: feature extraction, feature matching, and regularization using a specified transformation model (e. g. rigid, affine, diffeomorphic models, etc.). Such complex appearance changing of neural images brings enormous challenges for this pipeline's first two steps (feature extraction and matching). The key to solving serial section neural image registration is finding robust and reliable correspondences under such changeable appearance circumstances.

There are several conventional nonlinear registration methods using hand-crafted features for aligning serial section images, such as scale-invariant feature transform (SIFT) (Lowe, 2004) and block-matching. Some open-source ssEM image registration tools use these types of hand-crafted features, such as elastic alignment (Saalfeld et al., 2012) (available via TrakEM2 Albert et al., 2012). However, these types of methods could not guarantee that correct correspondence is detected in the adjacent sections with considerable appearance changing. Most of them (Wang et al., 2015) choose one of the serial images as the reference image and then sequentially perform forward or backward pairwise image registration. As a result, these methods may introduce error accumulation and propagation, which leads to artificial deformation for images that are far away from the reference image. Recently, many state-of-the-art results in deep learning have been proposed for computer vision problems. Some researchers have drawn from natural image registration methods and apply to ssEM image registration. They used the data-driven deep learning method to measure the similarity between sections and minimize the square error between the target image and the registered source image to obtain the deformation field. Yoo et al. (2017) proposed an end-to-end trained 2D convolutional autoencoder (Hinton and Salakhutdinov, 2006) to generate feature maps and then combined them with a spatial transformer network (STN) (Jaderberg et al., 2015) to obtain the vector field. Mitchell et al. (2019) leveraged a siamese convolutional network to encode images and extract features, and a coarse-to-fine recursion trained alignment module to transform the source image. Despite the above significant progress, serial section neural image registering methods do not take full advantage of neuronal morphological characteristics. Due to the lack of utilization of characteristics of image content and the variation adaptability, they can hardly capture the reliable corresponding relationship in adjacent images of neural tissues well, which may lead to registering failures. Note that there exists a strong spatial morphological relationship between the appearances of neural structures in consecutive sections. Therefore, we propose that serial image registering can benefit from the neuronal morphological modeling of neuronal structures in neural images.

In addition, instead of using image information, some methods employ the straightforward way to build the corresponding relation, in which implants are utilized as registration landmarks in the biological tissues (Pauchard et al., 2004; Yavariabdi et al., 2013). Such landmarks are steady, regular, and unvarying structures in the adjacent images. They are directly selected as correspondences to register biological tissues, reducing the difficulties of seeking robust and reliable correspondence in significant appearance-changing neural images. However, the method of using implants as correspondences for registration is limited by additional experimental operations. Therefore, we consider whether there are such relatively stable structures in neural images.

Motivated by the above observations, in this work, we model the nerve morphological structure to investigate the stable structure as the "landmarks" for serial neural image registration. In previous work (Bohao et al., 2022), the authors used the spherical deformation model to simulate neurons to investigate the effects of section thickness and neuronal structure size on the ssEM registration accuracy. Inspired by this work, we use the spherical deformation model to simulate the local neuronal structure, generate synthetic pre-aligned image data, and mathematically analyze the relationship between registration accuracy and neuronal structure shapes in two adjacent sections. These findings presented that registration accuracy is positively correlated with neuronal structure roundness, or in other words, the more regular the shape of neuronal structure, the more accurate it is to register adjacent sections. Based on the conclusion of neuronal morphological modeling analysis, we go a step further by proposing a novel framework for serial neural image registration, consisting of a neuronal morphological model-driven region selection module, correspondence extraction and global optimization module, and image deformation module. Specifically, a mask is obtained in the neuronal morphological model-driven region selection module as the adaption guidance by a deep segmentation network and neural morphological physical selection model. With this mask, the feature can be adaptively extracted by focusing on the stable region and paying less attention to changing region (unstable structure). In correspondence extraction and global optimization module, the reliable correspondences between the adjacent sections are extracted by the local feature description method SIFT-flow (Liu et al., 2011). Then, the corresponding points are adjusted globally (the whole sequence) based on an energy function to satisfy the position consistency of the extracted correspondences, reducing error propagation in computing multiple images, and preserving section continuities. Finally, each section image is deformed according to the position adjustment result of the correspondences in the image deformation module. We validate the effectiveness and efficiency of our approach on two datasets of serial EM sections.

To sum up, the main contributions of this work are as follows:

(1) We model the neuronal morphology by spherical deformation model. The relationship between registration accuracy and neuronal structure shapes in two adjacent sections is mathematically analyzed. The relevant analysis proves that more regular circular structures in the section images should be selected as "landmarks."

(2) We design a new serial section image registration framework driven by neuronal morphological model, fully utilizing the characteristics of the anatomical structure of nerve tissue and obtaining more reasonable corresponding relationship. Our proposed framework consists of a neuronal morphological model-driven region selection module, correspondence extraction and global optimization module, and image deformation module.

(3) Experimental results on different datasets show that the proposed register performs significantly better than the state-of-the-art algorithms, achieving more genuine deformation results.

The remainder of this paper is organized as follows. Section 2 briefly describes the datasets used in this paper, mathematically analyses neural image registration accuracy by neuronal morphology model, and presents the proposed registration framework. Section 3 provides experiment results on a real and synthetic dataset. Finally, we conclude the paper in Section 4 with a discussion of our findings and future research directions.

2. Materials and Methods

2.1. Data Preparation

In this subsection, we briefly describe datasets used in this paper. We use three biological datasets and one synthetic dataset to validate our serial section image registration framework and verify the registration accuracy theory. We first introduce the biological datasets, then describe how to model a tubular neuron as local sphere-like structures, and generate synthetic pre-aligned image data.

2.1.1. Serial Section Neural Images Acquisition

Drosophila ssTEM (DST) dataset: These data were released by Gerhard et al. (2013), containing 20 sections of a Drosophila brain from serial section transmission electron microscopy (ssTEM). The section thickness is 40-50 nm. The resolution is 1, 024 × 1, 024 pixels with a 4.6 nm pixel size in the xy plane.

Drosophila FIB-SEM (DFS) dataset: These data include 31 serial sections of a Drosophila brain acquired using FIB-SEM (Knott et al., 2008). The solution is 6, 684 × 6, 516 pixels with a 9.15 nm pixel size in the xy plane. Because the serial EM sections are highly anisotropic in resolution, the section thickness is set to 100 nm. The sections of this dataset are imaged in situ, while due to the imaging parameter drift, a little shift has existed. It could be corrected by simple linear alignment. The finely adjusted FIB-SEM images are often regarded as the ground truth in serial section image registration and can introduce synthetic deformation to simulate the deformation received by real slices. We can quantitatively evaluate registration methods by comparing the effect of aligning these artificially deformed images.

To generate the artificially deformed images, each section is first rotated by a random angle and shifted in a random displacement. The random angle is uniformly distributed between –90 and 90 degrees, and the random displacement is uniformly distributed between –100 and 100 pixels. Then, in order to simulate the nonlinear distortion during the section cutting, we further distort every five sections artificially using MLS (Schaefer et al., 2006) by displacing four control points. The image is evenly divided into four parts by cross line, each control point is randomly selected at random positions within the range of more than 50 pixels from the edge of each part.

Drosophila ssTEM dataset from CREMI challenge (DST-CREMI): CREMI dataset is from an adult Drosophila melanogaster brain and imaged with ssTEM. It consists of 125 serial EM images with the voxel resolution size of 4 × 4 × 40nm and pixel resolution of 1, 250 × 1, 250. CREMI datasets have undergone strict manual registration and neurite membrane annotation, which can be treated as the ground truth of the EM image to be registered. We follow the process of Yoo et al. (2017) to generate synthetic deformation data: deformed using the TPS method (Bookstein, 1989) by several random vectors on random positions. The random vectors were sampled from the normal distribution with a zero mean value, and the random positions were uniformly distributed in space.

2.1.2. Synthetic Dataset

Inspired by previous work modeled vessels as the envelope of a family of spheres in Wang et al. (2020), we model the tubular neuronal structure in biological tissue as a series of local sphere-like structures, as shown in Figure 1A. Each local structure is described by the spherical deformation model, which is analyzed in detail in Hobolth (2003). The spherical deformation model has been successfully used to fit the surfaces of neurons in the human hippocampus. Further, we intend to simulate the registration of local structure K in adjacent sections in neuronal circuit reconstruction. The local structure K is cut into two parallel planes with distance d, regarded as adjacent sections. d is regarded as section thickness. The projections of K on the cutting planes are regarded as image content. After the projecting process, i.e., the simulated imaging process, we add a translation to one image, regarded as image deformation during the slicing and imaging process. Then, we register these two images and analyze the relationship between registration accuracy and neuronal structure shape. Figure 1B illustrates the above process.

FIGURE 1
www.frontiersin.org

Figure 1. (A) A tubular-like neuron is modeled as a sequence of local sphere-like structures described by a spherical deformation model with different parameters. (B) Illustration of simulated imaging, image deformation, registration, and accuracy analysis process.

We generated a synthetic dataset based on the spherical deformation model with the above description of modeling local structure and simulated registration process. Now, we briefly introduce the spherical deformation model for the following accuracy analysis and synthetic data generation. For a local structure K, the spherical deformation model selects a point zK as origin and describes the surface of K by spherical coordinates as {z + r(θ, ϕ)ω(θ, ϕ) : 0 ≤ θ < 2π, 0 ≤ ϕ ≤ π}. The unit direction vector ω(θ, ϕ) = (cosθ sinϕ, sinθ sinϕ, cosϕ) is the vector on the unit sphere with polar longitude θ and polar latitude ϕ, and r(θ, ϕ) is the distance from z to the surface of K. For simplicity, we consider the normalized radius r(θ,ϕ)/r̄, where r̄ is the mean radius length. The normalized radius function is written as follows

r(θ,ϕ)=1+n=2m=-nnanmφnm(θ,ϕ),    (1)

where φnm(θ,ϕ) is the spherical harmonic function. The spherical harmonics are given by the following equations

φnm(θ,ϕ)={kn|m|Pn|m|(cosϕ)cosmθ,m=n,,1kn0Pn0(cosϕ),m=0knmPnm(cosϕ)sinmθ,m=1,,n,

where knm is normalizing constant and Pnm is the associated Legendre function of the first kind. In the normalized radius function, each coefficient anm of spherical harmonic function is modeled as Gaussian random variables subject to N(0,λnm). The deformation model supposes there is stationarity on K, which is obtained by assuming λnm=λn, n  2, m = -n,,n. Besides, the covariance of vectors r1, ϕ1) and r2, ϕ2) on K is as follows:

Cov(r(θ1,ϕ1),r(θ2,ϕ2))=n=2λnm=-nnφnm(θ1,ϕ1)φnm(θ2,ϕ2)          =n=2λn(kn0)2Pn(cosψ),    (2)

where cosψ = ω(θ1, ϕ1)·ω(θ2, ϕ2), and ψ represents the spatial angle between r1, ϕ1) and r2, ϕ2). The spherical deformation model also supposes λn decrease according to

1/λn=α+β(np2p),n  2,p > 2,α,β > 0.    (3)

The pre-defined parameter p determines the smoothness of K, while the other two pre-defined parameters α and β determine the global and local shape, respectively. More information about the spherical deformation model can be found in Hobolth (2003).

As described above, we generated synthetic structures containing displacements to validate registration accuracy theory. The synthetic dataset includes 121 subsets with different α and β. Since in the deformation model (3), the exponential term p non-linearly controls the convergence of λn to 0. To better demonstrate the influence of different α and β, we fix p to a proper value in all subsets. In Hobolth (2003), the author used this deformation model to fit neurons in the human brain, calculating p at [3.6, 4.4]. Here, we choose p = 4 for the following experiments. Each subset constitutes 1,000 pairs of generated shapes with r̄=150 pixels and thickness d = 40 pixels. We assume that the length of a pixel in the synthetic dataset corresponds to 1 nm in the imaging process. Thus, the above sizes in experiments correspond to a regular neuronal structure size and normal section thickness in ssTEM. The translation between images is (4.5 pixels, 4.5 pixels). Figure 2A illustrates the influence of α and β on the synthetic shapes.

FIGURE 2
www.frontiersin.org

Figure 2. (A) Simulated projections of the spherical deformation model with p = 4 and changing α and β. The background color indicated R value. (B,C) Comparison of theoretical variance and experiment results of estimated translation in X-axis on synthetic datasets. (D) The theoretical variance of estimated translation with p = 4 and changing α and β.

2.2. Registration Accuracy Analysis

In order to explore the influence of different morphological structures in nerve slices on image registration, we use the spherical deformation model mentioned above to simulate the neural morphological structure. Then, we analyze the registration accuracy between two images containing projections of a local structure.

Consider source image S and template image T in registration. These two images are different cross-sections of local structure K with a shift introduced by the slicing and imaging process, (Δt, Δs). In registration, we first extract N corresponding landmarks lik = (xik, yik), k = {S, T} in image S and image T. Landmarks are sampled at an equal angular interval on contours relative to the selected origin. Next, we estimate the translation based on the given cost function. The cost function in registration is as follows

i=1N(xiT-(xiS+u))2+(yiT-(yiS+v))2,    (4)

where u and v are the translation on X-axis and Y-axis, respectively. We align S and T by estimating the translation (Δt^,Δs^), which minimizes the cost function. Since we consider translation-only registration, landmark liT is corresponding to landmark liS, which are both sampled at the same longitude θi. We use ρik to represent the distance from the origin to the ith landmark lik. ρi is actually the projection of unnormalized radius-vector on the plane, which means ρik=r(θi,ϕk)r̄sinϕk. As mentioned above, researchers generally use ultra-thin sections of 40−50nm for volume reconstruction. As a result, the size of neuronal structures in adjacent sections are approximately identical, i.e., the mean distance ρ̄T from the origin to the neuronal structure contour in T roughly equals the mean distance ρ̄S in S. This relation can be written as EiT) ≈ EiS), and can be simplified to sinϕT ≈ sinϕS. In our analysis and experiments, we suppose sinϕT = sinϕS = sinϕ.

By computing the partial derivative of the cost function concerning u and making the derivative function be zero, we get

Δt^=1Ni=1N(xiT-xiS)=1Ni=1N(x̄T+ρiTcosθi-x̄S-ρiScosθi)         =Δt+r̄sinϕNi=1Ncosθi(riT-riS) ,

where x̄k is the X-axis coordinate of the origin in image k, and x̄T=x̄S+Δt. The second-order moment of (Δt^-Δt) is as follows

E((Δt^Δt)2)=r¯2(sinϕ)2N2E(i=1Nj=1Ncosθicosθj(riTriS)(rjTrjS))                                 =r¯2(sinϕ)2N2i=1,j=1N2cosθicosθj(CoviTjTCoviSjT)                                 =r¯2(sinϕ)2N2n=2λn(kn0)2n                                     i=1,j=1N2cosθicosθj(Pn(cosψiTjT)Pn(cosψiSjT))    (5)

where Covik1jk2 = Cov(ri, ϕk1), rj, ϕk2)). The above analysis is the same as Δs^ in the Y-axis. Equation (5) shows that the critical factor affecting the estimation accuracy is the covariance λn when the sampling interval and section latitude are fixed.

Consider the relationship between λn and neuronal structure shape. Now, we analyze why a circle-like projection is more likely to correspond to a more accurate registration result. According to the method of estimating λn based on the Fourier coefficients of a structure projection at polar latitude π/2 in Hobolth (2003), we apply a Fourier transform to the radius-vector function of a projection at polar latitude ϕ0, the radius-vector function in terms of the Fourier basis is written as

r(θ,ϕ0)=b02π+n=1(bncπcosnθ+bnsπsinnθ) ,

and the Fourier coefficients bns of sin are given by

bns=sinϕ0π02πr(θ,ϕ0)sinnθ dθ .

Under the spherical deformation model, the Fourier coefficient can be modeled as a random variable:

bns=sinϕ0π02πr(θ,ϕ0)sinnθ dθ      =sinϕ0π02πl=2m=-llalmφlm(θ,ϕ0)sinnθ dθ      =πsinϕ0l=nklnPln(cosϕ0)aln~N(0,κn),      κn=πsin2ϕ0l=n(klnPln(cosϕ0))2λl, n2.

The above analysis is also applicable to the Fourier coefficient of cos. Because of the characteristics of variance (3), λl decreases rapidly. In Hobolth (2003), the author found that κn is almost only related to the first item λn. We simplify κn as follows:

                     κnClnλn,Cln=πsin2ϕ0(klnPln(0))2,n2.    (6)

Since bns and bnc obey the same gaussian distribution, we have

bn=(bnc)2+(bns)2~κnχ2(2), n2.    (7)

We can get a good enough interval of λn with the interval estimation method by selecting a reasonable confidence level α. The estimated bound of λn is bn multiplied by a constant Cn, α determined by α. Therefore, a larger bn often corresponds to a larger λn. The accuracy of the bound of λn is determined by α. By Equations (6) and (7), a larger bn is more likely to correspond to a larger λn. We construct an index R to reflect the theoretical registration accuracy,

R=-log(n=1Nbn) .    (8)

On the one hand, based on the above analysis, we know that the larger R is, the larger λn is likely to be, and with Equation (5), the registration accuracy is expected to be worse. On the other hand, since R is an index related to the Fourier coefficients, R reflects whether a shape is close to a circle or not. We validate the mean of R on the synthetic dataset, as Figure 2A shows. The results show R increases with the shape roundness. According to Figure 2A, we select Rthres = 6 as a threshold to judge whether a local structure is a "sphere-like" object in the following process. The other three results in Figure 2 demonstrate that registration accuracy is positively correlated with shape roundness, just as above analysis.

2.3. Serial Neural Image Registration

Based on the above neuronal morphological modeling analysis, we further propose a novel framework for serial neural image registration, which utilizes the characteristics of the anatomical structure of nerve tissue to obtain a more reasonable corresponding relationship. The proposed image registration method for serial EM sections is presented in three parts: neuronal morphological model-driven region selection, correspondence extraction and global optimization, and image deformation.

2.3.1. Neuronal Morphological Model-Driven Region Selection

On the basis of Section 2.2, more regular circular regions in the section image of neural tissue should be selected to extract correspondences. Hence, in this part, we try to obtain a mask of neural images that focus on the stable structures and pay less attention to the variant structures. First, we utilize the deep membrane segmentation network "FusionNet (Quan et al., 2016)" and the traditional region segmentation algorithm watershed (Bieniek and Moga, 2000) to segment all neural structures. We then use the neural morphological physical selection model to pick the stable structures. Through the segmentation method and neural morphological physical selection model, a mask is learned to focus on the stable region and pay less attention to the changing regions (unstable structure).

FusionNet is widely used in cell membrane segmentation. Similar to the traditional U-Net framework, Fusionet includes a contraction path to extract features and a symmetrical expansion path to better locate. Furthermore, it embeds a residual module between lower sampling and upper sampling in each layer and retains the layer hopping connection of U-Net. As a result, it has more robust advantages than U-Net in neuron boundary extraction (Quan et al., 2016). Therefore, in this paper, we adopt the Fusionnet structure. The number of output feature channels of the first residual block of the network is 32, and the number of feature channels increases two times every subsequent sampling. In addition, there are four lower sampling layers. Similarly, there are four upper sampling layers. The number of feature channels is reduced by two times every upper sampling. The last layer of the network is the sigmoid function to ensure that the range of the final output probability diagram is between 0 and 1. Our loss function adopts the single-pixel MSE function, which can ensure that the width of the output cell membrane will not be thicker or thinner than the ground truth. The training dataset is from the ISBI two-dimensional electron microscopy segmentation challenge (Arganda-Carreras et al., 2015), which included a 2 × 2 × 1.5 um3 volume imaged from 30 sections and publicly available manual segmentations. During instance testing, the corresponding binary image of the cell membrane can be obtained. The watershed algorithm (Bieniek and Moga, 2000) is used for different segment cells in the process, which takes the similarity between adjacent pixels as a reference so that the pixels with similar spatial positions and similar gray values are connected to form a closed contour. Finally, the desired more rounded regions are selected by the morphological threshold Rthres. Algorithm 1 describes an implementation of the complete model-driven region selection process.

ALGORITHM 1
www.frontiersin.org

Algorithm 1: Region selection.

2.3.2. Correspondence Extraction and Global Optimization

After the specific structural regions were selected, we adopt SIFT-flow (Liu et al., 2011) to extract the reliable correspondences in these areas, assuming that the structure of the biological tissue changes independently, which means that the angle of the structure to the section plane is random. As an optical flow method, SIFT-flow searches the correspondence for every pixel. Pixelwise SIFT features between the adjacent section i and i + 1 are matched as a discrete optimization problem,

E(wi)=pimin(si(pi)-si+1(pi+wi)1,t)             +piη(|u(pi)|+|v(pi)|)             +(pi,qi)εmin(α|u(pi)-u(qi)|,d)             +min(α|v(pi)-v(qi)|,d)    (9)

where wi(pi) = (u(pi), v(pi)), and si(pi) are the displacement vector and the SIFT descriptor at location pi in section i, respectively, ε contains all the spatial neighborhoods. The third term in function 9 is used as a smoothness constraint, which constrains the adjacent pixels to have similar displacement. With the assumption that the neighboring structures of the biological tissue change independently, the extracted correspondences between the adjacent sections are more reliable compared with SIFT and block-matching methods.

Although all of the pixels in the selected regions of the adjacent sections can be used as correspondences, the vertexes of a grid placed in the section are selected due to computational burden. The positions of the extracted correspondences need to be adjusted so that the points of each correspondence at the adjacent sections have the same positions in the xy plane. To avoid error accumulation, the correspondences through all of the sections are adjusted simultaneously. In addition to the position consistency of the correspondences, the displacements of these correspondences are constrained to be smooth and small, which restricts the nonlinear deformation of the original images. An energy function with an equality constraint for the sections (i = 1, ⋯ , n) is utilized to calculate the displacement of the correspondences.

E(wi,,wi,,wn)=ipi(u(pi)2+v(pi)2)+i(pi,qi)λdist(pi,qi)((u(pi)-u(qi))2+(v(pi)-v(qi))2)s.t. pi+wi=pi+1+wi+1 for every correspondence (pi,pi+1)    (10)

where wi(pi) = (u(pi), v(pi)) is the displacement vector at location pi in section i, and dist(pi, qi) is the Euclidean distance between pi and qi. The equality constraint keeps the position consistency of the correspondences between the adjacent sections and avoids the value of the displacement vector wi in the energy function to be zero.

2.3.3. Image Deformation

With the displacement vector of the extracted correspondences, we have the positions of the points in the original section image and their positions in the aligned image for each section. Any image deformation method based on the control points could be used to transform the original section into the aligned section. Here, the moving-least squares (MLS) method (Schaefer et al., 2006) is used to warp each section image. The deformation result produced by the MLS method is globally smooth, and as a result of using rigid transformation, rigidity and scale are maintained locally so that the biological tissue can retain its local shape as rigidly as possible.

3. Experiments and Results

In this section, to illustrate the proprieties of our proposed framework, we tested the effect of sub-modules: neuronal morphological model-driven region selection module and correspondence extraction and global optimization module. Furthermore, we performed volume reconstruction experiments on the serial section neural image datasets to demonstrate the capability of the whole proposed framework.

3.1. Neuronal Morphological Model-Driven Region Selection

We assessed the efficiency of the neuronal morphological model-driven region selection. We implemented the segmentation network using Keras and used a GPU workstation equipped with a NVIDIA 2080 Ti GPU. The neural cell membrane segmentation network is trained and verified on the ISBI2015 dataset (Arganda-Carreras et al., 2015), which is from Drosophila brain imaging by FIB-SEM and has the ground truth of membrane segmentation results. We got 95% membrane segmentation accuracy in the test dataset. Then, we apply the network parameters to DFS dataset, which has the same image style and could directly utilize our pre-trained model. The images from DFS data are cut to a smaller size of 512 × 512 to suit the network input. The membrane segmentation result is shown in Figure 3B. After that, the watershed algorithm is utilized to separate the regions of neural structure. Then, through region selection guided by the neuronal morphological model, we could get our desired stable region mask as shown in Figure 3C. We could see that the selected regions are relatively regular and concentrated in the nerve fiber bundle areas.

FIGURE 3
www.frontiersin.org

Figure 3. The results of extracting membrane boundaries and the selection of the regular circular region: (A) original image; (B) the results of membrane segmentation; (C) the results of region selection.

To further verify that the selected areas are indeed the more stable areas, we performed experiments using DFS dataset to measure the registration accuracy. DFS dataset are imaged in situ and the corresponding point is the direct correspondence of the pixel position. We calculated the position deviation of the corresponding points of the mask areas, the areas outside the mask, and the overall image. The calculation results showed that the mean and variance in displacement calculated in our mask regions are smaller than other regions, as shown in Figure 4. This proves that the regular circular areas change slowly and verifies our hypothesis that the stable regions extracted are more effective for improving registration accuracy.

FIGURE 4
www.frontiersin.org

Figure 4. The mean and the variance in the displacement between the regular circular areas, other areas and the whole image in two consecutive layers.

3.2. Correspondence Extraction and Global Optimization

We quantitatively evaluated the accuracy of the correspondence detection methods on adjacent sections acquired by DFS dataset. As the sections are imaged in situ, we can measure the absolute distances of the extracted correspondences as the ground truth. Among the 30 pairs of adjacent images in DFS dataset, we randomly select five starting coordinates for cutting image pairs. A total of 150 image pairs with 512 × 512 pixels are acquired for evaluating the accuracy of the correspondences. Our method was compared with two classical methods on correspondence detection, block matching and SIFT. Block matching method uses intensity information to search the correspondences. In some previous reconstruction works, block matching method was often used for slice image alignment (Saalfeld et al., 2012). SIFT is a stable local feature descriptor, maintaining invariant to rotation, scale scaling, brightness change (Lowe, 2004).

We tested block-matching methods with block size 21 × 21 and 41 × 41, SIFT and SIFT-flow. SIFT detects correspondences, and the points in the reference section are selected to detect their correspondences in the other section by block matching and SIFT-flow. Figure 5 displays the comparison of the corresponding relationship extraction of image pairs. The mean and the variance in the distances between the extracted correspondences and the ground truth of 150 image pairs are shown in Figure 6. As expected, most correspondences lie in the selected regions, and the SIFT-flow method achieves the best result because the extracted correspondences are calculated considering the neighborhood information.

FIGURE 5
www.frontiersin.org

Figure 5. Left: the points in the reference section. Middle: the correspondences extracted in the adjacent section. Right: the correspondences extracted in the selected region map.

FIGURE 6
www.frontiersin.org

Figure 6. The mean and the variance in the distances between the extracted correspondences and the ground truth.

3.3. Evaluation of Serial Section Electron Microscope Image Dataset

The registration of the serial section recovers three-dimensional continuity and reconstructs the 3D structure of biological tissue. To test the effect of restoring the z-axis continuity of neural structure, we employed our proposed method on the serial neural section images datasets: DST and DFS. Both the subjective and objective evaluation criteria are considered as methods to evaluate the results of serial section image registration methods. We also employed our proposed method on the DST-CREMI dataset for reasonably comparing with the learning-based method ssEMNet (Yoo et al., 2017).

Subjective evaluation is usually conducted by experienced technicians according to certain rules to score. In serial neural image registration, this scoring evaluation depends on the restoration of the continuity of the nerve structure in the Z direction after registration. To visualize the continuity performance of serial registration methods, we stack the serial 2D images into volume and get the Z-axis side view of the alignment results. The side views of the reconstructed anatomical object of the DFS dataset by the ground truth and the compared methods are presented in Figure 7. We compared our proposed method with the state-of-the-art serial section registration methods: Wang CW's (Wang et al., 2015) and Elastic methods (Saalfeld et al., 2012). It is shown that the proposed method and Wang CW's method achieve solid reconstructed objects with less discontinuity. While block matching-based elastic method has more shaky alignment results than others. These are mainly caused by the differences in the ability to extract the corresponding relationships. The DFS datasets with 100 nm slice thickness has large changes in the content of adjacent slice images. Due to the limited feature expression ability of the block matching method, the elastic method easily gets the incorrect correspondences when dealing with intensity changing images, leading performance degradation. On the contrary, Wang CW's and our proposed method utilizing SIFT with stronger presentation ability can better handle this. Compared with Wang CW's and our proposed method, the results of our proposed method can better maintain the morphology of neural structure. The correspondence constraint item of our proposed method based on the flow field model can punish those abnormal matching relationships, so as to reduce the matching error and improve the accuracy of feature matching. What is more, the correspondences extraction in line with the selection of neural morphological structure can cope with the changes in image contents from biological tissue. Hence, better registration results are acquired by our proposed method.

FIGURE 7
www.frontiersin.org

Figure 7. Registration results in the yz plane on the Drosophila FIB-SEM (DFS) dataset. (A,E) Ground truth. (B,F) Our results. (C,G) Wang CW results. (D,H) Elastic results.

We also counted the evaluation of registration results by several experienced technicians. We produced a questionnaire that contains the registration results of different methods. The data of subjective evaluation are prepared as follows: the aligned sections are placed sequentially on the z-axis, and the anatomical object of the biological tissue is reconstructed by interpolating along the z-axis to match the pixel size in the xy plane. The reconstructed anatomical object should have a smooth neuron membrane and mitochondria with a solid border, which is evaluated subjectively on the randomly selected side view plane by five experienced technicians for labeling electron microscopy data. Each person has 20 images of each dataset, and the total number of the to be evaluated images of each dataset is 100. Each person can only choose the best method from the results of various registration methods for each dataset, and the method score of the dataset is increased by 1. The statistical scoring results by experienced technicians are shown in Table 1. It is obvious that our proposed method has the highest score. To some extent, it can also show that our registration effect has better visual continuity.

TABLE 1
www.frontiersin.org

Table 1. The scores of the subjective evaluation of Drosophila ssTEM (DST) and Drosophila FIB-SEM (DFS) dataset.

In order to further assess the effectiveness of our method on other datasets, we evaluated our method, Wang CW's method (Wang et al., 2015) and Elastic method (Saalfeld et al., 2012) on DST dataset. Compared with the experimental settings on the DFS dataset, the Elastic method had different parameter settings due to the different slice thicknesses (100 nm in the DFS dataset, 40−50 nm in the DST dataset). The hyperparameters of elastic method are block matching search radius and block radius. After testing, we set 100 pixels block matching search radius, 20 pixels block radius in DST dataset, while 200 pixels,40 pixels in DFS dataset. The side views of the reconstructed anatomical object of the DST dataset by these methods are presented in Figure 8. The proposed method still achieves the best visual effect considering continuity. We also notice that compared with the results in the DFS dataset, the result of the Elastic method in Figure 8 looks significantly less shaky. Wang CW's method and our proposed method do not suffer from inconsistent alignment results on the used dataset. We find that the difference in alignment resulting in a different dataset may be related to the slice thickness. The difference between adjacent images increases as the slices thickness increases. As mentioned earlier, based on the feature selection and extraction of neural morphological structure, combined with the constraints of global optimization, our proposed method has a more stable alignment performance when dealing with such intensity changing images. This also proves that our method has better adaptability to different datasets.

FIGURE 8
www.frontiersin.org

Figure 8. Registration results in yz plane on the Drosophila ssTEM (DST) dataset. (A) Wang CW results. (B) Elastic results. (C) Our results.

The structural similarity index (SSIM) (Wang et al., 2004) is usually adopted to quantitatively measure the registration accuracy as an objective evaluation criterion. We also computed this objective evaluation criterion for the DFS dataset to quantitatively evaluate the serial registration methods. Figure 9 presents the quantitative evaluation results of individual approaches, the vertical axis represents the structural similarity between each transformed source image and ground truth after serial image registration. We annotated our proposed neuronal morphological model-driven serial image registration method as "Region SF" for better plotting figures. As the section index increases, the score of the Wang CW's method declines obviously because of error accumulation and propagation, which means that artificial deformation of the images far away from the reference image is aggravated. We notice that elastic method shows higher consistency than other methods, for which we thank the effective constraint of the elastic mesh model in serial registration. However, elastic mesh model constraint also restricts the ability to process abnormal slices. When large deformation occurs (slice damage or wrong correspondence), the elastic method will propagate this abnormal condition under the transmission of spring, leading to the poor performance of the elastic method in Figure 9. Our proposed neuronal morphological model-driven serial image registration method has better SSIM value and retains the morphology of the original section images better than the comparing methods, except for the sections that suffered from the artificial generated distortion. Although some distortions exist in the sections, most of the section images still reflect the real 3D structure of the sectioned biological tissue, which proved that our proposed method could retain as much as the morphology of the acquired serial EM section images as possible.

FIGURE 9
www.frontiersin.org

Figure 9. Structural similarity index (SSIM) between the aligned sections and the ground truth.

In addition, we used the open-source dataset DST-CREMI for a reasonable comparison of our proposed method and the learning-based serial section image registration methods ssEMNet (Yoo et al., 2017), which receives extensive attention from the serial image registration community. We reproduced ssEMNet according to the parameters mentioned in the original paper and tested it on a workstation equipped with Tesla V100 Graphics Processing Units (GPUs). SSIM index is utilized to quantitatively measure the registration accuracy between the aligned results and the ground truth image and between the aligned results and the adjacent image. The quantitative results in the Table 2 shows that our proposed method achieve the comparable performance compared with ssEMNet and better than elastic and Wang CW's methods on DST-CREMI dataset. The generalization of our method is proved again, which is suitable for section image datasets from different imaging modes and different samples. We could find that our proposed method has a superior SSIM value (with adjacent). Although the feature information output from the deep network has stronger representation ability than the traditional feature descriptor, our method performs better in maintaining the continuity of biological structure by using biological prior information screening the more robust corresponding points.

TABLE 2
www.frontiersin.org

Table 2. Quantitative results comparison of structural similarity index (SSIM) on Drosophila ssTEM (DST)-CREM datasets.

4. Conclusion

In this paper, we consider whether there are such stable "markers" in the neural images to alleviate registration difficulty. We employ the spherical deformation model to simulate the local neuron structure and mathematically analyze the relationship between registration accuracy and neuronal structure shape in two adjacent sections by the second-order moment of estimated translation. Through modeling and analysis, we can prove that the registration accuracy is positively correlated with neuronal structure roundness, which is to say, the more regular the structure, the more stable it is, and it is more conducive to registration.

Based on the analysis results of neural morphological modeling, we designed a new serial image registration framework with a structure selection module, utilizing the characteristics of the anatomical structure of nerve tissue to obtain a more reasonable corresponding relationship. The proposed framework consists of a neuronal morphological model-driven region selection module, correspondence extraction and global optimization module, and image deformation module. Our proposed method leverages the deep membrane segmentation network and neural morphological physical selection model to select the changeless rounded regions in neural structure. The output region selected masks combined with feature extraction and global optimization of correspondence position to obtain the deformation field of multiple images.

Compared with the state-of-the-art approaches, the proposed method achieves better genuine deformation results and better quantitative results on serial EM section image datasets. However, we still have some limitations in that utilizing a deep segmentation network selected region and SIFT-flow method extracted correspondences to obtain the deformation field, not a totally end-to-end method. The lack of joint optimization in separate steps may lead to some reduction in registration accuracy. The proposed method is not robust enough when the images have large distortions. In the future, we will explore how to use an end-to-end method to utilize the structural information of ssEM images for better registration.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author Contributions

FZ and BC designed the experiments and wrote the article. XC and HH conceived the ideas for all the methods and provided suggestions based on the experimental results. All authors contributed to the article and approved the submitted version.

Funding

This research is supported by the Strategic Priority Research Program of Chinese Academy of Science (No. XDB32030208 to HH), International Partnership Program of Chinese Academy of Science (No. 153D31KYSB20170059 to HH), Program of Beijing Municipal Science & Technology Commission (No. Z201100008420004 to HH), CAS Key Technology Talent Program (No. 292019000126 to XC), Instrument Function Development Innovation Program of Chinese Academy of Sciences (No. E0S92308 to XC).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Albert, C., Stephan, S., Johannes, S., Ignacio, A. C., Stephan, P., Mark, L., et al. (2012). Trakem2 software for neural circuit reconstruction. PLoS ONE 7, e38011. doi: 10.1371/journal.pone.0038011

PubMed Abstract | CrossRef Full Text | Google Scholar

Arganda-Carreras, I., Turaga, S. C., Berger, D. R., Dan, C., Giusti, A., Gambardella, L. M., et al. (2015). Crowdsourcing the creation of image segmentation algorithms for connectomics. Front. Neuroanat. 9, 142. doi: 10.3389/fnana.2015.00142

PubMed Abstract | CrossRef Full Text | Google Scholar

Bieniek, A., and Moga, A. (2000). An efficient watershed algorithm based on connected components. Pattern Recognit. 33, 907–916. doi: 10.1016/S0031-3203(99)00154-5

CrossRef Full Text | Google Scholar

Bohao, C., Tong, X., Hua, H., and Xi, C. (2022). “Performance analysis in serial-section electron microscopy image registration of neuronal tissue,” in SPIE Medical Imaging 2022 (ACCEPT) (International Society for Optics and Photonics) (San Diego, CA).

Bookstein, F. L. (1989). Principal warps: thin-plate splines and the decomposition of deformations. IEEE Trans. Pattern. Anal. Mach. Intell. 11, 567–585. doi: 10.1109/34.24792

CrossRef Full Text | Google Scholar

Briggman, K. L., and Bock, D. D. (2012). Volume electron microscopy for neuronal circuit reconstruction. Curr. Opin. Neurobiol. 22, 154–161. doi: 10.1016/j.conb.2011.10.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerhard, S., Funke, J., Martel, J., Cardona, A., and Fetter, R. (2013). Segmented anisotropic sstem dataset of neural tissue. Figshare.

Harris, K. M., Perry, E., Bourne, J., Feinberg, M., Ostroff, L., and Hurlburt, J. (2006). Uniform serial sectioning for transmission electron microscopy. J. Neurosci. 26, 12101–12103. doi: 10.1523/JNEUROSCI.3994-06.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Hayworth, K. J., Morgan, J. L., Schalek, R., Berger, D. R., Hildebrand, D. G., and Lichtman, J. W. (2014). Imaging atum ultrathin section libraries with wafermapper: a multi-scale approach to em reconstruction of neural circuits. Front. Neural Circ. 8, 68. doi: 10.3389/fncir.2014.00068

PubMed Abstract | CrossRef Full Text | Google Scholar

Hinton, G. E., and Salakhutdinov, R. R. (2006). Reducing the dimensionality of data with neural networks. Science 313, 504–507. doi: 10.1126/science.1127647

PubMed Abstract | CrossRef Full Text | Google Scholar

Hobolth, A. (2003). The spherical deformation model. Biostatistics 4, 583–595. doi: 10.1093/biostatistics/4.4.583

PubMed Abstract | CrossRef Full Text | Google Scholar

Jaderberg, M., Simonyan, K., and Zisserman, A. (2015). “Spatial transformer networks,” in Advances in Neural Information Processing Systems (Montreal, QC), 2017–2025.

Kasthuri, N., Hayworth, K., Lichtman, J., Erdman, N., and Ackerley, C. (2007). New technique for ultra-thin serial brain section imaging using scanning electron microscopy. Microsc. Microanal. 13, 26–27. doi: 10.1017/S1431927607078002

CrossRef Full Text | Google Scholar

Knott, G., Marchman, H., Wall, D., and Lich, B. (2008). Serial section scanning electron microscopy of adult brain tissue using focused ion beam milling. Journal of Neuroscience 28, 2959–2964. doi,: 10.1523/JNEUROSC. I.3189-07.2008

PubMed Abstract | Google Scholar

Liu, C., Yuen, J., and Torralba, A. (2011). Sift flow: dense correspondence across scenes and its applications. IEEE Trans. Pattern Anal. Mach. Intell. 33, 978–994. doi: 10.1109/TPAMI.2010.147

PubMed Abstract | CrossRef Full Text | Google Scholar

Lowe, D. G. (2004). Distinctive image features from scale-invariant keypoints. Int. J. Comput. Vis. 60, 91–110. doi: 10.1023/B:VIS.I.0000029664.99615.94

CrossRef Full Text | Google Scholar

Luo, Q. (2017). Brainsmatics-bridging the brain science and brain-inspired artificial intelligence. Scientia Sinica Vitae. 47, 1015–1024.

Mitchell, E., Keselj, S., Popovych, S., Buniatyan, D., and Seung, H. S. (2019). Siamese encoding and alignment by multiscale learning with self-supervision. arXiv preprint arXiv:1904.02643.

Pauchard, Y., Smith, M., and Mintchev, M. (2004). Modeling susceptibility difference artifacts produced by metallic implants in magnetic resonance imaging with point-based thin-plate spline image registration. Conf. Proc. IEEE Eng. Med. Biol. Soc. 2004, 1766–1769. doi: 10.1109/IEMBS.2004.1403529

PubMed Abstract | CrossRef Full Text | Google Scholar

Paul Alivisatos, A., Miyoung, C., Church, G. M., Greenspan, R. J., Roukes, M. L., and Rafael, Y. (2012). The brain activity map project and the challenge of functional connectomics. Neuron 74, 970–974. doi: 10.1016/j.neuron.2012.06.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Quan, T. M., Hildebrand, D. G., and Jeong, W. K. (2016). Fusionnet: A deep fully residual convolutional neural network for image segmentation in connectomics. arXiv[Preprint]. arXiv: 1612.05360. Available online at: https://arxiv.org/ftp/arxiv/papers/1612/1612.05360.pdf

Saalfeld, S., Fetter, R., Cardona, A., and Tomancak, P. (2012). Elastic volume reconstruction from series of ultra-thin microscopy sections. Nat. Methods 9, 717. doi: 10.1038/nmeth.2072

PubMed Abstract | CrossRef Full Text | Google Scholar

Schaefer, S., McPhail, T., and Warren, J. (2006). Image deformation using moving least squares. ACM Trans. Graphics 25, 533–540. doi: 10.1145/1141911.1141920

CrossRef Full Text | Google Scholar

Sun, L., and Du, J. (2018). Progress in brain neural connectomics. Scientia Sinica Vitae. 48.

Wang, C. W., Gosno, E. B., and Li, Y. S. (2015). Fully automatic and robust 3d registration of serial-section microscopic images. Sci. Rep. 5, 15051. doi: 10.1038/srep15051

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Wei, X., Liu, F., Chen, J., Zhou, Y., Shen, W., et al. (2020). “Deep distance transform for tubular structure segmentation in ct scans,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (Seattle, WA: IEEE), 3833–3842.

Wang, Z., Bovik, A. C., Sheikh, H. R., Simoncelli, E. P., et al. (2004). Image quality assessment: from error visibility to structural similarity. IEEE Trans. Image Process. 13, 600–612. doi: 10.1109/TIP.2003.819861

PubMed Abstract | CrossRef Full Text | Google Scholar

Wanner, A. A., Kirschmann, M. A., and Genoud, C. (2015). Challenges of microtome-based serial block-face scanning electron microscopy in neuroscience. J. Microsc. 259, 137–142. doi: 10.1111/jmi.12244

PubMed Abstract | CrossRef Full Text | Google Scholar

Yavariabdi, A., Samir, C., Bartoli, A., Da Ines, D., and Bourdel, N. (2013). “Contour-based tvus-mr image registration for mapping small endometrial implants,” in International MICCAI Workshop on Computational and Clinical Challenges in Abdominal Imaging (Springer), 145–154.

Yoo, I., Hildebrand, D. G. C., Tobin, W. F., Lee, W. C. A., and Jeong, W. K. (2017). ssEMnet: Serial-section electron microscopy image registration using a spatial transformer network with learned features. arXiv:1707.07833 [cs.CV]. doi: 10.1007/978-3-319-67558-9_29

CrossRef Full Text | Google Scholar

Keywords: neuronal structure, spherical deformation model, image registration, registration accuracy, serial section electron microscopy

Citation: Zhou F, Chen B, Chen X and Han H (2022) Neuronal Morphological Model-Driven Image Registration for Serial Electron Microscopy Sections. Front. Hum. Neurosci. 16:846599. doi: 10.3389/fnhum.2022.846599

Received: 31 December 2021; Accepted: 04 April 2022;
Published: 05 May 2022.

Edited by:

Hasti Shabani, Shahid Beheshti University, Iran

Reviewed by:

Sergiy Popovych, Princeton University, United States
Min Liu, Hunan University, China

Copyright © 2022 Zhou, Chen, Chen and Han. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Xi Chen, xi.chen@ia.ac.cn; Hua Han, hua.han@ia.ac.cn

These authors have contributed equally to this work and share first authorship

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.