Pixel labeling in a second-order Markov mesh

Signal Processing 12 (1987) 59-82 North- Holland

PIXEL LABELING IN A S E C O N D - O R D E R

59

MARKOV MESH

Vinciane L A C R O I X Philips Research Laboratory, Avenue E. Van Becelaere 2, Box 8, B-1170 Brussels, Belgium Received 6 January 1986 Revised 15 May 1986

Abstract. In this paper, we address the pixel labeling problem under the assumption that contextual information in the picture is encoded in the transition probabilities of a second-order Markov mesh random field. We show that the problem can be solved by a linear algorithm generalizing Devijver's F - G - H algorithm. In order to label the current pixel, all the dependencies inside a triangular past neighborhood are taken into account. However, the computation of the probability of the labels lying at the border of the neighborhood is based upon an hypothesis. The complexity of the algorithm increases with the size of the considered neighborhood. Some experiments illustrating the theory proposed show how to recover an artificial binary image from its noisy version and how the size of the considered neighborhood influences the recovering. The very good results justify the use of this algorithm and moreover show that a best compromise between computer-time and precision leads to the F - G - H algorithm.

Zusammenfassung. In dieser Arbeit untersuchen wir das Problem der Bildpunktetikettierung unter der Voraussetzung, dab die Kontextinformation des Bildes in Form der Ubergangswahrscbeinlichkeiten eines zweidimensionalen Markov-Feldes vorliegt. Wir zeigen, dab das Problem dutch den folgenden Algorithmus gel6st werden kann. Betrachtet wird eine kausale dreieckf6rmige Umgebung jedes Bildpunktes. Wir definieren einen Satz von Funktionen, yon denen jede mit einer Diagonale der betrachteten Umgebung verkniipft ist; die Funktion liefert die Wahrscheinlichkeit des Auftretens eines gegebenen Etiketts entlang dieser Diagonale. In dieser Umgebung ist es m6glich, Funktionen, die zu benachbarten Diagonalen geh6ren, mit Hilfe einer 'lokalen Zerlegung' zueinander in Beziehung zu bringen. Eine weitere Beziehung wird ben6tigt fiir die Funktion, die zu der Diagonalen an der Grenze der betrachteten Umgebung geh6rt. Eine geeignete Hypothese erlaubt, diese Funktion aus vorher gewonnenem Wissen zu berechnen; d.h., aus Funktionen, die bereits in der Umgebung benachbarter bekannter Bildpunkte berechnet wurden. Diese Beziehung bezeichnen wit als die 'Kreuzglied-Rekursionsbeziehung'.AnschlieBend wird an jedem Bildpunkt die folgende Prozedur gestart. Zuniichst wird mit Hilfe der Kreuzglied-Rekursionsbeziehungdie Funktion berechnet, die der Diagonalen an der Grenze zugeordnet ist. Dann werden mit Hilfe des Verfahrens der "lokalen Zerlegung' Schritt fiir Schritt die Funktionen berechnet, die zu den anderen Diagonalen geh6ren. In h Schritten ist dann die Wahrscheinlichkeit jeder der m6glichen Etikettierungen am beobachteten Bildpunkt bekannt; h ist die Zahl der Diagonalen in der betrachteten Umgebung. Schliel31ich wird die Etikettierung mit der gr613ten Wahrscheinlichkeit als beste Sch/itzung der wirklichen Etikettierung gew/ihlt, die zu dem Bildpunkt geh6rt. Zu jeder Gr613e der betrachteten Umgebung existiert ein Algorithmus, dessen Komplexit/it mit der gew/ihlten Gr613e w/ichst. Dieser Algorithmus ist eine Verallgemeinerung des F-G-H-Algorithmus yon Devijver, bei dem die betrachtete Umgebung zwei Diagonalen enth~ilt; G und H sind hierbei die beiden Funktionen, die der ersten und zweiten Diagonalen zugeordnet sind. Basierend auf dieser Theorie zeigen einige Experimente, auf welche Weise ein Bild aus einer verrauschten Version rekonstruiert werden kann, und in welcher Weise die Gr68e der betrachteten Umgebung den Prozel3 der Rekonstruktion beeinfluBt. Die sehr guten Ergebnisse rechtfertigen den Einsatz dieser Algorithmen; dariiber hinaus zeigen sie, dab ein optimaler KompromiB zwischen Rechenzeit und Genaugkeit zum F-G-H-Algorithmus fiibrt. R6sum& Dans ce papier, nous adressons le probl~me de 1'6tiquetage de pixel sous l'hypoth6se que l'information contextuelle de l'image r6side dans des probabilit6s de transition d'un Markov mesh random field of second order. Nous montrons que le probl6me peut ~tre r6solu par un algorithme lin6aire qui g6ndralise l'algorithme F - G - H de Devijver. On consid~re toutes les d6pendances ~ l'int6rieur d'un voisinage triangulaire appartenant au pass6 du pixel afin de lui assigner une 6tiquette. Cependant le calcul de la probabilit6 des 6tiquettes aux fronti6res de ce voisinage repose sur une hypoth6se. La complexit6 de l'algorithme augmente avec la taille du voisinage consid6r6. 0165-1684/87/$3.50 O 1987, Elsevier Science Publishers B.V. (North-Holland)

V. Lacroix / Pixel labeling in a second-order Markov mesh

60

Quelques exp6riences illustrant la th6orie propos6e montrent comment restaurer une image binaire artificielle ~ partir de sa version bruit6e et comment la taille du voisinage consid6r6 influence la restauration. Les tr~s bons r6sultats obtenus justifient l'usage de cet algorithme et de plus, montrent qu'un meilleur compromis entre le temps ordinateur et la pr6cision conduit ~ l'algorithme F - G - H . Keywords. Context, image modeling, image restoration, Markov mesh, pixel labeling.

I. Introduction

We shall consider a two-dimensional rectangular lattice. To each pixel (or node) (m, n), is associated an unknown label A,,,n and a known feature vector xm,, containing imperfect information about the corresponding label. Our problem will be to label the pixels given the information comprised in the feature vectors. In the literature there exist several methods in order to label the pixels. Among them the contextual methods--so-called because they are based on spatial interaction models--take the neighborhood of a pixel into account in order to label it. The order of the model reflects the size of the considered neighborhood. A model is classified as either causal or noncausal depending on the nature of the neighborhood. In causal models, the concept of the past of a pixel is introduced and only the past neighborhood influences the label of the current pixel. Among the spatial interaction models, the Markov Random Field (MRF) models are often used (see [1, 3-8, 11-15, 17, 18]). The equivalence between MRFs and Gibbs distributions leads to interesting results in relaxation labeling [11, 17]. MRFs include both causal and non causal definitions. However, causal MRFs are generally called Markov mesh random fields or simply Markov meshes. The MRF literature offers basically two different approaches depending on the MRFs definition. According to Rosenfeld [16], MRFs may be defined as either wide sense or strict sense. The wide sense Markov fields are defined directly in terms of least squares estimates (see [5, 6, 18]) and the strict sense Markov fields--which we shall consider in this paper--are defined in terms of conditional distribution functions (see [1, 3, 4, 7, 8, 11-15]). More precisely, the strict sense Markov field models assume that the probability of a label )t, given all other labels, is equal to the probability of the label A, given the label6 of its neighborhood only (see [3, 4, 7, 11, 14, 15]). The assumption is exactly the same in the strict sense Markov mesh models, except that the given labels and the neighborhood are restricted to the past (see [1, 8, 12, 13, 14]). The latter models are thus seen as the two-dimensional generalization of the Markov chains. The spatial interaction models have already found various applications in image processing. Noncausal autoregressive models assimilate labels to grey-levels, and seem particularly suited to image restoration [5]. Moreover, they are successfully used in texture modeling [6]. Noncausal MRF models on the other hand are applied with a restricted number of labels viewed as identification of regions. Therefore, they are more appropriate to image segmentation [ 11 ] or to restoration of binary images [ 17]. The latter models are also quite effective in texture analysis and synthesis [7]. The success of these models in real image applications is in the hands of the learning phase, that is, the process which fits the parameters of the model to the real image [17]. Markov meshes could as well be applied to image segmentation, texture analysis and synthesis but the learning phase is still under investigation and will be discussed in another paper [9]. In this paper we shall be concerned with the developing phase, that is, the elaboration of a technique to solve the pixel labeling problem in a second-order Markov mesh. In order to overcome the exponential complexity of the problem, we shall need a simplifying hypothesis. In the proposed solution, we consider Signal Processing

V. Lacroix / Pixel labeling in a second-order Markov mesh

61

a past triangular neighborhood of fixed size at each pixel and we divide it into an ordered set of past diagonals. Then we define a set of functions each of which is associated to a diagonal and corresponding to the probability of occurrence of given labels along that diagonal. Within the considered neighborhood, a function associated to a diagonal can be computed from the function associated to the previous diagonal. We call that property the 'local decomposition relationship'. The adjective local underlines the fact that only the current pixel and its associated past neighborhood are involved in the relationship. But we need a hypothesis at the boundary diagonal of the neighborhood. The hypothesis assumes a form of independence between image elements. These elements are located at each extremity of the boundary diagonal so that the larger the neighborhood is, the weaker the hypothesis becomes. Under this hypothesis, it is possible to deduce the function associated to the boundary diagonal from previous knowledge, that is, from functions computed in the respective neighborhoods of the past adjacent pixels. This relationship is called the 'lattice recurrence relationship'. Now the term lattice makes explicit the interdependence of functions defined inside neighborhoods associated to different pixels, more precisely to adjacent pixels. The complete algorithm is generated by using the lattice recurrence relationship, the local decomposition relationships and the boundary conditions as summarized below. At each pixel, the function associated to the boundary diagonal is first computed by the lattice recurrence. Then, step by step, using the local decomposition, the functions associated to next diagonals are computed until the function associated to the current pixel. The last function gives the probability of occurrence of a given label at the current pixel. The final step consists in choosing the label with the maximum probability. The complexity of the algorithm increases with the size of the considered neighborhood. The F - G - H algorithm [8], from which the proposed generalization comes, considered a past neighborhood containing two diagonals. 1. I. T h e M a r k o v f i e l d m o d e l

Following [8], let: • VM.N = {(m, n ) : 0 ~< m ~< M rows, 0 ~< n <~ N columns} be a finite integer lattice, • (a, b) be a pixel at the intersection of row a and column b, • Aa,b and Xa.b respectively be the label and the feature vector associated with the pixel (a, b), and • va.b"(~)and --a, bv~x) denote the rectangular arrays of labels and feature vectors, respectively (see ~(~)-~,bin Fig. 1). As we shall consider a Markov mesh, the concept of the past of a pixel is to be defined. We adopt the definition first introduced by Abend et al. [1]: the past of the pixel (a, b) is the set of pixels {(m, n): m < a or n < b}. In this paper we consider --M,N ~/(~) as a homogeneous second-order Markov random mesh [1, 8].

Ao,o

Ao,l

• • •

Ao,b i

Ao,b

Al,o

Al,l

•

Al,b

Al,b

Aa-l,O

Aa 1,1

" " "

Aa-l,b I

Aa-l.b

A,,,O

A.,,

• • •

A.,~,_,

A~,,b

•

•

1

Fig. 1. Definition of v(~) --a,bVoI. 12, N o . l , January 1987

62

V. Lacroix/ Pixel labeling in a second-orderMarkov mesh

So we have, by definition:

Second-order Markov Mesh: P[Aa, b [{A,.,. : m < a or n < b}]

=

P[Aa, b lAa-l,b ; Aa, b-1].

Homogeneous: Let q, r, s be some pixel labels; then, for 0 < a <~ M, 0 < b ~< N, Pql~,s- P[A~,b = q lA~-Lb = r, Aa,b-1

= S]

is independent of a, b.

The boundary conditions may be introduced in two different ways. The first one is to consider the image as a periodic random field on a toroidal lattice. The second one is to consider the upper left pixel, the upper row and the left most column separately [8]. From an image processing point of view, the latter seems more natural. Therefore, we have the following boundary conditions:

Pq "-- P[A~ b = q]

if a = b = 0 (upper left pixel),

Pqfr,~= ~PqlG,s -- P[Aa, b = q IA~,b-1 = s] / [Pqr,.D -- P[A~,b = q[Aa-~,b = r]

if N t> b > a = 0 (upper row), if M i> a > b = 0 (left most column).

As a result of the second-order Markov mesh definitions we have the following property: For any

(a, b) ~ VM.N, b

P[ V(x)'-~,bj= I~I H P[Xm.lXm-l..; x . . . . . ,1. m:O

n:O

1.2. The noise model As in [8], we assume that the noise-free image is corrupted by conditionally independent noise. It follows that pq(X~,b)"---p(Xa,b lAa,b----q) depends on Z only. This implies that b

(x)

P[Va.b[ V~b)] = 1~I I] p(x,,,,,,lX,.,,,). m=O

n=O

1.3. Statement of the problem In this p a p e r we assume that both the Markov and the noise models are known to us. Our problem is to recover the original array V~!N from the measurements V~!N. In order to decide the label on each pixel, we shall find the ' m a x i m u m a posteriori estimate' ~ for the labels according to the following rule: h,,b

=

q

for P[Aa, b

--

q l V(a,~] max P[Aa,b =

r

=

rl

v(Zg].

1.4. Organization of the paper In Section 2, we give the motivations for this work and its origin: we describe the F - G ' - H algorithm and an intuitive feeling of a possible generalization. In Section 3, we first set up a new frame needed for the development of the general algorithm, then we show what F, (3, and H become in this new frame, and finally we give the generalization. In Section 4, the complete algorithm is presented, including boundary conditions and the special case of one-diagonal neighborhood. In the last subsection, a summary of the derived formulae is given in the form of an algorithm written in a pseudo-computer language. Signal

Processing

V. Lacroix / Pixel labeling in a second-order Markov mesh

63

Section 5 contains experimental results demonstrating the effectiveness of the proposed technique and showing the influence of the size of the neighborhood considered at each pixel. Section 6 proposes a brief discussion. Finally, Section 7 gives the conclusions of the work.

2. Motivation

2.1. Previous work: The F - G - H

algorithm

In [8], a solution to the problem is proposed with the F - G - H algorithm. This algorithm is expressed in terms of joint likelihoods and needs a normalization in order to avoid underflow when simulated on a computer. The normalization is introduced by the use of a posteriori probabilities. Therefore, the F - G - H algorithm translated in terms of a posteriori probabilities will briefly be described below. For the purpose of this work, we shall name the a posteriori form of ~, ~, N, respectively ~-o, ~ , ;T2 as shown in Fig. 2. ~-° b(q)-- p[A. b = q[ ,/(~)I V a , bJ~ ,

,

~la,b( r, S) -- P [ A a _ l , b = r ; }~a,b-I = S I V ( 2 ) \ { X a , b } ] , ~ 2 . b ( t , U, Z)) ":- e [ l ~ a _ 2 , b = t'~ A a

0

.........

b

l.b-I = U

0

; h a , b - 2 = 1)I V~g)b\{Xa.b, Xa ,,b, Xa,b-l}].

.........

b

0

0 0

G

"'"

b

v 2\ (xab, Xa-l,b~ Xa,b-1}

¢Z

8

(a)

IZ

(b)

! (c)

Fig. 2. Arguments in the definitions of (a) ~.b, ~ (b) ~-~~,b, (c) 5~, 2 b. It is easily seen that these functions are related by 1

.

b(q) = - - Y ~Tl.b( r, s)Pqlr.spq(Xa.b), N O r,s

whereY°=~

(1)

~ ~ , bl ( r , , S ,)Pq,l<,,pq,(X~,b),

q' r',s'

ff;la,b(r , S) = l-~i~l E J;],b(t, U, v)P~lt,~P~l.,opr(Xa_,,b)p~(x~,b_,).

(2)

t, U,1)

w h e r e ? ( ' = }2

~

r',s' t',u',v'

2 , U,, V ,)Prl,...P.l.,,~,p~(xa_l.b)p.(x..b_~). ~;.,b(t, Vol. 12, No. 1. January 1987

64

V. Lacroix / Pixel labeling in a second-order M a r k o v mesh

0

. . . . . .

b

~'a,h

Va(1_l,h_l

t

v I

s

q

R'-{~) Fig. 3. The basic assumption. We shall call the above equalities 'local decomposition relationships'. The adjective local is used in order to underline the fact that all j~-o, ~ , o~2 involved in (1) and in (2) refer to one fixed pixel (a, b). If O~2b(t, U, V) were known for (a, b), one could compute ~-~,b(r, S) and 5~°~,b(q). Devijver introduced a hypothesis under which it is possible to relate :T2,,b to ~,-,,b, gea,b , and ~°~_~,b_b SO that if o~° and ~-1 are known at the borders, the whole process is completely determined. Let R2b and C2,b respectively denote a truncated row and truncated column of V,.b as shown in Fig. 3. The hypothesis assumes that (l . D2(x)'l f,2(x) ~ x) .... b-2,,-~,b J and (A~-2.b ;'-a,b J are stochastically independent conditionally upon (Aa 1,b , ; V~-~,b-~). Formally, if t,'~2(x).

2

P(A. 2.b;Xa,b-2;'~o.b ,Ro, blX.-,,b-,; I¢22,,~-,) =P(A~

2,b;g'-'2tx)ll~a_l,b *'-~a,b

l"~ If(x) I b-2~l'ta,b .IT~2(x)[l~a_l,b_l ; ~/4£~_)1,b_1) ' • a l ' , b - l lhl D [klta,

then the following relationship holds: ff2,h(t, u, v)oc ~ l ~,b(t, U)~,b-l(U, V)/~'O~-~,b-I(U).

(3)

We shall call (3) the 'lattice recurrence relation' for the 'distance' 2. The term lattice is now used because (3) relates quantities referring to adjacent pixels. The distance refers to the distance (in pixel) between the truncated row and truncated column of Va, ,,(x)b involved in the hypothesis. In summary, the computation of ~a.b(q) for all q's on each pixel (a, b) is completely determined by (1), (2), (3) and the boundary conditions for ~ and o~1.

2.2. Towards a generalization We could think of two different ways of transforming the F - G - H simpler algorithm. Consider (1) alone:

~-~.b(q) = ~

1

algorithm. The first one leads to a

Y. o~.b(r, s)Pqlr.,Pq(X~,b). r.s

If ,~,b(r, S) were known for (a, b), the local problem would be solved as well. We would therefore need a--probably stronger--hypothesis in order to relate ,~.b(t:, s) to ,°,~_l.b(r) and ~-°~.b_I(S) SO that the Signal Processing

65

V. Lacroix / Pixel labeling in a second-order Markov mesh

0

. . . . . .

)\

b

,,,,,

27a 2,b~ X a - l , b - l ~ 3;a,b-2}

/3

O.

Fig. 4. Arguments in the definition of ~3. hypothesis equation (1) and the boundary conditions on 5~' completely determine the process. In fact, here we propose to use only one local decomposition relationship and to initialize the process one pixel ahead, by means of a lattice recurrence relation for the distance 1. ~3 The second one leads to a more complex algorithm. Let ,~a,b(a, fl, % 6) be a function as shown in Fig. 4, defined by

ff3.,,(o~,,B,"/, 6)=P[Aa-3.b= o~,A,~ 2.~, ,=,8, A,,-Lb-2=Y,A,~.b 3=61 (x)

Va, b\{Xa.b, Xa l.b ; Xa, b-1 ;Xa-2,b ; X a - l . b - 1 ; Xa.b-2}]"

Readily we have

,~,b( t, u, v ) = ~ where 3 ;2 =

~

a,~,%8

~

~3,b(a, fl, y, 6)P,l~.t~Pult3.~Pvl~.apt(x~-2.b)p~(x~-Lb)p~(xa.h-2),

y~

t',u',v' ot',,B',~/' 6'

(4)

~3 . , fl',. y , 6)Pt'l,~'trpu'lo'y'P~'l~/,6 . . ~a,~,ta pc(x,~-2,b)p,,,(Xa-j,h)p~'(X,,,b-2).

Like equations (1) and (2), expression (4) is a local decomposition relationship that we shall generalize in the next section. In the same manner, if we introduced a lattice recurrence one pixel further, i.e., for the distance 3, ~3a,b(Ol, ~, ~/, 6) would be known for (a, b) and the problem would be solved. We should introduce oZ-.3 0722 a - - p r o b a b l y weaker--hypothesis in order to relate ~,,h to ~o~2 - - L b , ~,,b-~, and ~O721 , ~,b ~- In such a case, (1), (2), (4), the new hypothesis, and the boundary conditions on ~o, ~ and ,~2 would completely determine the process. Intuitively, we see that we could similarly define further functions which involve labels along further diagonals of the past of (a, b) and introduce weaker hypotheses so that we could approach the exact solution as close as we want. The functions ~v°, , ~ , ~:2, ~3, and their arguments as they are, do not provide a general form of the local decomposition relations; therefore, our first aim will be to define a new notation more appropriate to the proposed generalization. This will be done by changing from pixels to diagonal arguments. In this section, the hypotheses are not given, and our second aim will be to specify them and to show that they are weaker as the algorithm thus generated gets more complex. Vol. 12, No. 1, January 1987

66

V. Lacroix / Pixel labeling in a second-order M a r k o v mesh

truncated rectangle V:~ 'x}

b (col)

diagonal nk(~X) ~'a,b

a

I II I II

I

(row) Fig. 5. Definitions of diagonal and truncated rectangle.

3. Linear time algorithm

3.1. Generalized frame All ~o, ~1, if2, ~3 functions introduced in Section 2 involved labels of v(a) --a,b along a given diagonal, and feature vectors of various truncations of ~/(x) r a,b. We shall then introduce some geometric objects called diagonals and truncated rectangles associated with the pixel (a, b). We assume in this section that k < min{a, b}; the boundary conditions, i.e., a or b ~< k, will be treated in Section 4. Let

Dka,b "--{(a - k + j , b-j):O<~j<~ k} denote the kth diagonal of Va,b and let

Vkb-- {(m, n)e Va,b:m+n<~a+b-k}, denote the truncated rectangle (see Fig. 5). In the previous section we referred to a truncated row and column of V.,b in order to establish the lattice recurrence. We shall now label the various truncated row and column of Vab, k b, respectively (see Fig. 6). When necessary, an upper script (A) or (x) , Ra,k b, and Ca, will be added as before in order to distinguish between associated labels and feature-vectors. From now on we shall deal with diagonals instead of pixels. Therefore, we shall first adapt the definitions concerning pixels already presented in Section 1. Then, we shall introduce the decomposition relationships and the lattice recurrence relation. When writing labels and feature vectors, we adopt the convention of using lower case letters for pixels and capitals for diagonals. Let Qk denote some list of pixel labels along the diagonal Okb:

Qk-- {qo . . . . . q i , ' " , qk}. Signal Processing

V. Lacroix / Pixel labeling in a second-order Markov mesh

67

b (col) X0,b

,b Xa-k,b

(row)

o,0

a

o,b-k

I

Fig. 6. Definitions of truncated row and column. We define the M a r k o v transition probabilities from one diagonal to the next one as

pOkiRk+,--p[Dka!~ )

Qklmk+'(a) ltJ

a, b

~

R k+l]

.

With the introduction o f these new definitions, the sequence of diagonals Dkb appears to form a M a r k o v chain. O f course, the transition from one diagonal to another can be c o m p u t e d from the pixels transition probabilities as follows:

PokiR "' --l-I Pq,l,,,r,+,, J where the rj's are the elements o f R k+~. Finally, we define po ~(Dak(~,)) - p [ D ak!~,)i ~k(~) ,~,o,~ = Q~]. Likewise, po~(Dk!~ )) can be c o m p u t e d from the pixels probabilities as follows:

P& (Dk:~ ))= 1-I Pq,(X~-k+:,b-:). J

3.2. F-G-H in the generalizedframe The local decomposition relationship involved only pixel (a, b), so that we shall omit the indices as long as no confusion is possible. W h e n changing the pixels arguments to diagonals, ~o, ~1, ~2 b e c o m e respectively: ~-O(q) ~ 5~O(QO) = p[Do(~) = QO I vO(x)]

where QO= {qo}, qo = q,

~:'(r, s ) ~ ~ ' ( R ' ) = P[D'(A) = R ' I V '(x>]

where R l = {ro, rl}, ro = r, rl = s,

~.~2(I, U, U) "-)' ~ ; 2 ( $ 2 ) = P[D2(~) = $2[ V 2(x)]

where S 2 = {So, Sl, s2}, So = t, s~ = u, s2 = v. Vol. 12, No. 1, January 1987

v. Lacroix / Pixel labeling in a second-orderMarkov mesh

68

We can rewrite (1) and (2) as follows: 1

(1)-* :~O(QO) =~-6 ~ :~(RI)PQ°IR'Po°(D°(X)),

(5)

Rt

where N O= ~ ~ ~;I(R'I)PQ,OIR,'pQ,°(D°(X)), Q,O R,t

(2) -> ~I(R1) = ~ ~ ~;2(S2)pR'Is:PR'(D'(X)),

(6)

vv S 2

where Azl = E E R ,1 S, 2

~I;2(S'2)pR"Is'2PR"(DI(X))"

In the lattice recurrence relation (3) the indices are different, so we specify them. (3)--> :~2,b($2) = S :~-l'b(SE(2)):~'b-~(S2(O)),

(7)

~a_,,~_,(S~(O, 2)) where N =

P[ V~._,,b]P[V~,b-,] P[ V~,b]P[ V°~_I,b_I]

and (see Fig. 7) S 2 ( 2 ) = {SO, s i t ,

S 2 ( 0 ) = {Sl, S2} ,

S2(0, 2) = {sit

are truncations of S2.,b= {So, s~, s2}. In (7) we wrote explicitly the proportionality factor N which is a normalization factor; we shall see in Section 4 that it is not necessary to compute it. In this new frame the problem becomes: at each pixel find 5~° for all labels QO. The solution becomes: at each pixel consider a triangular neighborhood containing two diagonals; compute first, for all possible labels, the function ~2 associated to the boundary diagonal by means of the lattice recurrence (7), then, using the decomposition formulae, compute ~ and finally 5~°.

3.3. Generalization of the F- G - H algorithm We consider at each pixel (a, b) a triangular past neighborhood containing h diagonals. Within this )--" P[O~,b = Q'I V']. neighborhood we define a list of functions associated to the ith diagonal: ~,.b(Q ~

b-1

b

b-1

b

.,=.

s=(o)

a-1 82 = Y']

a

(a) s2(2)

a

g

(b) S2(O) Fig. 7. S's involved in the hypothesis.

Signal Processing

(¢) s=(o,2)

b

V. L a c r o i x / P i x e l l a b e l i n g in a s e c o n d - o r d e r M a r k o v m e s h

69

Each function indexed by i is related to the function indexed by i + 1 according to a local decomposition relationship: i 1 ~'(Q )=~-7 ~

~;'+l(gi+l)Po'lR'+'Po'(Di(X)),

(8)

R i+1

whereyfi=

~ ~ ~;i+l(R,i+l)pQ,,[R,,÷,pQ,,(Di(X)). Q,i R,i+t

If ~h were known, one could compute all ~ for i < h, so that, in h steps, one reaches the goal, that is, the computation of ~-o. We shall prove (8) for i = 0. The generalization to other cases is straightforward.

Proof of equation (8) P [ D ° = Q°I V 1] = Z

P[D 1= R 1, D O= Q°I

V 1] = Y

R1

P[D ~= R~l V']poota,,

Rt

p[Do= QOI Vo] _- PEDO= QO[ v l k..)DO(X)] _

P[ D ° = Q°; D°(~)I v l ]

p[ D°(X) I V 1]

P[D°= Q°I V']pQo(D °(~)) Y~o,o P [ D ° =

Q'° I V']po,o( D °(x))

~R' P[ D1 = R'I V']PQ°IR'PQ°(D°(X)) =~o,O~g,, n[ P' = R" I V']Po,olR,,po,o(P°(X))"

[]

As for the F - G - H algorithm, we shall introduce a hypothesis which enables us to compute ~h, the function associated to the diagonal located at the boundary of the considered neighborhood, from the h-1 h-1 h--2 knowledge of ~a-l.b, ~:a.b-1, and ~a-l,b-1, SO that if all ~ i for i < h are known at the borders, the process is completely determined. For the time being, we shall be considering the case h I> 2 only; h = 1 will be treated as a special case of neighborhood in Section 4.2. The hypothesis assumes a form of independence between elements of the image which are at a characteristic distance from one another. More precisely, the hypothesis assumes that ( ) t a , b _ h ," Dl~a, h ( xb) ~] and ( • a - h , b ; f~a,b ' , h ( x ) ~, are stochastically independent, conditionally upon /~'-'a-l,b-1, /~h--2(X) . 1/h-2(x) --a-l.b-~J. Formally, if p[Aa_hb;Aab_h, , , =PEAa_h

f,h(x), oh(x) , " J a . b , "l'~a,b

lr~h-2(A) . l z h - 2 ( x ) q LJa--l,b--I , Va-l,b-lJ

b . f ' h ( x ) ll-~h--2(A) • | z h - 2 ( x ) 1 D F i , , "-~a,b [ z " ' a - l , b - 1 , C a - l , b - l J l l r t a ,

. D h ( x ) l l ' ~ h - 2 ( A ) • Lr h - 2 ( x ) 1 b - h ~ l"~a,b [ L S a - - l , b - - I , V a - l , b - l J ,

then the following equality holds: nl'r~h(A)=

rL

,o.b

QhI l/-h(x)l --o,b J =

I --o,b J =

where N -

N

proh-l(x) llh-l(x)lDFl')h--l(a) __ |Th-l(x)l t ~-,,b = Qh(h)l --~-l,b J--t'-'a,b-1 -- Qh(0)l--~,b-1 a prDh-2(A) Ll'~'a

I , b - 1 --~

Qh(0, h)[

~7h-2(x) 1 Va_l,b_lJ

(9)

p[ r~zh-l(x)l pr I z h - - l ( x ) " I a-l,b I L Va, b-I J P[ Vh.b]P[ vhz2!~),]

and (see Fig. 8)

Qh(h)

= (qo, • • •,

are truncations of

qh-1),

Qh = (qo, ql,''',

Qh(O)= qh-1,

(q~, • • •,

qh),

Qh(0, h) = (ql . . . . .

qh-l)

qh). Vol. 12, No.l, January 1987

70

V. Lacroix / Pixel labeling in a second-order Markov mesh

b

b

vh-*(~) a-l,b

vh-2(~) a-l,b-1 q0

FIJ'

Dh-I(X) a- l,b =:2'>

a

~h-l()~)

t ' ~ l

I

Dh-2(~)

~

a- l,b- l ~ J

, ...

ih-,ll

qh

a

F IV

qh

/-)h-l(A) -- Qh(h )

D h ,,,~-1 -l('x)

(b)

_

-__-

Qh(O)

(¢)

n h-2('x) = Q h ( o , h )

Fig. 8. Q's in the hypothesis.

Proof of equation (9) ( u p p e r scripts (A) and (x) will be omited below). We consider 2 ~< h ~ min{a, b}. Note that (see Fig. 9) h

h

h-2

h

V~,b = Ra, b co Va_l,b_ l k..) C a , b. Using this d e c o m p o s i t i o n we can write p [ D h , b = Qob, h V h~b] h-2 =P[Aa-h.b = qo, Da-l.b-1 =

= P[Aa-h,b

=

qo,

Qh (0, h),Ao, b_h = qhh, Ra.b, h

C ah, b ll~a,b-h =

-Lb-, = Qh(h), qh, D ah-2

x P[Xa,b-h=q~, Da-l,b h-2 I = Qh(0, h), h P [ A a - h ,b = q o , C a , b, A a , b - h = P[Aab

V ah-2 -l,b-1,

q hh , ~n.h, b l UJ ~ah- -l 2, b - 1

= qh, Oh~t~'a,bllDh-2*jal,b-1

h--2 Va-l,b-1,

V ah-2 _l,b_l~

C ah, b ]

R h~.b]

R ha,b] =

= Qh(O,h),

Qh(o,h),

V ah- l ,:b - l ]

Va_l,bh-2 1]

x p[Dh~b~_, = Qh(O), Vafibl_l] =P[Aa-h,b --

-~

I"~h I r l h - 2 q0,'~a,bl'-'a-,,b-1

=

Qa,h b(1, • .

. ,

h-l),

1 V ah-2 - , , b - , ] P [ D a , h b-,

=

Qh(O),

h 1 Va,b_l]

h 1 h-I , b-1 = Qh(O), V~,b_,] P [ D . _h l- ,Ib = Q h ( h ) , V~-~,b] x P [ D a h-1 h-2 P[ Da_l,b_ 1 = Qh(O, h), vah-21,b_l]

The translation in terms o f a posteriori probabilities leads to (9).

[]

The lattice recurrence (9) can be rewritten in terms of j:h as follows: ~;a,b[Qh]h

h-I

h

o~2h-1

h

= s ~ ; a _ l , b [ Q (h)]3Ua,b_l[Q (0)] h-2 h 3ra-,,b I[Q (0, h)]

(10)

Like in (7), N is a normalization factor which can be ignored as we shall see in Section 4. The p a r a m e t e r and (h h; Ca,h b) are assumed to be h can be seen as a characteristic distance at which (Aa,b-h,"Ra.b) h h--2 . h--2 stochastically i n d e p e n d e n t conditionally u p o n (Da_l.b_l, Va-l,b-1). Therefore, we expect that the higher Signal Processing

71

V. Lacroix / Pixel labeling in a second-order Markov mesh

b ( ol)

v h

a,b "~

h-2

Va-l,b-I

(row)

II

I

Fig. 9. D e c o m p o s i t i o n o f V~. h.

h is, the weaker the hypothesis becomes. On the other hand, we see from (8) that the higher h is, not only do the steps leading to the final computation increase, but also the complexity of the sums involved at each step. More precisely, the complexity of the algorithm per pixel is of order O(c 2h+1) as one can see from the decomposition relation (8) needed to compute ~h-1, where c denotes the number of label classes. Further insight [8] into the case h = 2 enables a reduction of complexity from O(c 5) to O(c3). This reduction of complexity seems to occur for h = 2 only. In summary, the complexity of the algorithms generated with h = 1, h = 2, and h -- 3, will be respectively O(c3), O(c3), and O(c7). Therefore, the passage from h = 1 to h = 2 is done without any additional cost; on the contrary, the passage from h -- 2 to h = 3 is very expensive in terms of complexity. This argument suggests to use the algorithm with h = 2 for a best performance versus complexity trade-off as it will be confirmed with the experimental results (see Section 5). Though costly and complex, but still linear, we have a theoretical algorithm in order to compute bJ that we might expect to approach its real value when h increases. p[ha~. = q] --a, vlx)l In the next section, we shall specify the complete algorithm, including the boundary conditions which were intentionally omitted so far.

4. The complete algorithm

4.1. Boundary conditions It this subsection we shall give the boundary conditions at the 'h-borders', that is, at each pixel (a, b) belonging to row zero up to row h - 1 or to a column zero up to column h - 1. For these pixels, it is not possible to consider a triangular neighborhood containing h diagonals and only the functions ffi such that i<~ min{a, b} are defined (see Fig. 10). Therefore, the local decomposition relationships (8) can be used only up to function ~:i where i = min{a, b } - 1 and another lattice recurrence is needed in order to initialize the process. The initialization is slightly different whether a < b, a = b, or a > b. We shall label these different cases using their geometric position in the array (see Fig. 11); the boundary pixels such that a < b, a = b, and a > b are respectively called row-border pixels, diagonal-border pixels, and columnborder pixels. Vol. 12, No. 1, January 1987

V. Lacroix /

72

b

° ° °

.

.

Pixel labeling in a second-order Markov mesh

0

...

n

.

b

0

...h-1

b

..°

•-':"

. . . .

°

i

i

:

i:!ZiiiiEi i:i:i DI

i...!

E

* * i

(b) at x: existence of 50 and 5rl only

(a) at ×: existence of :v0 only

F i g . 10. E x i s t e n c e

of functions

i i i

i:i!i 57 ,°°

at x: existence of (c) ~r0,..., 7h-1 only

at the borders.

Along a row-border, the initialization is then made at ~ = ~ . Noting that Da+l a,b : Da, b-1

and

a

lla+l --,,,b

=

a Va, b-l,

we use (8) in order to initialize the recurrence: - l - i - R~+, ~:,b-I( O'*)Po°IR"÷'PQ~(D ' ~ ) ) , o

(11)

~::,b( Q~) - Xa

where 2¢"a= •

~+1 :~a,b-t(Q ~ ,a )PQ'°IQ'°+'PQ'"(Da,,,(,0 b ).

O,aQ ,

0

1

..h-I

o

[r

,

D[r

h-1

r

--+

b

N

r

r

r

r

r

r

r

r

r

r

r

r

r

r

r

r

r

r

r

r

r

r

r

r

r

r

r

r

¢= h - b o r d e r s

r: row-border pixel d: diagonal-border pixel c: column-border pixel

<= Inside array

M F i g . 11. h - B o r d e r s Signal

Processing

and inside array.

V. Lacroix / Ih'xel labeling in a second-order Markov mesh

73

Along a diagonal-border, the initialization is m a d e at .~-,,=b. In this case, we note that o.b = D aa- l , b - I Da+l

and

V~.~l = Va-j.h-1, a

and we use (8) once m o r e for the initialization: 1

J;".,b( O a) = - ~ ~+, ~;]_,,b_,( o " ) eQ.iRO+'poo( O~a!~)), where ./¢'~

=

(12)

Y. Y~ ~¢ea-l.b-l(Q ~. . . )Po,ots . .... pQ,°(Do.b(x) ). Q,a R,a+l

Finally, along a column-border, the initialization m a d e at ~ = b leads to the equation b

b

1

~;..b(Q )=~--g Y..~]-I.b(Q°)Po~IRh+'P&(Db!~)),

(13)

Rb +1

where N b = E

E

Q,b R,b+I

b ,b )Po'blR .... Po'~(Da,b(~) J;~-~,b(Q b ).

4.2. Special case The hypothesis for h = 1 is slightly different because of the decomposition of V~.b (see Fig. 12),

V'~,b = Ra,1 bk9 vOa-l.b-I ~ Cla.b. We introduce another hypothesis in order to find a relationship similar to (9) for h = 1: if (a~-,,b;R'~,b) and (A._~.b; C~..b) are stochastically i n d e p e n d e n t conditionally to V°~_t,b_~, then the following equality holds:

P[O'..b =

Q'I v'.,~] =

N P f D ° = Q'(O)iV°~-"b]P[D°'b ' = O'(1)l v°~'b-']

(14)

P[ V°._,.~_,] where N - P[ V~-"b]P[ V~ b-,]

P[ v~,d

The p r o o f of equation (14) is similar to that o f (9).

a,b

Va-l,b-1 o

Cal,b

Aa-l,b

R~,b

[,Xa,b-,

Fig. 12. Decomposition of Vta.band hypothesis for h = 1. Vol. 12, No. 1, January 1987

V. Laeroix/ Pixellabeling in a second-orderMarkov mesh

74

4.3. Pseudo-computer form of the algorithm In Section 3, we assert that the normalization factor N in (7) or in (9) is not necessary. In order to prove this assertion, let us define

~;*h( Rh) ~--~ ,~h( Rh).

(15)

The knowledge of $:h is not necessary to compute ~o; it turns out that sr.h is sufficient. Indeed substituting (15) into (8) yields

~7;h_l( Qh_l) = - ~ 1

~R Nl~*h(Rh)pQh 'IRhpQ~-t(xh-1)'

where M h-'= o,~, R~ N.~'*h(R'h)po,"-'I.,"po,"-'(xh-'), so that the simplification by N leads to

g~h-,(Qh-,)= where

1

~ ~'*hPQ"-'IRhpQh-'(xh-')'

~ch-, = 0,~_' ~

(16)

~-*h(Rh)pQ,h_,lRhpQ,h_,(xh-,).

The following algorithm summarizes how all derived formulae are used in order to compute S ~ ( q ) = P[;ta.b = ql V°a.b] for all q's and for all a, b of the VM.N array. The array VM.N is scanned as follows: • Initialize: start at the upper left pixel, at row 0, column 0. • Go along row 0: scan {(a, b) : a = 0, 0 < b ~< N}. • Go along column 0: scan {(a, b):b=O,O
-

if row-border pixel: min{a, b} = a # b; use (11), if diagonal-border pixel: min{a, b} = a = b; use (12), if column-border pixel: min{a, b} = b # a; use (13), - f o r j = 1 to min{a, b} use (8) in order to compute ~min{a,b}-j(Qmin{a,b}-j) f o r all Q's, so that 5~'°(Q °) •

•

•

is now known at the 'h-borders'. • If the pixel (a, b) is such that min{a, b}~h: compute ~,-.h for all Q's; use (15) and ((14) or (10)), - compute ~;h-l(Qh-l) for all Q's; use (16), - for j = 2 to h use (8) in order to compute ~;h-j(Qh-j) for all Q's, so that ~O(QO) is now known inside the VN,N array. To each choice of h corresponds an algorithm. If h = 2, the F - G - H algorithm is generated. -

Signal Processing

V. Lacroix / Pixel labeling in a second-orderMarkov mesh

75

The final step is a classification process where the pixel label is chosen according to the following rule: /)~,b = Q0 for P[D°,b =

Q°I

= mRaoXP[D°,b= R°I V°~,b]),

where/~o denotes the estimate of the true D °. (Note that D o is the zero diagonal, so that it contains the pixel (a, b) only, and R ° represents only one value of the label. The rule above is therefore exactly the decision rule presented in Section 1.)

5. Experimental results The application of Markov meshes to real images necessitates two phases: the developing and the learning phases. The developing phase consists in elaborating a technique able to accurately handle images exactly corresponding to the model. The learning phase which makes uses of the developing phase, comprises the process which fits the parameters of the model to the real image. This paper is concerned with the developing phase only because the learning phase forms a major topic in its own right and will be the subject of a subsequent paper [9]. The experimental results have two aims. The first one is to show how well indeed the proposed algorithm recovers an image simulated by a Markov mesh model from its noisy version. The second aim is to show the supremacy of the F - G - H algorithm. This means that not much could be gained in choosing a neighborhood larger than two diagonals, i.e., h > 2. In order to realize these aims, we need a practical application for which the complexity of the algorithm for h = 3 is not prohibitive. This is why the algorithm is applied to the restoration of artificial binary images corrupted by white noise, like in [17]. The images used in our experiments were produced by computer simulation of a first order two-state Markov mesh. Figs. 13(a) and 14(a) show both a simulated Markov mesh; each one corresponds to specific transition probabilities. White noise was added to the simulated Markov mesh and the image thus generated served as input data to several algorithms. The algorithms have to assign a label to each pixel. The label of a pixel should correspond to its 'color', black or white, in the noncorrupted image. The feature vector, which has as only component the intensity of the pixel in the noisy image, contains all the available information at the pixel. The following algorithms were implemented. First, the 'block constraint' method [14] based on local Markov dependence has been used as a reference for comparison because it appeared to be the best among various approaches investigated in [14]. This algorithm maximizes the likelihood in small regions within the image. Specifically, it chooses the label A according to the following rule: Aa.b q for K¢b(q) =

=

maxKa.b(r),

where

Ka, b( q) = ~, P[ r, s]p~(x.-l.b)p~(Xa,b-l)Pql~,spq(X~.b), r,$

and P[r, s] is the probability of occurrence of the event {Aa-l,b = r ; ) t a . b - i = S} over the entire image. The next algorithms are three applications of the algorithm described in this paper, each one corresponding to a neighborhood of h diagonals, with h = 1, h = 2, and h = 3, respectively. Vol. 12, No. 1, January 1987

V. Lacroix / Pixel labeling in a second-order Markov mesh

76

5.1. Input parameters The input parameters in our experiments are the boundary conditions (initial state of the upper left pixel, the transition probabilities for the first row and column), the transition probability matrix and, finally, the signal and the noise variances. The upper left pixel was arbitrarily decided to be black. The transition probabilities for the first row and column were estimated so that they fit the assumed model as well as possible. According to Besag [4], in a Markov mesh, the boundary conditions have little influence on the inside image as Markov meshes do not suffer from long-range order, and therefore they will not be considered as critical parameters in the experiments. The transition matrix of a first order two-state Markov mesh is a four by four matrix and may be represented as a point in a four-dimensional unit cube centered about (0.5, 0.5, 0.5, 0.5).

5.2. Results as function of extracted parameters Following Devijver [8], we characterized an experiment by the SNR and entropy of the input image. We used the between-class convariance matrix given in [10], applied to the two-class problem as a definition of the SNR: S N R = PoPl( (tZo- lZ,)/ o') 2, where P0 and P~ are, respectively, the probability of occurrence of black and white pixel, while/Zo and /z~ correspond to the means of the signals, and tr is the noise variance. As a second parameter, we used the entropy (see [12, 2] for the definition). The lower the entropy is, the more the image contains Markovian information. Therefore, experiences at low entropy better correspond to the supposed model than experiences at high entropy; but the latter might give an idea of the behavior of our algorithms in other contexts. A typical experiment is shown in Fig. 13. The entropy of the markov mesh in Fig. 13(a) is equal to 0.05, its noisy version in Fig. 13(b) corresponds to an SNR equal to 0.75, and the recovered images can be found in Figs 13(d), (e), (f). Fig. 14 shows a similar experiment on another simulated Markov mesh of 0.13 entropy. These two figures are quite illustrative of the power of the proposed algorithm but they do not provide quantitative results. Therefore, we designed two sets of experiments in which the error rate is computed in function of a parameter: the SNR in the first, and the entropy in the second. The noise-free image of Fig. 13(a) was used in the first set of experiments and several noisy versions corresponding to different values of SNR served as input images. Fig. 15 shows how well our algorithms recover the noise-free image for low SNR, compared to the block-constraint method. As announced earlier, we also see that there is much advantage in passing from h = 1 to h = 2 because of improvement obtained without any cost of complexity. On the other hand, there is less sense in passing from h = 2 to h = 3 because the complexity of the algorithm increases from O(c 3) to O(c 7) while little improvement is obtained. In the next set of experiments, all parameters were kept constant except for the transition matrix, visualized as a point inside a four-dimensional cube. Noting that all vertices of this cube correspond to minimum entropy (E = 0), while its center corresponds to maximum entropy (E = 1), the variation of the matrix was made by choosing points along a line joining a vertex to the center. The Markov meshes so generated are quite different one from another and they are not displayed except for the first one shown in Fig. 13(a). The results of these experiments are summarized in Fig. 16 and show that our algorithm is quite effective, especially at low entropy. In these experiments, the SNR was kept to 0.75. Once more, Signal Processing

V. Lacroix / Pixel labeling in a second-order Markov mesh

77

Fig. 13. (a) Input image, (b) noisy image, (c), (d), (e), (f) recovered images.

these results suggest to use the F - G - H algorithm (h = 2) for a g o o d compromise between complexity and performance. It has to be noted that, in the first set o f experiments, a different transition matrix gave similar qualitative results but rather different quantitative ones. By qualitative results we mean a better recovering at a low Vol. 12, No. 1, January 1987

78

V. Lacroix / Pixel labeling in a second-order Markov mesh

Fig. 14. (a) Input image, (b) noisy image, (c), (d), (e), (f) recovered images.

SNR than at high S N R and a better recovering of our algorithms compared to the block-constraint algorithm. In the second set of experiments the same remark is to be made when a different starting transition matrix is chosen: our algorithms are performing less when the entropy increases and, for all values of the entropy ~0.5, are performing more than the block-constraint algorithm. SignalProcessing

V. Lacroix / Pixel labeling in a second-order Markov mesh

79

0 0

[]

I

0

0 i'N 0

=

o=

~= ÷=

block-constraint h=l h=2 h=3

o rye, L

0

o

I. I. U.l 0 0

0 0

0 0 0

0.0

0.5

1.0

1.5

2.0

2.5

SNR Fig. 15. E r r o r r a t e as f u n c t i o n o f S N R .

6. Discussion

The algorithm has two parts. The first part is the computation of the function ~h,b associated to the diagonal Dh.b located at the boundary of the considered past neighborhood of the pixel (a, b). The second part consists of the computations of all functions associated to the next diagonals Da.bh-1,..., Da.b.° The first part is based on a well-defined hypothesis. Is this hypothesis reasonable? It is well known that, in the Ising Model, a noncausal two-state Markov field, pixels far apart are indeed correlated [15]; therefore, in such a model, our hypothesis would not be correct. However, according to Kanefsky [14], the Markov mesh model results, at least in some cases, in a spatial correlation that is a function of Euclidian distance. Besag [4] goes further: Markov meshes do not suffer from long-range order, that is, the colors of pixels arbitrarily far apart are not positively correlated. So, we can say that afortiori the conditional independence assumed in our hypothesis will be reasonable as h increases. The second part is perfectly rigorous as it was shown in the proof of the local decomposition formula. However, from a computer point of view, one could argue that an error propagation can easily occur in this type of algorithm where the computation of each function is based on the previous one. But, the functions in the decomposition formula are weighted by factors depending on local information so that Vol. 12, NO, 1, January 1987

80

V. Lacroix / lh'xel labeling in a second-order Markov mesh 0 0

d

L-d 0

I. t.

uJ

~/

~

hb~ck-c°nstraint

~ = h:2 += h=3

0 0

0 0

d 0.0

I

0.1

I

i

0.2 0.3 Entropy

I

0.4

0.5

Fig. 16. Error rate as function of entropy. error propagation has not much meaning in this context. The important thing is that a past neighborhood is taken into account in order to label the pixel (a, b), so that a large set of feature vectors is considered for each pixel. Another computational problem may occur as h increases. The diagonal Dha,b contains more and more labels so that the computation of ~ , h for all Qh's may cause underflow especially if the *h set of possible labels is large. In such a case, ~a,b should also be normalized. This argument along with the argument of complexity, suggests that, for practical purpose, one uses the algorithm with h = 2, which is the F - G - H algorithm. We may conclude the discussion by noting that in the literature we already saw pretty good results though the considered neighborhood was small [14], and we might expect our algorithm to give even better results in most cases.

7. Conclusion and comments

We have developed a linear-time algorithm using a second order Markov mesh model in order to label the pixels of a given image. The method is nonrecursive. The complexity of the algorithm depends on the size of the neighborhood considered in order to label the pixel. The algorithm has two parts. The first Signal Processing

V. Lacroix / Pixel labeling in a second-order Markov mesh

81

one is based u p o n an independence hypothesis, for which we gave several arguments, and makes use o f a lattice recurrence relation. The second part purely involves exact calculations summarized in the local d e c o m p o s i t i o n formula. The very encouraging experimental results suggest to use the algorithm with a n e i g h b o r h o o d containing two diagonals which corresponds to Devijver's F - G - H algorithm. Further work can be done in the developing phase. The algorithm assumes that the transition matrix is given; therefore, it would be interesting to analyze the sensitivity o f the algorithm to some changes in the matrix. A n o t h e r suggestion would be to change the first part o f the algorithm, that is, to use another formula for the c o m p u t a t i o n o f ~ h b or to set it arbitrarily. This algorithm can also be generalized to higher order M a r k o v meshes. M u c h attention is actually devoted to the learning phase, a difficult topic o f critical importance for the application o f M a r k o v meshes to real world image processing. We expect the M a r k o v mesh models to perform well in image segmentation, texture analysis, and synthesis. The labels to assign in an image to a segment should correspond to identification numbers o f region. As the complexity o f the algorithm depends on the n u m b e r o f available labels, the input image should present a reasonable n u m b e r o f regions. The feature vector of a pixel will naturally be associated to its current grey level. Therefore, the input image should be c o m p o s e d o f rather h o m o g e n e o u s regions. Works similar to [7] and [6] could be realized with M a r k o v meshes for texture analysis and texture synthesis.

Acknowledgment I wish to express all m y gratitude to P. Devijver who introduced me to the field, gave me m u c h theoretical help and e n c o u r a g e m e n t s while supervising this work. Thanks are also due to M. Dekesel for his p r o g r a m m i n g assistance and to P. Delsarte for his suggestions and c o m m e n t s on the structure and the notations in a first draft.

References [1] K. Abend, T.J. Harley and L.N. Kanal, "Classification of binary random patterns", IEEE Trans. Inform. Theory, Vol. IT-11, October 1965, pp. 538-544. [2] N. Abramson, Information Theory and Coding, McGrawHill, New York, 1963. [3] J. Besag, "Spatial interaction and the statistical analysis of lattice systems", J. Roy. Statist. Soc., Set. B, Vol. 36, 1974, pp. 192-326. [4] J. Besag, "On the statistical analysis of dirty pictures", Paper read at the SERC Research Workshop on Statistics and Pattern Recognition, Edinburgh, July 1985. [5] R. Chellappa and R.L. Kashyap, "Digital restoration using spatial interaction models", IEEE Trans. Acoust., Speech, Signal Process., Vol. ASSP-30, No. 3, 1982, pp. 461-472. [6] R. Chellappa and R.L. Kashyap, "Texture synthesis using 2-D non-causal autoregressive models", IEEE Trans. Acoust., Speech, Signal "Process.,Vol. ASSP-33,No. 1, 1985, pp. 194-203. [7] G.R. Cross and A.K. Jain, "Markov random field texture models", IEEE Trans. Pattern Analysis Machine D~telligence, Vol. PAMI-5, 1983, pp. 25-39.

[8] P.A. Devijver, "Probabilistic labeling in a hidden second order Markov mesh", in: E.S. Gelsema and UN. Kanal, eds., Pattern Recognition in Practice, II, North-Holland, Amsterdam, 1986. [9] P.A. Devijver, "Hidden Markov models for speech and images", in P.A. Devijver and J. Kittler, eds., Pattern Recognition Theory and Applications, Springer, Berlin, to appear. [101 P.A. Devijver and J. Kittler, Pattern Recognition: A Statistical Approach, Prentice-Hall, Englewood-Cliffs, NJ, 1982. [ll] S. Geman and D. Geman, "Stochastic relaxation, Gibbs distributions and Bayesian restoration of images", IEEE Trans. Pattern Analysis Machine Intelligence, Vol. PAMI-6, November 1984, pp. 721-741. [12] L.N. Kanal, "Markov mesh models", Comput. Graphics Image Process., Vol. 12, 1980, pp. 371-375. [13] M. Kanefsky and C. Fong, "Predictive source coding techniques using maximum likelihood prediction for compression of digitized images", IEEE Trans. Inform. Theory, Vol. IT-30, No. 5, September 1984, pp. 722-727. [14] M. Kanefsky and M.G. Strintzis, "A decision theory approach to picture smoothing", IEEE Trans. Comput., Vol. C-27, January 1978, 32-36. Vol. 12, No. 1, January 1987

82

V. Lacroix / Pixel labeling in a second-order Markov mesh

[15] R. Kinderman and J.L. Snell, Markov Random Fields and their Applications, Amer. Math. Soc., Providence, RI, 1980. [16] A. Rosenfeld and A.C. Kak, Digitalt~'ctureProcessing, Vol. 1, 2nd ed., Academic Press, New York, 1982, Chap. 7, p. 312.

Signal Processing

[17] O. Wolberg and T. Pavlidis, "Restoration of binary images using stochastic relaxation with annealing", Pattern Recognition Lett., Vol. 3, December 1985, pp. 375-388. [18] J.W. Woods, "Two-dimensional discrete Markovian fields", IEEE Trans. Inform. Theory, Vol. IT-18, March 1972, pp. 232-240.