Neighbourhood search feature selection method for content-based mammogram retrieval

Abstract

Content-based image retrieval plays an increasing role in the clinical process for supporting diagnosis. This paper proposes a neighbourhood search method to select the near-optimal feature subsets for the retrieval of mammograms from the Mammographic Image Analysis Society (MIAS) database. The features based on grey level cooccurrence matrix, Daubechies-4 wavelet, Gabor, Cohen–Daubechies–Feauveau 9/7 wavelet and Zernike moments are extracted from mammograms available in the MIAS database to form the combined or fused feature set for testing various feature selection methods. The performance of feature selection methods is evaluated using precision, storage requirement and retrieval time measures. Using the proposed method, a significant improvement is achieved in mean precision rate and feature dimension. The results show that the proposed method outperforms the state-of-the-art feature selection methods.

Appendices

Appendix 1

1.1 GLCM features

The GLCM of an $N \times N$ pixel image I contains the probabilities $p_{d,\theta }(i,j)$ of the change over from a grey level i to a grey level j in a given direction $\theta$ separated by a pixel distance d which are expressed as:

$$\begin{aligned} p_{d, \theta }(i, j) = \frac{C_{d, \theta }(i, j)}{\sum _{i=1}^{N_g}\sum _{j=1}^{N_g}C_{d, \theta }(i, j)}, \end{aligned}$$

(3)

where

$$\begin{aligned} C_{d, \theta }(i,j)&= {} \# \{(x_1, y_1), (x_2, y_2) \in N \times N : f(x_1, y_1)=i, f(x_2, y_2)=j,\\&|(x_1, y_1) - (x_2, y_2)| = d, ((x_1, y_1), (x_2, y_2)) = \theta \}. \end{aligned}$$

Here, $\#$ denotes the cardinality of the set, $f(\cdot )$ denotes the grey level of the pixels at locations $(x_1,y_1)$ and $(x_2, y_2)$, and $N_g$ is the total number of grey levels in the image. The common choices for $\theta$ are $0^{\circ }, 45^{\circ }$ , $90^{\circ }$ and $135^{\circ }$ [22]. The expressions and the significance of fourteen statistical measures, called Haralick’s features [26], that is extracted from the GLCMs are given in [16].

1.2 Gabor features

A set of 24 Gabor filters [25] are obtained by appropriate dilations and rotations of the Gabor function,

$$\begin{aligned} g(x, y) = \frac{1}{2\pi \sigma _x \sigma _y} \exp \left[ - \frac{1}{2} \left( \frac{x^2}{\sigma _x^2}+\frac{y^2}{\sigma _y^2}\right) +2\pi jU_h x \right] \end{aligned}$$

(4)

through the generating function

$$\begin{aligned} g_{mn}(x,y) = a^{-m}g(x', y'). \end{aligned}$$

(5)

Here,

$$\begin{aligned} \sigma _x&= {} \frac{(a+1)(2\ln 2)^{\frac{1}{2}}}{2\pi (a-1) U_h },\\ \sigma _y&= {} \displaystyle \frac{1}{2\pi \tan \left( \frac{\pi }{2K}\right) \left[ U_h - 2 \ln \left( \frac{\sigma _u^2}{U_h}\right) \right] \left[ 2 \ln 2 - \frac{(2 \ln 2)^2\sigma _u^2}{U_h^2}\right] ^{-\frac{1}{2}}}, \end{aligned}$$

$j=\sqrt{-1}, x' = a^{-m}\left( x\cos \phi + y \sin \phi \right)$ and $y' = \left( -x \cos \phi + y \sin \phi \right)$; $\phi$ and a are taken as $\frac{n\pi }{K}$ and $\left( \frac{U_h}{U_l}\right) ^{-\frac{1}{s-1}}$, respectively. The following values are used in our implementation as reported in [25]: m and n are integers selected in the interval [1, S] and [1, K], respectively, $U_h = 0.04, U_l = 0.05$, and $K=6$ and $S=4$ are the total number of orientations and the number of scales, respectively. Six features, [37] are extracted from the outputs of Gabour filters given by

$$\begin{aligned} G_{mn}(x,y)= \text{ Re }(I(x, y) *\bar{g}_{mn}(x, y)), \end{aligned}$$

(6)

where $\text{ Re }(\cdot )$ represents the real part of the argument, $*$ is the convolution operator and $\bar{g}_{mn}(x,y)$ is the complex conjugate of the generating function.

1.3 Daubechies (Db4) wavelet features

For an integer r, Daubechies wavelet can be defined as [14],

$$\begin{aligned} \phi _{i,j,k}(x) = 2^{j/2} \phi _r (2^j x - k),\quad j, k \in Z \end{aligned}$$

(7)

where, j is a scale, k is a translation and r is a filter. The coefficients of the orthonormal decomposition pairs in the case of Daubechies 4 wavelet function can be expressed as follows:

$$\begin{aligned} G(n)= & {} \left[ \frac{1-\sqrt{3}}{4\sqrt{2}}, - \frac{3-\sqrt{3}}{4\sqrt{2}}, \frac{3+\sqrt{3}}{4\sqrt{2}}, - \frac{1+\sqrt{3}}{4\sqrt{2}} \right] \end{aligned}$$

(8)

$$\begin{aligned} H(n)= & {} \left[ \frac{1+\sqrt{3}}{4\sqrt{2}}, \frac{3+\sqrt{3}}{4\sqrt{2}}, \frac{3-\sqrt{3}}{4\sqrt{2}}, \frac{1-\sqrt{3}}{4\sqrt{2}} \right] \end{aligned}$$

(9)

As per the parameter estimation method mentioned in [7], the estimated GGD model parameters of the Daubechies (Db4) wavelet coefficient distribution in a $M \times N$ subband $X = {x_{i,j}, i = 1, \ldots ,M, j = 1,\ldots ,N}$, namely $\widehat{\alpha }$ and $\widehat{\beta }$, are obtained as follows,

$$\begin{aligned} \widehat{\alpha } = \left( \frac{\beta }{MN}\sum _{i=1}^M \sum _{j=1}^N \left| x_{i,j}\right| ^{\beta }\right) ^{1/\beta } \end{aligned}$$

(10)

where $\beta$ is an approximation of $\widehat{\beta }$ determined using the Newton–Raphson iterative procedure [7]. In this work, the GGD model parameters are estimated from the details subbands of each level resulted in 18 feature components forming the feature vector for a given mammogram.

1.4 CDF 9/7 wavelet features

The lifting scheme for the CDF 9/7 wavelet includes two set of predict ($p_1$ and $p_2$) and update ($u_1$ and $u_2$) filters. The values used for a given operation in the CDF 9/7 wavelet implementation are as follows: $\alpha \approx -1.58613432$ and $\delta \approx 0.8829110762$ for predicting; $\beta \approx -0.05298011854$ and $\gamma \approx 0.4435068522$ for updating; $\zeta \approx 1.149604398$ for scaling. The prediction errors, updated odd sets and normalized outputs are computed as given in Eqs. 11, 12 and 13, respectively [4].

$$\begin{aligned} x_e^1(n)&= {} x_e(n) - \alpha (x_o(n+1) + x_o(n))\nonumber \\ x_e^2(n)&= {} x_e^1(n) - \delta (x_o^1(n+1) + x_o^1(n)) \end{aligned}$$

(11)

$$\begin{aligned} x_o^1(n)&= {} x_o(n) + \beta (x_e^1(n) + x_e^1(n-1)) \nonumber \\ x_o^2(n)&= {} x_o^1(n) + \gamma (x_e^2(n) + x_e^2(n-1))\end{aligned}$$

(12)

$$\begin{aligned} x_o^3(n)&= {} \zeta x_o^2(n); \quad x_e^3 = \frac{x_e^2(n)}{\zeta } \end{aligned}$$

(13)

The above procedure is repeated with $x = x_o^3$ for the next level decomposition. As mentioned in [28], a maximum of three levels of decomposition is chosen for this work. Next, the GGD [7] model parameters from the details subbands of each level and 32-bin histogram from the third level approximation coefficients are estimated to generate a feature vector of size 50 for a given mammogram.

1.5 Zernike moments features

The mammogram region of interest is pre-processed [10] and 256 ZMs are extracted without the need of previous image segmentation. For a digital image the respective ZMs of $A_{pq}$ of order p with repetition q are computed as follows [24]:

$$\begin{aligned} A_{pq} = \frac{p+1}{\pi }\sum _{i} f(x_i, y_i)V_{pq}^*(\rho , \theta );\quad x_i^2+y_i^2\le 1 \end{aligned}$$

(14)

where i runs over all the image pixels, p is a non-negative integer, q is an integer subject to the constraint $p - |q| = even, |q| \le p$, ‘$*$’ denotes complex conjugation.

$$\begin{aligned} V_{pq}(\rho , \theta )= {} R_{pq}(\rho )\exp (jq\theta ); \quad \rho = x^2+y^2 ;\quad \theta = tan^{-1}(y/x)\end{aligned}$$

(15)

$$\begin{aligned} R_{pq}(\rho )= & {} \sum _{s=0}^{(p-|q|)/2} \frac{(-1)^s \left[ (p-s)!\right] \rho ^{p-2s}}{s!\left( \frac{p+|q|}{2}-s\right) !\left( \frac{p-|q|}{2}-s\right) !} \end{aligned}$$

(16)

Appendix 2

1.1 Genetic algorithm (GA)

In genetic algorithm, a population is initially generated by randomly creating a group of individuals (or feature subsets). Then, these individuals are evaluated using the fitness function, i.e. the total precision value given in Eq. 2. Next, two individuals are selected based on their fitness. There is a higher chance of being selected if the fitness is higher. These individuals then reproduce through crossover to create offspring, which are then mutated randomly. This process continues until an optimal or near-optimal feature subset is found or a certain number of generations is over. Each chromosome representing a feature subset is encoded as a binary string of size equal to the particular feature set size. For example, the chromosome length is nine and 10 for Db4 and CDF 9/7 wavelet based features, respectively. Each gene or element of the chromosome refers to a particular feature and represented as 1 or 0 when it is present or absent, respectively. In our work, the population size considered is 20. The control parameters are tournament selection, single point crossover with 0.7 (probability of crossover) and 0.001 (probability of mutation) [32]. The stopping condition is either the number of generations (e.g. 100) or the maximum fitness value, which is equal to the total number of queries.

1.2 t test

Statistical multivariate t test select the most discriminant features by assessing the significance of the difference between the means of two sample set A and B. The value of the t test ‘t’ is obtained as follows [29],

$$\begin{aligned} D_a= {} \sum {(A_i - \mu _a)^2} \end{aligned}$$

(17)

$$\begin{aligned} D_b= {} \sum {(B_i - \mu _b)^2} \end{aligned}$$

(18)

$$\begin{aligned} V= {} \frac{D_a+D_b}{(n_a-1)+(n_b-1)} \end{aligned}$$

(19)

$$\begin{aligned} \sigma= {} \sqrt{\left( \frac{V}{n_a}+\frac{V}{n_b}\right) } \end{aligned}$$

(20)

$$\begin{aligned} t&= {} \frac{\mu _a-\mu _b}{\sigma } \end{aligned}$$

(21)

where $A_i$ and $B_i$ are the ith elements of the set A and B, $\mu _a$ and $\mu _b$ are its means, $D_a$ and $D_b$ are the sum of squared deviates of the set A and B, V is the estimated variance of the source population and $\sigma$ is the standard deviation of the sampling distribution of sample mean differences. The degree of freedom (d.f.) is $(n_a - 1) + (n_b - 1)$. In our case study, with 20 normal images (set A) and 19 abnormal images (set B), [34] the d.f. is 37. The t value according to the Table of Critical Values for 37 degrees of freedom is 1.305. In the experiment, for a given feature if the ‘t’ value is greater than 1.305 implies that there is a significant mean difference between the normal and abnormal mammograms.

1.3 StARMiner and FD-ASE

The goal of StARMiner algorithm is to find statistical association rules to select a minimal set of features that preserves the ability of discerning image according to their types. Let $x_j$ be a category of an image and $f_i$, an image feature (attribute). The rules returned by the StARMiner algorithm have the format $x_j \rightarrow f_i$. The threshold values used are as follows: $\Delta \mu _{min}$ is the minimum allowed difference between the average of the feature $f_i$ in images from category $x_j$ and the average of $f_i$ in the remaining data set; $\sigma _{max}$ is the maximum standard deviation of $f_i$ values allowed in a category; $\gamma _{min}$ is the minimum confidence to reject the hypothesis $H_0$. StARMiner mines rules of the form $x_j \rightarrow f_i$, if the conditions given in Eqs. (22–24) are satisfied [39].

$$\begin{aligned} \mu _{f_i}(V)= {} \frac{\sum _{k\in V} f_{ik}}{|V|} \end{aligned}$$

(22)

$$\begin{aligned} \sigma _{f_i}(V)= {} \sqrt{\frac{\sum _{k\in V}\left( f_{i_k} - \mu _{f_i}(V)\right) ^2 }{|V|}} \end{aligned}$$

(23)

$$\begin{aligned} \mu _{f_i}(T_{x_j})-\mu _{f_i}(T-T_{x_j})\ge {} \Delta \mu _{min} \end{aligned}$$

(24)

$$\begin{aligned} \sigma _{f_i}(T_{x_j})\le {} \sigma _{max}\end{aligned}$$

(25)

$$\begin{aligned} H_0 : \mu _{f_i}(T_{x_j})= {} \mu _{f_i}(T-T_{x_j}) \end{aligned}$$

(26)

FD-ASE algorithm performs the dimensionality reduction in the feature vector. Considering the contribution of each feature to increment the fractal dimension of the data set, the algorithm finds dependence relationship between attributes and determines a set of independent ones, discarding the others. The approach of attribute forward inclusion is used to recreate the data set. The fundamental idea is to calculate the partial correlation fractal dimension of data set projections, integrating more and more attributes until its calculated correlation fractal dimension is achieved. More details of this algorithm is found in [5].

1.4 Sequential selection methods

SFS [17] algorithm starts from an empty set. It sequentially adds the feature that results in the highest fitness value when combined with the features that have already been selected and moves towards. The search continues until full feature set is obtained. In SFS, it is unable to remove features that become obsolete after the addition of other features. SBS [17] starts from the full set and sequentially removes the feature that results in the best possible fitness value. SBS works best when the optimal feature subset has a large number of features. Reevaluation of the usefulness of a feature after it has been discarded is not possible in SBS.