Open access peer-reviewed chapter

New Model and Applications of Cellular Neural Networks in Image Processing

By Radu Matei

Published: October 1st 2009

DOI: 10.5772/8223

Downloaded: 1547

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):


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


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:


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 (xj), inputs (uk) and the bias term (Ii). 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 byzithe argument of the piece-wise linear functionfin (3):


where the inputs and bias are constant and were put together into the termki.

Multiplying both members of (3) byAj, summing overiNi(Nidenotes the current cell neighborhood) and shifting again the indices i and j, we finally obtain (Matei, 2001):


having a form similar to (1), but in the new variablezi. 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:


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


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 variabley(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 variablexi, which is a voltage, and gives the output voltageyi=f(xi). In the modified CNN equation (3), the non-linearity applies to the amount denotedzij, summing up the entire information from the states and the inputs of the neighboring cells. The variablezijis a current (given by several linear controlled current sources), andf(zi)is a current as well. For the piecewise-linear current-controlled current source we can writeIyx=f(zi). The resistor denotedRyhas 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 resistorRx. On the cell circuit, considering normalized values for C andRxwe 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,xi[1,1]. This domain defines the linear region of the state space. For the modified CNN, the linear region of the state space in variablexidepends 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 havef(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 derivativedxi/dtvs. state variablexi(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 variablezi, we can first represent the driving-point plot (dzi/dt,zi) in the new variablezias 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 variablexi, we use the linear mapping (4) in the form:


In (8),p=A0is the central element of the feedback template. Supposing for the moment a fixed, constant value forgi(t), denotedgi, we obtain forp0:


The casep=0should be analyzed separately. Through a simple graphical construction shown in Figure 2 (a) we can draw the corresponding driving point plot in the original variablexi. This can be done using two auxiliary plots: (a) the plot(xi,zi), given by (10), represented by a straight line with slope of value1/pand offset valuegi; (b) the plot(dzi/dt,dxi/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 valuesp=2,gi=0.5.

The segmentA1B1with slope (p1), given by:


has an identical form inxi:


The pointsA1,B1, given by the coordinates:A1(1,p1+gi),B1(1,p+1+gi)become:


The linesA1D1andB1C1, described by the equations:

dzi/dt=zi+p+gi        dzi/dt=zip+giE15

becomeADandBCrespectively, given by:

dxi/dt=xi+1               dxi/dt=xi1E16

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(dxi/dt,xi)the points C and D result from the intersection of the straight lines BC and AD with axisxi, imposing the condition:dxi/dt=0.

We finally obtain:xi(C)=1andxi(D)=1. PointsC(1,0)andD(1,0)will therefore be fixed equilibrium points of the modified CNN, independent of the values of p andgi. For different values of the amountgi, the diagram(dzi/dt,zi)shifts up or down. Correspondingly, the(dxi/dt,xi)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 andgiare 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 statesxiin a CNN are bounded fort0and the boundvmaxis given by:


if we consider the normalized resistance valueRx=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 isxi=f(zi)and consequently its variation is bounded to the range[1,1]. Moreover, the output variable is identical to the state variable at any time:yi(t)=xi(t).

Referring again to the circuit structure, the argumentzi(t)of the nonlinear functionf(zi)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 variablezi, equation (17) may be formally applied to evaluate a maximum boundImaxfor 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 ifu,yBMN,x(0)BMN, whereB={-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 cellC(i,j)we will have:x˙ij(0)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:


Suppose that a given pixel is initially black, i.e.xi(0)=1. If it is to remain black, this impliesx˙i(0)=0and from (18) we must have:f(zi(0))=1, thereforezi(0)1. If the black pixel must switch to white, we impose:zi(0)1, sof(zi(0))=1; the initial state derivative is:


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

xi()=1          x˙i()=xi()1=0E20

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:


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 a3×3frame. 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:


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 assumings0; writing the local rules, we reach the following system:


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=[sssspssss];  B=[rrrrprrrr];  r=0.5sE24

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


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 its3×3neighborhood. 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=[rsrspsrsr];   B=[0000a0000]E26

and the following system results:


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π/4is 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=[rsrspsrsr];   B=[0.5r0.5s0.5r0.5sp0.5s0.5r0.5s0.5r]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 haveU=Pand the initial stateX(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. Withs0,p0, we finally obtain the system:


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 imagesI1andI2containing 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)=I1;U=I2. Interchanging the two images (X(0)=I2;U=I1) yields the same result. The feedback and control templates are identical,A=B, of the form (21). For this task we get the system:


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(z1,z2)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:|xij(t)|1,i,j=1...N.

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


whereA(ω1,ω2),B(ω1,ω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 than3×3or5×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:


The numerator and denominator of (32) can be factorized into first and second order polynomials incosω. For instance, the numerator can be decomposed as follows:


withn+2m=N(filter order). Correspondingly, the template B (1×N), can be decomposed into1×3or1×5templates.

A similar expression is valid for the denominatorA(ω). Coupling conveniently the factors ofA(ω)andB(ω), the filter transfer function (32) can be always written as a product of elementary functions of order 1 or 2, realized with pairs of1×3or1×5templates:


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 filterH(ω1,ω2), the final state can be written:


The imageX1(ω1,ω2)=U(ω1,ω2)H1(ω1,ω2)resulted in the first filtering step is re-applied to CNN input, giving the second output image:X2(ω1,ω2)=X1(ω1,ω2)H2(ω1,ω2)and so on, until the whole filtering is achieved. The frequency responseH(ω1,ω2)can be written:


which implies a cascade interconnection of IIR filters. The separable templates A and B can be written as convolution products, whereBi,Ajare elementary (3×3or5×5) templates:

B=B1B2Bn          A=A1A2AmE37

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 functionHp(ω), the corresponding 2D filterH(ω1,ω2)results through the frequency transformationωω12+ω22:


A currently-used approximation of the functioncosω12+ω22is given by:


and corresponds to the3×3template:


Let us consider a 1D filterHP(ω)described by the symmetric templatesBP,APof radius R:


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


SubstitutingcosωbyC(ω1,ω2)given by (39) in the prototype functionAP(ω), we obtain:


Therefore, the design of the 2D circular filter consists in substitutingcosωwithC(ω1,ω2)in each factor ofHP(ω);A(ω1,ω2)will have the form, corresponding to (42):


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 toA(ω1,ω2)results directly decomposed into elementary (3×3or5×5) templates:


If we use the3×3template C given in (40) to approximate the 2D cosine function, each templateCifrom (45) is obtained from C by addingaito the central element. Each5×5templateDjresults asDj=CC+b1jC1+b2jC0, whereC0is a5×5zero template, with central element one;C1is a5×5template 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 functionH(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[π,π]. We consider the magnitude|H(ω)|ofH(z)=H(ejω). In order to find a zero-phase 1D filter transfer function which approximates|H(ω)|, we can make the change of frequency variable:


Substitutingωin|H(ω)|byarccosx, we first get a real functionG(x)in variable x; next we find a convenient rational approximation forG(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 ofG(x)as a ratio of polynomials in variable x, we return to the original variableωby substituting backx=cosω, and finally we get the rational approximation of|H(ω)|in the variablecosωon the rangeω[π,π], similar to (32):


From the numeratorB(ω)and denominatorA(ω)we directly identify the corresponding control and feedback templates, preferably of equal size. The polynomialsB(ω)andA(ω)can be factorized according to relation (33).


withn+2m=N(filter order). A similar expression is valid forB(ω). 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: orderN=6, stop-band rippleR=40dB, stop-band edge frequencyωS=0.3(where generallyωS[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:

[bn]=[0.0176 0.0068 -0.0164 0.0038][am]=[-0.4240 1.2552 -1.1918 0.3724]E49

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


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


whereA1i(i=1...n)are3×3templates andA2j(j=1...m)are5×5templates, given by:


whereA01is a3×3zero template andA02is a5×5zero template with central element one;C0is a5×5template obtained by borderingC(3×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ωS=0.3, while the second has a larger cutoff frequency, being derived from a Chebyshev prototype withN=6,R=40dB,ω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ω2axis 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ωEφ(ω1,ω2):


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


where E and F are the stretching coefficients (EF),(ω1,ω2)are the current coordinates and(ω1',ω2')are the former coordinates. Thus, the unit circle is stretched along the axesω1andω2with factors E and F, then counter-clockwise rotated with an angleφ, 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φwhich give the shape and orientation, by using the frequency transformation:


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

Replacing the real frequency variablesω1andω2by the complex variabless1=jω1ands2=jω2, the functionEφ(ω1,ω2)can be written in the 2D Laplace domain as:


The next step is to find the discrete approximationEφ(z1,z2)forEφ(s1,s2)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:s1=z11,s2=z21. The discrete approximations fors12,s22ands1s2will result as (Matei, 2004):

s12=z1+z112s22=z2+z212          2s1s2=z1+z11+z2+z212z1z21z11z2E59

Replacing the above operators in (58),Eφ(s1,s2)takes the form:


After re-arranging terms and identifying coefficients of the 2D Z transform, we obtain the matrix operatorEφ:


The matrixEφdepends on the orientation angleφand the stretching coefficients E and F. We found the frequency transformationEφ:2,ω2Eφ(z1,z2)in the matrix form:


As anumerical example, taking the parameter valuesφ=π/3,E=2.4,F=0.6we get:


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


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


SinceEφ(z1,z2)corresponds to the matrixEφin (61), we determine the filter templates as:

BEφ=b0E0+b1Eφ1+b2EφEφ           AEφ=E0+a1Eφ1+a2EφEφE66
Bywe denoted matrix convolution.E0is a5×5zero matrix with the central element of value 1. The5×5matrixEφ1is obtained by bordering with zeros the3×3matrixEφin order to be summed withE0andEφEφ.

In general, for a zero-phase prototypeHP(ω)of order2N, we obtain following the method presented above a 2D filter transfer functionHEφ(z1,z2)whose numerator and denominator correspond to templatesBEφandAEφof size(2N+1)×(2N+1).

For a design example, we take as 1D prototype a type-2 Chebyshev digital filter withN=4,Rs=40dB and passband-edge frequencyωp=0.5. Its transfer function in z is:


The real-valued transfer function approximating the magnitude of the filterHp(z)is:


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) φ = π / 3 , E = 2.8 , F = 1 ; (c), (d) φ = π / 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φwith respect toω1-axis and is defined by the 1D to 2D spatial frequency transformation:ωω1cosφ+ω2sinφ. Considering a 1D zero-phase prototype filterH(ω)and using this substitution we obtain the transfer function of the oriented filter,Hφ(ω1,ω2):


The 2D oriented filterHφ(ω1,ω2)has the magnitude along the lineω1cosφ+ω2sinφ=0, identical with the prototypeH(ω), and is constant along the lineω1sinφω2cosφ=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φ:


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ω[π2,π2]is:


The range of approximation[π2,π2]was required due to the fact that in the frequency planeπω1π,πω2πthe largest value for the frequencyωof the prototype filter isπ2. Next we have to find a discrete approximation for the expression (70). Denotingfφ(ω1,ω2)=(ω1cosφ+ω2sinφ)2, we obtain using (70) and (71):

fφ(ω1,ω2)can be written in 2D Laplace domain in the complex variabless1=jω1,s2=jω2:

Using the forward Euler method as in section 4.3, we obtain the discrete approximations (59) fors12,s22ands1s2; substituting these expressions into (73),fφ(s1,s2)finally becomes:


corresponding to an "orientation matrix" operatorOφ, resulted after re-arranging terms and identifying the coefficients of 2D Z transform:


The matrixOφdepends only on the orientation angleφ. The 2D rational functionCφ(ω1,ω2), derived from relations (72) and (74), can be regarded itself as a recursive filter:


and can be realized using the pair of templates:


also depending only on the orientation angleφ. By "" we denoted the operation of matrix convolution;O0is a5×5zero matrix with the central element of value 1; in (77) and (78) the5×5matrixOφ'is derived by bordering with zeros the3×3matrixOφin order to be summed withO0andOφOφ.

For a design example we use the 1D prototype filter from section 4.2, whereA(ω)andB(ω)in the transfer function are factorized as in (50) and (51); we simply have to substitutecosωby the oriented cosine function (76). Corresponding toA(ω)factorized as in (51), the resulting functionA(ω1,ω2)is realized by a template A, decomposed into elementary templates as in (52), whereA1iare3×3templates andA2jare5×5templates, given by:


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φ=π/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:


where the value of p sets the filter selectivity. If in (80)cosωis substituted byCφ(ω1,ω2)given by (76), forp=25,φ=0.22π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 ratioE/F, and are designed as shown in section 4.3; for instance, the first filter given in Figure 13 (a), (b) has the parametersE=3,F=0.1andφ=0.22π.

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 (ω1,ω2) at an angle ofπ/2with 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φ=π/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 ( φ = π / 6 ); (c) contour plot; (d) 2D directional frequency response for φ = 0.22 π ; (e) contour plot

Figure 13.

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

Figure 14.

(a) Test image; (b),(c),(d) directionally filtered images with parameters: (a) φ = 0.8 π , p = 30 , E = 7 , F = 0.5 ; (b) φ = π / 7 , p = 70 , E = 7 , F = 0.3 ; (c) φ = − 0.2 π , p = 50 , E = 7 , F = 0.3

Figure 15.

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

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ρwhich is a periodic function of the current angleφ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:ρ=ρ(φ)whereφis the angle formed by the radius OP with theω1- axis, as shown in Figure 16 (a) for a four-lobe filter. Thereforeρ(φ)is a periodic function of the angleφin the rangeφ[0,2π].

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ω. In general this filter will be described by:


whereMNand 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:2,  ω2F(z1,z2)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ωto the complex plane(z1,z2). However, it will be defined initially through an intermediate frequency mapping of the form:


Here the functionρ(ω1,ω2)plays the role of a radial compressing function and is initially determined in the angular variableφasρ(φ). In the frequency plane(ω1,ω2)we have:


If the radial functionρ(φ)can be expressed in the variablecosφ, using (84) we obtain by substitution the functionρ(ω1,ω2). Generally the functionρ(φ)will have to be determined as a polynomial or a ratio of polynomials in variablecosφ. For instance, the four-lobe filter whose contour plot is shown in Figure 16 (a) corresponds to a function:


which is plotted in Figure 16 (b) in the rangeφ[0,2π].

As a first 1D prototype we consider a type-2 Chebyshev digital filter with the parameters: orderN=4, stopband attenuationRS=40dB and passband-edge frequencyωp=0.7, where 1.0 is half the sampling frequency. The transfer function in z for this filter is:


The magnitude of its transfer function in the range[π,π]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:


Let us now consider the radial function with the form:


whereB˜(φ)is a periodic function; letB˜(ω)=cos(4φ). We will use this function to design a 2D filter with four narrow lobes in the frequency plane(ω1,ω2). Using trigonometric identities, it is straightforward to obtain:


plotted forφ[π,π]in Figure 16 (d). This is a periodic function with periodΦ=π/4and 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ρ(φ)=kHr(φ). We obtain using (83):


Sinceω12=s12andω22=s22we derive the functionF1(s1,s2)in the plane(s1,s2)as:


Finally we determine a transfer function of the 2D filterH(z1,z2)in the complex plane(z1,z2). This can be achieved if we find a discrete counterpart of the functionρ(ω1,ω2).

A possible method is to express the functionρ(ω1,ω2)in the complex plane(s1,s2)and then to find the appropriate mapping to the plane(z1,z2). 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 ρ ( φ ) ; (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 valueT=1so the bilinear transform for the two variabless1ands2in the complex plane(s1,s2)has the form:

s1=2(z11)/(z1+1)        s2=2(z21)/(z2+1)E92

Substitutings1,s2in (91), we find the frequency transformation inz1,z2, in matrix form:


whereZ1andZ2are the vectors:


where×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 ofB(z1,z2),A(z1,z2)are the matrices of size5×5:

B=8[12p04p2012p080804p208p2004p20808012p04p2012p]       A=k[101040101][121242121]E95

We next use the maximally flat filter prototype (87). We only have to make the substitution of (93) inHa1(ω)given in (87) and we obtain the desired 2D transfer function:


whereb2=0.9403,b1=0.5756,b0=0.0947,a1=2.0677,a0=4.6631are the coefficients of the prototype (87). In our caseH(z1,z2)results of order 8. For a chosen prototype of higher order, we get a similar rational function in powers ofA(z1,z2)andB(z1,z2). Since the 2D transfer function (96) can be also described in terms of templatesBfandAfcorresponding toBf(z1,z2),Af(z1,z2)we have equivalently:


wheredenotes two-dimensional convolution. For our filterBfandAfare of size9×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(ω1,ω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(ω1,ω2)at an angleπ/2with 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 valuesp=50andk=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.


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”.

© 2009 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike-3.0 License, which permits use, distribution and reproduction for non-commercial purposes, provided the original is properly cited and derivative works building on this content are distributed under the same license.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Radu Matei (October 1st 2009). New Model and Applications of Cellular Neural Networks in Image Processing, Advanced Technologies, Kankesu Jayanthakumaran, IntechOpen, DOI: 10.5772/8223. Available from:

chapter statistics

1547total chapter downloads

1Crossref citations

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Concentration of Heterogeneous Road Traffic

By Thamizh Arasan Venkatachalam and Dhivya Gnanavelu

Related Book

First chapter

A Survey of Decentralized Adaptive Control

By Karel Perutka

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More about us