Keywords

1 Introduction

About 20% of the European population suffers from allergic reactions accompanied by symptoms such as swelling of the eyes, the urge to sneeze, or shortness of breath [16]. One of the main carriers of allergens causing such allergic reactions is the male gametophytes of seed plants, which is also known as pollen. Pollen with allergens trigger an over-sensitive reaction of the immune system by stimulating the production of the antibody immunoglobulin E. Symptoms appear when the allergen gets in contact with the antibody. Those symptoms are aggravated by air pollutants such as nitrogen dioxide, ozone, or fine dust (e.g. diesel soot particles) [5, 9, 10]. As the most common chronic disease, allergies not only reduce the quality of life, but also have a negative impact on the socio-economy. For example, additional costs of 56 billion US dollars were incurred due to allergies in the United States from 2002 to 2007 [3].

In Germany, the pollen count has been systematically recorded by the German Pollen Information Service Foundation (Stiftung Deutscher Polleninformationsdienst) since 1983. Allergologically relevant as well irrelevant pollen are collected at stationary measurement stations where experts identify and quantify them manually by using microscopes [2]. In cooperation with the German Weather Service, the collected data is used to determine and forecast the pollen count for the following days. Likewise, the pollen load for those areas where there are no monitoring stations can be predicted [20, 25].

Existing automated pollen classification systems use images that are generated by the means of mass spectrometry [21], vibration spectroscopy [11, 14], light microscopy, or electron microscopy [4]. While the former two methods require chemical fingerprints for the identification of the pollen, the latter two classify them on the basis of their texture [15] or their external morphology [23]. In this regard, the segmentation of the individual pollen and the subsequent manual generation of the characteristics constitute the biggest challenge while very noisy backgrounds and irrelevant particles make the classification process more difficult. In [12], these problems were circumvented by using a convolutional neural network for the classification.

Against this background, the existing pollen classification systems require a high amount of information (e.g. electron microscopy images [4]) and, thus, depend on measurement techniques which are able to provide the required data. Inevitably, the existing methods do not only involve high costs, but also need more space for the measurement devices. At the same time, in order to develop and conduct effective treatment and prevention measures in a targeted manner, it is necessary that the information about the local exposure profiles for the individual patients is as detailed and as differentiated as possible. Hence, there is a high demand for reliable, compact, and inexpensive systems that allow for real-time detection and classification of pollen grains at the individual level. For this purpose, a new measurement technique is developed which uses binary 2D projections to identify and classify individual pollen grains. Given that binary 2D projections require relatively little information, the present paper explores the question of to what extent it is possible to detect and classify pollen grains using this new approach. To verify the measurement method, artificial pollen is created in advance and the measurement methodology is simulated. The simulated data are used for later classification.

The present paper proceeds as follows: In Sect. 2, after introducing the allergologically relevant pollen, a simulation methodology is presented with which the respective pollen can be created. In the subsequent part, the new measurement methodology and the simulation are described. Section 4 again, presents the characteristics that are used while Sect. 5 presents the results. Finally, the concluding section sums up the findings with respect to the research question and provides an outlook on the future use of binary 2D projections for the detection and classification of individual pollen grains.

2 Used Pollen

In this section, the allergologically relevant pollen are presented, followed by a method to simulate them artificially.

2.1 Allergologically Relevant Pollen

In this study, seven allergologically relevant pollen, which have been registered by the Foundation German Pollen Information Service since 1983, are examined. They include birch pollen (Betula), which accounts for around 50 percent of the total distribution of all allergologically relevant pollen, alder pollen (Alnus), which makes up almost 20 percent, grass pollen (Poaceae) with about 15 percent, hazel (Corylus) with almost 10 percent, and mugwort (Artemisia), rye (Secale) and ragweed (Ambrosia), whose combined share amounts to 3 percent [25].

2.2 Artificial Pollen

The allergically relevant pollen all have a spherical shape [1, 6]. Thus, for their mathematical description, spherical coordinates (SC) are used. Thereby, a point of such a spherically shaped object is specified as the distance r from the centroid of the object to the two angles ϕ and θ [24]:

$$\displaystyle \begin{aligned} P = (r,\phi,\theta)_{SC}, \quad \! \begin{array}{ll} r &\geq 0, \\ {} 0^\circ &\geq \phi \geq 360^\circ, \\ {} 0^\circ &\geq \theta \geq 180^\circ. \end{array} \end{aligned}$$

A constant distance r and the variable angles ϕ and θ result in an ideal spherical surface (sphere). The allergologically relevant pollen, however, do not have an ideal spherical surface, which means that the distance r is not constant. Instead, depending on the two angles ϕ and θ, the distance r can be described as follows:

$$\displaystyle \begin{aligned} r &= f(\phi,\theta) = \hat{r}(\theta) + g(\phi,\theta),\\ \hat{r}(\theta) &= \sqrt{\frac{1}{\frac{cos^2(\theta)}{b^2}+\frac{sin^2(\phi)}{a^2}}}. \end{aligned} $$

The term \(\hat {r}(\theta )\) describes the rotationally symmetrical flattening and elongation of the ideal sphere so that with the two parameters a and b, a flattened or elongated rotational ellipsoid is created. The parameters a and b correspond to half the length of the polar axis DPol (see Table 1) and the corresponding transverse axis DLateral. The length of the transverse axis can be derived from the following relationship:

$$\displaystyle \begin{aligned} D_{\text{Lateral}} &= \frac{D_{\text{Pol}}}{P_{\text{Forml}}}.\\ \text{Thereby:} \quad a &= \frac{D_{\text{Pol}}}{2},\quad b = \frac{D_{\text{Lateral}}}{2}. \end{aligned} $$
Table 1 Quantitative description of the allergologically relevant pollen with a comparative depiction of natural and artificially created pollen
Full size table

Hereby, the PForml is the pollen shape index which describes the ratio of the length of the polar axis to the largest transverse axis of a pollen grain. All other deviations of the pollen under examination from the rotational ellipsoid, such as the shape of the aperture or the surface structure, are expressed with the term g(ϕ, θ). Each of the allergologically relevant pollen has a different surface variation. For this reason, this term is derived on the basis of the geometric measurement data in the existing literature [6] through an empirical analysis. The artificially created pollen are shown in Table 1.

3 Measurement Principle and Simulation

The used measuring principle is presented in this section. Afterwards the boundary conditions for the simulation are explained.

3.1 Measurement Principle

For the detection of pollen and fine dust particles, a particle sensor based on silicon photonics is used. Thereby, through an opening, the particle flow is run into the photonic chip. On its way, the particle flow passes by an optical sensor with which the particles are detected and analysed. The analysis is carried out in the first iteration with the discrete change in the optical signal. That is, the detection and analysis take place if the particles interrupt the light signal of the sensor when passing by (principle of a light barrier). By accumulating this signal over time, a binary two-dimensional projection of the particle is produced. Thereby, the resolution of the projection can be adjusted with any technically feasible finite number of receivers and transmitters. The measurement principle is depicted in Fig. 1.

Fig. 1
figure 1

Measurement principle of the light barrier: Measurement section without particles (left) and with particles (right)

Full size image

3.2 Simulation

To create the simulation environment, some assumptions need to be made in order to keep the model as simple as possible. Firstly, the simulation grounds on the assumption that the particle flow is a laminar flow. This, again, implies the following: The particles run through the measuring section with a constant velocity v0 (Eq. 1) and that the velocity vvertical(t) of the particles which is vertical to the flow velocity equals zero (Eq. 2). Moreover, the rotational direction (Eq. 3) does not change during the flow.

$$\displaystyle \begin{aligned} v(t) &= \textit{constant} = v_0,& t > t_0,{} \end{aligned} $$
(1)
$$\displaystyle \begin{aligned} v_{\text{vertical}}(t) &= 0,& t > t_0,{} \end{aligned} $$
(2)
$$\displaystyle \begin{aligned} \dot{\alpha}_{x}(t) &= \dot{\beta}_{y}(t) = \dot{\gamma}_{z}(t) = 0,& t > t_0.{} \end{aligned} $$
(3)

Hereby, the parameters \(\dot {\alpha }_{x}(t)\), \(\dot {\beta }_{y}(t)\) and \(\dot {\gamma }_{z}(t)\) specify the angular velocity around the axes x, y and z in the three-dimensional space.

Secondly, it is assumed that the light source is an infinitesimally small light source. That is, the receivers are illuminated with either 100 or 0 percent of the light. This, again, means that the recipients have only two modes: “recipient on” and “recipient off”.

Finally, it is assumed that all other interactions, such as the reflection of the optical signal, are irrelevant and, hence, neglected. Based on these assumptions the equations of particle motion can be simplified into the following equations:

$$\displaystyle \begin{aligned} P_x(t,r,\phi,\theta) &= P_{x,0}(r,\phi,\theta) = \textit{constant},\\ P_y(t,r,\phi,\theta) &= P_{y,0}(r,\phi,\theta) = \textit{constant}, \\ P_z(t,r,\phi,\theta) &= P_{z,0}(r,\phi,\theta)-v_0\cdot t. \end{aligned} $$

The parameters Px, Py and Pz represent a point of the particle in the Cartesian coordinate system. The parameters with the additional index 0 describe the points at the starting position at time t0. The derivation of the starting points in the Cartesian coordinate system from spherical coordinates is defined as follows:

$$\displaystyle \begin{aligned} P_{x,0}(r,\phi,\theta) &= r \cdot cos(\theta) \cdot cos(\phi) + x_{C,\text{Start}},\\ P_{y,0}(r,\phi,\theta) &= r \cdot cos(\theta) \cdot sin(\phi) + y_{C,\text{Start}},\\ P_{z,0}(r,\phi,\theta) &= r \cdot sin(\theta) + z_{C,\text{Start}}. \end{aligned} $$

In these equations, the terms xC,Start, yC,Start and zC,Start specify the starting position of the center of gravity at the time.

Figure 2 exemplarily illustrates the simulation for two points in time by using the measurement principle described above.

Fig. 2
figure 2

Visual depiction of the measurement method by using artificial pollen Betula in the starting position t = t0 (a) and while running through the measurement section at time t = t1 (b)

Full size image

4 Feature Extraction

The Zernike moments, topological features and the Fourier descriptors are used for feature extraction. The procedures are described in the following section.

4.1 Zernike Moments

Zernike moments belong to the orthogonal moments and are based on mapping an image onto a set of orthogonal, complex polynomials. Due to orthogonality, the image is represented with a minimum of information redundancy and the Zernike moments are both rotation and flip invariant [22]. In order to make the Zernike moments also invariant with respect to translation and scaling, they need to be normalized to their basic polynomials in a pre-processing step [13]. Another property of the Zernike moments includes information compression. Information compression means that low frequencies of a pattern are encoded into lower order of moments. Given these properties of Zernike moments, already a small number of descriptors is sufficient to avoid noise and deformation [19].

Because Zernike moments are specified in polar coordinates (r, ψ) while the image intensities are defined in Cartesian coordinates (x, y), the image intensity function on the unit disk needs to be transformed into polar coordinates. The center of the image corresponds to the coordinate origin of the new domain. For a digital image I of size N × M, the transformed coordinates can be noted as follows:

$$\displaystyle \begin{aligned} r_{x,y} &= \sqrt{\left(\frac{2 \cdot x-N+1}{N-1}\right)^2+\left(\frac{2 \cdot y-M+1}{M-1}\right)^2},\\ \psi_{x,y} &= tan^{-1} \left( \frac{\left(2 \cdot y-M+1 \right) \cdot \left(N-1\right)}{\left(2 \cdot x-N+1 \right) \cdot \left( M-1\right)} \right). \end{aligned} $$

Orthogonal Zernike moments (Eq. 4) of the order p with the repetition q are defined by the complex, conjugated Zernike polynomials \(V^*_{pq}\) (Eq. 5). One has order p ∈{0, …, ∞} and the repetition q attains a positive integer value or zero if, and only if, the repetition is smaller than the order (p ≥|q|) and the difference amount is supposed to be even (|q − p| = even). The complex, conjugated Zernike polynomials are defined by the radial Zernike polynomials (Eq. 6).

$$\displaystyle \begin{aligned} Z_{pq}&=\frac{p+1}{n} \cdot \sum_{x=0}^{N-1} \sum_{y=0}^{M-1} V^*_{pq}\left(r_{x,y}, \psi_{x,y} \right) \cdot I\left(x, y \right),{} \end{aligned} $$
(4)
$$\displaystyle \begin{aligned} V^*_{pq}\left(r_{x,y}, \psi_{x,y} \right)&=S_{pq}\left(r_{x,y}\right) \cdot e^{-i\cdot q \cdot \psi_{x,y}},{} \end{aligned} $$
(5)
$$\displaystyle \begin{aligned} S_{pq}\left(r_{x,y}\right)& = \sum_{k=0}^{\frac{p+q}{2}} \left( -1\right)^k \cdot \frac{\left(p+q\right)! \cdot r_{x,y}^{p-2\cdot k}}{k! \cdot \left(\frac{p+q}{2}-k\right)! \cdot \left(\frac{p-q}{2}-k\right)!}.{} \end{aligned} $$
(6)

Table 2 lists the calculated Zernike moments up to the 10th order. The amplitude and phases of these complex moments serve as features for the pollen classification.

Table 2 Calculated Zernike moments up to the 10th order
Full size table

4.2 Topological Features

Topological features refer to the global structural properties of a binary image. As such, they are inherently invariant to both rotation, translation and shear. Some topological features, such as the Euler number – which indicates the number of connected regions minus the number of their holes – or the form factor are also scaling-invariant. Due to this peculiar properties, the type of features remains unchanged even with significant image deformations [8].

In order to determine the global topological features of a binary image, the image is searched for local pixel configurations. Pixel configurations are 2 × 2-sized filter masks, which are also called bit quads. Sixteen different bit quads can occur for binary images. The bit squads, to each of which an index k is assigned, are listed in Table 3 [17].

Table 3 2 × 2 pixel configurations for binary images and the respective index k [17]
Full size table

In order to calculate the topological features, the number of the occurring bit quads Qk in the image is calculated. To this end, a mask M with the power of two is defined:

$$\displaystyle \begin{aligned} M = \begin{bmatrix} 2 & 1 \\ 8 & 4 \\ \end{bmatrix}. \end{aligned}$$

The image I is then convolved with the mask M resulting in

$$\displaystyle \begin{aligned} F = I*M. \end{aligned}$$

This calculation process yields a 4-bit grayscale image F, whereby its gray values fij corresponds to the bit quads

$$\displaystyle \begin{aligned} Q_k = \sum_{i=0}\sum_{j=0} \begin{cases} 1 \quad f_{ij} = k\\ 0 \quad f_{ij} \neq 0 \end{cases}. \end{aligned}$$

For the pollen classification, the bit quads are categorised into congruence classes TAc (c ∈{0, …, 4}). Congruence classes contain only those bit quads that can be converted into one another by rotation:

$$\displaystyle \begin{aligned} TA_0 &= Q_0, \\ TA_1 &= Q_1+Q_2+Q_4+Q_8,\\ TA_2&=Q_3+Q_5+Q_{10}+Q_{12},\\ TA_3&=Q_7+Q_{11}+Q_{13}+Q_{14},\\ TA_d&=Q_6+Q_9,\\ TA_4& = Q_{15}. \end{aligned} $$

In addition to the congruence classes, the topological features surface area A, perimeter P and circularity C are calculated by deriving them from the congruence classes as follows:

$$\displaystyle \begin{aligned} A&=\frac{1}{4}\cdot \left(TA_1+2\cdot TA_2+3\cdot TA_3+4\cdot TA_4+2 \cdot TA_d \right),\\ P&=TA_1+TA_2+TA_3+TA_4+2\cdot TA_d, \\ C &= \frac{P^2}{4\cdot \pi \cdot A}. \end{aligned} $$

4.3 Fourier Descriptors

Fourier descriptors are global shape features that represent a two-dimensional closed contour. Each feature characterises a peculiarity of the entire form. If one of these features is changed or omitted, the entire form changes. Fourier descriptors, which basically reproduce the shape in the frequency space, have the advantage that a small number of low frequencies are sufficient to reproduce shape properties [26]. By conducting certain modifications, the shape features can be made invariant to translation, rotation, and scaling [7].

To calculate Fourier descriptors, the two-dimensional closed contour B is mapped into a one-dimensional space. By choosing the scalar representation of the distance from the center of gravity C, invariance with respect to translation can be achieved. For the contour coordinates (xn, yn) ∈ B and the center of gravity (xc, yc) ∈ C, the equation takes the following form:

$$\displaystyle \begin{aligned} r_n = \sqrt{\left(x_n-x_c\right)^2+\left(y_n-y_c\right)^2} \quad \! \begin{array}{l} n= 0 \dots N-1 \end{array}. \end{aligned}$$

The frequency coefficients are calculated by using the discrete Fourier transformation as follows:

$$\displaystyle \begin{aligned} c_k = \frac{1}{N} \cdot \sum_{n=0}^{N-1} r_n \cdot e^{-i\cdot2\cdot \pi \cdot \frac{n}{N}} \quad \! \begin{array}{l} k = 0 \dots N-1 \end{array}. \end{aligned}$$

In order to determine the invariance towards rotation and scaling, the amplitudes of the coefficients are divided by the first coefficient c0. For the classification of pollen, it suffices (see in Sect. 5 Fig. 6) to determine the first 50 Fourier descriptors:

$$\displaystyle \begin{aligned} FD = \left[\frac{|c_1|}{|c_0|},\ldots, \frac{|c_{50}|}{|c_0|}\right]^T. \end{aligned}$$

5 Experiment and Results

Two different data sets, each of which contained 500 × 7 artificial pollen of the seven allergologically relevant pollen types, are simulated. In the first data set (name: data set without rotation), random pollen sizes in accordance with [6] are used. Moreover, all pollen are subjected to a fixed rotational direction. In the second data set (name: data set with rotation), the size of the pollen is changed based on [6] and the rotational direction of the pollen is varied.

By applying the measurement principle described, the artificial pollen are used to create binary 2D projections for different receiver distances (see Fig. 3). The receiver distance is varied from 0.1 to 5 μm in an interval of 0.1 μm.

Fig. 3
figure 3

Binary 2D projections of an Alnus pollen for different receiver distances: (a) 0.1 μm, (b) 5 μm

Full size image

In order to classify the pollen, a support vector machine with a Gaussian kernel is used. Thereby, for each type of pollen, 334 × 7 labeled pollen are used for training and 166 × 7 unlabeled pollen for testing. The Zernike moments, the topological features, and the Fourier descriptors are initially used individually, but the latter two are combined later.

Figure 4 shows the percentage of incorrectly classified pollen (error rate) for the data set without rotation in relation to different Zernike orders. As the graph illustrates, with a high resolution (receiver interval of 0.1 μm), the error rate decreases, the higher the order is. At the 10th order, while the error rate is about 4% at high resolution, it increases to over 20% at low resolution (receiver interval 2 and 5 μm).

Fig. 4
figure 4

Percentage of incorrectly classified pollen (error rate) for Zernike orders up to the 10th order and for the receiver intervals 0.1, 2 and 5 μm in the data set without rotation

Full size image

The findings for the topological features are shown in Fig. 5. As shown in the graphs, for the data set without rotation, the error rate in relation to different orders of Zernike moments is significantly lower than for the data set with rotation. For the former, the error rate remains below 5% even at low resolutions. In fact, with a receiver interval of 0.1 μm, 100% of the pollen are classified correctly. In contrast, with an error rate of up to about 40%, the data set with rotation is significantly more prone to errors. In comparison, the error rate for the data set without rotation alternates by approximately 10% for larger receiver intervals.

Fig. 5
figure 5

Percentage of incorrectly classified pollen (error rate) for receiver intervals of up to 5 μm for the data set with and without rotation by using topological features

Full size image

In Fig. 6a and b, it is depicted how the error rate proceeds when using Fourier descriptors as a discriminatory feature for the pollen classification. As Fig. 6a shows, the first six Fourier descriptors are already sufficient to correctly classify 100% of the pollen without rotation at a receiver interval of 0.1 μm. The error rate increases only minimally up to about 3% even if the receiver interval is high. For the data set with rotation, by using ten Fourier descriptors, an error rate of 0% is achieved with a resolution of 0.1 μm receiver interval. However, also for this data set, the error rate increases, the bigger the receiver interval is.

Fig. 6
figure 6

Error rate of pollen classification in relation to the number of Fourier descriptors for the receiver intervals 0.1, 2 and 5 μm for both data set without (a) and with rotation (b). In (c), the error rate is illustrated for the data set with rotation while using Fourier descriptors and topological features combined

Full size image

Given that the data set with rotation generally involves higher error rates, the topological features of the data set were combined with the Fourier descriptors. The findings depicted in Fig. 6c demonstrate that this combination significantly lowered the error rate so that it remains below 10% for up to 30 Fourier descriptors. Though, the error rate slightly increases if the number of Fourier descriptors increases. Apart from that, bigger receiver distances usually correspond with significantly higher error rates.

6 Conclusions

Existing pollen classification systems use the texture of digital pollen recordings for the discrimination of pollen species. The present paper, however, explores the question to what extent the contour of the pollen suffices to classify them. For this purpose, a new measuring principle is developed which creates binary 2D projections of the pollen grains by using light barriers.

In order to test the measurement method in advance, binary 2D projections of the pollen grains are created in a simulation by using this method. For this purpose, two data sets of artificial pollen are used: In the first data set, the shape and size of the pollen are varied while the rotational direction is the same for all pollen grains. By doing so, this data set serves the purpose of examining whether it is possible to classify the pollen grains on the basis of the chosen properties by using binary 2D projections only. In the second data set, in contrast, all properties, including the rotational direction of the pollen, are varied. Thereby, the classification of the pollen is conducted by using a support vector machine with a Gaussian kernel. The simulated binary 2D projections of the pollen are represented by Zernike moments, topological features, and the Fourier descriptors.

The findings suggest that Zernike moments are not suitable for classifying pollen, for large receiver distances, they involve a high error rate of more than 20% even for the data set without rotation. The topological features yield better results for the data set without rotation. However, for the data set with rotation, an error rate of more than 30% occurs for large receiver intervals. In contrast, the Fourier descriptors provided the best results for both the data set with and without rotation. With a high resolution, i.e. a small receiver interval, even 100% of the pollen can be correctly classified by using just a small number of Fourier descriptors. Moreover, the performance of the shape-based Fourier descriptors can be improved by combining them with the structure-based topological features. By doing so, the error rate of classification can be reduced to less than 10% even for large receiver intervals.

Although the results are promising, no clear conclusions can be derived from the findings of this study as to whether it is possible to classify allergologically relevant pollen solely on the basis of their binary 2D projection. To obtain a higher validity of the results, it is necessary that studies are conducted which include both allergologically relevant and non-relevant pollen grains. Apart from that, the present study is based on simulated pollen. In order to examine to what extent the simulated pollen provide a good approximation of real pollen, further research are required.