c18-7 (779614), страница 3
Текст из файла (страница 3)
MEMreconstructions of test images using any of these entropy forms are virtuallyindistinguishable [5].By all available evidence, MEM seems to be neither more nor less than oneusefully nonlinear version of the general regularization scheme A + λB that we haveby now considered in many forms. Its peculiarities become strengths when appliedto the reconstruction from incomplete Fourier data of images that are expectedto be dominated by very bright point sources, but which also contain interestinglow-intensity, extended sources.
For images of some other character, there is noreason to suppose that MEM methods will generally dominate other regularizationschemes, either ones already known or yet to be invented.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).In the problem of incomplete Fourier image reconstruction, the typical Rjµhas the form exp(−2πikj · xµ ), where xµ is a two-dimensional vector in the imagespace and kµ is a two-dimensional wave-vector.
If an image contains strong pointsources, then the effect of setting unmeasured cj ’s to zero is to produce sideloberipples throughout the image plane. These ripples can mask any actual extended,low-intensity image features lying between the point sources. If, however, the slopeof G is smaller for small values of its argument, larger for large values, then ripplesin low-intensity portions of the image are relatively suppressed, while strong pointsources will be relatively sharpened (“superresolution”).
This behavior on the slopeof G is equivalent to requiring f 000 (u) < 0. For f(u) = u ln u, we in fact havef 000 (u) = −1/u2 < 0.In more picturesque language, the nonlinearity acts to “create” nonzero valuesfor the unmeasured ci ’s, so as to suppress the low-intensity ripple and sharpen thepoint sources.824Chapter 18.Integral Equations and Inverse TheoryAlgorithms for MEMThe goal is to find the vector bu that minimizes A + λB where in the notationof equations (18.5.5), (18.5.6), and (18.7.13),XA = |b − A · bu|2B=f(buµ )(18.7.18)Compared with a “general” minimization problem, we have the advantage thatwe can compute the gradients and the second partial derivative matrices (Hessianmatrices) explicitly,∇A = 2(AT · A · bu − AT · b)0[∇B]µ = f (buµ )∂2A∂buµ ∂buρ∂2B∂buµ ∂buρ= [2AT · A]µρ(18.7.19)00= δµρ f (buµ )It is important to note that while A’s second partial derivative matrix cannot bestored (its size is the square of the number of pixels), it can be applied to any vectorby first applying A, then AT .
In the case of reconstruction from incomplete Fourierdata, or in the case of convolution with a translation invariant point spread function,these applications will typically involve several FFTs. Likewise, the calculation ofthe gradient ∇A will involve FFTs in the application of A and AT .While some success has been achieved with the classical conjugate gradientmethod (§10.6), it is often found that the nonlinearity in f(u) = u ln u causesproblems. Attempted steps that give bu with even one negative value must be cut inmagnitude, sometimes so severely as to slow the solution to a crawl. The underlyingproblem is that the conjugate gradient method develops its information about theinverse of the Hessian matrix a bit at a time, while changing its location in the searchspace.
When a nonlinear function is quite different from a pure quadratic form, theold information becomes obsolete before it gets usefully exploited.Skilling and collaborators [6,7,10,11] developed a complicated but highly successful scheme, wherein a minimum is repeatedly sought not along a single searchdirection, but in a small- (typically three-) dimensional subspace, spanned by vectorsthat are calculated anew at each landing point.
The subspace basis vectors arechosen in such a way as to avoid directions leading to negative values. One of themost successful choices is the three-dimensional subspace spanned by the vectorswith components given bybµ [∇A]µe(1)µ = ue(2)bµ [∇B]µµ = uPPuµ ∂buρ )buρ [∇B]ρ uuµ ∂buρ )buρ [∇A]ρubµ ρ (∂ 2 A/∂bbµ ρ (∂ 2 A/∂b(3)qPqPeµ =−22bρ ([∇B]ρ )bρ ([∇A]ρ)ρuρu(18.7.20)(In these equations there is no sum over µ.) The form of the e(3) has some justificationif one views dot products as occurring in a space with the metric gµν = δµν /uµ,chosen to make zero values “far away”; see [6].Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).µ18.7 Maximum Entropy Image Restoration825µBecause the gradient directions ∇A and ∇B are separately available, it is possibleto combine the minimum search with a simultaneous adjustment of λ so as finally tosatisfy the desired constraint.
There are various further tricks employed.A less general, but in practice often equally satisfactory, approach is due toCornwell and Evans [12]. Here, noting that B’s Hessian (second partial derivative)matrix is diagonal, one asks whether there is a useful diagonal approximation toA’s Hessian, namely 2AT · A. If Λµ denotes the diagonal components of such anapproximation, then a useful step in bu would be1(∇A + λ∇B)(18.7.22)∆buµ = −Λµ + λf 00 (buµ )(again compare equation 10.7.4). Even more extreme, one might seek an approximation with constant diagonal elements, Λµ = Λ, so that∆buµ = −1(∇A + λ∇B)Λ + λf 00 (buµ )(18.7.23)Since AT · A has something of the nature of a doubly convolved point spreadfunction, and since in real cases one often has a point spread function with a sharpcentral peak, even the more extreme of these approximations is often fruitful. Onestarts with a rough estimate of Λ obtained from the Aiµ ’s, e.g.,+*X2(18.7.24)[Aiµ ]Λ∼iAn accurate value is not important, since in practice Λ is adjusted adaptively: If Λis too large, then equation (18.7.23)’s steps will be too small (that is, larger steps inthe same direction will produce even greater decrease in A + λB).
If Λ is too small,then attempted steps will land in an unfeasible region (negative values of ubµ ), or willresult in an increased A + λB. There is an obvious similarity between the adjustmentof Λ here and the Levenberg-Marquardt method of §15.5; this should not be toosurprising, since MEM is closely akin to the problem of nonlinear least-squaresuµ ) can be used tofitting. Reference [12] also discusses how the value of Λ + λf 00 (badjust the Lagrange multiplier λ so as to converge to the desired value of χ2 .All practical MEM algorithms are found to require on the order of 30 to 50iterations to converge.
This convergence behavior is not now understood in anyfundamental way.“Bayesian” versus “Historic” Maximum EntropySeveral more recent developments in maximum entropy image restorationgo under the rubric “Bayesian” to distinguish them from the previous “historic”methods.
See [13] for details and references.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Within the three-dimensional subspace, the three-component gradient and ninecomponent Hessian matrix are computed by projection from the large space, andthe minimum in the subspace is estimated by (trivially) solving three simultaneouslinear equations, as in §10.7, equation (10.7.4). The size of a step ∆bu is requiredto be limited by the inequalityX(∆buµ )2 /buµ < (0.1 to 0.5)U(18.7.21)826Chapter 18.Integral Equations and Inverse TheoryCITED REFERENCES AND FURTHER READING:Jaynes, E.T. 1976, in Foundations of Probability Theory, Statistical Inference, and StatisticalTheories of Science, W.L.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
