1. Introduction

To gain information about the dynamics of a many particle system, which can in general be any system containing three or more particles, is of great interest in a wide range of research fields and especially in statistical physics. One way to obtain information about the dynamics is by studying the kinematics of the system.

An example of such a many particle system is a complex plasma, where a set of dust particles is introduced into a charge-neutral plasma [1]. Due to the repulsive interaction of the quickly negatively charged dust particles in a confinement, the particles organize themselves e.g. in crystalline or liquid structures. Therefore, the compelling property of this system is that it can be used as a model system for different states of matter. See [1] for a description of the transition of a plasma crystal from a solid into a liquid phase. A complex plasma is a system where the time evolution of the particles can be recorded by an optical device to estimate the kinematics [2, 3], independent of the state the complex plasma is in.

To study physical processes, such as for example melting transitions [1, 4] one has to track all particles for all time steps in order to obtain information about the kinematics. In this work we focus on the estimation of the particle positions, on which all particle tracking algorithms heavily rely on. The particle tracking methods themselves are not subject of this work.

The conditions under which optical images are recorded are often not ideal and can complicate the detection of a particle significantly. Bright field images may lack contrast, if the features in the image have a low light absorption and therefore appear invisible in the image [5]. Target detection in infrared images can be difficult due to background clutter and low contrast [6, 7]. Living cell imaging, where illumination intensities are reduced to prevent photobleaching and photodamage, is another example of an imaging process with poor contrast [8, 9].

Some of the prominent particle detection methods involve techniques like local background subtraction [10], smoothing kernels [11, 12], the combination of both [13, 14], morphological image filtering [6, 7, 15] and wavelet based filtering [16, 17]. An overview of different spot detection methods can be found in [18].

A quite different approach is the technique of template matching for detecting structures in optical images. The technique of template matching can be considered as a part of image registration methods. Image registration is the task of aligning images from the same scene but taken for example from different viewpoints or by different sensors. However [19] found that besides the traditional application of image alignment template matching can also be used to detect particles in optical images. Structures are detected by exploiting the statistical similarity of a subregion in an image and an equally sized template containing the structure one is searching for. An overview over different image registration techniques, including template matching, is given in [20, 21].

A crucial part of the template matching algorithm is the choice of the similarity measure, which quantifies the similarity between two images. As [19] have shown the correlation coefficient (CC), can be used to detect particles in the presence of strong signal independent noise. However, one drawback of the CC is that it only measures linear dependencies between two variables. In contrast the mutual information (MI), as a nonlinear measure, does not suffer from this limitation.

Both, the MI and the CC are well known similarity measures in a wide range of fields. While the development of the CC can be traced back to 1895 (see [22, 23] for reviews and interpretation) the foundations of MI were developed by Shannon in 1949 [24]. The use of MI for the purpose of image registration was developed both by Collignon [25] and Viola [26]. MI is applied for example in [27] for detecting relationships between variables and in [28] for image registration. For a survey of the use of MI especially for medical image processing please refer to [29].

While [19] show that the CC is very effective in detecting particles in the presence of high Gaussian white noise, the assumption of such a noise distribution is not always valid. In fact there exist several imaging scenarios where the noise is not distributed according to a zero mean Gaussian distribution with a constant standard deviation, but is a function of the signal itself. Some examples, where the noise in optical images is signal dependent are speckle noise in Synthetic Aperture Radar [30] and Ultrasonic Scanner [31] imaging systems. The noise in x-ray CT images is distributed according to a signal dependent Gaussian distribution [32]. The noise in magnetic resonance images follows a riccian distribution [33]. In the case of a CMOS sensor the intensity dependent photon noise contributes most to the overall noise variance [34]. A signal dependent noise model of a CCD sensor is presented in [35]. We therefore replace in this paper the noise model of [19] by a signal dependent noise model according to [36].

To access the kinematics of a many particle system, both the positions and the velocities of the particles have to be estimated. The common way to achieve this is by tracking the particles in consecutive frames and calculating the velocity simply by the quotient of the distance traveled and the time elapsed [37, 38]. This procedure involves the problem of identifying the same particle in consecutive images [10]. This identification problem, commonly referred to as single particle tracking, is present in many scientific fields and there exist therefore many methods to achieve this task. Examples to tackle the single particle tracking problem are the nearest neighbor approach used in [13, 14], the combination of the nearest neighbor method with an additional Kalman filter [15] and the seemingly most preferred method called multiple hypothesis tracking [17, 39]. A comparison of performance about the mentioned methods and other popular particle tracking algorithms can be found in [40].

In the end, all the algorithms have to deal with the problem of the particle identification in consecutive images. If one could estimate the velocity of a particle out of a single image one could overcome the problem related to the identification process or at least supply the tracking algorithm with additional information to improve its results.

At this point our second extension of the template matching technique of [19], which lies in the template itself, comes into play. The imaged shape of a particle can strongly depend on its velocity. Isotropic particles with a high velocity appear elongated along the direction of motion due to movement across several pixels during exposure time.

While this elongation could lead to serious complications in other particle detection algorithms, the template matching algorithm is not affected by the change of particle shapes. In fact we can take advantage of the information about the elongation of a particle to estimate its position, velocity and direction of motion simultaneously. To achieve this simultaneous estimation we use a new set of anisotropic templates, with which we are able to measure the distance a particle traveled during exposure time.

The problem of tracking particles in different frames to obtain their velocity can then be replaced by estimating the velocity of a particle out of a single image. The information about the velocities for all particles in a sequence of images, without any particle tracking, can reveal important physical properties of the system. For example estimating the velocities of the particles for all time steps out of a series of images enables one to estimate important parameters such as the viscosity in a complex plasma [41, 42].

2. Methods

In this section we introduce the methods we use to detect particles and their velocities in optical images via template matching. We start by introducing the templates and how to construct them. Next we describe the two similarity measures we use and how to calculate them for two gray value images. We show how to use the similarity measures to detect the particles in a gray value image. Finally, we introduce the method how to estimate the position, velocity and direction of motion of a single particle simultaneously by template matching.

2.1. Isotropic template profile

In the case of an isotropic particle shape we model the intensity profile T_x,y of the template to detect it according to [19] with a Gaussian intensity profile. We calculate the values of the Gaussian profile inside a window of size 2 · w_x + 1 in x and 2 · w_y + 1 in y direction, where the parameters w_x and w_y are integer valued. We always take the Gaussian profile to be centered at the center of the window. The number of rows and columns of the window is taken to be odd, in order to be able to define a center of the rectangle within pixel accuracy. The intensity values T_x,y for the template with standard deviations of σ_x and σ_y and maximum value T_max can then be calculated by

2.2. Anisotropic template profile

To generate a template mimicking a particle with high velocity one has to take a closer look at the underlying imaging process. Instead of being imaged at a constant position, the particle travels across several pixels during exposure time. This movement creates an elongated particle shape which we will call a trace in the following. To create a template with a similar shape as such a trace we simulate the particle movement during exposure time by adding up isotropic templates, shifted by equidistant steps Δx and Δy in x and y direction with respect to the preceding template, until we reach the required length of the trace. So for a particle traveling N · Δx equidistant steps in x and N · Δy equidistant steps in y direction we construct the corresponding template T_elong by following formula

For every combination of the two parameters v_x and ϕ we can create a unique template. Assuming a velocity range for v_x from 1 to 30 $\frac{pixels}{frame}$ and 5 rotation angle values for ϕ we have to build a bank of 150 different templates to test for every combination of v_x and ϕ. The resulting velocity in x-direction v_res,x and in y-direction v_res,y of a particle detected by a template with parameters v_x = N · Δx and v_y = N · tan(ϕ) · Δx can then be derived with the known quantities of the pixel size µ and the exposure time τ

In Fig. 1 we show examples of templates for a stationary and a fast moving particle.

 

Fig. 1 Examples of templates both with a standard deviation of 5 pixels: stationary particle (left), fast moving particle with v_x = 30 and ϕ = 0.2 * π(right). Colour coding indicates pixel brightness.

Download Full Size | PDF

2.3. Mutual information

The mutual information is in general defined for two random variables, which we will call A and B. For these random variables A and B, their marginal probability distributions p_A(a) and p_B(b) and their joint probability distribution p_AB(a, b) are first calculated. $\widehat{MI}$ (A, B) is then defined by [29]

Two random variables are statistically independent, if the joint probability distribution p_AB(a, b) can be rewritten as p_A(a) · p_B(b). In this case the mutual information becomes zero. Hence the mutual information is measuring the distance between the joint distribution p_AB(a, b) and the distribution in the case of independence p_A(a) · p_B(b).

If we think of the elements of two equally sized gray value images A and B as random variables we can calculate the $\widehat{MI}$ for these values. We approximate the marginal probability distributions p_A(a) and p_B(b) and the joint probability density p_AB(a, b) by binning gray values of the two images and calculate their histograms normalized by the number of elements N in the images. If using $\widehat{MI}$, as a similarity measure, one has to be aware of the fact that the result depends on the size of the image. Therefore, we use a normalized form of the $\widehat{MI}$, denoted by MI(A, B), which takes a maximum value of one independent of the image size:

Here H(A) and H(B) denote the marginal entropies of the random variables A and B [43].

2.4. Pearsons correlation coefficient

Again thinking of the elements of two gray value images A and B as realizations of two random variables, the correlation coefficient ρ of the two gray value images of equal size N is defined by

2.5. Particle detection algorithm

Detecting particles in a gray value image by template matching can be divided into two main steps:

The first is to calculate the similarity measure under consideration for the whole gray value image. We obtain the similarity measure for every pixel in the image by calculating the similarity value of a sliding window around it with the template and save the resulting value in a new matrix. We call this new matrix the similarity measure map in the following. The sliding window has to have the same size as the template. So every pixel in the similarity measure map corresponds to the similarity value between a sub-region in the image around the pixel and the template. The edges of the image are not taken into account and particles there will not be detected.

Particle detection is then performed by finding the local maxima in the similarity measure map. Here, we adapt the method from [13]. We search for the largest maximum in the similarity measure map, set a region around this local maximum to zero and search again for the next local maximum. This iteration is repeated until the new local maximum falls below a threshold we choose. If not stated otherwise the size of the region we set to zero is determined by the window size of the template we use to detect the particle. For example for a local maximum of the similarity measure map calculated with a template with a window size of 21 × 21 we set a region of 21 × 21 pixels around this maximum to zero. These local maxima correspond to the coordinates with the highest similarity to the template and therefore to the most probable particle centers according to the template matching algorithm.

In Fig. 2 we illustrate the whole process. For this illustration we used a synthetic image with 100 particles simulated by a isotropic Gaussian intensity distribution with σ_x = σ_y = 5. We superimposed signal independent noise and a signal to noise ratio of ten (see Sec. 3.1).

 

Fig. 2 Illustration of the particle detection algorithm: left: original noisy image with template in the inset; middle: similarity measure map calculated by MI; right: original noisy image with estimated centres plotted as green dots.

Download Full Size | PDF

2.6. Position, velocity and direction of motion estimation

In the case of estimating the position, speed and direction of motion of a particle out of an image, the dimensionality of the problem increases. In the case of the pure position estimation one has to calculate the similarity measure map once for every pixel in the image. Now one has to calculate the similarity measure map for each template with a different velocity and direction of motion for every pixel in the image. So for a image of R rows and C columns the resulting similarity measure map has the dimension of R×C×K×M for K different velocities and M different directions of motion.

To find the maximum values according to all three parameters we first determine the maximum similarity measure map value with respect to the parameters v_x and ϕ for every position in the image, independent of whether it is the estimated center of a particle or not. We thereby reduce the dimension of the similarity measure map to R×C and call it the reduced similarity measure map in the following. Afterwards we calculate the positions of the particles out of the reduced similarity measure map in the same way as before for the pure position estimation.

We demonstrate the whole process in Fig. 3. We generate a noisy synthetic image containing a single particle with $vx=30\frac{pixels}{frame}$ and ϕ = 0.1 · π and estimate the position velocity and direction of motion of this particle simultaneously.

 

Fig. 3 Illustration of the position, velocity and direction of motion estimation algorithm: Top row: reduced similarity measure maps with estimated particle position plotted as a green dot; center: original noisy image containing one noisy particle with estimated particle position plotted as green dot (both MI and CC estimate the position of the particle at the same location); lower row: similarity measure values for different velocities and angles of motion at the estimated position of the particle (maximum value is outlined by a black frame)

Download Full Size | PDF

From Fig. 3 we can see that both the MI and the CC estimate the position, the velocity v_x and the direction of motion measured by ϕ precisely. We further observe for different values of v_x the differences in similarity are far less pronounced than for different value for ϕ for both similarity measures. In addition we observe that the MI decays much faster for diverging values of ϕ from the true value for ϕ. This indicates are more accurate estimation of the true values for the MI.

3. Validation of the algorithm

To test the particle detection algorithms we use synthetic images, where we know the position, velocity and direction of motion for every particle in the image.

We apply our particle detection algorithm to the synthetic images under different noise conditions and different particle densities to measure the performance of both similarity measures. As a last part in this section we show a typical application of the velocity estimation procedure where we estimate the averaged velocity profile out of images simulating laminar flow.

3.1. Synthetic image generation

First we choose random positions for the particle centers, taking into account a minimal distance of the particles. Then at each random position we add a “particle” modeled by a Gaussian intensity profile to the image I_orig. The center of each particle corresponds to the maximum value I_max of I_orig. After adding the particles to the image we corrupt the image with white additive Gaussian noise denoted by η multiplied by Σ characterizes the strength of the noise.

Then we normalize the noisy image to the range from 0 to I_max, where I_max is the maximum value of the uncorrupted image I_orig. The strength of the noise Σ is defined by the signal to noise ratio (SNR), which is a free parameter in the image generation process.

For the performance evaluation of the similarity measures we compare the estimated particle locations to the true particle locations, which we know for synthetic images. We will consider a particle as truly detected with accuracy δ_c if the difference of the estimated position and the true position δ is smaller or equal δ_c. We define δ by the euclidean distance

Here x_true and y_true denote the known real x an y coordinate of the particle in the synthetic image and x_est and y_est the x and y coordinates estimated by the particle detection algorithm.

3.2. Noise model

As we mentioned in the introduction, we are considering both signal independent and signal dependent image noise. To model the noise in the synthetic images we adapt the linear noise model of [36], which is a function of the noiseless image I_orig:

The term $\sqrt{I_{orig}}\cdot {\xi }_{dep}$ is the signal dependent part of the noise model. It models the contribution of the so-called shot noise. The term ξ_ind models signal independent background noise in the image. To test the dependence of the detection results on the different noise components we introduce two additional parameters α and β, both either 0 or 1. The expression for $\widehat{\eta }$ then reads:

With α and β we can control the contributions of the signal dependent and independent noise components. In 3.1 we multiplied η by Σ and identified Σ with the strength of the noise. This is only valid if the noise distribution of η has a standard deviation of one. This is not the case in general if α = 1. We therefore normalize $\widehat{\eta }$ as a last step by its standard deviation $\sigma (\widehat{\eta })$:

3.3. Particle detection in the presence of noise

We compare MI and CC for different noise strengths and noise compositions. The performance of each similarity measure is assessed by the number of true positive detections compared with the number of false positive detections. In addition we compare the detection accuracy in terms of the mean detection error and the standard deviation of the detection error as a function of the SNR for the different noise compositions.

3.3.1. Signal independent noise (α = 0, β = 1)

Here we compare the number of true positive and false positive detections for different SNR values in the signal independent noise regime. We create the synthetic images for the performance evaluation with the following parameters. 100 isotropic particles modelled by templates with σ_x = σ_y = 5 are distributed randomly in the images. The noise strengths vary from a SNR of 10 to a SNR of 0.5. We calculate the performance of the similarity measures for an accuracy δ_c of 2 pixels.

The results for the MI and the CC are shown in Fig. 4.

 

Fig. 4 Number of true positive detections (tp) in the upper row and number of false positive detections (fp) in the lower row for signal independent noise at a detection accuracy of δ_c = 2. Both quantities are plotted against the threshold θ normalized by the maximum value of the respective similarity measure maps θ_max. Color coding indicated in the lower right plot is identical for all plots.

Download Full Size | PDF

In the signal independent noise regime the MI and the CC perform very similar until a SNR value of 4. Both similarity measures detect all of the particles correctly while in both cases the number of false positive detections is 0. For a SNR of 3 the performance of the MI is slightly lower, but still nearly all particles are correctly detected. At a SNR values of 2 and 1 the performance of the MI starts to break down. The number of true positive detections decreases to slightly above 80 while the number of false positive detections is still practically 0 for a SNR of 2. But for a SNR of 1 also the number of false positive increases significantly and is of the same order as the number of true positive detections (note the different scale of the tp and fp plot). The CC detects all 100 particles correctly with 0 false detections for a SNR of 2 and about 80 particles for a SNR of 1 while detecting only 14 particles by mistake. For a SNR of 0.5 the number of tp detections is again comparable for the MI and the CC and is in both cases very poor. While the MI detects a large amount of particles wrongly for a SNR of 0.5 the number of false detections are still comparably low even at these high noise levels for the CC. Overall one can say that in the signal independent regime the CC performs better for low SNR values as the MI.

We further calculate the mean detection error $\overline{\delta }$ as a function of the SNR value for signal independent noise in Fig. 5. To calculate $\overline{\delta }$ we position a single particle in an image and corrupt the image with signal independent noise. Then we estimate the position of the particle and calculate the detection error by Eq. (10). This procedure is repeated independently 100 times. Then $\overline{\delta }$ corresponds to the mean value of the 100 values for δ. We indicate the spread of the detection errors with vertical error bars calculated by the standard deviation the distribution of δ. It can happen that a similarity measure “sees” only noise and detects a noisy subregion of the image by chance as a particle. Therefore we disregard detections with a distance larger as half of the size of the particle from the true particle positions in the calculation of $\overline{\delta }$. We separately plot the number of detections purely due to noise, which we denote by ε.

 

Fig. 5 Mean detection accuracy $\overline{\delta }$ and number of noise-induced detections ε for signal independent noise.

Download Full Size | PDF

In Fig. 5 we plot $\overline{\delta }$ and ε both as a function of SNR value for pure signal independent noise.

Here we observe a similar performance of the MI and CC as in Fig. 4. The error in the position estimation is in general smaller in the case of the CC as compared to the MI. The MI also detects many particles wrongly for a SNR value of 0.5.

3.3.2. Signal independent and signal dependent noise (α = 1, β = 1)

We perform again the same evaluation as in the pure signal independent case for the MI and CC. The results are shown in Fig. 6.

 

Fig. 6 Same as Fig. 4 but for signal dependent noise.

Download Full Size | PDF

Here the MI is little affected by decreasing SNR values in contrast to the signal independent noise regime. For this reason we choose a very low minimum threshold to make the differences especially for the number of false positive detections visible. The MI can detect nearly all of the particles correctly until θ = 0.4 · θ_max while the number of false positive detections decreases very rapidly and are all very close to 0 for θ ≥ 0.2 · θ_max. The CC shows a very similar performance to the signal independent regime until a SNR value of 4. But in contrast to the signal independent noise regime the performance of the CC starts to break down at a SNR value of 3. For a SNR value of 2 and below the MI clearly outperforms the CC both in a larger amount of true positive detections and with nearly zero false negative detections.

In this noise regime we show that MI is able to detect nearly all of the particles correctly practically independently of the SNR if the minimum detection accuracy of two pixels is sufficient. In general, we can say that CC is not suited for particle detection for signal dependent noise at low SNR values. As for the signal independent noise we plot $\overline{\delta }$ and ε for signal dependent noise as a function of the SNR in Fig. 7.

 

Fig. 7 Same as Fig. 5 but for signal dependent noise.

Download Full Size | PDF

Here we observe again the superior performance of the MI as compared to the CC. In case of the MI $\overline{\delta }$ is nearly constant for all the SNR values. For the CC $\overline{\delta }$ is about three times larger than for the MI at a SNR of 0.5. For a SNR of 4 and above $\overline{\delta }$ is very similar for both the MI and the CC, which is in good agreement with the results we obtained in Fig. 6.

3.4. Detecting particles with low interparticle distances

In case of dense particle clouds, detecting the particles can be difficult due to overlapping image shapes of nearby particles. We therefore investigate the dependence of the particle detection performance of the MI and CC in the case of high density particle clouds using synthetic images containing five particles. The particles in the synthetic images are all modeled by the same isotropic Gaussian intensity profile with σ_x = σ_y = 5. Four particles are arranged in the corners of a square. The fifth particle is positioned in the middle of the square. The size of the square is then gradually decreased while the size of the five particles is left constant. The minimum distance of two particles in one square is then the distance from the middle particle to one of the outer particles. Although the impact of the noise is not in the focus of this investigation, we add some noise with a SNR of 20 to the images to create a more realistic environment.

In Figs. 8 and 9 we show the response of the MI and the CC to the five particles with two different minimal interparticle distances of 18 and 12 pixels.

 

Fig. 8 Upper row: noisy image in the first, horizontal slice through the image in the second and diagonal slice in the third column, middle row: response of the MI to the five particles in the first, horizontal in the second and diagonal slice through MI in the third column, lower row: response of the CC to the five particles in the first, horizontal in the second and diagonal slice through CC in the third column. Slices through the noiseless image are plotted in blue together with the slices through the noisy images and the similarity measures, both in green, for better visualization. The minimum peak to peak distance is 18 pixels.

Download Full Size | PDF

 

Fig. 9 Same as Fig. 8 but with a minimum peak to peak distance of 12 pixels.

Download Full Size | PDF

One can observe that in the case of the CC the peak corresponding to the center particle is significantly below the peaks corresponding to the outer particle centers in the diagonal slice. In contrast the MI exhibits three sharp and distinct peaks in the diagonal slice.

All five particles, which can still be separated by the MI, can no longer be separated any more by the CC.

We further analyze the performance of the MI and the CC similar to Sec. 3.3 by calculating the number of true positive and false positive detections for a larger set of close-by randomly distributed particles. We focus on the influence of the particle density on the particle detection performance by MI and CC. To investigate this we use synthetic images containing isotropic Gaussian shaped particles with σ_x = σ_y = 5. We vary the particle density by using different minimum interparticle distances Δ. We distribute the particles in the image by iteratively choosing a random position for a particle and check that there is no other particle within the chosen minimum interparticle distance Δ. If there is already a particle within Δ we try again and choose a different random location. This process is repeated until we cannot find a new position for a particle which has a minimum distance to all other existing particles within one million iterations. Because the particles are distributed randomly in the image the number of particles for a given value of Δ is also a random number. In contrast to Sec. 3.3 we corrupt the image with a moderate amount of signal independent noise with a SNR of 20 and therefor impose a higher detection accuracy δ_c of 1 pixel. We only use signal independent noise, because at this noise level the type of noise has no impact on the performance of the similarity measures as seen in Sec. 3.3. In Fig. 10 we show the performance of the MI and the CC for two different values of Δ = 20,15.

 

Fig. 10 Number of true positive (tp) and false positive (fp) detections of the MI and the CC for Δ = 20 with 385 particles (first row) and Δ = 15 with 693 particles (second row). The true particle centers are depicted by red dots in the noisy image in the first column.

Download Full Size | PDF

We observe that the CC detects many particles wrongly despite the low noise while the MI performs significantly better for Δ = 20. If we decrease Δ to 15 the performance of the CC breaks down. The number of false positive detections exceeds considerably the number of true detections at low thresholds. MI shows a much better performance for Δ = 15 in very good agreement with the results we presented in Figs. 8 and 9, although one has to note that a non negligible amount of misclassification occur, too.

The capability of the MI to separate particles even in the case of a low interparticle distance becomes crucial if one tries to locate particles in a dense three dimensional cloud. Observed with one camera, all of the particle shapes are projected on to the imaging plane of the camera. Particles which are originally behind each other along the axis perpendicular to the imaging plane can reduce the interparticle distance substantially and even overlap in the camera image. In such a scenario one has to be aware of the fact that particles which are out of the imaging plane of the camera are imaged differently than particles on the imaging plane. As we will also see in Sec. 4 they appear much darker in the resulting 2D image. While this could cause a problem because the noise in the image has a much larger impact on these dark particles, it also opens a new opportunity to access information about the z component of the position of a particle perpendicular to the imaging plane of the camera. By imposing different template shapes the particle cloud could be segmented into several layers with respect to their distance to the camera. As long as the imaged shape of the out of plane particle is not fully covered by noise the particle will get detected by our method, as it detects particles by structural similarity independent of the brightness.

3.5. Velocity profile estimation

One of the typical application of our velocity estimation method is the computation of the velocity profile of a fluid. In order to validate our algorithm we model such a fluid motion by distributing particles in a sequence of images, with velocities according to a laminar pipe flow velocity profile. Our goal is to calculate the known velocity profile out of the sequence of images. For a detailed discussion on the velocity profiles please refer to [44].

In the laminar case the velocity can be described as a function of the distance r from the center of the pipe

Here v_max denotes the maximum velocity of the flow and R the radius of the pipe. To simulate the laminar flow we create 30 synthetic images each containing 100 particles modeled by Gaussian intensity profiles with σ_x = σ_y = 3 and velocities according to the noisy laminar velocity profile described by Eq. (15). We choose a maximum velocity of 20 $\frac{pixels}{frame}$. The center of the horizontal flow is set to the image center. To simulate the fluctuation around the time averaged velocity profile we add noise to the velocities of the particles in the single images according to a zero mean Gaussian distribution with a standard deviation of one. The final velocity v_final of a particle at position r is then simulated by

Here we set negative values of v_final due to the noise to zero because in a laminar flow the velocity of the flow has to always point in the same direction for all values of r. As a last step the image is degraded by signal independent image noise with a SNR of 50. We choose this high SNR value because faster particles get much more affected by the noise than slow particles which appear much brighter in the image. Also we did not vary the direction of motion for the particles, because in laminar pipe flow the movement is parallel to the long axis of the pipe. To estimate the original velocity profile we first estimate the velocities of the individual particles in every frame separately using 26 templates ranging from v_x = 0 to v_x = 25. Then we calculate the mean velocity of the particles $\overline{v}$ as a function of the displacement r to the center of the flow using the previously estimated velocities of the single frames. As a last step we fit the calculated mean velocities by Eq. (15) with free parameters v_max and R. The results are shown in Fig. 11.

 

Fig. 11 Velocity profile estimation by MI (top row) and CC (bottom row) for signal independent image noise: on the left we show an example image to demonstrate the position and velocity estimation by the MI and the CC in a single frame, on the right we show the estimated velocities of the particles in green, the estimated velocity profile of the flow in blue and the true velocity profile in red. In the left images we can see that all particles are detected by both the MI and the CC (coloured dots) and that the estimated velocities (indicated by colour coding of the dots) matches the elongation of the particles very well. The estimated velocity profile by either MI or CC (blue crosses) matches the true velocity profile (red line) precisely.

Download Full Size | PDF

We can observe that the calculated velocities scatter around the original velocity profile shown in red for both similarity measures. This results in a nearly perfect estimation of the velocity profile for both similarity measures and is an excellent proof of concept for our velocity estimation method applied to synthetic data.

4. Application to experimental data

In this section we apply our methods to experimental data from the PK-4 experiment, which is currently operated on the Columbus module on the International Space Station (ISS). The PK-4 experiment is designed for the investigation of a complex plasma in a dc discharge. For a detailed description of the experimental setup please refer to [45–47]. The particles of the complex plasma are recorded by the scattering light of the microparticles illuminated by a laser. Besides the illumination laser the experimental setup also involves a manipulation laser which makes it possible to accelerate regions of the particle cloud. The laser therefore creates a velocity gradient of the particles perpendicular to the axis along the laser. An example of such a laser driven velocity gradient is shown in Fig. 12.

 

Fig. 12 Example image of a dust particle cloud accelerated by a six watt laser. The laser is centered at the vertical center of the image.

Download Full Size | PDF

The manipulation laser can be operated at different laser powers hence inducing different particle velocities. Furthermore it is possible to modulate the laser power over time. We use this feature to test our algorithms in Sec. 4.2 for experimental data.

4.1. Position estimation of dust particles

It is naturally less straightforward to compare results for position estimation for experimental data than in the case of synthetic data, where the ground truth about position and size of the particles is known. However, we can compare results of both similarity measures by eye. Particles which are at the edges of the laser sheet appear much darker in the image and are almost invisible for the eye. Nevertheless, these particles are detected by the similarity measures. To make them visible we show the images on a logarithmic scale, reducing the dynamic range of the brightness values in the image. We can therefore judge more easily the performance of the different similarity measures also for the particles at the edges of the laser sheet.

In Fig. 13 and 14 we show the particle detection performance for both similarity measures. Both similarity measure maps are calculated with an isotropic template with standard deviation of one pixel and a window size of 11 pixels. The threshold for the maximum detection is 30 percent of the maximum value in the respective similarity measure maps and the region set to zero around the local maxima is of size 5 × 5 pixels. All parameters were chosen by hand to maximize performance.

 

Fig. 13 Position estimation of dust particles: original images each with zoomed in section with positions estimated by CC (left) and MI (right) both plotted as red dots in the top row and similarity measure maps to estimate the positions for CC (left) and MI (right) on the lower row. The CC only detects one particle in the left particle pair in the zoomed in sections whereas the MI detects both.

Download Full Size | PDF

 

Fig. 14 Same as in Fig. 13. Again, the CC only detects one of the particles in the middle particle pair whereas the MI detects both.

Download Full Size | PDF

In the middle of the zoomed in sections of Fig. 13 we observe two pairs of two nearby particles each. While both particles of the pair on the right are detected correctly by both similarity measures, the pair on the left is only detected correctly by the MI. If we look at the corresponding zoomed in sections of the similarity measure maps we can observe that for the left pair the response of the CC is a highly bright spot for the left particle and rather dark for the right one. In the case of the MI we observe two distinct spots with less difference in brightness.

In Fig. 14 we show a similar example and find that again the MI detects two particles where the the CC only detects one particle. In the case of the MI, the response to the two particles is practically of equal brightness.

To summarize, in the case of the CC one can observe that if two particles are close together, one of the spots in the similarity map can be significantly darker which can cause the maxima detection algorithm to disregard it because it falls below the chosen threshold. In contrast the spots of two close particles in the similarity map calculated by MI are more similar in brightness. The maximum detection algorithm is therefore much more likely to detect both particles. The results of the examples shown give clear evidence that the MI outperforms the CC in case of dense particle clouds and are in good agreement with the results for synthetic data we presented in Sec. 3.4.

4.2. Velocity wave analysis

Here we analyze experimental data where a part of the particle cloud is accelerated by the manipulation laser with a modulated laser power of (6 ± 2) W. The modulation frequencies of the laser intensity, known to be 3 and 10 Hz, ideally create a modulation of the particle velocities with the same frequency. This setup gives us the perfect opportunity to test our velocity estimation algorithm on real world data where part of the ground truth, here the modulation frequencies, is known. In Fig. 15 we show an example of the accelerated region of the particle cloud for maximum and minimum laser power. Only sections of the whole images are shown for demonstration purposes.

 

Fig. 15 Example images illustrating the effect of the laser power modulation of (6 ± 2)W and 3Hz. For minimum laser power (upper row) the elongation of the particles is much smaller than for the maximum laser power (lower row). The difference between the two images is 5 frames, which corresponds to approximately 0.143 s.

Download Full Size | PDF

To determine the velocity modulation frequency we calculate the mean velocity $\overline{v}$ of the particles in a part of the laser accelerated region for each of 300 consecutive frames. By this we obtain the time evolution of $\overline{v}$ in the accelerated region of the particle cloud. To estimate the velocities of the single particles in each frame we use templates of σ_x = σ_y = 1. We choose a velocity range for v_x from 0 to 30 pixels in steps of 1 pixel and the values for ϕ from −0.2 π to 0.2 · π in steps of 0.1 · π. Thus, the images are processed with a bank of 150 anisotropic templates and one isotropic template. From the time evolution of $\overline{v}$ we calculate its power spectrum, which should have a distinct peak at the modulation frequency of the laser. To calculate the power spectrum of a discrete signal x_t one has to first calculate the discrete Fourier transform of the signal X_f. The discrete Fourier transform can be calculated from the signal values in the time domain x_t of length N by following formula:

The power spectrum is then given by the square of the absolute value of every value of every value of X_f. In Figs. 16 and 17 we show the time evolution of the mean velocity and its power spectrum. We fit a linear function of the form v(t) = a · t + b to the time series to test for temporal trends. We begin the fit at the 20th frame to safely disregard the frames of the initial acceleration of the particles. From the power spectrum we can determine the frequency $f_{max}=\arg \max \left({\left|X_f\right|}^2\right)$. Due to the finite frame rate and length of the image sequence the frequency resolution of the power spectrum is limited to a finite frequency step width of Δ_f.

 

Fig. 16 Frequency estimation of velocity wave driven by laser of 6 ± 2W laser with modulation frequency of 3 Hz: time evolution of the mean velocity with fitting parameters in the insets (left), power spectrum of the time evolution of the mean velocity with frequency $f_{max}=\arg \max \left({\left|X_f\right|}^2\right)$ and frequency step Δf in the inset (right). Results for the MI are shown in the upper row, for the CC in the lower row.

Download Full Size | PDF

 

Fig. 17 same as Fig. 16 but with a laser modulation frequency of 10 Hz.

Download Full Size | PDF

We observe a practically constant fit of the time evolution of $\overline{v}$, ensuring that the system is in a steady state. The maximum value of the power spectrum appears at f_max = 2.9538 Hz for both similarity measures.

In the case of the modulation frequency of 10 Hz the slope of the fits m slightly increases for both similarity measures. Still they are both on the order of 10⁻³ which justifies the assumption of a steady state. We calculate f_max = 9.9554 Hz for both MI and CC.

For both modulation frequencies, we show that the power spectra of the mean particle velocities show a distinct peak matching the respective modulation frequencies of the laser. Taking into account the uncertainty in frequency due to the finite signal length these results represent an excellent proof of concept for our velocity estimation method applied to real world data.

5. Conclusions

In this work we extended the method of particle detection by template matching introduced by [19] by the nonlinear similarity measure MI and a new template design. The newly developed anisotropic templates enable us to estimate not only the position of a particle, but also its velocity and direction of motion out of a single image. Thereby, we highly increase the information available in a many particle systems, whose time evolution is recorded by an optical device. We show the difference in performance of the MI and the CC for particle detection depending on the type of noise applied to synthetic images. In the signal independent noise regime the CC shows better performance for strong noise. However, for signal dependent noise the MI outperforms the CC by far for low SNR values. We further demonstrate that the MI can separate nearby particles to a much greater extent than the CC and show how this affects the ability to accurately detect particles in dense particle clouds. Again we demonstrate the superior performance of the MI compared to the CC.

As a typical application of our velocity estimation algorithm we calculate the velocity profile of a laminar pipe flow, by estimating the particle velocities out of a sequence of synthetic images. The very accurate reconstruction of the velocity profile represents an excellent proof of concept for synthetic data.

Applying our methods to experimental data of the dusty plasma microgravity experiment (PK-4) on the International Space Station (ISS) we found strong evidence for the superior performance of the MI for dense particle clouds. Furthermore, we show that we can correctly estimate the frequencies of velocity waves in the experimental data induced by a laser with modulated power.

Both for synthetic images as for the experimental data in this work we investigate the position and velocity estimation of isotropic shaped particles. But particle detection by template matching can also be applied for anisotropic shaped particles. Given an adequate template the particle position estimation can be done the same way as for isotropic particles. If the anisotropy in the image of the particles is caused by both the original shape and the velocity of the particles one can in principle also infer the velocity of the particles. The procedure will be more complicated: More calibration measurements and a larger set of (less trivial) templates will be required to disentangle the shape and velocity information of the moving anisotropic particles.