Particle Detection Algorithms for Complex Plasmas

TL;DR

This paper evaluates various particle detection algorithms in complex plasmas, combining image processing techniques to improve accuracy, and provides insights applicable to other particle tracking fields like colloids and granular matter.

Contribution

It systematically compares algorithms and processing methods, including automatic threshold detection, to enhance particle position accuracy in complex plasma imaging.

Findings

Different algorithms have unique advantages depending on the scenario.

Automatic threshold detection improves robustness of particle identification.

Results are applicable to other fields involving image-based particle tracking.

Abstract

In complex plasmas, the behavior of freely floating micrometer sized particles is studied. The particles can be directly visualized and recorded by digital video cameras. To analyze the dynamics of single particles, reliable algorithms are required to accurately determine their positions to sub-pixel accuracy from the recorded images. Typically, straightforward algorithms are used for this task. Here, we combine the algorithms with common techniques for image processing. We study several algorithms and pre- and post-processing methods, and we investigate the impact of the choice of threshold parameters, including an automatic threshold detection. The results quantitatively show that each algorithm and method has its own advantage, often depending on the problem at hand. This knowledge is applicable not only to complex plasmas, but useful for any kind of comparable image-based particle tracking, e.g. in the field of colloids or granular matter.

Equations

$$\displaystyle w(r)=\left\{\begin{array}[]{ll}\frac{1}{2}-\frac{1}{2}\cos{\left(\frac{2\pi\left(r-\frac{N-1}{2}\right)}{N-1}\right)}&\mbox{if }\left|r\right|\leq\frac{N-1}{2}\\ 0&\mbox{else.}\\ \end{array}\right.$$

$$P_i\subset$$

$$\left(\left|P_i\right|=1\right)\vee$$

$$(\forall p_j\in P_i:$$

$$\exists p_k\in P_i\setminus \left\{p_j\right\}:d(p_j,p_k)<r)$$

$$\forall p_j\in P_i:$$

$$\forall p_k\in \left\{p:I(p)>I_{threshold}\right\}\setminus P_i:$$

$$d(p_j,p_k)\ge r$$

$$\left|P_i\right|\ge m_p$$

$$x_d\left(P_i\right)\ge m_x$$

$$y_d\left(P_i\right)\ge m_y$$

$$\frac{\left|P_i\right|}{x_d\left(P_i\right)\times y_d\left(P_i\right)}\ge m_d$$

$$\frac{1}{\left|P_i\right|}\sum _{p\in P_i}I(p)\ge m_{bd}$$

$$x_d(P_i):=\max \{x\left(p_j\right)-x\left(p_k\right):$$

$$p_j,p_k\in P_i\}$$

$$y_d(P_i):=\max \{y\left(p_j\right)-y\left(p_k\right):$$

$$p_j,p_k\in P_i\}$$

$$x_c(P_i)=\frac{{\sum }_{p\in P_i}x(p)\left(I(p)-I_{base}\right)}{{\sum }_{p\in P_i}\left(I(p)-I_{base}\right)}$$

$$y_c(P_i)=\frac{{\sum }_{p\in P_i}y(p)\left(I(p)-I_{base}\right)}{{\sum }_{p\in P_i}\left(I(p)-I_{base}\right)}$$

$$\frac{1}{2\pi {\sigma }_x{\sigma }_y}\exp \left(\frac{-1}{2}\left(\frac{(x-{\mu }_x{)}^2}{{\sigma }_x^2}+\frac{(y-{\mu }_y{)}^2}{{\sigma }_y^2}\right)\right)$$

$${\int }_{x_1}^{x_2}{\int }_{y_1}^{y_2}$$

$$\exp \left(\frac{-1}{2}\left(\frac{(x-{\mu }_x{)}^2}{{\sigma }_x^2}+\frac{(y-{\mu }_y{)}^2}{{\sigma }_y^2}\right)\right)dydx$$

$$=$$

$$\times \left({\int }_{y_1}^{y_2}\exp \left(\frac{-(y-{\mu }_y{)}^2}{2{\sigma }_y^2}\right)dy\right)$$

$$=$$

$$\times \left(\mathrm{erf}\left(\frac{y_1-{\mu }_y}{\sqrt{2}{\sigma }_y}\right)-\mathrm{erf}\left(\frac{y_2-{\mu }_y}{\sqrt{2}{\sigma }_y}\right)\right)$$

$$B_{noisy}=\max \left\{0,\min \left\{1,B+\frac{1}{\mathrm{SNR}}{B}_{noise}\right\}\right\}$$

$$\frac{p^{1-\frac{1}{p}}}{2\sigma \Gamma \left(\frac{1}{p}\right)}\times \exp \left(\frac{-{\sqrt{x^2+y^2}}^p}{p{\sigma }^p}\right)$$

$$\displaystyle d\left(\left(\begin{array}[]{cc}x_{i}\\ y_{i}\\ \end{array}\right),\left(\begin{array}[]{cc}x_{j}\\ y_{j}\\ \end{array}\right)\right):=$$

$$\max \left\{|x_i-x_j|,|y_i-y_j|\right\}$$

$$\displaystyle\forall i:d\left(\left(\begin{array}[]{cc}x_{i}\\ y_{i}\\ \end{array}\right),p_{0}\right)\geq d\left(p_{1},p_{0}\right)$$

$$\left\{\lceil x_0-\delta \rceil ,\lfloor x_0+\delta \rfloor ,\lceil y_0-\delta \rceil ,\lfloor y_1+\delta \rfloor \right\}$$

$$\underset{x,\sigma ,p,b}{\min }\sum _{i,j}{\left(I_{i,j}-A_{i,j}(x,\sigma ,p)-b\right)}^2$$

Full text

Particle Detection Algorithms for Complex Plasmas

Daniel P. Mohr Deutsches Zentrum für Luft- und Raumfahrt e. V., Institut für Materialphysik im Weltraum, 82234 Wessling, Germany. complex-plasmas.dlr.de   

Christina A. Knapek11footnotemark: 1

  

Peter Huber11footnotemark: 1

  

Erich Zaehringer11footnotemark: 1

Abstract

In complex plasmas, the behavior of freely floating micrometer sized particles is studied. The particles can be directly visualized and recorded by digital video cameras. To analyze the dynamics of single particles, reliable algorithms are required to accurately determine their positions to sub-pixel accuracy from the recorded images. Typically, straightforward algorithms are used for this task. Here, we combine the algorithms with common techniques for image processing. We study several algorithms and pre- and post-processing methods, and we investigate the impact of the choice of threshold parameters, including an automatic threshold detection. The results quantitatively show that each algorithm and method has its own advantage, often depending on the problem at hand. This knowledge is applicable not only to complex plasmas, but useful for any kind of comparable image-based particle tracking, e.g. in the field of colloids or granular matter.

Introduction

Particle detection in digital images is a crucial first step in the analysis of many-particle systems in the case that individual particles can be detected by direct optical measurements. Efforts to optimize particle detection can be found in a wide range of fields: in biophysics, single particle tracking is used to study the motion of particles (e.g. proteins, molecules or viruses) involved in cell membrane and intracellular activities (Saxton and Jacobson, 1997; Sbalzarini and Koumoutsakos, 2005; Chenouard et al., 2014). Particle detection and tracking from optical measurements is utilized in granular matter research (Tsai and Gollub, 2004; Harrington et al., 2014), and in colloidal physics, where the dynamics of systems of nano- to micrometer sized particles can be investigated by analyzing single particle motion from direct video microscopy (Crocker and Grier, 1996; Leocmach and Tanaka, 2013).

Complex plasmas (Fortov et al., 2005; Morfill and Ivlev, 2009; Ivlev et al., 2012) consist of micrometer sized particles injected into a low-temperature plasma composed of electrons, ions and neutral gas. These particles are large enough to be visible to digital cameras with an appropriate optics, and provide an excellent opportunity to study fundamental dynamics on the kinetic level of individual particles. In contrast to colloids, where particles are embedded in a viscous medium and therefore over-damped, complex plasmas are virtually undamped, and the time scale of dynamical processes is short and therefore easily accessible. The particles are typically illuminated with a sheet of laser light, and the reflected light can be observed with digital cameras. Since the particle distances are large (with a magnitude of several hundreds of micrometers), individual particles can be observed directly as mostly disjunct small groups of illuminated pixels on the camera sensor.

From those “blobs” of pixels, particle positions can be determined to sub-pixel accuracy – a necessity for the study of dynamics of single particles – with an adequately chosen algorithm.

By detecting individual particles, and tracing them through consecutive images (this is possible if the particle displacement between two images is small enough to allow for a unique assignment), velocities can be obtained. This method is called Particle Tracking Velocimetry (PTV), and has the advantage of more precise velocity measurements (Pereira et al., 2006) in contrast to Particle Image Velocimetry (PIV)(Williams, 2016), where only spatially averaged velocity vectors are obtained, especially in particle clouds too dense for single particle detection.

Complex plasmas are three-dimensional systems, and currently the interest in 3D optical particle diagnostics is growing (Jambor et al., 2016; Melzer et al., 2016). To triangulate the real position of a particle in 3D space, additional requirements are imposed on particle detection algorithms. Hence, we are also looking for algorithms for the detection of particles which are nearby each other on the image plane due to their overlapping motion in different layers. These algorithms can also be useful for particle tracking in systems with a high packing density.

With the methods presented in this paper we show that the commonly acquired accuracy can be exceeded without unnecessary increasing the complexity of the procedure. This is accomplished by combining simple image pre- and post-processing procedures with an improved version of the commonly used algorithm for blob detection, and to some extent by applying automatic threshold detection.

Usually, straightforward algorithms are used for blob detection (Crocker and Grier, 1996; Feng et al., 2007; Ivanov and Melzer, 2007), which is justified by the simple search feature and the typically low image noise. We show that this algorithm can be improved by generalizing it to blobs being not necessarily simply connected sets of pixels. Other more complex blob detection algorithms, such as SimpleBlobDetector (Bradski, 2000, SimpleBlobDetector) or MSER (Matas et al., 2002), did not turn out to be satisfactory in our case.

Though some of the techniques are well-known, a combination of them as well as an investigation of their individual influence on the accuracy of particle detection was not performed elsewhere to such an extent, especially for the typical particle shapes obtained in complex plasma experiments (for example, Feng et al. (2007) is mostly involved in examining one particular core algorithm without pre- and post-processing, while Ivanov and Melzer (2007) investigated some methods for pre- and post-processing, but does not combine the methods in the result).

Here, we not only investigate pre-processing, particle detection and post-processing in combination, but also take into account particle sizes and several kinds of image noise in our results. Additionally, we introduce Otsu binarization as an automated procedure. We also show that the choice of methods strongly depends on the image features (e.g. noise).

The paper is organized as follows: After a description of the general approach in section 1, section 2 shows how the artificial images are generated to test the quality of the algorithms. In sections 3, 4 and 5 the different steps of image processing and particle detection are presented in more detail, followed by some examples in section 6. Finally, the conclusion summarizes the paper in section 7.

1 General Approach

The process of obtaining particle coordinates from experimental data (images) can be divided into the following necessary steps:

Image acquisition

Get the image from real world.

Image processing

Prepare/enhance the image (e.g. by filtering).

Blob detection

Identify particle positions.

Postprocessing

Enhance found positions of particles (e.g. fitting).

Each step is an own field of research. In the following, they are explained in the depth necessary for this work.

1.1 Image Acquisition

Image acquisition is part of the experiment, and is only mentioned here for completeness, since the details of the experimental procedures go beyond the scope of this paper. To get good images, we need a proper illumination of the particles, a matched optical system, and an appropriate camera with an applicable storage system.

In this step, image noise is introduced. The sources are manifold, e.g. thermal behavior of the camera chip, noise of the involved electronics, defect pixels or radiation influencing the complete system. The noise causes uncertainties, which can be abstracted as additive white Gaussian noise and salt and pepper noise, superimposed on the pixels.

In the following, we assume a camera giving $8$ Bit gray scale images.

1.2 Image Processing

Preparing the image is a task extremely dependent on the blob detection algorithm to be chosen for the next step, e.g. an algorithm using edge detection will not work well if the edges of the blob are destroyed by applying a smoothing filter. In that case, a sharpening filter would be preferable.

One particle can be seen approximately as a point source of light, and the point spread function describes what we can expect to see on the image sensor. The point spread function defines how an ideal point source will be mapped by a system of optical components. In the case of point-like particles, the Airy disc (Airy, 1835) gives a good approximation of this mapping.

Optical side lobes of the point spread function can be reduced by a Hanning amplitude filter (a convolution with the Hann function) (Kumar et al., 2013). The Hann function, visualized in Figure 2, with the parameter $N$ for a point $r$ is given by:

[TABLE]

The parameter $N$ influences the width of the window. The Hanning amplitude filter is in principal a low-pass filter. This kind of filter passes signals with a spatial frequency less than the (user chosen) cutoff frequency, and can therefore reduce high-frequency image noise, e.g. Gaussian white noise.

This filter can easily be implemented by using template matching from the library opencv (Bradski, 2000).

In Figure 1, an example shows the effect of a Hanning amplitude filter.

Of course it is in general a good idea to use combined low-pass and high-pass filtering. High-pass filtering does the opposite to a low-pass filter: pass spatial frequencies below a cutoff and thus reduce image noise such as large-scale intensity gradients. But a high-pass filter can mask the behavior of a low-pass filter – e.g. the blurring of a low-pass filter would be reduced by a high-pass filter. Since we want to investigate the effect of specific filters, we do not want this masking. Usually, we do not observe low frequency noise in our images, and therefore omit high-pass filtering in this paper.

In general it should be mentioned, that Crocker and Grier (1996) and Ivanov and Melzer (2007) use a simple but effective filter, which behaves similar to a high-pass filter. They subtract a background calculated by a convolution with a boxcar kernel (moving average) from the image after low-pass filtering by a convolution with a Gaussian kernel.

1.3 Blob Detection

In the complex plasma community, a typical approach for blob detection is the moment method (Feng et al., 2007; Ivanov and Melzer, 2007), which is a simplified version of the approach by Crocker and Grier (1996) 111You can find many implementations in the internet, e. g.: www.physics.emory.edu/faculty/weeks/idl:

Find connected pixels brighter than a threshold (a particle). 2. 2.

Calculate center of every particle (position of a particle).

In the literature (Feng et al., 2007; Ivanov and Melzer, 2007), connected pixels are assumed to be a set, which is simply connected. More general, we now define a set of pixels $P_{i}$ belonging to a particle as:

[TABLE]

[TABLE]

[TABLE]

Here, $P_{i}$ is a set of pixels and represents the particle with the number $i$, $I(p)$ is the intensity of the pixel $p$, $I_{threshold}$ is the intensity of the threshold, $\left|P_{i}\right|$ denotes the cardinal number of $P_{i}$, $d(p_{j},p_{k})$ is the distance of the two pixels $p_{j}$ and $p_{k}$, $r$ is a search radius, $m_{p}$ is the minimum number of pixels a particle needs to be composed of, $x(p)$ is the $x$ coordinate of $p$, $y(p)$ is the $y$ coordinate of $p$, $m_{x}$ is the minimum length in $x$ direction in pixel, $m_{y}$ is the minimum length in $y$ direction in pixel, $m_{d}$ is the minimum density of a particle (density being the total number of pixels weighted by the area of the smallest rectangle envelope of $P_{i}$), $m_{bd}$ is the minimum brightness density. The brightness density is defined as the sum of all intensity values of the pixels in $P_{i}$, weighted by the total number of pixels in the set $P_{i}$.

The parameter $r$ allows to identify sets of pixels as a particle $P_{i}$ even if those pixels are not directly connected. For example, setting $r=1$ leads to a simply connected set as used in the mentioned literature, while setting $r=1.5>\sqrt{1^{2}+1^{2}}$ (assuming quadratic pixels with side length $1$) allows pixels in $P_{i}$ to be connected only by a corner. For larger values of $r$, the pixels in the set $P_{i}$ do not need to be simply connected at all. This can be used for compensation of pepper noise or intensity jitter. In addition, to be recognized as separate particles, the shortest distance between the particle contours of two neighboring particles must be $\geq r$.

The center $\left(x_{c}(P_{i}),y_{c}(P_{i})\right)$ can be calculated using the pixel positions and — as often done in the mentioned literature — the brightness of them:

[TABLE]

Here, $I_{base}$ gives an offset. In Feng et al. (2007) this offset is discussed and it is recommended to use $I_{base}=I_{threshold}$ to reduce the error.

Other blob detection algorithms (Bradski (2000, SimpleBlobDetector and MSER)) were tested, but proved to be unreliable and could only detect some of our largest particles. Since those algorithms increase the complexity and computation time without reaching the quality of our proposed blob detection method for the small particle images prevalent in complex plasmas, they were not investigated further.

1.4 Postprocessing

Since the blob detection is not an exact deconvolution, we are bound to have errors. To overcome this, we can fit a function to the approximate particle coordinates obtained from the blob detection. We now use the concept of a particle as an approximate point source of light, and the subsequent description of $P_{i}$ as a point spread function similar to the Airy disc (Airy, 1835). The latter can be approximated by a Gaussian or a generalized Gaussian point spread function (Claxton and Staunton, 2008) (see (1) in section 5), visualized in Figure 2.

In our procedure, we choose a generalized Gaussian point spread function and fit it to the approximate coordinates from the blob detection.

2 Simulated Images

To test our implementation we need well-defined, artificial images of particles. Here, the use of artificial images with well-defined particle positions is crucial to be able to calculate the deviation of tracked position to real position and thus to quantify the quality of our algorithms. The images are modelled after real-world experimental images of complex plasmas recorded by optical cameras.

The particles are represented by a bivariate normal distribution with a correlation of [math]:

[TABLE]

Since $(2\pi\sigma_{x}\sigma_{y})^{-1}$ is just a constant factor, it can be ignored. Furthermore, the image is scaled to values between [math] and $1$.

In a real-life camera image, the brightness of a pixel is the integration over time and space. Therefore, we integrate the intensity values over one pixel:

[TABLE]

[TABLE]

[TABLE]

Again, the constant factor $(\pi\sigma_{x}\sigma_{y})/2$ can be ignored, because the image is rescaled in the end. This procedure is repeated for each pixel.

Examples for artificial particle images are given in Figure 3. The figure also illustrates the impact of the given sub-pixel location of the particle center on the intensity distribution.

To be able to describe the strength of the noise by one single parameter, we create an additive white Gaussian noise (AWGN) with a mean of [math] and a standard deviation of $1$. We can scale the noise to the image by a signal to noise ratio $\operatorname{SNR}\,$ with $B$ a matrix representing the noise free image, $B_{noise}$ a matrix representing the noise and $B_{noisy}$ a matrix representing the image with noise:

[TABLE]

With this widely-used, simple noise model (e.g. it is often used in information theory Shannon (1948); Proakis (2001)), we can create a noise which behaves roughly similar to the thermal noise of camera sensors. In Pitas (2000, pages 43–44), this approach of cutting values is used to generate additive Laplacian noise222Because of a simple pseudo random number generator it was necessary to use a Laplacian instead of a Gaussian distribution in (Pitas, 2000, pages 43–44)..

Our simple $\operatorname{SNR}\,$ is consistent with the well-known Rose criterion (Rose (2013, page 97)), which states that a $\operatorname{SNR}\,$ of at least $5$ is necessary for a reliable detection. Due to this fact, Figure 4 does not show bars for $\operatorname{SNR}\,=5$.

By setting pixel intensities to [math] or $1$ with a given probability, we can add salt and pepper noise. This kind of noise simulates defective pixels usually present on typical camera sensors. Pixels can appear dark (“pepper”) or bright (“salt”), regardless of the exposure, e.g. due to errors in the pixel electronics. Bright pixels are easy to detect by taking dark-frame images (e.g. an image taken with covered lens), dark pixels can be detected with more effort by taking gray images. If a list of defective pixels is available, some cameras are able to correct these listed pixels by averaging the surrounding ones.

Though the occurence of excessive salt and pepper noise in an experimental setup should normally lead to an exchange of hardware, there are situations in which this is not an option. Good examples are experimental instruments in remote locations not accessible to technicians, e.g. satellites or sealed experimental devices on a space station, such as complex plasma microgravity research facilities (PK-3 Plus (Thomas et al., 2008), PK-4 (Pustylnik et al., 2016)). Here, cameras are exposed to higher levels of radiation, and pixel deterioration, causing salt and pepper noise, becomes an issue. To still obtain good scientific results over an extended period of time, one needs to handle such noise sources adequately during data analysis as long as it is feasible.

3 Preprocessing

(Image Processing)

Image preprocessing is not restricted to the use of general filters preserving the brightness distribution of particles, but can be extended to procedures for e.g. threshold detection, especially with regard to the requirements of the moment methods.

In the first step, the moment method needs a separation of the pixels belonging to particles, and pixels composing the background. Since our images represent particles illuminated with a laser, we can assume to have a bi-modal histogram.

This can be clustered for example by Otsu’s method (Otsu, 1979). This method separates the histogram of the image in $2$ classes – below and above the threshold – with the goal to minimize the variances of both classes. This leads to a maximal interclass variance. The image is then binarized according to the classes — pixels of one class are usually shown as white, and those belonging to the other class as black.

There are other thresholding techniques available (for an overview see e.g. Sezgin and Sankur (2004)). We use Otsu’s method since it is in the top ten ranking of Sezgin and Sankur (2004), one of the most referenced (therefore well-known), and implemented e.g. in opencv (Bradski, 2000).

Furthermore, a quick visual check of our example images with the tool ImageJ (Schneider et al., 2012) shows, that most available other techniques lead to erroneous binarizations, with background pixels becoming falsely detected as signals and set to white.

We analyzed one of the more promising methods further: the intermodes thresholding of Prewitt and Mendelsohn (1966) (e.g. implemented in ImageJ (Schneider et al., 2012)) shows a detection rate similar to that of Otsu’s method. It smoothes the histogram by running averages of size 3 until there are exactly two local maxima. The threshold is then the arithmetic mean of those. But only in the example of particle separation (subsection 6.2) the intermodes thresholding performs superior to Otsu’s method, because the intermodes thresholding chooses a higher threshold than Otsu’s method. The higher threshold is chosen in all our examples, and the reason for this is simple and shows also the drawback of intermodes thresholding: dominant peaks in the histogram — such as the peak at 1 in our perfectly “illuminated” artificial images — are detected as one maximum and shift the average towards the maximum brightness value. Nonetheless, the performance of such a simple approach is excellent.

In Figure 1, an image is shown representing the clustering by Otsu’s method. For all further steps of calculating the particle center only the threshold value detected by Otsu’s method is used, not the binarized image itself. This means that in the first step of the moment method the threshold is used to identify “white” pixels belonging to a particle, while in the second step the position is calculated with the brightness values of the original image.

Different algorithms are compared with respect to different signal to noise ratios in Figure 4:

alg01

moment method (subsection 1.3)

alg02

Preprocessing (Image Processing) with $I_{base}=I_{threshold}$

alg03

Preprocessing (Image Processing) preprocessed by a Hanning filter ($N=5$)

alg04

Preprocessing (Image Processing) with the threshold automatically detected by Otsu’s method

alg05

Preprocessing (Image Processing) with gamma adjustment with $\gamma=3$

alg06

Preprocessing (Image Processing) preprocessed by a Hanning filter ($N=5$) and fitted by a generalized Gaussian (section 5)

alg07

Preprocessing (Image Processing) with the threshold automatically detected by Otsu’s method

Examples of single noisy particle images are shown in Figure 5. For a $\operatorname{SNR}\,$ of $5$, not all of the $10000$ particles could be detected by any of the algorithms. While particles in the example images Figure 5 are easy to identify for human eyes, the algorithms are more sensitive to the noise.

Comparing Preprocessing (Image Processing) and Preprocessing (Image Processing) in Figure 4, we can see that the Hanning filter used in Preprocessing (Image Processing) leads to a better detection rate in the case of high noise.

Feng et al. (2007) recommend using $I_{base}=I_{threshold}$ in the moment method to reduce uncertainties in the found particle positions. Preprocessing (Image Processing) uses this method, and Figure 4 shows that indeed the error can be reduced in comparison with the pure moment method in Preprocessing (Image Processing).

However, this is not true for small particles ($\sigma_{\{x,y\}}\in\{0.1,0.5\}$), as shown in Figure 6. While Feng et al. (2007) explain, why an inappropriately chosen threshold leads to pixel locking, here we see that another reason for pixel locking can be missing information, such as particles consisting of not enough pixels, as seen in Figure 7. This is not an error of the algorithm, but of the measurement itself. Figure 8 illustrates the influence of the particle size on Preprocessing (Image Processing). For particle sizes $\sigma_{\{x,y\}}\in\{0.1,0.5\}$, the positions calculated by Preprocessing (Image Processing) are not statistically fluctuating around the real position. Instead, there is a systematic deviation depending on the real position – a similar behavior can be observed for all presented algorithms. For example, a particle with $\sigma_{\{x,y\}}=0.1$ consists of more than one pixel only, if the absolute value of the chosen sub-pixel coordinate is greater than $0.25$ (cf. Figure 7). Therefore, for a coordinate with an absolute value of the chosen sub-pixel coordinate of less than $0.25$, any algorithm can find just that one pixel, and consequently only detect the exact coordinate of it, which yields [math] as the sub-pixel coordinate.

The clustering by Otsu’s method used in Preprocessing (Image Processing) and Preprocessing (Image Processing) performs well. Only for very small particles, in the example given by Figure 6 and visualized in Figure 7, a stable detection is not possible. Increasing gamma (Preprocessing (Image Processing)) does slightly improve the accuracy of Preprocessing (Image Processing), but not all particles can be detected any more. Comparing Preprocessing (Image Processing) and Preprocessing (Image Processing) shows that Otsu’s method does not choose the best threshold. But as an automatic procedure processing all available pixel values, it can reduce human errors in the process of choosing the threshold.

4 Moment Method

(Blob Detection)

In Figure 9, different algorithms are compared with respect to different probabilities of pepper noise:

alg08

Preprocessing (Image Processing) with a search radius $r=1$ (moment method with particles being only single connected sets, similar to (Feng et al., 2007; Ivanov and Melzer, 2007))

alg09

Preprocessing (Image Processing) with a search radius $r>1$

alg10

Preprocessing (Image Processing) preprocessed by a Hanning filter ($N=5$)

We can see that for high pepper noise Moment Method (Blob Detection) is not able to detect all particles correctly — it finds too many, because some particles are split in two by the pepper noise.

Using Moment Method (Blob Detection), the generalized moment method described in subsection 1.3, we are able to detect all particles correctly. The same holds for the Hanning filter in Moment Method (Blob Detection). The quality of the latter is comparable to the generalized moment method. The only draw-back is the larger computing time of Moment Method (Blob Detection) (see Figure 10, comparison of used processor times of Preprocessing (Image Processing) and Preprocessing (Image Processing)).

5 Fitting (Postprocessing)

Given approximate coordinates from the blob detection of the last section 4, we can try to enhance them by fitting a generalized Gaussian point spread function, which is visualized in Figure 2 and given as:

[TABLE]

The fit is performed locally to every single particle. Therefore, we split the given image in non-overlapping squares with an approximated particle coordinate located in the center of the square. Every square is chosen with a maximal side length under the given restrictions.

Initially, the distance $d$ of two particles $i$ and $j$ is defined as:

[TABLE]

Then, for a given particle coordinate $p_{0}:=(x_{0},y_{0})^{T}$ the closest particle $p_{1}:=(x_{1},y_{1})^{T}$ is found as:

[TABLE]

With $\delta:=\frac{1}{2}d((x_{1},y_{1})^{T},(x_{0},y_{0})^{T})$ and $\left\lfloor{.}\right\rfloor$, $\left\lceil{.}\right\rceil$ the floor and ceiling functions333In practice, this should not be a mapping to integers $\mathbb{Z}$, but to image coordinates — a subset of non-negative integers $N_{0}$. we get the vertices of the square as:

[TABLE]

Now we generate separate problems for every square or particle. Let the given image be a matrix $I=\left(I_{i,j}\right)$. Here, we use the original image and not the prefiltered one.

An artificial image can be created with (1) as $A(x,\sigma,p)=\left(A_{i,j}(x,\sigma,p)\right)$ with the particle coordinate $x\in\mathbb{R}^{2}$ and $\sigma$, $p$ the parameters for the generalized Gaussian point spread function. With $b$ the averaged brightness of the background, this results in the optimization problem:

[TABLE]

For solving this optimization problem we use the algorithm L-BFGS-B (Byrd et al., 1995; Zhu et al., 1997; Morales and Nocedal, 2011), implemented in the python module/library SciPy (Jones et al., 2001–).

The gradient of the objective function is calculated numerically by a symmetric difference quotient if possible (e.g. on the boundary of the feasible solutions we cannot calculate a symmetric difference quotient).

In Figure 4, different algorithms were compared with respect to different signal to noise ratios, including those with fitting. Additionally, in Figure 10 the processor time used for the simulation is given. It is obvious that a small improvement by fitting a generalized Gaussian (Preprocessing (Image Processing) and Preprocessing (Image Processing)) leads to a large calculation time (Figure 10 shows a factor of $60$ to $144$).

The improved detection rate by Hanning filtering (e.g. Preprocessing (Image Processing)), and the automatically chosen threshold by Otsu’s method (e.g. Preprocessing (Image Processing)) each lead to a larger error, as demonstrated in Figure 4 and Figure 6. This can be corrected by fitting (e.g. Preprocessing (Image Processing) and Preprocessing (Image Processing)).

Instead of successively fitting every individual particle in one image, one can try to fit all particles in one image simultaneously. This assumption leads to a high dimensional optimization problem. In our implementation with the algorithm L-BFGS-B, this problem could not always be solved successfully. If it was successful, the result was sometimes slightly better than fitting every individual particle, but at the cost of a considerably increased computing time: it was about $100$ times higher for $100$ particles, and about $700$ times higher for $1000$ particles.

6 Examples

6.1 Velocity

In this section we will regard the velocity in images; this is the velocity of a particle in the image plane – e. g. a 2D mapping of a real 3D motion.

Let us assume a sequence of images with a temporal distance of $dt=10\mbox{ ms}$ between $2$ consecutive images (equivalent to a frame rate of $100$ images per seconds), with particles modeled as a Gaussian with $\sigma=1\mbox{ pixel}$ and a $\operatorname{SNR}\,=100$.

In our simulations (see Figure 4 and Figure 6), the presented algorithms Preprocessing (Image Processing) and Preprocessing (Image Processing) yield a root mean square error of about $0.014\mbox{ pixels}$ (or better). Assuming a distribution around [math], the root mean square is the standard deviation.

When calculating a particle velocity from $2$ consecutive particle positions, each subject to the same uncertainties, error propagation leads to an error of $(2\times 0.014)/dt$ $[\mbox{px}/\mbox{s}]$ in the velocity.

As an example, for a $4$ megapixel camera ($2048$x$2048$ $\mbox{px}^{2}$) with a field of view of $4\mbox{ cm}$ by $4\mbox{ cm}$ this leads to an uncertainty in the velocity of $0.056\mbox{ mm}/\mbox{s}$.

In the following table a few examples are given with the resolution in $\mbox{mm}/\mbox{px}$, $dt$ in ms and the velocity in $\mbox{mm}/\mbox{s}$:

[TABLE]

We neglect here that changing the resolution also changes the size of a particle on the image sensor. Otherwise, the error for the resolution of $0.005\mbox{ mm}/\mbox{px}$ would be reduced dramatically.

With the last example representing a $4$ megapixel camera ($2048$x$2048$ $\mbox{px}^{2}$) at 80 frames per seconds with a field of view of $11\mbox{ mm}$ by $11\mbox{ mm}$ we would be able to measure the velocity of particles (\varnothing=9.19\mbox{ \mum}, $\rho=1.51\mbox{ g/cm}^{3}$) at room temperature, which would be about $0.08\mbox{ mm/s}$. Experiments with a crystalline 2D complex plasma, and a comparable spatial camera resolution, were analyzed with the presented Preprocessing (Image Processing) by Knapek et al. (2007a), yielding reasonable kinetic energies of the particles.

Knowing the uncertainties, especially for particle velocity calculation, should not be underestimated: Gaussian noise can easily mask a Maxwellian velocity distribution, but it is not possible to separate the two distributions (see e.g. Knapek (2011, chapter 7); here, as well as in other publications (Knapek et al., 2007a, b), the applicability of our algorithm to real-world data is demonstrated in more detail). Therefore, it is of high importance to know the limit of resolvable particle motion (depending on particle size, $\operatorname{SNR}\,$ and the algorithms) for a specific experiment before interpreting the results.

6.2 Particle Separation

Here, we assume two particles which are nearby each other on the image plane (see Figure 11), e.g. due to their overlapping motion in different layers. Further, we assume both particles have the same size on the image, e.g. an uniform illumination, a good enough depth of field and same particle size and texture.

We neglect the Gaussian beam profile of the laser in the simulation. This profile might give additional information about the depth of a particle: the pixel intensity values would reflect the (ambiguous) position of the particle within the spatial extent of the laser beam, and could be used for a relative depth evaluation between particles, but we do not use this kind of information in the presented algorithms.

Now we can use all algorithms introduced in section 3. Since we want to separate both particles, we do not want to detect particle pictures like Figure 9 (right) as one particle. Also, we cannot use a too large Hanning amplitude filter, because it would wash-out distinctive edges of nearby particles. Therefore, we have to demand a good $\operatorname{SNR}\,$ to avoid the necessity of preprocessing. This is usually available in typical images of complex plasmas obtained with a laser filter which suppresses background illumination (e.g. from the plasma glow).

We now introduce additional algorithms similar to Preprocessing (Image Processing) and Preprocessing (Image Processing):

alg11

Preprocessing (Image Processing) with a large threshold of $190$ and a search radius $r=1$

alg12

Particle Separation with fitting a generalized Gaussian (section 5)

In Figure 12, different particle distances, as visualized in Figure 11, are compared. Here, for all presented algorithms a search radius of $r=1$ was used. Particles with distances of $3$ pixels or less could not be separated. The moment method (Particle Separation) described in subsection 1.3 with a large threshold of $130$ is able to separate particles with $\sigma_{\{x,y\}}=1$ with a distance of only $4$ pixels. Choosing the threshold automatically (for minimal particle distances of $3$, $4$, $5$ and $10$ pixels, the threshold was $78$, $71$, $73$ and $73$, respectively) with Preprocessing (Image Processing) separates these particles down to a distance of $5$ pixels. In both cases, postprocessing the images by fitting (Particle Separation and Preprocessing (Image Processing)) reduced the uncertainties.

The possibility to separate close-by particles can prove helpful in 3D diagnostics, e.g. for the analysis of data taken with a stereoscopic setup (several cameras viewing the same volume from different angles). Particles located close to each other on the image plane are typical features for this kind of diagnostics, and algorithms are needed to reliably detect particles in each of the camera views as the basis for a subsequent triangulation(Alpers et al., 2015).

7 Conclusion

In this paper we presented a comparison of several methods and algorithms for particle tracking from images. The methods and algorithms were tested on artificial images simulating data as they are obtained in complex plasma experiments, including realistic image noise (additive white Gaussian noise, salt and pepper noise). To increase the statistical significance, images with a large number of particles ($10000$) were analyzed. The proposed procedure for particle tracking consists of three major steps: image processing, blob detection and postprocessing.

In section 3, we show that using a Hanning filter to remove Gaussian noise during image processing results in a better detection rate in the presence of high noise, whereas the accuracy of the found positions is slightly reduced (Figure 4).

For images consisting of features (the particles) and a background (noise), the choice of a good threshold is important during image processing. With Otsu’s method (used in Preprocessing (Image Processing), Preprocessing (Image Processing)), we introduce this concept of automatic thresholding for particle detection in complex plasma for the first time (section 3). Other automatic thresholding techniques were tested, but did not prove to be suitable. The clustering by Otsu’s method performs very well (Figure 4), yielding almost the same results as the manually chosen threshold for all but the smallest particle sizes (Figure 6). On the one hand, choosing the right threshold value is not an easy task, and an automatic method can dramatically reduce human errors. On the other hand, an automatism prohibits using expert knowledge of the user in special circumstances, e.g. for the task of particle separation (subsection 6.2).

In subsection 1.3 we introduce an improved algorithm for blob detection: we generalized the set used for the moment method to a not necessarily simply connected set, and show that we can considerably improve particle detection in the presence of certain kinds of noise (e.g. salt and pepper noise, Figure 9) with this generalization.

We present a postprocessing method in section 5 to further enhance the accuracy of the detected particle positions by fitting a generalized Gaussian function to the intensity profiles of the particles. This is in particular interesting if prefiltering is necessary due to noisy images. Then, the postprocessing can reduce errors introduced by the prefilter (Figure 6, Figure 4). Also, it can increase the sub-pixel resolution of particle positions. This is especially interesting for applications where small particle velocities, e.g. thermal velocities, are calculated from the positions (subsection 6.1).

Another application is shown in subsection 6.2: Particles which are close-by each other on the image plane can be separated by either manual or automatic threshold detection, and position accuracy was improved by the above postprocessing method. This kind of situation typically appears in the individual camera images of a stereoscopic imaging system.

In summary, image processing with a Hanning filter (Preprocessing (Image Processing)), and a subsequent blob detection with the moment method detects in the most cases all particles in our simulations, but needs a manually chosen threshold. Automatic threshold detection (Preprocessing (Image Processing)) results in a slightly reduced accuracy and a reduced detection rate, but has the advantage of the automatism. In both cases, postprocessing the acquired positions by fitting (Preprocessing (Image Processing) and Preprocessing (Image Processing)) reduced uncertainties in the particle coordinates at the cost of a large calculation time (Figure 10 shows a factor of $60$ to $144$), but for specific experiments with the requirement of a good sub-pixel resolution this can be very useful and worth the effort.

8 Acknowledgments

The authors are funded by DLR/BMWi under grant FKZ 50WM1441, and by STMWi (Bayerisches Staatsministerium für Wirtschaft und Medien, Energie und Technologie).