A Variational Model For Longitudinal Brain Tissue Segmentation

Abstract

Longitudinal quantification of brain changes due to development, aging or disease plays an important role in the filed of personalized-medicine applications. However, due to the temporal variability in shape and different imaging equipment and parameters, estimating anatomical changes in longitudinal studies is significantly challenging. In this paper, a longitudinal Magnetic Resonance(MR) brain image segmentation algorithm proposed by combining intensity information and anisotropic smoothness term which contain a spatial smoothness constraint and longitudinal consistent constraint into a variational framework. The minimization of the proposed energy functional is strictly and effectively derived from a fast optimization algorithm. A large number of experimental results show that the proposed method can guarantee segmentation accuracy and longitudinal consistency in both simulated and real longitudinal MR brain images for analysis of anatomical changes over time.

1. Introduction

Longitudinal quantification of anatomical changes provides a unique opportunity for human brain study and analysis in many clinical scenarios [1-3]. Estimation of anatomy changes are essential in clinical applications for disease onset prediction, disease progression evaluation, recovery and treatment efficacy quantization and common patterns of growth or atrophy determination [4-6]. In these studies, consistent tissue segmentation on Magnetic Resonance (MR) images plays an indispensible role for measuring subject changes with time [7-9]. On series images segmentation, longitudinal stability is critical since the quantification accuracy will be affected by measurement error [10]. However, due to different imaging parameters, different scanner calibration, and joint processing of multiple 3-D images, existing 3-D segmentation algorithms may not be able to provide adequate longitudinal stability. In the past two decades, lots of image techniques have been used for image segmentation, like thresholding method [11], the region growing [12], statistical models [13, 14], clustering [15, 16] and active contour models [17, 18]. However, these intensity-based methods cannot be directly applied for MR image segmentation. A commonly encountered problem is the intensity inhomogeneity or bias field of brain MR image [19]. The intensity inhomogeneity due to the device limitations and patient-induced electrodynamics interactions can overlap the strength range of different organizations. So the voxels in different tissues can not be distinguished by intensity. To solve this problem, one kind of the methods which integrate the bias field with segmentation model was widely studied [18, 20-24]. In these methods, the tasks of bias field correction and segmentation are merged and mutually benefical to yield better results. Another tissue segmentation accuracy impact is the noise. The heavy noise in the image can seriously affect the spatial consistent of tissue segmentation. Therefore, a spatial constraint was widely imposed to the energy function for spatial smooth segmentation results [21-23]. However, these existing brain MR image segmentation methods are designed for 3D brain image segmentation at single time point even though longitudinal scans are available. The segmentation results of these methods might result in inconsistent tissue variation in temporal dimension when they are applied to a series of scans of the same subject. In the previous works[10], a time-consistent segmentation algorithm for longitudinal images proposed by Xue et al., using iterated registration and segmentation. Whereas, this method which uses the local clustering means in the spatiotemporal is extremely sensitive to the initialization. The convergence of algorithm and the parameters, such as number of iterations, were also not analyzed in this paper. Furthermore, the framework iteratively performed the image registration and segmentation is very complicated and time consuming.

In this paper, we propose a novel longitudinally guided variational model for consistent tissue segmentation of longitudinal brain MR images. In our model, the data fitting term and the anisotropic smoothness term are combined together to generate a variational energy function. The data fitting term in our variational model is the longitudinal version of 3D intensity fitting model that incorporates the bias field into the segmentation model to deal with intensity inhomogeneities in the MR data. Our anisotropic smoothness term is a linear combination of the spatial and temporal consistent constrain terms. These two consistent constrain terms are generated by the total variation (TV) regularization which can simultaneously preserve the structural information whilst smooth away the noises in flat regions. With this anisotropic smoothness term, we can overcome the influence of the noise and keep the longitudinal consistent variation of the tissue segmentation results. To minimize the proposed energy function, a fast optimization algorithm is applied to quickly obtain the minimization results. Experimental results on real and virtual data show that our longitudinal segmentation results can ensure segmentation accuracy and longitudinal consistency for analysis of anatomical changes over time.

2. Background and Preprocess

2.1 3D brain image segmentation and bias correction method

Most frequently, the intensity inhomogeneity is assumed to be a multiplicative field as it is consistent with the inhomogeneous sensitivity of the reception coil. According to this assumption, the acquired image I(x) can be expressed as:

I(x) = b(x)J(x) + n(x),       (1)

where b(x) is the bias field, J(x) is the inhomogeneity-free image, n(x) is the additive noise respectively. Noise is usually ignored and logarithmic transform is performed on both sides, in order to simplify the calculation:

log(I) = log(bJ) = log(J) + log(b)       (2)

We name log(I) with \(\begin{aligned}\check{I}\end{aligned}\), log(J) with \(\begin{aligned}\check{J}\end{aligned}\) and log(b) with \(\begin{aligned}\check{B}\end{aligned}\), respectively. Therefore the equation (2) is writen as \(\begin{aligned}\check{I}=\widetilde{B}+\check{J}\\\end{aligned}\). Based on this multiplicative model, the J(x) can be approximated as a constant c_i [22-25]. And a linear combination of basis functionsg₁, g₂, … , g_m.is helpful to the bias field b(x).

According to these two assumptions, an intensity fitting model for 3-D brain segmentation and bias correction as follows: \(\begin{aligned}\mathrm{E}=\sum_{i=1}^{\mathrm{N}}\left|\check{I}(\mathrm{x})-\omega^{\mathrm{T}} \mathrm{g}(\mathrm{x})-\mathrm{c}_{\mathrm{i}}\right|^{2} \mathrm{dx}\\\end{aligned}\), whereω = (ω₁, ω₂, ⋯ , ω_m)^Tare the combination coefficients, and g = (g₁, g₂₁, … , g_m)^T are the basis functions.

This method can obtain the accurate segmentation results and smoothness bias field of 3-D MR image. However, without considering the spatial and temporal consistent of each tissue, this method is sensitive to the noise, and only can be applied for 3-D image segmentation. In this subsection, we extend this model as our data-fitting term and add to the spatial and temporal smoothness constrain terms. 

 

2.2 Preprocess

A set of 3-D brain MR images I(t) t ∈ T = {t₁, t₂, ⋯ , t_n} scanned from the same patient at different time point constitute the longitudinal brain MR image I in this article. Due to different imaging parameters, different imaging position and different scanner calibration of each 3-D brain MR image, we need to do some preprocesses for consistent 4-D brain segmentation. First, we use the rigidly registration method [26, 27] to integrate the subsequent image features into the first image features. Then, a robust skull-stripping method [28] was used to remove the non-brain tissue parts from the image. At last, the image intensities were globally normalized that the intensity interval of each image is much closer.

3. The longitudinal brain MR image segmentation model

In this section, a novel longitudinally guided variational model is proposed for consistent tissue segmentation of longitudinal brain MR images. Our global energy function is built by the data fitting term and the smoothness term. The data fitting term is generated by extending the intensity fitting term in [6] to longitudinal version. Also, the anisotropic smoothness term is the linear combination of two consistent constrain terms in spatial dimension and time dimension. With this anisotropic smoothness term, our project is noise-free and can keep the longitudinal consistent variation of the tissue segmentation results.

3.1 Data-fitting term

To simultaneous estimate the bias field and segment the longitudinal brain MR image, we adapted (2) to longitudinal version as our data-fitting term:

\(\begin{aligned}\mathrm{E}_{\text {data }}(\check{B}, \mathrm{c})=\sum_{\mathrm{i}=1}^{4} \iint_{\Omega_{\mathrm{i}}}\left(\check{\mathrm{I}}(\mathrm{x}, \mathrm{t})-\widetilde{\mathrm{B}}(\mathrm{x}, \mathrm{t})-\mathrm{c}_{\mathrm{i}}(\mathrm{t})\right)^{2} \mathrm{dxdt}\\\end{aligned}\)       (3)

where \(\begin{aligned}\widetilde{\mathrm{B}}(\mathrm{x}, \mathrm{t})=\omega(\mathrm{t})^{\mathrm{T}} \mathrm{g}(\mathrm{x})\\\end{aligned}\) and c_i(t) is the bias field value of pixel and the mean value of i-th class in time t, c = (c₁, c₂, c₃, c₄)^T is the vector of mean values in different brain tissue. By using two fuzzy membership functions u₁, u₂ to generate four regions, we rewrite the (3) as follows:

\(\begin{aligned}\left(\omega, c, \mathrm{u}_{1}, \mathrm{u}_{2}\right)=\sum_{\mathrm{i}=1}^{4} \iint\left(\check{\mathrm{I}}(\mathrm{x}, \mathrm{t})-\omega(\mathrm{t})^{\mathrm{T}} \mathrm{g}(\mathrm{x})-\mathrm{c}_{\mathrm{i}}(\mathrm{t})\right)^{2} \mathrm{M}_{\mathrm{i}}(\mathrm{x}, \mathrm{t}) \mathrm{dxdt}\\\end{aligned}\)       (4)

where M₁ = u₁u₂, M₂ = u₁(1 − u₂), M₃ = (1 − u₁)u₂, M₄ = (1 − u₁)(1 − u₂) . There are three advantages of this data-fitting term: (1) the bias field and the means in each time point are independent which can simplify the energy minimization for longitudinal image segmentation; (2) this data-fitting term can estimate the bias field and segment the brain tissue simultaneously, the estimated bias field can be used for bias correction; (3) the energy of each variables is convex .

3.2 Anisotropic smoothness term

Anatomical structures are consistent in spatial and temporal dimension throughout the developmental stages [3]. Therefore, we propose an anisotropic smoothness term to better guide the segmentation. Our anisotropic smoothness term can be written as follows:

E_smooth = E_spatial + β ⋅ E_temporal       (5)

where E_spatial and E_temporal are spatial and temporal consistent terms respectively. β is a parameter to control the balance of these two smooth terms. The characteristic of this smooth constrain term is that the whole smooth term is an anisotropic problem. It is very convenient to adjust the contribution of the temporal smooth term to strengthen/weaken the temporal smooth constrain according to the different purpose.

The spatial and temporal constraint term in this paper are generated by total variation (TV) regularization [29]. Total variation method, which is widely used in image processing, can effectively maintain edges whilst eliminating noise in flat regions, even in case of low signal-to-noise ratios

Using total variation regularization to generate our spatial and longitudinal regulation term, our method is able to ensure the smoothness segmentation results both in spatial and temporal flat regions and keep the anatomical structures very well. According to this analysis, the spatial and temporal constraint terms can be written as follows:

E_spatial(u₁, u₂) = ∑_i=1²∬|𝛻u_𝑖(x, t)|dxdt       (6)

_temporal(u₁, u₂) = ∑_i=1²|d_t(u_i(x, t))|dxdt       (7)

3.3 The unified energy function for longitudinal brain tissue segmentation

According to the analysis above, we can define the 4-D brain MR image segmentation energy function, which combines data fitting term and spatial-temporal smooth constraint term, as

\(\begin{aligned} \mathrm{E}\left(\omega, \mathrm{c}, \mathrm{u}_{1}, \mathrm{u}_{2}\right)= & \alpha \sum_{\mathrm{i}=1}^{4} \iint\left(\check{\mathrm{I}}(\mathrm{x}, \mathrm{t})-\omega(\mathrm{t})^{\mathrm{T}} \mathrm{g}(\mathrm{x})-\mathrm{c}_{\mathrm{i}}(\mathrm{t})\right)^{2} \mathrm{M}_{\mathrm{i}}(\mathrm{x}, \mathrm{t}) \mathrm{dxdt} \\ & +\beta \sum_{i=1}^{2} \iint\left|d_{t}\left(u_{i}(x, t)\right)\right| d x d t+\sum_{i=1}^{2} \iint\left|\nabla u_{i}(x, t)\right| d x d t\end{aligned}\)       (8)

where 𝛼 is the regulating parameter to control the contribution of data fitting term. The three terms in our energy function reflect the intensity variations, the spatial constrain and the temporal constraint respectively. By minimizing this energy functional, we can obtain accurate, spatial and temporal consistent segmentation results.

4. Energy minimization

Every variable of energy function E(ω, c, u₁, u₂) is convex . The energy of each variable ought to be minimized in virtue of interleaved minimization which is oprated iteratively. The minimization with respect to ω and c can be easy solved by gradient descent algorithm. According to (4), the value of b and c in each time point is mutually independent, therefore we can fix tand solve the minimization with respect to ω(t) and c(t) respectively.

For fixedu₁(t), u₂(t) and c(t) , the optimization method in [23] is used to minimize E(ω, c, u₁, u₂) with respect to ω. The expression of ω(t) is written as:

ω(t) = (A(t))⁻¹v(t)       (9)

where A(t) = ∑_i=1⁴∫g(x)g^T(x)M_i(x, t)dx and \(\mathrm{v}(\mathrm{t})=\sum_{\mathrm{i}=1}^{4} \int(\check{\mathrm{I}}(\mathrm{x})-\left.c_{i}(t)\right) g(x) M_{i}(x, t) d x\\\). Keeping u₁(t), u₂(t) and ω(t) fixed and setting the partial derivative of E(ω, c, u₁, u₂) with respect to c_i(t) equal to zero, the equation to update the c_i(t) can be given by:

\(\begin{aligned}\mathrm{c}_{\mathrm{i}}=\frac{\int(\check{\mathrm{I}}(\mathrm{x}, \mathrm{t})-\mathrm{B}(\mathrm{x}, \mathrm{t})) \mathrm{M}_{\mathrm{i}}(\mathrm{x}, \mathrm{t}) \mathrm{dx}}{\int \mathrm{M}_{\mathrm{i}}(\mathrm{x}, \mathrm{t}) \mathrm{dx}}\\\end{aligned}\)       (10)

The difficulty of the energy minimization is the minimization respect with u₁, u₂ due to the TV term. Because of the non-differentiability and non-linearity of the TV term, this problem is simple in form but computationally challenging. To solve this problem, we apply the Split Bregman algorithm to fast and robust minimizes the energy respect to u₁and u₂. Compared with other methods, the Split Bregman method is independent of regularization, continuation, or the enforcement of inequality constraints, which earns the leading edge for it. When it comes L1 regularized issues, it takes little time for the proposed method to converge [18, 25, 30, 31].

Using u₁ as an example, we fixed c, b and u₂, the energy function can be rewritten as follows:

\(\begin{aligned}\min _{u_{1} \in[0,1]}\left|\nabla u_{1}\right|+\beta\left|d_{t}\left(u_{1}\right)\right|+\alpha\left\langle r_{1}, u_{1}\right\rangle\end{aligned}\)       (11)

where r₁ = (e₁ − e₃)u₂ + (e₂ − e₄)(1 − u₂) and \(\begin{aligned}\mathrm{e}_{\mathrm{i}}(\mathrm{x}, \mathrm{t})=\left(\check{\mathrm{I}}(\mathrm{x}, \mathrm{t})-\widetilde{\mathrm{B}}(\mathrm{x}, \mathrm{t})-\mathrm{c}_{\mathrm{I}}(\mathrm{t})\right)^{2}\\\end{aligned}\). Due to containing two L1 regularized terms, we introduce two auxiliary variables \(\begin{aligned}\overrightarrow{d_{1}} \leftarrow \nabla u_{1}\\\end{aligned}\) and p₁ ← d_tu₁ and apply Bregman iteration [31] to strictly enforce the constraint variables \(\begin{aligned}\overrightarrow{\mathrm{d}_{1}}=\nabla \mathrm{u}_{1}\\\end{aligned}\) and p₁ = d_tu₁. The resulting optimization problem is given by:

\(\begin{aligned}\begin{array}{l}\quad\left(u_{1}^{n+1}, d_{1}^{n+1}, p_{1}^{n+1}\right)=\underset{u_{1} \in[0,1], \overrightarrow{d_{1}}, p_{1}}{\operatorname{argmin}}\left|\overrightarrow{d_{1}}\right|+\beta\left|p_{1}\right|+\alpha\left\langle r_{1}, u_{1}\right\rangle+\frac{\mu}{2}\left\|\overrightarrow{d_{1}}-\nabla u_{1}-{\overrightarrow{b_{1}}}^{n}\right\|^{2}+ \\ \frac{\mu \beta}{2}\left\|p_{1}-d_{t}\left(u_{1}\right)-q_{1}^{n}\right\|^{2}\end{array}\\\end{aligned}\)       (12)

\(\begin{aligned}{\overrightarrow{b_{1}}}^{n+1}={\overrightarrow{b_{1}}}^{n}+\nabla u_{1}^{n}-{\overrightarrow{d_{1}}}^{n}\\\end{aligned}\)       (13)

\(\begin{aligned}q_{1}^{n+1}=q_{1}^{n}+d_{t}\left(u_{1}^{n}\right)-p_{1}^{n}\\\end{aligned}\)       (14)

The Euler-Lagrange equation of (12) with respect to u₁ is given by:

\(\begin{aligned}\Delta u_{1}+\beta \cdot d_{t t} u_{1}=\frac{\alpha}{\mu} r_{1}+\nabla \cdot\left(\overrightarrow{d_{1}}-\overrightarrow{b_{1}}\right)+\beta\left(d_{t}\left(p_{1}-q_{1}\right)\right)\\\end{aligned}\)       (15)

A fast approximated solution of (15) is provided by Gauss-Seidel iterative scheme and projection algorithm:

\(\begin{aligned} \begin{cases} \begin{array}{l}g_{i, j, k, t}^{n}=\quad d_{1_{i-1, j, k, t}}^{x, n}-d_{1_{i, j, k, t}}^{x, n}-b_{1_{i-1, j, k, t}}^{x, n}+b_{1_{i, j, k, t}}^{x, n}+d_{1_{i, j-1, k, t}}^{y, n}-d_{1_{i, j, j, t}}^{y, n} \\ -b_{1_{i, j-1, k, t}}^{y, n}+b_{1_{i, j-1, k, t}}^{y, n}+d_{1_{i, j, k-1, t}}^{z, n}-d_{1_{i, j, k, t}}^{z, n}-b_{1_{i, j, k-1, t}}^{z, n}+b_{1_{i, j, k, t}}^{z, n} \\ h_{i, j, k, t}^{n}=\quad p_{1, i, j, k, t-1}^{t, n}-p_{1, i, j, k, t}^{t, n}-q_{1, i, j, k, t-1}^{t, n}+q_{1, i, j, k, t}^{t, n} \\ l_{i, j, k, t}^{n}=\frac{1}{6+2 \beta}\left(u_{1_{i-1, j, k, t}}^{n}+u_{1_{i+1, j, k, t}}^{n}+u_{1_{i, j-1, k, t}}^{n}+u_{1_{i, j+1, k, t}}^{n}+u_{1_{i, j, k-1, t}}^{n}+u_{1_{i, j, k+1, t}}^{n}\right. \\ \left.+2 \beta\left(u_{1_{i, j, k, t-1}^{n}}^{n}+u_{1_{i, j, k, t+1}^{n}}^{n}\right)-\frac{\alpha}{\mu} r_{1_{i, j, k, t}}+g_{i, j, k, t}^{n}+2 \beta h_{i, j, k, t}^{n}\right) \\ u_{i, j, k, t}^{n+1}=\quad \max \left\{\min \left\{l_{i, j, k, t}^{n}, 1\right\}, 0\right\} \\\end{array}\end{cases} \end{aligned}\)       (16)

Minimization with respect to \(\begin{aligned}\overrightarrow{\mathrm{d}_{1}}\\\end{aligned}\) and p₁ is performed using the following formula:

\(\begin{aligned}{\overrightarrow{d_{1}}}^{n+1}=\operatorname{shrink}\left({\overrightarrow{b_{1}}}^{n}+\nabla u_{1}^{n}, \frac{1}{\mu}\right)\\\end{aligned}\)       (17)

\(\begin{aligned}p_{1}^{n+1}=\operatorname{shrink}\left(q_{1}^{n}+d_{t}\left(u_{1}^{n}\right), \frac{1}{\mu}\right)\\\end{aligned}\).       (18)

The u₂ also can be solved by using the same optimization strategy. The whole iterative algorithm of global energy minimization is given as follows:

Algorithm

Initializing u₁⁰,u₂⁰ ;

For k = 1 : iternum

Updating c_i^k(t) and r₁^k, r₂^k ;

Update u₁^k and u₂^k ;

Updating ω^k(t) ;

if |E^k – E^k-1| < total

break;

end

End

5. Experiments and Validation

In this section, we conduct a large number of experiments on simulated and clinical MR brain images to verify the performance of longitudinal image segmentation method. Quantitative comparison results show the superiority of the proposed method.

5.1 Simulated Data without longitudinal deformation

We generate synthetic image series without longitudinal deformation. Fig. 1 shows three brain images. Each of them has same white matter, gray matter and cerebrospinal fluid but different bias field, intensity contrast and noise. The second row shows the 3-D segmentation results of parametric method [22]. The fourth and fifth images in second row show the difference of two adjacent segmentation results. Due to the noise and low intensity contrast, the segmentation results obtained by [22] are big changes and cannot maintain longitudinal smooth. The first three images in third row demonstrate the segmentation results obtained by the proposed method. The last two images in third row show the difference of the adjacent segmentation results obtained by our method. With the penalty of TV norm in both longitudinal and spatial dimensions, the results are longitudinally smooth, and the volume change of adjacent results is more stable.

Fig. 1. Examples of the segmentation results for simulated longitudinal images without longitudinal deformation. Top: simulated serial images, middle: the results obtained by 3-D parametric method, bottom: results obtained by the proposed method.

5.2 Simulated Data with longitudinal deformation

The ability of proposed method is evaluated for simulated image series with deformations segmentation in this chapter. To simulate the longitudinal deformation, five images were simulated for atrophy and intensity/contrast decrease. The images in first row show the simulated MR brain image with longitudinal deformation. The second and third rows show the segmentation results obtained by the 3-D method and our method, respectively. The 3-D method without longitudinal regulation cannot obtain accurate results. Compared to 3-D method, the proposed method considering longitudinal and spatial regulations is able to segment the images with longitudinal deformation. Fig. 2 (a) and Fig. 2(b) show the number of White Matter(WM) and Gray Matter (GM) pixels in each time point obtained by 3-D method and the proposed method, respectively. These two figures demonstrate that the proposed method with longitudinal and spatial regulations can obtain longitudinally consistent segmentation results.

Fig. 2. Comparison of temporal consistency of GM and WM obtained by 3-D parametric method and the proposed method, respectively. Left: the number of WM and GM pixels in each time point segmented by 3-D parametric method, right: the number of WM and GM pixels in each time point segmented by the proposed method.

5.3 Clinical longitudinal MR brain images

In this experiment, our method is used to segment 6 sets of clinical longitudinal MR brain images given by OASIS database (http://www.oasis-brains.org/). Images in each set of the data are obtained from a same health old adult during a period of four consecutive years. The image in each time point is 3-D T1 MR image, and its size is 256×256×128.

Fig. 3 shows a typical the segmentation result of 3-D parametric method and the proposed method respectively. The first row shows four longitudinal images of a same object. The second and the third row show the segmentation results and the difference of two neighboring point data segmentation results from 3-D parametric method and the proposed method, respectively. The fourth and last row show the corresponding results obtained by our method. From these results, we can see that the changes of longitudinal images not only located in atrophy regions but also other longitudinal smooth regions. Compared with 3-D segmentation results, the most changes of our segmentation results located in atrophy regions and the results of other longitudinal smooth regions are longitudinal smoothness. This demonstrates that the proposed method is a promising tool for longitudinal consist segmentation on longitudinal brain MR images.

Fig. 3. Comparison of temporal consistency of GM and WM in real clinical data obtained by 3-D parametric method and the proposed method, respectively.

In order to quantitatively analyze the segmentation results, the numbers of WM and GM pixels of entire brain are calculated and shown in Fig. 4. From the figures, as we can see, the curves representing the results are fairly smooth. According to our calculation results, the pixels of GM and WM gradually decrease over time. Figure 5 shows the mean and Standard Deviation (SD) of segmentation differences for 8 subjects at different time points. The smaller the SD value is, the smoother the longitudinal change is and the better the segmentation effect is. It can be seen that the changes obtained by our method are more stable than those using 3-D parametric method.

Fig. 4. GM and WM volumes in entire brain. Left: results of 3-D segmentation method. Wright: results of our longitudinal segmentation method.

Fig. 5. The average and standard deviation (SD) of segmentation difference calculated at different time points from 6 subjects. Left: results of 3-D segmentation method. Wright: results of our longitudinal segmentation method.

6. Conclusion and Discussion

In this paper, we proposed a longitudinal MR brain image segmentation algorithm by combining data term, spatial and longitudinal smoothness constraints terms. The data term combining the bias field, membership functions and the clustering centers reflects the characteristic of the image in each time point. The spatial and longitudinal smooth constraint terms are used for consistent results. The advantage of using total variation norm for longitudinal smooth term is that it is able to remove the noise in the flat regions and keep the anatomical structures very well.

There are two parameters in our energy function. We select these two parameters according to the experience. The value of α can be chose in [0.03, 0.08]. It can be shown that the spatial smoothing segmentation will be decreased with α. Analogously, the value of β can be chose in [0, 10], and the longitudinal smoothing segmentation will be increased with β. If we use β=0, our energy function is same as that of 3-D segmentation method. In this respect, our energy function is a generalized form for 3-D and 4-D brain tissue segmentation. In this paper, we use α=0.05, β=6.

In conclusion, we proposed a longitudinal brain tissue segmentation method which yields spatially adaptive and longitudinal consistent segmentation results. The advantage of our method is that it is able to yields spatially and longitudinally consistent segmentation results while adapting to the longitudinal changes of anatomical structures. Experimental results on simulated and real longitudinal MR brain images show that our longitudinal segmentation results can allow both segmentation accuracy and longitudinal consistency for analysis of anatomical changes over time.