© E. Shahamatnia et al., Published by EDP Sciences 2016

1. Introduction

To develop an automatic solar feature tracking system, a hybrid approach combining an image processing technique and an evolutionary computation algorithm is used. The evolutionary computation algorithm, Particle Swarm Optimization (PSO), is a search algorithm for finding the global optimum, inspired by the social behavior of bird flocking (Kennedy & Eberhart 1995). The original idea of PSO is to simulate how birds interact with each other and with their environment in search for food. That makes PSO a population-based evolutionary computation algorithm consisting of relatively simple individuals, called particles, where each particle represents a potential solution to the optimization problem at hand. PSO tries to find the optimum by iteratively improving the candidate solutions in regard to a given measure of quality. In the search for optimal solutions, particles keep track of their best values and communicate the best solutions, per iteration, back to the population. Instead of exhaustively searching the whole search space, PSO evaluates only a small portion of the search space. The particle swarm searches more efficiently by paying attention to the more promising areas of the search space. The search is conducted in a combination of stochastic manner with some guiding mechanisms. Since the algorithm uses primitive mathematical operations it is computationally inexpensive to implement. The simplicity, flexibility, and good performance of PSO have made it a popular choice as a global problem solver in a wide range of real-world applications such as human tremor analysis (Eberhart & Hu 1999), tracking dynamic systems (Eberhart & Shi 2001), RNA molecule structure prediction (Agrawal & Agrawal 2015), and synthesis of antenna arrays (Ram et al. 2014). In many applications where PSO has been used, it has shown consistently good performance (Hu et al. 2004; Poli 2008; Yang 2015). Moreover, thanks to its speed, simplicity, and flexibility in formulating problems, PSO has been successfully used in many hybrid algorithms to solve specific problems such as antenna optimization, classification of biological data, and vehicle routing (Robinson et al. 2002; Holden & Freitas 2005; Marinakis et al. 2010).

PSO is a member of the family of stochastic optimization methods, along with genetic algorithms and simulated annealing, two other popular algorithms in this field. That type of optimization algorithm does not guarantee finding an optimal solution, but these algorithms are good at finding near optimum solutions, by searching very large spaces and making few or no assumptions about the problem. Although the choice of algorithm is usually affected by the problem specifications, several studies in the literature show that PSO outperforms genetic algorithms and simulated annealing in many test problems investigated (Habib & Al-kazemi 2005; Hassan et al. 2005; Panda & Padhy 2008; Yang et al. 2008; Ethni et al. 2009). Further, it is reported (Hassan et al. 2005) that the convergence rate of PSO is steadier than genetic algorithms and that its computational efficiency and execution time are better because of the smaller number of function evaluations needed for convergence. Unlike classic optimization methods, such as gradient descent and quasi-Newton methods, PSO does not use the gradient, therefore it can also be used for optimization problems that are not differentiable, are partially irregular, noisy, or change over time.

A specialized image processing technique is included in the hybrid PSO-Snake algorithm. Image segmentation is a frequently addressed problem in digital image processing and deformable contours are well-known algorithms for object boundary delineation. Introduced by Kass et al. (1988), the active contour model is a class of deformable contours for finding features of interest in an image, by formalizing it as an optimization problem. Starting by a rough approximation of the object boundary, contours evolve to find the precise boundary of that object (Ballerini 1999; Ballerini & Bocchi 2003). Evolving the contour takes into account the low-level image data such as image gradient as well as high-level image information such as contour continuity, shape, texture, color, etc. The active contour model is also known as the “snake model” and the contours in this model are called snakes. Due to their flexibility, snakes are widely used in several applications such as image segmentation, shape modeling, stereo matching, and object tracking (Ballerini & Bocchi 2003; Karlsson et al. 2003; Niu 2006; Wildenauer et al. 2006). For every snake there is an energy function associated with it. Specifically, the snake model is an optimization algorithm that works as an energy minimization procedure where the snake with the lowest energy is considered the best match for representing the object. The snake’s evolution process starts by setting up an initial curve (snake) and moving it toward the target object. The snake movement is governed both by internal forces, within the curve, and external forces from the image. The former maintains the snake shape while the latter steers the snake toward the image feature, defined by the external forces. Snake energy is calculated as a combination of those two forces. In the standard snake model, this energy is minimized in a numerical way by iteratively solving a pair of Euler equations. This is usually a computationally expensive operation.

Other problems associated with the snake model are sensitivity in placing the initial snake and the poor convergence at concave object boundaries (Davatzikos & Prince 1994; Bresson et al. 2007). Many researchers have tried to ameliorate this problem by improving the capture range of image forces (Cohen & Cohen 1993; Leroy et al. 1996; Prince 1997; Park et al. 2001; B. Li & Acton 2007). Several works (Amini et al. 1988; Mun et al. 2004; Bresson et al. 2007) have tried to address the snake model limitations in regard to noise tolerance, local minima sensitivity, and stability. Very few have succeeded in proposing an alternative that would solve the snake model’s limitations without compromising its performance, flexibility, and simplicity.

One of the successful variations of the snake model is the geometric active contour model, particularly popular for medical image analysis (McInerney & Terzopoulos 1996). Geometric snakes, sometimes called geodesic snakes (Caselles et al. 1997; Paragios & Deriche 2000; Xu et al. 2000; He et al. 2008), are mainly implemented based on the level-set theory to capture boundaries through a continuous curve, and thus are able to adapt to topology changes while evolving contours. With their intrinsic capability to manage splitting and merging contours, they can be used to overcome the sensitivity to the initial snake configuration. In its standard implementation, the level-set function of all pixels in the image is updated at each iteration, which makes the geometric active contour model computationally intensive. There have been suggestions for lowering the computational complexity of this method, but issues related to noise and ill-defined boundaries (both present in the solar images) remain problematic. Thus, geometric snakes are not further discussed in the scope of the current work.

In the standard snake model, the contour is defined by a set of finite points called control points. For evolving the initial contour toward the final object boundaries, snake energy is minimized, i.e. snake control points are iteratively updated by solving a pair of Euler equations. Other methods to minimize the snake energy have been suggested such as: dynamic programming (Amini et al. 1988), greedy algorithms (Lam & Yan 1994), genetic algorithms (Ballerini 1999; Ballerini & Bocchi 2003; Mun et al. 2004), and swarm-based optimization algorithms (Asl & Seyedin 2006; Zeng & Zhou 2008; Nebti & Meshoul 2009; R. Li et al. 2009; Tseng et al. 2009; Shahamatnia & Ebadzadeh 2011). In recent years, several researchers have used PSO to optimize the snake energy, mostly by constraining the search space. Some works have used multi-population PSO (R. Li et al. 2009; Tseng et al. 2009) where every control point is confined to a sub-swarm spatially distinct from other sub-swarms. Nebti & Meshoul (2009) restrict the search space of each contour control point using polar coordinates. Zeng & Zhou (2008) use an iterative method to rank the best set of particles’ position at each epoch, preventing particles from intersecting.

Most of those methods are used as a general optimization technique in solving the snake model equations without modifying the snake structure. They formulate the snake energy calculations as a minimization problem according to each method’s specifications. In this paper, we take the hybrid PSO-Snake approach introduced in Shahamatnia & Ebadzadeh (2011) and explore its versatility by further extending it to solve a real-world problem from the astrophysics domain. The hybrid PSO-Snake model is no longer a general problem solver, but it is a specialized image processing technique for searching the image space. The method presented here customizes the previous PSO algorithm (Shahamatnia & Ebadzadeh 2011) to overcome some snake model drawbacks including snake initialization, concave boundaries, sensitivity to noise, and local minima. The performance of the PSO-Snake algorithm regarding those issues was measured using synthetic images, specifically generated to isolate each issue, and the results are discussed in detail by Shahamatnia & Ebadzadeh (2011).

In the PSO-Snake model the simple structure of the PSO is preserved yielding an algorithm with a low order of complexity and hence good processing time. These factors are of utmost importance for precisely calculating the differential rotation of solar features. Determining the exact nature of the differential rotation of both the solar surface and the solar interior is still one of the most serious open issues of solar physics (Thompson et al. 1996; Scherrer et al. 2011; Hanasoge et al. 2015). The solar surface rotates differentially, lower latitudes rotate faster than higher latitudes (Howard 1984). However, the differential rotation mechanism, most likely caused by the interactions between the convection and the overall rotation, is not exactly known. Differential rotation plays an important role in solar activity generation – at least the large-scale manifestations of solar activity are related to changes in the local magnetic field which may have their roots in variations in the differential rotation. Rotational irregularities may also serve as indicators of hypothetical processes beneath the solar surface. One example could be the location of a layer where rotational speed changes abruptly (a jet stream). Sometimes it is called a layer of torsional oscillation (Howard & LaBonte 1980). The current paper focuses on reporting the results of applying the proposed tool on coronal bright points (CBPs). In this work we trace the location of bright points in a series of selected images obtained by the Atmospheric Imaging Assembly (AIA; Lemen et al. 2012), an instrument on board the Solar Dynamics Observatory (SDO; Pesnell et al. 2012).

Coronal bright points (CBPs) or simply bright points are small and bright structures observed in the extreme ultraviolet (EUV) and in the X-ray frequencies of the solar spectrum (Habbal & Withbroe 1981; Brajša et al. 2001). They are known to have a mean lifetime of about 8 h, a typical maximum area of 2 × 10⁸ km², but still they look like a small dynamic loop structure on the solar images (Brown et al. 2001) as shown in Figure 1. Bright points are associated with small, bipolar magnetic features in the photosphere. A large quantity of bright points (several thousand) emerges over the surface of the Sun each day.

Fig. 1. A sample solar image from the 5th October 2010 with several CBPs marked. (Image courtesy of NASA/SDO).

Tracking the coronal bright points with high precision, over extended periods of time, will help solar physicists and space weather scientists to better understand this important solar feature. Such automatic tools will allow solar researchers to precisely process large amounts of solar data and hence improve solar physics models. Hence, the main aim of this paper is to assess the results of applying a hybrid PSO-Snake algorithm to tracking of coronal bright points. Due to the dynamic nature of PSO-Snake hybrid algorithm, detected contours are flexible and can conform to changes in shape and size of the deformable objects like CBPs. The tracking result is then used for calculating the differential rotation of coronal bright points. The result of PSO-Snake hybrid algorithm is cross-referenced and compared with a state-of-the-art study, which entails a manual procedure done by an expert (Lorenc et al. 2012).

Small CBPs were chosen for tracking because they were well detectable and suitable for the PSO-Snake (or any other) tracking algorithm. Compared to other solar features such as sunspots, CBPs have simpler shapes and smaller structures. CBPs are particularly suitable for the PSO-Snake algorithm, because of their small size and their shape that make them visually distinguishable. This facilitates the convergence of the initial snake to the boundaries of the CBP in the initial phase of the PSO-Snake algorithm. In this study, for calculating the velocity, the CBP is represented by a single reference point. In the case of CBPs this point can be considered as the center of mass of the structure, while for complex shapes such as some sunspots such assumption is not intuitive. Another characteristic of the CBPs that makes them suitable for the PSO-Snake algorithm is that they are usually spatially confined, contrary to some other solar features such as filaments. The PSO-Snake algorithm is a contour-based tracking scheme and while it is well capable of adapting to deformities, changing shapes, and recovering from partial edge information, it cannot overcome too many faint borders, completely disconnected object segments, and fine details of the complex shapes.

CBPs are also very good tracers since they extend to much higher latitudes than sunspots (Sudar et al. 2015). They are also one of the solar corona’s most ubiquitous features (McIntosh & Gurman 2005), quite numerous in all phases of the solar cycle while, for example, sunspots are often absent in the minimum of the cycle (Sudar et al. 2015). McIntosh & Gurman (2005) developed a method of automatically detecting EUV bright points and applied it to the archive of EIT data from the launch of SOHO in 2005. As a result, they produced a database of all detected bright points that can be used to extract numerous diagnostics of the solar corona over the 23rd solar cycle. CBPs are also suitable tracers for the determination of the solar differential rotation, because they are localized objects which are very well distributed over the solar disk (Brajša et al. 2014). The rest of this paper is organized as follows: Snake model, PSO, and PSO-Snake algorithms are reviewed in Section 2. Section 3 provides the experimental results and discussions. Finally, conclusions are provided in Section 4.

2. PSO-Snake hybrid algorithms

The first version of the hybrid algorithm is a merger of the snake model and PSO (Shahamatnia & Ebadzadeh 2011). It integrates the snake model’s contour evolving paradigms with particle dynamics from PSO. After that the region of interest is approximated and the snake model will be able to find the precise boundary of that object.

As mentioned before, in the snake model, the contour of snake has an energy associated with it that is related to the location of the snake on the image (external force) and with the geometrical characteristics of the snake (internal force). The idea is to minimize the integral measure that represents the snake’s total energy, by evolving the snake over time. The original snake model achieves this minimization by iteratively solving a pair of Euler equations on the discrete grid, resulting in a computationally expensive algorithm (Karlsson et al. 2003). Two main approaches for snake representation are geometric and parametric representations. Geometric models use an implicit representation based on the curve evolution theory and are usually implemented with level-set techniques. Effectively handling multiple objects and topology alteration is the advantage of this approach, with the cost of being computationally more complex. On the other hand, the parametric approach is computationally efficient and easy to interact with users (Horng et al. 2010). In the parametric implementations, the snake is defined as a curve p(s) = (x(s), y(s)), having arc length s. As it is shown by Eq. (1), a number of discrete points, called control points or snaxels, characterize the snake at each given time, t (Kass et al. 1988). The PSO-Snake hybrid uses this representation since it matches the snaxels to the PSO particles well. The total snake energy is the sum of its internal (spatial) and external (geometrical) integrals as shown in Eq. (2).(1) (2)

The snake model is considered to be a controlled continuity spline under the influence of internal and external forces that induce the snake energy. Internal energy, shown in Eq. (3), consists of two terms that are the first and the second derivatives of the snake with respect to s. The first term coerces the spline to act like a membrane while the second term makes the snake to act like a thin plate (Kass et al. 1988). The external energy determines the snake’s relationship to the image. It is formulated in a way that its local minima correspond to the image features of interest. Various external energies can be employed such as image intensity, image gradient, object size or shape. One common definition used for gray-level images, I(x, y), is the gradient of Gaussian, as shown in Eq. (4), where γ is the tuning coefficient and ∇ is the operator to calculate the gradient of the resulting image. G_σ (x, y) is a two-dimensional Gaussian function with the standard deviation σ which is an important parameter of the model. Depending on the image noise, complexity of the object, and thickness of the boundaries, a proper value for σ can be chosen to control the blurring effect of the Gaussian function.(3) (4)

Equations (3) and (4) represent the internal and external energies of the snake and the total snake energy is the sum of these two energies. Solving the classical snake model means minimizing the total snake energy, or in other words the sum of the integral of the two equations. In its original form, minimizing the energy function of the total snake energy gives rise to two independent Euler equations which are explained in detail in the original reference (Kass et al. 1988). The classical snake model optimization approach, which uses the Euler equations for minimizing the snake energy, is a complex method because it requires calculating the higher-order derivatives of the contour.

The leading part of the PSO-Snake hybrid algorithm is its PSO component. PSO is a population-based evolutionary optimization algorithm. Population in the PSO is called a swarm and consists of a number of particles; each particle is a potential solution to the optimization problem depending on its location in the search space. Each particle, i, has a position, x_i and a speed, v_i, which are initialized with random values. Over a set of iterations, each particle’s position on the search space is updated by revising its velocity according to its best experience, y, and its neighbors’ experiences, . Particle position and its corresponding fitness value are stored as personal best experience and form the cognitive aspect of particle evolution. Another aspect of the particle position update is called the particle social behavior, which shows particles’ influence on their neighbors. The PSO neighborhood can be defined with various topologies such as ring, star, von Neumann, and random. If all particles of the population are in the same neighborhood it is called global best (gbest) PSO, while if the neighborhood is restricted to a subset of the swarm it is called local best (lbest) PSO (used in this study). Figure 2 shows how the position of each PSO particle in the search space is updated under the influence of its velocity vectors. As it is shown in Figure 2, updating the particle velocity is also affected by the inertia velocity, i.e. the current velocity of the particle. Inertia velocity is important for controlling the particle velocity and preventing particle movement from radical changes. The following equations display the dynamics of the canonical PSO algorithm for updating particle velocity and position:(5) (6)where x_i(t) and v_i(t) are the position and the velocity of the ith particle at time t, y_i(t) and denote the best positions discovered by the ith particle and its neighborhood up to the time t. The inertia weight, ω(t), controls the previous velocity impact. Usually the inertia weight is decreased dynamically during the run time, to balance between exploration in the early iterations and exploitation in the later iterations. Coefficients r₁ and r₂ are random numbers. Weights of cognitive and social aspects of the algorithm are represented by acceleration factors c₁ and c₂, respectively. As shown by Van den Bergh (2002) regulated values for inertia and acceleration weights can be used to achieve guaranteed convergence.

Fig. 2. Search mechanism of PSO particle in a two-dimensional space. (a) Shows particle x at time t under influence of its inertia, cognitive and social velocity components, (b) new position of the particle at time t + 1 is the algebraic sum of the velocities acting upon it.

The PSO-Snake hybrid algorithm integrates the snake model mechanisms with PSO dynamics. While most swarm intelligence approaches in the literature that are used in conjunction with the snake model try to optimize the snake model equations, the PSO-Snake hybrid does not employ PSO algorithm only as a general problem solver to optimize snake energy minimization, but it also customizes the standard PSO to better solve this specific type of image processing problem. Early experiments on medical image segmentation (Shahamatnia & Ebadzadeh 2011) and sunspot tracking (Shahamatnia et al. 2012) reported promising results. It should be noted that running the standard snake model with the balloon force and the improved snake model with the Gradient Vector Flow forces (Xu & Prince 1998) on the CBP input images used in this paper did not yield any usable results. Due to the high noise in the image and weak CBP boundaries, the standard snake model initialized around a CBP could not converge to the CBP and evolved arbitrarily according to the sensitivity set by the balloon force. Testing the enhanced snake model improved the capture range of the image force, but still for the case of input images used in this paper it failed to be of practical use, since it requires meticulous parameter tuning for each individual CBP. A detailed discussion of the advantages of the PSO-Snake model over the traditional snake model is presented in (Shahamatnia & Ebadzadeh 2011).

In the Hybrid PSO-Snake model we use a population of particles where each particle is a snaxel of the contour. All particles together form the contour and hence the population is the final solution. As the algorithm runs, each particle updates its position and its velocity according to its personal best experience, local best experience, and according to the internal force of the snake and external force of the image. This gives the PSO-Snake dynamics a wider range of informative guides to update the particle position so that it converges to the region of interest.

In relation to the PSO algorithm (and the PSO-Snake for that matter), the term convergence refers to a stable condition where particles have found an optimum position in the search space (which might be a local or a global optimum). In general, the initial particle state is not at equilibrium, and convergence of the PSO is defined as the eventual settling of particles at the equilibrium (Trelea 2003). In PSO-Snake the equilibrium is reached when the snaxels latch to the boundaries of the object of interest, i.e. the CBPs within the initial contour.

PSO-Snake hybrid explores the search space according to the PSO trajectory disciplines. This eliminates the need to have a separate searching window around each particle as many swarm-based snake optimization algorithms do (Nebti & Meshoul 2009; Tseng et al. 2009; Horng et al. 2010). These methods consider a searching window around each particle and evaluate every position inside that window to determine the snaxels’ next position. Since this local search is performed for each particle per iteration, it is a computationally expensive operation that is avoided in the PSO-Snake hybrid model. In contrast, the PSO trajectory disciplines do not require performing an exhaustive search for each particle move and hence alleviate the need for a local search window. According to the PSO kinematics, particles move toward the locations that they think will be promising by considering their own and their neighbors’ previous experiences. The PSO-Snake adopts the same concept and pushes snaxels to move toward promising locations according to some guiding terms. This is implemented via the PSO-Snake’s velocity update equation:(7)where pbest_i(t) and lbest_i(t) denote the best values experienced up to time t by particle i (personal best) and by its neighbors (local best), respectively. In the standard PSO, the fitness or the best experience of each particle is a location in the search space visited by the particle, for which the value of a specified objective function becomes minimum. In the PSO-Snake algorithm every particle constitutes a portion of the solution, hence the fitness of each particle should be calculated in regard to other particles. The fitness of the population corresponds to the snake energy and is calculated by Eq. (2). In the PSO-Snake algorithm the personal best experience of each particle is defined as the velocity of the particle in the best-experienced position (i.e. the objective function with the minimum value). Moreover, to update the pbest the current fitness is not compared with previously calculated fitness for the particle’s best experienced position, but the fitness is recalculated for that position in the new swarm. is the average of positions at time step t, approximating the center of mass of particles. This term pushes the snake to contract or expand with respect to the sign of its weighting factor, r₃, speeds up the algorithm, and is particularly useful when the snake is stagnated and there is no other compelling force. f.Image_i is the normalized image force corresponding to the external energy from snake model principles. For particle i, f.Image_i gives the image force at the position specified by that particle. Image force can be any arbitrary function depending on the application, but generally external energies such as the image gradient and the gradient of a Gaussian functional are enough for satisfactory performance. It must also be noted that the image force does not vary by time and it is calculated only once for an image pixel. c₄ is the weighting factor to control the effect of image force. Inertia weight, ω, is taken to be a relatively small constant and r₁, r₂, and r₃ denote random numbers. In PSO-Snake, as well as traditional PSO algorithm, velocity values are influenced by a uniformly distributed random weight. That’s because by its very nature, this algorithm is a stochastic search technique. Introducing some randomness in the search process can help to improve the algorithm speed in finding good results (Hoos & Stützle 2004) and several approximate global optimization algorithms such as simulated annealing, genetic algorithms, evolutionary strategies, PSO, and the bee algorithm employ mechanisms to inject randomness in their search process. Having a touch of randomness along with other forces guiding the particles helps the PSO-Snake algorithm to be more flexible in escaping from local minima and handling noise. Coefficients c₁, c₂, and c₃ are determined dynamically as a negative logarithmic function of f.Image_i and coefficient c₄ is set as a constant coefficient. This ensures that if there is high image force, i.e. direct image information exists for guiding the particle evolution path, the image force will have a higher impact on the velocity update. Otherwise if the particle is out of the image force capture range, coefficients c₁, c₂, and c₃ will take higher values and the particle evolution will be guided by other velocity update terms. The PSO-Snake hybrid algorithm uses lbest with ring structure and a radius of 3. The whole process can be summarized in the following steps:

Step 1. Initialization. If the input images need to be preprocessed it is done in this step. The algorithm requires all images in a tracking sequence to have the same resolution and image size. The RMS contrast measurement method explained in Peli (1990) can be used to check the consistency of image contrast in successive frames. Moreover for the coordinate logging system to work correctly, the orientation of the solar images should be corrected. The SDO-AIA images used in the current study are already preprocessed and do not require any further corrections.

Step 2. Initial Contour. The region of interest is chosen by the operator. This is the initial snake. For most cases a rough estimation of the initial contour is enough. This step is done only once when the coronal bright point appears.

Step 3. Internal parameters setup. The weight parameters for the PSO-Snake hybrid algorithm are initialized.

Step 4. Snake force calculation. The external force (image force) is calculated, once for every image. Depending on the application, noise, and complexity of the image, different image forces can be employed. In this work we use the gradient of the Gaussian function presented in Eq. (4) for calculating the image force. Figure 3 shows an illustration of sample image force calculated as a negative image.

Fig. 3. Illustration of sample image force as a basis for calculating external snake energy.

Step 5. Calculation of social and cognitive parts. We update the pbest value (the best velocity the snaxel ever experienced) and the lbest value as the average of velocities of neighboring particles.

Step 6. Moving snaxels. Each snaxel’s velocity is evaluated and then each snaxel’s velocity and position are updated.

Step 7. Snake detection. The convergence of the snake contour to the coronal bright point outline is checked, i.e. choosing the snake with the lowest total energy calculated. If the results are not satisfactory, the algorithm goes back to step 5. The outcome of this step is the CBP contour for an image frame.

Step 8. Tracking CBPs. The same CBP in the next image is tracked by feeding the subsequent image frame to the system as input. The algorithm loops back to step 4 and passes the specifications of the detected CBP.

Step 9. Stopping tracking. Tracking a CBP stops when it reaches the solar limb and disappears into the other side of the Sun, or when the CBP shrinks to a size smaller than a predefined threshold, according to the size and resolution of image.

3. Results and discussion

Our benchmark data are coronal images at 9.4 nm. This line is emitted by the Fe XVIII ion with a formation temperature of 7 × 10⁶ K, that can dominate the emission in CBPs. This allows us to detect and track relatively bright structures. Images containing the CBPs with appropriate features such as proper cadence are selected from JPEG images taken between 14 September 2010 and 20 October 2010, downloaded from the SDO archive (http://sdo.gsfc.nasa.gov/data/aiahmi/). Altogether we have observed the motion of 24 CBP structures with small and simple shapes so that the data gathered for the evolution profile of the CBPs will be more reliable. Furthermore, it helps to avoid complex CBP evolution such as splitting CBPs. The CBPs used for tracking were chosen by an expert (coauthor) and the dataset used in this study is a subset of the dataset used by Lorenc et al. (2012). Tracking the 24 CBPs over their lifetime resulted in 1577 measurements. In a previous study (Lorenc et al. 2012), through a manual procedure, the CBP structures were observed directly on a PC monitor in an interactive session. Figure 4 shows latitudinal dependence of sidereal angular speed of coronal rotation obtained by the manual procedure in comparison with previous works. Further details can be found in Lorenc et al. (2012). In that paper, the position of each of the 4998 CBPs is determined by an expert operator manually on JPEG images and converted to heliographic coordinates. We employ the same method described in Lorenc et al. (2012) to compute the latitude b from the pixel coordinates of the tracked CBPs on the input images. In brief, in that method first the Cartesian coordinates of each CBP in solar image are converted to the spherical coordinate, knowing the radius of the Sun and center of the solar disk in the image. Then the spherical coordinate can be converted to the heliographic coordinate by adjusting to the tilt angle of the solar axis for the observation date of the CBP. This value is looked up from an ephemeris maintained by the BASS2000 archive (http://bass2000.obspm.fr/ephem.php).

Fig. 4. Comparison of the derived values of the rotational speed. Points in the chart represent the measurements for the center of point-like structures with error bars showing the 95% confidence level intervals. The dotted curve shows the fit to the mean ω(b) values as a function of latitude b calculated by Eq. (8). Overplotted are the results of Howard & Harvey (1970) in solid line and Hara (2009) and Brajša et al. (2004) both in the dashed-dotted curve because they are almost identical. Image courtesy of Lorenc et al. (2012).

In order to choose which CBP to track we use the initial location of the CBPs from the previous work. Initial snakes are defined as a circle encompassing the CBP. Then, we run our PSO-Snake hybrid algorithm to track those CBPs and calculate measurements. Input images are converted to grayscale color map with 256 gray levels per pixel and image force is calculated by a gradient of Gaussian functional with σ = 3 pixels. The value for σ is chosen as an empirical parameter. Input images have 1024 × 1024 resolution. For test purposes we have chosen the same CBPs for which we have the benchmark data available from the expert’s manual CBP positioning. It should be noted that in the automated process, after choosing the CBP to be tracked (only once), the tracking process is automatic during the life span of that CBP.

Figure 5 shows a screenshot of the PSO-Snake hybrid algorithm tracking tool for a test image. The red circle shows the initial snake roughly locating the region of interest containing a CBP. Figure 6 shows how the initial snake is evolved under the PSO-Snake algorithm and the CBP boundary is detected. After the CBP is identified, its characteristics including the heliographic coordinates of its center of mass are calculated and stored. Then the next frame in the sequence is fed into the system. The detected CBP contour from the previous frame is used as a baseline to automatically track the CBP in the new frame. The temporal resolution of two successive images is calculated based on the image time stamp, then the new CBP initialization contour is started in the expected CBP location. Figure 7 shows a closer look at a tracked CBP along 37 frames. The results show the close matching between the automatic tracked positions and the manual CBP marking by an expert.

Fig. 5. A screenshot of the CBP Tracking system showing the Sun on 2010-10-05. (a) First snake (red circle) is initiated around the CBP we want to track, (b) zoomed-in view of the initial snake shows that the initial snake is a rough estimation of the location of the CBP. The algorithm will evolve the snake to find the exact CBP contour.

Fig. 6. Detection process of a selected CBP (2010-10-05). Initial contour around the region of interest is evolved to precisely delineate the CBP boundaries. The outer contour is the initial snake and the inner contours are the transitional contours of every 10 iterations of the PSO-Snake detection algorithm. In the final stage contours converge and do not change much through iterations.

Fig. 7. (a) Red contour represents the initial snake on first image (the Sun at 2010-10-05), (b) and (c) cropped view of the tracking process of the selected CBP during time, showing the CBP in frame 1 and frame 37, respectively. The cyan contour is the boundary of the tracked CBP, the red square is the expert’s manual CBP positioning result, and the yellow circle is the PSO-Snake hybrid algorithm’s automated tracking result for CBP’s center of mass.

To compare the precision of the algorithm, we calculate several parameters proposed in Lorenc et al. (2012). In that paper, after an expert manually determined positions of CBPs on the solar images, the following measurements were calculated (reported in Table 1 of the referenced paper): angular rotation velocity denoted by ω and measurement error at 95% confidence level denoted by Δω. The Sun’s surface rotational speed is calculated by Eq. (8) (Howard & Harvey 1970), where ω is the sidereal angular velocity of the rotation of the solar corona at latitude b and A, B, and C are the solar differential rotation parameters. Equation (8) is a common representation of the solar differential rotation (see references of Table 3). ω is expressed in [deg/day], b is expressed in [deg] and coefficients A, B, and C are in [deg/day]. Here we try to determine coefficients A, B, and C in order to best approximate the fitting of the heliographic latitude to the solar angular rotational velocity.(8)

Figure 8 illustrates the difference between our calculated values and the benchmark values for the test CBPs. It should be noted that part of this deviation (10%–15%) is due to code implementation differences, which in precise calculations impose minute variations, called calculation error. Tables 1 and 2 show the results obtained with manual CBP tracking and results obtained by the PSO-Snake hybrid algorithm for some structures. In these tables, structure is the tracking ID given to each CBP as it can be seen in Figure 1, n is the number of frames the structure is tracked in, and b is the detected heliographic latitude of the structure. ω_E is the orbital angular rotation velocity of the Earth, which can be looked up from solar almanacs. For the input images used in this study the almanac data is available only with a coarse temporal resolution, i.e. for the year 2010 the data is available for every 5 days. The exact ω_E for each CBP observation time is calculated by an interpolating technique.

Fig. 8. Deviation of the obtained results for b, ω, Δω on test data from the manual CBP marking of Lorenc et al. (2012). Every point in the horizontal axis represents a test CBP structure and its position in the vertical axis represents its error from the manual marking.

In addition, the PSO-Snake results were assessed by an expert. Table 2 and Figure 9 show that the obtained results are very close to the result of manual CBP tracking. Computed angular rotation velocity is within ±0.2 [deg/day] of the benchmark data most of the time. It also worth mentioning that in several cases, where results were showing bigger differences, with further investigation we found out that PSO-Snake hybrid algorithm behaves consistently and the user-error was the main cause of error. Table 3 compares the values obtained for the solar differential rotation constants A, B, and C of Eq. (8) in this method with results previously reported in the literature. As shown in Figure 9 left column, these coefficients are calculated by using a non-linear least squares curve-fitting function, with Levenberg-Marquardt algorithm (LMA) and least absolute residuals (LAR) estimation method and with no constraint applied on the coefficients. As part of the curve-fitting algorithm the confidence bounds on coefficients are set to be 95%. Figure 9 right column presents a visual comparison of Table 3 and illustrates how the calculated values for these coefficients determine the sidereal rotational velocity throughout all the solar latitudes. Early in our study, we had access to part of the manual database, and we did run our algorithm on 650 measurements. Later on, we got access to more data (with manual CBP markings for the comparison) and we ran the algorithm with more measurements (1577 measurements including the previous dataset). In Table 3, we decided to present the algorithm results with both the number of CBP measurements to assess the scalability and stability of the method. It should be noted that the number of measurements from the comparison is not included in the table, because it is very different in terms of implementation, method used for the coefficient estimation, number of measurements used, representing date and solar cycle, dataset size, type of input images and channels, filtering methods, etc. This table only aims to compare the final results.

Fig. 9. Left: The result of Levenberg-Marquardt fitting algorithm (LMA) with least absolute residuals (LAR), on different number of tracked bright point structures. Right: Comparison of the curves based on Eq. (8) from different works. The black line represents this work, and it should be noticed that the similarity between the obtained results from this method with different numbers of tracked structures ((a) 650 CBP measurements and (b) 1577 CBP measurements) indicates the stability of this method.

4. Conclusions and future works

In this paper the PSO-Snake hybrid algorithm has been used to address a real-world problem from the solar physics domain. The PSO-Snake hybrid model inherits a particle navigation scheme from PSO and modifies it to include snake model algorithm concepts. Combining PSO dynamics with snake model kinematics enables us to successfully overcome active contour difficulties, while preserving the simplicity of PSO. A detailed discussion of the advantages of PSO-Snake model over traditional snake model is presented in Shahamatnia & Ebadzadeh (2011). Comparing Eq. (5), the standard PSO velocity update equation, with Eq. (7) of the PSO-Snake algorithm, shows that two new terms have been added to the PSO-Snake algorithm, namely the that refers to the average of positions at time t and (f.Image_i), which is the normalized image force corresponding to external energy from snake model principles. By adding two new terms to the PSO velocity update equations, the PSO-Snake model still can evolve even if some of the components are missing or misleading.

Other aspects of the PSO were also changed in the PSO-Snake hybrid model. For example, the population in the PSO-Snake represents the positions of snaxels and hence, the whole population forms one viable solution to the problem. The PSO-Snake model is similar to the usual snake formulations in the definition of snaxels and external forces. However, unlike the usual snake formulation, which optimizes the snake by solving the Euler-Lagrange equations or represents it as a level-set formulation, PSO-Snake inherits PSO kinematics by moving the snaxels in the search space toward the global optimum (i.e. the final contour). This is fundamentally different from the snake’s elastic boundary constraints because the PSO particle evolving mechanisms are based on particles’ experience in the search space, while the usual snake formulations emphasize regulating the shape and physical characteristics of the contour. Elastic boundary constraints can be set to restrain the movement of a snaxel in a way that snaxels would demonstrate a desired physical characteristic. In contrast, PSO particles move in the search space irrespective of such constraints, i.e. they evaluate the position that they are visiting (underlying pixel) and tend to move to the promising areas by sharing information with their neighbors. In other words, while the elastic boundary constraint emphasizes the shape of the contour, the PSO-Snake emphasizes the search direction.

This algorithm uses velocity vectors to guide particles in the promising direction based on the information obtained from the search space that the particle itself or its neighbors have visited before. This is all implemented with the notion of consolidated velocity vectors, including the inertia of particles. PSO-Snake model can be used for different applications in image processing for object detection, image segmentation and tracking. Since the particle/snaxels have embedded velocity information, which can adapt to the object movement over time, this method is able to naturally incorporate and exploit that information.

By tracking the CBPs over time, the angular rotational velocity of the Sun is calculated. Lorenc et al. (2012) used a manual method to track and mark the location of each individual CBP in a dataset, however, the method is laborious and with a large number of images becomes unworkable for practical reasons. In this work we devised an adaptation of the PSO-Snake to automate the manual method of Lorenc et al. (2012) in tracking CBPs. Further, the human-defined dataset of Lorenc et al. (2012) was used as the ground truth for validating our algorithm. The obtained CBP tracking results were used to estimate the sidereal rotational angular velocity of the rotation of the solar corona. Further, we compared our results for the angular rotational velocity of the Sun with those reported in the literature. The results obtained in this study by running the PSO-Snake algorithm to track CBPs, on a sample dataset, showed good conformity with previous works. Despite the limited number of test cases used in this study, the proposed algorithm proved to be stable and gives consistent results. Moreover, we expect the parameter errors to reduce by running the algorithm on larger datasets. Higher number of structures would allow easier detection and removal of the outlier data and better sampling of data in all latitudes, thus reducing the parameter error in the curve-fitting function. Outliers could be due to human error, wrong database entry, calculation error, or eventually the algorithm failure. Detection and removal of these values will help improve the quality of results for the sidereal angular rotational velocity calculation. With larger datasets and sampling of data from all latitudes, the detection and removal of outlier data would be easier. We expect that an automatic outlier detection algorithm based on Peirce’s criterion would be suitable as future work.

In summary, this paper’s objective was to introduce the PSO-Snake algorithm to the solar physics community and report early results on CBP tracking problems. The algorithm showed the feasibility of using it for solar image analysis by addressing a current solar physics question. At the current state, the algorithm is suitable for tracking solar features, pre-marked by a user or indicated as an outcome from another system. As future work, authors are extending the PSO-Snake algorithm to make it fully automated, including initial detection of the CBPs. This will allow to analyze CBPs automatically, for extensive number of solar images and compare with the results of McIntosh & Gurman (2005), obtained by analyzing nine years of SOHO EUV images. Also it should be noticed that JPEG images were used in this study for their ease of use and as a proof of concept for the algorithm, but we will extend this work to use science-grade FITS images, for improved quality and precision.

Acknowledgments

We thank the SDO (NASA) and AIA science team for the provided observational material. This work was partially supported by Grants SFRH/BPD/44018/2008 (I.D.) and SFRH/BD/62249/2009 (E.S.) from Fundação para a Ciência e a Tecnologia (FCT), MCTES, Portugal, and by FCT Strategic Program UID/EEA/00066/203 of UNINOVA, CTS. The editor thanks three anonymous referees for their assistance in evaluating this paper.

References

All Tables

All Figures

	Fig. 1. A sample solar image from the 5th October 2010 with several CBPs marked. (Image courtesy of NASA/SDO).
In the text

	Fig. 2. Search mechanism of PSO particle in a two-dimensional space. (a) Shows particle x at time t under influence of its inertia, cognitive and social velocity components, (b) new position of the particle at time t + 1 is the algebraic sum of the velocities acting upon it.
In the text

	Fig. 3. Illustration of sample image force as a basis for calculating external snake energy.
In the text

Fig. 4. Comparison of the derived values of the rotational speed. Points in the chart represent the measurements for the center of point-like structures with error bars showing the 95% confidence level intervals. The dotted curve shows the fit to the mean ω(b) values as a function of latitude b calculated by Eq. (8). Overplotted are the results of Howard & Harvey (1970) in solid line and Hara (2009) and Brajša et al. (2004) both in the dashed-dotted curve because they are almost identical. Image courtesy of Lorenc et al. (2012).

In the text

	Fig. 5. A screenshot of the CBP Tracking system showing the Sun on 2010-10-05. (a) First snake (red circle) is initiated around the CBP we want to track, (b) zoomed-in view of the initial snake shows that the initial snake is a rough estimation of the location of the CBP. The algorithm will evolve the snake to find the exact CBP contour.
In the text

	Fig. 6. Detection process of a selected CBP (2010-10-05). Initial contour around the region of interest is evolved to precisely delineate the CBP boundaries. The outer contour is the initial snake and the inner contours are the transitional contours of every 10 iterations of the PSO-Snake detection algorithm. In the final stage contours converge and do not change much through iterations.
In the text

Fig. 7. (a) Red contour represents the initial snake on first image (the Sun at 2010-10-05), (b) and (c) cropped view of the tracking process of the selected CBP during time, showing the CBP in frame 1 and frame 37, respectively. The cyan contour is the boundary of the tracked CBP, the red square is the expert’s manual CBP positioning result, and the yellow circle is the PSO-Snake hybrid algorithm’s automated tracking result for CBP’s center of mass.

In the text

	Fig. 8. Deviation of the obtained results for b, ω, Δω on test data from the manual CBP marking of Lorenc et al. (2012). Every point in the horizontal axis represents a test CBP structure and its position in the vertical axis represents its error from the manual marking.
In the text

Fig. 9. Left: The result of Levenberg-Marquardt fitting algorithm (LMA) with least absolute residuals (LAR), on different number of tracked bright point structures. Right: Comparison of the curves based on Eq. (8) from different works. The black line represents this work, and it should be noticed that the similarity between the obtained results from this method with different numbers of tracked structures ((a) 650 CBP measurements and (b) 1577 CBP measurements) indicates the stability of this method.

In the text