### Bilateral filter

Bilateral filter was developed by Tomasi and Manduchi [

13]. The basic idea of bilateral filter is to replace a pixel value in an image by a weighted mean of its neighbors, which the weights depend on both the spatial distance and the intensity distance [

14][

15]. There are many types of bilateral filters depending on the choice of weighting functions. What we develop in this paper is based on the Gaussian bilateral filter [

13][

16]. For Gaussian bilateral filter, it can be expressed mathematically as [

13][

17]

where

is the output pixel value,

*J* (

*Y*) is the input pixel values,

*X* and

*Y* are the coordinates vectors,

*σ*
_{
d
}
^{2} and

*σ*
_{
r
}
^{2} are the parameters controlling the fall-off of weights in spatial and intensity domains, respectively,

*N* (

*X*) is a spatial neighborhood of pixel

*J* (

*X*), || || is Euclidean distance, C is used for the normalization and is expressed as [

13][

17]

In the above equation, when *X* and *Y* are 2-D vectors, the bilateral filter is called 2-D bilateral filter, which can be used to reduce the noise in 2-D images.

Bilateral filter is a good choice for image de-noising because it is stable and simple. The effectiveness of bilateral filter lies in the combined use of the domain filter, which is used to enforce spatial closeness by weighting pixel values with coefficients that fall off with distance [18], and the range filter, which assigns greater coefficients to those neighbouring pixels with light intensity that is more similar to the centre pixel value [18]. In bilateral filter, the choice of the parameters *σ*
_{
d
}
^{2} and *σ*
_{
r
}
^{2} is very important. If their values are set too high, the filter will act as a smoothing filter and will blur the edges. If their values are set too low, the noise cannot be removed. Generally speaking, the choice of *σ*
_{
d
}
^{2} and *σ*
_{
r
}
^{2} depends on the variance of the noise. Based on the research in [17], the optimal *σ*
_{
d
}
^{2} is relatively insensitive to noise variance while the optimal *σ*
_{
r
}
^{2} changes significantly as the noise standard deviation changes [17]. [17] also demonstrates that *σ*
_{
r
}
^{2} and noise variance are linearly related to a large degree. The research in [17] is based on additive noise model. If bilateral filter is applied to speckle noise, the relationship between *σ*
_{
r
}
^{2} and noise variance will not be established because speckle noise is multiplicative noise. In order to reduce the speckles in ultrasound images effectively, we develop speckle reducing bilateral filter.

### Speckle reducing bilateral filter

Generally speaking, noise can be modelled by an additive model or a multiplicative model. Additive noise model is the simpler case of the two and can be described by a linear model

where *J* (·) is the noised image, *I* (·) is the original image and *n* (·) is the noise. Multiplicative noise is generally expressed by a multiplicative model

It is well known that multiplicative noise appears much worse in bright image regions than dark regions since it multiplies the gray intensities.

Speckle noise is generally treated as multiplicative noise and can be modelled using equation (4). Thus, compared with other types of noise, speckle noise is generally difficult to be removed. Our research below shows that the conventional bilateral filter described in equation (1) and (2) generally gets bad results when it is applied to speckle reduction directly. Thus, the bilateral filter described in (1) and (2) needs improvement or enhancement so that it can be applied to reduce the speckles in images effectively. In order to do this, we will first analyze the behavior of
of the bilateral filter in equation (1) in a homogenous region for both additive noise and multiplicative noise, then we will propose an adaptive bilateral filter for speckle reduction.

Let *J* (*Y*) and *J* (*X*) be two different pixels from image *J.* If *J* is corrupted by additive noise, then we can use equation (3) to compute the difference between these two pixels

|| *J* (*Y*) - *J* (*X)* || *= ||I* (*Y*) + *n(Y*) - *I* (*X*) - *n* (*X*)|| (5)

If both *J* (*Y*) and *J* (*X*) are from the same homogenous region, then we have *I* (*Y*) = *I* (*X*), thus

Equation (6) means that the difference between any two pixels from the same homogenous region is only related to the difference of the noise. If *J* is corrupted by multiplicative noise, then we can use equation (4) to compute the difference between these two pixels. From equation (4), we have

|| *J* (*Y*) - *J* (*X)|| =* || *I* (*Y*) * *n* (*Y*) - *I(X*) * *n* (*X*)|| (7)

Similarly, if both *J* (*Y*) and *J* (*X*) are from the same homogenous region, then we have *I* (*Y*) = *I* (*X*), thus

From equation (8), we can understand that the difference between two pixels in the same homogenous region(in the image corrupted by multiplicative noise) is not only related to the difference of the noise. It also depends on the intensity of the region. As is seen in equation (8), the difference is big when the intensity of the region is big while the difference is small when the intensity of the image is small.

The above analysis shows the bilateral filter described in (1) and (2) is not suitable for removing speckle noise, which is multiplicative noise. The reason lies in the difference of the corrupted image has different distributions in different homogenous regions. For example, if

*σ*
_{
r
}
^{2} is fixed in the processing, when

*σ*
_{
r
}
^{2} is set to be big, the edge in lower intensity regions will be removed, while the noise can’t be removed in the higher-intensity regions when

*σ*
_{
r
}
^{2} is set to be small. Thus, in order to develop an effective bilateral filter to remove the speckle, we need to develop a new representation of the difference. Dividing each side of equation (

8) by ||

*J* (

*X*)||, in the homogenous regions, we have

Equation (

9) shows that the normalized difference is only related to the noise and doesn’t depend on the intensities of the region. Thus, the proposed adaptive bilateral filter can be expressed as follows

Bilateral filter is famous because it is non-iterative, however, the non-iterative bilateral filter doesn’t yield good results. In order to improve its effectiveness, we use iterative bilateral filter. The basic idea of iterative bilateral filter is to use the filtered image obtained by equation (

10) as the input of equation (

1) and implement it many times, the mathematical expression is as follows:

Where
. Experiments show that iterative bilateral filter gives much better results than the non-iterative bilateral filter.

### Cattle follicle segmentation

In order to analyze and monitor the reproduction of cattle, the acquisition of some quantitative parameters is very important. These parameters include diameters, areas and perimeters of the follicles. These parameters can be used to monitor the development and maturity of follicles. In order to get these parameters, we need to segment the follicles.

Many image segmentation methods have been proposed, which includes histogram based methods, edge detection based methods, region based methods, active contour model based methods, etc. Active contour model based methods have drawn a lot of attention in the past decade because of their significant advantages. In this paper, we adopt active contour model based method for the segmentation of the follicles. An active contour or a snake [

19] is defined as a controlled continuity contour that is attracted to salient image features. However, there are some disadvantages related to the original model. Thus, many improved active models have been proposed based on the original model. The gradient vector flow (GVF) model is one of them [

20][

21]. GVF model is designed to overcome one of the disadvantages of original model, i.e. the original model is sensitive to the initialization of the snake. In GVF model, GVF fields are computed by another diffusion process, which can be implemented by minimizing the following energy function [

20][

21]:

where

*g* is a decreasing function of the edge-force magnitude and is defined as follows:

Here *k* is a non-negative smoothing parameter for the field (*u*, *v*) *.* The functional described by equation (15) smoothes the force field (*u*, *v*) only when the edge strength is low. Solving the energy functional optimization problem in (14), we can obtain the generalized gradient vector flow, which can be used as external forces that attract the snake to the follicle boundary [20][21].

GVF provides external forces for a snake model, we also need internal forces to smooth the contour. In this paper, we use B-spline to represent the contour instead of the real internal forces. B-spline has been used in snake model in several applications and get pretty good results [22][23][24]. Let the control points be denoted by *P*
_{0} through *P*
_{
m
}. The knot-value sequence is a non-decreasing sequence of knot values *t*
_{0} through *t*
_{
m
}
_{+4}, and *Q*
_{
i
} is a curve segment defined by control points *P*
_{
i-
}
_{3}, *P*
_{
i-
}
_{2}, *P*
_{
i-
}
_{1}, *P*
_{
i
} and blending functions *B*
_{
i-
}
_{3,4}, *B*
_{
i
}
_{-2,4}, *B*
_{
i
}
_{-1, 4}, *B*
_{
i
}
_{, 4} (*t*) as follows [22]:

*Q*
_{
i
} (*t*) = *P*
_{
i-
}
_{3} · *B*
_{i-3, 4} + *P*
_{i-2} · *B*
_{
i
}
_{-2,4} + *P*
_{
i
}
_{-1} · *B*
_{
i
}
_{-1,4} + *P*
_{
i
} · *B*
_{
i
}
_{,4} (*t*)(16)

where 3

*≤ i ≤ m* and

*t*
_{
i
}
*≤ t ≤ t*
_{
i+
}
_{1}. The blending functions can be obtained using recursion as follows [

22]:

When *p* =4, we obtained the blending function of cubic splines. The GVF snake with B-spline is called B-spline GVF snake [23][24][25].

For the segmentation of the follicles, we initialize the B-spline GVF snake using a circle inside each follicle. The circle is represented by B-spline and the number of control points is set to 48 in this paper. Then, starting from the initial contour, the GVF is used to drive the contour to the boundary of the follicle. The evolution of the contours is the same as that in the B-spline GVF snake in single scale proposed by [24].