Diatom identification including life cycle stages through morphological and texture descriptors

Elliptical Fourier descriptors

We obtain EFD of the contour using the method described by Kuhl & Giardina (1982). An implementation of EFD is available in BielStela (2017). Taking a contour image as the starting point, we calculate the Freeman chain code. Then being a_i the ith element in the Freeman chain code we obtain: (1)${\displaystyle \Delta x_i=sgn\left(6-a_i\right)sgn\left(2-a_i\right)}$ (2)${\displaystyle \Delta y_i=sgn\left(4-a_i\right)sgn\left(a_i\right)}$ (3)${\displaystyle \Delta t_i=1+\left(\frac{\sqrt[]{2}-1}{2}\right)\left(1-{\left(-1\right)}^{a_i}\right)}$

Δx_i and Δy_i are the changes in the x, y projections of the chain code. Δt_i is the modulus of the segment between two points i and i + 1. View Fig. 4 for more information.

Then we calculate the Fourier coefficients of the x and y projections of the Freeman chain code. a_n, b_n, c_n, d_n represent the n harmonic Fourier coefficients and T the perimeter. (4)${\displaystyle a_n=\frac{T}{2n^2{\pi }^2}\sum _{p=1}^K\frac{\Delta x_p}{\Delta t_p}\left[\right.cos\frac{2n\pi t_p}{T}-cos\frac{2n\pi t_{p-1}}{T}\left]\right.}$ (5)${\displaystyle b_n=\frac{T}{2n^2{\pi }^2}\sum _{p=1}^K\frac{\Delta x_p}{\Delta t_p}\left[\right.sin\frac{2n\pi t_p}{T}-sin\frac{2n\pi t_{p-1}}{T}\left]\right.}$ (6)${\displaystyle c_n=\frac{T}{2n^2{\pi }^2}\sum _{p=1}^K\frac{\Delta y_p}{\Delta t_p}\left[\right.cos\frac{2n\pi t_p}{T}-cos\frac{2n\pi t_{p-1}}{T}\left]\right.}$ (7)${\displaystyle d_n=\frac{T}{2n^2{\pi }^2}\sum _{p=1}^K\frac{\Delta y_p}{\Delta t_p}\left[\right.sin\frac{2n\pi t_p}{T}-sin\frac{2n\pi t_{p-1}}{T}\left]\right..}$

Finally, it was empirically found that the first 30 coefficients provide an accurate approximation to the contour.

The amplitude of the nth harmonic can be calculated as: (8)${\displaystyle am{p}_n=\frac{1}{2}\sqrt{a_n^2+b_n^2+c_n^2+d_n^2}}$

Figure 5 shows an example of the reconstruction of a contour using EFD when a different number of coefficients are used.

Phase congruency descriptors

The method of phase congruency (PC) is based on the concept that all Fourier components are in phase in the areas where signal changes occur. In the case of the images, these zones correspond to the edges, corners and textures of the objects. Therefore, the method seeks to obtain the maximum phase components in the Fourier domain. The main advantage of this method has to do with the fact that it is very robust to changes in lighting and contrast. This is due to the fact that the method works with the phases of the Fourier components but not with their amplitude.

Phase congruency was previously used as a preprocessing stage to contour extraction. Sosik & Olson (2007) found that a simple threshold-based edge detection applied to phase congruency provides excellent results to an ample set of phytoplankton images. Verikas et al. (2012) applied preprocessing through phase congruency-based methods for the enhancement of image edges. Phase congruency, as well as the gradient operator, is sensitive to noise. Li et al. (2006) performed a comparison of phase congruency with Canny edge detector and observed that the former is able to extract more detailed information than the later. They applied a contour length filter to remove the noise. In this paper, PC was not applied for contour improvement but for the extraction of discriminant descriptors.

Phase congruency (PC) based descriptors are calculated as in Verikas et al. (2012) by computing mean and standard deviation from phase congruency maximum (M) and minimum (m) momentum images (described in Kovesi, 2003 and available in Kovesi, 2000). Those images combine the phase congruency information of each orientation. At the end, we have four phase congruency descriptors.

The phase congruency is obtained according to the following expression: (9)${\displaystyle PC=\frac{\sum _{\theta }\sum _n{w}_{\theta }\left(x\right)\lfloor A_{n\theta }\left(x\right)\Delta {\Phi }_{n\theta }\left(x\right)-T_{\theta }\rfloor }{\sum _{\theta }\sum _n{A}_{n\theta }\left(x\right)+ϵ}}$ (10)${\displaystyle \Delta {\Phi }_{n\theta }\left(x\right)=cos\left({\varphi }_{n\theta }\left(x\right)-{\overline{\varphi }}_{n\theta }\left(x\right)\right)-\mid sin\left({\varphi }_{n\theta }\left(x\right)-{\overline{\varphi }}_{\theta }\left(x\right)\right)\mid }$ where θ indicates the orientation, A_nθ(x) and ϕ_nθ(x) the amplitude and phase angle respectively used with the frequency component n, orientation θ and location $x.{\overline{\varphi }}_{\theta }$ is the amplitude of the average phase angle in the orientation θ, w_θ a frequency parameter in the orientation θ and ϵ a constant to prevent division by zero. T_θ is the noise estimated in the orientation θ that must be suppressed.

Maximum and minimum momentum images are calculated as in Eqs. (11) and (12) as a function of the PC. (11)${\displaystyle M\left(x\right)=\frac{1}{2}\left[c+a+\sqrt{b^2+{\left(a-c\right)}^2}\right]}$ (12)${\displaystyle m\left(x\right)=\frac{1}{2}\left[c+a-\sqrt{b^2+{\left(a-c\right)}^2}\right]}$ where a, b, c are: (13)${\displaystyle a=\sum _{\theta }{\left[P{C}_{\theta }\left(x\right)cos\left(\theta \right)\right]}^2}$ (14)${\displaystyle b=2\sum _{\theta }\left[P{C}_{\theta }\left(x\right)cos\left(\theta \right)\right]\left[P{C}_{\theta }\left(x\right)sin\left(\theta \right)\right]}$ (15)${\displaystyle c=\sum _{\theta }{\left[P{C}_{\theta }\left(x\right)sin\left(\theta \right)\right]}^2}$

where PC_θ(x) is the phase congruency value at orientation θ: (16)${\displaystyle P{C}_{\theta }\left(x\right)=\frac{\sum _n{w}_{\theta }\left(x\right)\lfloor A_{n\theta }\left(x\right)\Delta {\Phi }_{n\theta }\left(x\right)-T_{\theta }\rfloor }{\sum _{\theta }\sum _n{A}_{n\theta }\left(x\right)+ϵ}}$

Figure 6 shows an example of a diatom and its corresponding m and M images.

Gabor filters

Gabor based descriptors are calculated by the same method as in Bueno et al. (2017) and originally described in Fischer et al. (2007). The implementation can be found in Cristobal, Fischer & Redondo (2016). First we calculate the log-Gabor filters as Gaussians shifted from the origin at different scales, s, and orientations, t, and they are applied to the input images. The formulation of the log-Gabor filters is: (17)${\displaystyle G_{\left(s,t\right)}\left(\rho ,\theta \right)=exp\left(-\frac{1}{2}{\left(\frac{\rho -{\rho }_s}{{\sigma }_{\rho }}\right)}^2\right)exp\left(-\frac{1}{2}{\left(\frac{\theta -{\theta }_{\left(s,t\right)}}{{\sigma }_{\theta }}\right)}^2\right)}$

where (ρ, θ) are the log-polar coordinates, the number of scales is S = 4 and the number of orientations is O = 6. Thus, s ∈ {1, …, S} and t ∈ {1, …, O} indexes the scale and the orientation of the filter, respectively. And (ρ_s, θ_(s,t)) are the coordinates of the center of the filter; (σ_ρ, σ_θ) are the angular and radial bandwidths in ρ and θ (see Fischer et al., 2007) for more details.

Then first and second order statistics are acquired for every sub-band. Gabor feature reduction was obtained using correlation for removing redundant information. With this reduction procedure we finally obtain a 177 Gabor feature vector from the original 1460 one. Figure 7 shows an example of log-Gabor filtering (four bands. G1, G2, G3 and G4) applied to a diatom example.

During the feature selection process other alternative descriptors were considered. We did some experiments with SURF features (Bay, Tuytelaars & Van Gool, 2006) but the results did not improve the final performance. Figure 8 presents the results of applying the Relieff algorithm to the full set of features selected, including SURF. Due to the low scores of SURF such features were discarded.

Supervised

An extensive range of supervised classifiers like nearest-neighbor, Supported Vector Machines (SVM), Random Forest and Bagging Trees, the quadratic Bayes normal classifier and the Fisher classifier were compared. The best results were obtained with KNN and SVM. Both algorithms need previous training. This training was carried out by selecting a small subset from the image dataset as training data and the rest as test data. Thus, a 10-fold cross-validation (10fcv) scheme was followed, where 10 image samples are used as the validation set and the remaining data as the training set. This was repeated $C_{10}^n$ times, where n is the total number of images in each dataset, and $C_{10}^n$ is the binomial coefficient. This is done to divide the original sample on a validation set of 10 samples and a training set. Finally, the arithmetic mean of the results of each iteration was performed to provide a single and final result.

• KNN: refers to K-Nearest Neighbors (Friedman, Bentley & Finkel, 1976). This classification algorithm assigns a class to each sample, choosing between the class of the K nearest neighbors, giving a confidence value of the assigned class. The distance between different elements can be computed as a simple Euclidean distance. KNN is one of the most common and straightforward classification methods. The procedure is highly dependent on the value of K which is usually determined empirically.

• SVM: refers to Support Vector Machine. It was first proposed by Vapnik (1999). This method determines the hyperplane that best divides the data into the different classes. The algorithm increases the dimensionality of the features space, so a non-linear classification problem can be solved linearly. The distance between the hyperplane and the training data is called the functional marging, which is used as a confidence interval of classification results (Wang et al., 2017).

Non supervised clustering

Three different clustering algorithms were selected for non-supervised classification: K-means, Hierarchical Agglomerative Clustering and BIRCH. Although the methods are labeled as unsupervised, they can actually be considered as semi-unsupervised because the number of clusters was identified with the number of classes.

• K-means: K-means algorithm (Lloyd, 1982) separates all the data into K clusters where the distance between the data and the cluster centroid (mean of all the data in the cluster) is minimized. The centroids are initialized at random points and they are updated after a feature is assigned to a cluster. As a result of the algorithm, the data space is partitioned into Voronoi cells, i.e., a diagram that divides the space in a given number of regions. A related method to K-means is the so-called K-medoids (Park & Jun, 2009). Unlike the K-means, K-medoids chooses datapoints (medoids) and uses squared Euclidian distance to define distance between data points.

• Hierarchical Agglomerative Clustering (Manning, Raghavan & Schutze, 2008): This algorithm starts with a cluster for each observation. These clusters are successively merged together, minimizing a distance function between the clusters. The process ends when the predefined number of clusters has been reached.

• BIRCH: refers to Balanced Iterative Reduced Clustering using Hierarchies (Zhang, Ramakrishnan & Livny, 1996). This algorithm is divided into four different phases. In phase 1 a Clustering Feature Tree (CF tree) is built with all the data. The second phase scans the CF tree to remove outliers and group crowded subclusters into larger ones to get a new smaller CF tree. This step is optional. Phase 3 is a global clustering procedure applied to the reduced CF tree. The last step is optional and it refines the clustering results. It obtains new clusters using the centroids of phase 3 as seeds for the clustering algorithm.

Clustering validation metrics

In order to validate clustering results five different metrics were considered (Vinh, Epps & Bailey, 2010; Kassambara, 2017).

• Adjusted RAND index: This index measures the similarity between the clustering labels assignment and the given ground truth. Its value varies in the range [−1,1], where 1 is perfect score.

• Silhouette: This metric evaluates the similarity between an element and others members of the same cluster compared to members of other clusters. This metric can indicate how compact a cluster is and how it is separated from other clusters. It can take on a value between −1 and 1. Measures close to 1 mean well defined clusters.

• Adjusted Mutual Information (AMI): This metric compares the agreement between clustering class assignments and the ground truth classes. It can take a value between −1 and 1 and values close to 1 indicate significant agreement.

• Homogeneity: It measures if a cluster consists only of member of the same class or not.

• Completeness: It measures if all the members of a class are assigned to the same cluster.