SciELO - Scientific Electronic Library Online

vol.34 issue3New insight on the relationship between lethal electrical fields versus cardiomyocyte orientationAutomatic respiratory phase detection for functional electrical stimulation synchronization author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand




Related links


Research on Biomedical Engineering

Print version ISSN 2446-4732On-line version ISSN 2446-4740

Res. Biomed. Eng. vol.34 no.3 Rio de Janeiro July./Sept. 2018  Epub Sep 13, 2018 

Original article

Robust pulmonary segmentation for chest radiography, combining enhancement, adaptive morphology and innovative active contours

Daniel Aparecido Vital1  * 

Barbara Teixeira Sais1 

Matheus Cardoso Moraes1 

1 Laboratory of Image and Signal Processing, Institute of Science and Technology, Federal University of São Paulo, São José dos Campos, SP, Brazil.



Statistical data reveal that approximately 140 million radiological exams are performed annually in Brazil. These exams are designed to detect and to analyze fractures, caused by different types of trauma; as well as, to diagnose pathologies such as pulmonary diseases. For better visualization of those lesions or abnormalities, methods of image segmentation can be implemented. Such methods lead to the separation of the region of interest, which allows extracting the characteristics and anomalies of the desired tissue. However, the methods developed by researchers in this area still have restrictions. Consequently, we present an automatic pulmonary segmentation approach that overcomes these constraints.


This method is composed of a combination of Discrete Wavelet Packet Frame (DWPF), morphological operations and Gradient Vector Flow (GVF). The methodology is divided into four steps: Pre-processing - the original image is enhanced by discrete wavelet; Processing - where occurs a combination of the Otsu threshold with a series of morphological operations in order to identify the pulmonary object; Post-processing - an innovative form of using GVF improves the binary information of pulmonary tissue, and; Evaluation – the segmented images were evaluated for accuracy of detection the pulmonary region and border.


The evaluation was carried out by segmenting 247 digital X-ray challenging images of the thorax human. The results show high for values of Overlap (97,63% ± 3.34%), and Average Contour Distance (0.69mm ± 0.95mm).


The results allow verifying that the proposed technique is robust and more accurate than other methods of lung segmentation, besides being a fully automatic method of lung segmentation.

Keywords Lung segmentation; Chest radiographs; Discrete wavelet packet frame; Gradient vector flow; Binary morphology


According to the Brazilian College of Radiology and Diagnostic Imaging (Colégio Brasileiro de Radiologia e Diagnóstico por Imagem - CBR), about 140 million x- ray examinations are performed annually in Brazil (Departamento…, 2015). Approximately 30% of these procedures aim at the pulmonary region, due to trauma or some pathology. Among the different types of pulmonary pathologies, Tuberculosis, Pneumonia and Lung Cancer are the leading causes of death in the world (Departamento…, 2015; Priya et al., 2017 ). Tuberculosis is one of the leader causes of adult mortality, since it is responsible for the death of 2 million people annually ( Costa et al., 2015 ). In addition to that, it is manifested in eight million people around the world yearly where 50% of untreated individuals die because of this disease ( Kumar et al., 2010 ). Pneumonia is a common pathology that affects approximately 450 million people yearly all over the world ( Ruuskanen et al., 2011 ). It is one of the leading causes of death among all age groups, resulting in 4 million deaths, 7% of the world's annual death toll ( Lorgelly et al., 2010 ; Ruuskanen et al., 2011 ). Lung cancer is the most common cancer in terms of incidence and mortality worldwide; additionally, it is the one with the highest mortality between men. Approximately 1.6 million people are estimated to be affected by lung cancer annually, causing 1.38 million of deaths ( Ferlay et al., 2010 ). The abovementioned abnormalities lead to changes in morphological patterns, detectable in any exam of a radiological imaging modality. Those patterns can be variation of the organs, borderlines, size or structural shape and textures ( Zhu et al., 2008 ). However, in many cases, these patterns are hardly detected or seen by radiologists or doctors. Therefore, there is a great demand for fast and accurate tools and methods that support diagnoses; hence, increasing the chances of patient recovery ( Zhu et al., 2008 ).

Radiological examination is the oldest and most used medical imaging technique to assist diagnosis and medical treatment ( Chen and Suzuki, 2014 ). X-ray based image modalities discriminate tissues by representing the intensity of the attenuation coefficient of each tissue ( Chaudhary and Patil, 2013 ). Specifically, chest radiography is the most recommended technique for the diagnosis of pulmonary abnormalities, and it is widely used ( Chen and Suzuki, 2014 ). However, the image per se, without any processing, may lack information for accurate analysis, which may lead to subjective diagnosis and thus, poor pathological discrimination. This results in a seriously compromised diagnosis. And the use of techniques that allow the improvement of the evaluation from a single image of x-ray allows reduce radiation index in human that need to be evaluated ( Michel-González et al., 2011 ). Consequently, image segmentation of the pulmonary region is of fundamental importance to extract characteristics, bringing detailed objective information for an accurate quantification and pathology classification.

Image segmentation is the process of delimiting the image in related regions. This separation of the region of interest allows extracting the characteristics and anomalies of the desired tissue (Udupa et al., 2006). Consequently, it is possible to quantify and to analyze features, which are difficult to detect or to distinguish with the naked eye ( Sokashe, 2013 ). Therefore, the improvement and creation of new methods of pulmonary segmentation for chest radiographies is of extreme importance for better diagnosis, therapies and interventions ( Chen and Suzuki, 2014 ). Specifically, it provides quantitative and objective anatomical information, such as shape, area morphology, etc. Precisely and objectively infer about morphological anomalies are very important to support anatomically related diagnostic, therapeutic procedures and treatments; such as Chagas disease, and cancer removal or treatment. In an indirect way, a lung segmentation is an important step to identify inner pulmonary pathologies. Once the lung region is isolated, a reduced amount of information is left to be searched; hence, it can improve the seek for diseases in visual inspection; as well as being used as an initial step to build another computational method devoted to segment or extract a specific disease.

The state of the art presents different segmentation methods based on threshold intensity, deformable models, edge detection or pixel classification. In Armato et al. (1998) , pulmonary segmentation was based on combining local and global intensity level threshold and contour smoothing. In the work by McNitt-Gray et al. (1995) common pixel characteristics for different regions such as heart, lung and clavicle are the bases for the segmented region. These methods may be useful, but they present loss of detail in the pulmonary borders that may be extremely necessary depending on the application. Other methods have been proposed to solve such problem and to increase the effectiveness of pulmonary segmentation. Van Ginneken et al. (2006) make use of supervised methods for chest anatomical structures segmentation. Deformable flexible models are applied by Stegmann et al. (2003) . Shi et al. (2008) presented a hierarchical and global Principal component analysis (PCA) technique method for lung segmentation. In Rebouças et al. (2013) , active contour techniques, such as Gradient Vector Flow (GVF), performed the 2D lung segmentation in computed tomography (CT) images of the thorax. They observed the ability of this method to identify the edges of the lung tissue in detail, since when performing the gradient of the images, the border locations are the regions that present less energy. However, GVF carry out limitations by being semi-automatic, and without a proper contour initialization, it may end up having a high computational cost, given the number of interactions required to make the contour go towards the border. Wan Ahmad et al. (2015) performed lung segmentation in chest radiographs using oriented Gaussian derivatives filter. Dai et al. (2017) recognized standards for segmenting organs in thorax-rays through network training and correction. Also, in the work by Xiong et al. (2017) , the lung segmentation occurred by a weighted sparse shape composition with robust initialization.

These methods present elegant solutions in which accuracy, efficiency and practicability have gradually evolved with each work. However, some limitations such as the lack of details at organs borders and dependence on more than one exam or image for improving accuracy could be found. The lack of border information may be critical, since the frontier of a lesion and a vital part of an organ may be a concern for a surgery; additionally, as shape is totally related to vital information, a small lack of border information may dangerously lead to a totally different diagnosis. In the case of dependence on many acquisitions to provide enough information for the diagnosis, it is totally against the image processing techniques duties. As the radiological exams are quite unhealthy to the patient, due to the use of radiation, image-processing methods are meant to decrease the need of high doses. In other words, image-processing methods may be able to support hardware limitation, providing high quality of relevant information from heterogeneous acquisitions. Consequently, while thinking and developing methods, different manufactures and/or hospitals are welcome to provide images to evaluate not only accuracy, but also precision, practicability and method robustness, regarding manufacturing. Therefore, the search for more complete and ideal methods is still ongoing, towards gathering practicability, efficiency and robustness.

Seeking to satisfy such capabilities, we performed a combination of different image processing techniques to overcome current limitations, obtaining superior effectiveness and robustness. The enhancement provided by the discrete wavelet transform, during pre-processing, increased the discrimination between lung regions and background ( Vidaurrazaga et al., 2000 ). Next, we performed the Otsu technique ( Otsu, 1979 ), to optimally binarize the highly bimodal histogram image produced by the previous discrete wavelet enhancement ( Katouzian et al., 2008 ); with information binarized, a series of morphological operations took place to polish the desired information ( Moraes and Furuie, 2010 ), the pulmonary object. Finally, so as to produce a fine adjustment on the final segmentation, we carried out an adapted and automatic version of the Gradient Vector Flow (GVF) method (Xu and Prince, 1998). In this adapted version, we take advantage of the previously obtained pulmonary objects to initialize the GVF contour; hence, reducing the computational cost, since the objects are in lung limit, ( Rebouças et al., 2013 ). Such combination of methods aims to increase practicability, efficiency and reduction of computational cost. Consequently, contributing with an alternative lung segmentation methodology.

Finally, the general goal of this paper is to present an alternative fully automatic method of lung segmentation in chest radiography consisting of the combination of techniques presented to overcome the limits of the effectiveness of current methods.


The objective of this work is to present and to validate an alternative method that combines medical image processing techniques for highly accurate lung segmentation of chest radiography.

The material used consists of chest X-ray images from the online database of the Japanese Society of Radiological Technology (JSRT) ( Shiraishi et al., 2000 ). This is composed of 247 digital X-ray images of the thorax, being 154 images with nodules and 93 non-nodular images. These images are available for research and education purposes ( Shiraishi et al., 2000 ). Additionally, a Desktop computer with an Intel® CORE i7, and microprocessor of 3.5 GHz, 16 GB of RAM, Windows 8.1 Pro 64 bits, with MATLAB® (R2017a) (MathWorks, Inc., Natick, MA, USA) were used.

The methodology is divided into 4 steps: Pre-processing, Processing , Post-processing and Evaluation ( Figure 1 ). In the pre-processing block ( Figure 1 in yellow), the original image (IO) ( Figure 1 a) went through an image enhancement technique: the Discrete Wavelet Packet Frame (DWPF) ( Figure 1 b) ( Coifman and Wickerhauser, 1992 ). DWPF increased the discrimination between lung region and background by over filtering image data. During the processing ( Figure 1 in blue), the image information highlighted by the DWPF method (IRW ) was binarized using the Otsu method, resulting in (IRWB ) ( Figure 1 c) ( Otsu, 1979 ). Accordingly, a series of morphological operations were performed to correct the missing information ( Figure 1 d); as a result, two binary images, one of the left (IBL) ( Figure 1 e) and another of the right (IBR) lungs were obtained ( Figure 1 f). Next, in the post-processing ( Figure 1 in red), the edge information of the binary images of the left and right lungs (I BL and IBR) were used to initialize the GVF method ( Xu and Prince, 1998 ) ( Figure 1 h). The gradient image of IRC, obtained by the enhancement technique of Contrast-limited adaptive histogram equalization (CLAHE) ( Ikhsan et al., 2014 ) ( Figure 1 g), used to move the contour in GVF. After this process, two segmented images referring to the left lung (IACL) and right lung (IACR ) were obtained ( Figure 1 i and j). Finally, in the evaluation block, the segmented images were evaluated concerning their segmentation accuracy of the pulmonary region and border detection ( Figure 1 k) in gray). This evaluation was performed by comparing the segmented images with their respective gold standards (IG) ( Figure 1 l) in two parts: first, the parameters of segmentation accuracy proposed by Udupa et al. (2006) were calculated; secondly, the Average Contour Distance ( Candemir et al., 2014 ) between the edges of segmented and gold standards images were computed.

Figure 1 Combination of binary operations for obtaining the binary objects of the left and right lungs. (a) original image (IO); (b) Image enhancement technique of Discrete Wavelet Packet Frame (DWPF); (c) Otsu method; (d) Morphological operations; (e) Left (IBL) and (f) Right (IBR) images of binary lungs obtained; (g) Enhancement technique of Contrast-limited adaptive histogram equalization (CLAHE); (h) Method of active contour Gradient Vector Flow (GVF); (i) Left (IACL ) and (j) Right (IBCR) segmented images of the lungs obtained; (k) Evaluation of the accuracy of the segmentation method; and (l) Image gold standard ( IG).  

Wavelet enhancement

Wavelet decomposition is an excellent operation to discriminate information, and it has been successfully used in other studies involving radiological images ( Vidaurrazaga et al., 2000 ; Katouzian et al., 2008 ). Therefore, we used DWPF to separate lung tissue from other image information. Specifically, for this project, we performed two DWPF decompositions from the original image ( Io) and approximation coefficient image (cA1), respectively ( Figure 2 ) ( Coifman and Wickerhauser, 1992 ). The wavelet chosen to perform the DWPF was the 4th derivative of Symlet ( Symlet4), because it was the wavelet that produced the highest correlation with the desired information ( Figure 2 b and c). Once we had the two levels of coefficients, we combined them to produce an image with good level of discriminated information (IRW). Therefore, the resultant image (IRW) can be obtained by:

IRW=cA1+cDh1+cDv1+cDd1+cDh2+cDv2+cDd2 (1)

where, cA1 is the approximation coefficient of the first decomposition, and cDh1, cDv1, cDd1, cDh2, cDv2 and cDd2 are all the detail coefficients from the first and second decomposition, respectively ( Figure 2 b, c and f). The mother wavelet, number of decompositions, and coefficient combination were defined during a parameter calibration stage, aiming at the best lung tissue discrimination as well as the highest final segmentation accuracy. In spite of having acquired a substantial amount of information, the result provided by the sum of the coefficients still has to be refined to obtain complete information of the pulmonary structure. One way is by binary morphology; hence, it becomes necessary to binarize the previous resultant image (IRW ).

Figure 2 Combination using the DWPF to obtain the wavelet enhancement technique, using the Symlet4 wavelet. (a) Original image; (b) 1st decomposition using the wavelet Symlet4; (c) Sub-images obtained from the approximation coefficients: cA1, and from the horizontal, vertical and diagonal detailing coefficients: cDh1, cDv1, and cDd1, referring to the 1st decomposition.; (d) 2nd decomposition using the wavelet Symlet4; (e) Sub-images obtained from the approximation and detailing coefficients referring to the 2nd decomposition; (f) The resulting image (IRW ) is obtained from the sum of sub-images cA1, cDh1, cDv1, cDd1, cDh2, cDv2 and cDd2.  


The Otsu method is an automatic threshold method ideal for bimodal histogram images ( Otsu, 1979 ). A highly bimodal histogram is produced by the resultant combination of wavelet coefficients, since the wavelet decompositions extract information by the highest correlation between the texture of the desired region and the mother wavelet ( Katouzian et al., 2008 ). Thus, the image highlighted by discrete wavelet (IRW) ( Figure 2 f), is binarized using Otsu (IRWB) ( Figure 3 a). Specifically, a binary image is obtained, with white for the pulmonary region and black for others. Once we have the binary image (IRWB) ( Figure 3 a), better accuracy of lung tissue detection can be obtained by applying a binary morphology to it.

Figure 3 Combination of binary operations for obtaining the binary objects of the left and right lungs. (a) Image enhanced by the DWPF method binarized by Otsu (IRWB ); (b) Selection phase; (c) Selection of the binary objects of the left and right lungs from points PL and PR, respectively; (d) Left (ISelL) and (e) Right (ISelR ) images of selected binary objects; (f) Opening operation applied to the pulmonary objects using variable dimensions of circular structure element; (g) Checking parameter of the opening operation; (h) Left (IOpenL) and (i) Right (I IOpenR) images of binary objects after opening; (j) closing operation with a circular structure element with radius of 30 pixels; and (k) Left (ISL) and (l) Right (ISR ) images of binary objects after closing.  

Binary morphology

Binary morphology is a set of morphological operations designed to improve the segmentation process, thus allowing a better separation between the region of interest, lung tissue, and the background of the image. At this stage, a binary image goes through a combination of binary morphological operations that aims to estimate and to enhance missing and/or incorrect information. Therefore, this phase is divided into three steps: selection, opening and closing. In the selection block ( Figure 3 in yellow), a selection of the representative objects of the left and right lungs in the binary image (IRWB) ( Figure 3 a) was performed to obtain separated images (ISelL and ISelR) ( Figure 3 d and e). Next, in the opening block ( Figure 3 in blue), each image selected ISelL and ISelR underwent a loop of the opening operation ( Figure 3 f), which has a check parameter to evaluate the operation ( Figure 3 g). This parameter performs the change in the size of the structure element until it is satisfied by some defined criteria, obtaining the IOpenL and IOpenR images for the left and right lungs ( Figure 3 h and i). Finally, in the closing block ( Figure 3 in red), the image went through a closure operation with fixed size structure element ( Figure 3 j), resulting the ISE and ISD images ( Figure 3 k and l).

The selection phase ( Figure 3 b) is to obtain the binary objects representing the left (ISelL ) and right (ISelR) lungs ( Figure 3 d and e). This step occurs to separate the two lungs in each image, so that they go through the rest of the process separately. This separation prevents the parameters adapted for improving one lung segmentation from harming the other lung segmentation. The selection of lung objects was performed automatically, selecting objects from seeds in pre-established locations. The location of the seeds was defined in the vertical position representing half of the image height, and the horizontal position by one third of the image width for the left lung ( PL) ( Figure 3 c), and two thirds of the image width for the right lung (PR) ( Figure 3 c) - position with high probability of being only lung tissue region.

Even after Otsu threshold and selection, the resultant binary image may present no desired components, with different size, connected to the representative lung region, which may decrease the final segmentation accuracy. Therefore, the second phase ( Figure 3 f) consisted in a loop of the opening operation applied in the pulmonary objects using variable dimensions of the circular structure element. Therefore, the following operation occurs for the left (ISelL) and right lungs (ISelR ) using a circular structure element according to:

IOpen=ISelropen (2)

where IOpen is the image after opening operation, ISel is the selected image representing ISelL and ISelR ( Figure 3 d), ropen is a variable radius size of the circular structure element. Initially, the radius size of the circular structure element was 1 pixel, but it may increase due to the requirement of each image. This requirement is determined in a check procedure, in which it is verified if the pulmonary object is connected to the edges of the image ( Figure 3 g). If the check procedure finds that the object is connected to the border, the image is re-opened with an additional unit in the radius size of the circular structure element. This loop occurs until the object is no longer connected to regions at the edge of the image, causing each selected image to have a different size of circular structure element, hence allowing the best aperture for each image. Finally, after the end of this loop process, the disengagement of these regions from the binarized lung tissue is performed, and the IOpenL images for the left lung and IOpenR for the right lung are obtained ( Figure 3 h and i).

The above operation isolated the pulmonary region from any background noise; nonetheless, missing information ( Figure 4 ) still needed to be estimated to construct a better approximation of the desired object. Therefore, the last phase consisted of a closing operation on each of the IOpenL and IOpenR images with a circular structure element with radius of 30 pixels ( Figure 3 j and Figure 4 ). This operation occurs for the left and right lungs ( Figure 4 ) according to:

IS=IOpenrclose (3)

where IS is the image that underwent the closing operation, IOpen is the previously opened image, and r close is the radius size of the closing circular structure element, rclose = 30. The possible discontinuities or missing information ( Figure 4 ), reason to carry out this phase, occur for two reasons. The first is that the opening operation in the previous phase may unnecessarily remove some possible border information. The second is due to the limitation of the Otsu thresholding method, which may cause the resultant image not to have the complete lung information. Consequently, after this opening operation, more complete binary images are obtained for the left ISL lung and for the right lung ISR ( Figure 3 k and l, Figure 4 ). However, these images still not have the highest accuracy in the lung tissue detection, mainly regarding the fine adjustments of border. Therefore, it is necessary to use some technique that allows increasing such a detail. In this work, we used the gradient vector flow technique, as a post processing phase, to improve the detailed information of the pulmonary border, and as a consequence obtained a high accurate pulmonary segmentation.

Figure 4 Closing procedure with example of the pulmonary missing information followed by the closing procedure to correct them.  

Active contour

Active Contour models were introduced by Kass, Witkin and Terzoupoulos in 1987 ( Candemir et al., 2014 ). The general idea of this method is to use a curve of energy minimization to extract important characteristics of an image. This energy associated with the curve is defined so that it is minimal when the curve lies on a border region with the characteristics to be extracted. In this way, the energy function starts to work as an objective function in an optimization problem, thus having different techniques to achieve this objective ( Cavalcante et al., 2013 ). One of these techniques is the Gradient Vector Flow (GVF) ( Xu and Prince, 1998 ). This method consists in searching the region of interest by detecting the lower energy region of the image. The method is carried out by an initial segment, following the image gradient. This segment must be selected near to the region of interest; hence, it will follow the vector pointing to the edge of this segment. Thus, in each interaction, the contour analyzes the image gradient, according to its positioning on image gradient; thus, each contour point is directed to the point of less energy, changing the contour coordinates.

In this work, the contour of the objects, obtained after the step of binary morphology ( ISL and ISR) ( Figure 3 k and l) were used as initial contour of GVF (segment), taking advantage of the morphological information previously obtained. Initializing the contour with the morphological information reduces the computational cost, since it decreases the number of interactions by having the initial contour very close to the final object. Consequently, this increases the accuracy of lung tissue detection from the detailed detection of the pulmonary borders. Accordingly, this phase was divided into 5 steps: border detection, gradient acquisition, contour initialization, contour processing and filling ( Figure 5 ). In the border detection block ( Figure 5 in yellow), the border was obtained to the resulting objects of the checked for possible displacement towards regions of lower energy. When all the points find the region of lower energy, they form a more accurate border, forming the new contour of objects of segmentation for the left and right lungs (IACL and IACR ) ( Figure 5 m and n). Finally, in the filling block ( Figure 5 in green), the contours of images IACL and I ACR are filled up ( Figure 5 o) so that the objects and area of segmentation of the left and right lungs (I FL and IFR) ( Figure 5 p and q) are obtained.

Figure 5 Methodology of the GVF Active Contour method used in this project for detailing the edge of the chest x-ray segmentation. (a) Left (ISL) and (b) Right (ISR) images of binary objects obtained after the binary morphology; (c) Canny edge detection technique; (d) Left (IedgeL) and (e) Right (IedgeR ) lung border of binary object obtained by Canny to be used as initial contour in the GVF; (f) Image highlighted by the technique of Contrast-limited adaptive histogram equalization (CLAHE); (g) Obtaining the gradient of the image; (h) Image gradient enhanced by the CLAHE method; (i) Initialization of the GVF method; (j) Left (IiniL ) and (k) Right (IiniR) image gradient enhanced by CLAHE initialized by the edge of the binary object; (l) Processing of the GVF method; (m) Left (IACL) and (n) Right (I ACR) images of the final contour obtained by the GVF method; (o) Contour fills; and; (p) Left (IFL) and (q) Right (I FR) images of the segmented lungs.  

In the border detection phase, the edge of the binary objects related to the left and right lungs were obtained. The edges of each lung were acquired by the Canny technique ( Canny, 1987 ) ( Figure 5 c), which located every border pixel position for the left (IedgeL ) and right lungs (IedgeR) ( Figure 5 d and e). This step is performed due to the need of the GVF method contour initialization; thus, the vector of border pixel position replaces the default initial manual contour of deformable models. Consequently, since the initial contour is automatically inserted very close to the desired final border, it removes the need of human interaction, increasing the final accuracy; moreover, it reduces the computational cost, because it requires a smaller number of interactions. In addition to the initial vector position, the GVF method requires the image gradient to direct the deformable contour through the desired final contour. Accordingly, in the gradient acquisition phase, we used the gradient of the image enhanced by the CLAHE method ( Ikhsan et al., 2014 ) (IRC) ( Figure 5 f). The CLAHE method is inspired by the histogram equalization (HE), which works by adjusting the image contrast, improving the appearance of the image by scaling the intensity interval by redistributing pixel values ( Ikhsan et al., 2014 ). For this, the histogram is extended to produce a broader uniform distribution, obtaining several histograms, one in each distinct part of the image, using an exponential interpolation technique between each sub histogram, resulting in the enhanced image (IRC ). We chose the image enhanced by the CLAHE method to be the source of the gradient matrix creation, as its enhanced pixel intensities produces gradient vectors with greater magnitudes, comparing to the original intensities; hence, enforcing the active contours shifting, in addition to preventing getting stacked in local minimums. Additionally, the use of CLAHE also presented highest final segmentation accuracy. Thus, the image gradient enhanced by the CLAHE method ( Figure 5 g) was obtained according to:

IRCG(i,j)=[dIRC(i,j)dx,dIRC(i,j)dy] (4)

where IRCG(i,j) is the gradient vector in the position (i,j), dIRC(i,j)dx is the derivative with respect to x (gradient in horizontal), dIRC(i,j)dy is the derivative with respect to y (gradient in vertical), and (i,j) is the respective position of each point in the vertical (i) and horizontal (j) position of the image. Thus, the gradient image (IRCG) ( Figure 5 h), which will be used as the basis for the GVF method process, was obtained.

Next, in the contour initialization phase ( Figure 5 i), the binary border obtained previously had its positions transformed to the deformable contour. The contours stand over the gradient images of the left and right lungs, IiniL and IiniR ( Figure 5 j and k and Figure 6 a and b); hence, the deformable process is ready to begin. In this phase, the points of the deformable contours over the images IedgeL and IedgeR are stored in a row vector (VR) ( Figure 6 c), so that each point can be analyzed for the possibility of displacement as a function of the energy of the gradient.

Figure 6 Methodology of the process of GVF Active Contour. (a) Left (IiniL ) and (b) Right (IiniR) image gradient enhanced by CLAHE initialized by the edge of the binary object; (c) Row vector (V R) whereby the points of edge of the binary object are ordered clockwise; (d) Displacement Procedure, whereby the points of the row vector are analyzed and displaced towards lower energy pixels; (e) Vector for future analyses (VFA) whereby the points that will still be analyzed are stored; (f) Interpolation procedure among the points present in the row vector; (g) GVF process shutdown checkpoint (Stop Condition ); and (h) Left (IACL) and (i) Right (I ACR) images of the final contour obtained by the GVF method.  

Afterwards, the contour processing phase is divided into 3 sub phases, Displacement Procedure, Interpolation and Stop Condition ( Figure 6 ). The first interaction of the Displacement Procedure occurs from the edge positions of the obtained contour, results of the previous step (IiniL and IiniR) ( Figure 6 a and b) and stored in the row vector (VR) ( Figure 6 c). In the Displacement Procedure, the points of the row vector are analyzed and displaced towards lower energy pixels ( Figure 6 d). A shift of 1 pixel per interaction was chosen to increase accuracy. During each interaction, the position of displaced points are stored in a vector for future analyses (V FA) ( Figure 6 e) and the value position of the same points are changed in VR , so that they are the only ones to be verified in the next interaction, together with the points added in the interpolation phase, maximizing the computational cost reduction. After that, an Interpolation and Reorder procedure among the points presented in the row vector (VR) was performed ( Figure 6 f). During this interpolation procedure, new points may be added in (VR ) to satisfy the interpolation function. The interpolation procedure used was Spline, since it presented the best result among the techniques tested. In addition, as the points reach the minimum energies, VFA continuously decreases its size at each interaction, showing fewer points as its positions keeps being displaced, and then, when it reaches a total size of zero, a Stop Condition is activated. ( Figure 6 g). Finally, when the condition is eventually satisfied, the final segmented border of the left and right lungs, IACL and IACR , respectively ( Figure 5 m and n and Figure 6 h and i) are obtained.

Lastly, in the filling phase ( Figure 5 o), the borders obtained in the previous phase underwent a filling procedure, obtaining the left (IFL) and right (IFR ) segmented lungs ( Figure 5 p and q). This filling is performed so that the segmented images can be evaluated for the detection of lung tissue. These are the final images of the pulmonary segmentation to be evaluated.


The proposed approach was evaluated by segmenting and computing the parameters of the lung tissue region and border detection. A number of 247 chest X-ray images from JSRT with different intensities and anatomical variations of thorax structures were used. The images are in digital format, and have resolution of 2048 x 2048 pixels (resolution of pixel: 0.175 mm/pixel) ( Shiraishi et al., 2000 ). Additionally, experts from the mentioned database kindly produced and made available the gold standards of the lung region for evaluating image segmentation ( Shiraishi et al., 2000 ). The averaged time of segmentation of the lung region using the software and computer described above was 7.03 s per image. This processing time is much faster than manual segmentation, which is well known to be more than one minute per image. It is also more practical and much less exhaustive, avoiding the possibility of error and variation of inter and intra operator. Additionally, it was also possible to verify an appreciable high correlation between the images segmented by the method presented in this work from the calculation of the following metrics:

Assessing the accuracy of the segmented area

This evaluation was performed to evaluate the accuracy of the segmentation method regarding the distinction of the lung tissue region. This comparison was carried out computing the averaged and standard deviation of the measurements of accuracy proposed by Udupa et al. (2006): “True Positive” (TP), “False Positive” (FP) and “False Negative” (FN), and the parameters of Overlap Dice (OD) and Overlap (OR ) indices ( Dice, 1945 ; Kupinski and Giger, 1998 ). The averaged, standard deviation (std) and maximum (max) and minimum (min) of TP, FP, FN, OD and OR indices obtained between the comparisons of the 247 images segmented with their corresponding gold standards are presented in Table 1 .

Table 1 Values of average, standard deviation and maximum and minimum of accuracy of area segmentation and edge detection of the proposed segmentation method, which are the True Positive (TP), False Positive (FP), False Negative (FN), Overlap Dice (OD), Overlap (OR), Average Contour Distance (ACD) and maximum error (Emax).  

TP (%) FP (%) FN (%) OD (%) OR (%) ACD (mm) Emax (mm)
Average 97.83 0.007 2.17 98.87 97.63 0.69 3.27
std ±3.54 ±0.0016 ±3.54 ±1.91 ±3.34 ±0.95 ±2.8
min 72.66 0 0.12 84.14 72.46 0.08 0.55
max 99.88 0.35 27.34 99.94 99.65 6.05 10.46

Table 1 shows high TP average (approximately 98%), with a low FP in the order of 7x10-3%. Additionally, when the combination of three parameters are concerning, represented by the Overlap (OR), the outcome of accuracy is still amazingly high, also very close to 98%. The very small standard deviation demonstrates the method precision, since it shows how accuracy changes very little with all the different 247 images. Moreover, only the number of three images have OR smaller than 80%; finally, to conclude the accuracy and precision analyses of the lung area, an appreciable number of 80% of the images have OR greater than 95%, definitely corroborating the method huge accuracy.

Assessing the accuracy of the corresponding edge

This evaluation was performed to compute the accuracy of the corresponding edge of the lung area. In this evaluation, the positions of the edges of the segmented objects (β(S)) and gold standards (β(G)) were compared, where both edges of each of the 247 images were performed by the Canny method ( Ikhsan et al., 2014 ). In the evaluation, the averaged and standard deviation of the Average Contour Distance (ACD), as calculated by Xiong et al. ( Candemir et al., 2014 ; Xiong et al., 2017 ), and the maximum error ( Emax) between the edges of the segmented images and their Gold Standards were calculated. The values of averaged, standard deviation (std), maximum (max) and minimum (min) of the ACD and Emax indices obtained between the comparisons of edges of the 247 images segmented with their edge gold standards are presented in Table 1 (where the accuracy, precision and robustness f the method are confirmed). The averaged of maximum false positive error ( EmaxFP ), slightly higher than 3mm, produces a visually undetectable maximum error. In addition to that, the averaged value of ACD is close to 0.7 mm, which represents an insignificant impact of less than 0.2% of average error in the final segmented object, considering the image size of (35.84 x 35.84 cm) ( Shiraishi et al., 2000 ).


As presented in the introduction, there is a large demand for tools and methods that support the diagnosis of pulmonary abnormalities, increasing the probability of patient recovery ( Zhu et al., 2008 ). These diagnoses must be fast, accurate, precise and practical, given the urgency, and high possibility of death caused by the numerous pathologies associated with the lung ( Priya et al., 2017 ). A rapid and non-invasive diagnostic solution for initiating treatment is chest radiography ( Michel-González et al., 2011 ). However, techniques that can quantitatively identify the different abnormalities of the lung region that can hardly be seen with the naked eye are needed. In addition to that, another demand to obtain quantitative results is the possibility of numerically comparing the pathology evolution or treatment progression, overcoming many subjective interpretations that could result image quality of different equipment and operator expertise. Thus, the use of segmentation techniques becomes a fundamental step for performing more accurate diagnoses regarding the abnormalities present in the image ( Sokashe, 2013 ; Udupa et al., 2006 ). Moreover, the use of automatic segmentation techniques brings advantages by not requiring specialized operators, time optimization compared to manual segmentation, and zero vulnerability of intra and inter-operators interpretations. The few automatic segmentation methods in the literature, completely or partially presented limitations to be overcome. They have not achieved a full potential of high accuracy or they depend on consecutive exams or images of the same patient to increase accuracy. In both cases, it could lead to a higher patient risk, by either the misinterpretation of a low accurate method or a higher incidence of patient exposure, as the method depends on many acquisitions; thus, there is stillroom for the growth of the efficacy of such methods.

The methodology presented herein consists of a combination of enhancement techniques and binary morphology to obtain the lung binary object; hence, the object contour is used to initialize the active contour method of GVF, allowing the final fine contour adjustment. Accordingly, this specialized methodology has overcome the limitations of the literature practicability and accuracy, presenting high accuracy and practicability. The accuracy presented by our methodology could be corroborated by the analysis of the overlap (OR), which was close to 98% ( Table 1 ). This result is almost 6% above the standard automatic segmentation methods that used the same database in their studies ( Dai et al., 2017 ), and 2% higher than the best methods found in the literature, using these or other images ( Shi et al., 2008 ; Xiong et al., 2017 ). Moreover, in order to reach the overlap of 96%, the highest found in literature ( Shi et al., 2008 ), Shi and colleagues presented a methodology, which depends on a series of examinations of the same patient region. Accordingly, a method which exposes the patients to a greater amount of radiation and expend more time and cost to obtain the diagnosis. On the other hand, our proposed method requires a single image of the patient to perform the segmentation with superior accuracy. Another recent study, using the same database for evaluation, was presented by Xiong et al. ( Dai et al., 2017 ; Xiong et al., 2017 ). Xiong study presented an average OR of 95.5%, while our proposed method presented almost 98%, with a margin of accuracy 2% higher. In addition, the parameter of edge accuracy, ACD ( Table 1 ), of our proposed approach, in which averaged, deviation, minimum and maximum values obtained were 0.69, 0.95, 0.08 and 6.05, respectively, are better or similar to the best results of the outcome presented by Xiong et al. ( Dai et al., 2017 ; Xiong et al., 2017 ), (1.43, 0.74, 0.68 and 5.96, respectively). In addition, it is possible to evaluate with more detail the overcoming of the presented method by comparing the averaged values of OR and ACD of others lung segmentation methods of literature using the same database as shown in Table 2 . The values in Table 2 confirm the high accuracy of the proposed method in the identification of the area and the pulmonary borders. Thus, the obtainment of a new automatic methodology of pulmonary tissue segmentation that surpasses the accuracy present in the current literature is fundamental to serve as a basis for the development of techniques that improve the identification and quantification of pulmonary pathologies, helping to obtain the diagnosis and speeding up initiation of treatment.

Table 2 Comparison between our method and other methods present in the literature in terms of averaged values of overlap (OR) and average contour distance (ACD).  

Method OR ± std (%) ACD ± std (mm)
Hybrid voting ( Van Ginneken et al., 2006 ) 94.90±2.00 1.620±0.660
PC post-processed ( Van Ginneken et al., 2006 ) 94.50±2.20 1.61±0.80
Hybrid ASM/PC ( Van Ginneken et al., 2006 ) 93.40±3.70 2.08±1.40
Hybrid AAM/PC ( Van Ginneken et al., 2006 ) 93.30±2.60 2.06±0.84
ASM tuned ( Van Ginneken et al., 2006 ) 92.70±3.20 2.30±1.03
ASM default ( Van Ginneken et al., 2006 ) 90.30±5.70 3.23±2.21
Mean shape ( Van Ginneken et al., 2006 ) 71.30±7.50 10.06±3.18
Shi et al. (2008) 92.00±3.10 1.78±0.78
Ibragimov et al. (2012) 94.60±2.20 1.59±0.68
Shao et al. (2014) 94.60±1.90 1.67±0.76
Candemir et al. (2014) 95.40±1.50 1.32±0.32
Wu et al. (2015) 94.70±2.10 NA *
Xiong et al. (2017) 94.40±2.10 1.43±0.74
Dai et al. (2017) 94.7±0.40 NA *
Our Method 97.63±3.34 0.69±0.95

*NA: Not Available.

The need for a great number of adaptive blocks, such as the adaptive construction of structuring elements ( Figure 4 (f) and (g)) and the use of GVF ( Figure 6 ), may be seen as a limitation of the method regarding computational cost. Nonetheless, they were thought to start very close to the final desired parameter value, which required very few adjustments, therefore keeping the method with low computational cost. However, futures works may focus on alternatives for these operations as well as including new images for evaluations and adapting new versions of this method to other applications, such as ribs and heart segmentation. The values here presented emphasize the outstanding accuracy of the method proposed in this study, given the efficiency of discrimination of edge lung tissue with loss below 0.2% of average error in the final segmented image.

In Conclusion, the proposed approach presented significant values regarding accuracy and practicability for an automatic chest X-ray segmentation. The methodology of the proposed approach is based on a unique combination of operations devoted to lung region segmentation. As well as the outcome, this methodology presented many other specific contributions: (i) The investigation of an enhancement method, which successfully permits Otsu’s threshold to discriminate lung tissue ( Figure 1 b and d); (ii) A combination of adaptive binary morphological operations devoted to estimating and to correcting the morphological misinterpretation of lung anatomy ( Figure 1 e); (iii) The use of an adapted automatic version of GVF as post-processing with specific features for fine final contour adjustments ( Figure 1 h); (iv) The automatic initialization of GVF with a very close lung object obtained previously, excluding the need of an operator and allowing fine adjustment and low computational cost, since the initialized contour is close to the fine object ( Figure 5 d and e); (v) The use of a gradient vector from an enhanced image, increasing the absolute values of vectors, amplifying the contour shifting and preventing GVF from getting stacked in local minimums ( Figure 5 g, h and i); (vi) A specialized routine to manage the necessary number of points used for interpolation ( Figure 6 ) - hence, increasing the efficiency of GVF, while reducing computational cost. Consequently, with the combination of the operations mentioned, we obtained an automatic segmentation method with exceptional accuracy and practicability. As a result, we increase the precision and practicability of the lung tissue segmentation, improving the diagnostic for pulmonary tissue abnormalities; hence, increasing the patient's chance of survival and recovery.

Finally, for future works, we will seek new collaborators to increase our database to expand the analyses of the potential applications. In addition to seek for methodological improvements for segmentation accuracy, techniques for discriminating pulmonary abnormalities, such as cancer, tuberculosis and pneumonia will be also in the scope of future works.


The authors are grateful to the Coordination for the Improvement of Higher Education Personnel (CAPES-Brazil) and the Laboratory of Image and Signal Processing of the Institute of Science and Technology of UNIFESP (LaPIS-ICT-UNIFESP).

How to cite this article: Vital DA, Sais BT, Moraes MC. Robust pulmonary segmentation for chest radiography, combining enhancement, adaptive morphology and innovative active contours Res Biomed Eng. 2018; 34(3):. DOI: 10.1590/2446-4740.180035.


Armato SG 3rd, Giger ML, MacMahon H. Automated lung segmentation in digitized posteroanterior chest radiographs. Acad Radiol. 1998; 5(4):245-55. . PMid:9561257. [ Links ]

Candemir S, Jaeger S, Palaniappan K, Musco JP, Singh RK, Zhiyun Xue, Karargyris A, Antani S, Thoma G, McDonald CJ. Lung segmentation in chest radiographs using anatomical atlases with nonrigid registration. IEEE Trans Med Imaging. 2014; 33(2):577-90. . PMid:24239990. [ Links ]

Canny J. A computational approach to edge detection. USA: Readings Comp Vision, 1987. p. 184-203. [ Links ]

Cavalcante TS, Cortez PC, Almeida TM, Felix JHS, Holanda MA. Segmentação automática 2D de vias aéreas em imagens de tomografia computadorizada do tórax. Rev Bras Eng Bioméd. 2013; 29(4):389-403. . [ Links ]

Chaudhary C, Patil MK. Review of image enhancement techniques using histogram equalization. Int J Applic Innovation Eng Managem. 2013; 2(5):343-9. [ Links ]

Chen S, Suzuki K. Separation of bones from chest radiographs by means of anatomically specific multiple massive-training ANNs combined with total variation minimization smoothing. IEEE Trans Med Imaging. 2014; 33(2):246-57. . PMid:24132005. [ Links ]

Coifman RR, Wickerhauser MV. Entropy-based algorithms for best basis selection. IEEE Trans Inf Theory. 1992; 38(2):713-8. . [ Links ]

Costa CFF Fo, Levy PC, Xavier CDM, Fujimoto LBM, Costa MGF. Automatic identification of tuberculosis mycobacterium. Res Biomed Eng. 2015; 31(1):33-43. . [ Links ]

Dai W, Doyle J, Liang X, Zhang H, Dong N, Li Y, Xing EP. SCAN: structure correcting adversarial network for chest x-rays organ segmentation. ArXiv. 2017:1-10. [ Links ]

Departamento de Informática do Sistema Único de Saúde – DATASUS. Informações de Saúde [Internet]. Brasília: Ministério da Saúde; 2015. [cited 2017 Feb 10]. Available from: [ Links ]

Dice LR. Measures of the amount of ecologic association between species. Ecology. 1945; 26(3):297-302. . [ Links ]

Ferlay J, Shin HR, Bray F, Forman D, Mathers C, Parkin DM. Estimates of worldwide burden of cancer in 2008: GLOBOCAN 2008. Int J Cancer. 2010; 127(12):2893-917. . PMid:21351269. [ Links ]

Ibragimov B, Likar B, Pernus F, Vrtovec T. A game-theoretic framework for landmark-based image segmentation. IEEE Trans Med Imaging. 2012; 31(9):1761-76. . PMid:22692901. [ Links ]

Ikhsan IAM, Hussain A, Zulkifley MA, Tahir NM, Mustapha A. An analysis of x-ray image enhancement methods for vertebral bone segmentation. In: Proceedings of the IEEE 10th International Colloquium on Signal Processing and its Applications; 2014 Mar 7-9; Kuala Lumpur, Malaysia. USA: IEEE; 2014. p. 208-11. [ Links ]

Katouzian A, Baseri B, Konofagou EE, Laine AF. Texture-driven coronary artery plaque characterization using wavelet packet signatures. In: Proceedings of the 5th IEEE International Symposium on Biomedical Imaging: From Nano to Macro; 2008 May 14-17; Paris, France. USA: IEEE; 2008. p. 197-200. . [ Links ]

Kumar SV, Deka MK, Bagga M, Kala MS, Gauthaman K. A systematic review of different type of tuberculosis. Eur Rev Med Pharmacol Sci. 2010; 14(10):831-43. PMid:21222369. [ Links ]

Kupinski MA, Giger ML. Automated seeded lesion segmentation on digital mammograms. IEEE Trans Med Imaging. 1998; 17(4):510-7. . PMid:9845307. [ Links ]

Lorgelly PK, Atkinson M, Lakhanpaul M, Smyth AR, Vyas H, Weston V, Stephenson T. Oral versus iv antibiotics for community-acquired pneumonia in children: a cost-minimisation analysis. Eur Respir J. 2010; 35(4):858-64. . PMid:19717479. [ Links ]

McNitt-Gray MF, Huang HK, Sayre JW. Feature selection in the pattern classification problem of digital chest radiograph segmentation. IEEE Trans Med Imaging. 1995; 14(3):537-47. . PMid:18215858. [ Links ]

Michel-González E, Cho MH, Lee SY. Geometric nonlinear diffusion filter and its application to X-ray imaging. Biomed Eng Online. 2011; 10(1):47. . PMid:21639933. [ Links ]

Moraes MC, Furuie SS. An approach to automatically segment the media-adventitia borders in IVUS. Rev Bras Eng Bioméd. 2010; 26(3):219-33. . [ Links ]

Otsu N. A threshold selection method from gray-level histograms. IEEE Trans Syst Man Cybern. 1979; 9(1):62-6. . [ Links ]

Priya CL, Gowthami D, Poonguzhali S. Lung pattern classification for interstitial lung diseases using an ANN-back propagation network. In: Proceedings of the 2017 International Conference on Communication and Signal Processing (ICCSP); 2017 Apr 6-8; Chennai, India. USA: IEEE; 2017. p. 1917-22. . [ Links ]

Rebouças PP Fo, Cortez PC, Félix JHS, Cavalcante TS, Holanda MA. Modelo de Contorno Ativo Crisp Adaptativo 2D aplicado na segmentação dos pulmões em imagens de TC do tórax de voluntários sadios e pacientes com enfisema pulmonar. Rev Bras Eng Bioméd. 2013; 29(4):363-76. . [ Links ]

Ruuskanen O, Lahti E, Jennings LC, Murdoch DR. Viral pneumonia. Lancet. 2011; 377(9773):1264-75. . PMid:21435708. [ Links ]

Shao Y, Gao Y, Guo Y, Shi Y, Yang X, Shen D. Hierarchical lung field segmentation with joint shape and appearance sparse learning. IEEE Trans Med Imaging. 2014; 33(9):1761-80. . PMid:25181734. [ Links ]

Shi Y, Qi F, Xue Z, Chen L, Ito K, Matsuo H, Shen D. Segmenting lung fields in serial chest radiographs using both population-based and patient-specific shape statistics. IEEE Trans Med Imaging. 2008; 27(4):481-94. . PMid:18390345. [ Links ]

Shiraishi J, Katsuragawa S, Ikezoe J, Matsumoto T, Kobayashi T, Komatsu KI, Matsui M, Fujita H, Kodera Y, Doi K. Development of a digital image database for chest radiographs with and without a lung nodule: receiver operating characteristic analysis of radiologists’ detection of pulmonary nodules. AJR Am J Roentgenol. 2000; 174(1):71-4. . PMid:10628457. [ Links ]

Sokashe MSV. Computer assisted method for cervical vertebrae segmentation from x-ray images. Computer. 2013; 2(11):4387-8. [ Links ]

Stegmann MB, Ersboll BK, Larsen R. FAME-a flexible appearance modeling environment. IEEE Trans Med Imaging. 2003; 22(10):1319-31. . PMid:14552585. [ Links ]

Udupa JK, Leblanc VR, Zhuge Y, Imielinska C, Schmidt H, Currie LM, Hirsch BE, Woodburn J. A framework for evaluating image segmentation algorithms. Comput Med Imaging Graph. 2006; 30(2):75-87. . PMid:16584976. [ Links ]

Van Ginneken B, Stegmann MB, Loog M. Segmentation of anatomical structures in chest radiographs using supervised methods: a comparative study on a public database. Med Image Anal. 2006; 10(1):19-40. . PMid:15919232. [ Links ]

Vidaurrazaga M, Diago LA, Cruz A. Contrast enhancement with wavelet transform in radiological images. In: Proceedings of the 22nd Annual International Conference of the IEEE Engineering in Medicine and Biology Society; 2000 July 23-28; Chicago, IL. USA: IEEE; 2000. p. 1760-3. . [ Links ]

Wan Ahmad WSHM, W Zaki WMD, Ahmad Fauzi MF. Lung segmentation on standard and mobile chest radiographs using oriented Gaussian derivatives filter. Biomed Eng Online. 2015; 14(1):20. . PMid:25889188. [ Links ]

Wu G, Zhang X, Luo S, Hu Q. Lung segmentation based on customized active shape model from digital radiography chest images. J Med Imaging Health Inf. 2015; 5(2):184-91. . [ Links ]

Xiong J, Shao Y, Ma J, Ren Y, Wang Q, Zhao J. Lung field segmentation using weighted sparse shape composition with robust initialization. Med Phys. 2017; 44(11):5916-29. . PMid:28875551. [ Links ]

Xu C, Prince JL. Snakes, shapes, and gradient vector flow. IEEE Trans Image Process. 1998; 7(3):359-69. . PMid:18276256. [ Links ]

Zhu H, Sun W, Wu M, Guan G, Guan Y. Pre-processing of X-ray medical image based on improved temporal recursive self-adaptive filter. In: Proceedings of the 2008 The 9th International Conference for Young Computer Scientists; 2008 Nov 18-21; Hunan, China. USA: IEEE; 2008. p. 758-63. . [ Links ]

Received: May 14, 2018; Accepted: August 14, 2018

*Corresponding author: Daniel Aparecido Vital, Laboratory of Image and Signal Processing, Institute of Science and Technology, Federal University of São Paulo, Rua Talim, 330, Sala 108, Jardim Aeroporto, CEP 12231-280, São José dos Campos, SP, Brazil. E-mail:

Creative Commons License  This is an Open Access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.