Molecular Image Segmentation Based on Improved Fuzzy Clustering

Segmentation of molecular images is a difficult task due to the low signal-to-noise ratio of images. A novel two-dimensional fuzzy C-means (2DFCM) algorithm is proposed for the molecular image segmentation. The 2DFCM algorithm is composed of three stages. The first stage is the noise suppression by utilizing a method combining a Gaussian noise filter and anisotropic diffusion techniques. The second stage is the texture energy characterization using a Gabor wavelet method. The third stage is introducing spatial constraints provided by the denoising data and the textural information into the two-dimensional fuzzy clustering. The incorporation of intensity and textural information allows the 2DFCM algorithm to produce satisfactory segmentation results for images corrupted by noise (outliers) and intensity variations. The 2DFCM can achieve 0.96 ± 0.03 segmentation accuracy for synthetic images under different imaging conditions. Experimental results on a real molecular image also show the effectiveness of the proposed algorithm.


INTRODUCTION
Molecular imaging techniques such as positron emission imaging, fluorescent imaging, and isotope radiation imaging have undergone explosive growth over the past few decades. It will allow clinicians not only to measure concentrations of interesting molecules quantitatively, but also to visualize the interactions of molecular markers in vivo, thus extending the emphasis of radiological imaging beyond the anatomical and functional levels [1]. Integrations of molecular information specific to each patient with anatomical information obtained by conventional imaging methods such as the magnetic resonance imaging (MRI), X-ray computed tomography (CT), and ultrasound (US) will undoubtedly enhance the ability to fight disease. Image segmentation is a preliminary and crucial step for subsequent image applications such as quantification of molecular concentration, image registration, and integration. However, molecular images often suffer from a low signal-to-noise ratio (SNR); this will lead to difficulties with its segmentation.
The fuzzy clustering algorithm, more widely used as fuzzy C-means algorithm (FCM) [2], has been successfully utilized in medical image segmentation [3][4][5][6]. The most important feature of the FCM is that it allows each pixel to belong to multiple clusters according to its degree of member-ship in each cluster, which makes the clustering methods able to retain more information from the original image as compared to the case of hard segmentation. FCM works well on images with low levels of noise, but there are two disadvantages of the FCM used in segmentation of noise-corrupted images. One is that the FCM does not incorporate the information about the spatial context, which makes it sensitive to the noise and other imaging artifacts. The other is that the cluster assignment is based solely on the distribution of the pixel intensity, which makes it sensitive to intensity variations due to the illumination or the object geometry [7]. In order to improve the robustness of conventional FCM, many algorithms have been presented in the literatures. These methods can be divided into two main groups: imposing spatial constraints to clustering algorithms [3,[5][6][7] and introducing other features or dissimilarity index that is insensitive to intensity variations in the objective function of FCM [5,7]. This paper presents a novel algorithm based on fuzzy logic for molecular image segmentation. In this algorithm, two factors to improve the robustness of conventional FCM are considered. Due to the low SNR of molecular images, image denoising is taken for a prelude to the segmentation. A denoising method which combines a Gaussian noise filter with an anisotropic diffusion (AD) technique is presented to alleviate noise in molecular images. Since the Gabor wavelet 2 International Journal of Biomedical Imaging representation of molecular images is relatively robust to intensity variations, a texture characterization method derived from Gabor filters bank is presented to extract texture information from images. Spatial constraints provided by the denoising data and texture information provided by the Gabor wavelet are embedded in the objective function of a novel two-dimensional fuzzy clustering (2DFCM) algorithm.
The remainder of this paper is organized as follows. Section 2 proposes the denoising method. Section 3 introduces the multichannel Gabor filters and the texture feature characterization. In Section 4, we present in detail the new two-dimensional FCM algorithm (2DFCM) which integrates both intensity information and texture information. The experimental comparisons are presented in Section 5. Section 6 concludes the paper.

MOLECULAR IMAGE DENOISING
Gaussian noise is the most common noise broadly existed in signal processing sciences. Ling and Bovik [8] proposed a method to smooth molecular images by assuming that the noise follows an additive Gaussian model. Following Ling and Bovik's notion, we also assume that molecular images are corrupted by a zero-mean Gaussian white noise.
The FIR filter is well known for its ability to remove Gaussian noise from signals but it does not work very well in the image processing since it blurs edges within the image. The Gaussian noise filter (GNF) [9], combining a nonlinear algorithm and a technique for automatic parameter tuning, is a valid method for estimation and filtering of Gaussian noise. The GNF used in this paper can be summarized as follows. Let X = {x 1 , x 2 , . . . , x n } be a set of n data points in the noisy image. The output Y = {y 1 , y 2 , . . . , y n } is defined as where N i stands for the neighborhood configuration with respect to a center pixel x i , and N R is the cardinality of N i . The automatic tuning of the parameter p is a key step in GNF. Let MSE(k) denote the mean square error between the noisy image filtered with p = k and the same image filtered with p = k − 1. A heuristic estimate of the optimal parameter value is where MSE k m = MAX MSE(k) .
The GNF can remove intensity spikes due to the Gaussian noise. However, it has limited effect on suppressing little intensity variations caused by the neighboring smoothing.
Since the conventional FCM is a method based on the statistical feature of the image intensity, a piecewise-smooth intensity distribution will be greatly beneficial to it. We pursue a more desirable denoising result by following the GNF with an anisotropic diffusion filter. Yu and Acton [10] provided an improved anisotropic diffusion filter called speckle reducing anisotropic diffusion (SRAD) which outperforms the traditional Perona-Malik nonlinear diffusion [11]. Although SRAD is proposed for the speckle reduction in synthetic aperture radar (SAR) or ultrasound images, its advantages in mean preservation, variance reduction, and edge localization are also preferable for molecular images. The SRAD used in this paper can be formulated as a diffusive process: where c(q) represents the diffusion coefficient, q(x, y; t) is the instantaneous coefficient of variation served as the edge detector in the noise image. q(x, y; t) combines a normalized gradient magnitude operator and a normalized Laplacian operator: where ∇ is the gradient operator and Y is the image filtered by GNF. q 0 (t) is the scale function serving as the diffusion threshold which can be approximated by using a heuristic constant q 0 with the exponential decay function Here ρ is a constant typically set to 1/6. Suppose that the output of the SARD with Y = {y 1 , y 2 , . . . , y n } as the input can be represented by X * = {x 1 * , x 2 * , . . . , x n * }. To clearly illustrate the denoising effect of GNF plus SRAD, Figure 1 shows a group of filtering results of GNF alone, SRAD alone, GNF plus SRAD, and the anisotropic median-diffusion (AMD) [8] on a synthetic molecular image. From the filtering results comparison, it is seen that the denoising method of integrating GNF with SRAD can overcome the intensity fluctuation effect of GNF and the "blocky" effect of SRAD and MAD.

TEXTURE CHARACTERIZATION
A molecular image illustrates the distribution of a certain molecule [8]. Since the photon has different transportation characteristics in different turbid tissues, a molecular image can be divided into several separate regions with each region showing similar intensity (implying similar molecular concentration) and certain kind of textural pattern. Because the photon distribution in a turbid tissue is not usually uniform, the intensity within a region usually changes gradually. This intensity variation can cause errors when attempting to segment images using intensity-based classification methods. Intuitively, if a feature insensitive to the slowly varying intensity can be introduced into the classification, the performance of the image segmentation could be improved. Here, a texture characterization method based on the Gabor wavelet is utilized to obtain this desirable feature. A large number of texture classification techniques have been proposed in the past two decades [12]. Gabor wavelet has been a popular method because it can capture the local structure corresponding to spatial frequency, spatial localization, and orientation selectivity. As a result, Gabor wavelet representation of an image should be robust to intensity variations [13,14]. A Gabor function in the spatial domain is a sinusoidal modulated Gaussian. The real impulse response of Gabor filter is given by where x = x cosθ + y sinθ, y = −x sinθ + y cosθ, (x, y) represent rotated spatial-domain rectilinear coordinates, u is the frequency of the sinusoidal wave along the direction θ from the x-axis, σ x and σ y define the size of the Gaussian envelope along xand y-axes, respectively, which determine the bandwidth of the Gabor filter. The frequency response of the filter is given by where σ u = 1/2πσ x , σ v = 1/2πσ y . By tuning u and θ, multiple filters that cover the spatial frequency domain can be obtained. In our study, Gabor wavelets with four different scales, μ ∈ {π/4 √ 2, π/4, π/2 √ 2, π/2}, and eight orientations, θ ∈ {0π/8, 1π/8, . . . , 7π/8}, are used. Let X(x, y) be the intensity level of an image. The Gabor wavelet representation is the convolution of X(x, y) with a family of Gabor kernels: where * denotes the convolution operator, and G u,θ is the convolution result corresponding to the Gabor kernel at the scale μ and the orientation θ. The next step is to compute the textural energy in G u,θ . The textural energy is a measure widely used to characterize the image texture. The en-ergy that corresponds to a square window of the image G u,θ centered at x and y is defined as where M 2 is the total number of pixels in the window, and F(.) is a nonlinear, sigmoid function of the form where α equals 0.25. The texture feature image is finally given by As an example, Figure 2(a) shows a synthetic image with the intensity inhomogeneity. Figure 2(b) gives the texture energy bank (E μ,θ ) illustration. Figure 2(c) shows the texture feature image. From this example, it is seen that the texture feature characterization using Gabor wavelet is insensitive to the intensity inhomogeneity.

FCM
. , x n } be a set of n data points, and let c be the total number of clusters. The objective function of the FCM [2] for partitioning X into c clusters is given by where m j , j = 1, 2, . . . , c, represent the cluster prototypes and μ i j gives the membership of pixel x i in the jth cluster m j . The parameter b is the fuzzy index that satisfies b ∈ (1, ∞) and controls the degree of "fuzziness" in the resulting classification. The fuzzy partition matrix satisfies Under the constraints condition of (15), taking the first derivations of (14) with respect to μ i j and m j and setting those equations to zero yield necessary conditions for (14) to be minimized. Performing iteration through these two necessary conditions leads to an iterative scheme for minimizing the objective function. The objective function (14) is minimized when high membership values are assigned to pixels whose intensities are close to the centroid of its particular class, and low membership values are assigned when the pixel data is far from the centroid [2]. After FCM clustering, a segmentation of the image can be obtained by assigning each pixel solely to the class that has the highest membership value for that pixel.
Although the membership allows a pixel to deviate from multiple cluster prototypes, the spatial correlation between adjacent pixels is not considered.

FCM with spatial constraints
A popular method to introduce the local spatial context into the pixel classification is the spatial constraint. The spatial constraint is to let the spatial information influence the classification of the pixel of interest [5,6]. Let N i denote the configuration of neighbors that exists in a window around x i . According to the assumption that real-world images usually have strong correlation among neighboring pixels, if the pixel x i belongs to the cluster with the prototype m j , then pixels in N i and the center pixel x i should have similar and high membership values in m j . This original idea of incorporating local spatial constraints in the FCM is formulized as [15] where N i stands for the neighborhood configuration with respect to a center pixel x i , N R is the cardinality of N i , α controls the effect of the neighboring penalty. The second term on the right side of (16) allows the labeling of a pixel to be influenced by the labels in its immediate eight neighborhoods and aims at keeping continuity in the neighboring window. The problem with (16) is that computing the neighborhood terms will cost much more time than clustering. In order to reduce the complexity of computing the neighborhood terms, the dissimilarity measurements between the whole neighborhood configuration and the prototype m j can be replaced by a distance from a feature data of N i to m j . The feature data of the neighborhood configuration can be obtained by several kinds of neighboring window filters, such as the linear filter or the median filter. This approach is expressed in the following objective function [5]: where x ∧ i is a mean or median of neighboring pixels lying within a window around x i . Here, we modify (17) by substituting x ∧ i with denoising molecular image data x * i . The objective function for the FCM with spatial constraints (called FCM S later) is given by

2DFCM
Equation (18) introduces spatial constraints into the clustering procedure. However, the classification result of (18) still solely depends on the intensity distribution of the image, which makes it sensitive to intensity variations within a turbid tissue. With the texture information obtained by the Gabor wavelet bank, the two-dimensional fuzzy C-Means (2DFCM) algorithm is constructed by integrating both the intensity and the texture information. Suppose that the texture feature image is T = {t 1 , t 2 , . . . , t n }, the objective function of 2DFCM can be expressed as The influence of the texture characterization imposed on the clustering procedure can be controlled by a constant vector β i (i = 1, . . . , n); the prototype of texture image data is represented by v j ( j = 1, . . . , c). The choice of β i is based on the following principle. If t i is large, implying the texture energy is dominant, and β i should be large; if t i is small, implying the texture energy is weak, and β i should be also small. The β i is determined by β i = (Bt i )/max(T), where B is a constant and its optimized value is determined by "trial-and-error" technique (see Section 5 for details). The optimization problem under the constraint of U as stated in (15) can be solved using one Lagrange multiplier: Taking the derivative of F with respect to μ i j and setting the result to zero, we can obtain an equation for μ i j with unknown Utilizing the constraint of U can be solved as Substituting (22) into (21), a necessary condition for (19) to be at a local minimum will be obtained: Similarly, zeroing the derivative of F with respect to m j and v j , we have

Implementation of 2DFCM
For the 2DFCM, the number of prototypes (c) and the initial centroids (M = {(m j , v j ) | j = 1, . . . , c}) ought to be known at the beginning of iterative procedures. A maximum likelihood approach by processing and analyzing the twodimensional (2D) histogram of X and T is used to estimate c and M. The number of prototypes (c) and initial prototypes (M ) is estimated by following steps.
(1) Count the number of peaks in the 2D histogram and record it as PeakNum prev . (2) Filter the histogram using a five-by-five Gaussian filter with zero mean and a standard deviation of 0.6. The procedure of 2DFCM can be summarized in the following steps.
(1) Filter the image using GNF followed by SRAD to generate the denoising data X * . (2) Filter the image using Gabor wavelet band and compute the texture feature image T. (3) Formulate the 2D histogram using the denoising data X * and the texture feature image T. Estimate the number of clusters (c) and initial prototypes (M). (4) Repeat the following steps until the centroids variation is less than 0.001.
(a) Update the membership function matrix using (23). (b) Update the centroids using (24). (c) Calculate the centroids variation between before updating and after updating.

EXPERIMENTAL RESULTS AND DISCUSSIONS
We perform experiments on a PC with 2.0 GHz Pentium processor using Visual C++ 6.0. To illustrate the performance of the 2DFCM, we first test it using simulated molecular images from which the ground truth data is available. Simulated molecular images are obtained by using MOSE (Monte Carlo optical simulation environment) [16][17][18] developed by Bioluminescence Tomography Lab, Department of Radiology and Department of Biomedical Engineering, University of Iowa (http://radiology.uiowa.edu/). MOSE is based on Monte Carlo method to simulate bioluminescent phenomena in the mouse imaging and to predict bioluminescent signals around the mouse. The optimized α and B in the 2DFCM should be obtained by "trial-and-error" technique. We first choose an appropriate value for α based on the segmentation performance of the FCM with spatial constraints (FCM S) (the objective 6 International Journal of Biomedical Imaging function is formularized as (18). We take a group of values for α ranging from 0.25 to 6 to test the misclassification rate. With the increasing of α, the number of misclassified pixels reduces. However, after α exceeds 3, the segmentation performance of the FCM S has no apparent changes. Therefore, we set α = 3.5 in our study, which is a value that can produce steady and good results. Then, we choose an appropriate value for B based on the segmentation performance of the 2DFCM. We also take a group of values for B ranging from 5 to 80 to test the misclassification rate. After B exceeds 36, the segmentation performance of the 2DFCM has no apparent changes. Therefore, we set B = 36 in our work, which gives steady and satisfactory results. The computation time of the proposed algorithm on an image of 128 × 128 is approximately 12 seconds. About two thirds of total time are consumed in texture characterization based on Gabor wavelet.
The first example is applying algorithms to a synthetic cellular image and comparing the 2DFCM with other three algorithms, including the FCM on the original image, the FCM with spatial constraints, and the FCM on the texture feature image. The model to generate synthetic molecular images is illustrated in Figure 3(a). The simulated molecular image (128 × 128) corresponding to Figure 3(a) is shown in Figure 3(b). Then Figure 3(b) is corrupted by the intensity inhomogeneity (as shown in Figure 3(c)) to generate the final synthetic image (as shown in Figure 3(d)). Figure 3(e) shows the image filtered by the GNF plus SRAD. Figure 3(f) shows the texture feature image obtained by the Gabor wavelet bank. Figures 3(g)-3(j) give the segmentation results of the FCM on the original image (Figure 3(d)), the FCM with spatial constraints, the FCM on the texture feature image (Figure 3(f)), and the 2DFCM, respectively. We quantify the algorithm performance in terms of three parameters defined as follows: SA represents the total segmentation accuracy; US is the under segmentation rate; OS denotes the over segmentation rate. N CORRECT is the number of correctly classified pixels; N TOTAL is the total number of pixels; N f p is the number of pixels that do not belong to the class and are segmented into this class; N f n is the number of pixels that belong to the class and are not segmented into the class; N p is the number of all pixels that belong to the class; N n is the number of all pixels that do not belong to the class. There are totally four algorithms that are compared in our experiments. Table 1 gives the SA, US, and OS comparisons among the four algorithms, correspondingly.
To further test the segmentation performance of the proposed method, a group of synthetic images under different imaging conditions are utilized. Nine synthetic images are shown in Figure 4. These images are organized into the form with different photons density along the vertical direction and different types of intensity inhomogeneity along the horizontal direction. Table 2 summarizes the segmentation accuracy of the FCM on the original image, the FCM with  spatial constraints, the FCM on the texture feature image, and the 2DFCM, respectively. The second example is applying the algorithms to a realmolecular image (256 × 256) (as shown in Figure 5 Figure 5(a) that the middle of the tissue appears homogeneously bright. However, the molecular concentration decreases in the boundary area, which leads to the intensity variation near the boundary. The conventional FCM on the original image and the FCM with spatial constraints produce undersegmentation results and the FCM on the texture feature image shows oversegmentation.
From the experimental results, we can see that the denoising effects of the GNF plus SRAD are satisfactory. The Gabor wavelet bank can represent the texture information in the molecular image without being disturbed by the intensity variation. The FCM produces the worst result due to the fact that no spatial constraints are used in it. The FCM with spatial constraints produces more smoothed segmentation results than the FCM. However the intensity inhomogeneity makes the segmentation result degenerate. Since the 2DFCM utilizes both the intensity and texture information simultaneously, it produces more satisfactory results than other methods.

CONCLUSIONS
In this paper, we have developed a novel algorithm based on the fuzzy clustering for the molecular image segmentation. Considering that there are two disadvantages for the conventional FCM in the image segmentation, its successful employment in the molecular image segmentation requires overcoming nonrobust factors by introducing spatial constraints and the texture feature of images into the clustering. To alleviate noises in molecular images, a denoising method combining GNF plus SRAD is proposed. We use the denoising data obtained by GNF plus SRAD to compose spatial constraints for the new 2DFCM. By utilizing the Gabor wavelet representation and the texture energy characterization, the texture feature that is insensitive to intensity variations is introduced into the 2DFCM. Quantitative evaluation demonstrates the superiority of the 2DFCM over the conventional FCM in the molecular image segmentation.