Algorithm for brain extraction on Magnetic Resonance
Images T1 using Morphological 3D Transformations
J. D. Mendiola-Santibañez* I. M Santillán Méndez ** C. Paredes Orta *** I. R Terol Villalobos **** *,***,**** CIDETEQ Centro de Investigación y Desarrollo Tecnológico en Electroquímica *,** División de Investigación y Posgrado de la Facultad de Ingeniería, Universidad Autónoma de Querétaro |
Keywords: brain segmentation, connected transformations, filtering, viscous transformations. |
Correspondencia: |
Palabras clave: segmentación del cerebro, transformaciones conexas, filtrado, transformaciones viscosas. |
Introduction
Mathematical Morphology (MM) is a technique where the images can be analyzed from a structural point of view. A special field of application of the MM is the segmentation of medical imaging. In particular, in this paper a composition of two morphological transformations defined for the 3D case and capable of separating brain tissue from skull is presented as proposal. The composed transformation takes into consideration the lower leveling[1] and the viscous alternated sequential filters[2]. The lower leveling propagates a marker (the marker is a part of the original image obtained from an operator) at the interior of the original volume in a controlled manner. Hence, the output image is contained into the original one. While the viscous alternated sequential filter allows disconnecting and smoothing chained components[2]. The viscous name appears because when chained components are separated, the output image visually resembles with a viscous effect.
The proposal is tested on 20 volumes of MRI T1 taken from the IBSR database[3] to separate brain. In order to compare our results, the Brain Extraction Tool (BET) -which is a popular skull-stripping algorithm- will be used. An important characteristic of BET is its quickness to separate the skull from MRI volumes of the brain. Other algorithms reported in the current literature used to carry out the same task can be found in references [4,5,6,7,8,9,10,11,12]. The composed operator provided in this paper is similar to that presented in [12]. However, there are two significant differences, i) here, the composition is obtained from the application of the lower leveling and the alternated sequential filters; while in [12], the lower leveling and the viscous opening are utilized; and ii) the execution time to segment the brain is improved in this paper.
The computed segmentations will be compared using the Jaccard[13] and Dice[14] indices. It is noteworthy to mention that, in the literature, two problems are associated during the deskulling process, 1) time used to separate brain, and 2) the quality of the segmentation. For example, in [12], the Jaccard and Dice indices have values around 0.93 and the proposal given to separate brain consumes around 15 minutes. In this paper it is introduced a morphological transformation capable of processing a MRI T1 in 3D with 60 slices in 83 seconds with Jaccard indices around 0.93.
This paper is organized as follows, in section Methodology some basic morphological transformations are defined which will be useful to introduce our proposal. In section Results and Discussion, a comparison is made between the segmented brains using our proposal and those obtained from: i) manual segmentations provided by the IBSR database[3]; and ii) BET algorithm implemented in the MRicro software[15]. In Discussion section, the advantages and disadvantages of our proposal are presented. The conclusions correspond to the last section.
METHODOLOGY
Opening by reconstruction
In MM, the basic transformations are the erosion εμB(f) and dilation δμB(f) [16], where B represents the 3D structuring element which contains its origin in the center. Figure 1 illustrates the shape of the structuring element used in this paper. denotes the transposed set of B respect its origin, = {-x : x ∈ B}, μ is a size parameter, f is the input image defined on Z3, Z represents the integer set and x is a point on the definition domain.
The morphological erosion and dilation are expressed as follows:
where ∧ and ∨ are the inf and sup operators. The morphological erosion and dilation permit to build other type of transformations; these include the morphological opening γμB(f) and closing φμB(f):
Figure 1: Shape of the 3D structuring element.
In order to reduce the notation, consider that μB = μ. The transformations having the characteristics of modifying the regional maxima or minima, without considerably affecting the remaining components, are the opening and closing by reconstruction. These operators use the geodesic transformations[17]. The geodesic dilation δf1(g) and erosion εf1(g) are expressed as δf1(g) = f ∧ δB(g) with g ≤ f, and εf1(g) = f ∨ εB(g) considering g ≥ f, respectively. When function g is equal to the morphological erosion or dilation, and the basic geodesic transformations are iterated until stability, the opening μ(f) or closing μ(f) by reconstruction are obtained,
Viscous opening and closing
These transformations allow dealing with overlapped or chained components and provide a solution to disconnect them, formally[2]:
(λ,μ)(f) = δλμ-λελ(f) with λ ≤ μ | (1) | |
(λ,μ)(f) = ελ μ-λδλ(f) with λ ≤ μ | (2) |
Figure 2: Viscous opening illustration. (a) Original image; (b) erosion size λ = 5; (c) opening by reconstruction μ - λ = 3 of the image in Fig. 2(b); (d) dilated size λ = 5 of the image in Fig. 2(c).
Viscous opening and closing utilize the parameters λ and μ. The ελ(f) and δλ(f) are used to detect those components that can be separated. To illustrate this situation consider the image in Fig. 2. Original image is presented in Fig. 2(a). In Fig. 2(b), the erosion size λ = 5 disconnects chained components; posteriorly in Fig. 2(c) the opening by reconstruction eliminates the components of size minor than μ - λ = 3 while the remaining components are maintained unchanged. Finally, in Fig. 2(d) the dilated image with λ = 5 recover the size of the original components but now they are disconnected. The components that do not support the erosion and the opening by reconstruction are merged with the background.
Viscous Alternating Sequential Filters (VASFs)
Serra[18] defined and characterized four operators where the size of the structuring element μ is indexed over a size distribution with 1 ≤ μ ≤ k.
Figure 3: Viscous Alternating Sequentail Filters illustration. (a) Original volume; (b) VASFs with λ1 = 1, λ2 = 2 and μ = 5; the filter applied is Nλ2,μ = nλ2,μnλ1,μ; (c) a threshold between 90-255 sections is computed for the image in Fig. 3(b) and lately those pixels different of 0 are recovered from the original volume.
Let us take one of these operators defined as follows in terms of connected viscous openings and closings[2]:
For fixed parameters λ1, λ2, with λ1 ≤ λ2 ≤ μ,
From last equation and for a family {λi}, with λj ≤ λk if j < k, the following VASF is defined[2]:
| (3) |
with the condition λn ≤ μ.
The VASFs will be used to segment some regions of interest. The segmentation process consist in flat the maxima and minima iteratively through the viscous opening and closing. At the end with a threshold is posible to recover the flatened zones under study. Fig. 3 illustrates the performance of Nλn,μ. Notice that, the brain has been completely separated.
Figure 4: Lower leveling illustration. (a) Original image; (b) marker by gradient, i.e., gradiμ(f) = f -εμ(f); (c) lower leveling with α = 3; (d) arithmetic difference between images in Figs. 4(a) and 4(c).
Nevertheless, the problem of Eq. 3 is that not always produces the same result and parts of the skull can be obtained in the final segmentation. Due to this situation, the VASFs will be used with the lower leveling1 to improve the final result. Lower leveling transformation is presented below.
Lower Leveling
The extended lower leveling is expressed as follows[1]:
| (4) |
where f is the reference image, g is a marker, α ∈ [0,255] is a positive scalar called slope, and μ is the size of the structuring element. The marker g used in this paper corresponds to,
| (5) |
The Eq. 5 is iterated until stability with the purpose of reconstructing the marker g at the interior of the original mask f, i.e.:
| (6) |
The image in Fig. 4, ilustrates the performance of Eq. 6. In Fig. 4(b) a marker obtained from the internal gradient is shown. The internal gradient is defined as the arithmetic difference between the original image and the eroded one, i.e., gradiμ(f) = f - εμ(f). The internal gradient uniquely detects regions on the skull. The idea of applying Eq. 6 using as marker the internal gradient is to reconstruct precisely the skull. The output image after applying Eq. 6 with α = 3 is displayed in Fig. 4(c). In this image the skull has been separated and reconstructed. Hence the arithmetic difference f - Ψμ,αn(f,g), allows us to recover the remaining components, i.e., brain, dura matter, among others, as is illustrated in Fig. 4(d). The problem of separating brain in this way is that the intensity levels of the white matter are similar to the skull bringing as consequence that the marker cannot be obtained uniquely on the bone.
On the other hand, the appropriated size of the structuring element to detect the marker utilized in this paper is found through a granulometric study. The granulometry definition is introduced below.
Granulometry
Granulometry is the distribution by sizes of particles that constitute an aggregate; it is employed in diverse areas to describe the qualities of size and shape of individual grains within a product. The concept of granulometry was introduced by G. Matheron at the end of the sixties and is presented as follows[19].
Definition (Granulometry). Let (ψλ≥0)λ∈Z+ be a family of transformations depending on an unique positive parameter λ. This family constitutes a granulometry if and only if the next three properties are verified: (i)∀ positive λ,ψλ is increasing; (ii) ∀ positive λ,ψλ is antiextensive, and (iii) ∀ positive λ and μ, ψμ = ψμψλ = ψmax(λ,μ).
Figure 5: Granulometry. The Eq. 7 is computed on a volume of MRI T1 segmented in this paper to detect the size μ corresponding to discover the brain.
The family of morphological openings and closings for the numerical case {γμ}, {φμ} with μ = {1,…,n} fulfills this last definition. The volume distribution χ is computed with Eq. 7:
| (7) |
where vol represents the volume of f, i.e., the sum of all intensity levels into the image. Figure 5 presents a granulometric study applied on a volume of MRI T1 segmented in this paper. Notice in this figure that, for sizes of μ between 1-5 the openings detect a component corresponding to the skull while openings between sizes μ in the interval 5-21 detect an important presence of the brain. From this analysis, appropiated sizes of μ to detect the marker of the brain using the morphological opening are within the interval μ ∈ [6,22].
The operators presented so far will be used to introduce a methodology to separate brain in volumes of MRI T1 for the 3D case.
Methodology to segment brain
The idea to separate brain and skull consists in combining Eq. 3 and Eq. 6 using as marker the Eq. 5. The separation of the brain consists in: i) to compute a marker of the brain using the morphological opening; ii) posteriorly, the lower leveling will reconstruct the portion of the brain into the head. However, the slope α will avoid the complete reconstruction of the marker at the interior of the volume; and iii) to clean those regions adhered to the reconstructed brain, Eq. 3 will be applied together with a threshold operation.
Formally, next filter expresses the aforementioned in steps i, ii and iii,
ξλn,μz,μx,α,μy(f)(x) = Nλn,μz[Ψμx,αn(f,γ μy(f))](x) | (8) |
Equation 8 is in term of the parameters λn, μz, μx, α and μy. Now, an explanation of how selecting the values for these parameters is provided.
- μy
- : This parameter represents the size of the structuring element of the morphological opening γμy(f) used to obtain an initial marker of the volume f. According to the granulometric study presented in Fig. 5 (see section Granulometry), adequate sizes for μy are contained in the interval μy ∈ [6,22]. In particular, in this paper, μy = 12 is utilized.
- μx
- : In this paper, μx takes the value of 1 to reconstruct the marker γμy (f) as close to the original volume f.
- α
- : This parameter is used in Eq. 6 and its purpose is to control the propagation of the marker at the interior of the original volume f. When α increases, less volume is reconstructed. In Fig. 6 this effect is illustrated. Adequated values for α according to the image in Fig 6 are in the interval 6 to 12. With α in this interval, the skull component is almost eliminated. In our case α = 10 is considered.
- λn,μz
- : These parameters correspond to the VASFs. The task of the VASFs is smoothing the remaining components of the skull after the application of the lower leveling.
Figure 6: Parameter α determination. The input volume f is processed for the 3D case; however to visualize the effect of α parameter, slices in 2D are displayed. (a) Slice taken from volume f; (b) Marker g = γμ (f) considering a structuring element size 12; (c) Eq. 6 with μ = 1 and α = 1; (d) Eq. 6 with μ = 1 and α = 1; (e) Eq. 6 with μ = 1 and α = 4; (f) Eq. 6 with μ = 1 and α = 6; (g) Eq. 6 with μ = 1 and α = 9; (h) Eq. 6 with μ = 1 and α = 12.
- An example of this situation can be seen in Fig. 6(g) where the lower leveling is applied considering α = 9. This example is important because it illustrates the approximate size of the remaining components of the skull that should be eliminated. As described before when the VASFs were introduced, once has been considered the dimension of the bone to be removed, the complete skull is supressed by applying the VASFs with λn = 2 and μz = 5 (see Fig. 3). However, due to the volumes processed with the VASFs will contain less information after applying the lower leveling, then the μz parameter can be reduced to 3 instead of 5 with the purpose of reducing the execution time. If λn = 2, the minimum value taken by μz is 3 because of next relation must be fulfilled, λ1 ≤ λ2 ≤ μz, see section of Viscous Alternating Sequential Filters.
RESULTS AND DISCUSSION
In this paper, 20 MRIs T1 (denoted as IBSR) of the head taken from the Centre for Morphometric Analysis (CMA) at the Massachusetts General Hospital (http://www.nitrc.org)[3] are processed with Eq. 8. The parameters utilized to separate the brain are shown in Table 1. From Table 1 next comments are presented:
- i.
- The parameters λn, μz, μx, α, μy associated to Eq. 8 are similar for all segmented brains.
- ii.
- The parameter Th1 represents a threshold used to eliminate cerebrospinal fluid with the purpose of disconnecting several regions between the brain and skull. According to the results presented in 20, the cerebrospinal fluid is eliminated with threshold values among 20-255 sections. Fig. 7 illustrates this situation. Original image is located in Fig. 7(a), whilst the threshold among 20-255 sections is presented in Fig. 7(b). Notice that, several components among the brain and the skull have been merged with the background, while the remaining regions are maintained unchanged. In Fig. 7(c) it is presented a mask between the original image (7(a)) and the image in Fig. 7(b). The image in Fig. 7(c) is convenient during the processing because when the leveling transformation reaches the regions where cerebrospinal fluid was located (dark regions); the propagation of the marker will be slower, obtaining in this way almost the separated brain with the leveling transformation.
- iii.
- The selected size of the opening to obtain a marker is μy = 12. This value belongs to the interval 6 to 22 according to the graph in Fig. 5.
- iv.
- The parameter λn = 2 indicates that the following filter will be applied:
with μ fixed to 3.
- v.
- The parameter Th2 allows recovering the brain region. The threshold value is according to the information provided in [20], to eliminate dark regions.
The algorithm utilized to separate the brain is given as follows. The notation .* , * and > are similar that in Matlab. .* represents the arithmetic multiply of the elements of the matrices in the same positions, * the arithmetic multiply, while > produces a logical image, true(1) for all values fulfilling the condition, and false(0) in other case.
Algorithm:
- Load the 3D image and assign to f
- threshold1 =1*( f >Th1)
- h = f .* threshold11
- μy = 12; α = 10, λn = 2; μz = 3
- Compute the marker g = γμy (h)
- Iterate Eq. 6 using g = γμy (h), the result is kept in J that represents Ψμx,αn
- Apply the filter Nλn,μ to the image J and the resultant image is kept in R that represents ξλn,μz,μx,α,μy.
- threshold2 =1*(R >Th2)
- brain=h.* threshold 2
Figure 7: Threshold parameter. (a) Original image; (b) threshold application among 20-255 sections to separate cerebrospinal fluid; (c) Mask between images in Figs. 7(a) and 7(b).
Figure 8: Illustration of the algorithm to segment brain. (a) Original volume f; (b) Mask between the original volume and the thresholded image computed from Fig. 7(a) considering sections 80-255; (c) Morphological opening size 12, γμ(y=12) (f); (d) Lower leveling considering μ = 1 and α = 8 using as marker the morphological opening, i.e., Ψμx=1,α=8n (f, γμy=12 (f)); (e) The Nλn=2,μz=3 filter is applied on the output image obtained in (d); (f) The segmented brain is obtained from a mask between the original volume and the thresholded image computed from Fig. 7(e) considering sections 90-255.
Figure 9: Brain segmentation of the volume IBSR_7. (a) Original images; (b) Manual segmentations obtained from the Centre for Morphometric Analysis (CMA) at the Massachusetts General Hospital (http://www.nitrc.org); (c) segmented brain obtained by applying the proposed algorithm; (d) segmentations obtained from BET algorithm.
Volume | Th1 | μy | α | μx | λn | μz | Th2 |
IBSR_001 | 50 | 12 | 10 | 1 | 2 | 3 | 90 |
IBSR_002 | 50 | 12 | 10 | 1 | 2 | 3 | 90 |
IBSR_004 | 50 | 12 | 10 | 1 | 2 | 3 | 90 |
IBSR_005 | 40 | 12 | 10 | 1 | 2 | 5 | 100 |
IBSR_006 | 50 | 12 | 10 | 1 | 2 | 3 | 90 |
IBSR_007 | 50 | 12 | 10 | 1 | 2 | 3 | 90 |
IBSR_008 | 50 | 12 | 10 | 1 | 2 | 3 | 90 |
IBSR_011 | 70 | 12 | 10 | 1 | 2 | 3 | 100 |
IBSR_012 | 70 | 12 | 10 | 1 | 2 | 3 | 100 |
IBSR_013 | 70 | 12 | 10 | 1 | 2 | 3 | 100 |
IBSR_015 | 40 | 12 | 10 | 1 | 2 | 3 | 80 |
IBSR_016 | 50 | 12 | 10 | 1 | 2 | 3 | 90 |
IBSR_017 | 50 | 12 | 10 | 1 | 2 | 3 | 90 |
IBSR_100 | 50 | 12 | 10 | 1 | 2 | 3 | 90 |
IBSR_110 | 70 | 12 | 10 | 1 | 2 | 3 | 100 |
IBSR_111 | 70 | 12 | 10 | 1 | 2 | 3 | 100 |
IBSR_112 | 50 | 12 | 10 | 1 | 2 | 3 | 90 |
IBSR_191 | 50 | 12 | 10 | 1 | 2 | 3 | 90 |
IBSR_202 | 50 | 12 | 10 | 1 | 2 | 3 | 90 |
IBSR_205 | 50 | 12 | 10 | 1 | 2 | 3 | 90 |
In Fig. 8 several of the steps previously mentioned are illustrated. For implementing the step f of the algorithm, the condition used as stop criterion is the volume, i.e., when the volume of the processed image does not change then the idempotence or stability has been reached.
In Fig. 9 some output images belonging to the volume IBSR_7 (see Fig. 9(a)) and processed with the algorithm aforementioned are exhibited, see Fig. 9(c). In order to compare our results, the Jaccard and Dice indices are computed on: i) manual segmentations obtained from the IBSR volumes, and ii) segmentations obtained from the BET algorithm implemented in the MRicro software with default parameters.
Some slices of the manual segmentation corresponding to the volume IBSR_7 are presented in Fig. 9(b), while a set of slices belonging to the volume segmented with the MRicro software are displayed in Fig. 9(d). The Jaccard and Dice indices associated to the aforementioned segmentations are presented in Table 2.
IBSR1 | Eq. 8 | BET
| ||
Volume | Jaccard | Dice | Jaccard | Dice |
IBSR1_001 | 0.9373 | 0.9676 | 0.7949 | 0.8857 |
IBSR1_002 | 0.9431 | 0.9707 | 0.9091 | 0.9524 |
IBSR1_004 | 0.9220 | 0.9594 | 0.8539 | 0.9212 |
IBSR1_005 | 0.8415 | 0.9140 | 0.4721 | 0.6414 |
IBSR1_006 | 0.9117 | 0.9538 | 0.5335 | 0.6958 |
IBSR1_007 | 0.9544 | 0.9767 | 0.8790 | 0.9356 |
IBSR1_008 | 0.9475 | 0.9731 | 0.7587 | 0.8628 |
IBSR1_011 | 0.9222 | 0.9595 | 0.8444 | 0.9157 |
IBSR1_012 | 0.9038 | 0.9495 | 0.8130 | 0.8968 |
IBSR1_013 | 0.9012 | 0.9480 | 0.8873 | 0.9403 |
IBSR1_015 | 0.8947 | 0.9444 | 0.3976 | 0.5690 |
IBSR1_016 | 0.9318 | 0.9647 | 0.6575 | 0.7933 |
IBSR1_017 | 0.9376 | 0.9678 | 0.6730 | 0.8045 |
IBSR1_100 | 0.9445 | 0.9714 | 0.9085 | 0.9520 |
IBSR1_110 | 0.9162 | 0.9562 | 0.9085 | 0.9520 |
IBSR1_111 | 0.9240 | 0.9605 | 0.8233 | 0.9031 |
IBSR1_112 | 0.9282 | 0.9628 | 0.8347 | 0.9099 |
IBSR1_191 | 0.9484 | 0.9735 | 0.9243 | 0.9607 |
IBSR1_202 | 0.9441 | 0.9712 | 0.9082 | 0.9519 |
IBSR1_205 | 0.9451 | 0.9718 | 0.9085 | 0.9520 |
According to the computed indices, the Eq. 8 outperforms the BET algorithm. However there are some points to mention:
- BET algorithm is faster and only requires 2 seconds to segment the volume. Nevertheless, during the segmentation of the brain, this is cropped like a circle, due to this, the Jaccard and Dice indices have minor values than those computed with Eq. 8.
- The proposed algorithm utilizes 83 seconds to segment a brain with 60 slices, and although the time is superior to the BET algorithm, the segmentations are better. In 12, the maximum time reported for other methodologies to separate brain is 35 minutes. Also in [12], two expressions to segment brain are proposed, the Eq. 9 and Eq. 12. The results obtained in this paper outperform those obtained with Eq. 12 and the reduction of the execution time is considerable. Our algorithm utilizes 83 seconds while Eq. 9 and Eq. 12 employ 174 and 354 seconds to process a volume of 60 slices respectively.
Conclusions
In this paper a morphological transformation is presented to segment brain from MRI T1 of the head. The Jaccard and Dice indices indicate that our proposal provides better segmentations than those obtained from the BET algorithm which is very popular software among the scientific community. The time used to separate the brain with our proposal is adequate (83 seconds for a volume of 60 slices), in addition the algorithm can be implemented with Matlab because several of the required transformations are already implemented in such software.
References
- Mendiola-Santibañez JD, Terol-Villalobos IR, Jimenez Sanchez AR, Gallegos-Duarte M, Rodriguez-Reséndiz J, Santillán I. Application of morphological connected openings and levelings on magnetic resonance images of the brain. Int J Imaging Syst Technol, 2011; 21(4): 336-348.
- Santillán I, Herrera-Navarro AM, Mendiola-Santibáñez JD, Terol-Villalobos IR. Morphological Connected Filtering on Viscous Lattices. Journal of Mathematical Imaging and Vision, 2010; 36(3): 254-269.
- Center for Morphometric Analysis. Massachu setts General Hospital, The Internet Brain Segmentation Repository (http://www.nitrc.org/projects/ibsr), IBSR 2014.
- Liu JX, Chen YS, Chen LF. Accurate and robust extraction of brain regions using a deformable model based on radial basis functions. Journal of Neuroscience Methods, 2009; 183(2):255-266.
- Smith SM. Fast robust automated brain extraction. Hum Brain Mapp, 2002; 17(3):143-155.
- Hahn H, Peitgen H-O. The skull stripping problem in MRI solved by a single 3D watershed transform. In Proc MICCAI, LNCS, 2000; 1935:134-143.
- Beare R, Chen J, Christopher AL, Silk T, Thompson DK, Yang JYM, et al. Brain extraction using the watershed transform from markers. Frontiers in Neuroinformatics, 2013; 7(32):1-15.
- Dogdas B, Shattuck DW, Leahy RM. Segmentation of skull and scalp in 3-D human MRI using mathematical morphology. Hum Brain Mapp, 2005; 26(4):273-285.
- Sandor S, Leahy R. Surface-based labeling of cortical anatomy using a deformable database. IEEE Trans Med Imaging, 1997; 16(1):41-54.
- Ségonne F, Dale AM, Busa E, Glessner M, Salat D, Hahn HK, Fischl B. A hybrid approach to the skull stripping problem in MRI. NeuroImage, 2004, 22(3):1060-1075.
- Iglesias JE, Liu CY, Thompson PM, Tu Z. Robust Brain Extraction Across Datasets and Comparison With Publicly Available Methods. IEEE Transactions on Medical Imaging, 2011; 30(9): 1617-1634.
- Mendiola-Santibañez JD, Gallegos-Duarte M, Arias-Estrada MO, Santillán-Méndez I, Rodríguez-Reséndiz J, Terol-Villalobos IT. Sequential application of viscous opening and lower leveling for three-dimensional brain extraction on magnetic resonance imaging t1. J. Electron. Imaging. 0001; 23(3):033010. doi: 10.1117/1.JEI. 23.3. 033010, 2014.
- Jaccard P. The distribution of the flora in the alpine zone. New Phytologist, 1912;11(2):37-50.
- Dice LR. Measures of the amount of ecologic association between species. Journal of Ecology, 1945; 26(3): 297-302.
- Rorden C, Brett M. Stereotaxic display of brain lesions. In Behavioural Neurology, 2000; 12(4): 191-200.
- Heijmans H. Morphological Image Operators. Academic Press (USA), 1994.
- Vincent L. Morphological grayscale reconstruction in image analysis: applications and efficient algorithms. IEEE Transactions on Medical Imaging, 1993; 2(2):176-201.
- Serra J. Image Analysis and Mathematical Morphology. Theoretical Advances, vol. 2. Academic Press (San Diego), 1988.
- Matheron G. Random sets and integral geometry. En: John Wiley and Sons (New York),1975.
- Mendiola-Santibáñez JD, Terol-Villalobos IR, Herrera-Ruiz G, Fernández-Bouzas A. Morphological contrast measure and contrast enhancement: One application to the segmentation of brain MRI. Signal Processing, 2007, 87(9): 2125-2150.