INTERNATIONAL JOURNAL OF CIRCUIT THEORY AND APPLICATIONS Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) CIRCUIT REALIZATION OF MARKOV RANDOM FIELDS FOR ANALOG IMAGE PROCESSING MAURO PARODI,* MARCO STORACE AND CARLO REGAZZONI Biophysical and Electronic Engineering Department, University of Genoa, Via Opera Pia 11a, I-16145 Genova, Italy SUMMARY A methodological approach to the definition of nonlinear circuits for real-time image processing is presented. The image processing problem is formulated in terms of the minimization of a functional based on the Markov random fields (MRFs) theory. The terms of such a functional are related to the co-contents of proper nonlinear multiterminal resistors, thus reporting the minimization process to the achievement of an equilibrium solution in a circuit made up of these multiterminal resistors and of linear capacitors. Coupled image restoration and edge extraction in the presence of additive Gaussian noise are the specific problems addressed in this paper. Simulations shap the ability of the proposed circuits to yield real-time solutions for these problems. ( 1998 John Wiley & Sons, Ltd. KEY WORDS: nonlinear circuits; minimization; parallel processing 1. INTRODUCTION A wide number of image processing methods developed in the last decades is based on the minimization of nonlinear functionals. The digital nature of the computational resources oriented the research in image processing to conceive algorithms performing sequential operations in a discrete state space. Frequent drawbacks inherent to such a traditional paradigm concern the need of a compromise between computation times and accuracy of the solutions. On the other hand, these drawbacks are compensated by the algorithms’ robustness. One of the emerging ideas for performing low-level vision tasks in applications with tight real-time requirements is to obtain analog image processing by dedicated nonlinear analog circuits.1 In this case, the traditional paradigm is replaced with a properly defined dynamical system operating in a continuous space. Significant examples of this image processing strategy are found in the literature concerning the cellular neural networks (CNN).2—6 In particular, in Reference 7 a CNN framework for various adaptive nonlinear filters is developed, in order to show that a wide class of image processing problems can be efficiently solved by simple analogic algorithms in the environment of the CNN Universal Machine.8 The theoretical efforts regarding these topics are also justified by the advances in VLSI technology, that allow the integration of complex circuit architectures on a single chip. This paper deals with the synthesis of nonlinear circuits solving a particular class of image processing problems. As shown in Reference 9, many ill-posed early vision problems (e.g. optical flow computation, edge detection) can be solved by minimizing functionals H which are a weighted sum of two terms, i.e., H"H #jH , according to the regularization theory. 1 2 The term H , called the regularization term, restricts the class of admissible solutions by introducing 1 suitable a priori knowledge. *Correspondence to: M. Parodi, Biophysical and Electronic Engineering Department, University of Genoa, Via Opera Pia 11a, I-16145 Genova, Italy. Email: parodi@dibe.unige.it CCC 0098—9886/98/050477—22$17.50 ( 1998 John Wiley & Sons, Ltd. Received 4 February 1997 Revised 2 April 1998 478 M. PARODI, M. STORACE AND C. REGAZZONI The term H , called the data term, constraints the solution to remain close to the observed data. 2 The weighting coefficient j controls the compromise between the degree of regularization of the solution and its closeness to the data. The literature numbers some circuits that minimize, in an analogic and parallel way, functionals structured as H. The general aspects concerning analogue networks for solving ill-posed variational problems in early vision were extensively discussed in References 9, 10 starting from Maxwell’s minimum heat theorem. The resistive part of the pertinent circuits is made up of linear resistors, together with proper voltage and/or current sources. When capacitances are added, the steady-state current or voltage distribution in the network corresponds to the solution of the regularization problem. The functionals considered here are built up on the basis of the Markov random fields (MRFs) theory. In these functionals, the regularization term takes into account the statistical properties of the solution space and the term jH is formulated starting from the probability distribution of the noise affecting the data.9 2 The specific class of MRF image processing models addressed in this paper is oriented to coupled image restoration and edge detection. A detailed description of the general ideas (Gibbs distributions and Bayesian restoration of images) concerning the structure of the corresponding functionals can be found as a part of a rather general paper.11 In such a paper, a stochastic relaxation method is adopted to minimize the functional H , which depends on two sets of variables. Simplified methods for the functional minimization 1 are based on the so-called deterministic relaxation and require less computational time than the stochastic relaxation. On this basis, Sziranyi and Zerubia recently implemented a probabilistic based MRF approach in a CNN Universal Machine.12 Another deterministic approach based on the mean field approximation is proposed in Reference 13. A further method for simplifying the structure of the functional considered in Reference 11 and for minimizing it according to deterministic criteria is described in Reference 14. In such a Reference, that constitutes the starting point of our paper, a method is considered which allows to reconduce the original functional H to a functional that depends on one set of variables only. Starting from 1 this basis, we propose a method to synthesize a nonlinear circuit solving the minimization problem. Each term of the corresponding functional is first put into relation with a potential function associated to a proper nonlinear multiterminal resistor. Then, it is shown that the minimum of the specific MRF functional can be interpreted as an equilibrium state of a circuit made up of both the multiterminal resistors and linear capacitors. A brief outline of the general elements of the considered MRF model is given in the next section. Then, the general structure of the nonlinear multiterminal resistors and of the complete circuit is obtained in Section 3. Sections 4—6 deal with specific aspects concerning the relations between the multiterminal resistors characteristics and the regularization and data terms. In particular, it is shown that the building block for the multiterminal resistor associated to the regularization term is a nonlinear two-terminal element, whose characteristic can be obtained directly from the functional by applying a well-defined method. As a support to the theoretical discussion, some examples of image processing obtained by simulating the circuit behaviour are offered. The main dynamic properties of the circuit are briefly investigated in Section 7. 2. THE CONSIDERED MRF MODEL A well-known method from the image processing theory enables to perform simultaneously image recovery and boundary detection through a coupled MRF model.11,14 2.1. Basic elements In most MRF models, the image to be processed is represented through a set of M]N pixels organized over a regular lattice S"M( j, k): j"1, 2 , M; k"1, 2 , NN. Each pixel jk of S is associated with proper random variables (e.g. the grey level x ). As a general assumption, a random variable at the jk site directly interacts jk with a limited number of neighbour sites. For each pixel jk, the neighbourhood set can be defined as N "Mpq : E( j, k)!(p, q)E2)rN, where r is an integer number representing the order of the neighbourhood set. jk ( 1998 John Wiley & Sons, Ltd. Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) 479 MARKOV RANDOM FIELDS Each kind of variables individuates a field over S. The coupled MRF model considered in this paper needs the definition of two fields. The first (field X) collects the grey levels x . The second (field B) collects functions jk b of the gradient terms (x !x ). jkpq jk pq Taking a first order (r"1) neighbourhood set, the image recovery and boundary detection process in the presence of zero-mean additive Gaussian noise is put in correspondence with the minimization of a functional whose general structure is14: with H*(X, B)"H*(X, B)#jH (X) 1 2 C A (1) B D x !x 2 jk pq #t(b ) H* (X, B)" + + b (2) 1 jkpq jkpq * 2, N jk j/1, 2 p, q|N k/1, , M H (X)" + (y !x )2 (3) 2 jk jk 2, N j/1, 2 k/1, , M where * is a positive constant that will be identified later, j is the regularization parameter, and y denotes jk the noisy grey level of the pixel jk before image processing. The functional H* (X, B) is the regularization term, 1 and the functional H (X) is the data term. The minimization process is made with respect to the terms x and 2 jk b . jkpq For the sake of compactness, we define the dimensionless variable x !x pq u " jk jkpq * (4) that gives a measure of the gradient of the field X between two points jk and pq. The term b depends on jkpq u and assumes values over an interval [0, M ], with M '0. In particular, b "0 (resp., b "M ) jkpq b b jkpq jkpq b when the field gradient term u between the two points is sufficiently large (resp., is null), i.e. jkpq when a discontinuity exists. Thus, since X is the field representative of the grey levels of the image, b is jkpq null when there is an edge between jk and pq. A more detailed discussion, not needed here, is reported in Reference 14. 2.1.1. From H* (X, B) to H (x). Under proper conditions, that will be presented later, the minimization of 1 1 the functional H*(X, B) is completely equivalent to that of the simpler functional with H(X)"H (X)#jH (X) 1 2 A (5) B x !x pq H (X)" + + ' jk 1 * 2, N j/1, 2 jk p,q|N k/1, , M Under the hypothesis that each term b is independent from the others, the equivalence jkpq min H*(X, B)"min H(X) (6) (7) holds provided that the function ' fulfils the condition14 where '(u )" inf h(u , b ) jkpq jkpq jkpq 0)b )M (8) h(u, b)"bu2#t(b) (9) jkpq b The use of H(X) instead of H*(X, B) allowed by equivalence (7) delineates a simplified way to solve the original problem. ( 1998 John Wiley & Sons, Ltd. Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) 480 M. PARODI, M. STORACE AND C. REGAZZONI 2.1.2. From H (X) to H*(X, B). In a dual perspective, a given h(x) can be put into relation with a functional 1 1 H*(X, B) such that equivalence (7) holds, provided that the conditions on '(u) reported in the following are fulfilled. For u3[0,#R), Theorem 1 in Reference 14 provides the following three conditions on '(u) (which is supposed to be an even function) sufficient for the existence of the function t(b) (see expressions (2) and (9)) '(0)"!1 '(Ju) concave (10) lim '(u)"0 u?`= Under these requirements, t(b) exists and exhibits the following properties: t(0)"0 t(b) strictly decreasing for b3[0, M ]. b (11) t(M )"!1 b 2.2. Preliminary steps towards the circuit definition The minimization of H(X) can be related to the achievement of a stationary solution in an electric circuit consisting of two proper resistive nonlinear multiterminal elements and of linear capacitors. The most significant part of the resistive structure is a grid of identical nonlinear two-terminal elements. This grid is related to the structure of the regularization term. For the first-order neighbourhood set, each pixel jk interacts directly only with the four closest ones in the lattice S. Similarly, in a rectangular grid of two-terminal elements, each node is connected directly only to the four closest ones. The nonlinear characteristics of the resistors are related to the nonlinear nature of the function ' in H (X). 1 In the next sections, it is shown that the behaviour of the two resistive multiterminals can be described through potential functions that are homologous to H (X) and jH (X). In the circuit, the stationary values of 1 2 the capacitors voltages individuate the position of a minimum of the functional H(X). The outline of the remaining part of this paper is as follows: (1) The structure of the multiterminal resistor associated to H (X) is defined by taking as reference the grid 1 of (identical) nonlinear resistors. Then, the second multiterminal resistor is obtained in such a way that its potential function has the same structure of jH (X). The convergence of the complete circuit to 2 stationarity points is briefly discussed. (2) The relation between H*(X, B) and the nonlinear characteristic of the grid resistors is discussed. First, 1 a specific shape is taken for the nonlinear characteristic of the grid resistors and a corresponding regularization term H*(X, B) is obtained following the method described in Reference 14. Then, 1 a method is shown for the synthesis of the characteristic when h(u, b) is known. (3) The regularization parameter j for the case of null-mean additive Gaussian white noise is considered, showing its relation with the parameters of the circuit. Two examples of image restoration are given. (4) The limit case of noise absence is considered, showing that the corresponding circuit contains only the multiterminal resistor pertinent to H (X). A general property following from the particular topological 1 structure and two examples of edge extraction are presented. (5) Some general properties concerning the dynamic behaviour of the complete circuit are presented. Finally, the circuit is briefly compared with the classical CNN structure given by Chua and Yang in Reference 2. ( 1998 John Wiley & Sons, Ltd. Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) 481 MARKOV RANDOM FIELDS Figure 1. The M]N-terminal resistive network R 1 3. DEFINITION OF THE CIRCUIT STRUCTURE 3.1. Definition of the resistive multiterminal R 1 According to the first-order neighbourhood set defined in the previous section, we consider the structure shown in Figure 1, made up of (nonlinear) resistors connected so as to form a rectangular grid. The M]N nodes of the grid are connected to terminals. The resulting M]N-terminal resistive network is denoted by R . All the resistors are identical and voltage-controlled. The behaviour of R can be described in terms of 1 1 the M]N node voltages x with respect to a common node P, as shown in Figure 1. These voltages jk individuate a complete set15 of variables. For this circuit the global co-content function G is given by the 1 sum of the co-contents of the nonlinear resistors.15 Denoting by i"f(v) the voltage-controlled characteristic of each resistor, the corresponding co-content term 'ª (v) is v P0f(o) do 'ª (v)" (12) By means of the grid resistors, each inner node jk interacts directly with the four closest ones, that individuate the neighbourhood set N "M( j, k!1), ( j, k#1), ( j!1, k), ( j#1, k)N. jk Denoting the terminal nodes of one of these resistors by jk and pq, the pertinent voltage v is given by x !x and the global co-content G can be thought of as a function of the vector x of the node jk pq 1 potentials G (x)" + + 'ª (x !x ) (13) 1 jk pq 2, N j/1, 2 p,q|Njk k/1, , M The structure of G (x) is identical to that of the term H (X) in expression (6). G (x) corresponds to a given 1 1 1 H (X) when a proper characteristic f(v) for the nonlinear resistors is individuated. Though the voltages 1 x are the variables of both the co-content G (x) and the functional H (X), we maintain distinct symbols for jk 1 1 the field X of random variables Mx N and for the vector x collecting the circuit variables Mx N, for the sake of jk jk definiteness. ( 1998 John Wiley & Sons, Ltd. Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) 482 M. PARODI, M. STORACE AND C. REGAZZONI Figure 2. (a) The two terminal resistor R . (b) The (M]N#1)-terminal resistive network R jk 2 3.2. Definition of the resistive multiterminal R 2 The term jH (X) in expressions (1) and (5) can now be considered. As previously stated, this term forces 2 the x variables to remain close to the observed data y representing the image to be filtered. Being jk jk the y terms fixed, then, also jH (X) can be interpreted as the co-content of a proper resistive network jk 2 R . 2 Elementary considerations lead to represent the term j(y !x )2 through the two-terminal voltagejk jk controlled resistor R represented in Figure 2(a). The co-content of R can be written as jk jk 1 )ª (x )" [(y !x )2!y2 ] jk jk jk jk 2R jk (14) The network R is made up of M]N resistors R , as shown in Figure 2(b). Each of them is connected 2 jk between the reference node P and the pertinent terminal jk of R . 1 The co-content G (x) of R is given by the sum of the M]N terms )K 2 2 jk 1 G (x)" 2 2R 1 + (y !x )2! + y2 jk jk jk 2R , 2 N 2, N j/1, 2 j/1, 2 , , k/1, M k/1, M (15) The last term of this expression is constant, so that it does not take part in any minimization process. For this reason, the structure of G (x) can be considered equivalent to that of H (X). 2 2 3.3. The complete circuit The (M]N#1)-terminal R resulting from the connection of R and R is shown in Figure 3(a). Its 1 2 co-content G is given by G (x)#G (x) and exhibits the same structure of H(X) (see expression (5)). 1 2 The minimization of G (x) can be accomplished by connecting R to M]N capacitors as shown in Figure 3(b). For simplicity, all the capacitors are equal. The KCL equations for the M]N nodes yield ![C] dx "+G(x) dt (16) where [C]"diag MCN and C +G(x)" ( 1998 John Wiley & Sons, Ltd. D LG(x) LG(x) , , . Lx 2 Lx 11 MN Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) 483 MARKOV RANDOM FIELDS Figure 3. (a) The (M]N#1)-terminal resistive network R. (b) The complete circuit Let x(0) denote the initial state vector for the capacitors. Starting from this state, the vector x will change making the function G(x(t)) decrease with time until a stationary value is reached, as it follows from the expression of dG/dt A B A B dG dx T dx dx [C] "(+G(x)) ) "! )0 dt dt dt dt (17) A point x giving dG/dt"0 is a stationarity point for G and corresponds to an equilibrium point of equation (16). As usual in the circuit-based optimization algorithms (see, e.g. Reference 16), the stationarity point is assumed to be a minimum point for G(x), thus individuating a minimum also for H(X). 4. THE ELEMENTS OF R1 AND THE REGULARIZATION TERM 4.1. From the resistor characteristic to H* 1 Starting from the co-content G (x) of the grid of resistors with a given characteristic i"f(v), a method to 1 obtain a corresponding functional H*(X, B) such that equation (7) is fulfilled consists (i) in defining the 1 relation between H (X) and G (x). Then, (ii) a functional H*(X, B) can be obtained from H (X) following the 1 1 1 1 standard procedure given in Reference 14. We choose to exemplify the elements of the procedure making reference to the piecewise linear (PWL) resistor described below. 4.1.1. The reference resistor. The characteristic i"f(v) of the reference resistor is shown in Figure 4(a), together with its parameters I , » and » . The corresponding co-content function 'ª (v) is represented in 0 0 T Figure 4(b). Its analytical expression is reported in Appendix I. The stationary values of 'ª (v) correspond to the null-current points of f(v). Then, 'ª (v) is characterized by two regions of stationarity (individuated by DvD'» ) and by a stationarity point at v"0. T When the resistor connects the nodes jk and pq, and the difference between the node potentials x and jk x is sufficiently large, it works within one of the stationarity regions of its characteristic. This occurs when pq the two corresponding points of the image are separated by an edge. As previously stated, this case ( 1998 John Wiley & Sons, Ltd. Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) 484 M. PARODI, M. STORACE AND C. REGAZZONI Figure 4. (a) The PWL characteristic f(v) together with its parameters I , » and » . (b) The corresponding co-content function 'ª (v) 0 0 T corresponds to b "0. Since all the resistors have the same characteristic, from now on we will use the jkpq symbol b instead of b . jkpq The stationary point v"0 deals with the case of identical node potentials (identical grey levels in the image) and corresponds to the upper limit M of b. b The global co-content G (x) of the multi-terminal resistor R can be calculated by expression (13) taking 1 1 for 'ª (v) the expression (40). The PWL characteristic of Figure 4(a) could be viewed as a limit case of the fading resistor cited in Reference 17. 4.1.2. Application of the procedure. Owing to the electrical nature of the parameters and variables in G (x), 1 the correspondence with the dimensionless terms of H (X) must be defined through a proper normalization 1 process. In Appendix II it is shown that the relation between the co-content term 'ª and the term '(u) in equation (6) can be formulated as follows: 2 '(u)" 'ª (u» )!1 (18) T I » 0 T The variable u"v/» is a dimensionless term defined by taking *"» (see expression (4)). The function '(u) T T defined in equation (18) fulfils the three requirements (10) needed to apply the theorem given in Reference 14. Then, as shown in Appendix II, it can be related to C A BD A B AA B B » » T u2!1, b" T » » 0 0 b » h(u, b)" bu2! , b3 0, T 1#(1!» /» ) b » 0 T 0 0, b"0 u3 0, u3 » 0 » T » 0 ,1 » T (19) u3[1,#R) being » (20) M" T b » 0 The function h(u, b) yields a regularization term H*(X, B) with the same minima of the co-content of R . 1 1 ( 1998 John Wiley & Sons, Ltd. Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) 485 MARKOV RANDOM FIELDS 4.2. From H* to the resistor characteristic 1 The derivation of the characteristic i"f(v) from a given H*(X, B) can be obtained as follows. 1 We first observe that expression (2) for H*(X, B) yields the function h(u, b). The next step consists in 1 obtaining the even function '(u) according to condition (8). From a geometrical point of view, '(u) can be interpreted as the infimum of a b-parameterized family of quadratic functions. Then, following Reference 14, the function f (u)"'(Ju)" inf [bu#t(b)] u3[0,#R) (21) 0)b)M b is thought as the lower envelope of a one-parameter affine family. For any given b3[0, M ], and for u3[0,#R), the conditions for the function f (u) to be tangent to the b b-parameterized family of straight lines bu#t(b) are f (u)"bu#t(b) (22) df "b du (23) By combining these two conditions, a Clairaut’s differential equation18 for f (u) is obtained f (u)"u A B df df #t du du (24) Instead of considering the general solution (family of straight lines) of equation (24), we refer to the particular solution (envelope of the family of straight lines) that results by eliminating the parameter b from the system of equations f (u)"ub#t(b) (25) dt(b) u# "0 db Thus, according to definition (21), the even function '(u) is obtained from f (u2). As stated in the previous section, the function '(u) can be interpreted as the normalized co-content of each nonlinear resistor inside R . From expression (18), then, we find 1 d'ª d I » v 0 T' " #1 (26) i" dv dv 2 » T As an example, this procedure can be applied to the function h(u, b) reported in the previous subsection (see expression (19)) to find the PWL characteristic of Figure 4(a) (see Appendix III for the details). A B C A B D 4.3. A limit case for the f(v) characteristic: » "» 0 T For » "» , the f(v) characteristic and the related co-content function 'ª (v) are shown in Figure 5. 0 T Taking into account expression (19), the functional H*(X, B) can be written as + + h(u , b )D» "» #j + (y !x )2 (27) jkpq jkpq 0 T jk jk 2, N 2, N jk j/1, 2 p, q|N j/1, k/1, , M k/1, 2 , M Because of the ‘limit’ character of this case, that makes the f(v) characteristic discontinuous at v"» , the T infimum in expression (8) is degenerate, and it is achieved at b"1 whenever D x !x D(» and at b"0 jk pq T elsewhere. Its expression is u2 !1 forDu D)1 jkpq '(u )" jkpq (28) jkpq 0 elsewhere H*(X, B)D» "» " 0 T G ( 1998 John Wiley & Sons, Ltd. Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) 486 M. PARODI, M. STORACE AND C. REGAZZONI Figure 5. (a) The PWL characteristic f(v) of Figure 4(a), in the limit case » "» . (b) The corresponding co-content function 'ª (v) T 0 This is the so-called ‘truncated quadratic function’ considered in Reference 19, often used for image restoration11,14. The positions of the minima are the same as for the general case » '» . T 0 5. THE ELEMENTS OF R2, THE DATA TERM AND j The term jH in expressions (1) and (5) can be built up by using the co-content terms )ª (see expres2 jk sion (14)). The comparison between jH (X) and the corresponding sum of the dimensionless functions 2 2 ) (x)" )ª (x) (29) jk I » jk 0 T leads to identify the relative weight j as 1 j" R» I T 0 (30) In the case of null-mean Gaussian white noise with variance p2, the maximum a posteriori criterion14 leads to assign the regularization parameter j the value 1/2p2. Then, the second term in equations (1) and (5) can be written as 1 2p2 + (y !x )2 jk jk 2, N j/1, 2 k/1, , M (31) Comparing this sum with that originating by the ) and ignoring the constant term (see expressions (14), (15) jk and (29)), we obtain 2p2 R" »I T 0 (32) that directly relates the circuit parameter R to the variance of the noise model through the reference power » I . T 0 It could be easily shown that, by replacing R with a proper PWL resistor, it is possible to model jk asymmetrical noise distributions. This case is not addressed in this paper for brevity. ( 1998 John Wiley & Sons, Ltd. Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) 487 MARKOV RANDOM FIELDS Figure 6. (a) A noiseless piecewise constant image. (b) The image of Figure 6(a) corrupted by Gaussian noise G(k, p), with k"0 and p"35. (c) The steady-state image Table I t/q (q"500 ns) Signal/noise ratio (dB) 0 2 4 6 8 15)062 17)594 18)268 18)494 18)582 To simulate the behaviour of the complete nonlinear circuit, we define the lower and upper limits for the x ’s voltages as x and x , respectively. Then we take the scaling relation jk .*/ .!9 x !x .*/#x x "f .!9 jk jk .*/ 255 (33) between the grey levels f 3[0, 255] and the corresponding voltages x . jk jk The slope » /I of the central region in Figure 4(a) multiplied by C can be arbitrarily taken as a rough 0 0 reference term q for the definition of the time scale at which the circuit works. As a first example, a circuit defined by N"M"128, » "2» "1)8 V, I "1)8 mA, C"1 nF, x "0 T 0 0 .*/ and x "10 V is simulated. .!9 Figure 6(a) shows a locally constant image. Figure 6(b) shows the image to be processed, obtained by corrupting the image of Figure 6(a) by Gaussian noise G(k, p), with k"0 and p"35. Starting from the initial capacitor voltages x (0) generated through expression (33), the circuit simulation gives, after a few multiples jk of q, the steady-state values corresponding to the image in Figure 6(c). The signal/noise ratios measured at different times grow monotonically and are summarized in Table I. The example reported in Figure 7 concerns a natural image with more structural details. The original image shown in Figure 7(a) was corrupted by Gaussian noise G(k, p), with k"0 and p"10, thus obtaining the input picture represented in Figure 7(b). This image was processed by the circuit defined by N"M"256, » "2» "0)8 V, I "0)8 mA, C"1 nF, x "0 V and x "10 V. The circuit T 0 0 .*/ .!9 simulation gives, after a few multiples of q, the steady-state values corresponding to the image in Figure 7(c). The result of the simulation with »@ "2»@ "1)1 V and I@ "1)1 mA is shown in 0 0 T Figure 7(d). Some significant details still present in Figure 7(c) are lost (see, for instance, the background buildings). ( 1998 John Wiley & Sons, Ltd. Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) 488 M. PARODI, M. STORACE AND C. REGAZZONI Figure 7. (a) A noiseless natural image. (b) The image of Figure 7(a) corrupted by Gaussian noise G(k, p), with k"0 and p"10. (c) The steady-state image obtained with » "2» "0)8 V and I "0)8 mA. (d) The steady-state image obtained with »@ "2»@ "1)1 V and T 0 0 T 0 and I@ "1)1 mA 0 6. THE PARTICULAR CASE R,R 1 When pP#R, the j coefficient in expressions (1) and (5) annihilates. Owing to expression (32), each resistor R and the multiterminal R can be viewed as open circuits. 2 In this case, the node P individuates a cutset formed by all the capacitors C only, as shown in Figure 8. As a consequence, the overall charge of the cutset remains constant with time. ( 1998 John Wiley & Sons, Ltd. Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) 489 MARKOV RANDOM FIELDS Figure 8. The circuit in the case R,R and the C-cutset (dashed line) 1 The x vector of the state variables x evolves from the initial value x(0) (that represents, in this case, the jk image to be processed) so as to drive the co-content function towards a stationary value (see expression (17) with G"G ) under the constraint 1 + x (t)"C + x (0) jk jk 2, N 2, N j/1, 2 j/1, k/1, , M k/1, 2 , M C for any t'0 (34) Since the voltages x represent the grey levels of the image, the above constraint means that the average jk grey level of the image remains constant with time. It can now be considered a particular case when the grey levels are such that a portion P of the image never interacts with the rest of it during the dynamical evolution of the circuit (i.e. the nonlinear resistors connecting the boundary pixels of P to the outer region behave as open circuits). In this case the capacitors at the nodes of P form a cutset, so a relation completely analogous to (34) can be written for this set of capacitors and the average grey level inside P remains constant with time. This average level is the uniform grey level for all pixels of P (homogeneous area) in a further particular case when all these pixels interact with their neighbours inside P during the whole process. None of the above properties remains valid when R also is connected, because cutsets of capacitors are no 2 longer definable. In Reference 20, Perona and Malik showed that their anisotropic diffusion-based method for edge enhancement may be seen as a gradient descent on some energy function corresponding to H . In 1 Section 4, we showed a method to obtain a multiterminal resistor R whose co-content function G 1 1 is completely equivalent to H . When R,R , then, the circuit can also be viewed as the analogic 1 1 implementation of the anisotropic diffusion algorithm,20 and the discussion developed in Section 4 evidences the relationships between the structure of H and the characteristic i"f(v) of the nonlinear 1 resistors of R . 1 The input figure represented in Figure 9(a) was processed by a circuit with N"M"256, » "2» " T 0 0)8 », I "0)8 mA, C"1 nF, x "0 and x "10 V. At the end of the process, null currents flow through 0 .*/ .!9 each element of R , because each nonlinear resistor 1 (a) either works in the origin of the (v, i) plane of Figure 4(a) (b) or it behaves as an open circuit at voltage E , with DE D'» (extreme regions in Figure 4(a)). jkpq jkpq T ( 1998 John Wiley & Sons, Ltd. Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) 490 M. PARODI, M. STORACE AND C. REGAZZONI Figure 9. (a) The input image. (b) The edge extraction result obtained with » "2» "0)8 V and I "0)8 mA; (c) »@ "2»@ "1 V and T 0 0 T 0 I "1 mA; (d) »A "2»A "0)5 V and I "0)5 mA 0 T 0 0 Distinguishing the elements according to these two possibilities, the resulting distribution map yielding the edges of the image is obtained by associating a binary variable d "M0, 255N with each resistor, jkpq i.e. d "0 in the case (b) jkpq d "255 in the case (a). jkpq ( 1998 John Wiley & Sons, Ltd. Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) MARKOV RANDOM FIELDS 491 Figure 10. (a) The input image. (b) The edge extraction result obtained with » "2» "0)7 V and I "0)7 mA; (c) »@ "2»@ "0)6 V T 0 0 T 0 and I@ "0)6 mA; (d) and »A "2»A "0)8 V and IA"0)8 mA 0 T 0 0 The d distribution map gives the image shown in Figure 9(b). In this process a key role is played by the jkpq parameter » . The higher the value of » , the higher the gradient needed to individuate an edge between two T T pixels. Figure 9(c) shows the edge extraction obtained with »@ "2»@ "1 V and I "1 mA, and Figure 9(d) 0 0 T shows that obtained with »A "2»A "0)5 V and I "0)5 mA. As expected, the reduction in » results in an 0 0 T T increased sensitivity of the circuit to the presence of edges. ( 1998 John Wiley & Sons, Ltd. Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) 492 M. PARODI, M. STORACE AND C. REGAZZONI Another example with more structural image details is shown in Figure 10. The original image (Figure 10(a)) is first processed with N"M"256, » "2» "0)7 V, I "0)7 mA, C"1 nF, x "0 and T 0 0 .*/ x "10 V (Figure 10(b)). A lower threshold value (»@ "2»@ "0)6 V and I@ "0)6 mA) allows to extract 0 0 T .!9 more edges (Figure 10(c)), while a higher threshold value (»A "2»A "0)8 V and IA "0)8 mA) causes a loss 0 0 T of details (Figure 10(d)). 7. ON THE DYNAMIC BEHAVIOUR OF THE CIRCUIT AND ON ITS STRUCTURE 7.1. Time evolution of the vector x(t) Expression (17) clearly evidences that the vector x changes so as to make the co-content function G decrease with time. This ensures that the circuit reaches a stationary state, but does not give indications on the behaviour of x(t). The stationary state depends on both the passive resistors in R and the R elements of 1 jk R . 2 We first observe that lim G(x)"#R (35) ExEPR Following the notations of the Theorem B1 in Reference 21, we choose a Lyapunov function V(x),G(x). According to expression (16), we can write dx "!f(x) dt (36) LV(x) f(x)"(+G(x))T f(x)"f(x)T[C]f(x)*0 L(x) (37) with f(x)"[C]~1+G(x). Therefore as each solution x(t) is bounded. The calculations reported in Appendix IV show that, for any t and any j, k Dx (t)D)S jk jk where S is a function of the parameters of R and of the initial values x (0). jk jk (38) 7.2. Interpretation of the circuit in terms of the CNN architecture As shown in Figure 1, each node jk directly interacts with the four closest neighbour nodes of the grid in R . With this caveat in mind, it is rather simple to reformulate the circuit structure as a lattice of cells. As 1 shown in Figure 11, each cell can be viewed as a modified version of the CNN cell proposed in Reference 2. A first difference lies in the coincidence of the output variable with the state variable. Moreover, the overall nonlinearity of the cell is confined within the four current sources controlled by the pertinent terms x !x . jk pq Figure 11. The structure of the elementary cell in the CNN interpretation of the circuit. ( 1998 John Wiley & Sons, Ltd. Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) 493 MARKOV RANDOM FIELDS 8. CONCLUSIONS A nonlinear circuit solving a coupled edge-extraction and image restoration problem has been presented. It has been shown that well-defined relationships hold between the MRF functional defining the image processing method and the co-content of the resistive part of the nonlinear circuit. The analysis of such relationships leads to synthetize proper multiterminal resistors as parts of a circuit able to solve a specific image processing problem. The analogic image processing performed by the circuit has the following two main advantages: (1) the circuit yields a noise-robust solution for the considered problem, thanks to the intrinsic properties of the MRF approach; (2) the circuit structure implies parallel processing. Nonlinear circuits related to other classes of image processing problems could be identified by properly generalizing the described method. This could be useful for problems with strong real-time requirements (e.g. obstacle detection for autonomous driving,1 and dense optical flow field computation for motion compensation in image coding).22 Following Reference 14, the approach proposed in this paper is based on the widely diffused assumption that each term b is independent from the others (see Subsection 2.1.1). Removing this hypothesis could jkpq allow the definition of nonlinear circuits to perform more general image processes, at the cost of an increased complexity in the structure of the circuit’s resistive part. ACKNOWLEDGEMENTS Research supported by M.U.R.S.T. and by University of Genova. APPENDIX I The driving point characteristic shown in Figure 4(a) can be analitically expressed as follows: i"f(o)" 0 o)!» T I 0 (o#» ) !» )o)!» T T 0 » !» 0 T I 0o !» )o)» 0 0 » 0 I 0 (o!» ) » )o)» T 0 T » !» 0 T 0 (39) o*» T Then, according to equation (12), the co-content of each resistor of R is 1 I v2 0 DvD)» 0 » 2 0 v I I » 0 (DvD!» )2 » )DvD)» 'ª (v)" f(o) do" 0 T! T 0 T 2(» !» ) 2 0 T 0 I » 0 T DvD*» T 2 P ( 1998 John Wiley & Sons, Ltd. (40) Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) 494 M. PARODI, M. STORACE AND C. REGAZZONI APPENDIX II The normalized function '(u) is subject to the requirements (10). With this in mind, ' can be obtained from the co-content function 'ª (v) defined in expression (40) as follows. First, the parameter » is chosen so to relate the dimensionless variable u to the voltage v: T v"u» T thus » "* owing to expression (4). T Then, observing that 'ª (0)"0 I » lim 'ª (v)" 0 T 2 V?`= and recalling that we must have '(0)"!1 lim '(u)"0 u?`= we define 'ª (u» ) T !1 '(u)" I » 0 T 2 (41) which in turn yields » » T u2!1 DuD) 0 » » T 0 » » T 0 '(u)" ! (DuD!1)2 )DuD)1 » !» » T 0 T 0 DuD*1 (42) It is easy to verify that, for u3[0,#R) and defining z"Ju, the function f (u)"'(z)"inf h(z, b)" b » T u!1 » 0 » T (Ju!1)2 ! » !» T 0 0 C A BD AA B B » 2 0 » T » 2 0 ,1 » T u("z2)3 0, u("z2)3 (43) u("z2)3[1,#R) is concave. Then all the requirements (10) necessary to apply Theorem 1 in Reference 14 are fulfilled. We point out that the constant term in expression (41) is influent for the position of the minima of the functional H (X). As a consequence, the minimization of H (X) (related to ') in the 1 1 normalized domain is completely equivalent to the minimization of G (x) (related to 'ª ), pertinent to the 1 circuit. By following the methods described by Geman and Reynolds,14 a function h(u, b) that fulfils the condition (8) can be easily determined. Such a function yields a possible functional H*(X, B). Following Reference 14, the unknown term t(b)"f (u)!ub ( 1998 John Wiley & Sons, Ltd. (44) Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) with MARKOV RANDOM FIELDS 495 df b(u)" du (45) must be determined in order to obtain h(u, b) (see expression (9)). Taking into account expressions (42), (43) and (45), the function b(u) turns out to be » T » 0 » 1 T b(u)" !1 » !» T 0 Ju 0 A A B 0)u) B A B » 2 0 » T » 2 0 )u)1 » T (46) u*1 Thus, the inverse function u(b) can be defined over the interval [0, » /» ] as T 0 1 u(b)" [1#(1!» /» ) b]2 0 T Owing to the limit nature of the considered f(v) characteristic, the values of u outside the interval (47) CA B D I" » 2 0 ,1 » T are compatible either with b"0 (u'1) or with » 2 » b" T u3 0, 0 » » T 0 Then, recalling that z"Ju, the expressions (43) and (47) lead to C A A B BD f (u(b))" Taking into account » » T u!1 b" T » » 0 0 1 2 » » T b3 0, T ! !1 » !» 1#(1!(» /» )) b » 0 T T 0 0 C 0 D A B b"0 C A BD AA B B » 2 0 » T » 2 0 ,1 » T u("z2)3 0, u("z2)3 (48) u("z2)3[1,#R) expression (44), we finally obtain b (49) t(b)"! 1#(1!(» /» )) b 0 T over the whole interval [0, » /» ], as can be easily verified. T 0 The function h(z, b) reported in expression (19) can be directly obtained from definition (9) and from expression (48) » » » T z2!1 b" T z3 0, 0 » » » 0 T 0 » b » 0 ,1 h(z, b)" bz2! b3 0, T z3 (50) » » 1#(1!(» /» )) b T 0 0 T b"0 z3[1,#R) 0 C A BD A B AA B B ( 1998 John Wiley & Sons, Ltd. Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) 496 M. PARODI, M. STORACE AND C. REGAZZONI APPENDIX III From h(z, b) defined as in expression (50), we have (see expression (9)) b t(b)"! 1#(1!(» /» )) b 0 T (51) over the whole interval [0, » /» ]. T 0 Then, defining z"Ju (see expression (21)), we obtain the first equation of the system (25) » » b" T u T!1 » » 0 0 b(u) » b3 0, T f (u)" ub(u)! » 1#(1!(» /» )) b(u) 0 0 T b"0 0 C A BD A B AA B B » 2 0 » T » 2 0 ,1 » T u3 0, u3 (52) u3[1,#R) The second equation of (25) yields 1 u(b)" [1#(1!(» /» )) b]2 0 T Then, by eliminating b from (52) and (53), we obtain » u T!1 »0 (1!Ju)2 f (u)" ! 1!(» /» ) 0 T 0 (53) C A BD AA B B » 2 0 » T » 2 0 ,1 » T u3 0, u3 (54) u3[1,#R) It can be directly verified that the even function '(u)"f (u2) is given by expression (42), which yields the PWL characteristic of Figure 4(a). APPENDIX IV In order to evaluate the dynamic range of x (t), we first refer to Figure 12 and write jk dx 1 y i (t) jk"! x # jk # jk jk dt RC RC C i " + f(x !x ) jk pq jk pq|Njk Thus, the time evolution of x can be written as jk t i (q) y x (t)"x (0) e~t@RC# e~(t~q)@RC jk # jk dq jk jk RC C 0 and its absolute value is subject to the following inequality: P KP Dx (t)D)Dx (0)De~t@RC# jk jk ( 1998 John Wiley & Sons, Ltd. t 0 C D C D K i (q) y e~(t~q)@RC jk # jk dq C RC (55) (56) (57) Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) 497 MARKOV RANDOM FIELDS Figure 12. The second term at the r.h.s. of this inequality can be handled further, observing that Di D44I jk 0 owing to both the f(v) characteristic (see Figure 4(a)) and expression (55), and that t P0e~(t~q)@RC dq)RC (58) (59) We thus obtain the chain of inequalities KP t C D K A i (q) y 4I Dy D e~(t~q)@RC jk # jk dq ) 0# jk C RC C RC BP t e~(t~q)@RCdq)4RI #Dy D 0 jk 0 0 If the initial voltage value x (0) ranges between the values x and x , we can write jk .*/ .!9 Dx (0)D)xN jk with xN "maxMDx D, Dx DN .*/ .!9 From expressions (57), (60), and (61), we finally obtain where (60) (61) Dx (t)D)S jk jk (62) S "xN #4RI #Dy D jk 0 jk (63) REFERENCES 1. C. Koch and B. Mathur, ‘Neuromorphic vision chips’, IEEE Spectrum, 33(5), 38—46 (1996). 2. L.O. Chua and L. Yang, ‘Cellular neural networks: theory’, IEEE ¹rans. Circuits Systems, CAS-35(10), 1257—1272 (1988). 3. Int. J. Circuit. ¹heory. Appl., Special Issue on Cellular Neural Networks, CTA-20(5), (1992). ( 1998 John Wiley & Sons, Ltd. Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998) 498 M. PARODI, M. STORACE AND C. REGAZZONI 4. IEEE ¹rans. Circuits Systems, Special Issue on Nonlinear Waves, Patterns and Spatio-Temporal Chaos in Dynamic Arrays, CAS-42 (10), (1995). 5. Int. J. Circuit ¹heory Appl., Special Issue on Cellular Neural Networks II: Part I, CTA-24(1), (1996). 6. Int. J. Circuit ¹heory Appl., Special Issue on Cellular Neural Networks II: Part II, CTA-24(3), (1996). 7. C. Rekeczky, T. Roska and A. Ushida, ‘CNN based self adjusting nonlinear filters’, Proc. CNNA+96, Seville, Spain, 1996, pp. 369—373. 8. T. Roska and L. O. Chua, ‘The CNN universal machine’, IEEE ¹rans. Circuits Systems, CAS-40(3), 163—173 (1993). 9. T. Poggio, V Torre and C. Koch, ‘Computational vision and regularization theory’, Nature, 317, 314—319 (1985). 10. T. Poggio and C. Koch, ‘Ill-posed problems in early-vision: from computational theory to analogue networks’, Proc. Roy. Soc. ¸ond. B, 226, 303—323 (1985). 11. D. Geman and S. Geman, ‘Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images’, IEEE ¹rans. PAMI, PAMI-6(6), 721—741 (1984). 12. T. Sziranyi and J. Zerubia, ‘Markov random fields image segmentation using cellular neural network’, IEEE ¹rans. Circuits Systems, CAS-44(1), 86—89 (1997). 13. D. Geiger and F. Girosi, ‘Parallel and deterministic algorithms for MRF’s: surface reconstruction’, IEEE ¹rans. PAMI, PAMI13(5), 401—412 (1991). 14. D. Geman and G. Reynolds, ‘Constrained restoration and the recovery of discontinuities’, IEEE ¹rans. PAMI, PAMI-14(3), 367—383 (1992). 15. L. O. Chua, ‘Stationary principles and potential functions for non-linear networks’, J. Franklin Inst., 296(2), 91—114 (1973). 16. M. P. Kennedy and L. O. Chua, ‘Neural networks for nonlinear programming’, IEEE ¹rans. Circuits Systems CAS-35(5), 554—562 (1988). 17. T. Roska and L. O. Chua, ‘Cellular neural networks with non-linear and delay-type template elements and non-uniform grids’, Int. J. Circuit ¹heory Appl., 20(5), 469—481 (1992). 18. N. Piskunov, Differential and Integral Calculus, MIR Publishers, Moskow, 1977. 19. A. Blake and A. Zissermann, »isual reconstruction, MIT press, Cambridge, MA, 1987. 20. P. Perona and J. Malik, ‘Scale-space and edge detection using anisotropic diffusion’, IEEE ¹rans. PAMI, PAMI-12(7), 629—639 (1990). 21. L. O. Chua and D. N. Green, ‘A qualitative analysis of the behavior of dynamic nonlinear networks: stability of autonomous networks’, IEEE ¹rans. Circuits Systems., CAS-23(6), 355—379 (1976). 22. J. C. Brailean and A. Katsaggelos, ‘Simultaneous recursive displacement estimation and restoration of noisy-blurred image sequences’, IEEE ¹rans. Image Process., 4 (9), 1236—1251 (1995). ( 1998 John Wiley & Sons, Ltd. Int. J. Circ. ¹heor. Appl., 26, 477—498 (1998)

1/--страниц