New Model and Applications of Cellular Neural Networks in Image Processing

1. Introduction

The cellular neural network (CNN) in its standard form as defined by Chua and Yang has become a paradigm of analog, high-speed computation, with various applications in parallel signal processing, especially in image processing tasks. Within two decades this domain has extensively developed both theoretically and from the applications point of view. Although the mainstream research basically followed the standard CNN model, many other CNN models and architectures were proposed, envisaging new and more complex classes of applications.

We propose here an alternative CNN model, inspired from Hopfield-type recurrent networks. This new model has the same general architecture as the standard model, while the cell structure, preserving the basic elements, is modified from a dynamical point of view. The cell non-linearity, having a major role in CNN dynamics, no longer applies only to the current state variable as it does in the standard model, but to the whole information fed to the current cell, thereby leading to a different dynamic behavior. The new cell structure and state equations are presented and the dynamics and equilibrium points are analyzed.

The advantages of this model, which we call “CNN with overall non-linearity”, are a simpler design of the templates, and the fact that all voltages in the cells are limited to a fixed dynamic range regardless of template parameters; thus the boundedness problem regards the current range instead of voltage, which makes it less critical from an implementation point of view. The image processing capabilities of the proposed CNN model will be approached. Some new image processing tasks are proposed. We first analyze the behavior of the proposed CNN in processing binary images, then using so-called local rules we design the templates for several specific tasks, which are: convex corner detection, erasing thin lines, detection of closed curves, finding all paths between two selected points through a labyrinth, rotation detector, marked object reconstruction and finding the intersection points of one-pixel thin lines from two superposed binary images. Many other processing tasks both on binary or grayscale images are possible. Further research may focus on extending the range of applications of the proposed model.

Next some applications of CNNs in image linear filtering will be approached. In this range of applications, the cells of the CNN operate in the linear domain, and the system stability must be ensured. The design of CNN linear filters with high selectivity may lead in general to large-size templates. Due to the fact that local connectivity in CNNs is an essential restriction, the decomposition of a large-size template into small-size templates is often necessary. Fast linear filtering using parallel computing systems like CNNs may be useful in early vision aplications and pre-processing tasks. We will present some original and efficient design methods for several types of two-dimensional filters: circular, maximally-flat, elliptically-shaped, orientation-selective filters etc. We also propose a class of 2D linear filters which are designed starting from a desired shape in the frequency plane, described in polar coordinates. Using this method, some interesting filters can be designed, such as concave-shaped filters presenting narrow lobes in the frequency plane. A particular class are the zero-phase 2D filters, which are extensively used in image processing due to the absence of phase distortions. For all the proposed filters, the design starts from an imposed 1D prototype and uses rational approximations and various frequency transformations to obtain the 2D frequency response. The resulted filters are efficient, at the same time of low complexity and relatively high selectivity. The design methods are entirely analytical, but they could be combined with numerical optimization techniques. The filters can be implemented either digitally or on analog processing systems like CNNs. Further research also envisages an efficient and reliable implementation of this class of filters.

2. Cellular Neural Networks with Overall Nonlinearity

2.1. Equations and Circuit Model

The state equation of the current cell i for a standard CNN has the form (Chua & Yang, 1988):

$$\frac{d{x}_i}{dt}=-x_i+{\displaystyle \sum _{j\in N_i}{A}_j\cdot f(x_j)}+{\displaystyle \sum _{k\in N_i}{B}_k\cdot u_k}+I_i$$ E1

and the output is given by the piece-wise linear function:

$$y_i=f(x_i)=0.5\left(\left|x_i+1\right|-\left|x_i-1\right|\right)$$ E2

Here a new version of CNN is proposed (Matei, 2001), with a modified structure as compared to the standard one. The difference lies in the form of the state equation, which has the form:

$$\frac{d{x}_i}{dt}=-x_i+f\left({\displaystyle \sum _{j\in N_i}{A}_j\cdot x_j}+{\displaystyle \sum _{k\in N_i}{B}_k\cdot u_k}+I_i\right)$$ E3

For simplicity we used a single index for the cells, as in the case of 1D CNNs. Unlike the standard CNN, in this case the nonlinear function applies to the entire information gathered from the neighboring cells, i.e. state variables ( $$x_j$$ ), inputs ( $$u_k$$ ) and the bias term ( $$I_i$$ ). For this reason, we have called this a CNN with "overall non-linearity". This CNN version was inspired by the well-known dynamical model of a neuron which is the elementary unit of so-called recurrent neural networks or Hopfield networks. Within this framework, the standard CNN can be regarded as a particular Hopfield network with local connectivity. The neuron's output is related to the current state ("activation potential") through a nonlinear activation function with various shapes (logistic, piece-wise linear,etc.). There are two basic dynamical neuron models, related by a linear and generally invertible transformation (Haykin, 1994). The proposed CNN corresponds to the second recurrent model. Therefore, the proposed CNN model and the Chua-Yang model are equivalent, being described by similar equations if a linear transformation is applied.

Let us denote by $$z_i$$ the argument of the piece-wise linear function $$f$$ in (3):

$$z_i={\displaystyle \sum _{j\in N_i}{A}_j\cdot x_j}+{\displaystyle \sum _{k\in N_i}{B}_k\cdot u_k}+I_i={\displaystyle \sum _{j\in N_i}{A}_j\cdot x_j}+k_i$$ E4

where the inputs and bias are constant and were put together into the term $$k_i$$ .

Multiplying both members of (3) by $$A_j$$ , summing over $$i\in N_i$$ ( $$N_i$$ denotes the current cell neighborhood) and shifting again the indices i and j, we finally obtain (Matei, 2001):

$$\frac{d{z}_i}{dt}=-z_i+{\displaystyle \sum _{j\in N_i}{A}_j\cdot f(z_j)}+{\displaystyle \sum _{k\in N_i}{B}_k\cdot u_k}+I_i$$ E5

having a form similar to (1), but in the new variable $$z_i$$ . Equation (4) defines the linear transformation which relates the two CNN models.

2.2. System and Circuit Description

The modified CNN will be analytically described from dynamical point of view, as well as at circuit level. Figure 1 shows a block diagram representation of the dynamic behavior of the modified CNN, which can be compared to the standard CNN block diagram (Matei, 2001).

In this block diagram all the variables are vectors or matrices. The correlation sums in the state equation can be written as convolutions. The variable z at the input of the nonlinear block is:

$$z=A\ast x+B\ast u+I$$ E6

and the matrix state equation of the modified CNN will be:

$$dx/dt=-x+f\left(A\ast x+B\ast u+I\right)$$ E7

Unlike the standard CNN case (where the nonlinear block is placed in the feedback loop), in

Figure 1.

(a) Block diagram of proposed CNN model; (b) Circuit structure of a cell.

the modified CNN it operates in the forward path. Here we must emphasize an important feature of the proposed CNN: the output is identical to the state variable $y(t)=x(t)$ throughout the operating domain – a significant advantage in the CNN implementation.

In Figure 1(b) we present the cell circuit of the modified CNN, which preserves the same elements as the standard CNN, but has a different structure. The major difference between the two structures lies in the location of the nonlinear element. In the standard CNN case, the nonlinear function applies to the state variable $x_i$ , which is a voltage, and gives the output voltage $y_i=f(x_i)$ . In the modified CNN equation (3), the non-linearity applies to the amount denoted $z_{ij}$ , summing up the entire information from the states and the inputs of the neighboring cells. The variable $z_{ij}$ is a current (given by several linear controlled current sources), and $f(z_i)$ is a current as well. For the piecewise-linear current-controlled current source we can write $I_{yx}=f(z_i)$ . The resistor denoted $R_y$ has the single role of summing the currents given by the linear sources, while in the standard CNN cell it converts voltage into current. Unlike the Chua-Yang model, in this case the nonlinear current source feeds the current directly into the capacitor C and the resistor $R_x$ . On the cell circuit, considering normalized values for C and $R_x$ we can also derive the state equation (3).

2.3. Dynamics and Equilibrium Points

Using the Lyapunov method and the linear mapping (4) it can be proven that the modified CNN also has the property of convergence to a constant steady state following any transient. The stability properties are also established in Hopfield network theory (Haykin, 1994).

For standard CNNs, the state equation (1) is linear as long as all the state variables are less than one in absolute value, $x_i\in [-1,1]$ . This domain defines the linear region of the state space. For the modified CNN, the linear region of the state space in variable $x_i$ depends on the feedback state information, the input information and the bias term, as shows (3). In the linear region of the state space, as long as no cell reached saturation, equations (1) and (3) have a similar form, since in this domain we have $f(x)=x$ , which implies that the linear dynamics of the standard CNN and modified CNN are identical. The difference lies in the way in which the two types of CNNs enter the nonlinear region and reach equilibrium.

In the following paragraph we proceed to an analysis of the dynamic behavior and a study of the equilibrium points of the modified CNN. We will use the so-called driving-point plot, which gives the state derivative $d{x}_i/dt$ vs. state variable $x_i$ (Chua & Yang, 1988). We will emphasize the dynamic routes and the possible equilibrium points the system can reach.

Using equation (5) which corresponds to a standard CNN in the state variable $z_i$ , we can first represent the driving-point plot ( $d{z}_i/dt,z_i$ ) in the new variable $z_i$ as in Figure 2 (a). The drawing of this plot is shown in detail in (Matei, 2001). In order to derive the driving point plot in the original variable $x_i$ , we use the linear mapping (4) in the form:

$$z_i(t)=p\cdot x_i(t)+g_i(t)$$ E8

$$g_i(t)={\displaystyle \sum _{j\in N_i,j\ne i}{A}_j\cdot x_j}+{\displaystyle \sum _{k\in N_i}{B}_k\cdot u_k}+I_i={\displaystyle \sum _{j\in N_i,j\ne i}{A}_j\cdot x_j}+k_i$$ E9

In (8), $p=A_0$ is the central element of the feedback template. Supposing for the moment a fixed, constant value for $g_i(t)$ , denoted $g_i$ , we obtain for $p\ne 0$ :

$$x_i(t)=(1/p)\cdot \left(z_i(t)-g_i\right)$$ E10

$$d{z}_i/dt=p\cdot d{x}_i/dt$$ E11

The case $p=0$ should be analyzed separately. Through a simple graphical construction shown in Figure 2 (a) we can draw the corresponding driving point plot in the original variable $x_i$ . This can be done using two auxiliary plots: (a) the plot $(x_i,z_i)$ , given by (10), represented by a straight line with slope of value $1/p$ and offset value $g_i$ ; (b) the plot $(d{z}_i/dt,d{x}_i/dt)$ , given by (11), consisting in a straight line with slope p crossing the origin.

The plots depicted in Figure 2 (b) were drawn for the particular values $p=2$ , $g_i=0.5$ .

The segment $A_1{B}_1$ with slope ( $p-1$ ), given by:

$$d{z}_i/dt=(p-1)\cdot z_i+g_i$$ E12

has an identical form in $x_i$ :

$$d{x}_i/dt=(p-1)\cdot x_i+g_i$$ E13

The points $A_1$ , $B_1$ , given by the coordinates: $A_1(1,p-1+g_i)$ , $B_1(-1,-p+1+g_i)$ become:

$$A\left(\frac{1-g_i}{p},\frac{p-1+g_i}{p}\right),B\left(\frac{-1-g_i}{p},\frac{-p+1+g_i}{p}\right)$$ E14

The lines $A_1{D}_1$ and $B_1{C}_1$ , described by the equations:

$$d{z}_i/dt=-z_i+p+g_i d{z}_i/dt=-z_i-p+g_i$$ E15

become $AD$ and $BC$ respectively, given by:

$$d{x}_i/dt=-x_i+1 d{x}_i/dt=-x_i-1$$ E16

Figure 2.

(a) Graphical construction for determining the driving point plot; (b) Driving point plot, dynamic routes and equilibrium points for different values of $g_i$ .

In the diagram $\left(d{x}_i/dt,x_i\right)$ the points C and D result from the intersection of the straight lines BC and AD with axis $x_i$ , imposing the condition: $d{x}_i/dt=0$ .

We finally obtain: $x_i(C)=-1$ and $x_i(D)=1$ . Points $C(-1,0)$ and $D(1,0)$ will therefore be fixed equilibrium points of the modified CNN, independent of the values of p and $g_i$ . For different values of the amount $g_i$ , the diagram $(d{z}_i/dt,z_i)$ shifts up or down. Correspondingly, the $(d{x}_i/dt,x_i)$ plot slides along the two parallel dashed lines crossing the two fixed points C and D. In Figure 2 (b) one typical driving point plot is shown with continuous line and the dynamic routes are indicated by arrows. The number of equilibrium points and their nature for different values of p and $g_i$ are analyzed in (Matei, 2001).

2.4. Dynamic Range of Modified CNNs

For standard CNNs the following result has been proven in (Chua & Yang, 1988):

Theorem: All states $x_i$ in a CNN are bounded for $t0$ and the bound $v_{\max }$ is given by:

$$v_{\max }=1+\left|I\right|+\underset{i}{\max }{\displaystyle \sum _i\left(\left|A_i\right|+\left|B_i\right|\right)}$$ E17

if we consider the normalized resistance value $R_x=1$ .

The possibility of estimating the maximum dynamic range is very important from the implementation point of view, because it allows one to optimally design the circuit parameters in order to keep the voltage values within reasonable bounds.

In (Rodriguez-Vasquez et al., 1993) a new CNN model, called "full-range" model was introduced, adding a nonlinear term to the state equation, which confines the state value to the range [-1,1]. This simplifies the design, allowing for reduced area and power consumption in VLSI implementation.

Our approach gives in fact another solution to this problem, forcing the state and output values to be equal. Therefore, the proposed CNN with "overall non-linearity" is in fact a "full-range" model as well. In the modified CNN, this problem is solved by the form of the state equation itself. Indeed, at any equilibrium point, the current cell state variable is $x_i=f(z_i)$ and consequently its variation is bounded to the range $[-1,1]$ . Moreover, the output variable is identical to the state variable at any time: $y_i(t)=x_i(t)$ .

Referring again to the circuit structure, the argument $z_i(t)$ of the nonlinear function $f(z_i)$ given by (4) is physically a current. However, since equation (5) is similar to the state equation of a standard CNN cell, but written in variable $z_i$ , equation (17) may be formally applied to evaluate a maximum bound $I_{\max }$ for the cell currents.

Therefore in the case of proposed CNN with overall nonlinearity, the boundedness problem regards the current range instead of voltage, which makes it less critical from an implementation point of view. All voltage variations in the modified CNN are limited to a fixed range corresponding to the extreme pixel values, regardless of template parameters.

3. Applications of CNNs with Overall Nonlinearity

In this section we approach the image processing capabilities of the proposed CNN model. Some processing tasks will be proposed on binary (black and white) images. We first analyze the behavior of the proposed CNN model in processing binary images, then we design the templates for a few specific tasks, namely: convex corner detection, erasing thin lines, detection of closed curves, finding all paths between two selected points through a labyrinth, rotation detector, marked object reconstruction and finding the intersection points of one-pixel thin lines from two superposed binary images. Many other processing tasks on binary or grayscale images are possible.

A vast reference for different classes of standard CNN tasks is the CNN Software Library, a template catalog for a large range of analog processing operations on binary and grayscale images. These also include more complex analog algorithms which consist of sequences of elementary tasks which are fulfilled under digital control.

We use here the convention that the pixel value -1 corresponds to white, and the pixel value +1 corresponds to black. Generally, one binary image is loaded into the CNN as initial state, while another binary image may be applied to the array input. In this range of applications the system must behave as a bipolar CNN. A standard CNN is said to be bipolar if $u,y\in B^{MN}$ , $x(0)\in B^{MN}$ , where $B=\{-1,1\}$ , and M, N are the number of rows and columns of the CNN array (Hänggi & Moschytz, 2000).

Based on the information contained in the two images, the CNN must accomplish a given task upon the state image, according to a set of local rules for which the feedback and control templates are designed. If a given binary image is loaded as initial state, this implies that all cells are initially saturated. Depending on the chosen template parameters, the initial state of a given cell may be stable or unstable. If the state is stable, the corresponding pixel will remain unchanged (black or white); if on the contrary the state is unstable, the pixel may switch to the opposite state, following a transient evolution. Taking into account that we want to obtain a binary output image in response to the processed image, we impose that the cells which are initially unstable switch into the opposite state, therefore for a generic cell $C(i,j)$ we will have: ${\dot{x}}_{ij}(0)\ne 0$ . These cells will switch into the opposite saturated state which is supposed to be stable according to local rules. The switching of a cell between two saturated states implies a dynamic evolution through the middle linear region. As stated from the beginning in (3), the dynamics of the current cell is described by:

$$d{x}_i(t)/dt=-x_i(t)+f\left(z_i(t)\right)$$ E18

Suppose that a given pixel is initially black, i.e. $x_i(0)=1$ . If it is to remain black, this implies $${\dot{x}}_i(0)=0$$ and from (18) we must have: $$f\left(z_i(0)\right)=1$$ , therefore $z_i(0)1$ . If the black pixel must switch to white, we impose: $z_i(0)-1$ , so $$f(z_{{}_i}(0))=-1$$ ; the initial state derivative is:

$${\dot{x}}_i(0)=-x_i(0)-1=-20$$ E19

The cell state value decreases exponentially from 1 to -1:

$$x_i(\infty )=-1 {\dot{x}}_i(\infty )=-x_i(\infty )-1=0$$ E20

therefore the cell finally reaches a stable equilibrium point.

In the following paragraphs we present several tasks on binary images in order to show the image processing capabilities of the proposed CNN model with overall nonlinearity. For each of these tasks, the templates are designed and simulation results are provided.

3.1. Convex Corner Detection

Although the convex corner detection on binary images is implemented as well on standard CNNs, we show here the template design method for the proposed model. We assume a feedback template A of the general symmetric form:

$$A=\left[\begin{array}{ccc}s& s& s\\ s& p& s\\ s& s& s\end{array}\right]$$ E21

The first step is to find the local rules according to the desired result. We have to take into account all the possible combinations of black and white pixels within a $3\times 3$ frame. Although the total number of combinations is quite large, due to template A symmetry the number of distinct combinations is largely reduced. All inner black pixels must turn white, while any edge pixel must remain black if it locates a convex corner or turn white otherwise. In Figure 3 (a), the combinations (a)–(d) represent convex corners, while cases (e)-(h) do not. Assuming for s a negative value, we obtain a system of inequalities which reduces to the set:

$$p+2s+I-1p+I1p+8s-I1$$ E22

Finding a solution to the above system is a linear programming problem which can be solved graphically or numerically. We finally reach the set of parameters: $p=7.5$ ; $s=-1.5$ ; $I=-6$ . A simulation of this task is shown in Figure 3 (b), (c). The result of convex corner detection is shown, resulted from the test image (b) loaded as initial state. We remark that all convex corners are correctly detected, while the concave corners do not appear in the output image (c). A thin line emerging from a convex corner is also detected.

Figure 3.

(a) Possible pixel combinations in convex corner detection; (b) Image loaded as initial state; (c) Output images containing detected corners.

3.2. Erasing Thin Lines

We intend to design a template which erases every thin line (one-pixel thick) from a binary image, leaving all other objects unchanged. For example, the test image in Figure 4(a) is loaded into the CNN as initial condition. We look for the feedback template A in the form (21), this time assuming $s0$ ; writing the local rules, we reach the following system:

$$p-4s+I-1p-2s+I1p-I1$$ E23

and we choose a convenient set of parameter values: $p=8$ ; $s=2$ ; $I=-2$ . The CNN input is not used. The final result of the processing is shown in Figure 4(c). We notice that the two thin curves from the left, as well as the spiral and the lines emerging from the black square have been erased. The lines perpendicular to the square edges may leave a single residual pixel. In Figure 4(b) a snapshot of an intermediate state is shown, in which the progressive deletion of thin lines is noticeable. The compact objects present in the image remain unchanged.

3.3. Detection of Closed Curves

In some applications it might be useful to discern between closed and open curves. These curves in turn may result from a previous processing, such as edge detection etc.

We design a CNN task which detects the one-pixel thin closed curves and deletes the open curves from a binary image. The absence of at least one pixel from a closed curve turns it into an open curve. Both the input and the initial state are used. The same binary image P is loaded as initial state and applied to input: $U(t)=P$ ; $X(0)=P$ . The templates are chosen as:

$$A=\left[\begin{array}{ccc}s& s& s\\ s& p& s\\ s& s& s\end{array}\right]; B=\left[\begin{array}{ccc}r& r& r\\ r& p& r\\ r& r& r\end{array}\right]; r=-0.5s$$ E24

Applying the corresponding local rules, the set of inequalities is derived:

$$2p-2s+I12p-3s-I12p-2.5s+I-1$$ E25

A convenient choice of parameters is: $p=9$ ; $s=6$ ; $I=-4.5$ . In this task we observed the 8-pixel connectivity, i.e. the current black pixel is considered connected to any other black pixel lying in its $3\times 3$ neighborhood. A simulation is shown in Figure 5. The test image (a) contains six thin curves, which feature all types of possible edges and corners. There is also a compact concave object in the image. The image (b) shows the CNN state at an intermediate moment, and the final image (c) contains just the closed curves and the compact object, which was preserved. The process of erasing open curves relies on a dynamic propagation starting from the free ends. The black pixels from the ends turn white according to the imposed local rules until the entire curve is erased.

Figure 4.

(a) Binary image; (b) Transient state; (c) Steady state output image.

Figure 5.

Closed curve detection: (a) binary image; (b) transient state; (c) output image.

3.4. Finding All Paths Between Two Selected Points Through a Labyrinth

A processing task related to the previous one is presented below. We consider a labyrinth made up of white one-pixel thin curves against a black background. In this case we suppose valid the convention that one pixel is possibly connected only to four neighbors (up, down, left, right). The target of the proposed task is to find all the existing paths between two given points situated in the labyrinth, which are marked by white pixels in the input image. A point is represented by one pixel, either on the input or state image.

The binary image P representing the labyrinth is loaded into the CNN as initial state ( $X(0)=P$ ). The input image may be zero everywhere (average grey level) except for two white pixels (of value -1) which mark the points selected in the labyrinth.

The dynamic process which finds the routes between the chosen points is similar to the one used in application 3.3. The template parameters are designed such that all labyrinth routes with unmarked ends are gradually erased starting from end pixels which turn from black to white (cell outputs switch from -1 to +1). The templates are searched in the form:

$$A=\left[\begin{array}{ccc}r& s& r\\ s& p& s\\ r& s& r\end{array}\right]; B=\left[\begin{array}{ccc}0& 0& 0\\ 0& a& 0\\ 0& 0& 0\end{array}\right]$$ E26

and the following system results:

$$-p+4r+I-1-p+2s-4r+I1-p+2s+2r+I1-p+2s+4r+I-a-1$$ E27

with convenient parameter values: $p=12$ ; $s=4$ ; $r=0.5$ ; $a=8$ ; $I=8$ .

In Figure 6 a simulation example is shown. The input image (b) contains two white pixels on a uniform gray background (zero pixel value). The two pixels locate on the labyrinth image (a) (loaded as initial state) the points between which the routes are selected. At the end of dynamic propagation shown in (c), the output image (d) results.

More generally, routes can be found among several points in the labyrinth, or a single point is chosen and we find all the closed routes containing that point. The classical problem of finding the shortest path through a labyrinth is more complex and it can be solved through a CNN algorithm using several nonlinear templates, which are run sequentially under digital control.

Figure 6.

(a) Labyrinth image as initial state; (b) input image; (c) transient state; (d) output image containing detected routes between selected points.

3.5. Erasing Objects with Tilted Edges (Rotation Detector)

We will design a CNN task which detects in a binary image the compact convex or concave objects having only horizontal and vertical edges and removes all tilted objects or objects having at least one tilted edge. In particular the CNN can detect the rotation of such objects.

This is based on the discrete representation on an array of a straight line, in particular the edge of a compact object, with a certain slope. Due to the limited spatial resolution of the CNN array, any tilted edge will be mapped as a stair-like edge. An edge with the minimum angle will map into a staircase with only two stairs, while an edge with a slope of $\pi /4$ is mapped into a staircase with a maximum number of one-pixel stairs. The binary image P is loaded into the CNN as initial state ( $$X(0)=P$$ ) and applied at the CNN input as well, $U(t)=P$ . We look for the templates in the general forms:

$$A=\left[\begin{array}{ccc}r& s& r\\ s& p& s\\ r& s& r\end{array}\right]; B=\left[\begin{array}{ccc}0.5r& -0.5s& 0.5r\\ -0.5s& p& -0.5s\\ 0.5r& -0.5s& 0.5r\end{array}\right]$$ E28

and finally the following convenient values are found: $p=5$ ; $s=5$ ; $r=-0.8$ ; $I=-11.2$ .

Figure 7.

(a) Binary image; (b) Transient state; (c) Steady state output image.

A simulation is presented in Figure 7. Image (a) contains several compact objects: a rectangle and a square placed horizontally and their rotated versions etc. The objects with at least one tilted edge are gradually eroded through a dynamic propagation, as shows the transient snapshot (b); the final image (c) preserves only objects with horizontal and vertical edges. Therefore, this task can be used to detect the rotation of simple objects (rectangles etc.).

A simulation is presented in Figure 7. Image (a) contains several compact objects: a rectangle and a square placed horizontally and their rotated versions etc. The objects with at least one tilted edge are gradually eroded through a dynamic propagation, as shows the transient snapshot (b); the final image (c) preserves only objects with horizontal and vertical edges. Therefore, this task can be used to detect the rotation of simple objects (rectangles etc.).

3.6. Marked Object Reconstruction

A binary image P containing compact shapes or curves is applied to the CNN input. From this image we intend to select certain objects which we want to preserve, and delete the others. The selection of objects is done by marking them with a black pixel, in the image loaded as initial state. Therefore at the input we have $U=P$ and the initial state $X(0)$ consists in black marking pixels against a white background. These can be placed anywhere within the area of selected object (in the interior or on the edge). The marked objects are reconstructed in the state image through a gradual process. Starting from the marking pixel, the white pixels within the marked object turn black. This propagation process stops on the edge of the object in the input image. The unmarked objects will not appear in the final image. We assume the feedback template of the form (21) and the control template zero, except the central element equal to a. With $s0$ , $p0$ , we finally obtain the system:

$$p-8s+I+a1p+8s-I-a1p+6s-I-a-1$$ E29

A convenient set of parameters results as: $p=3$ ; $s=1.5$ ; $a=12$ ; $I=1.5$ . An example of selecting marked objects is shown in Figure 8. The image (a) contains six compact objects and curves. Using the marking points from (b) the objects from (d) were selected, while using points in (e), we selected the objects as in image (f). The image (c) is an intermediate snapshot showing the reconstruction of the selected objects around the marking points.

Figure 8.

(a) Test image; (b) initial state (marking points); (c) transient state; (d) output image.

3.7. Finding the Intersection Points of One-Pixel Thin Lines from Two Binary Images

We consider two binary images $I_1$ and $I_2$ containing one-pixel thin lines or curves, among other compact objects. Superposing the two images, the objects will have overlapping areas. We propose to design the feedback and control templates such that the CNN detects only the intersection points (marked by single black pixels) of the thin lines present in the two images. All other overlapping areas will be ignored and will not appear in the output image. This task may be regarded as logic selective AND operation, since it only applies to some objects of particular shape from the images to be processed. One of the images is loaded as initial state while the other is applied to the input: $X(0)=I_1$ ; $U=I_2$ . Interchanging the two images ( $X(0)=I_2$ ; $U=I_1$ ) yields the same result. The feedback and control templates are identical, $A=B$ , of the form (21). For this task we get the system:

$$2p-8s+I12p-2s+I-114s-I1$$ E30

and a set of suitable parameter values is: $p=3$ ; $s=-0.5$ ; $I=-8.5$ . The described task is shown in Figure 9. The binary images (a) and (b) contain thin lines and various other objects. The first image is loaded as initial state, the second is applied to the input. The output result

is the image (c) which contains only the intersection points of thin lines.

Figure 9.

(a) Image I1 as initial state; (b) Input image I2; (c) Resulted output image.

4. Applications of CNNs in Image Linear Filtering

Although cellular neural networks are basically designed for nonlinear operation, an important class of applications of CNNs in image processing is also linear filtering (Crounse & Chua, 1995). Generally, the two-dimensional spatial filters are most commonly designed in the frequency domain, by imposing a set of specifications, or starting from various 1D prototype filters with given characteristics and applying frequency transformations which finally lead to the desired 2D filter. For each specific method the filter stability must also be ensured. We will approach here only the IIR spatial filters. Generally the existing design methods of 2D IIR filters rely to a large extent on 1D analog filter prototypes, using spectral transformations from s to z plane via bilinear or Euler transformations followed by z to $$(z_1,z_2)$$ transformations. The most investigated were circular, elliptic-shaped, fan filters etc.

When CNNs are used as stable linear filters, the image I to be filtered is usually applied at the input (U=I) and the initial state is usually zero (X=0). The CNN cells must not reach saturation during operation, which implies the restriction for the cell state: $$\left|x_{ij}(t)\right|1$$ , $$i,j=\mathrm{1...}N$$ .

The linear CNN filter is described by the spatial transfer function (Crounse & Chua, 1995):

$$H({\omega }_1,{\omega }_2)=-B({\omega }_1,{\omega }_2)/A({\omega }_1,{\omega }_2)$$ E31

where $$A({\omega }_1,{\omega }_2)$$ , $$B({\omega }_1,{\omega }_2)$$ are the 2D Discrete Space Fourier Transforms of templates A, B.

An essential constraint for CNNs is local connectivity, specific to these parallel systems. Design specifications may lead to high-order 2D filters, requiring large-size templates, that cannot be directly implemented. The templates currently implemented are no larger than $$3\times 3$$ or $$5\times 5$$ . Larger templates can only be implemented by decomposing them into a set of elementary templates. The most convenient are the separable templates, which can be written as a convolution of small-size templates. As we will show, the templates of the 2D filters resulted from 1D prototypes can always be written as convolution products of small-size templates, therefore the filtering can be realized in several steps. The filtering function will be a product of elementary rational functions.

For an 1D recursive CNN filter of order N, its transfer function has the following expression, where equal degrees are taken for numerator and denominator:

$$H(\omega )={\displaystyle \sum _{n=1}^N{b}_n{\cos }^n\omega }/{\displaystyle \sum _{m=1}^N{a}_m{\cos }^m\omega }=B(\omega )/A(\omega )$$ E32

The numerator and denominator of (32) can be factorized into first and second order polynomials in $$\cos \omega$$ . For instance, the numerator can be decomposed as follows:

$$B(\omega )=k\cdot {\displaystyle \prod _{i=1}^n(\cos \omega +b_i)}\cdot {\displaystyle \prod _{j=1}^m({\cos }^2\omega +b_{1j}\cos \omega +b_{2j})}$$ E33

with $$n+2m=N$$ (filter order). Correspondingly, the template B ( $$1\times N$$ ), can be decomposed into $$1\times 3$$ or $$1\times 5$$ templates.

A similar expression is valid for the denominator $$A(\omega )$$ . Coupling conveniently the factors of $$A(\omega )$$ and $$B(\omega )$$ , the filter transfer function (32) can be always written as a product of elementary functions of order 1 or 2, realized with pairs of $$1\times 3$$ or $$1\times 5$$ templates:

$$H(\omega )=H_1(\omega )\cdot H_2(\omega )\cdot \mathrm{...}\cdot H_k(\omega )$$ E34

Here we will discuss zero-phase IIR CNN filters, i.e. with real-valued transfer functions, which correspond to symmetric control and feedback templates. In IIR spatial filtering with a 2D filter $$H({\omega }_1,{\omega }_2)$$ , the final state can be written:

$$\text{X(}{\text{ω}}_{\text{1}}\text{,}{\text{ω}}_{\text{2}}\text{)=U(}{\text{ω}}_{\text{1}}\text{,}{\text{ω}}_{\text{2}}\text{)×H(}{\text{ω}}_{\text{1}}\text{,}{\text{ω}}_{\text{2}}\text{)=U(}{\text{ω}}_{\text{1}}\text{,}{\text{ω}}_{\text{2}}\text{)×}{\text{H}}_{\text{1}}\text{(}{\text{ω}}_{\text{1}}\text{,}{\text{ω}}_{\text{2}}\text{)×}{\text{H}}_{\text{2}}\text{(}{\text{ω}}_{\text{1}}\text{,}{\text{ω}}_{\text{2}}\text{)…}{\text{H}}_{\text{n}}\text{(}{\text{ω}}_{\text{1}}\text{,}{\text{ω}}_{\text{2}}\text{)}$$ E35

The image $$X_1({\omega }_1,{\omega }_2)=U({\omega }_1,{\omega }_2)\cdot H_1({\omega }_1,{\omega }_2)$$ resulted in the first filtering step is re-applied to CNN input, giving the second output image: $$X_2({\omega }_1,{\omega }_2)=X_1({\omega }_1,{\omega }_2)\cdot H_2({\omega }_1,{\omega }_2)$$ and so on, until the whole filtering is achieved. The frequency response $$H({\omega }_1,{\omega }_2)$$ can be written:

$$H({\omega }_1,{\omega }_2)={\displaystyle \prod _k{H}_k({\omega }_1,{\omega }_2)}={\displaystyle \prod _k{B}_k({\omega }_1,{\omega }_2})/{\displaystyle \prod _k{A}_k({\omega }_1,{\omega }_2})$$ E36

which implies a cascade interconnection of IIR filters. The separable templates A and B can be written as convolution products, where $$B_i$$ , $$A_j$$ are elementary ( $$3\times 3$$ or $$5\times 5$$ ) templates:

$$B=B_1\ast B_2\ast \dots \ast B_n A=A_1\ast A_2\ast \dots \ast A_m$$ E37

4.1. Design Methods for Circularly-Symmetric Filters

In this section we propose an efficient design technique for 2D circularly-symmetric filters, based on 1D prototype filters. Circular filters are currently used in many image processing applications. Given a 1D filter function $$H_p(\omega )$$ , the corresponding 2D filter $$H({\omega }_1,{\omega }_2)$$ results through the frequency transformation $$\omega \to \sqrt{{\omega }_1^2+{\omega }_2^2}$$ :

$$H({\omega }_1,{\omega }_2)=H_p\left(\sqrt{{\omega }_1^2+{\omega }_2^2}\right)$$ E38

A currently-used approximation of the function $$\cos \sqrt{{\omega }_1^2+{\omega }_2^2}$$ is given by:

$$\cos \sqrt{{\omega }_1^2+{\omega }_2^2}\cong C({\omega }_1,{\omega }_2)=-0.5+0.5(\cos {\omega }_1+\cos {\omega }_2)+0.5\cos {\omega }_1\cos {\omega }_2$$ E39

and corresponds to the $$3\times 3$$ template:

$$C=\left[\begin{array}{ccc}0.125& 0.25& 0.125\\ 0.25& -0.5& 0.25\\ 0.125& 0.25& 0.125\end{array}\right]$$ E40

Let us consider a 1D filter $$H_P(\omega )$$ described by the symmetric templates $$B_P$$ , $$A_P$$ of radius R:

$$\begin{array}{l}{B}_P=[\begin{array}{ccccccc}\dots & b_2& b_1& b_0& b_1& b_2& \dots \end{array}]\\ A_P=[\begin{array}{ccccccc}\dots & a_2& a_1& a_0& a_1& a_2& \dots \end{array}]\end{array}$$ E41

Applying DSFT and using trigonometric identities we finally get the transfer function as a ratio of polynomials in $$\cos \omega$$ . Its denominator has the factorized form (with, the filter order):

$$A_p(\omega )=c_0+{\displaystyle \sum _{k=1}^R{c}_k(\cos \omega }{)}^k=k\cdot {\displaystyle \prod _{i=1}^n(\cos \omega +a_i)}\cdot {\displaystyle \prod _{j=1}^m({\cos }^2\omega +a_{1j}\cos \omega +a_{2j})}$$ E42

Substituting $$\cos \omega$$ by $$C({\omega }_1,{\omega }_2)$$ given by (39) in the prototype function $$A_P(\omega )$$ , we obtain:

$$A({\omega }_1,{\omega }_2)=A_p\left(\sqrt{{\omega }_1^2+{\omega }_2^2}\right)=c_0+{\displaystyle \sum _{k=1}^R{c}_k\cdot C^k({\omega }_1,{\omega }_2)}$$ E43

Therefore, the design of the 2D circular filter consists in substituting $$\cos \omega$$ with $$C({\omega }_1,{\omega }_2)$$ in each factor of $$H_P(\omega )$$ ; $$A({\omega }_1,{\omega }_2)$$ will have the form, corresponding to (42):

$$A({\omega }_1,{\omega }_2)=k\cdot {\displaystyle \prod _{i=1}^n\left(C({\omega }_1,{\omega }_2)+a_i\right)}\cdot {\displaystyle \prod _{j=1}^m\left(C^2({\omega }_1,{\omega }_2)+a_{1j}\cdot C({\omega }_1,{\omega }_2)+a_{2j}\right)}$$ E44

Since the 1D prototype filter function can be factorized, the 2D circularly-symmetric filter is separable, a major advantage in implementation. The large-size template A corresponding to $$A({\omega }_1,{\omega }_2)$$ results directly decomposed into elementary ( $$3\times 3$$ or $$5\times 5$$ ) templates:

$$A=k\cdot (C_1\ast \mathrm{...}\ast C_i\ast \mathrm{...}{C}_n)\ast (D_1\ast \mathrm{...}\ast D_j\ast \mathrm{...}{D}_m)$$ E45

If we use the $$3\times 3$$ template C given in (40) to approximate the 2D cosine function, each template $$C_i$$ from (45) is obtained from C by adding $$a_i$$ to the central element. Each $$5\times 5$$ template $$D_j$$ results as $$D_j=C\ast C+b_{1j}\cdot C_1+b_{2j}\cdot C_0$$ , where $$C_0$$ is a $$5\times 5$$ zero template, with central element one; $$C_1$$ is a $$5\times 5$$ template obtained by bordering B with zeros.

4.2. Maximally-Flat Zero-Phase Filters

We propose a method of designing 2D spatial filters of different types starting from approximations commonly used in temporal filters. However, we will design spatial filters which will preserve only the magnitude characteristics of their temporal prototypes, while the designed filters will be zero-phase; their frequency response will be real-valued. The design of the desired 2D filters will be achieved in several steps, detailed as follows.

The starting point is a discrete prototype filter of a given type and desired specifications, with a transfer function $H(z)$ . Using this discrete filter we first derive an 1D spatial filter. This 1D filter will inherit only the magnitude characteristics of its temporal counterpart, while its phase will be zero throughout the frequency domain, namely the range $$[-\pi ,\pi ]$$ . We consider the magnitude $$\left|H(\omega )\right|$$ of $$H(z)=H(e^{j\omega })$$ . In order to find a zero-phase 1D filter transfer function which approximates $$\left|H(\omega )\right|$$ , we can make the change of frequency variable:

$$\omega =\arccos x\iff x=\cos \omega$$ E46

Substituting $$\omega$$ in $$\left|H(\omega )\right|$$ by $$\arccos x$$ , we first get a real function $$G(x)$$ in variable x; next we find a convenient rational approximation for $$G(x)$$ . One of the most efficient (best tradeoff between accuracy and approximation order) is the Chebyshev-Padé approximation. The main drawback of this method is that the coefficients of the rational approximation can be only determined numerically using for instance a symbolic computation software. After finding the approximation of $$G(x)$$ as a ratio of polynomials in variable x, we return to the original variable $\omega$ by substituting back $$x=\cos \omega$$ , and finally we get the rational approximation of $$\left|H(\omega )\right|$$ in the variable $$\cos \omega$$ on the range $$\omega \in [-\pi ,\pi ]$$ , similar to (32):

$$\left|H(\omega )\right|\cong {\displaystyle \sum _{n=1}^N{b}_n{\cos }^n\omega }/{\displaystyle \sum _{m=1}^N{a}_m{\cos }^m\omega }=B(\omega )/A(\omega )$$ E47

From the numerator $B(\omega )$ and denominator $A(\omega )$ we directly identify the corresponding control and feedback templates, preferably of equal size. The polynomials $B(\omega )$ and $A(\omega )$ can be factorized according to relation (33).

$$A(\omega )=k\cdot {\displaystyle \prod _{i=1}^n(\cos \omega +a_i)}\cdot {\displaystyle \prod _{j=1}^m({\cos }^2\omega +a_{1j}\cos \omega +a_{2j})}$$ E48

with $$n+2m=N$$ (filter order). A similar expression is valid for $$B(\omega )$$ . In the following sections we will use this technique to design several types of 2D filters: circularly-symmetric filters, elliptically-shaped filters, directional filters etc.

Let us consider as prototype a type 2 Chebyshev low-pass digital filter with the following specifications: order $$N=6$$ , stop-band ripple $$R=40dB$$ , stop-band edge frequency $${\omega }_S=0.3$$ (where generally $${\omega }_S\in [0.0,1.0]$$ , with 1.0 corresponding to half the sample rate).

Following the procedure described above, we finally get an approximation of the form (47). The numerator and denominator coefficients are given by the vectors:

$$\begin{array}{l}[b_n]=[\text{0}\text{.0176 }\text{0}\text{.0068 }\text{-0}\text{.0164}\text{ 0}\text{.0038}]\\ [a_m]=[\text{-0}\text{.4240}\text{ 1}\text{.2552 }\text{-1}\text{.1918}\text{ 0}\text{.3724}]\end{array}$$ E49

and the factorization of the numerator and denominator of the 1D prototype function is :

$$B(\omega )=0.\text{0176}\cdot \text{(cos}\omega +\text{1}\text{.261)}\cdot \text{(cos}\omega \text{-0}\text{.2957)}\cdot \text{(cos}\omega \text{-0}\text{.5789)}$$ E50

$$A(\omega )=\text{-0}\text{.424}\cdot \text{(cos}\omega \text{-1}\text{.4035)}\cdot {\text{(cos}}^{\text{2}}\omega \text{-1}\text{.5568}\cdot \text{cos}\omega +0.\text{6257)}$$ E51

To obtain a circularly-symmetric filter from the 1D factorized prototype function, we simply replace in (50) and (51) $$\cos \omega$$ with the circular cosine function (39). For instance, the feedback template A results as the discrete convolution:

$$A=k\cdot A_{11}\ast A_{12}\ast \mathrm{...}\ast A_{1n}\ast A_{21}\ast A_{22}\ast \mathrm{...}\ast A_{2m}$$ E52

where $$A_{1i}$$ $$(i=\mathrm{1...}n)$$ are $$3\times 3$$ templates and $$A_{2j}$$ $$(j=\mathrm{1...}m)$$ are $$5\times 5$$ templates, given by:

$$A_{1i}=C+a_i\cdot A_{01}$$ E53

$$A_{2j}=C\ast C+a_{1j}\cdot C_0+a_{2j}\cdot A_{02}$$ E54

where $$A_{01}$$ is a $$3\times 3$$ zero template and $$A_{02}$$ is a $$5\times 5$$ zero template with central element one; $$C_0$$ is a $$5\times 5$$ template obtained by bordering $$C(3\times 3)$$ with zeros. The above expressions correspond to the factors in (44). In Figure 10 the frequency response and the contour plot are shown for two 2D circularly-symmetric maximally-flat filters. The first one has a narrower characteristic, corresponding to $${\omega }_S=0.3$$ , while the second has a larger cutoff frequency, being derived from a Chebyshev prototype with $$N=6$$ , $$R=40dB$$ , $${\omega }_S=0.7$$ .

Figure 10.

(a), (c) maximally-flat circular filter frequency responses; (b), (d) contour plots.

4.3. Elliptically-Shaped Oriented Filters

Next we study a class of 2D low-pass filters having an elliptically-shaped horizontal section. These filters will be specified by imposing the values of the semi-axes of the ellipse, and the orientation is given by the angle of the large axis with respect to ${\omega }_2-$ axis in the frequency plane. Starting from the frequency response for a 1D zero-phase filter given by (47), we can derive a 2D elliptically-shaped filter using the frequency transformation $\omega \to \sqrt{E_{\phi }({\omega }_1,{\omega }_2)}$ :

$$\begin{array}{l}{E}_{\phi }({\omega }_1,{\omega }_2)={\omega }_1^2\left(\frac{{\cos }^2\phi }{E^2}+\frac{{\sin }^2\phi }{F^2}\right)+{\omega }_2^2\left(\frac{{\sin }^2\phi }{E^2}+\frac{{\cos }^2\phi }{F^2}\right)\\ +{\omega }_1{\omega }_2\sin (2\phi )\left(\frac{1}{F^2}-\frac{1}{E^2}\right)=a\cdot {\omega }_1^2+b\cdot {\omega }_2^2+c\cdot {\omega }_1{\omega }_2\end{array}$$ E55

An elliptically-shaped filter results from a circular one through the linear transformation:

$$\left[\begin{array}{c}{\omega }_1\\ {\omega }_2\end{array}\right]=\left[\begin{array}{cc}E& 0\\ 0& F\end{array}\right]\cdot \left[\begin{array}{cc}\cos \phi & -\sin \phi \\ \sin \phi & \cos \phi \end{array}\right]\cdot \left[\begin{array}{c}{\omega }_1^{\text{'}}\\ {\omega }_2^{\text{'}}\end{array}\right]$$ E56

where E and F are the stretching coefficients ( $EF$ ), $({\omega }_1,{\omega }_2)$ are the current coordinates and $({\omega }_1^{\text{'}},{\omega }_2^{\text{'}})$ are the former coordinates. Thus, the unit circle is stretched along the axes ${\omega }_1$ and ${\omega }_2$ with factors E and F, then counter-clockwise rotated with an angle $\phi$ , becoming an oriented ellipse, which is the shape of the proposed filter in horizontal section. Consequently, given a 1D prototype filter of the general form (47), we can obtain a 2D elliptically-shaped filter, specified by the parameters E, F and $\phi$ which give the shape and orientation, by using the frequency transformation:

$${\omega }^2\to E_{\phi }({\omega }_1,{\omega }_2)=a\cdot {\omega }_1^2+b\cdot {\omega }_2^2+c\cdot {\omega }_1{\omega }_2$$ E57

where the coefficients a, b and c are identified from relation (55).

Replacing the real frequency variables ${\omega }_1$ and ${\omega }_2$ by the complex variables $s_1=j{\omega }_1$ and $s_2=j{\omega }_2$ , the function $E_{\phi }({\omega }_1,{\omega }_2)$ can be written in the 2D Laplace domain as:

$$E_{\phi }(s_1,s_2)=-(a\cdot s_1^2+b\cdot s_2^2+c\cdot s_1{s}_2)$$ E58

The next step is to find the discrete approximation $E_{\phi }(z_1,z_2)$ for $E_{\phi }(s_1,s_2)$ in (58), using the forward or backward Euler approximations, or the bilinear transform. We will use here the forward Euler method, therefore along the two axes of the plane we have: $s_1=z_1-1$ , $s_2=z_2-1$ . The discrete approximations for $$s_1^2$$ , $$s_2^2$$ and $$s_1{s}_2$$ will result as (Matei, 2004):

$$\begin{array}{l}{s}_1^2=z_1+z_1^{-1}-2\\ s_2^2=z_2+z_2^{-1}-2\\ 2s_1{s}_2=z_1+z_1^{-1}+z_2+z_2^{-1}-2-z_1{z}_2^{-1}-z_1^{-1}{z}_2\end{array}$$ E59

Replacing the above operators in (58), $E_{\phi }(s_1,s_2)$ takes the form:

$$E_{\phi }(z_1,z_2)=-(a+c/2)\cdot (z_1^{-1}+z_1)-(b+c/2)\cdot (z_2^{-1}+z_2)+(c/2)\cdot (z_1{z}_2^{-1}+z_1^{-1}{z}_2)+(2a+2b+c)$$ E60

After re-arranging terms and identifying coefficients of the 2D Z transform, we obtain the matrix operator $E_{\phi }$ :

$$E_{\phi }=\left[\begin{array}{ccc}0& a+c/2& -c/2\\ b+c/2& -(2a+2b+c)& b+c/2\\ -c/2& a+c/2& 0\end{array}\right]$$ E61

The matrix $E_{\phi }$ depends on the orientation angle $\phi$ and the stretching coefficients E and F. We found the frequency transformation $E_{\phi }:ℝ\to {ℂ}^2$ , ${\omega }^2\to E_{\phi }(z_1,z_2)$ in the matrix form:

$${\omega }^2\to \left[\begin{array}{ccc}{z}_1^{-1}& 1& z_1\end{array}\right]\times \left[E_{\phi }\right]\times {\left[\begin{array}{ccc}{z}_2^{-1}& 1& z_2\end{array}\right]}^T$$ E62

As anumerical example, taking the parameter values $\phi =\pi /3$ , $E=2.4$ , $F=0.6$ we get:

$$E_{\phi }=\left[\begin{array}{ccc}0& 3.2544& -1.1277\\ 1.9523& -8.1581& 1.9523\\ -1.1277& 3.2544& 0\end{array}\right]$$ E63

Let us consider the simple 1D zero-phase prototype filter:

$$H_P(\omega )=\left(b_0+b_1{\omega }^2+b_2{\omega }^4\right)/\left(1+a_1{\omega }^2+a_2{\omega }^4\right)$$ E64

Applying the frequency transform (62), we obtain the 2D elliptically-shaped filter:

$$H_{E\phi }(z_1,z_2)=\frac{b_0+b_1\cdot E_{\phi }(z_1,z_2)+b_2\cdot E_{\phi }^2(z_1,z_2)}{1+a_1\cdot E_{\phi }(z_1,z_2)+a_2\cdot E_{\phi }^2(z_1,z_2)}=\frac{B_{E\phi }(z_1,z_2)}{A_{E\phi }(z_1,z_2)}$$ E65

Since $E_{\phi }(z_1,z_2)$ corresponds to the matrix $E_{\phi }$ in (61), we determine the filter templates as:

$$B_{E\phi }=b_0\cdot E_0+b_1\cdot E_{\phi 1}+b_2\cdot E_{\phi }\ast E_{\phi }{A}_{E\phi }=E_0+a_1\cdot E_{\phi 1}+a_2\cdot E_{\phi }\ast E_{\phi }$$ E66

By $\ast$ we denoted matrix convolution. $E_0$ is a $5\times 5$ zero matrix with the central element of value 1. The $5\times 5$ matrix $E_{\phi 1}$ is obtained by bordering with zeros the $3\times 3$ matrix $E_{\phi }$ in order to be summed with $E_0$ and $E_{\phi }\ast E_{\phi }$ .

In general, for a zero-phase prototype $H_P(\omega )$ of order $2N$ , we obtain following the method presented above a 2D filter transfer function $H_{E\phi }(z_1,z_2)$ whose numerator and denominator correspond to templates $B_{E\phi }$ and $A_{E\phi }$ of size $(2N+1)\times (2N+1)$ .

For a design example, we take as 1D prototype a type-2 Chebyshev digital filter with $N=4$ , $R_s=40$ dB and passband-edge frequency ${\omega }_p=0.5$ . Its transfer function in z is:

$$H_p(z)=\frac{0.012277\cdot z^2-0.012525\cdot z+0.012277}{z^2-1.850147\cdot z+0.862316}$$ E67

The real-valued transfer function approximating the magnitude of the filter $H_p(z)$ is:

$$\left|H_p(e^{j\omega })\right|\cong H_{a1}(\omega )=\frac{0.940306-0.57565\cdot {\omega }^2+0.0947\cdot {\omega }^4}{1-2.067753\cdot {\omega }^2+4.663147\cdot {\omega }^4}$$ E68

Figure 11 shows two elliptically-shaped filters based on (68), having a flat characteristic in the pass-band and a very low stop-band ripple. An advantage of this design method is that the filter characteristics in the vertical and horizontal plane are specified independently.

Figure 11.

Frequency response and contour plots of two elliptically-shaped filters: (a), (b) $\phi =\pi /3$ , $E=2.8$ , $F=1$ ; (c), (d) $\phi =\pi /3$ , $E=2.4$ , $F=0.6$

4.4. Directional Filters

There are several classes of 2D filters with orientation-selective frequency response, which find applications in many image processing tasks, such as edge detection, motion analysis, texture analysis etc. For instance, steerable filters which are synthesized as a linear combination of a set of basis filters (Freeman & Adelson, 1991). Another important class are the Gabor filters, efficiently implemented also on CNNs (Shi, 1998), with useful applications in image processing and computer vision. A particular class of filters, which can be also implemented on CNNs, are directional (orientation-selective) filters, which select narrow domains along specified directions in the frequency plane. Several design examples of 2D recursive directional filters with an arbitrary orientation were given in (Matei, 2004). They can be used in selecting lines with a given orientation in the input image.

The spatial orientation is specified by an angle $$\phi$$ with respect to $${\omega }_1$$ -axis and is defined by the 1D to 2D spatial frequency transformation: $$\omega \to {\omega }_1\cos \phi +{\omega }_2\sin \phi$$ . Considering a 1D zero-phase prototype filter $$H(\omega )$$ and using this substitution we obtain the transfer function of the oriented filter, $$H_{\phi }({\omega }_1,{\omega }_2)$$ :

$$H_{\phi }({\omega }_1,{\omega }_2)=H({\omega }_1\cos \phi +{\omega }_2\sin \phi )$$ E69

The 2D oriented filter $$H_{\phi }({\omega }_1,{\omega }_2)$$ has the magnitude along the line $${\omega }_1\cos \phi +{\omega }_2\sin \phi =0$$ , identical with the prototype $$H(\omega )$$ , and is constant along the line $${\omega }_1\sin \phi -{\omega }_2\cos \phi =0$$ (longitudinal axis). The design of an oriented filter based on a 1D prototype relies on finding a realization of the 2D oriented cosine function, which depends on the orientation angle $$\phi$$ :

$$C_{\phi }({\omega }_1,{\omega }_2)=\cos ({\omega }_1\cos \phi +{\omega }_2\sin \phi )$$ E70

Next we will determine a convenient 1D to 2D complex transformation which allows an oriented 2D filter to be obtained from a 1D prototype filter. We will treat here the zero-phase filters, which sometimes are desirable in image processing tasks.

An accurate rational approximation (Chebyshev- Padé) of the cosine for $$\omega \in [-\pi \sqrt{2},\pi \sqrt{2}]$$ is:

$$\cos \omega \cong \frac{1-0.447754{\omega }^2+0.018248{\omega }^4}{1+0.041694{\omega }^2+0.002416{\omega }^4}$$ E71

The range of approximation $$[-\pi \sqrt{2},\pi \sqrt{2}]$$ was required due to the fact that in the frequency plane $$-\pi \le {\omega }_1\le \pi$$ , $$-\pi \le {\omega }_2\le \pi$$ the largest value for the frequency $$\omega$$ of the prototype filter is $$\pi \sqrt{2}$$ . Next we have to find a discrete approximation for the expression (70). Denoting $$f_{\phi }({\omega }_1,{\omega }_2)={({\omega }_1\cos \phi +{\omega }_2\sin \phi )}^2$$ , we obtain using (70) and (71):

$$C_{\phi }({\omega }_1,{\omega }_2)\cong \frac{1-0.447754f_{\phi }({\omega }_1,{\omega }_2)+0.018248f_{\phi }^2({\omega }_1,{\omega }_2)}{1+0.041694f_{\phi }({\omega }_1,{\omega }_2)+0.002416f_{\phi }^2({\omega }_1,{\omega }_2)}$$ E72

$$f_{\phi }({\omega }_1,{\omega }_2)$$ can be written in 2D Laplace domain in the complex variables $$s_1=j{\omega }_1$$ , $$s_2=j{\omega }_2$$ :

$$f_{\phi }(s_1,s_2)=-{(s_1\cos \phi +s_2\sin \phi )}^2=-(s_1^2{\cos }^2\phi +s_2^2{\sin }^2\phi +2s_1{s}_2\sin \phi \cos \phi )$$ E73

Using the forward Euler method as in section 4.3, we obtain the discrete approximations (59) for $$s_1^2$$ , $$s_2^2$$ and $$s_1{s}_2$$ ; substituting these expressions into (73), $$f_{\phi }(s_1,s_2)$$ finally becomes:

$$\begin{array}{l}{f}_{\phi }(z_1,z_2)=-[{\cos }^2\phi (z_1+z_1^{-1}-2)+{\sin }^2\phi (z_2+z_2^{-1}-2)\\ +\sin \phi \cos \phi (z_1+z_1^{-1}+z_2+z_2^{-1}-2-z_1{z}_2^{-1}-z_1^{-1}{z}_2)]\end{array}$$ E74

corresponding to an "orientation matrix" operator $$O_{\phi }$$ , resulted after re-arranging terms and identifying the coefficients of 2D Z transform:

$$O_{\phi }=-\frac{1}{2}\left(\sin 2\phi \cdot \left[\begin{array}{ccc}0& 1& -1\\ 1& -2& 1\\ -1& 1& 0\end{array}\right]+\cos 2\phi \cdot \left[\begin{array}{ccc}\text{}0& 1& \text{}0\\ -1& 0& -1\\ \text{}0& 1& \text{}0\end{array}\right]+\left[\begin{array}{ccc}0& 1& 0\\ 1& -4& 1\\ 0& 1& 0\end{array}\right]\right)$$ E75

The matrix $$O_{\phi }$$ depends only on the orientation angle $$\phi$$ . The 2D rational function $$C_{\phi }({\omega }_1,{\omega }_2)$$ , derived from relations (72) and (74), can be regarded itself as a recursive filter:

$$C_{\phi }({\omega }_1,{\omega }_2)=B_{\phi }({\omega }_1,{\omega }_2)/A_{\phi }({\omega }_1,{\omega }_2)$$ E76

and can be realized using the pair of templates:

$$B_{\phi }=O_0-0.447754\cdot O_{\phi }^{\text{'}}+0.018248\cdot O_{\phi }\ast O_{\phi }$$ E77

$$A_{\phi }=O_0+0.041694\cdot O_{\phi }^{\text{'}}+0.002416\cdot O_{\phi }\ast O_{\phi }$$ E78

also depending only on the orientation angle $$\phi$$ . By " $$\ast$$ " we denoted the operation of matrix convolution; $$O_0$$ is a $$5\times 5$$ zero matrix with the central element of value 1; in (77) and (78) the $$5\times 5$$ matrix $$O_{\phi }^{\text{'}}$$ is derived by bordering with zeros the $$3\times 3$$ matrix $$O_{\phi }$$ in order to be summed with $$O_0$$ and $$O_{\phi }\ast O_{\phi }$$ .

For a design example we use the 1D prototype filter from section 4.2, where $$A(\omega )$$ and $$B(\omega )$$ in the transfer function are factorized as in (50) and (51); we simply have to substitute $$\cos \omega$$ by the oriented cosine function (76). Corresponding to $$A(\omega )$$ factorized as in (51), the resulting function $$A({\omega }_1,{\omega }_2)$$ is realized by a template A, decomposed into elementary templates as in (52), where $$A_{1i}$$ are $$3\times 3$$ templates and $$A_{2j}$$ are $$5\times 5$$ templates, given by:

$$\begin{array}{l}{A}_{1i}=B_{\phi }+a_i\cdot A_{\phi }\\ A_{2j}=B_{\phi }\ast B_{\phi }+a_{1j}\cdot A_{\phi }\ast B_{\phi }+a_{2j}\cdot A_{\phi }\ast A_{\phi }\end{array}$$ E79

Using this method, some design examples are presented as follows.

1) With a zero-phase elliptic filter (Figure 12 (a)), we design a directional filter oriented at an angle $\phi =\pi /6$ ; its frequency response and contour plot are shown in Figure 12 (b) and (c).

2) A simple and efficient low-pass zero-phase filter prototype has the frequency response:

$$H_p(\omega )=1/(p+1-p\cos \omega )$$ E80

where the value of p sets the filter selectivity. If in (80) $$\cos \omega$$ is substituted by $$C_{\phi }({\omega }_1,{\omega }_2)$$ given by (76), for $p=25$ , $\phi =0.22\pi$ a very selective directional filter results (Figure 12 (d), (e)). Another two very selective directional filters are shown in Figure 13; these are in fact elliptically-shaped filters with a large ratio $E/F$ , and are designed as shown in section 4.3; for instance, the first filter given in Figure 13 (a), (b) has the parameters $E=3$ , $F=0.1$ and $\phi =0.22\pi$ .

Directional filters may be used in detecting straight lines with a given orientation in an image (Matei, 2004). In order to demonstrate the directional selectivity of the proposed filters, we present an example of detecting lines with a given inclination from an image, by means of filtering with a 2D IIR oriented low-pass filter. The spectrum of a straight line is oriented in the plane ( $${\omega }_1$$ , $${\omega }_2$$ ) at an angle of $$\pi /2$$ with respect to the line direction. The binary image in Figure 15 (a) contains straight lines with different orientations and some ellipses and is filtered with an oriented elliptically-shaped filter with $$\phi =\pi /6$$ , designed using the method presented in section III. Depending on the filter selectivity, in the filtered image in Figure 15 (b), only the lines which have the spectrum oriented more or less along the filter characteristic, remain practically unchanged, while all the other lines appear more or less blurred, due to low-pass filtering. For a good orientation resolution, to separate lines with very close orientation

Figure 12.

(a) Elliptic 1D LP prototype filter; (b) 2D directional filter response ( $\phi =\pi /6$ ); (c) contour plot; (d) 2D directional frequency response for $\phi =0.22\pi$ ; (e) contour plot

Figure 13.

Directional elliptically-shaped filters: (a),(b): $p=5$ , $E=3$ , $F=0.1$ , $\phi =0.22\pi$ ; (c), (d) $p=5$ , $E=7$ , $F=0.5$ , $\phi =0.22\pi$

Figure 14.

(a) Test image; (b),(c),(d) directionally filtered images with parameters: (a) $\phi =0.8\pi$ , $p=30$ , $E=7$ , $F=0.5$ ; (b) $$\phi =\pi /7$$ , $p=70$ , $E=7$ , $F=0.3$ ; (c) $$\phi =-0.2\pi$$ , $p=50$ , $E=7$ , $F=0.3$

Figure 15.

(a) Binary test image; (b) Directionally filtered image; (c) Texture image;

(d) Directionally filtered image with $$\phi =\pi /12$$ , $p=70$ , $E=7$ , $F=0.3$

angles, the filter selectivity has to be properly chosen. Another example is shown in Figure 15 (c), (d). The texture-type grayscale image (c) representing a portion of a straw-chair upper face is directionally filtered with the given parameter values and the result is image (d), which suggests a possible use in texture segmentation.

4.5. Filters Designed in Polar Coordinates

In this section we study a family of 2D zero-phase filters which can be designed starting from a specified shape in polar coordinates in the frequency plane. The contour plots of their frequency response can be defined as closed curves which can be described in terms of a variable radius $$\rho$$ which is a periodic function of the current angle $$\phi$$ formed with one of the axes. This periodic radius can be expressed using a rational approximation. Then, using some desired 1D prototype filters with factorized transfer functions, spatial filters can be obtained by a 1D to 2D frequency transformation. Their frequency response is symmetric about the origin and has at the same time an angular periodicity. If we consider any contour resulted from the intersection of the frequency response with a horizontal plane, the contour has to be a closed curve which can be described in polar coordinates by: $\rho =\rho (\phi )$ where $$\phi$$ is the angle formed by the radius OP with the ${\omega }_1$ - axis, as shown in Figure 16 (a) for a four-lobe filter. Therefore $\rho (\phi )$ is a periodic function of the angle $$\phi$$ in the range $\phi \in [0,2\pi ]$ .

The proposed design method is based on zero-phase filter prototypes whose transfer function is real-valued and can be expressed as a ratio of polynomials in even powers of the frequency $\omega$ . In general this filter will be described by:

$$H_p(\omega )={\displaystyle \sum _{j=0}^M{b}_j\cdot {\omega }^{2j}}/{\displaystyle \sum _{k=0}^N{a}_k\cdot {\omega }^{2k}}$$ E81

where $M\le N$ and N is the filter order. This function may be obtained through a rational approximation of the magnitude of an usual prototype (Chebyshev, elliptic etc.). The proposed design method for this class of 2D filters is based on a frequency transformation of the form:

$$F:ℝ\to {ℂ}^2, {\omega }^2\to F(z_1,z_2)$$ E82

Through this transformation we will be able to obtain low-pass type filters, symmetric about the origin, as in fact are most 2D filters currently used in image processing. The frequency transformation (82) is a mapping from the real frequency axis $\omega$ to the complex plane $(z_1,z_2)$ . However, it will be defined initially through an intermediate frequency mapping of the form:

$F_1:ℝ\to {ℝ}^2$ ,

$${\omega }^2\to F_1({\omega }_1,{\omega }_2)=({\omega }_1^2+{\omega }_1^2)/\rho ({\omega }_1,{\omega }_2)$$ E83

Here the function $$\rho ({\omega }_1,{\omega }_2)$$ plays the role of a radial compressing function and is initially determined in the angular variable $\phi$ as $\rho (\phi )$ . In the frequency plane $({\omega }_1,{\omega }_2)$ we have:

$$\cos \phi ={\omega }_1/\sqrt{{\omega }_1^2+{\omega }_2^2}$$ E84

If the radial function $\rho (\phi )$ can be expressed in the variable $\cos \phi$ , using (84) we obtain by substitution the function $$\rho ({\omega }_1,{\omega }_2)$$ . Generally the function $\rho (\phi )$ will have to be determined as a polynomial or a ratio of polynomials in variable $\cos \phi$ . For instance, the four-lobe filter whose contour plot is shown in Figure 16 (a) corresponds to a function:

$$\rho (\phi )=a+b\cos 4\phi =a+b-8b{\cos }^2\phi +8b{\cos }^4\phi$$ E85

which is plotted in Figure 16 (b) in the range $\phi \in [0,2\pi ]$ .

As a first 1D prototype we consider a type-2 Chebyshev digital filter with the parameters: order $N=4$ , stopband attenuation $R_S=40$ dB and passband-edge frequency ${\omega }_p=0.7$ , where 1.0 is half the sampling frequency. The transfer function in z for this filter is:

$$H_p(z)=\frac{0.012277\cdot z^2-0.012525\cdot z+0.012277}{z^2-1.850147\cdot z+0.862316}$$ E86

The magnitude of its transfer function in the range $[-\pi ,\pi ]$ is shown in Figure 16 (c). Using a Chebyshev-Padé approximation, we determine the real-valued zero-phase frequency response which approximates quite accurately the magnitude of the original filter function:

$$\left|H_p(e^{j\omega })\right|\cong H_{a1}(\omega )=\frac{0.940306-0.57565\cdot {\omega }^2+0.0947\cdot {\omega }^4}{1-2.067753\cdot {\omega }^2+4.663147\cdot {\omega }^4}$$ E87

Let us now consider the radial function with the form:

$$H_r(\phi )=1/\left(p\cdot \tilde{B}(\phi )-p+1\right)$$ E88

where $\tilde{B}(\phi )$ is a periodic function; let $\tilde{B}(\omega )=\cos (4\phi )$ . We will use this function to design a 2D filter with four narrow lobes in the frequency plane $({\omega }_1,{\omega }_2)$ . Using trigonometric identities, it is straightforward to obtain:

$$H_r(\phi )=1/\left(1+8p\cdot {(\cos \phi )}^2-8p\cdot {(\cos \phi )}^4\right)$$ E89

plotted for $\phi \in [-\pi ,\pi ]$ in Figure 16 (d). This is a periodic function with period $\Phi =\pi /4$ and has the shape of a “comb” filter. In order to control the shape of this function, we introduce a parameter k, such that the radial function will be $\rho (\phi )=k\cdot H_r(\phi )$ . We obtain using (83):

$${\omega }^2\to F({\omega }_1,{\omega }_2)=\frac{{\omega }_1^4+(2+8p){\omega }_1^2{\omega }_2^2+{\omega }_2^4}{k({\omega }_1^2+{\omega }_2^2)}$$ E90

Since ${\omega }_1^2=-s_1^2$ and ${\omega }_2^2=-s_2^2$ we derive the function $$F_1(s_1,s_2)$$ in the plane $$(s_1,s_2)$$ as:

$$F_2(s_1,s_2)=-\left(s_1^4+(2+8p)s_1^2{s}_2^2+s_2^4\right)/\left(k(s_1^2+s_2^2)\right)$$ E91

Finally we determine a transfer function of the 2D filter $H(z_1,z_2)$ in the complex plane $(z_1,z_2)$ . This can be achieved if we find a discrete counterpart of the function $$\rho ({\omega }_1,{\omega }_2)$$ .

A possible method is to express the function $$\rho ({\omega }_1,{\omega }_2)$$ in the complex plane $(s_1,s_2)$ and then to find the appropriate mapping to the plane $(z_1,z_2)$ . This can be achieved either using the forward or backward Euler approximations, or otherwise the bilinear transform, which gives better accuracy. The bilinear transform is a first-order approximation of the natural logarithm

Figure 16.

(a) Contour plot of a four-lobe filter; (b) Variation of the periodic function $$\rho (\phi )$$ ; (c) Maximally-flat low-pass prototype; (d) Very selective radial function

function, which is an exact mapping of the z-plane to the s-plane. For our purposes the sample interval takes the value $T=1$ so the bilinear transform for the two variables $s_1$ and $s_2$ in the complex plane $(s_1,s_2)$ has the form:

$$s_1=2(z_1-1)/(z_1+1) s_2=2(z_2-1)/(z_2+1)$$ E92

Substituting $s_1$ , $s_2$ in (91), we find the frequency transformation in $z_1$ , $z_2$ , in matrix form:

$${\omega }^2\to F(z_1,z_2)=\frac{B(z_1,z_2)}{A(z_1,z_2)}=\frac{Z_1\times \left[B\right]\times Z_2^T}{Z_1\times \left[A\right]\times Z_2^T}$$ E93

where $Z_1$ and $Z_2$ are the vectors:

$$Z_1=\left[\begin{array}{ccccc}{z}_1^{-2}& z_1^{-1}& 1& z_1& z_1^2\end{array}\right]Z_2=\left[\begin{array}{ccccc}{z}_2^{-2}& z_2^{-1}& 1& z_2& z_2^2\end{array}\right]$$ E94

where $\times$ denotes matrix/vector product. The expression of the frequency transformation (93) can be considered in itself a 2D recursive filter. The templates B, A giving the coefficients of $B(z_1,z_2)$ , $A(z_1,z_2)$ are the matrices of size $5\times 5$ :

$$B=8\cdot \left[\begin{array}{ccccc}-1-2p& 0& 4p-2& 0& -1-2p\\ 0& 8& 0& 8& 0\\ 4p-2& 0& -8p-20& 0& 4p-2\\ 0& 8& 0& 8& 0\\ -1-2p& 0& 4p-2& 0& -1-2p\end{array}\right] A=k\cdot \left[\begin{array}{ccc}1& 0& 1\\ 0& -4& 0\\ 1& 0& 1\end{array}\right]\ast \left[\begin{array}{ccc}1& 2& 1\\ 2& 4& 2\\ 1& 2& 1\end{array}\right]$$ E95

We next use the maximally flat filter prototype (87). We only have to make the substitution of (93) in $H_{a1}(\omega )$ given in (87) and we obtain the desired 2D transfer function:

$$H(z_1,z_2)=\frac{b_2{B}^2(z_1,z_2)+b_{1}A(z_1,z_2)B(z_1,z_2)+b_0{A}^2(z_1,z_2)}{B^2(z_1,z_2)+a_{1}A(z_1,z_2)B(z_1,z_2)+a_0{A}^2(z_1,z_2)}=B_f(z_1,z_2)/A_f(z_1,z_2)$$ E96

where $b_2=0.9403$ , $b_1=-0.5756$ , $b_0=0.0947$ , $a_1=-2.0677$ , $a_0=4.6631$ are the coefficients of the prototype (87). In our case $H(z_1,z_2)$ results of order 8. For a chosen prototype of higher order, we get a similar rational function in powers of $A(z_1,z_2)$ and $B(z_1,z_2)$ . Since the 2D transfer function (96) can be also described in terms of templates $B_f$ and $A_f$ corresponding to $B_f(z_1,z_2)$ , $A_f(z_1,z_2)$ we have equivalently:

$$\begin{array}{l}{B}_f=b_2\cdot B\ast B+b_1\cdot A\ast B+b_0\cdot A\ast A\\ A_f=B\ast B+a_1\cdot A\ast B+a_0\cdot A\ast A\end{array}$$ E97

where $\ast$ denotes two-dimensional convolution. For our filter $B_f$ and $A_f$ are of size $9\times 9$ .

Figure 17.

(a), (b) Frequency response and contour plot for the 4-lobe filter; (c) Test image; (d) Filtered image.

The designed filter has the frequency response and contour plot as in Figure 17 (a), (b). We remark that the filter is very selective simultaneously along both axes of the plane $({\omega }_1,{\omega }_2)$ .

We now present an example of image filtering using a selective four-lobe filter. This type of filter can be used in detecting lines with a given inclination from an image, as the directional filters discussed in section 4.4. The spectrum of a straight line is oriented in the plane $({\omega }_1,{\omega }_2)$ at an angle $$\pi /2$$ with respect to the line direction. The test image in Figure 17 (c) contains straight lines with different orientations and lengths, and also a few curves. We consider the filter shown in Figure 17 (a), with parameter values $$p=50$$ and $$k=0.8$$ . Depending on filter selectivity, only the lines with the spectrum oriented more or less along the filter pass-bands will be preserved in the filtered image. We notice that in the output image shown in Fig.17 (d), the lines roughly oriented horizontally and vertically are preserved, while the others are filtered out or appear very blurred. The intersections of the detected lines appear as darker pixels; this allows also these intersections to be detected, if after filtering a proper threshold is applied.

5. Conclusion

An alternative CNN model was proposed, in which the state equations are modified such that the piece-wise linear function operates on the entire information fed to the current cell, including input, feedback and offset term. The circuit structure, dynamic behavior and equilibrium points were analyzed. Through a linear transformation, the proposed model is shown to be equivalent to the standard CNN model. We aimed at making a connection between these two models and compare them from the recurrent neural network point of view. Certain advantages regarding implementation result from the limited voltage range. Several processing tasks on binary images were proposed, templates were designed and simulation results on test images are provided. The template design is relatively easy and can be optimized in order to obtain high robustness to parameter variations. A deeper insight into the proposed CNN dynamic behavior may lead to further image processing capabilities.

As regards the applications of CNNs in linear filtering, some original and efficient design methods were approached, for several types of two-dimensional filters: circular and elliptically-shaped filters, maximally-flat, orientation-selective filters etc. The design starts from an imposed 1D prototype filter and uses rational approximations and various frequency transformations to obtain the 2D frequency response. These analytical methods could be combined with numerical optimization techniques in order to realize an efficient design.

Even if theoretically all the presented spatial linear filters can be implemented in analog processing systems like CNNs, we have to take into account the fact that these analog filters may be very sensitive to template elements variations due to technological dispersion which is inherently present in any VLSI process. A systematic theoretical study on the effect of these dispersions on filter characteristics and stability can be found in (David, 2004). It has been found both analytically and through simulations that the sensitivity of the frequency response to process variations increases with filter selectivity. Since we have designed also some very selective filters, a sensitivity analysis has to be carried out in each case, to find out if a given filter can be reliably implemented. The filter robustness and stability requirements may imply strong constraints on the realization of filters with high selectivity. Further research also envisages an efficient and reliable implementation of this class of filters.

Acknowledgments

This work was supported by the National University Research Council under Grant PN2 – ID_310 “Algorithms and Parallel Architectures for Signal Acquisition, Compression and Processing”.