Unilateral approximation of Gibbs random field images

Unilateral approximation of Gibbs random field images

CVGIP: GRAPHICAL MODELS Vol. 53, No. 3, May, pp. AND 240-257, Unilateral IMAGE PROCESSING 1991 Approximation of Gibbs Random Field Images ...

2MB Sizes 0 Downloads 45 Views

CVGIP:

GRAPHICAL

MODELS

Vol. 53, No. 3, May,

pp.

AND

240-257,

Unilateral

IMAGE

PROCESSING

1991

Approximation

of Gibbs Random Field Images JOHN GOUTSIAS

Deparlment

of Electricul

and Computer The Johns Hopkins

Received September

Engineering, University, 12,

1989;

important, since, in some cases, various problems related to the of images via Gibbs random fields may be simplified

by

restricting our interest to a special and convenient class of Gibbs random fields, the class of mutually compatible Gibbs random fields. A general analysis is presented, which results in some

interesting conclusions about the desired relationship. We consider two approaches to the approximation problem. The first approach, which is more accurate, can be easily applied to the case of binary images. For nonbinary images this approach becomes cumbersome, since it requires the solution of a large system of nonlinear equations. The second approach, which is simpler but less accurate, is based on the direct relationship between the pa-

rameters of a general Gibbs random field and the parameters of a mutually compatible Gibbs random field. Simulation examples demonstrate VariOUS aSpeCb Of Our ZiUalySiS. 0 IY!N Academic PW, IW.

I.

INTRODUCTION

An important problem in image analysis and computer vision applications is the development of mathematical, or statistical, models for comprehensive image description. A popular statistical solution to this problem is via random fields, and specifically, via Gibbs (or Markou) random jields (GRF) [l-3]. The popularity of these models is due to the fact that they can effectively describe contextual characteristics of images in terms of a few parameters. This approach to image modeling has been adopted by many researchers for solving a variety of problems such as image smoothing and segmentation [4-81, image restoration [9, 101, texture synthesis [I 1, 121, motion analysis [13, 141, and other computer vision tasks 115, 161. The wide applicability of these models is usually offset by a number of theoretical and computational drawbacks. For example, in order to analyze a Gibbs random field by means of a digital computer, we need a procedure which produces a sample from a general Gibbs random field. 240 1049-9652191

Copyright All rights

$3.00

0 1991 by Academic Press, Inc. of reproduction in any form reserved.

Cornmunicutions 21218

accepted

4, 1990

September

Luhorutory.

The direct solution to this problem is not possible in general, and computationally inefficient Monte Carlo simulation algorithms have been deviced to deal with this problem [17]. Additionally, in many important application areas (e.g., ecology, image analysis, computer vision, neural networks) the Gibbs random field has potential for modeling two-dimensional spatial data and, consequently, is worth studying in a purely statistical context. It is well known that a Gibbs random field is parametrized via its neighborhood structure, yielding exponential families in canonical form. Fitting a Gibbs random field model to real data requires the estimation of its parameters from the given data. The solution to this problem via maximum-likelihood parameter estimation seems, at a first glance, straightforward, but unfortunately, it is quite difficult since there is no known technique for the computation of the likelihood function. Moreover, Gibbs random fields may exhibit phase transitions and long-range interactions, thereby creating identifiability problems [lS]. Although a general solution to these problems is under investigation, a specialized solution can be obtained, in many important situations, by restricting our interest to a special class of Gibbs random fields, the class of mutually compatible Gibbs random fields [ 19, 201. It is known that a mutually compatible Gibbs random field is equivalent to a Markov mesh [20-221 or unilateral Markov random $e/d [23, 241. These models have been extensively used as crystal growth models in crystallography [25, 261, for image smoothing and segmentation [27-301, and image coding [31], and as we shall show in this paper, they are directly related to general Gibbs random fields. The purpose of this paper is to exploit the mathematical relationship between a general and a mutually compatible Gibbs random field image, to investigate the possibility of approximating a general Gibbs random field by a mutually compatible Gibbs random field, and to study possible consequences of such an approximation. Our goal is not to replace a general Gibbs random field by a mutually compatible Gibbs random field, but rather to exploit and demonstrate their relationship. We show that, in many interesting cases, a general Gibbs random field can be well approximated by a simpler mutually

We exploit the relationship between a general Gibbs random field and a mutually compatible Gibbs random field image, we investigate the possibility of approximating a general Gibbs random field by a mutually compatible Gibbs random field, and we study the consequences of such an approximation. This study is modeling

Image Anulysis und Bultimore, Maryland

GIBBS RANDOM

compatible Gibbs random field (in which cases a mutually compatible Gibbs random field can be used directly for image modeling), whereas, in other cases such an approximation is poor (in which cases a general Markov random field model is more appropriate). The approximation is shown to be parameter value-dependent and the parameters responsible for the quality of such an approximation are identified. We believe that this study is important, since, in many practical situations, various problems related to the modeling of images via Gibbs random fields can be simplified substantially by restricting our interest to the special and convenient class of unilateral Gibbs random fields. Our results allow us to establish a clear picture of the advantages and disadvantages of using unilateral versus general Gibbs random field models. This paper is an extended version of the work in [32] and is structured as follows. In Section II we briefly introduce the elementary theory of Gibbs random fields, whereas, in Section III we briefly summarize basic facts about the behavior and properties of mutually compatible Gibbs random field images. Except for the fact that a mutually compatible Gibbs random field is a special case of a general Gibbs random field, no other information about their relationship is generally known.’ In Section IV a general analysis is presented which results in important conclusions about the desired relationship. We consider two approaches to the approximation problem. The first approach can be applied easily in the case of binary images. For nonbinary images this approach becomes cumbersome, since it requires the solution of a large system of nonlinear equations. The second approach, which in principle is a slight modification of the first approach, is an alternative simpler approach based on the direct relationship between the parameters of a general Gibbs random field and the parameters of a mutually compatible Gibbs random field. We introduce the concept of the relative entropy as a measure of the distance between any probability measure and the Gibbs measure. We show that in some important cases, the second approach results in a particular mutually compatible Gibbs random field whose probability measure minimizes the relative entropy, among all other probability measures from the class of mutually compatible Gibbs random fields. In Section V we show experimentally that, in many interesting cases, a general Gibbs random field can be well approximated by a mutually compatible Gibbs random field. We also show that, by initializing the Monte Carlo simulation algorithm with the computed approximation, we can speed up the simulation of a Gibbs random field by a factor of up to about 10. An experimental analysis is also performed in order to study the sensitivity of the proI A possible exception to this remark are Refs. [26, 33, 341, where a specialized analysis toward a better understanding of this relationship is presented. Our work constitutes a natural extension of this analysis.

241

FIELD IMAGES

posed approximation to the variation of certain modeling parameters. We conclude the paper by summarizing our results in Section VI. II.

GIBBS RANDOM

FIELDS

It is well known [1] that, a second-order2 GRF H, which is defined over a two-dimensional rectangular lattice AMN = {(i, j): 1 5 i 5 M, 1 5 j 5 N}, is statistically characterized by the joint probability measure Pr[H

= h] = G

1

exp

- + Wh)J,

0)

where

Z CPF

=

c

exp

-

i

f

(lb)

U(h)}.

h&,

is the partitionfunction. In Eq. (I), T is the temper&we, which controls the degree of “peaking” in the probability measure Pr[H = h], whereas, U(h) is the energyfunction, given by3 + Vij”‘(hij,

(J(h) = - 5 5 [Vfj’(h;j)

j=, j=,

+ V!j’2’(hij, hi,j-‘)+

hi-l,j)

V:j’3’(h;j, hi-‘,j-‘)

+ V!,2’4)(hi-‘,j, hi,j-‘) + Vi:“‘(hij,

hi-‘,j, hi-‘,j-‘)

+ Vij’2’(h[j, hi-‘,j, hi+1) + Vij’3’(hij, hi-l,i-l)

hi,,-,)

+ V::.4)(him’,j, hi-l,j-l, + Vj:‘(hij,

hi-‘,j,

hi,j-j)

hi-t,j-t 9hi,j-l)l,

(2)

in terms of the canonical potential functions (CPFS) Vij. The CPFs are characterized by the property that they equal zero if at least one of their independent variables equals zero [35]. If the random variable H(i, j), which is assigned at site (i, j), is discrete-valued and takes R (R 2 2) distinct values from a discrete ensemble EH = (4 ’ , $2, $R}, then cH in Eq. (lb) denotes the set of all possible RMN distinct realizations of H. It can be shown [1, 351 that the GRF H is a Murkou random jield (MRF) with neighborhood Nii, given by r The ideas presented in this paper can be generalized easily to include first-order, or higher-order, GRFs. For clarity of presentation, we shall base our analysis on second-order GRFs. 3 Capital letters denote random variables or random fields, whereas small letters denote their realization.

JOHN GOUTSIAS

242

= {(i + 1, j + I), (i, j + I), (i - 1, j + I),

Nij

(i + 1, j), (i - 1, j), (i + 1, j - I), (i, j - l), (i - 1, j - l)}.

(3)

In the following, we shall assume that the random field H is zero in {(i, j), i 4 0 or j 5 O}. This is clearly the case of a free boundary. An alternative definition for a GRF is given in terms of the local transfer function [ 19,201. Let us assume that the statistical contribution to the lattice system described above, due to various interactions among the random variables H(i, j), H(i - 1, j), H(i - 1, j - I), and H(i, j - I), is given by the homogeneous local transferfunction (LTF) v(hij, hi-r,j, hi-r,j-r, hi.j-1). In this case

where the partition

ZLTF

=

C IlEE,,

fi l=I

function Z,,

fi j=l

c(h;j,

hi-l,j-l,

hi+,).

III.

MUTUALLY

(4b)

It can be easily proved that the GRF H, whose joint probability measure is given by Eq. (4), is also a MRF with neighborhood Nij, given by Eq. (3). The following lemma, whose proof is reminiscent of the proof of the Hammersley-Clifford theorem [ 11, results in a useful representation for the LTF u. If 4, = 0, then, for any LTF a(x, y, z, w), there exists a unique expansion of the form LEMMA.

Pr[HA = hAI = k nn A

Y)XY

+ P4(Y9

Z)YZ

+

P2k

WhJ

+ PdY,

+

Ylk

Y, Z)XYZ

+

Y3b,

z, w)xzo

+

+

O)YW Y2(X,

+ Y4(Y,

+ a3(z)z + w(w)w

y,

P3(x,

z)xz

w(hij, h;-l,j, hi-l,j-1, hi+,),

(LJEA

where =

cfln hn

Plk

u(hij,

hi-l,j,

hi-l+,,

I

f

hi+,).

(i,jEA

The next theorem results in a necessary and sufficient condition for a GRF H to be a MC-GRF. Its proof can be found in [20]. THEOREM I. A necessary and sufficient condition for a GRF H to be a MC-GRF is that there exists a LTF Q-(X, y, z, w) such that

+ pb(z,o)zo w)xyo

z, o).Yzo

4x7 y, z, WI = cdx, y, z, 01, + 6(x, Y7 z, ~)XYZWl

FIELDS

DEFINITION 2. The GRF H defined over a rectangular lattice RwN is a mutually compatible Gibbs random jield (MC-GRF) if there exists a LTF (T satisfying Eq. (4), so that the restriction HA of H over any primary sublattice A of lattice RMN is also a GRF with joint probability measure

ZA

= exp i -) [r + a~(x)x + dy)y

GIBBS RANDOM

Let us now assume that the probability structure of a subset HA of random field H, restricted on a finite sublattice A of lattice AMN, is independent of the size or shape of A. Following [ 19, 201 we have the following two definitions.

r(x, Y, z, w)

+

COMPATIBLE

DEFINITION 1. A set of sites A is called the primary sublattice of lattice AMN ifA c RM,v, {(i, j): 1 I i 5 K, j = l} c A, {(i, j): i = I, 1 I j I L} c A, for some (K, L) with 1 5 K 5 M, I 5 L 5 N, and if (i, j) E A, then {(i, j), (i 1, j), (i - 1, j - I), (i, j - l)} C A.

is now given by

hi-l.j,

random variables H(i, j), H(i - 1, j), H(i - 1, j - I), and H(i, j - 1). For this reason, we shall call them the interaction parameters. The relationship between these parameters and the corresponding LTF is summarized in Appendix A (see Eq. (Al)), which also provides an expression of the CPFs Vij in tbrms of these parameters (see Eq. 642)).

(5)

The assumption that $r = 0 is not restrictive, since we can always modify the values of EH to obtain 4, = 0. Without loss of generality, we shall assume that 4, = 0 throughout. Parameters r, (Y, p, y, and 6 quantitatively characterize the amount of statistical interaction between

(64

where ,,& G,

Y, z, 0) = 1,

for every triplet (y, z, w) E Ei, independent of (x, y, z, w).

(6b)

with c being a constant

GIBBS

RANDOM

FIELD

+ P:“‘k zbx +

For a MC-GRF the computation of the partition function is trivial. Indeed, from Eqs. (4) and (6) it is easy to prove that

Z LT/T = fi

fi

i=l

J=I

C 2

T(U,

hi-l,j,

h;-l,j-l,

hi,j-1)

= CMN.

From Eqs. (4a), (6), and (7) we see that the joint probability measure of a MC-GRF H(O) is given by = h(o)]

= fi fi Pd”)[hi,O’ / h)!,,j, i-1 j-1

hi?,,j-,,

hjyj-,],

@a)

where PtJ”)[h;;’

1hjo’,.j, hj!!!,+,

, hi;;-, ]

= T(h!?’,J ’ h!O’ h!O’ __ h!‘!_ ) 1 ],J’ I

1.J

I’

I../

@b)

1

is the conditional probability of H(O)(i, j) = hly’, given that H(O)(i - 1, j) = hr,,j, H(O)(i - 1, j - 1) = hlO),+,, and ~(o)(i j - 1) = h!‘!- . The simulation of a%&-GRF is a relatively easy task due to the special form of the joint probability measure, given by Eq. (8). The values of Q-(X, y, z, w) are specified such that Eq. (6b) is satisfied and a point-by-point simulation (e.g., in a lexicographic ordering) is performed easily by using one of the algorithms in [36]. IV.

+ pp(z,

w)zw

+

y:o’(x,

z,

+

6 (O’(x

p:“‘(Y,

Z)YZ + pg”‘(Y,

+ yfyx,

w)xzw

+

UNILATERAL

APPROXIMATION

OF GRFs

We shall now examine the relationship between a general GRF and a MC-GRF in terms of their interaction parameters. Let us assume that a MC-GRF H(O) is statistically characterized by a joint probability measure Pr(O)[H(O) = h(O)], given by Eq. (8). Let us also assume that a general GRF H is statistically characterized by a joint probability measure Pr [H = h], given by Eq. (4). According to the lemma

y,

z)xyz

+

z,

w)yzw

y:“)(y,

3 y 7 z 3 w)xyzo]

O)YO

y:o’(x,

y,

0)xyo

(9)

17

NEE”

(7)

p,Co)[H’o

243

IMAGES

where the interaction parameters cr(‘), p(O), y(O), and 6”’ equal zero if at least one of their arguments equals zero. The relationship between the interaction parameters r(O), a(‘), /3(O),y”‘, and disco’and the corresponding LTF 7(x, y, z, o) is similar to the relationship between the interaction parameters r, (Y, p, y, and 6 and the corresponding LTF u(x, y, z, w) (see Eq. (Al)). From Eqs. (4) and (8) we conclude that a MC-GRF is a GRF with LTF 7(x, y, Z, w), a partition function which equals one, and CPFs V;:(x), vg;;‘c& Y), vf;:‘cx, w), Vfl:‘(x, z), vgy, o), V&1’<& y, z), V~~‘(X, y, w), Vs;:‘(x, z, o), Vf$y, z, w), and Vp:(x, y, z, w). The CPFs are related to the interaction parameters r(O), (Y(O),p(O), y(O), and 6(O) with equations similar to Eqs. (A2). We now consider the following approximation problem: Given the interaction parameters (Y,p, y, and 6 of a general GRF H, compute the interaction parameters r(o’, am , p(O), y(O), and a(” of a MC-GRF H”’ such that pr(o)[H’o’ = h] = Pr[H = h]. We shall consider two approaches to this problem. The first approach (Approach A) is easily applied in the case of a binary random field (i.e., when R = 2). For R Z- 3 Approach A becomes quite cumbersome, since it requires the solution of a large system of nonlinear equations. The second approach (Approach B) is an alternative simpler approach which is based on a direct relationship between the LTF of a general GRF H and the LTF of a MC-GRF H(O). A.

Binary Random

Fields (Approach

A)

In Appendix B we derive Eqs. (B2) and (B3), which form a set of necessary and sufficient conditions on the interaction parameters r(O), a(‘), p(O), y”‘, and 8”’ such that the LTF 7(x, y, z, w) satisfies Eq. (6b). If a MC-GRF H(O) is statistically equivalent to a general GRF H (i.e., if Pr[H”’ = h] = Pr[H = h]), then the following (R - 1)(R3 + R* - R) equations should be satisfied (see also Eqs. (4), (5), (8), (9), and (A2)): a(O)(x) + a f’(x) + do’(x) + a(O)(x) 4 I 3 =

a,(x)

+

(Yz(x)

+

(Y3(x)

+

(Y4(x), (104

= exp + [r(O) + a(O)(x)x ] + a(O)(y)y + a$“‘(z)z 2 i +

arp)(w),

+

p;O’(x,

y)xy

+

p’O’(x 2

’ w)xw

p’1”(&

y)

+

@“(y,

x)

=

p,(x,

Y>

+

p6(Y7

x),

(lob)

p:“‘(x,

w)

+

p:“‘(x,

w)

=

p&7

0)

+

P4k

WI,

(1Oc)

z)

=

/33(x,

z),

p:o’(x,

(104

244

JOHN

fl”‘tY, a) = PjtY, @I, 5 Y’p’k Y, z> = Ylb, Y, z),

GOUTSIAS

tloe) uw

y:"'(x,

Y,

WI

=

Y2b,

Y,

WI,

(log)

yytx,

z,

w)

=

y3(x,

z,

WI,

(10h)

rptu,

z,

w>

=

Y4(Y,

z, WI,

For the case of binary random field (i.e., for R = 2) Eqs. (12a) and (12b) become4

B1t1, 1) = [I + A(I)l[l

[I + A(l)Bdl,

(lOi)

x exp

and f3(O’(x3 y 3 z 3 w) = 6(x 1 y 3 z 2 w).

(loj)

The GRF approximation problem is now reduced to the problem of solving the system of (R - l)(R3 + 2R 2 + 1) + 1 equations (lo), (B2), and (B3) (i.e., the (R - l)(R3 + R2 - R) linear equations (10) and the (R - l)(R2 + R + 1) + 1 nonlinear equations (B2)), in terms of the (R - l)(R3 + R2 + R + 1) + 1 interaction parameters r(O), a!(‘), /3”‘, y”‘, and a(“‘. Although an exact solution to this over-determined system of equations is not possible, an approximate solution can be obtained by assuming that ff I”‘(x) = a,(x) 9 p:"'k Yy'k

y,

(114 tllb)

z)

=

/33(x,

z),

7.)

=

y,(x,

y,

z>,

(1 Ic)

y;"'k

y,

0)

=

y&G

y,

01,

(114

Y:O't&

23

WI

=

73(x,

z,

w),

(1 le)

and 8’O’(x, y, z, w) = 6(x, y, z, w).

(110

From Eqs. (11) and (B3) we see that parameters M-G z), TI(x,

r&,

Y, z),

Y, 4,

o) can be easily calculated. (B3b) we obtain

r&x,

z, ~1,

and

A(X),

&,

Y, z,

From Eqs. (lob), (B2g), and

&ej,A(u)%rwfi A(u) x

B1(x’ ‘) = &a,,

YP3G4,

BZb,

A(u)

x)r3(u,

y) &Et,

Y. xl

{~[~,(l,

1)

+

1) r,tl, 1, I)1 IIlL + A(lM1, 111 l)B3(1,

/36(l,

1;

[@1(x,

Y)

+

p6t.Y~

dlrq'].

B3(1,

1)

= 11+ A(lN1 + A(l)B,(l, 1)B2(1, 1) rltl, 1, 111 [I + AtI)B,tl, 1)1[1+ A(1M1, 111 x exp

{$32u1

1)

+

/34u,

(12a)

v:;T;j’(x, y) = vy’(x,

w) = by(x,

B,(x,

WI

wvx4

X,

x exp

A@MI(~,

i+

[P2b-,

x)

&a,,

WI

+

A(U)

P4k

4 We

0)

dlxw}.

(14b)

(14c)

4

= 2 uEEH

(14a)

w),

forl5i~M-l,l~j~N, x)B~(u,

(13b)

y),

forl5iiM,lSjSN-1, vg’(x,

~~(4

IH).

Parameters B, (1, 1) and B3( 1, 1) can now be computed by solving the nonlinear equations (13).5 These equations can be solved iteratively by employing a technique from [37]. In the case when R 2 3, the system of 2(R - 1)’ nonlinear equations (12) must be solved simultaneously in terms of the 2(R - 1)2 parameters B,(x, y) and B3(x, w). If a solution to this system exists, then its computation will be prohibitive, for large R; therefore, this approach is limited generally to the case of binary random fields. After the parameters B,(x, y) and B3(x, o) are calculated, the interaction parameters r(O), a\“(y), a:“(z), a:“(w), @y’(y, z), pi”(y, w), pd”(z, w), and yF’(y, z, w) can be calculated directly from Eq. (B2). Finally, the interaction parameters fiy’(x, y) and &“(x, W) can be calculated easily from Eqs. (lob) and (IOc), respectively. According to the previous analysis, Approach A results in a set of interaction parameters (r(O), p(O), y”‘, and 8(O), which satisfy the constraints (lob)-( lOd), (lOf)(IOh), and (1Oj). Constraints (lOa), (IOe), and (lOi) are not satisfied in general. From the previous remarks and from Eq. (A2) observe that, in general,

whereas, from Eqs. (lOc), (B2e), and (B3d) we have

x

(13a)

and

for i = M, j = N, X exp

I)])

VXJj(X, = Vjj’(x),

x)

A(N)B~(u,

+ A(IMl,

(12b)

shall assume, without loss of generality, that C& = 1. 5 Observe that by using Eq. (13b) in Eq. (13a) we obtain a nonlinear equation with respect to B,(l, 1). After this equation is solved, Eq. (13b) can be used to determine the value of BZ(l, 1).

GIBBS

RANDOM

V'2!3'(X 0.11 ' z) = V!2,3'(x r, 3 z) 3

for 1 I i 5 M, 1 5 j 5 N,

(144

V3!!‘(X 0,r.l ’ y ’ z) = v!?“(x ‘J 5 y 3 z) 9

for 1 5 i 5 M, 1 % j 5 N, v3:*yx O,lJ

)

y ) w) = v!?*‘(x LJ

y 3 0)



(14e)

3

for 1 5 i 5 M, 1 5 j 5 N,

(140

V3!3’(X O,lJ 7 z ’ 0) = VF3)(x IJ ‘9z w) 1

(14g)

and

245

IMAGES

GRF H is indeed a MC-GRF. This procedure may prove to be quite useful in modeling images via GRFs, since, given the CPFs, it provides a simple criterion for the type of GRF under consideration. To conclude this section we mention the following interesting property. Its proof can be found in Appendix C. PROPOSITION. A GRF with homogeneous CPFs can never be equivalent to a nontrivial MC-GRF.

B.

for 1 5 i I M, 1 5 j 5 N,

Multivalued

Random Fields (Approach B)

Let us now define the relative entropy (or, discrimination) D(nl]p) of the general Gibbs measure rr(h) = Pr[H = h], with respect to the probability Prc”)[Hco) = h] of a MC-GRF H(O), by

V4’..(x y z w) = V!?(x y z w) for 1 5 i % M, 1 I j I N, O,lJ

FIELD

)

)

IJ

7

7

9

3

measure p(h) =

3

(14h)

whereas Vpij(X,

# V!!‘(x) lJ ’ for (6 A E h-N - {CM, N)),

Vf;j’(x,

y) # vyyx,

y),

for 1 d i 5 M, j = N, vg;;yx,

w) # vy’(x,

(15b)

w),

for i = M, 1 I j 5 N, vfjy’(y,

(154

w) # vy4yy,

(15c)

w),

for 1 5 i I M, 1 5 j I N,

(154

and v3:q y z w) # v!334’(y z 0) o.lJ

)

)

,J

‘)

)

for 1 5 i 5 M, 1 I j P N.

WeI

This approximation is achieved at the expense of computing the solution of a system of 2(R - 1)2 nonlinear equations. Given the CPFs V, of a GRF H, we can easily check if the GRF H is a MC-GRF by using the following three steps: (a) given the CPFs Vij we calculate the parameters al(x), Plk Y), P2k WI, P3(x, 4, Y,(X, Y, 21, Y2k Y, WI, Y&X, z, w), and 6(x, y, z, o) from Eqs. (A2b)-(A2e), (A2g)-(A2i), and (A2k);6 (b) we compute the quantities A, B,, B2, B3, Fr, F2, F3, and A from Eq. (B3); and (c) we calculate the parameters a*(y), a3(z), cx4(w), p4(y, z), Ps(Y, 01, p6k d, and ydy, z, w) from Eq. W-9. If the calculated parameters satisfy Eq. (A2), then the given 6 From Eqs. (A2b), /3,(x, y) = (llxy)V~j(x,

where 7Tk = r(hk), pk = p(hk), and K = RMN [IT, 381. The relative entropy D(r](p) is a measure of the “distance” between the probability measures p(h) and n(h), for every h E xc, [38]. It is worthwhile to approximate a GRF H by a MC-GRF H(O) whose probability measure is given by Pr(O)[H(O) = h] = a(h), where B(h) satisfies the equation

(ARC), and (A2d) we have a,(x) = (l/x)V!&,(,(x), y), and pz(x, w) = (l/xw)V~~(x, co).

where PWC-oRFis the set of all probability measures p(h) which correspond to a MC-GRF. We have the following theorem. Its proof can be found in [17]. THEOREM 2. The probability measure e(h) which minimizes the relative entropy D(n(lp), among all possible probability measures p(h) which correspond to a MCGRF, is given by

B(h) = fi fi +(hij, h;-l,jT hi-l,,i-t, i=l j=l

hi.i-I),

(164

wjrere qx,

y, z, w) = +

8 ln &TF y7,o I aucf,

Yv z7 WI

1u(x,z,0)(16b) Y,

and

A yzw = ,&

[du;;;y~~w,]

u(u, y, z, w).

(16~)

Although the LTF +(x, y, z, w) given by Eqs. (16b) and (16~) cannot be calculated in general, for every LTF a(x,

JOHN GOUTSIAS

246

y, z, w) of a GRF H, there exists a LTF T&, y, z, o) of a MC-GRF Ho, such that 70(x, y, z, w) =

4x7

Y, z, w)

Cl, (T(u, y, z, WI.

v;.;.,i(x, = v;;‘(x),

Let us denote by P(),J = pO(hk), k = 1, 2, . . . , K, the probability measure of the MC-GRF Ho with LTF T&, y, z, 0). We have the following theorem. Its proof can also be found in [17]. THEOREM 3.

for (x, y, z, w) E Ei. In this case, conditions (lOd), (lOf)(lOh), and (1Oj) are satisfied, whereas conditions (lOa)(lOc), (lOe), and (lOi) are not satisfied in general. From Eqs. (18) and (A2) we have that

V$f’(x,

y) = ly(x,

y),

(19b)

for 1 5 i 5 M, j = N, v&;‘(x,

If

(19d

for i = M, j = N,

w) = q2yx,

w),

for i = M, 1 5 j % N, vgx,

for every (x, y, z, o) E E$, where U&Y, y, z, w) is the total number of elementary configurations (h;j = X, hi-1.j = y, /~-r.,~-~ = z, /~;.~-i = o) in a realization hk, then

z) = q3’(x,

dx,

(194

V3:!‘(x 0,lJ ’ y ’ z) = V!V)(x rJ 3”y z) for 1 5 i % M, 1 5 j 5 N,

We)

v’3!2’(x IJ 3 y ’ w) ’ O,l/ ’ y ’ w) = VF2)(x for 1 5 i 5 M, 1 5 j 5 N,

(19~1

for 1 5 i I M, 1 5 j I N,

According to the previous theorem, a general GRF with LTF (T can be approximated by a MC-GRF with LTF Q, given by Eq. (17b), provided that Eq. (17a) is satisfied. Although condition (17a) cannot be verified analytically, its validity, in some interesting cases, is demonstrated experimentally in the next section. From Eqs. (Al), the equations which relate the interaction parameters r(O), &v, pm, p, and 6(O) with the LTF Q-~(x,y, z, w) (which are similar to Eqs. (Al)), and from Eqs. (17b) and (A2), with 7(x, y, z, w) = T~(x, y, z, w), we obtain

and ‘g’ljCx,

Y) = Pl(X, Y),

(18b)

P:"'(&

WI = Pzk,

WI,

(18c)

z),

(1X4

z)

= Pdx,

Y'p'(X,

Y, z) = Yl(X,

Y, z),

(18e)

Y:o'(&

Y? 0)

Y, WI,

(1W

rl"'cx

= Yz(X,

3 z 7 WI = ydx

53z

w) 3

(1%)

and 8’O’(x 3 y 7 z 9 w) = 6(x ,179 y z w)

(18h)

(19h)

whereas

v;,‘,(x) # v;;‘(x), for vg;jyx,

y)

#

kj)

vyyx,

E AMN

-

W,

W,

(20d

y),

forlSi%M,15j
(18a)

/p(x,

Y, Z, w) = V1,4’(X, Y, Z, W),

for 1 5 i 5 M, 1 5: j 5 N,

Vff(x,

a(O)(x) = a,(x), 1

(1%)

YT z, w) (17b)

p:oYx,

z, w) = v;j~3’(x, z, o),

Y, Z? w)

lx ,,GL,,d~,

z),

for 1 9 i I M, 1 5 j I N,

v$;‘(x, f(x, y. z, 0) = 7()(X, y, z, WI =

(19c)

GObI

0) # v;;-2)(x, w), forlSiSM-l,lsj
V:T;j’(y, w) # vp4yy,

(2Oc)

w),

for 1 % i 5 M, 1 zs j 5 N,

(204

and Vf;i”‘(y, z, w) f v;jg4yy, z, w), for 1 5 i I M, 1 5 j I N,

WeI

in general. By comparing Eqs. (14) and (15) with Eqs. (19) and (20) we observe that the previous approaches to the approxi-

GIBBS

FIG. initialized

1. Nine realizations

mation problem result in a similar approximation of the CPFs, the difference being only in the approximation for the CpFs v!?i)(x, y) and V!?*‘(x, w). Finally, observe that in Apprgach A, differen;LTFs (+ with the same in-

(i.i.d.) image. The interaction parameters (Y, /3, y, and 6 are given in Table 1, which also depicts the values of

teraction

parameters

image

air

PI,

p2,

p3,

The

p4?

Gibbs

P6r

yl,

sampler

247

IMAGES

FIG. 2. Nine realizations of a MC-GRF, which are obtained by approximating the GRFs depicted in Fig. I, by using Approach A (Experiment I).

an i.i.d.

GRF. 1).

FIELD

is

with

of a general (Experiment

RANDOM

y2,

and 6, but with different interaction parameters CY*, (Ye, (~4, Ps, and y4, result in the same LTF 7. In the case of Approach B, different LTFs (T with the same parameters (~1, fir, j32, p3, yi, y2,y3, and 6, but with different interaction parameters 02, (Ye,ff4, p4, pS, &, and y4, result in the same LTF Q.’ V.

SIMULATION

sty,

y3,

EXPERIMENTS

We shall now demonstrate the importance of our results with some simulation experiments (in all cases T = 1). EXPERIMENT 1. The purpose of this experiment is to show some “good” unilateral approximation of various GRF images. Figure 1 depicts nine 64 x 64 realizations of binary GRFs. These realizations are obtained by employing the Gibbs sampler with random site updating [9, 171. The Gibbs sampler is initialized by a purely random

’ This remark demonstrates the fact that, a MC-GRF is unable to capture 45” diagonal interactions, especially when the 45” diagonal interaction is due to a large value of the parameter ps(y, w), a result of a similar flavor as the one derived in [39].

z, w) =

c

du,

y, z, WI.

IlEE”

The GRF realizations depicted in Fig. 1 are not MC-GRF realizations, because their LTF (T does not satisfy Eq. (6).8 Figures 2 and 3 depict the MC-GRF realizations obtained by using Approaches A and B, respectively. Observe that no significant improvement is made by using Approach A over the simpler Approach B. Figure 4 depicts the GRF realizations obtained by the Gibbs sampler when it is initialized by the MC-GRF realizations depicted in Fig. 3. Observe the striking similarity between these results. Various statistics of the Monte Carlo simulation algorithms, which are used to obtain the realizations depicted in Figs. 1 and 4, are summarized in Table 2. Observe that, by using the MC-GRF realizations depicted in Fig. 3 as the initial realizations, we are able to obtain a speedup in calculations by an average factor of 5.76. Finally, we should emphasize the fact that the calculation of the MC-GRF realizations depicted in Figs. 2 and 4 took an average of 2.2 s of CPU time. In order to compare the quality of the previous approximations, the quantities 4(x, y, Z, o) and 40(x, y, z, o) are 8 An “close”

exemption to this to a MC-GRF.

remark

is GRF32,

which

is taken

to be

248

JOHN

GOUTSIAS

TABLE 1 Interaction Parameters of the General GRFs Depicted in Fig. 1 (Experiment 1) Parameters

s(O,O,O) s(l,O,O) SK414)

s(1,l.O) s(O.O.1) s(1,O.l) s(0,l.l)

s(1.1.1)

GRFII

GRFI2

-0.2600 1.0000 -0.5000 -0.5000 2.1000 -2.0000 0.1300 0.0000 0.0150 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

-2.0000 1.OOoo -0.5000 ~0.5000 I .oooo 0.7000 0.2222 0.3000 - I .6666 0.0000 -0.1666 0.6666 -0.2333 0.6666 0.3333

1.7711 19.8340 1.1391 13.4712 0.6698 3.0998 0.4116 2.0002

1.1353 3.7183 0.7090 3.091 I 9.7718 0.7608 0.4670 I .3337

GRF22

GRF23

GRF3

0.5000 -2.0000 I .oooo 0.0000 0.1200 0.0000 0.1000 0.4000 0.1000 -0.6200 0.9000

0.2000 0.0000 0.0000 0.0000 -4.0000 -5.0000 -0.3000 0.0000 0.2000 0.0000 0.6000 -0.2500 0.0100 0.3000 0.0010

-0.6265 -0.7350 0.7350 0.0000 0.0492 0.5927 0.0000 0.0000 -0.2286 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

-3.2000 -0.9180 0.0300 ~0.4000 3.4754 3.0000 -1.8200 0.3380 -0.7250 -0.3500 1.1225 - 1.4380 2.5000 -0.8300 -0.9000

2.0000 2.6487 3.7183 5.9530 I.1353 I .5028 I .4066 2.2553

2.2214 I .0224 0.9048 1 .0302 1.0082 I.221 1.0062 I .6490

I .5345 0.7487 3.2000 I.5614 I .9668 0.7689 4.1017 1.6036

I .0408 5.8026 1.0373 5.9848 I.2191 5.9189 1.2734 5.9366

GRF13

GRF2

0.0000 0.0000 0.0000 0.0000

FIG. 3. Nine realizations of a MC-GRF, which approximating the GRFs depicted in Fig. 1, by using _. periment 1).

I

are obtained by Approach B (Ex-

FIG. 4. Nine realizations initialized with the MC-GRF

1

GRF32

GRF33

0.8954 0.0000 0.0000 0.0000 - 1.7652 - 1.8428 0.3345 0.0000 - I .2500 0.0000 -0.1307 4.9244 -0.0678 -0.5996 0.4038

-2.0000 0.0000 0.0000 -2.9000 0.0000 5.0000 - I .oooo 0.0000 0.0000 0.9000 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 - I .2500 - 1.2500 0.5950 o.oooo 0.5950 o.oooo 0.0000 0.0000 0.0000 0.0000 0.0000

3.4483 1.4190 4.4209 1.5138 I .3877 2.9028 I .5063 2.6216

1.1353 1.1353 1.0498 I .0498 1.1602 1.1602 1.1353 I.1353

2.0000 I .2865 2.8130 1.5194 1.2865 I.9619 1.5194 2.0829

of a general GRF. The Gibbs sampler images depicted in Fig. 3 (Experiment

is I).

GIBBS RANDOM

Statistics of the Monte Carlo Simulation Statistics i.i.d. initial realization

MC-GRF initial realization

GRFll

Iterations Attempted replacements Replacements CPU time (s) Iterations Attempted replacements Replacements CPU time (s)

Absolute

TABLE 2 of the General GRFs Depicted

GRF13

GRF2l

GRF22

in Fig. 1 (Experiment GRF23

GRF3 1

1) GRF32

GRF33

13

20

5

14

10

22

16

13

7

49664 10502 328

79488 23585 523

19584 6613 130

54272 12134 359

39936 18585 264

87808 15688 576

64000 12627 418

50944 8004 334

26752 7566 176

2

5

2

3

2

2

3

5

2

6528 1385 45

I9968 5419 135

8064 2505 56

9856 2183 67

5248 2451 37

7424 873 51

10240 1662 70

I7024 1354 114

4352 1325 31

7.29

Speedup

Maximum

GRF12

249

FIELD IMAGES

3.87

2.32

5.36

7.14

11.30

5.97

2.93

5.68

TABLE 3 Relative Error between the Estimated Parameters q of the General GRFs Depicted in Fig. 1, and Their Corresponding Approximations with Approach B, Depicted in Fig. 3 (Experiment 1)

4

GRFl I

GRFl2

GRFl3

GRF21

GRF22

GRF23

GRF31

GRF32

GRF33

Maximum absolute relative error

0.0666

0.1099

0.1025

0.1014

0.0561

0.0374

0.1072

0.0074

0.0860

computed by a Monte Carlo simulation also Eq. (17a)). z, q(xy y’ z’ w) = lx:=, [lx,,,,

1401, where (see

l&(x, y, z, o)?Tk Uk(U, y, z, fiJJ 7rA (21a)

and (21b)

FIG. 5. Four realizations of a general GRF. The Gibbs sampler is initialized with an i.i.d. image (Experiment 2).

The quantities 4 and 90, given by Eq. (21), are related to the fourth-order statistics of the GRF H and the MC-GRF Ho, respectively; therefore, they can be effectively used to quantify the similarity between the realizations depitted in Figs. 1 and 3 [41]. The maximum absolute relative error in q, obtained by approximating the GRFs of Fig. 1 by the MC-GRFs of Fig. 3, is depicted in Table 3. In this case, we are able to obtain an average error of 0.0750, which is quite small. Finally, observe that the results depicted in Table 3 approximately verify the validity of condition (17a).

FIG. 6. Four realizations of a MC-GRF, which are obtained by approximating the GRFs depicted in Fig. 5, by using Approach A (Experiment 2).

250

JOHN GOUTSIAS

TABLE 4 Interaction Parameters of the General GRFs Depicted in Fig. 5 (Experiment 2) Parameters

GRFll

- 1.oooo 0.0000 0.0000 0.0000 0.0500 -0.0500 2.0000 0.0000 0.1000 0.0000 0.1000 -0.9000 0.0000 0.0000 -1.0000

s(O,O,O) s (1 AO) s(O, 1,O) s(l,l,O) s(O,O,1) s(1,0,1) s(O,1,1) .v(l,l,l)

1.3679 1.3867 3.7183 4.1582 1.3499 1.2705 3.5857 1.6018

GRF12

GRF21

2.5000 0.0000 0.0000 0.0000 -1.2500 -1.2500 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 ~2.5000 -2.5000

0.0000 0.0000 0.0000 0.0000

0.2000 0.0000 0.0000 0.0000 0.5000 -2.0000 ~5.0000 0.0000 5.0000 0.0000 0.1000 0.4000 0.1000 0.2000 0.9000

2.0000 1.0821 4.2871 1.2698 1.0821 3.3092 1.2698 3.3599

2.2214 3.0138 1.0082 1.0150 1.1653 208.7534 1.0012 182.7648

1.1900 0.0000 I.1900 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 13.1825 4.4903 13.1825 4.4903 4.4903 2.0000 4.4903 2.0000

GRF22

EXPERIMENT 2. We now present some simulation examples in which the proposed techniques fail to provide a “good” unilateral approximation to a general GRF. Figure 5 depicts four 64 x 64 realizations of binary GRFs, obtained by using the Gibbs sampler which is initialized by an i.i.d. image. The interaction parameters of these GRFs are given in Table 4, which also depicts the values of s(y, z, w). Observe again that these GRFs are not MC-GRFs. The GRFs of Fig. 5 are approximated by using Approaches A and B. The results of this approxima-

FIG. 7. Four realizations of a MC-GRF, which are obtained by approximating the GRFs depicted in Fig. 5. by using Approach B (Experiment 2).

tion are depicted in Figs. 6 and 7. Observe that, although some reasonable approximation is obtained for the cases of GRFll and GRF12, the proposed techniques fail to provide “good” approximations in the cases of GRF21 and GRF22. The reason for this is the fact that Eqs. (lOa), (lOe), and (lOi) are not satisfied (even approximately, as is the case in the previous experiment). For example, for GRF22 we have that ps = 5.0, whereas pi”’ = 0.6, which violates Eq. (10e). Figure 8 depicts the results of the Gibbs sampler when it is initialized by the MC-GRFs depicted in Fig. 7, whereas, Table 5 depicts statistics of the Monte Carlo simulation of the GRF realizations depicted in Figs. 5 and 8. In this case an average speedup of about 1.26 is obtained. Observe also that, in the case of GRF22, the simulation algorithm initialized by an i.i.d. image is 1.1 times faster than the simulation algorithm initialized by a MC-GRF. This is a direct consequence of the fact that condition (17a) is not satisfied in this case. The results depicted in Table 6 demonstrate this fact (i.e., we have an average error of 0.4390).

TABLE5

Statistics of the Monte Carlo Simulation of the General GRFs Depicted in Fig. 5 (Bxheriment 2) Statistics

GRFll

GRFl2

GRF21

i.i.d. initial realizations

Iterations Attempted replacements Replacements CPU time (s)

10 37632 9583 249

6 21248 7633 141

15 57600 6570 379

MC-GRF initial realization

Iterations Attempted replacements Replacements CPU time (s)

2649: 6306 176

5 17664 6278 118

37632 3219 249

Speedup

1.42

1.20

10

1.52

GRF22 13 52224 5204 343 I5 57600 4851 380 0.91

GIBBS

RANDOM

FIELD

251

IMAGES

TABLE 6 Maximum Absolute Relative Error between the Estimated Parameters q of the General GRFs Depicted in Fig. 5, and Their Corresponding Approximations with Approach B, Depicted in Fig. 7 (Experiment 2)

4 Maximum relative

absolute error

GRFll

GRFl2

GRFZI

GRF22

0.2928

0.5343

0.3345

0.5945

EXPERIMENT 3. A very important aspect of our work is studying the sensitivity of our approximations to the value of the interaction parameters. For example, in the case of Approximation B, different LTFs (+ with the same interaction parameters (Y,, p, , p2, pj, y, , y2. y3, and 6, but different interaction parameters (Ye, (Ye, cz4, p4, ps,

FIG. 8. Four realizations initialized with the MC-GRF

of a general GRF. The Gibt 1s sampler is images depicted in Fig. 7 (Ex .periment 2).

(a)

Ps(l,

1) = 0.0150

Ps(l,

1) = 0.1500

&(I,

1) = 1.ccoo

Ps(l,

I) = 1.5ooo

y&(1, 1) = 0.5950

Ps(l.

1) = l.m

Ps(l,

1) = 1.5Mxl

Ps(l.

1) = 2.cooo

obtained by an i.i.d. FIG. 9. Realizations of a general GRF with different values of parameter /?Z5(y, 0). The first row depicts the realizations initial image, whereas the second row depicts the realizations obtained by a MC-GRF initial image. The values of the interaction parameters (except ps) are the same as (a) the parameters of the GRFI 1, and (b) the parameters of the GRF33, in Experiment 1 (Experiment 3).

252

JOHN

Statistics of the Monte Carlo Simulation

GOUTSIAS

TABLE 7 of the General GRFs Depicted in Fig. 9 (Experiment

GRFI I p = 0.0150

Statistics i.i.d. initial realization

Iterations Attempted replacements Replacements CPU time (s)

MC-GRF initial realization

Iterations Attempted replacements Replacements CPU time ts)

Speedup

Maximum

Absolute

GRF33

pc = 0.1500

p< = 1.0000

II

IO

42880 9730 281

39808 10561 262

57984 10331 381

2

I

I3

2

7

13

I3

1263

27776 6766 184

51200 9090 337

4352 1325 31

25728 5862 170

51456 8084 339

50816 739s 335

13 49664

10502 327 2 6528 1385 45

5888 41

7.27

6.85

p5 = I.5000

pr = 0.5950 7

IS

1.42

p< = 1.0000

26752 7566 176

I.13

Maximum relative

absolute

p< = 1.5000

9 33192 7943 222

5.67

p5 = 2.0000

I4

15

54784 9502 362

1.30

60672 8866 400

1.06

I.19

TABLE 8 Relative Error between the Estimated Parameters q of the GeneralGRFs Depicted in Fig. 9, and Their Corresponding Approximations with Approach B, Depicted in Fig. 10 (Experiment 3) GRFl1

q

3)

GRF33

ps = 0.0150

ps = 0.1500

ps = 1.0000

ps = 1.5000

PC = 0.5950

ps = 1.0000

ps = 1.5000

ps = 2.0000

0.1251

0.1560

0.3529

0.8279

0.1047

0.2223

0.3377

0.3940

error

p6, and y4, result in the same LTF TV, Therefore, the quality of the resulting approximation will be a function of the interaction parameters cy2, (Ye,a4, p4, /35, &, 74. In this experiment we demonstrate the sensitivity of the approximation to the value of parameter ps. A similar sensitivity analysis can be performed easily in the case of other interaction parameters. Figure 9 depicts five realizations of a general GRF with different values of parameter ps. In Fig. 9a all the interaction parameters (except parameter ps) are the same as the parameters of GRFl 1, depicted in Fig. 1, whereas, in Fig. 9b all the interaction parameters (except parameter ps) are the same as the parameters of GRF33, depicted in Fig. 1. The first row in Fig. 9 depicts the realizations obtained by an i.i.d. initial image, whereas, the second row depicts the realizations obtained by a MC-GRF initial image. Figure 10 depicts the corresponding approximations obtained by Approach A (first column) and Approach B (second column). As the value of parameter pS increases, the quality of approximation deteriorates (see also Table 8). Tables 7 and 8 depict the statistics of the Monte Carlo simulation algorithms used to obtain the results of Fig. 9, and the maximum absolute relative error, respectively. Observe that the speedup obtained by initializing the Monte Carlo simulation algorithm by using a MC-GRF realization strongly depends on the value of parameter ps, and this

speedup decreases as the value of parameter creases. VI.

ps in-

SUMMARY

We have exploited the relationship between a gen era1 Gibbs random field image and a mutually compal :ible

FIG. 10. approximating

Realizations of four MC-GRFs, which are obtained by the general GRFs depicted in Fig. 9. First column: approximation with Approach A. Second column: anoroximation with .. Approach B (Experiment 3).

GIBBS

RANDOM

Gibbs random field image, we investigated the possibility of approximating a general Gibbs random field by a mutually compatible Gibbs random field, and we studied the consequences of such an approximation. We presented a general analysis, which resulted in important conclusions about the developed relationship. We considered two approaches to the approximation problem. We showed that in some important cases the second approach results in a particular case of a mutually compatible Gibbs random field whose probability measure minimizes the relative entropy. We demonstrated experimentally that, in many interesting cases, a general Gibbs random field image can be approximated by a mutually compatible Gibbs random field image. We also showed that substantial computa-

FIELD

253

IMAGES

tional savings can be obtained if we initialize the Monte Carlo simulation algorithm, for the simulation of Gibbs random field images, with the resulting approximation.

APPENDIX

On the Relationship

A

between GRF Parameters

In this appendix we shall summarize the relationship between the interaction parameters of a GRF, its LTF, and its CPFs. These relationships allow us to switch freely among various GRF representations. From Eq. (5) it is quite straightforward to prove that

Y = T In ~(0, 0, 0, O),

(Ala)

T a(~,0,0,0) a1(x)= x Inu(0,0,0,0)’ C.UZ(Y)

@lb)

T ~(0, y, 0, 0) = ; In 40, 0, 0, 0)’

64’~)

T cr(0, 0, z, 0) a3(z) = 7 In a(0, 0, 0, 0) ’

(A’4

T ~(0, 0, 0, 0) a4(w) = ii In u(0, 0, 0, 0) ’

(Ale)

:

In

u(x, y, 0, W(0, dx, 0, 0, W(O,

0, 0, 0) y, 0, 0)’

(A’f)

& In @(X7 0, 0, wb(O, 0, 0, 0) 4x9 0, 0, okdo, 0, 0, 0) ’

(AM

s In

u(x, 0, z, W(O, @(XT 0, 0, W(O,

0, 0, 0) 0, z, 0) ’

(Al’4

5 In

aa y, z, W(O, a@, Y, 0, W(O,

0, 0, 0) 0, z, 0)’

(Ali)

gin

do, 40,

&In

do, 0, z, wb(O, 0, 0, 0) -to, 0, z, Ob(O, 0, 0, WI’

y, 0, w)u(O, 0, 0, 0) Y, 0, okdo, 0, 0, w) ’

Wj) (A’k)

$1,

u(x, y, Z? O)u(x, 0, 0, Ob(O, y, 0, okm dx, y, 0, Obk 0, 2, Ob(O, y, z, Okm

0, z, 0) 0, 0, 0)’

&In

u(x, y, 0, wkdx, 0, 0, W(O, y, 0, Ob(O, 070, 0) 4x, y, 0, Ob(x, 0, 0, wb(O, y, 0, w)u(O, 0, 0, 0)’

(AM

&‘”

ah, 0, z, wb(x, 0, 0, W(O, 0, z, O)dO, 0, 0, WI u(x, 0, z, O)u(x, 0, 0, wb(O, 0, z, wb(0, 0, 0, 0)’

(AIn)

&In

do, Y, z, wb(O, y, 0, Ob(O, 0, z, %(O, 0, 0, 0) UK4 y, z, Ob(0, y. 0, wb(O, 0, 2, wb(0, 0, 0, 0)’

(Alo)

(All)

254

JOHN GOUTSIAS

and

(Alp) From Eqs. (l), (2), (4), and (5) we can also prove that T ’ = m [(YI(X)

+

az(x)

+

ZLTF In z,,, ’

+ ~4Wl-G

a3(x)

Wa) - 1

forl5i5M--l,l5j5N forlSiSM-

l,j=N Wb)

for i = M, 1 5 j 5 N - 1 for V!?‘(x 1.l



y)

[pl(x,

y>

+

PI@,

YbY,

/36(y,

i

=

M,

j

=

IV,

forl%i5M,lrj5N-1

x)]xy,

=

642~)

for 1 5 i 5 M, j = N, forI5i5M-l,IljSN

(A24

for i = M, 1 5 j 5 N, for 1 5 i 5 M, 1 I j % N, v!2s4'(y fJ vy(x,

' y,

0) z)

=

&(y,

=

y,(x,

o)yo, y,

z)xyz,

for 1 5 i 5 M, I 5 j 5 N,

(A2e) WW (A&) (A2h)

for 1 5 i 5 M, 1 % j 5 N, for 1 5 i 5 M, 1 I j 5 N,

v!3,3'(x ,J

'3 z

0)

=

y3(x

2 z 5 o)xzw,

for 1 5 i 5 M, 1 % j I N,

(A2i)

v!!,4'(y r,

3) z

w)

=

y4(y

9 z 2 0)yzo:

for 1 I i % M, 1 I j 5 N,

Wj)

and

vqx I, ' y ' z ' 0)

= 6(x 3 y 3 z 5 w)xyzw,

Parameter r, given by Eq. (A2a), cannot be computed in general since it is given in terms of the partition functions Ze,, and ZLrF. The particular value of this parameter does not affect the statistical behavior of the GRF H, since the probability measure Pr[H = h] does not depend on this parameter (see Eqs. (4) and (5)). Equation (A2) is a set of necessary and sufficient conditions for representation (4) to be equivalent to representation (1). Observe also that, even if the LTF u is spatially invariant, the GRF H has not spatially invariant CPFs I’,, in general. For this to be true it is required that (Ye = (w3(x)= (Ye

for 1 5 i 5 M, 1 I j I N.

= 0, for every x E EH and every (x, y) E EL.

p4(x,

APPENDIX

On the Relationship

WW

y)

=

&(x,

y) = 0, for

B

between the Interaction of a MC-GRF

Parameters

In this appendix we shall derive the relationship between the interaction parameters of a MC-GRF. This relationship, which is a direct result of the fact that the LTF

GIBBS

RANDOM

FIELD

T(X, y, z, W) of a MC-GRF satisfies Eq. (6b), is used in the derivation of Approach A. From Eqs. (6b) and (9), we easily obtain R3 equations of the form Y(O)+ cr:“‘(y)y + =

P:“‘(Y? -Tin

+ a~“‘(z)z + a~“)(“)w + @‘(y, W)YO

+ PIp’(z,

o)zw

+ yi$“)(y,

Q(Y, I, ~1 = ,z, exp I-) [a;“)(u)u + p\O’(u, y)uy + p$O’(u, w)uw + P:“‘(u, z>uz + r$“‘(u, y, z)uyz + y:“‘(u, y, w)uyw

z)yz z, w)yzw

Q(Y, z, w),

255

IMAGES

+ y:O’(u,

z, w)uzo

+ 6’“‘(u,

y,

z, w)uyzo]

I

.

@lb)

0314

for (y, z, W) E EL, where

From Eq. (Bl) we can easily show that

do) = T In CX(~)( y) = T In 2

Y

1

c ua,, A(u) [ %aH

[ %,aH

L2

ueEH

w) = j$ In 2

w) = 2 In

z) I ’

~EEH

ae,,

trcs

~1 I ’

Y>

ACuP

I (~1,

A(u)

&EH

A(u)BI(u,

lx

A(u)BI(u,

r,e,,

A(LO

&EH

%a,

c

U-Q4

.4(uW3(~,

[ &ran

pr’(z,

Wb)



A(u)

c

p$“(y, z) = $ In

p$“(y,

1

032~)

A(u)B,(u,

&,,

CY(~)(OJ) = T In 4 0

Y)

A(u)B,(u,

&,EEH A (~1

a(O)(z) = _TIn 3 Z

WW

[ 2 ~EEHA(u) I ’

%aH

z)

zN’~(u, ~1

w)rz(u,

z) %a~,, A(u)B3(u, A@)B2(u,

z)&(u,

We)

Y, z)

A(uW3@4,

YW~(U,

A(ulB~(u,

&E,,

A(u)B~(u,

Y)Bz(u, Y)

AbW2(u,

A(u)

%I,,

Y, ~1 01

w)r&,

z,

o)

1 1

Wf)



W&9



and

Y$(Y,

2 ~EEH A(u) &a z,

w>

=

&

A(WI(U,

~)B2@4,

Y,

z)~I(u,

z)

ln lx urz~H

Y)

A(u)BI(u,

I2lrEE,,

X ~,,,~w-h(u,

~,,~E,,

z)B~(u,

Aw2k

YP32b,Zuw4,

AbW2b4,

axu,

y,

z>

w)r2(U, dr2(u,Y,

Z, 0) oP73(u,z,

w)A(u,

y,

z, w)

1 WW ’

where

A(x) = exp 1~cx~~)(x)x I ,

GW

BI(x,

Y) = ew

B2k

z) = exp {+ p:o’(x, zbz],

(tp;“‘k

YbY],

(B3b)

033~)

JOHN GOUTSIAS

256

APPENDIX

B3(x, 0) = exp

C

(BW

Proof of Proposition We)

(B3f)

In this appendix we shall show that a GRF with homogeneous CPFs can never become equivalent to a nontrivial MC-GRF. If the CPFs are homogeneous then (see also Appendix A) a?(x) = a&c) = a4(x) = 0,

(B3g)

for x E EH,

and and p4k

Ah

y, z, w) = exp + 8(“Yx, y, z, ohyzw

W-d

Y) =

p6t-h

y) = 0,

for (X, Y) ‘2 EL,

from which we obtain that (see also Eqs. (B2b)-(B2e) and (B2g))

for every y E EH,

,,& A(u) = ,,& A(~)BI(~, YL

for every z E EH,

for every w E EH,

&

A(u) = &

A(~)BI(u,

YWZ(U,

z)r3(u, Y, z>, for every (y, z) E EL,

c

A(U) = c

A(u)&(M,

z)B3(u,

w)r,(u,

and

If we take B, (u, L) = minyEE,, {BI(u, y)>, then from Eq. (Cla) we have that ,,z,, A(~)BI(u,

Y) =

,,z,, A(u)BI(u, Yh

or &

A(u)tB~(u,

Y)

Z,

w),

for every (z, W) E EL.

(Cle)

sinceA(u)zOandBl(u,y)rB,(u,y),foreveryyEEH, from which we conclude that p,(x, y) is not a function of y. From Eqs. (Cl), and by following a similar argument, we can also show that the parameters p2(x, w), &(x, z), y 1(x, y, z), and y3(x, z, o) do not depend on ( y, z, w), which, together with Eqs. (B2) and (B3), let us conclude that the resulting random field is the trivial (i.i.d.) MCGRF [20].

- Bl(U, 31 = 0, ACKNOWLEDGMENT

or Bl(U, Y) = B,(u, L),

for every y E EH,

This work was supported by the Office of Naval Research, Mathematical Sciences Division, under Grant N00014-90-J-1345.

GIBBS RANDOM

FIELD IMAGES

257

tion and identification, in Proceedings,

REFERENCES

ence on Acoustics,

1. .I. Besag, Spatial interaction and the statistical analvsis of lattice systems, J.-R. Statist. Sot. Ser. B 36, 1974, 192-238. 2. R. C. Dubes and A. K. Jain, Random field models in image analysis, J. Appl. Statist. 16, 1989, 131-164. 3. F. S. Cohen, Markov random fields for image modeling and analysis, in Modeling and Application of Stochastic Processes (U. B. Desai, Ed.), pp. 243-272, Kluwer Academic Publishers, Boston, MA, 1986. 4. F. R. Hansen and H. Elliott, Image segmentation using simple Markov field models, Compur. Graphics Image Process. 20, 1982, 101-132. 5. H. Derin, H. Elliott, R. Cristi, and D. Geman, Bayes smoothing algorithms for segmentation of binary images modeled by Markov random fields, IEEE Trans. Pattern Anal. Mach. Intelligence 6, 1984, 707-720. 6. H. Derin and W. S. Cole, Segmentation of textured images using Gibbs random fields, Comput. Vision Graphics Image Process. 35, 1986, 72-98. 7. H. Derin and H. Elliott, Modeling and segmentation of noisy and textured images using Gibbs random fields, IEEE Trans. Pattern Anal. Mach. Intelligence 9, 1987, 39-55. 8. P. A. Kelly, H. Derin, and K. D. Hartt, Adapative segmentation of speckled images using a hierarchical random field model, IEEE Trans. Acoust. Speech Signal Process. 36, 1988, 1628-1641. 9. S. Geman and D. Geman, Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images, IEEE Trans. Pattern Anal. Mach. Intelligence 6, 1984, 721-741. 10. J. Besag, On the statistical analysis of dirty pictures, J. R. Statist. Sot. Ser. B 48, 1986, 259-302. 11. M. Hassner and J. Sklansky, The use of Markov random fields as models of texture, Comput. Graphics Image Process. 12, 1980, 357-370.

12. G. R. Cross and A. K. Jain, Markov random field texture models, IEEE Trans. Pattern Anal. Mach. Intelligence 5, 1983, 25-39. 13. J. Konrad and E. Dubois, Estimation of image motion fields: Bayesian formulation and stochastic solution, in Proceedings, IEEE International cessing, April

Conference on Acoustics, 11-14, New York, NY,

Speech,

and

Signal

Pro-

pp. 1072-1075, 1988. 14. P. Bouthemy and P. Lalande, Motion detection in an image sequence using Gibbs distributions, in Proceedings, IEEE International Conference on Acoustics, May 23-26, Glasgow, Scotland,

Speech,

and

Signal

Processing,

pp. 1651-1654, 1989. 15. J. Marroquin, S. Mitter, and T. Poggio, Probabilistic solution of illposed problems in computational vision, J. Am. S/atist. Assoc. 82, 1987, 76-89. 16. J. W. Modestino and J. Zhang, A Markov random field modelbased approach to image interpretation, in Proceedings, IEEE Computer Vision San Diego, CA,

and Pattern

17. J. Goutsias, A Theoretical the Simulation

Recognition

Conference,

June

4-8,

pp. 458-465. 1989. of Gibbs

Analysis of Monte Carlo Algorithm.sfbr Random Field Images, Technical Report

JHU/ECE 89-07, Department of Electrical and Computer Engineering, The Johns Hopkins University, Baltimore, MD, 1989. 18. D. K. Pickard, Inference for general Ising models, J. Appl. Probab. A 19, 1982, 345-357.

19. J. Goutsias, Mutually compatible Gibbs images: Properties, simula-

Speech,

and

Signal

IEEE International ConferProcessing, April 11-14,

New York, NY, pp. 1020-1023, 1988. 20. J. Goutsias, Mutually compatible Gibbs random fields, IEEE Trans. Inform. Theory 35, 1989, 1233-1249. 21. K. Abend, T. J. Harley, and L. N. Kanal, Classification of binary random patterns, IEEE Trans. Inform. Theory 11, 1965, 538-544. 22. L. N. Kanal, Markov mesh models, Comput. Graphics Image Process.

12, 1980, 371-375.

23. D. K. Pickard, Unilateral Ising models, Suppl. Adu. Appl.

Probab.

10, 1978, 58-64.

24. D. K. Pickard, Unilateral Markov fields, Adu. Appl.

Probab.

12,

1980, 655-671.

25. T. R. Welberry and R. Galbraith, A two-dimensional model of crystal-growth disorder, J. Appl. Crystallogr. 6, 1973, 87-96. 26. I. G. Enting, Crystal Growth models and Ising models: Disorder points, J. Phys. C Solid Sfate Phys. 10, 1977, 1379-1388. 27. M. Kanefsky and M. G. Strintzis, A decision theory approach to picture smoothing, IEEE Trans. Comput. 27, 1978, 32-36. 28. J. Haslett, Maximum likelihood discriminant analysis on the plane using a Markovian model of spatial context, Pattern Recognit. 18, 1985, 287-296.

29. P. A. Devijver, Probabilistic labeling in a hidden second order Markov mesh, in Pattern Recognition in Practice 11 (E. S. Gelsema and L. N. Kanal, Eds.), pp. 113-123. Elsevier, Amsterdam, 1986. 30. V. Lacroix, Pixel labeling in a second-order Markov mesh, Signal Process.

12, 1987, 59-82.

31. M. Kanefsky and C-B. Fong, Predictive source coding techniques using maximum likelihood prediction for compression of digitized images, IEEE Trans. Inform. Theory 30, 1984, 722-727. 32. J. Goutsias, Unilateral approximation of Gibbs random field images, in Proceedings, IEEE International Conference on Acoustics, Speech, and Signal Processing, April 3-6, Albuquerque, NM, pp. 2049-2052, 1990. 33. H. W. J. Blote, A. Compagner, and A. Hoogland, The simple quadratic Ising model with crossing bonds, Physical A 141, 1987, 375402.

34. I. G. Enting and T. R. Welberry, Connections between Ising models and various probability distributions, Suppl. Adu. Appl. Probab.

10, 1978, 65-72.

35. R. Kindermann and J. L. Snell, Markou Random Fields and Their Applications, Contemporary Mathematics, Vol. 1, American Mathematical Society, Providence, RI, 1980. 36. B. D. Ripley, Stochastic Simulation, Wiley, New York, 1987. 37. G. E. Forsythe, M. A. Malcolm, and C. B. Moler, Computer Methods for Mathematical Computntions, Prentice-Hall, Englewood Cliffs, NJ, 1977. 38. R. E. Blahut, Principles and Pructice of Information Theory, Addison-Wesley, Reading, MA, 1987. 39. F. C. Jeng and J. W. Woods, On the relationship of the Markov mesh to the NSHP Markov chain, Patrern Recognit. Left. 5, 1987, 273-279. 40. M. H. Kalos and P. A. Whitlock, Monte Car/o Methods, Vol. I, Basics, Wiley, New York, 1986. 41. A. Gagalowicz, A new method for texture fields synthesis: Some application to the study of human vision, IEEE Trans. Pattern Anal.

Mach.

Intelligence

3, 1981, 520-533.