Fast Adapting Mixture Parameters Schemes for Probability Density Difference-Based Deformable Model

Abstract

This paper presents a new region-driven active contour using the pdf difference to evolve. The pdf estimation is done via a new and fast Gaussian mixture model (GMM) parameters updating scheme. The experiments performed on synthetic and X-ray images have shown not only an accurate contour delineation but also outstanding performance in terms of execution speed compared to the GMM estimation based on EM algorithm and to non-parametric pdf estimations.

1 Introduction

Image segmentation is of major importance in the field of computer vision. It is the process of dividing a digital image into something that is more meaningful and easier to analyze by the image analysis stages. Image segmentation can reveal to be hard for numerous reasons: Firstly, partitioning the image into non overlapping regions and extracting regions of interest require a tradeoff between the computational efficiency of the involved algorithm, its degree of automation and the accuracy of its outcomes. Secondly, image noise linked to the image acquisition, and poor contrast are very difficult to reckon with segmentation algorithms without the user interacting [2]. So, designing a robust and efficient segmentation method is still not trivial and difficult for practical applications. Among segmentation methods, the active contour is one of the most successful ones. Broadly speaking, an active contour is a curve that evolves from an initial position, under some constraints and energies to match the desired features. This curve can be edge-driven or region-driven depending on the nature of these energies. The region-driven ones tend to rely on global information to guide the contour evolution. Hence, the inner and the outer region defined by the model are considered, which leads to better handling of noise and smooth or vanished boundaries as well as sensitivity to initial conditions [8]. Global information can be established by considering the probability density function (pdf) of regions intensities and then, are based on the probability theory. Image segmentation techniques built on probability theory have been extensively used and can be considered as clustering or classification problems [5]. They must return necessary knowledge that would enable to assign a label to each pixel. Most of them merely obtain a label map that is inferior to active contours outcomes, which have smooth and regular boundaries of connected regions [5]. Many statistical modeling can be applied to distribute image features. We can cite among others, parametric modeling by mixtures. The most commonly used mixture model is the GMM [12, 17]. GMMs have been extensively employed to this aim, due to their simplicity and ease of implementation. Associating statistical approaches and active contours method bring several advantages as the merits of both are combined. Unfortunately, this association leads to prohibitive computation cost [11]. That is why a particular attention should be directed to how this combination is achieved to take advantage from it by finding strategies to speed up the contour evolution. In this paper, a novel region-driven parametric active contour is proposed, whose evolution is based on the differences between the estimated pdf inside and outside the model curve. The pdf estimations are carried out by the mean of GMMs, whose parameters are computed in a new recursive scheme to alleviate the computation load and accelerate the parameters estimations.

The structure of the remainder is as follows. In Sect. 2, we introduce the mathematical foundation of active contours and the adaptive mixtures. Our method for object extraction via the new pdf estimation is presented in Sect. 3. The experimental results are shown in Sect. 4. We draw the main conclusions in Sect. 5.

2 Background

2.1 Active Contours

The active contour, also called “snake”, is a curve $\mathbf c (s)= [x(s), y(s)]'$, $s =[0 \ 1]$ which evolves towards image features to minimize the following energy [9]

$$\begin{aligned} E(\mathbf c )=\int _0^1( E^{int}(\mathbf c (s))+E^{ext}(\mathbf c (s)))ds \end{aligned}$$

(1)

Where s is the curvilinear abscissa, $E^{int}$ the internal energy and $E^{ext}$ the image energy responsible for driving the contour toward edges. The minimization of Eq. (1) leads to the iterative solution of Eq. (2).

$$\begin{aligned} \left\{ \begin{array}{l} x_t=(A+\gamma I_d)^{-1}(\gamma x_{t-1}+\nabla E^{ext}_x(x_{t-1},y_{t-1}))\\ y_t=(A+\gamma I_d)^{-1}(\gamma y_{t-1}+\nabla E^{ext}_y(x_{t-1},y_{t-1})) \end{array} \right. \end{aligned}$$

(2)

If the model is made of N nodes, then, A is a $N*N$ matrix, $I_d$ an identity matrix sized as A, $\gamma$ an evolution coefficient, $x_t$ and $y_t$ are the model nodes coordinates at the iteration t. $\nabla E^{ext}_x(x_t,y_t)$ and $\nabla E^{ext}_y(x_t,y_t)$ are the external forces of the input image at the model nodes locations in the x and y direction respectively. In addition to the external force proposed in [9], a variety of external forces have been proposed to improve the snakes performance. They can be classified as dynamic and static [10]. The dynamic forces are those which depend on the snake position and, as a result, change as the snake deforms. In turn, the static ones are those that are computed once for all and remain unchanged.

Pressure Forces. The pressure force given by Eq. (3) [3] is an inflation/deflation dynamic force if the model is considered as a balloon.

$$\begin{aligned} F_B(\mathbf c (s^*))=k_p{.}{\varvec{n}}(\mathbf c (s^*)) \end{aligned}$$

(3)

${\varvec{n}}(\mathbf c (s^*))$ is the normal unit vector to the curve at $\mathbf c (s^*)$ and $k_p$ a weight. By introducing an individual pressure weight $k_p$ for each node as done in [6] for example, the model is strengthened, regarding to the initialization issue, since some parts of it can inflate/deflate independently from the other ones [7]. The individual $k_p$ may be released to various forms. An example is proposed in [1] (Eq. (4)), where $k_p$ is introduced as a difference between the pdfs estimated inside and out side the model. This difference is evaluated at the node position grey level value.

$$\begin{aligned} k_p=p(z_{s^*}/O)-p(z_{s^*}/B) \end{aligned}$$

(4)

$z_{s^*}$ is the node position grey value, whereas $p(z_{s^*}/O)$ and $p(z_{s^*}/B)$ are the conditional pdf of the object (O) and the background (B) respectively. The problem here is to accurately estimate the pdf inside and outside the active contour to achieve a successful progression to the boundaries.

2.2 Probability Density Function Estimation

To exploit the weight of Eq. (4), one should have a good estimation of the pdf. To this end, the GMM is used here. Expectation Maximization (EM) algorithm is usually employed to compute the mixture parameters [4]. However, numerical procedure with EM algorithm can reveal to be very expensive [16]. That is why an adaptive updating is chosen instead of the EM one. It relies on the works of [16] and later [14] to achieve such update. The following set of recursive equations for normal components presented in [14] are used to the aim of parameters vector computation $\theta$. For K components, $\theta$ consists of 3K elements (mixture proportions, means and standard deviations); $\theta _n =\lbrace \pi _n^1, \mu _n^1, \sigma _n^1, ..., \pi _n^K, \mu _n^K, \sigma _n^K\rbrace$.

$$\begin{aligned} \left\{ \begin{array}{l} \beta _{(n)}^{(i)}=1/n\\ \rho _{n+1}^{(i)}=\frac{\pi _{n}^{(i)}\phi ^{(i)}(x_{n+1})}{\hat{F}_{n}(x_{n+1}) }\\ \pi _{n+1}^{(i)}=\pi _{n}^{(i)}+\beta _{(n)}^{(i)}(\rho _{n+1}^{(i)}-\pi _{n}^{(i)})\\ \mu _{n+1}^{(i)}=\mu _{n}^{(i)}+(\pi _{n}^{(i)})^{-1}\beta _{(n)}^{(i)}\rho _{n+1}^{(i)}(x_{n+1}-\mu _{n}^{(i)})\\ \sigma _{n+1}^{2(i)}=\sigma _{n}^{2(i)}+(\pi _{n}^{(i)})^{-1}\beta _{(n)}^{(i)}\rho _{n+1}^{(i)}((x_{n+1}-\mu _{n}^{(i)})^2-\sigma _{n}^{2(i)}) \end{array} \right. \end{aligned}$$

(5)

$\rho$ denotes the membership function, $x_{n}$ the $n^{th}$ observation, $\theta _n$ the parameters estimates after n observations, $\phi ^{(i)}$ the normal distribution with parameters $\theta _n^{(i)}$ and $\hat{F}_n=\sum _{i=1}^K \pi ^{(i)}_n\phi ^{(i)}$ the estimated pdf after n observation. The superscript (i) indicates that the calculation is done for the $i^{th}$ mixture component.

3 Fast Mixture Updating Scheme for Contour Evolution

The proposed active contour uses probabilistic pressure forces as external forces, where $k_p$, as is [1], are computed as the differences between the estimated pdfs inside and outside contour curve. During the model progression, the pdf estimates change as the inner and the outer regions change with the deformations. When the inner region is considered, from an iteration t to $t+1$, new pixels get inside the curve while some others are left outside. The idea is to remove from the pdf estimate $\hat{F}_t$, represented by its parameters, the gray levels contribution of the pixels that have been left out. This is done in order to get new mixture parameters defining $\hat{F}^-_{t}$, the intermediate pdf estimate. These parameters are once again updated with the new gray levels values that have been added inside, to get $\hat{F}_{t+1}$, the new pdf estimate for $t+1$. In other words, if M pixels $y_j, j=1\dots M$ have been left out in $t+1$, and the inner region in t consists of n pixels (n observation), $\hat{F}^-_{t}$ is obtained by M successive updates achieved with Eq. (6).

$$\begin{aligned} \left\{ \begin{array}{l} \beta _{(n-j+1)}^{(i)}=1/(n-j+1) \\ \rho _{n-j}^{(i)}=\frac{\pi _{n-j+1}^{(i)}\phi ^{(i)}(y_j)}{\hat{F}_{n-j+1}(y_j) }\\ \pi _{n-j}^{(i)}=\pi _{n-j+1}^{(i)}-\beta _{(n-j+1)}^{(i)}(\rho _{n-j}^{(i)}-\pi _{n-j+1}^{(i)})\\ \mu _{n-j}^{(i)}=\mu _{n-j+1}^{(i)}-(\pi _{n-j+1}^{(i)})^{-1}\beta _{(n-j+1)}^{(i)}\rho _{n-j}^{(i)}(y_j-\mu _{n-j+1}^{(i)})\\ \sigma _{n-j}^{2(i)}=\sigma _{n-j+1}^{2(i)}-(\pi _{n-j+1}^{(i)})^{-1}\beta _{(n-j+1)}^{(i)}\rho _{n-j}^{(i)}((y_j-\mu _{n-j+1}^{(i)})^2-\sigma _{n-j+1}^{2(i)}) \end{array} \right. \end{aligned}$$

(6)

After the updating performed with Eq. (6) for $j=1\dots M$, the recursive operations of Eq. (5) are applied then on the pixels gray level values that have been added inside the curve, beginning with $\hat{F}^-_{t}$ ($n\longleftarrow n-M$, $\pi _{n-M}^{(i)}$, $\mu _{n-M}^{(i)}$, $\sigma _{n-M}^{(i)}$, $i=1\dots K$). Furthermore, the intersection of the two sets representing the added and the subtracted pixels gray levels values, could be not empty, which means that same gray levels are going to be processed twice, once by Eq. (6) and then by Eq. (5). This situation can be avoided by removing the set intersection to prevent unnecessary operations, as the pdfs are smooth-wise gray levels values counting. This can be explained by an example: if a gray level value g is found m times in the subtracted pixels and k times in the added ones, so this gray level value will be processed $m-k$ times by Eq. (5) or Eq. (6), for $k<m$ or $m<k$ respectively. Moreover, if L pixels of the same gray level value g, have to be processed in the same iteration, the updates are performed once for all the L pixels. Indeed, if we assume that for a big n, $\hat{F}$ does not significantly change when adding few gray level values, then $\rho$ will remain utmost the same and kept unchanged in the L iterations. The same assumption is done for $\beta$. As for L iterations, $\beta$ grows from 1/n to $1/(n+L-1)$, we choose to use the median value $\beta _{n+(L-1)/2} =1/(n+(L-1)/2)$. An example of the updating development for the parameter $\pi$ is given in Eq. (7) where $\beta _{n+(L-1)/2}$ is noted $\beta$

$$\begin{aligned} \left\{ \begin{array}{l} \pi _{n+1}^{(i)}(1)=\pi _{n}^{(i)}+\beta ^{(i)}(\rho _{n+1}^{(i)}-\pi _{n}^{(i)}) \quad { 1^{st}~\texttt {pixel}} \\ \pi _{n+1}^{(i)}(2)=\pi _{n}^{(i)}+2\beta ^{(i)}(\rho _{n+1}^{(i)}-\pi _{n}^{(i)})+(\beta ^{(i)})^2(\rho _{n+1}^{(i)}-\pi _{n}^{(i)}) \quad {2^{nd}~\texttt {pixel}}\\ \pi _{n+1}^{(i)}(3)=\pi _{n}^{(i)}+3\beta ^{(i)}(\rho _{n+1}^{(i)}-\pi _{n}^{(i)})-(\beta ^{(i)})^2(\rho _{n+1}^{(i)}-\pi _{n}^{(i)})- \\ (\beta ^{(i)})^3(\rho _{n+1}^{(i)}-\pi _{n}^{(i)}) \ \ {3^{rd}~\texttt {pixel}}\ \\ \pi _{n+1}^{(i)}(L)=\pi _{n}^{(i)}+L\beta ^{(i)}(\rho _{n+1}^{(i)}-\pi _{n}^{(i)})+O(\beta ) \quad {L^{th}}~\texttt {pixel} \end{array} \right. \end{aligned}$$

(7)

$O(\beta )$ is a polynomial of $\beta$ with degrees going from 2 to L. $O(\beta )$ is negligible as $\beta \ll 1$. The same results are obtained for the other mixtures parameters. Updates for the L pixels of the same gray level value g are then given by Eq. (8):

$$\begin{aligned} \left\{ \begin{array}{l} \beta ^{(i)}=1/(n+(L-1)/2) \\ \rho _{n+1}^{(i)}=\frac{\pi _{n}^{(i)}\phi ^{(i)}(g)}{\hat{F}_{n}(g)} \\ \pi _{n+1\rightarrow n+L}^{(i)}=\pi _{n}^{(i)}+L\beta ^{(i)}(\rho _{n+1}^{(i)}-\pi _{n}^{(i)})\\ \mu _{n+1\rightarrow n+L}^{(i)}=\mu _{n}^{(i)}+L(\pi _{n}^{(i)})^{-1}\beta ^{(i)}\rho _{n+1}^{(i)}(g-\mu _{n}^{(i)})\\ \sigma _{n+1\rightarrow n+L}^{2(i)}=\sigma _{n}^{2(i)}+L(\pi _{n}^{(i)})^{-1}\beta ^{(i)}\rho _{n+1}^{(i)}((g-\mu _{n}^{(i)})^2-\sigma _{n}^{2(i)}) \end{array} \right. \end{aligned}$$

(8)

The index ${n+1\rightarrow n+L}$ means that the update is computed as a contribution of all the L pixels. Thus, updating $\hat{F}$ no longer relies on the observations number (pixels number) as suggested by Eq. (5) but only on the gray level values of the pixels involved in the computation. It is worth to note that the same procedure is applied for the background (the outer region), to update the outer pdf estimate. These modified recursions give a convenient way for computing quickly the mixture parameters as the model curve progresses.

Once the pdf estimation is performed, the contour evolution is carried out by the following iterative equations where $k_p$ is the pdf estimation difference for each node position gray level value.

$$\begin{aligned} \left\{ \begin{array}{l} x_t=(A+\gamma I_d)^{-1}(\gamma x_{t-1}+ k_p{.}N_x)\\ y_t=(A+\gamma I_d)^{-1}(\gamma y_{t-1}+ k_p{.}N_y) \end{array} \right. \end{aligned}$$

(9)

$N_x(N_y)$ is the normal unit vector components in the x(y) direction.

During the active contour initialization, the first round of pdf estimation is performed with EM algorithm to get quickly the first pdf estimate. For the other rounds, the proposed updates are carried out as the active contour evolves in the image domain.

Algorithm Complexity. Assume that the number of the processed pixels in one iteration is equal to $N_p$ represented on $L_g$ gray levels, then, each of the pdf parameters $\rho$, $\pi$, $\mu$ and $\sigma$ updates requires $L_g\times K$ operations for the proposed pdf parameters updating, while it requires $N_p\times K$ operations for the recursive model provided in [14]. Then, a gain of execution time with our approach is obviously obtained since $L_g\le N_p$.

4 Experiments

As first experiment, we have compared the proposed pdf estimation method to the EM-based one. As illustrated in Fig. 1, the implementation results of Eq. (9) on a $250\times 250$ synthetic image corrupted by a Gaussian noise with $\sigma =20$, show a successful contour delineation for the two methods. However, the proposed pdf update have reduced the computation time by $50\%$ compared with the EM algorithm. One can note that in addition of a good performance in object contour extraction, thanks to a good pdf estimation, as shown in Fig. 1, the strength of the proposed method is related to the computation time. Indeed, for comparison purpose, the proposed model is faster than the non-parametric pdf estimation-base active contour with Parzen [13] and with the Averaged shifted histogram ASH [15], where the execution time was reduced by $50\%$ for ASH, and more than $95\%$ for Parzen. Another experiment shows capabilities of the proposed active contour in boundaries extraction in very noisy images (Fig. 2).

Fig. 1.

Active contours results and pdfs estimations.

Fig. 2.

Final contours on images corrupted with pronounced Gaussian noises. A: $\sigma =30/50$ for object/background. B: $\sigma =50/30$ for object/background

Fig. 3.

The proposed model application to an X-ray image. A: Defect histogram and defect pdf estimate. B: Background histogram and background pdf estimate

The last experiment consists in applying the proposed model on an X-ray image shown in Fig. 3 which contains a region of interest (ROI) of welded joint X-ray image to segment, the initial and the final contours, and the histograms of the extracted defect and background depicted with the pdf estimates. The background of the weld X-ray images can be modeled by more than one component because of the high illumination non-uniformity characterizing such images [11]. Here the number of background components is taken equal to 3. The bad quality of the image shown in Fig. 3 does not prevent our model to extract the weld defect indication successfully. In conclusion, the main advantage of our method consists in reducing the computation time. Indeed, applying an active model evolving with EM-based pdf estimate on the image of Fig. 3 has slowed the progression down 2 times compared to the proposed model progression speed.

5 Conclusion

In this paper, a new probabilistic active contour is presented. The contour progression based on a new pdf updating scheme of GMM parameters shows to be effective. Indeed, the contour delineation is achieved successfully despite the noise and the images bad quality. Moreover, the convergence was faster compared to others pdf-based methods, which ascertains the opportunity of exploiting the proposed pdf estimation via our new parameters updating scheme in such active contour model. For further works, we plan to extend the proposed updating to other mixtures of distributions related to normal one.