- Email: [email protected]
Microcanonical optimization applied to visual processing
_l!!il3
25 September 1995 PHYSICS
me
ELSEVTER
LETTERS
A
Physics Letters A 205 (1995) 377-382
Microcanonical optimization applied to visual processing JosC R.A. Torreiio l, Edward Roe 2 Departamenro de Informdtica. Universidade Federal de Pematnbuco, 50.732-970 Recife, PE, Brazil
Received 16 January 1995; revised manuscript received 26 January 1995; accepted for publication 2 August 1995 Communicated by A.&?Fordy
Abstract We introduce an optimization algorithm, based on Creutz’s microcanonical simulation technique, which has proven very efficient for non-convex optimization tasks associated with image-processing applications. Our algorithm should also constitute a useful heuristic for applications in other domains requiring combinatorial optimization searches. Keywords: Creutz’s algorithm; Image processing; king systems; Microcanonical simulations; Simulated annealing; Stochastic optimization
Since the 1980’s, statistical physics techniques have found their way into the fields of image processing and image analysis [ l-41. Such techniques are usually based on, or suitable modifications of, the simulated annealing approach, which is a stochastic heuristic appropriate for global optimization problems where a non-convex functional is to be minimized [ 51. The simulated annealing procedure is most commonly run on top of a Metropolis-type Monte Carlo process [ 61, and is therefore based on the canonical ensemble a parallel being drawn between the visual processing problem at hand and the evolution of a physical system in equilibrium at a given temperature. Lately, there have also appeared simulated annealing applications based on the microcanonical ensemble, where the system’s evolution is controlled not by its temperature, but by its internal energy [ 1,4]. Microcanonical annealing (MA) has some definite advantages over the canonical simulated annealing (SA) , since it does not require the generation of quality random numbers or ’ E-mail: [email protected]. 2 E-mail: [email protected]. Bkevier Science B.V. SSDI 0375-9601(95)00585-4
the evaluation of transcendental functions, thus allowing much faster implementations [ I]. Here we present some novel contributions to the use of microcanonical simulation techniques for image-related applications. We have developed a microcanonical optimization algorithm (the ~0 algorithm) which retains the positive features of MA, but which does not involve annealing, thus yielding further gains in processing times. ~0 has been applied to the problems of binary-image restoration and halftoning - which correspond to estimating the ground states of a ferromagnetic and an antiferromagnetic Ising system, respectively [2] - and it has been found to outperform the annealing approaches. Our algorithm consists of two procedures which are alternately applied, starting from an arbitrary initial configuration, until the system has stabilized in a lowenergy state: an initialization procedure and a sampling procedure. In the initialization phase, the system is placed in a state with a lower energy than its starting configuration through an iterative improvement process by which local changes in the system dynamical variables are
318
J. R.A. Torreiio, E. Roe/Physics
proposed, and accepted only when leading to smaller energies. For optimization problems involving binary variables, such as the ones discussed below, this procedure is equivalent to the method of the iterated conditional modes (ICM) introduced in Ref. [ 71, where the local configuration chosen is always that which minimizes the energy function. The initialization procedure in ~0 is made to stop when a local minimum in the energy landscape has been reached. From the state obtained in the initialization, the sampling phase proceeds to generate samples of the system’s configurations at that energy, through the Creutz microcanonical simulation [ 81. In Creutz’s technique, an extra degree of freedom, called a demon, travels around the system carrying a variable amount of (always positive) energy, ED, in its bag. The purpose of the demon is to exchange energy with the system, in such a way that the total energy, Etotat = Es + ED, is preserved, where ES is the system’s energy. In the Creutz procedure, the demon attempts a local change in the system’s configuration; all such changes which lower the system’s energy are accepted, and the difference -AEs > 0 is added to the demon’s bag. On the other hand, those transitions which bring about an increase in the system’s energy are accepted only if the energy difference AEs > 0 satisfies AEs < ED, so that it can be provided by the demon (recall that ED is constrained to be always positive). In this way, the total energy ES + ED is guaranteed to remain constant. For the implementation of MA, as suggested in Ref. [ 11, an annealing schedule is followed for the gradual decrease of the demon’s energy. In the ~0 algorithm, we impose an upper bound on the energy which the demon can hold. Thus, assuming that the demon’s initial energy is El, and that it carries a bag of capacity EF, the sampling phase in our algorithm corresponds to a microcanonical simulation of the system in the energy range ES - EF + EI < E < ES + EI, where ES, here, is the system’s final energy obtained in the initialization procedure. The sampling phase is made to stop when thermodynamical equilibrium has been reached, which can be rigorously established from an analysis of the distribution of the demon’s states [ 1,4], although a simpler heuristic has been followed in our experiments (see below). The final state resulting from the sampling phase will then be used as the initial configuration for another initialization run, and thus consecutively, un-
Letters A 205 (1995) 377-382
til no further reduction in the system’s energy can be obtained. When compared to the annealing-based approaches, we see that the ,uO algorithm, besides preserving the advantages of the microcanical simulations, brings the extra benefit of not requiring a time-consuming annealing schedule, with the associated trial-and-error process of determination of free parameters (initial value and update rate for the temperature/demon energy, and number of sweeps at each temperature/demon energy). In our algorithm, the prescriptions for halting the initialization and sampling phases are unequivocal, and essentially only the values for EI and EF, respectively, the demon’s initial energy and capacity, have to be chosen a priori. In the experiments discussed below, we chose EI and EF, at each sampling phase, as fractions of the energy obtained in the previous initialization run ( EI = Es/200 and EF = Es/lOO), but largely equivalent results have been obtained when those parameters were kept fixed throughout the whole implementation. Here we show examples of the application of ~0 to some image-based processes, and compare its performance to those of the canonical and the microcanonical annealing algorithms. The applications which we discuss have been formulated as optimization problems and solved via SA by Camevali et al. in Ref. [ 21. They are image-processing applications (binaryimage smoothing, and halftoning of a continuous-tone image), which have been demonstrated to be formally equivalent to ground state problems for 2-D Ising systems. We have arrived at solutions to both problems via SA and MA, and through the ,xO algorithm, which has yielded faster and usually higher-quality results. The first application (experiment 1) is to obtain a smoothed version of a given binary image which has been corrupted by random, also binary, noise. The solution is obtained as the image p = {&j} which minimizes the energy function (cost function)
E(P)=
C(Pij-~kt)2+hC(cuij--Pij)2, (1) (ij,kl)
ij
where (Y= {aij} is the degraded image (aij and pij can assume the values 0 or 1) , and the notation (ij, kl) indicates that the summation is over nearest-neighbor pixels, only. The first term in ( 1) encodes a smoothness condition for the reconstructed image, /?, while
1. R.A. Torreiio, E. Roe/Physics
Letters A 205 (I 995) 377-382
la
lb
1C
Id
379
Fig. 1. Binary image smoothing experiment: (a) original image; (b) degraded image; (c) restoration via microcanonical optimization (~0) ; (d) restoration via microcanonical annealing (MA).
the second term enforces its fidelity to the observed image, CL The parameter A controls the balance between the two terms, and has been chosen as A = 1.2. Fig. la shows the original image from which the degraded version CY,in Fig. lb, was obtained. Fig. lc shows the smooth image, /3, resulting from the minimization of ( 1) via the ,uO algorithm. For comparison, in Fig. Id we show the result obtained through the application of microcanonical annealing. As a second application (experiment 2), we consider the problem of halftoning, which is to construct a binary image whose average intensities emulate a given continuous-tone image. Calling ,u = {pLij} the binary image (intensities + 1 or - 1) and y = {rij} the
gray-level image (intensities between +l and - 1) , and defining the average intensities of p through the introduction of a filter (blurring function) V, as Pij =
c
Kj,klPkl,
(2)
kl
we can formulate the problem of halftoning as that of minimizing the cost function
ij
The images p obtained from the minimization of (3) via the ~0 algorithm, and via MA, are presented in Figs. 2b and 2c, respectively, with the input image
380
J.R.A. TorreZo. E. Roe/ Physics Letters A 205 (1995) 377-382
(a)
(b)
Fig. 2. Halftoning experiment: (a) continuous-tone image; (b) halftone image obtained via microcanonical optimization (,uO) ; (c) halftone image obtained via microcanonical annealing (MA).
1. R.A. Torreiio. E. Roe/Physics
381
Letters A 205 11995) 377-382
ENERGY 60000
"MICROCC\NONICAL_OPTIMIZATION" ___. "~IICROCIWNICAL-ANNEfiLING" "CFNONICAL_BNNEALING" "KM"
Y””
.. .. '. '. '.
'. '_ '. ._ '.., '. '.
'. '..
*..,,.....,,,, '...., c-.__________.---me-- "'."".........,,.,,,,,,,,,
1
10000
I
0
5
20
15
10
I
I
I
25
30
35
I
40
STEPS ENERGY I
I
I
I
i (b)
"MICROC~INONIC~L_OPTI~~IZ~TION" -
I
120
I I r:... .._ *'. .-._,, I .-'. . .-.. ._._ '. '... .
._..
"MICROCANONICAL_ANNEALING'--"CANONICAL_ANNEALING".... "ICM" t.mw ”
.\ .. i.
6(
1
-. ‘. ..
60 STEPS Fig.3. Convergence of microcanonical optimization (PO), as compared to simulated annealing (SA), microcanonical annealing (MA), and ICM optimization: (a) binary image smoothing experiment; (b) halftoning experiment.
J.R.A. Tome&
382
E. Roe/Physics
Table 1 Performance of microcanonical optimization (~0) as compared to simulated annealing (SA) and microcanonical annealing (MA), for image smoothing (experiment I ) and halftoning (experiment 2) experiments Experiment 1 final energy ~0 13255.91 MA 13446.50 SA 13558.70
Experiment 2 CPU time (s) final energy 73.80 455.88 1917.06
1.74 2.04 2.07
CPU time (s) 217.16 434.10 729.26
y being shown in Fig. 2a. We have assumed V as a uniform 3 x 3 filter, i.e., xj,k[ = Yi__kl,(j_[l = i, for 0 < Ii - kl, Ij - II < 1. Table 1 summarizes the performance of the ~0, SA and MA algorithms for the experiments discussed above. The CPU times were obtained with a Sun Spare-II station. The algorithms were run sequentially, with SA (respectively, MA) following a geometric schedule, T c 0.9T (respectively, ED 6 0.9E~), for the decrease of temperature (respectively, demon energy), and both MA and ~0 being implemented with a single demon per lattice, in contrast with the multiple-demon implementations of Refs. [ 1,4]. The initial values of the temperature and the demon energy were 100, for the image smoothing experiment, and 10 for the halftoning experiment. When scanning the image lattices for updating the system’s states, a maximum of 100 sweeps (raster-scan visits of the whole image) were performed at each temperature (respectively, demon energy) for the SA (respectively, MA) algorithm, and at each sampling stage of ~0. When less than 10% of the tested state transitions were accepted, the system was considered to be in equilibrium, and the temperature (or demon energy) was updated according to the annealing schedules in SA (or MA) ; for the ~0, a new initiaIization procedure was run. The algorithms stopped when a fraction
Letters A 205 f 19951377-382
of less than 1% of the proposed transitions was being accepted. Fig. 3 illustrates the convergence of ~0, as compared to SA and MA, for the two applications discussed. In the graphs, a step is counted after each change in temperature, for SA, or demon energy, for MA, and after each initialization stage for ~0. For comparison, we also show the curves corresponding to the ICM optimization, which is equivalent to performing just the initialization procedure of our algorithm. In summary, the microcanonical algorithm (,uO) has been introduced as an alternative to the annealing approaches to non-convex optimization problems. Such an algorithm follows the ideas introduced by Creutz [ 81 and first exploited by Barnard [ l] in the context of stochastic optimization. We have made an assessment of the performance of ~0 when employed in some simple image-processing applications, and we have shown its higher computational efficiency as compared to the annealing approaches. We would like to thank Professor F.G. Brady Moreira for introducing us to Creutz’s algorithm.
References [ 11 ST. Barnard, Int. J. Comput. Vis. 3 (1989) 17. [2] P. Carnevali, L. Coletti and S. Patamello, IBM J. Res. Dev. 29 (1985) 569. [ 3 J S. Geman and D. Geman, IEEE Trans. Pattern Anal. Mach. Intell. 6 (1984) 721. [4] L. Herault and R. Horaud, IEEE Trans. Pattern Anal. Mach. Intell. 15 (1993) 899. [5] S. Kirkpatrick, C.D. Gelatt and M.P. Vecchi, Science 220 (1983) 671. [6] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller and E. Teller, J. Chem. Phys. 21 (1953) 1087. [7] J. Besag, J. R. Statist. Sot. B 48(3) (1986) 259. 181 M. Creutz, Phys. Rev. Lett. 50 (1983) 1411.
25 September 1995 PHYSICS
me
ELSEVTER
LETTERS
A
Physics Letters A 205 (1995) 377-382
Microcanonical optimization applied to visual processing JosC R.A. Torreiio l, Edward Roe 2 Departamenro de Informdtica. Universidade Federal de Pematnbuco, 50.732-970 Recife, PE, Brazil
Received 16 January 1995; revised manuscript received 26 January 1995; accepted for publication 2 August 1995 Communicated by A.&?Fordy
Abstract We introduce an optimization algorithm, based on Creutz’s microcanonical simulation technique, which has proven very efficient for non-convex optimization tasks associated with image-processing applications. Our algorithm should also constitute a useful heuristic for applications in other domains requiring combinatorial optimization searches. Keywords: Creutz’s algorithm; Image processing; king systems; Microcanonical simulations; Simulated annealing; Stochastic optimization
Since the 1980’s, statistical physics techniques have found their way into the fields of image processing and image analysis [ l-41. Such techniques are usually based on, or suitable modifications of, the simulated annealing approach, which is a stochastic heuristic appropriate for global optimization problems where a non-convex functional is to be minimized [ 51. The simulated annealing procedure is most commonly run on top of a Metropolis-type Monte Carlo process [ 61, and is therefore based on the canonical ensemble a parallel being drawn between the visual processing problem at hand and the evolution of a physical system in equilibrium at a given temperature. Lately, there have also appeared simulated annealing applications based on the microcanonical ensemble, where the system’s evolution is controlled not by its temperature, but by its internal energy [ 1,4]. Microcanonical annealing (MA) has some definite advantages over the canonical simulated annealing (SA) , since it does not require the generation of quality random numbers or ’ E-mail: [email protected]. 2 E-mail: [email protected]. Bkevier Science B.V. SSDI 0375-9601(95)00585-4
the evaluation of transcendental functions, thus allowing much faster implementations [ I]. Here we present some novel contributions to the use of microcanonical simulation techniques for image-related applications. We have developed a microcanonical optimization algorithm (the ~0 algorithm) which retains the positive features of MA, but which does not involve annealing, thus yielding further gains in processing times. ~0 has been applied to the problems of binary-image restoration and halftoning - which correspond to estimating the ground states of a ferromagnetic and an antiferromagnetic Ising system, respectively [2] - and it has been found to outperform the annealing approaches. Our algorithm consists of two procedures which are alternately applied, starting from an arbitrary initial configuration, until the system has stabilized in a lowenergy state: an initialization procedure and a sampling procedure. In the initialization phase, the system is placed in a state with a lower energy than its starting configuration through an iterative improvement process by which local changes in the system dynamical variables are
318
J. R.A. Torreiio, E. Roe/Physics
proposed, and accepted only when leading to smaller energies. For optimization problems involving binary variables, such as the ones discussed below, this procedure is equivalent to the method of the iterated conditional modes (ICM) introduced in Ref. [ 71, where the local configuration chosen is always that which minimizes the energy function. The initialization procedure in ~0 is made to stop when a local minimum in the energy landscape has been reached. From the state obtained in the initialization, the sampling phase proceeds to generate samples of the system’s configurations at that energy, through the Creutz microcanonical simulation [ 81. In Creutz’s technique, an extra degree of freedom, called a demon, travels around the system carrying a variable amount of (always positive) energy, ED, in its bag. The purpose of the demon is to exchange energy with the system, in such a way that the total energy, Etotat = Es + ED, is preserved, where ES is the system’s energy. In the Creutz procedure, the demon attempts a local change in the system’s configuration; all such changes which lower the system’s energy are accepted, and the difference -AEs > 0 is added to the demon’s bag. On the other hand, those transitions which bring about an increase in the system’s energy are accepted only if the energy difference AEs > 0 satisfies AEs < ED, so that it can be provided by the demon (recall that ED is constrained to be always positive). In this way, the total energy ES + ED is guaranteed to remain constant. For the implementation of MA, as suggested in Ref. [ 11, an annealing schedule is followed for the gradual decrease of the demon’s energy. In the ~0 algorithm, we impose an upper bound on the energy which the demon can hold. Thus, assuming that the demon’s initial energy is El, and that it carries a bag of capacity EF, the sampling phase in our algorithm corresponds to a microcanonical simulation of the system in the energy range ES - EF + EI < E < ES + EI, where ES, here, is the system’s final energy obtained in the initialization procedure. The sampling phase is made to stop when thermodynamical equilibrium has been reached, which can be rigorously established from an analysis of the distribution of the demon’s states [ 1,4], although a simpler heuristic has been followed in our experiments (see below). The final state resulting from the sampling phase will then be used as the initial configuration for another initialization run, and thus consecutively, un-
Letters A 205 (1995) 377-382
til no further reduction in the system’s energy can be obtained. When compared to the annealing-based approaches, we see that the ,uO algorithm, besides preserving the advantages of the microcanical simulations, brings the extra benefit of not requiring a time-consuming annealing schedule, with the associated trial-and-error process of determination of free parameters (initial value and update rate for the temperature/demon energy, and number of sweeps at each temperature/demon energy). In our algorithm, the prescriptions for halting the initialization and sampling phases are unequivocal, and essentially only the values for EI and EF, respectively, the demon’s initial energy and capacity, have to be chosen a priori. In the experiments discussed below, we chose EI and EF, at each sampling phase, as fractions of the energy obtained in the previous initialization run ( EI = Es/200 and EF = Es/lOO), but largely equivalent results have been obtained when those parameters were kept fixed throughout the whole implementation. Here we show examples of the application of ~0 to some image-based processes, and compare its performance to those of the canonical and the microcanonical annealing algorithms. The applications which we discuss have been formulated as optimization problems and solved via SA by Camevali et al. in Ref. [ 21. They are image-processing applications (binaryimage smoothing, and halftoning of a continuous-tone image), which have been demonstrated to be formally equivalent to ground state problems for 2-D Ising systems. We have arrived at solutions to both problems via SA and MA, and through the ,xO algorithm, which has yielded faster and usually higher-quality results. The first application (experiment 1) is to obtain a smoothed version of a given binary image which has been corrupted by random, also binary, noise. The solution is obtained as the image p = {&j} which minimizes the energy function (cost function)
E(P)=
C(Pij-~kt)2+hC(cuij--Pij)2, (1) (ij,kl)
ij
where (Y= {aij} is the degraded image (aij and pij can assume the values 0 or 1) , and the notation (ij, kl) indicates that the summation is over nearest-neighbor pixels, only. The first term in ( 1) encodes a smoothness condition for the reconstructed image, /?, while
1. R.A. Torreiio, E. Roe/Physics
Letters A 205 (I 995) 377-382
la
lb
1C
Id
379
Fig. 1. Binary image smoothing experiment: (a) original image; (b) degraded image; (c) restoration via microcanonical optimization (~0) ; (d) restoration via microcanonical annealing (MA).
the second term enforces its fidelity to the observed image, CL The parameter A controls the balance between the two terms, and has been chosen as A = 1.2. Fig. la shows the original image from which the degraded version CY,in Fig. lb, was obtained. Fig. lc shows the smooth image, /3, resulting from the minimization of ( 1) via the ,uO algorithm. For comparison, in Fig. Id we show the result obtained through the application of microcanonical annealing. As a second application (experiment 2), we consider the problem of halftoning, which is to construct a binary image whose average intensities emulate a given continuous-tone image. Calling ,u = {pLij} the binary image (intensities + 1 or - 1) and y = {rij} the
gray-level image (intensities between +l and - 1) , and defining the average intensities of p through the introduction of a filter (blurring function) V, as Pij =
c
Kj,klPkl,
(2)
kl
we can formulate the problem of halftoning as that of minimizing the cost function
ij
The images p obtained from the minimization of (3) via the ~0 algorithm, and via MA, are presented in Figs. 2b and 2c, respectively, with the input image
380
J.R.A. TorreZo. E. Roe/ Physics Letters A 205 (1995) 377-382
(a)
(b)
Fig. 2. Halftoning experiment: (a) continuous-tone image; (b) halftone image obtained via microcanonical optimization (,uO) ; (c) halftone image obtained via microcanonical annealing (MA).
1. R.A. Torreiio. E. Roe/Physics
381
Letters A 205 11995) 377-382
ENERGY 60000
"MICROCC\NONICAL_OPTIMIZATION" ___. "~IICROCIWNICAL-ANNEfiLING" "CFNONICAL_BNNEALING" "KM"
Y””
.. .. '. '. '.
'. '_ '. ._ '.., '. '.
'. '..
*..,,.....,,,, '...., c-.__________.---me-- "'."".........,,.,,,,,,,,,
1
10000
I
0
5
20
15
10
I
I
I
25
30
35
I
40
STEPS ENERGY I
I
I
I
i (b)
"MICROC~INONIC~L_OPTI~~IZ~TION" -
I
120
I I r:... .._ *'. .-._,, I .-'. . .-.. ._._ '. '... .
._..
"MICROCANONICAL_ANNEALING'--"CANONICAL_ANNEALING".... "ICM" t.mw ”
.\ .. i.
6(
1
-. ‘. ..
60 STEPS Fig.3. Convergence of microcanonical optimization (PO), as compared to simulated annealing (SA), microcanonical annealing (MA), and ICM optimization: (a) binary image smoothing experiment; (b) halftoning experiment.
J.R.A. Tome&
382
E. Roe/Physics
Table 1 Performance of microcanonical optimization (~0) as compared to simulated annealing (SA) and microcanonical annealing (MA), for image smoothing (experiment I ) and halftoning (experiment 2) experiments Experiment 1 final energy ~0 13255.91 MA 13446.50 SA 13558.70
Experiment 2 CPU time (s) final energy 73.80 455.88 1917.06
1.74 2.04 2.07
CPU time (s) 217.16 434.10 729.26
y being shown in Fig. 2a. We have assumed V as a uniform 3 x 3 filter, i.e., xj,k[ = Yi__kl,(j_[l = i, for 0 < Ii - kl, Ij - II < 1. Table 1 summarizes the performance of the ~0, SA and MA algorithms for the experiments discussed above. The CPU times were obtained with a Sun Spare-II station. The algorithms were run sequentially, with SA (respectively, MA) following a geometric schedule, T c 0.9T (respectively, ED 6 0.9E~), for the decrease of temperature (respectively, demon energy), and both MA and ~0 being implemented with a single demon per lattice, in contrast with the multiple-demon implementations of Refs. [ 1,4]. The initial values of the temperature and the demon energy were 100, for the image smoothing experiment, and 10 for the halftoning experiment. When scanning the image lattices for updating the system’s states, a maximum of 100 sweeps (raster-scan visits of the whole image) were performed at each temperature (respectively, demon energy) for the SA (respectively, MA) algorithm, and at each sampling stage of ~0. When less than 10% of the tested state transitions were accepted, the system was considered to be in equilibrium, and the temperature (or demon energy) was updated according to the annealing schedules in SA (or MA) ; for the ~0, a new initiaIization procedure was run. The algorithms stopped when a fraction
Letters A 205 f 19951377-382
of less than 1% of the proposed transitions was being accepted. Fig. 3 illustrates the convergence of ~0, as compared to SA and MA, for the two applications discussed. In the graphs, a step is counted after each change in temperature, for SA, or demon energy, for MA, and after each initialization stage for ~0. For comparison, we also show the curves corresponding to the ICM optimization, which is equivalent to performing just the initialization procedure of our algorithm. In summary, the microcanonical algorithm (,uO) has been introduced as an alternative to the annealing approaches to non-convex optimization problems. Such an algorithm follows the ideas introduced by Creutz [ 81 and first exploited by Barnard [ l] in the context of stochastic optimization. We have made an assessment of the performance of ~0 when employed in some simple image-processing applications, and we have shown its higher computational efficiency as compared to the annealing approaches. We would like to thank Professor F.G. Brady Moreira for introducing us to Creutz’s algorithm.
References [ 11 ST. Barnard, Int. J. Comput. Vis. 3 (1989) 17. [2] P. Carnevali, L. Coletti and S. Patamello, IBM J. Res. Dev. 29 (1985) 569. [ 3 J S. Geman and D. Geman, IEEE Trans. Pattern Anal. Mach. Intell. 6 (1984) 721. [4] L. Herault and R. Horaud, IEEE Trans. Pattern Anal. Mach. Intell. 15 (1993) 899. [5] S. Kirkpatrick, C.D. Gelatt and M.P. Vecchi, Science 220 (1983) 671. [6] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller and E. Teller, J. Chem. Phys. 21 (1953) 1087. [7] J. Besag, J. R. Statist. Sot. B 48(3) (1986) 259. 181 M. Creutz, Phys. Rev. Lett. 50 (1983) 1411.