Abstract
3D volume segmentation is the process of partitioning voxels into 3D regions (subvolumes) that represent meaningful physical entities which are more meaningful and easier to analyze and usable in future applications. Multiresolution Analysis (MRA) enables the preservation of an image according to certain levels of resolution or blurring. Because of multiresolution quality, wavelets have been deployed in image compression, denoising, and classification. This paper focuses on the implementation of efficient medical volume segmentation techniques. Multiresolution analysis including 3D wavelet and ridgelet has been used for feature extraction which can be modeled using Hidden Markov Models (HMMs) to segment the volume slices. A comparison study has been carried out to evaluate 2D and 3D techniques which reveals that 3D methodologies can accurately detect the Region Of Interest (ROI). Automatic segmentation has been achieved using HMMs where the ROI is detected accurately but suffers a long computation time for its calculations.
1. Introduction
Volume segmentation allocates the voxels in 3D images into partitions or 3D regions that represent meaningful physical entities. The goal is to distinguish between different regions in the 3D volume and cover the extracted contours from the entire volume. Voxels' classification into regions is performed according to a certain region to which the voxels belong, and some shared, predefined properties. Those voxels comprise an isolated or segmented Object Of Interest (OOI) from the input volume.
There are many existing techniques used for medical image segmentation, including Multiresolution Analysis (MRA), statistical methods, and thresholding- and clustering-based techniques. Clustering technique aims to classify each pixel in an image into the proper cluster, and then these clusters are mapped to display the segmented images. A certain clustering criterion can be adopted to group each pixel into a specific number of clusters depending on the image histogram [1, 2]. Medical images can also be segmented using thresholding approaches by partitioning their intensities. When images contain different structures with contrasting intensities, thresholding provides a simple but effective means for obtaining segmentation. Generally, the thresholds are generated based on visual assessment of the resulting segmentation [3, 4].
MRA allows the preservation of an image according to certain levels of resolution. Consequently, wavelets have been useful in image compression, de-noising, and classification. Wavelet theory, which is built on solid mathematical foundations uses well-established tools such as quadrature mirror filtering, subband coding, and pyramidal image processing. Wavelet analysis enables the exploitation of signal or image characteristics associated with a particular resolution level, which may not be detected using other analysis techniques [5, 6].
Statistical modeling is a set of mathematical equations which describes the behavior of an object of study in terms of random variables and the associated probability distribution. Markov Random Field Model (MRFM) is a statistical approach which has been utilized within segmentation methodologies to model spatial interactions between neighbor pixels [1, 3, 7]. These local correlations between pixels are used as a mechanism for modeling various image properties. From a medical imaging perspective, this kind of analysis is useful, as most pixels can be categorized into the same class as their neighbors [1].
Statistical models using Hidden Markov Models (HMMs) observe a sequence of emissions with a hidden sequence of states that the model went through to generate the emissions [8, 9]. HMM states are not directly visible to the observer, but the variables influenced by the states are visible. Each state in HMM has a probability distribution over the other states which evaluate the state sequence. The challenge in HMM is to determine the hidden parameters from the observable parameters to be used in performing further analysis [8, 9].
This paper focuses on the implementation of efficient and robust medical volume segmentation techniques. MRA including wavelet and ridgelet transforms have been deployed for feature extraction, while statistical modeling using HMMs has been used for segmentation. The outline of this paper is as follows. In the following section, the proposed segmentation system is illustrated and discussed. In Section 3, the mathematical background of the proposed segmentation methods is presented with some test images. Section 4 presents the results and analysis of the proposed techniques for medical volumes segmentation. Conclusions and implications for future work are discussed in Section 5.
2. Proposed Segmentation System
In medical applications, the source of the 3D data set is the acquisition systems such as Positron Emission Tomography (PET), Computerized Tomography (CT), or Magnetic Resonance Imaging (MRI). Such devices are capable of slicing an object in a physical sectioning. 3D data set from those devices can be considered as parallel slices stacked to form a 3D volume. A segmented medical volume into sub-volumes which are more meaningful and easier to analyze and understand is the output the proposed system illustrated in Figure 1.
Hybrid multiresolution statistical approaches and other segmentation techniques are used to achieve accurate segmented volumes. System input is a 3D phantom or real volume from scanner acquisition. Acquisition systems produce a number of 2D slices resulted from the scanned volume of the body. These slices can be individually segmented using 2D segmentation methods such as thresholding, clustering, and HMMs followed by volume reconstruction or directly using 3D segmentation methods such as 3D that thresholding or 3D discrete wavelet transform (3D-DWT) after volume reconstruction process.
This paper explains the new application of wavelet transform directly on the 3D medical volumes from the acquisition systems using 3D-DWT with Haar wavelet filter. HMMs have been also applied on those volumes slice-by-slice to segment the Region Of Interest (ROI) into a number of classes based on the grey scale values of the original volume pixels.
3. Segmentation Methods and Their Mathematical Backgrounds
The mathematical background of the developed techniques for 3D medical volume segmentation system is presented in this section.
3.1. Thresholding
Scalar images can be segmented using thresholding approaches by partitioning image intensities. This methodology attempts to determine an intensity value that can separate the signal into a desired number of classes. Segmented images can be achieved by clustering all pixels with intensities larger than the threshold value into one class, and all others into another. In many applications, the threshold values selection can be done depending on the basis of histogram. Multithresholding occurs when more than one threshold value is determined [10].
The voxels of a certain object are not necessarily connected after thresholding because this technique does not consider the spatial characteristics of an image, thus causing it to be sensitive to noise and intensity fluctuations. For this reason it cannot be easily applied to many medical imaging modalities. These drawbacks essentially corrupt the histogram of the image-making partitioning via the selection of more problematic appropriate thresholds [11].
Hard thresholding technique is a boolean filter [10] which depends on pixel value and threshold value. As illustrated in Algorithm 1, it either makes the input zero or keeps it without any changes [12].
| ||||||
Hard thresholding process is less complex than soft thresholding (Algorithm 2), where pixels that have values greater than threshold value do not change [13]; soft thresholding replaces each pixel which has greater value than the threshold value by the difference between threshold value and pixel value [10]. This makes the process more complicated and increases the processing time for the algorithm. Figure 2 illustrates a segmented medical brain slice from an MRI scanner using hard and soft thresholding techniques. Many other thresholding types are widely used in different areas including adaptive thresholding.
| ||||||||||
(a)
(b)
(c)
It can be seen here that applying thresholding techniques is a very easy process and can be affected easily by surrounding noise, but it has been used as a preprocessing step and postprocessing step with other segmentation techniques. It is worth mentioning that thresholding can be replaced by the clustering technique which will be explained in Section 3.2.
3D thresholding method is similar to the 2D approaches where thresholding process is applied on all pixels in the volume instead of that in the plane. Algorithm 3 explains the pseudocode for 3D thresholding.
3.2. Clustering
Clustering technique is the process of classifying each group of pixels in an image into one class, each class has the same or similar properties which evaluate a specific part of an image. Each class is highlighted in the segmented image to illustrate the image as a number of separated regions, and the ROI may be one of those regions. Clustering technique is based on multithreshold values which can be set depending on the image histogram. Amira et al. used in [11] -means method for segmenting medical volumes; it is the most commonly used clustering technique. -means will be used in the fourth section of this research paper as a previously available work to validate the proposed techniques in the comparison tables.
-means clustering classifies voxels into clusters or classes ( less than ). This algorithm chooses the number of clusters () then randomly generates clusters and determine the cluster centers. The next step is assigning each point in the volume to the nearest cluster center and finally recompute the new cluster centers. The two previous steps are repeated until the minimum variance criterion is achieved. This approach is similar to the expectation-maximization algorithm for Gaussian mixture in which they both attempt to find the centers of clusters in the volume. Its main objective is to achieve a minimum intracluster variance : where is the number of clusters, and is the mean of all voxels in cluster . Figure 3 illustrates an example of segmenting a real CT slice by classifying the image into 6 different classes.
(a)
(b)
3.3. Wavelet Transform
In the last decade, wavelet transform has been recognized as a powerful tool in a wide range of applications, including image/video processing, numerical analysis, and telecommunication. The advantage of wavelet over existing transforms such as Discrete Fourier Transform (DFT) and Discrete Cosine Transform (DCT) is that wavelet performs a multiresolution analysis of a signal with localization in both time and frequency [14, 15]. In addition to this, functions with discontinuities and functions with sharp spikes require fewer wavelet basis vectors in the wavelet domain than sine-cosine basis vectors to achieve a comparable approximation. Wavelet operates by convolving the target function with wavelet kernels to obtain wavelet coefficients representing the contributions in the function at different scales and orientations. Wavelet or multiresolution theory can be used alongside segmentation approaches, creating new systems which can provide a segmentation of superior quality to those segmentation approaches computed exclusively within the spatial domain [16].
3.3.1. Discrete Wavelet Transform
Discrete wavelet transform (DWT) can be implemented as a set of filter banks, comprising a high-pass and low-pass filters. In standard wavelet decomposition, the output from the low-pass filter can then be decomposed further, with the process continuing recursively in this manner. According to [17], DWT can be mathematically expressed by(2) The coefficients and refer to approximation and detailed components in the signal at decomposition level respectively. The and represent the coefficients of low-pass and high-pass filters, respectively.
DWT decomposes the signal into a set of resolution-related views. The wavelet decomposition of an image creates at each scale a set of coefficient values , with an overall mean of zero. This set of coefficient values contains the same number of voxels as the original 3D volume, and therefore, this wavelet transform is redundant [18, 19]. A nondecimated or redundant wavelet transform is useful for the detection of fine features within the signal. For the case of images, the one-dimensional DWT can be readily extended to two-dimensions. In standard two dimensional wavelet decomposition, the image rows are fully decomposed, with the output being fully decomposed columnwise. In nonstandard wavelet decomposition, all the rows are decomposed by one decomposition level followed by one decomposition level of the columns. Figure 4 illustrates the process of extending 1D-DWT into 2D-DWT.
3.3.2. Discrete Wavelet Packet Transform
Wavelet Packet (WP) is a wavelet transform where the signal is passed through more filters compared to DWT-based approach. Applying DWT or WP on images generates four coefficients; three of them are the detail coefficients, and the remaining one is the average coefficient. It is worth mentioning that the first level of decomposition is the same for both DWT and WP, as illustrated in Figure 5. The differences start being noticed from the second level of decomposition.
(a)
(b)
The differences between DWT and WP can be seen in the detail coefficients where the next decomposition of DWT is applied on the average coefficients from the previous decomposition (Figure 5(a)). The next decomposition of WP is applied on all previous decomposition coefficients (Figure 5(b)). An example of applying DWT and WP on a phantom slice at different levels of decomposition is illustrated in Figure 6 where the average quadrants in both DWT and WP are the same, but the details quadrants are transformed in WP, unlike in DWT [20]. The number of quadrants in DWT is increased linearly by 3 as the decomposition level is increased by 1 (replace the LL-filter output by its four transforms), as illustrated in Details quadrants from previous decomposition levels of WP are transformed; the number of quadrants are increased exponentially by 4 as the level is increased by 1 (replace each filter output by its 4 transforms), as illustrated in
(a) Original image
(b) DWT level1
(c) WP level1
(d) DWT level2
(e) WP level2
(f) DWT level3
(g) WP level3
(h) DWT level4
(i) WP level4
(j) WP full decomposition
3.3.3. Discrete Wavelet Transform (3D-DWT)
Section 3.3.1 has demonstrated that 2D-DWT is a generalization of 1D-DWT applied on all rows and columns using either standard or non-standard decomposition. Applying 3D-DWT is not easy; the difference between 2D images and 3D volumes is the third dimension (depth or Z-axis). The expected transform after applying 3D-DWT is illustrated in Figure 7.
The original volume is transformed into 8 octants (features) in the wavelet domain. Mathematically, 3D-DWT is the process of applying 1D-DWT on each vector in Z-axis which has the same X-axis and Y-axis coordinates after applying 2D-DWT for all comprising frames. Algorithm 4 explains the pseudo code for applying 3D-DWT on 3D data set, and the filter structure of 3D Haar wavelet transform is illustrated in Figure 8.
3.4. Ridgelet Transform
Recently, ridgelet transform [21–23] has been generating a lot of interest due to its superior performance over wavelets. While wavelets have been very successful in applications such as denoising and compact approximations of images containing zero dimensional (point singularities), they do not isolate the smoothness along edges that occurs in images [24]. Wavelets are thus more appropriate for the reconstruction of sharp point-like singularities than lines or edges. These shortcomings of wavelets are well addressed by the ridgelet transform, as they extend the functionality of wavelets to higher dimensional singularities, and are effective tools to perform sparse directional analysis [25]. The basic building block of these transforms is the finite radon transform (FRAT), and HWT has been used to perform FRIT. Applying FRAT on image can be presented as a set of projections of the image taken at different angles to map the image space to projection space. Its computation is important in image processing and computer vision for problems such as pattern recognition and the reconstruction of medical images.
For discrete image data, a projection is computed by summation of all data points that lie within specified unit-width strips; those lines are defined in a finite geometry [26]. It can be obtained by applying 1D Inverse Fast Fourier Transform (1D-IFFT) on the 2D Fast Fourier Transform (2D-FFT) restricted to radial lines going through the origin.
FRAT of a real function on the finite grid is defined in (5), [16]: Here, denotes the set of points that make up a line on the lattice as in (6): To compute the th radon projection (i.e., the th row in the array), all pixels of the original image need to be passed once and use histogrammers: one for every pixel in the row [16]. At the end, all histogrammed values are divided by to get the average values.
Once the wavelet and radon transforms have been implemented, the ridgelet transform is straightforward. Each output of the radon projection is simply passed through the wavelet transform before it reaches the output multiplier.
As shown in Figure 9, ridgelets use FRAT as a basic building block. FRAT maps a line singularity into point singularity, and the wavelet transform has used to effectively handle the point singularity for the radon domain.
Analysing an object with curve singularity implies that ridgelet coefficient will be not sparse and object with curved singularities is still curved or linear after the radon transform where the wavelet transform cannot detect it properly because it is still not a point singularity [28]. Figure 10 shows a real chest slice from a CT scanner [29] in ridgelet domain at different block sizes.
(a)
(b)
(c)
(d)
(e)
(f)
3.5. Statistical Modeling Using Hidden Markov Models
For block-based segmentation using statistical classification, an image is divided into blocks, and a feature vector is formed for each block by grouping statistics of its pixel intensities [3, 30]. Conventional block-based segmentation algorithms classify each block separately, assuming the independence of feature vectors.
Segmentation in MRFM is achieved by maximizing a posteriori probability of the segmentation depending on a given image data [31]. MRFM is an earlier version of the HMM, where states in MRFM are directly visible to the observer, and then the state transition probabilities are the only parameters [32]. Markov chain is an edge labeled directed graph (Figure 11), where each node represents a state, and the edge-label have probabilities of moving the state to the end of the directed edge [9, 33].
HMMs represent a widespread approach to the modeling of sequences as they attempt to capture the underlying structure of a set symbol strings. The use of HMM for shape recognition has not been widely addressed. Only a few works have been found to have some similarities with the proposed approach. In the first, He and Kundu [34] utilized HMMs to model shape contours through autoregressive (AR) coefficients. The use of circular HMM for shape recognition improving scaling and deformation robustness is proposed at [35],[36].
HMM is basically a stochastic finite state automaton, formally defined by the following elements [37]: a set of states ; a state transition probability distribution matrix representing the probability to go from state to ; a set of observation symbols where used to be a -dimensional vector (in the case of discrete HMM); an observation symbol probability distribution or emission matrix indicating the probability of emission of symbol when system state is ; an initial state probability distribution representing probabilities of initial states. For convenience, HMM has been denoted as a triplet which uniquely determines the model.
When modeling a sequence of observation symbols it is usual to use a so-called “left to right HMM” [37], which has only partial state transition matrix, such that +step, where step is a constant usually equal to 1 or 2.
3.5.1. -HMM
The problem with 2D-HMM is the double dependency of on its two neighbors, and , as illustrated in Figure 12. This does not allow the factorization of computation as in 1D, where must only depend on one neighbor at a time. However, this neighbor may be the horizontal () or the vertical () [9, 38].
Each slice is a two-dimensional matrix which can be classified by an optimum set of states with maximum probability; these states are mapped into classes or segmented objects. The basic assumption of applying HMM on medical images is to use the embedded-HMM by defining a set of Markovian superstates, within each superstate there is a set of simple Markovian states. The superstate is first chosen using a first-order Markov state transition probability based on the previous superstate. A simple Markov chain is then used to generate observations in this superstate. Thus, superstates are related to rows (or any equal size blocks), and simple states are related to columns (or smaller blocks comprised the superstate).
3.5.2. Defining the Initial States
To generate an observation sequence using HMM, an initial state must be chosen according to the initial state distribution; then an observation sequence should be chosen according to the probability distribution in the initial state [9]. Many feature selection techniques have been tested in medical images, and the best results were achieved using the grey-scale values of the medical image pixels to generate the probability distribution matrix. After defining the initial states, transition to a new state is taking place according to the state transition probability matrix for the current state and depending on the sequence of observations.
3.5.3. Training HMM
HMM considers observations statistically dependent on neighboring observations through transition probabilities organized in a Markov mesh. Training HMM for images is achieved by dividing the image into nonoverlapping, equally sized blocks, from each of which a feature vector is extracted. Each block and its feature vector evaluate the observation which has its own transition probability matrix. Training HMM produces an estimated state transition probability matrix and estimated emission probability matrix. After building the observation sequence, the model parameters are estimated based on the blocks' statistics. These classes or states were determined using Viterbi algorithm, which depends on the sequence of observations; estimated state transition probability matrix; emission probability matrix [30, 40]. Finally, pixels which belong to the same class are grouped together to evaluate a segmented image.
3.5.4. Testing HMM
The feature vectors for a testing image are generated to find the set of classes with maximum posteriori according to the trained HMM. The feature vector for each block may be changed at every single state. Once the block state is known, the feature vector will be independent of the other blocks; any two blocks may be more likely to be in the same state if they have close intensities [8]. In other words, testing HMM will generate the state transition probability matrix and emission probability matrix for a given comparable data [8, 9].
4. Results and Analysis
The proposed approach has been tested on NEMA IEC body phantom [41] (DATA SET 1) and real chest images from a CT scanner [29] (DATA SET 2). DATA SET 1 consists of an elliptical water filled cavity with six spherical inserts suspended by plastic rods of inner diameters: 10, 13, 17, 22, 28, and 37 mm [41]; this phantom is illustrated in Figure 13(a). Figure 13(b) is a slice example of the CT outputs which stacked with all other slices to form the 3D volume which is illustrated in Figure 13(c). DATA SET 2 is a CT scanner output of the chest area in a human body [29]. This DATA SET includes blood, bones, muscles, tumor and, other organs in the body.
(a)
(b)
(c)
4.1. Segmentation Performance Metrics
Many techniques can be used for segmentation, and each technique has a different segmentation performance and quality. Listed below are some performance measurement methods that have been used in this paper to test and compare the segmentation techniques.
(i) Dice Similarity Coefficients
Dice Similarity Coefficients (DSCs) are a statistical validation metric used to evaluate the performance of both the reproducibility of manual segmentations and the spatial overlap accuracy of automated probabilistic fractional segmentation. The DSC value is a simple and useful summary measure of spatial overlap, which can be applied to studies of reproducibility and accuracy in image segmentation.
The value of a DSC ranges from 0 indicating no spatial overlap between two sets of binary segmentation results to 1 indicating complete overlap [42]. As illustrated in Figure 14, DSC measures the spatial overlap between two samples, A and B target regions, and is defined as where is the intersection.
(ii) Euclidean Distance
Euclidean Distance (ED) is the straight line distance between two points. It can be used with DATA SET 1 to compare the measured diameters with the original diameters provided with the phantom description [41].
(iii) Processing Time
The future work of this paper is to segment the medical images automatically in real time and get the results while the patient is waiting. the processing time of the segmentation methods are different, and the processing time may be used as a comparison factor for these methods.
Signal to Noise Ratio (SNR) is also used to differentiate between wavelet and ridgelet output quality. SNR is used in image processing as a physical measure of the sensitivity of an imaging system. Industry standards measure SNR in decibels (dB) of power, and therefore apply the rule to the peak SNR ratio. Industry standards measure and define sensitivity in terms of the ISO film speed equivalent; image quality is excellent when SNR is 32 dB and if SNR is 20 dB then, it means acceptable image quality.
4.2. Wavelet Versus Ridgelet
Parameters in the wavelet transform are points in the cartesian grid. As illustrated in the first part of Figure 15, each point performs a pixel in the image or an entry in a 2D matrix. However, in ridgelet transform, straight lines evaluate the image in the frequency domain rather than those points in wavelet domain. Parameters in ridgelet domain are in a polar domain (second part of Figure 15) where is the intercept, and is the angle [24, 25].
(a)
(b)
Table 1 illustrates the SNR values of extracted features from DATA SET 2 in spatial domain, wavelet domain at different levels of decomposition, and in ridgelet domain at different block sizes.
It can be seen from Table 1 that small values of SNR have been obtained for all techniques; this is due to the noise from the acquisition systems which will be a part of the medical image itself after the reconstruction of all the slices. Relatively, better SNR values can be achieved with the second level of wavelet decomposition and as the block size is getting bigger with the ridgelet transform where the transformed image is getting more similar to the original image.
4.3. DWT
Applying segmentation techniques on 2D slices requires more time compared to the 3D volumes-based approaches time. The time required to look for the best slice that includes the spheres in full diameters is not required in 3D volume segmentation processes. Segmented volume for DATA SET 1 using 3D-thresholding is illustrated in Figure 16(a), 3D-wavelet technique is applied on the same phantom data (DATA SET 1), and the selected spheres are detected as illustrated in Figure 16(b). Another example of applying 3D-DWT on real CT images for DATA SET 2 is illustrated in Figure 16(c). Table 2 compares the errors in spheres diameters using 3D segmentation techniques with existing measurements using 2D segmentation techniques for slice number 19. ED has been used to measure the spheres diameters and calculate the error percentages for each technique sphere diameter error percentages have been calculated according to the following:
(a) Segmented volume using 3D-thresholding
(b) 3D-DWT for DATA SET 1
(c) 3D-DWT for DATA SET 2
From Table 2, in the case of -means clustering, tumor volumes are underestimated by approximately 5-6% in most cases, however, for the two smaller spherical inserts, with diameter of 10 mm and 13 mm, respectively, these underestimations are significantly greater. For the smallest sphere, more than a 13% volume discrepancy is recorded, with the -means algorithm finding it difficult to quantify the tumor accurately. Sphere 2 similarly is massively underestimated (11 : 5%). Unlike -means clustering, MRFM tends to overestimate the volumes of the spherical inserts, with the exception of Spheres 1 and 2.
Outer diameters have been measured in the case of 3D segmentation, and the inner diameters errors have been calculated based on the thickness of spheres edges. All spheres diameters detected using 3D-thresholding and the errors were over estimated by (0.5–2.5%), they were increasing while the sphere diameter increasing.
Spheres diameters are reduced to the half with each decomposition level of wavelet transform. Three decomposition levels of DWT have been applied on NEMA phantom using two different filters (Haar, Daubechies), and the measured diameters were doubled at each level to produce a fair comparison with the other available techniques. It can be seen that most of the error percentages were decreasing while the spheres diameter increasing, it is worth mentioning that there is no upper bound of the spheres diameters to keep the errors decreasing because the ROI becomes clearer and easier to be detected and measured properly. But tumors in real life are usually very small in the early stage cancer, and the problem is to detect those turnouts as soon as possible.
By applying one decomposition level of 3D-DWT on spatial domain and using the LLL filter output, underestimated percentages have been achieved for the three small spheres (10, 13, and 17 mm) and overestimated percentages for the three big spheres (22, 28, and 37 mm). DWT proved efficient in detecting the big obstacles where the biggest sphere (37 mm) was detected with a very small error percentage (0.8%).
The two smallest spherical inserts are still underestimated in all techniques except the 3D-thresholding. The large volumetric errors encountered using this acquisition exist as a consequence of the poor slice thickness setting selected for the scan. The 4.25 mm slice thickness causes large fluctuations in transaxial tumor areas to occur between image slices. This problematic characteristic occurs most notably with the smallest spherical inserts, where single voxel reallocation causes a large deviation in percentage error. In Figure 16, the percentage error computed between the actual sphere volume and the volumes obtained using all methodologies for each of the six tumor inserts is plotted. It can be seen that all techniques are settled down according to the error percentages as the sphere diameters increased (Figure 17).
4.4. HMM Experimental Results
HMMs have been used for segmentation, which can be applied either in the spatial or multiresolution domain using wavelet or ridgelet transforms. The weakness of HMMs is its long processing time, compared with the other evaluated techniques.
DATA SET 1 has been utilized to perform the experimental study using HMMs and MRA for medical image segmentation, and the achieved results are illustrated in Table 3. Figures 18 and 19 illustrate the outputs of applying HMM on a slice from DATA SET 1 and real brain slice, respectively.
(a)
(b)
(a)
(b)
From Table 3, applying HMMs on NEMA phantom slice in ridgelet domain failed to segment ROI where most of the spheres diameters could not be detected. As well as the long computational time problem, it is five times more than the required time for applying HMMs on the slice in the wavelet domain. Transforming any image into ridgelet domain changes the comprising pixel specifications as well as the dimensions of the image. The padding row is added to each block leading to long computational time for HMM segmentation. Relatively, the most accurate results have been achieved with applying HMMs directly on the spatial domain without any transformation. But in wavelet domain, applying HMMs requires less time which make it very useful in real time segmentation systems.
It can be seen from Figure 18 that different obstacles were merged such as bones and brain; the same thing happened in Figure 19 where the liver is merged with the skin. This merge is due to the number of HMM classes (states) where HMMs decided that brain and bones in Figure 18 will be in the same class because just three HMM classes were used. But each obstacle will be in separate class if the number of classes increased to the number of all obstacles in the image, the same explanation for Figure 19. There are some available techniques that can be used to define the best number of classes in HMM such as BiC.
HMMs have been also applied on DATA SET 2 as illustrated in Figure 20, and the DSC values for the segmented slices have been measured and illustrated in Table 4 which evaluates the performance based on a manual segmented image where better results have been achieved using HMMs in the spatial domain.
(a)
(b)
(c)
Many techniques have been previously implemented for medical volumes segmentation; some of them were illustrated in this paper and compared with the proposed techniques. Promising results have been achieved using 3D segmentation techniques directly to medical volumes and the statistical models using HMMs; both techniques have the same problem which is the computational time. Many acceleration processing methods have been recently implemented such as FPGAs, Matlab accelerators, Feature Reduction (FR) techniques, GPUs and many more. HMM can be applied with different MRA transforms such as wavelet, ridgelet, and curvelet to achieve promising results. It is worth mentioning that based on the experiments carried out on these specific medical data in this paper, HMMs can be classified as one of the ideal medical volume segmentation techniques compared to the other proposed techniques.
5. Conclusions
Segmentation is very important for medical image processing, to detect tumors and regions of interest for radiotherapy planning and cancer diagnosis. A novel and sophisticated segmentation system was developed specifically for 3D data segmentation. The developed techniques within the system have been tested on phantom data [41]. The system used to quantify tumors within the data of predefined volumes, and these results were compared with those obtained from 2D approaches such as thresholding. Thresholding addressed many problems in multiresolution analysis including de-noising and edge detection, but data loss was the main problem caused. Statistical models using hidden Markov models have been also investigated for segmentation, which can be deployed in the spatial domain or multiresolution domain. The weakness of HMMs is its long computation time required for the calculation of their models which was much smaller compared with the other evaluated techniques.
The proposed system was also tested on other data set for real human chest images from a CT scanner and has shown promising results. Ongoing research is focusing on the implementation of other 3D novel feature extraction techniques for medical images based on 3D-ridgelet and 3D-curvelet transforms. In order to speed up the proposed techniques; a graphical processing unit (GPU) will be deployed with more focus on the implementation of higher dimensional HMMs for a more accurate and automatic volume segmentation system.