Automatic identification of species with neural networks

Images

Image data were taken from two sources: natural history museum records, and online (Data S1). Each collection was analyzed according to the country of origin. Ichthyology collections from Colombia were compiled from the Instituto de Investigaciones Marinas y Costeras (INVEMAR), the Colección de Referencia Biología Marina Universidad del Valle (CRBMUV), and the Coleccion Ictiologica Universidad de Antioquia (CIUA). Ichthyology collections from Brazil were found in the Museu de Zoologia da USP (MZUSP), the Instituto Nacional de Pesquisas da Amazônia Manaus (INPA), and the Museu Nacional Rio de Janeiro (MNRJ). Image data from Spain came from the Museo Nacional de Ciencias Naturales Madrid (MNCN). We tested a data set that included a total of 740 species and 11,198 individuals of fish, plants, and butterflies. Fish specimen images were taken using a Canon EOS 6dD one-use camera with a 1,280 × 960 pixel resolution. A total of 697 fish species previously identified by experts, were photographed (see Fig. 1 for a subset of photographed species). Images of 32 plant species were downloaded from the Flavia database (2009) (http://flavia.sourceforge.net/) (see Fig. 2). Image data for 11 species of butterflies were downloaded from the MorphBank database (Erickson et al., 2007) (see Fig. 3).

System development

Based on pattern recognition theory (De Sá, 2001) and basic computer-processing pathways used in typical automated species identification systems (Gaston & O’Neill, 2004), we designed a system for automatic individual identification at the species level (Fig. 4). In a novel way, our system shared preprocess and extraction components with both training and recognition processes. Features of training images are used to build a model of the classification progress pattern after feature extraction. These features and the trained model were then recorded in the database and incorporated in the analysis of subsequent photos. This process used two types of data to model image features and resulted in better species identification results. Features can be found in Supplemental Information.

Image preprocessing

Image heterogeneity in terms of orientation, size, brightness, and illumination was common (Fig. 5.1). The image background was removed with Grabcut’s algorithm (Rother, Kolmogorov & Blake, 2004) (Fig. 5.2) and converted to grayscale (Fig. 5.3). Different filters were applied to improve the image by removing image noise; the filters used were smooth and median (Figs. 5.4 and 5.5), and the image was then reduced to one of two possible levels, 0 or 1 (Fig. 5.6). Next, the processed image was brought to a contour (Fig. 5.7) and then a skeleton (Fig. 5.8). All of these processes were performed for each taxonomic group using the image processing in MATLAB R2009b.

Feature extraction

A series of 15 geometrical, morphological, and texture features, that could be efficiently extracted with image processing and were unique to species, were used in our automatic identification system; these features can be efficiently extracted with image processing (Table 1).

Geometrical

Geometric features contain information about form, position, size, and orientation of the region. The following six geometric features were commonly used in pattern recognition.

1-Area was the total number of pixels of the specimen area, and was defined as: ${\displaystyle A\left(s\right)={\int }_x{\int }_{y}I\left(x,y\right)dydx}$ I(x, y) depended on the limits of the shape (Fig. 5.7).

2-Perimeter was the number of pixels that belonged to the edge of the region (Fig. 5.8). In other words, it was the curve that enclosed a region S, defined as ${\displaystyle P\left(s\right)={\int }_t\sqrt{x^2\left(t\right)+y^2\left(t\right)}dt.}$ 3-Diameter was the value that represented the diameter of a circle with the same area as the region. 4-Compatibility was the efficiency of the contour or perimeter P(s) that enclosed the area A(s) ${\displaystyle C\left(s\right)=\frac{4\pi A\left(s\right)}{P^2\left(s\right)}.}$ 5-Compactness was the efficiency with which area A(s) enclosed a speciment and was determined by P(s) ${\displaystyle Co\left(s\right)=\frac{P^2\left(s\right)}{4\pi A\left(s\right)}.}$ 6-Solidity was the scalar specifying the proportion of the pixels in the convex hull that were also in the region. This property was supported only for 2-D input label matrices.

Texture

Textures are important visual patterns for homogeneous description of regions. Intuitive measures provide properties such as smoothing, roughness, and regularity (Glasbey, 1996). Textures depend on the resolution of the image and can follow two approaches: statistical and frequency. We used the statistical approximation in which statistical values are analyzed first order (on the histogram) and second order (on the co-occurrence matrix).

Statistical first order was obtained from the gray level histogram of the image. Each value was divided by the total number of pixels (area) and had a new histogram representing the probability that a determined gray level was displayed in the region of interest.

The properties obtained were:

7-Median ${\displaystyle \mu =\sum _{x=1}^{n}xh\left(x\right).}$ 8-Variance ${\displaystyle {\delta }^2=\sum _{x=1}^n\left(x-{\mu }^2\right)h\left(x\right).}$ The second order of statistics were the matrix of spatial dependence of gray levels or co-occurrence matrices. Given a vector of polar coordinates, δ = (r, θ) we calculated the conditional probability that two properties appeared separated by a given distance δ, P_δ using an angle θ of −45 and a distance r equal to one pixel. The features that were extracted from this matrix were:

9-Uniformity ${\displaystyle \sum _{x=1}^n\sum _{y=1}^n{P}_{\delta }{\left(x,y\right)}^2.}$ 10-Entropy co-occurrence ${\displaystyle -\sum _{x=1}^n\sum _{y=1}^n{P}_{\delta }\left(x,y\right)log{P}_{\delta }\left(x,y\right).}$ 11-Homogeneity ${\displaystyle \frac{\sum _{x=1}^n\sum _{y=1}^n{P}_{\delta }\left(x,y\right)}{1+|x-y|}.}$ 12-Inertia ${\displaystyle \sum _{x=1}^n\sum _{y=1}^n{P}_{\delta }\left(x,y\right){\left(x-y\right)}^2.}$

Morphological

The morphological features were those that concentrate on the organization of pixels. They performed a comprehensive description of the region of interest. They fell into two categories: two-dimensional Cartesian moments and normalized central moments.

The two-dimensional Cartesian moments were variable at minor order, and were initiated at zero at higher orders. The moment of order p and q of a function I(x, y) was defined as: ${\displaystyle m_{pq}={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{x}^p{y}^{q}I\left(x,y\right)dxdy.}$ m_pq statistical moments, the parameters p and q denoted the order of the moment. When p = 0 and q = 0, which determined the center of mass or gravity of the overall function in binary images, the center of mass or gravity of the region under study was: ${\displaystyle \bar{x}=\frac{m_{10}}{m_{00}}\bar{y}=\frac{m_{01}}{m_{00}}.}$ The center of mass or gravity can defined the central moments that were invariant to displacement or translation of the image’s region of interest defined as: ${\displaystyle u_{pq}=\sum _x\sum _y{\left(x-\bar{x}\right)}^p{\left(y-\bar{y}\right)}^{q}I\left(x,y\right)\Delta A.}$ Where ΔA was the area of a pixel.

The normalized central moments were invariant to scale which was defined as: ${\displaystyle n_{pq}=\frac{u_{pq}}{u_{00}^{\gamma }}}$ where $\gamma =\frac{p+q}{2}\hspace{0.17em}\forall p+q⩾2$.

The above equations were defined by seven moments that were invariant to rotation, translation, and scale changes, known as the Hu invariant set of moments (Hu, 1962). In this study, we used the first Hu moment defined as:

13-Hu1 ${\displaystyle {\phi }_1=m_{20}+m_{02}.}$ Normalized central moments were generated by related moment invariants “AMI” (Flusser & Suk, 1993), based on the theory of algebraic invariants and invariants under general affine transformation. We used two of the four invariants associated with discriminant character moments defined as:

14-Ami1 ${\displaystyle I_1=\frac{u_{20}{u}_{02}-u_{11}^2}{u_{00}^4}.}$ 15-Ami2 ${\displaystyle I_2=\frac{u_{30}^2{u}_{03}^2-6u_{30}{u}_{21}{u}_{12}+4u_{30}{u}_{12}^3+4u_{21}^3{u}_{03}-3u_{21}^2{u}_{12}^2}{u_{00}^{10}}.}$ These moments enable a high degree of insensitivity to noise that is not altered by rotation, translation, or staggering.

The use of the above 15 features (Table 1) characterised the structure of the individual’s body, which was important for the identification at species level. We designed and realized automatic extraction algorithms to compute the values of these features so that all variables and features were calculated automatically.

Neural network

A neural network is defined as a parallel computer model composed of a large number of adaptive processing (neural) units which communicate via interconnections with variables. A multiple layer network has one or more layers (neurons) that enable the learning of complex tasks by progressively extracting more meaningful features from the input image patterns (Wu, 1997). Compared to other machine learning methods, neural networks learn slower but predict faster and have very good models for nonlinear data. The simple perceptron was assigned multiple inputs but generates a single output, similar to different linear combinations that depend on input weights and generated a linear activation function (Rosenblatt, 1958). Mathematically, the neural network was described with the following equation: ${\displaystyle y=\phi \left(\sum _{i=1}^n{w}_i*x_i+b\right)}$ w_i: weight vector, x_i: input vector, b: bias activation function.

A multilayer perceptron consisted of a set of source nodes containing one or more input layer and a set of hidden-node outputs. The input signal propagated through the network layer by layer (Zhang, Patuwo & Hu, 1998) (Fig. 6).

The neural network structure was composed of N inputs N = [N₁, N₂, …, N_n], a hidden layer h and an output vector S = [S₁, S₂, …, S_m]. Each S_i was assessed by a single step that transformed the vector S binary signal [0, 1]. A supervised training phase, or sigmoid activation, is based on the back propagation algorithm in which the weights and biases were updated in the direction of the negative gradient of the performance and then updated in the opposite direction (Werbos, 1974; Rumelhart, Hinton & Williams, 1986; Parker, 1987; Smith & Brier, 1996). The sigmoid activation function for the hidden layer and output layer was determined by the following equation: ${\displaystyle f\left(x\right)=\frac{1}{1+{\mathrm{e}}^{-x}}.}$ In this study, the number of input neurons was determined by the number of descriptors that were available in each pattern, which in this case was 15 (see variables section). The number of neurons in the hidden layer, h, was experimentally determined from the error set by comparing with the general training data of the ANN. The number of output neurons was determined by the number of species classified in each database.

To determine the optimal number of neurons given a data image, the relationship between the identification success rate and the number of neurons was explored. Figure 7 shows this relationship for the different configurations considered. We finally established our networks with 200 neurons for FC-MZUSP, 180 neurons for FC-INPA, 60 neurons for FC-MNRJ, 250 neurons for FC-INVEMAR, 60 neurons for FC-CIUA, 300 neurons for FC-CRBMUV, 250 neurons for FC-MNCN, 60 neurons for FLAVIA, and 35 neurons for BUTTERFLIES (see Table 2). The number of generations (i.e., a finite set of input patterns presented sequentially) for training and testing the ANNs was variable in the different collections between 50,000 and 140,000. It is evident that a high number of neurons and generations is required to process the information of the images in each collection.

The data set was randomly divided into 60–70–80–90% training images, resulting in 40–30–20–10% test images (Table 3). The results with the highest average accuracy for species identification were networks using 80–90% training and 20–10% test images. For these tests, the declared success rate was related to the number of species identified correctly. Recognition became more difficult with increased species number, as observed in the results from collections from MZUSP, INPA, INVEMAR, CRBMUV, and MNCN which averaged below 90% recognition.