Border Figure Detection Using a Phase Oscillator Network with Dynamical Coupling

Abstract

Oscillator networks have been developed in order to perform specific tasks related to image processing. Here we analytically investigate the existence of synchronism in a pair of phase oscillators that are short-range dynamically coupled. Then, we use these analytical results to design a network able of detecting border of black-and-white figures. Each unit composing this network is a pair of such phase oscillators and is assigned to a pixel in the image. The couplings among the units forming the network are also dynamical. Border detection emerges from the network activity.

1. Introduction

The synchronous firing of neurons seems to be important for accomplishing cognitive tasks such as attention (e.g., [1]), comprehension (e.g., [2]), coordination (e.g., [3]), perception (e.g., [4]), and sensory segmentation (e.g., [5]). Such experimental findings have inspired works about image processing (e.g., recognition, segmentation, symmetry detection) based on the synchronism of coupled oscillators (e.g., [6–11]). In addition, several neural systems, such as spinal cord (e.g., [12, 13]), hippocampus (e.g., [14, 15]), and visual cortex (e.g., [16, 17]), have been modelled by phase oscillators—one of the simplest oscillator models. In spite of its simplicity, this approach is suitable because more complex models for neurons like pulse-coupled and Hodgkin-Huxley-type models can be transformed into a phase oscillator model through coordinate changes [18]. Networks of phase oscillators can be electronically built using phase-locked loops (PLLs) (e.g., [19–22]).

Here we analytically study a phase oscillator model with dynamical coupling, and then use it in order to form a network capable of successfully detecting border of black-and-white figures. Each oscillator of our network corresponds to a first-order PLL (e.g., [8, 11, 17, 23]), which is equivalent to an overdamped pendulum.

The aim of any border detection process is to capture the main structural properties of an image, which can be useful, for instance, for image compression or diagnosis in echocardiography (e.g., [24]). Usually, this kind of image processing involves the use of partial derivatives (gradient, Laplacian; e.g., [25]), which are discretized in the space domain and in the time domain. Here we propose a scheme based on a network where each unit is formed by a pair of dynamically coupled PLLs corresponding to a pixel of the image. Border detection is obtained from the network activity. Notice that such a model can be transformed into a dedicated hardware for executing this image processing task. Also it naturally presents space discretization, because it is composed of a finite number of PLLs. This feature plus the fact that the time domain does need to be discretized can reduce the stability problems related to the algorithms employed for numerically calculating partial derivatives.

2. Model of a Single Unit

Consider that the temporal activities of two first-order PLLs with dynamical coupling are described by the following equations:𝑑𝜃₁𝑑𝑡=𝜔₁+𝑘₁₂sin(𝜃₂−𝜃₁),𝑑𝜃₂𝑑𝑡=𝜔₂+𝑘₂₁sin(𝜃₁−𝜃₂),𝑑𝑘₁₂𝑑𝑡=𝜇𝛼₁cos(𝜃₁−𝜃₂)−𝜇𝑘₁₂,𝑑𝑘₂₁𝑑𝑡=𝜇𝛼₂cos(𝜃₂−𝜃₁)−𝜇𝑘₂₁, (2.1)where 𝜃_𝑗 and 𝜔_𝑗(𝑗=1,2) are the phase and the natural frequency of the oscillator 𝑗, respectively; 𝑘_𝑖𝑗 (𝑖,𝑗=1,2 with 𝑖≠𝑗) is the connection strength from the oscillator 𝑗 to the oscillator 𝑖; 𝛼_𝑗 (𝑗=1,2) is the coefficient related to the Hebbian connection modification; and 𝜇 represents a natural (exponential) decay. The parameters 𝜔_𝑗, 𝛼_𝑗, and 𝜇 are positive numbers. Thus, the connection strengths are enhanced when the oscillators are in phase, and they are weakened when they are out of phase. This kind of (synaptic) modification between two oscillators (neurons) was suggested by Hebb [26].

By defining 𝑞≡𝜃₂−𝜃₁; 𝜔≡𝜔₂−𝜔₁; 𝛼₀≡(𝛼₁+𝛼₂)/2; 𝑘≡(𝑘₁₂+𝑘₂₁)/2, system (2.1) can be rewritten as𝑑𝑞𝑑𝑡=𝑓₁(𝑞,𝑘)=𝜔−2𝑘sin𝑞,𝑑𝑘𝑑𝑡=𝑓₂(𝑞,𝑘)=𝜇𝛼₀cos𝑞−𝜇𝑘. (2.2)The variables 𝑞 and 𝑘 represent the phase difference and the average connection strength, respectively. Notice that 𝑓_𝑗(𝑞,𝑘)=𝑓_𝑗(𝑞+2𝜋,𝑘) (𝑗=1,2). The formation of synchronized clusters in networks of phase oscillators described by system (2.2) was investigated by Seliger et al. [27]. In this section we study the asymptotic behaviors of this model.

Synchronism occurs when both nodes oscillate in a common frequency, which means that 𝑑𝜃₂/𝑑𝑡=𝑑𝜃₁/𝑑𝑡 or 𝜃₂−𝜃₁=𝑞=constant. Thus, synchronism can happen if there are constants (𝑞^∗,𝑘^∗) so that if (𝑞(𝑡),𝑘(𝑡))=(𝑞^∗,𝑘^∗) then 𝑑𝑞/𝑑𝑡=0 and 𝑑𝑘/𝑑𝑡=0 for all time 𝑡. Notice that a synchronous solution is an equilibrium point (𝑞^∗,𝑘^∗) of the dynamical system (2.2). Such a solution corresponds to an intersection point of the nullclines 𝑓₁=0 (𝑘=𝜔/(2sin𝑞)) and 𝑓₂=0 (𝑘=𝛼₀cos𝑞).

The stability of an equilibrium point (𝑞^∗,𝑘^∗) can usually be determined by calculating the eigenvalues 𝜆 of the Jacobian matrix corresponding to system (2.2) linearized around this point. It is asymptotically stable when all eigenvalues have negative real parts (e.g., [28]). For system (2.2), the eigenvalues 𝜆 are the roots of the polynomial𝜆²+(𝜇+2𝑘^∗cos𝑞^∗)𝜆+2𝜇𝑘^∗cos𝑞^∗−2𝜇𝛼₀(sin𝑞^∗)²=0. (2.3)According to the Routh-Hurwitz criterion, both roots of the polynomial 𝜆²+𝑎₁𝜆+𝑎₂=0 have negative real parts if 𝑎₁>0 and 𝑎₂>0 (e.g., [29]).

When 𝜔=0, the equilibrium points are 𝑃₁=(0,𝛼₀), 𝑃₂=(𝜋/2,0), 𝑃₃=(𝜋,−𝛼₀), 𝑃₄=(3𝜋/2,0), where 𝑃₁ and 𝑃₃ are asymptotically stable, and 𝑃₂ and 𝑃₄ are unstable.

When 𝜔≠0, there are also four equilibrium solutions if 𝜔<𝛼₀. The points with𝑘^∗𝛼=±⁰√21+𝜔1−(𝛼₀)² (2.4)are asymptotically stable; the ones with𝑘^∗𝛼=±⁰√21−𝜔1−(𝛼₀)² (2.5)are unstable. When 𝜔=𝛼₀, there are only two equilibrium points because two saddle-node bifurcations (e.g., [28]) occur. When 𝜔>𝛼₀, there is not synchronism; however, there is a limit cycle, which corresponds to a closed and isolated trajectory in the state space 𝑞×𝑘.

For an autonomous two-dimensional system, Poincaré-Bendixson theorem (e.g., [28]) ensures that there is an asymptotically stable limit cycle in a region of the state space if the vector field (𝑓₁,𝑓₂) points inward everywhere on the boundary of this region, which must not contain any equilibrium point. For system (2.2), it is easy to verify that if 𝜔>𝛼₀, then there is such an attracting trajectory inside the region 𝑅, where 𝑅 is the rectangle given by 𝑅={(𝑞,𝑘)∶0≤𝑞<2𝜋,𝛼₀+(𝜖/𝜇)<𝑘<−𝛼₀−(𝜖/𝜇)} for 𝜖→0₊. In this case, 𝑘 oscillates and 𝑞 grows as the time goes by, as shown in Figure 2. Figure 1 illustrates the case with synchronous (equilibrium) asymptotic solution.

Figure 1 

Temporal evolution of 𝑞(𝑡) and 𝑘(𝑡) obtained by numerically integrating system (2.2). Parameter values: 𝜔=1, 𝛼=2, 𝜇=1. Initial conditions: 𝑞(0)=0, 𝑘(0)=5. In this case, the asymptotic solution corresponds to the equilibrium point (𝑞∗,𝑘∗)=(0.26,1.93).

Figure 2 

Temporal evolution of 𝑞(𝑡) and 𝑘(𝑡) obtained by numerically integrating system (2.2). Parameter values: 𝜔=10, 𝛼=2, 𝜇=1. Initial conditions: 𝑞(0)=0, 𝑘(0)=5. In this case, the asymptotic solution corresponds to the limit cycle, where 𝑘(𝑡) oscillates with frequency 𝐵≃9.8 and amplitude 𝐴≃0.20, and 𝑞(𝑡) linearly grows with slope 𝐴.

In order to analytically characterize the limit cycle, suppose that the asymptotic behavior of 𝑞(𝑡) in this case is given by 𝑞_asymp(𝑡)≃𝐵𝑡. By inserting this approximated (linear) solution in the equation for 𝑑𝑘/𝑑𝑡 and integrating it, the asymptotic solution of 𝑘(𝑡) is𝑘_asymp(𝑡)≃𝐴cos(𝐵𝑡−𝜑), (2.6)where 𝐴=𝜇𝛼₀/𝜇²+𝐵² and sin𝜑=𝐵/𝜇²+𝐵². By substituting this expression in the equation for 𝑑𝑞/𝑑𝑡 in (2.2), an equation for calculating 𝐵 is found. In fact, 𝐵 is the root of the polynomial𝐵³−𝜔𝐵²+(𝜇²+𝜇𝛼₀)𝐵−𝜔𝜇²=0. (2.7)Observe that when 𝜔≫1, then the angular velocity 𝐵 can be written as 𝐵≃𝜔+Δ, with |Δ|≪1. The value of Δ is estimated byΔ≃−𝜇𝛼₀𝜔𝜔²+𝜇²+𝜇𝛼₀. (2.8)Consequently, the oscillation amplitude 𝐴 of 𝑘_asymp(𝑡) when 𝜔≫1 is𝐴≃𝜇𝛼₀𝜔. (2.9)

For instance, for 𝜔=10, 𝛼₀=2, and 𝜇=1, these approximated expressions give 𝐵≃9.8 and 𝐴≃0.20, which are in good agreement with the numerical solutions of system (2.2) shown in Figure 2.

3. Network for Detecting Figure Border

Our two-dimensional network for detecting figure border is composed of units consisting of a pair of phase oscillators forced by an external input 𝐼. In this network, the coefficients 𝛼₁ and 𝛼₂ can vary with the time. Thus, system (2.2) is rewritten as𝑑𝑞𝑑𝑡=𝐼+𝜔−2𝑘sin𝑞,𝑑𝑘𝑑𝑡=𝜇𝛼(𝑡)cos𝑞−𝜇𝑘. (3.1)

Assume that a black-and-white figure will be presented to this network, as illustrated in Figure 3(a). The variable 𝐼 represents the input for each unit according to the following rules: 𝐼=0 corresponds to a white part of the figure; 𝐼=9 corresponds to a black part. Notice that the input can be translated into a new natural frequency Ω by defining 𝐼+𝜔≡Ω. Thus, system (3.1) is reduced to system (2.2) if 𝛼(𝑡) is a constant.

(a)

(b)

Figure 3 

(a) At 𝑡=0, this is the image presented to the network with 9×9 units. (b) At 𝑡=40, the border of the image is dynamically determined.

Each unit is coupled with its four closest neighbors (left, right, up, and down), which is usually known as four-neighborhood in image processing literature (e.g., [30]) or von Neumann neighborhood in cellular automaton literature (e.g., [31]). As a consequence of this coupling, the value of 𝛼(𝑡) for each cell is given by𝛼𝛼(𝑡)=⁰if0≤𝑡≤𝑇,𝑝𝛼₀if𝑡>𝑇. (3.2)The parameter 𝑇 is a settling time; 𝑝≡1+𝑞𝑥 where 𝑞 is a positive number (here 𝑞=14); 𝑥=𝐺(𝑦) where 𝐺(𝑦)=0 if 𝑦=0; and 𝐺(𝑦)=1 if 𝑦≠0. The value of 𝑦 is obtained by𝑦=⁴_𝑖=1[𝑘_𝑖(𝑇)−𝑘_𝑖(𝑇−𝛿𝑇)], (3.3)where 𝛿𝑇 is the time step (here 𝛿𝑇=0.01) of the integration method (here fourth-order Runge-Kutta method) used for numerically solving the dynamical system. The variable 𝑘_𝑖(𝑡) corresponds to the average connection strength of the closest neighbor 𝑖(𝑖=1,2,3,4). If the neighbor 𝑗 does not exist, then 𝑘_𝑗(𝑡)≡0.

An image is presented to the network at 𝑡=0. For 0≤𝑡≤𝑇, all the units behave as if they were isolated, because the value of 𝛼(=𝛼₀) for each unit is independent of the neighbor activity. By taking for all units 𝜔=1, 𝛼₀=2, 𝜇=1, 𝑇=20, 𝑞(0)=0, and 𝑘(0)=5, the units with 𝐼=0 will tend to a stationary solution and the units with 𝐼=9 will tend to a limit cycle, as explained in the last section. For these parameter values, the permanent regime was already reached when 𝑡=𝑇, as shown in Figures 1 and 2.

At 𝑡=𝑇, the value of 𝛼 for each unit can be changed. Expressions (3.2) and (3.3) imply what follows. If all of its four neighbors are in a limit cycle (𝑦≠0), then the value of 𝛼 is increased from 2 to 30; if at least one neighbor is in a stationary state or if the unit does not have a complete von Neumann neighborhood (𝑦=0), then the value of 𝛼 is kept equal to 2. The limit cycle is characterized by 𝑘(𝑡) oscillating with frequency 𝐵 and amplitude 𝐴, and 𝑞(𝑡) linearly growing with slope 𝐴, as presented in the last section.

For 𝛼=30 and Ω=10, 𝑞(𝑡) and 𝑘(𝑡) tend to a stationary activity, which in our model corresponds to the white color. Thus, at 𝑡=2𝑇=40, the unique units remaining in oscillatory activity (the black color) are the ones in the figure border. Figure 3(b) shows the contour of Figure 3(a), after the units have reached the (new) permanent regime.

The higher the value of 𝜇 is, the shorter the transitory phase will be. Hence, border detection can be made faster by increasing the value of this parameter.

4. Conclusion

Networks governed by differential (e.g., [6–8, 10, 11]) or difference equations (e.g., [9, 32]) have been employed for image processing. Here we used a network with local and dynamical coupling for identifying contour of black-and-white images. The unit of such a network is the variant of the Kuramoto model [33] for two oscillators (neurons), which was proposed by Seliger et al. [27]. In our scheme, after applying the input corresponding to the figure, the value of 𝛼 for each unit is switched or not, depending on the neighborhood activity. This scheme could be implemented using first-order PLLs. This investigation could be extended for higher-order PLLs, which could be designed in order to shorten the transient.

Acknowledgment

L. Monteiro and J. Piqueira are partially supported by CNPq.