Expectation-maximization Algorithms For Image Processing .
Expectation-maximization algorithms for imageprocessing using multiscale models andmean- field theory, with applications to laser radarrange profiling and segmentationAndy TsaiLaboratory for Information and DecisionSystemsMassachusetts Institute of Technology77 Massachusetts AvenueCambridge, Massachusetts 02139E-mail: atsai@mit.eduJun ZhangDepartment of Electrical Engineering andComputer ScienceUniversity of Wisconsin–MilwaukeeMilwaukee, Wisconsin 53201E-mail: junzhang@csd.uwm.eduAlan S. WillskyLaboratory for Information and DecisionSystemsMassachusetts Institute of Technology77 Massachusetts AvenueCambridge, Massachusetts 02139E-mail: willsky@mit.edu1Abstract. We describe a new class of computationally efficient algorithms designed to solve incomplete-data problems frequently encountered in image processing and computer vision. The basis of this framework is the marriage of the expectation-maximization (EM) procedurewith two powerful methodologies. In particular, we have incorporatedoptimal multiscale estimators into the EM procedure to compute estimates and error statistics efficiently. In addition, mean-field theory (MFT)from statistical mechanics is incorporated into the EM procedure to helpsolve the computational problems that arise from our use of Markovrandom-field (MRF) modeling of the hidden data in the EM formulation.We have applied this algorithmic framework and shown that it is effectivein solving a wide variety of image-processing and computer-vision problems. We demonstrate the application of our algorithmic framework tosolve the problem of simultaneous anomaly detection, segmentation,and object profile estimation for noisy and speckled laser radar rangeimages. 2001 Society of Photo-Optical Instrumentation Engineers.[DOI: 10.1117/1.1385168]Subject terms: EM algorithm; mean-field theory; multiscale models; laser radar;image processing.Paper 200230 received June 8, 2000; accepted for publication Feb. 19, 2001.IntroductionEstimation problems for 2-D random fields arise in contextsranging from image processing to remote sensing. In somespecial cases—most notably spatially stationary statisticsover a rectangular grid with regularly spacedmeasurements—very efficient algorithms based on the fastFourier transform 共FFT兲 can be used to obtain optimal estimates and error statistics. However, when one strays fromthese cases, the complexity of optimal solutions for manypopular classes of stochastic models can become prohibitive with computational loads that do not scale well withproblem size.One approach that overcomes many of these problemsinvolves the use of multiscale stochastic models1,2 and associated multiscale estimation algorithms. These algorithmsdo not require stationarity, regularly shaped domains, orregularly spaced measurements. Moreover, these modelshave been shown to capture a rich class of random-fieldstatistical behavior, making them useful for a number ofimportant applications.In this paper we extend the applicability of this multiscale framework by marrying it with two other computationally powerful techniques, namely the expectationmaximization 共EM兲 algorithm and mean-field theory共MFT兲. The result is a methodology that can be applied toa greatly expanded range of applications, including texturesegmentation, simultaneous classification and gain correction of magnetic resonance imaging 共MRI兲 brain scans,3Opt. Eng. 40(7) 1287–1301 (July 2001)0091-3286/2001/ 15.00and a relaxed version of the Mumford-Shah variational approach to image segmentation. A number of these applications are described in some detail in Ref. 4, and Fig. 1depicts an example of the application to a texture segmentation problem.In particular, as described in Ref. 5, each of the individual textures shown in Fig. 1 can be well modeled withmultiscale models, allowing us to apply efficient multiscaleestimation procedures to problems such as restoring noisecorrupted versions of any one of these textures. However,the multiresolution algorithm by itself cannot solve theproblem of restoring the image shown in Fig. 1, as it consists of several different textures, nor can it directly solvethe problem of segmentation of such an image. However,by employing the EM concept of hidden data 共in this case,the segmentation map indicating which of the two woodtextures is visible in each pixel兲, we can indeed employmultiscale estimation as a core component of an algorithmthat produced the results 共restoration, segmentation, and associated error variances兲 shown in Fig. 1.Figure 2 depicts the conceptual structure of our methodology. Viewing the EM algorithm as the central structure ofan iterative procedure, we use multiscale estimation tosolve the so-called M step of the procedure and MFT toprovide an efficient approximation to the E step. In the nextsection we briefly describe the background topics that formthe components of the triad in Fig. 2. In Secs. 3 and 4 wethen develop our algorithm in the context of an important 2001 Society of Photo-Optical Instrumentation Engineers 1287
Tsai, Zhang, and Willsky: Expectation-maximization algorithms . . .in problems where the observation can be viewed as incomplete data. The MLE of X, denoted as X̂ ML , based on theincomplete observed data Y, is defined asX̂ ML arg max兵 log p 共 Y 兩 X 兲 其 ,共1兲XFig. 1 Motivational example: texture image segmentation and estimation.application, namely laser radar 共ladar兲 image processing,target segmentation, and range profiling. We have chosenthis application for several reasons. First, as we will see,the employment of our framework in this context involvesmore than just segmentation, and, as a result, using it as avehicle illustrates both the breadth of problems to which itcan be applied and how it is applied. Secondly, this particular application is of significant current interest, and, indeed,our results represent a substantial extension of previouswork on ladar anomaly rejection and background rangeplane profiling. Following the development of our approach, we present results in Sec. 5 demonstrating its efficacy in simultaneously rejecting anomalies, profiling nonplanar backgrounds, detecting and segmenting targets, andproviding range profiles of those targets. The paper thencloses with a brief discussion and conclusion.2 BackgroundAs we have indicated, our framework involves the synthesis of three distinct methodologies for statistical inference,and in this section we provide concise summaries of eachof these. More complete description of these topics can befound in some of the references 共e.g., Refs. 1,2, and 5–8兲.2.1 The EM ProcedureThe EM procedure is a powerful iterative technique suitedfor calculating the maximum-likelihood estimates 共MLEs兲where log p(Y 兩X) is the log likelihood of Y given X. Inmany applications, calculating the MLE is difficult becausethe log-likelihood function is highly nonlinear and not easily maximized. To overcome these difficulties, the EM algorithm introduces an auxiliary function Q 共along withsome auxiliary random variables兲 that has the same behavior as the log-likelihood function 共in that when the loglikelihood function increases, so does the auxiliary function兲 but is much easier to maximize.Central to the EM method is the judicious introductionof an auxiliary random quantity W with log likelihoodlog p(W兩X). The data W is referred to as the complete databecause it is more informative than Y. The complete data Wis not observed directly, but indirectly through Y via therelationship Y G(W), where G is a many-to-one mappingfunction. Those unobserved variables are referred to in theEM formulation as the hidden data and denoted by H. TheEM approach calculates X̂ ML through an iterative procedurein which the next iteration’s estimate of X is chosen tomaximize the expectation of log p(W兩X) given the incomplete data Y and the current iteration’s estimate of X. Theiteration consists of two steps: The E step computes the auxiliary functionQ共X兩X[k]兲 具log p共W兩X兲兩Y,X[k]典,共2兲where 具 典 represents expectation, and X [k] is the estimate of X from the k’th iteration. Often, this step reduces to calculating the expected value of H given Yand X [k] . The M step finds X [k 1] such thatX[k 1] arg max 兵Q共X兩X[k]兲其.XThe maximization here is with respect to X, the firstargument of the function Q. Intuitively, the M step isdesigned to use the expected value of the hidden dataH found in the E step as if it were measured data inorder to obtain the ML estimate of X. The EM algorithm can easily be extended to a maximum a posteriori 共MAP兲 estimator by imposing a prior model,p(X), on the estimated quantity X during this step.With this modification, the M step finds the X [k 1]such thatX[k 1] arg max 兵Q共X兩X[k]兲 log p共X兲其.共3兲XFig. 2 Conceptual representation of the algorithmic framework.1288 Optical Engineering, Vol. 40 No. 7, July 2001Since the EM algorithm is iterative, initial conditionsmust be given to start it. As it is guaranteed to convergeonly to a local maximum of the likelihood function, choos-
Tsai, Zhang, and Willsky: Expectation-maximization algorithms . . .ing initial conditions requires some care. The E and the Msteps are evaluated repeatedly until convergence, typicallyspecified as when 储 X [k 1] X [k] 储 for some small 0.2.2 Mean-Field TheoryMFT, which provides a computationally efficient procedurefor approximate calculation of expectations of large Markov random fields 共MRFs兲, has its origins in statistical mechanics and is concerned with finding the mean field energyof an MRF. In particular, let L denote a two-dimensionalrectangular lattice. Let u 兵 u i ,i苸L其 be a MRF taking onthe well-known Gibbs distribution 共see, e.g., Ref. 8兲 givenbyexp关 U 共 u 兲兴. p共 u 兲 The constant 0 is used to adjust the strength of theenergy function U(u) defined asU共 u 兲 兺ũ exp关 U 共 ũ 兲兴 .For simplicity of exposition, pixel interactions are assumedto be at most pairwise. The energy function can now bewritten asU共 u 兲 兺i苸L冋V clq共 u i 兲 12兺l苸Ni册V clq共 u i ,u l 兲 ,共4兲where the singleton clique V clq( ) and the doubleton cliqueV clq( , ) are the clique potentials for a single site and apair of neighboring sites 共horizontal or vertical兲, respectively. Finally, as a result of our simplifying assumption,Ni denotes the set of first-order neighbors of pixel i.MFT is used to compute an approximation to the meanof the field u, i.e., to find具 u i典 兺 u i p 共 u 兲 兺 u i 兺 p 共 u 兲uui sion of Eq. 共5兲. In particular, for each site i and each possible value of u i , we approximate this weight by replacingthe values of each of the u j , j i, by its mean value. Thatis, we define the mean field energy at site i asU mfi 共 u i 兲 U 共 u 兲 兩 u j 具 u j 典 ,V clq共 u 兲兺clqwhere the V clq(u) are the clique potentials used to capturethe interactions between pixels. The normalization constant is defined as Fig. 3 An example of a multiscale quadtree.兺u u iiL\ i兺 L\ i exp关 U 共 u 兲兴 共5兲for each i苸L 共where L\ i denotes the set of all pixels excluding i). It is clear from 共5兲 that a large number of configurations associated with the energy function need to beevaluated in order to calculate this expectation, thus making its precise calculation both complex and computationally impractical. To overcome this, MFT approximates theweight on each possible value of u i in the rightmost expres-mf mf j i U i 共 u i 兲 R i 共 具 u L i 典 兲 ,共6兲 where U mfi (u i ) includes all the clique potentials involvingthe value at site i and is given by U mfi 共 u i 兲 V clq共 u i 兲 兺l苸NiV clq共 u i , 具 u l 典 兲 ,共7兲 and R mfi ( 具 u L i 典 ) consists of terms involving only means atsites other than i. Further we use the mean-field approximation for the mean at each site:具u j典 兺 u juj mfj exp关 U mfj 共 u j 兲兴 mfj,共8兲 exp关 U mf兺j 共 ũ j 兲兴 .ũjThe MFT formulation supposes that for a particularpixel i苸L, the influence of the energy field from otherpixels can be approximated by that of their statisticalmeans. Random fluctuations from the other pixels are neglected. Note that since U mfi (u i ) depends on the means具 u j 典 , j i, the approximations to 具 u i 典 in Eq. 共8兲 dependson the means at other pixels. Thus the set of equationscorresponding to Eq. 共8兲 for all i苸L represents a set ofsimultaneous algebraic equations for all of the 具 u i 典 . Iterative algorithms can then be used to find the solution to thisset of equations.Multiscale Stochastic Models and OptimalEstimationIn this subsection, we briefly review the multiscale statistical framework for modeling and estimating 2-D randomfields.1,2,5 Multiscale stochastic models for image processing are defined on a pyramidal data structure such as thequadtree shown in Fig. 3. Each level of the tree corresponds2.3Optical Engineering, Vol. 40 No. 7, July 2001 1289
Tsai, Zhang, and Willsky: Expectation-maximization algorithms . . .to a different scale of resolution in the representation of therandom field of interest. Define the index s to specify thenodes on the tree, s to represent the parent of node s, ands i to denote the q offspring of node s with i 1, . . . ,q(q 4 for a quadtree兲. The scale of node s is denoted asm(s) and is numbered sequentially from coarse to fineresolution. A multiscale model with state X(s)苸Rn isspecified by a recursion of the formX共 s 兲 A共 s 兲 X共 s 兲 B共 s 兲 W共 s 兲共9兲with X(0) N(0,P o ), W(s) N(0,I), and W(s)苸R .关The notation N( , ) represents a Gaussian random vector x with mean and variance .兴 The matrices A(s) andB(s) are the parameters of the tree model, appropriatelychosen to represent the random field of interest. The rootnode of the tree X(0) with prior covariance P o provides theinitial condition for starting the coarse-to-fine recursion.The driving white noise W(s) is independent of the initialcondition X(0). The states at a given scale can be thoughtof as information about the process at that level of the tree.As X(s) evolves from coarse to fine scale, the termA(s)X(s ) predicts the states of the process for the nextfiner scale while the term B(s)W(s) adds new details tothese states.The measurement model associated with this multiscalestochastic model ismY共 s 兲 C共 s 兲 X共 s 兲 V共 s 兲共10兲with V(s) N(0,R(s)), R(s) the covariance of V(s), andC(s) the measurement matrix used to describe the nature ofthe measurement process. The general multiscale estimation framework allows easy and direct fusion of measurements at all scales 共see, for example, Ref. 9兲. For the application on which we focus here, all measurements will bescalar observations at the finest scale, representing pointmeasurements of the random field of interest.The algorithm for optimal multiscale estimation consistsof two steps: a fine-to-coarse sweep 共which resembles thetime-series Kalman filtering algorithm兲 followed by acoarse-to-fine sweep 共which is analogous to the RauchTung-Striebel smoothing algorithm兲. Among the most important properties of this algorithm is that it produces bothoptimal estimates and the covariance of the errors in theseestimates with total complexity that is O(k 3 n), where n isthe number of fine-scale pixels and k is the dimension ofthe state X(s). Thus the computation load for this algorithm grows linearly with problem size. This should be contrasted with the complexity of other estimation approaches共e.g., those based on MRF models兲, which have complexitythat is O(n ) for values of 1 for calculation of theestimates, and typically even greater complexity 共often prohibitively high兲 for the computation of error covariances.Of course, the utility of the multiscale approach also depends on having models with modest state dimension k.Demonstration that such models do exist for a wide varietyof random fields can be found in Refs. 5,10,11. Of specificrelevance to this paper are the multiscale models developed1290 Optical Engineering, Vol. 40 No. 7, July 2001Fig. 4 Down-looking geometry for collecting 3-D images of groundbased scenes.in Ref. 10 corresponding to so-called thin-plate and thinmembrane random surfaces. A brief summary of thesemodels is presented in Sec. 7.1.3Probabilistic Models for Laser Radar RangeProfiling and SegmentationAs illustrated in Fig. 4, and as described in detail in Refs.12,13, raster-scanned and range-resolved ladar imagingprovides the capability for producing 3-D images ofground-based objects from airborne platforms using adown-looking geometry. Utilization of such images 关e.g.,for automatic target recognition 共ATR兲兴 involves severalbasic functions. First, of course, is dealing with the uncertainties and sources of error in the ladar measurements. Asdescribed in Ref. 12, these errors consist of relatively smallperturbations due to noise and occasional large anomaliesdue to deep speckle fades. Since the frequency of anomalies can be significant, taking account of them is a seriousissue. Second, accurate localization and profiling of anytarget present in a ladar range image requires estimation ofthe unknown and spatially varying range of the ground.Detection and profiling of the target then requires segmenting out the region over which the range measurementsstand out from the background, indicating the presence ofthe target.As a first-generation preprocessor for an ATR system,Green and Shapiro12 developed an algorithm that estimatesthe background range profile for a ladar range image corrupted by range anomalies. Only the problem of estimatingthe background range profile is addressed in their algorithm. The idea is that pixels that do not belong to thebackground are grouped together as either range anomaliesor target range measurements. Their algorithm is based onan idealized single-pixel statistical model for the ladarrange measurements together with the assumption that thebackground range profile is a plane parameterized by theelevation and azimuth range slopes and by the range intercept. The EM procedure is employed in Green and Shapiro’s algorithm. The E step localizes and suppresses therange anomalies; the M step finds the ML estimate of thethree parameters that characterize the planar backgroundrange profile. Simulation studies have confirmed that Greenand Shapiro’s technique is computationally simple withgood range anomaly suppression capabilities.The application of our EM–multiscale-MFT estimationframework to the ladar range imaging problem extendsGreen and Shapiro’s work in several fundamental and prac-
Tsai, Zhang, and Willsky: Expectation-maximization algorithms . . .tically significant ways. First, rather than modeling thebackground surface as planar, we allow smooth spatialvariation of the background, using a multiresolution surfacemodel analogous to the so-called thin-membrane surfacemodel 共see Sec. 7.1兲. Employing this multiresolution modelallows us to use the very fast estimation procedure described in the previous section. In addition to estimating thebackground range profile, the approach we describe simultaneously performs two other significant functions, namelytarget segmentation and target range profile estimation. Toaccomplish this, we use a second multiscale model, in thiscase one corresponding to a so-called thin-plate prior, tomodel the target range profile within the image, togetherwith several important sets of hidden variables. The firsttwo of these sets, namely, the sets of pixels that correspondto anomalous background and target range measurements,are analogous to those used in Ref. 12. The other hiddenvariable, namely, the set of pixels segmenting the targetand background regions and explicitly identifying the lineof pixels corresponding to ground attachment 共i.e., wherethe target meets the ground兲, allows us to perform targetsegmentation and profiling.A uniformly spaced raster scan is used to collect a rangeimage consisting of n pixels. We denote the 2-D latticeconsisting of these n pixels as L. Let y 兵 y i ,i苸L其 denotethe measured ladar range values, whose background rangeprofile is denoted as x bg 兵 x bgi ,i苸L其 and whose targetrange pro
Expectation-maximization algorithms for image processing using multiscale models and . ML, based on the incomplete observed data Y, is defined as . only to a local maximum of the likelihood function, choos-Fig. 1 M
Bruksanvisning för bilstereo . Bruksanvisning for bilstereo . Instrukcja obsługi samochodowego odtwarzacza stereo . Operating Instructions for Car Stereo . 610-104 . SV . Bruksanvisning i original
Expectation Maximization We first focus on the analysis of the convergenceproperties of the Expectation-Maximization (EM) algorithm. Con-sider a probabilistic model of observed data x which uses latent variables z. The log-likelihood (objective function
Machine Learning for Computer Vision Expectation-Maximization EM is an elegant and powerful method for MLE problems with latent variables Main idea: model parameters and latent variables are estimated iteratively, where average over the latent variables (expectation) A typical exam
2. Expectation Maximization (EM) Iteration. A maximum likelihood (ML) method for image reconstruction based on Poisson data was introduced by Shepp and Vardi [42] in 1982 for image reconstruction in emission tomography. In fact, this algorithm was originally proposed by Richardson [37]
10 tips och tricks för att lyckas med ert sap-projekt 20 SAPSANYTT 2/2015 De flesta projektledare känner säkert till Cobb’s paradox. Martin Cobb verkade som CIO för sekretariatet för Treasury Board of Canada 1995 då han ställde frågan
service i Norge och Finland drivs inom ramen för ett enskilt företag (NRK. 1 och Yleisradio), fin ns det i Sverige tre: Ett för tv (Sveriges Television , SVT ), ett för radio (Sveriges Radio , SR ) och ett för utbildnings program (Sveriges Utbildningsradio, UR, vilket till följd av sin begränsade storlek inte återfinns bland de 25 största
Hotell För hotell anges de tre klasserna A/B, C och D. Det betyder att den "normala" standarden C är acceptabel men att motiven för en högre standard är starka. Ljudklass C motsvarar de tidigare normkraven för hotell, ljudklass A/B motsvarar kraven för moderna hotell med hög standard och ljudklass D kan användas vid
Zoo Animal Nutrition IV Zoo Animal Nutrition IV (2009) was edited by M. Clauss, A. Fidgett, G. Janssens, J.-M. Hatt, T. Huisman, J. Hummel, J. Nijboer, A. Plowman. Filander Verlag, Fürth ISBN-13: 978-3-930831-72-2 To obtain a copy of the book, contact Filander Verlag at info@filander.de Dierenfeld, E. S. Conservation collaborations: nutrition .
