c18-5 (779612), страница 2

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 2)

You don’t need to knowthe values of your boundary conditions, since B can have fewer rows than columns,as above; hopefully, your data will determine them. Of course, if you do know someboundary conditions, you can build these into B too.With all the proportionality signs above, you may have lost track of what actualvalue of λ to try first. A simple trick for at least getting “on the map” is to first try1 −33−100 −310 −126 −1019 −156−1 3 −126 −1520 −156 −16 −1520 −15 0 −1 ....H=. .0 −16 −15 0 ···00 −16 0 ··· 0 ···000−1 0 ···00000 ···0000000−16λ = Tr(AT · A)/Tr(H)0000−100000···············(18.5.16)where Tr is the trace of the matrix (sum of diagonal components). This choicewill tend to make the two parts of the minimization have comparable weights, andyou can adjust from there.As for what is the “correct” value of λ, an objective criterion, if you knowu − b|2 ) equalyour errors σi with reasonable accuracy, is to make χ2 (that is, |A · bto N , the number of measurements.

We remarked above on the twin acceptablechoices N ± (2N )1/2 . A subjective criterion is to pick any value that you like in theSample 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).µ=1812Chapter 18.Integral Equations and Inverse Theoryrange 0 < λ < ∞, depending on your relative degree of belief in the a priori and aposteriori evidence.

(Yes, people actually do that. Don’t blame us.)Two-Dimensional Problems and Iterative MethodseA(i − µ, j − ν) ⇐⇒ A(k)b(i,j) ⇐⇒ eb(k)ub(i,j) ⇐⇒ ue(k) (18.5.17)We also need a regularization or smoothing operator B and the derived H = BT · B.One popular choice for B is the five-point finite-difference approximation of theLaplacian operator, that is, the difference between the value of each point and theaverage of its four Cartesian neighbors. In Fourier space, this choice implies,eB(k)∝ sin2 (πk1 /M ) sin2 (πk2 /K)eH(k)∝ sin4 (πk1 /M ) sin4 (πk2 /K)(18.5.18)In Fourier space, equation (18.5.7) is merely algebraic, with solutionue(k) =eeb(k)A*(k)2ee|A(k)| + λH(k)(18.5.19)where asterisk denotes complex conjugation.

You can make use of the FFT routinesfor real data in §12.5.Turn now to the case where A is not translationally invariant. Direct solutionof (18.5.8) is now hopeless, since the matrix A is just too large. We need somekind of iterative scheme.One way to proceed is to use the full machinery of the conjugate gradientmethod in §10.6 to find the minimum of A + λB, equation (18.5.6).

Of the variousmethods in Chapter 10, conjugate gradient is the unique best choice because (i)it does not require storage of a Hessian matrix, which would be infeasible here,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).Up to now our notation has been indicative of a one-dimensional problem,b(xµ ).

However, all of the discussion easily generalizes to thefinding ub(x) or ubµ = uproblem of estimating a two-dimensional set of unknowns ubµκ , µ = 1, . . . , M, κ =1, . . . , K, corresponding, say, to the pixel intensities of a measured image. In thiscase, equation (18.5.8) is still the one we want to solve.In image processing, it is usual to have the same number of input pixels in ameasured “raw” or “dirty” image as desired “clean” pixels in the processed outputimage, so the matrices R and A (equation 18.5.5) are square and of size M K × M K.A is typically much too large to represent as a full matrix, but often it is either (i)sparse, with coefficients blurring an underlying pixel (i, j) only into measurements(i±few, j±few), or (ii) translationally invariant, so that A(i,j)(µ,ν) = A(i−µ, j−ν).Both of these situations lead to tractable problems.In the case of translational invariance, fast Fourier transforms (FFTs) are theobvious method of choice.

The general linear relation between underlying functionand measured values (18.4.7) now becomes a discrete convolution like equation(13.1.1). If k denotes a two-dimensional wave-vector, then the two-dimensionalFFT takes us back and forth between the transform pairs18.5 Linear Regularization Methods813and (ii) it does exploit gradient information, which we can readily compute: Thegradient of equation (18.5.6) is∇(A + λB) = 2[(AT · A + λH) · bu − AT · b](18.5.20)bu(k) + AT · bu(k+1) = [1 − (AT · A + λH)] · b(18.5.21)Here is a parameter that dictates how far to move in the downhill gradient direction.The method converges when is small enough, in particular satisfying0<<2max eigenvalue (AT · A + λH)(18.5.22)There exist complicated schemes for finding optimal values or sequences for ,see [7]; or, one can adopt an experimental approach, evaluating (18.5.6) to be surethat downhill steps are in fact being taken.In those image processing problems where the final measure of success issomewhat subjective (e.g., “how good does the picture look?”), iteration (18.5.21)sometimes produces significantly improved images long before convergence isachieved.

This probably accounts for much of its use, since its mathematicalconvergence is extremely slow. In fact, (18.5.21) can be used with H = 0, inwhich case the solution is not regularized at all, and full convergence would bedisastrous! This is called Van Cittert’s method and goes back to the 1930s. A numberof iterations the order of 1000 is not uncommon [7].Deterministic Constraints: Projections onto Convex SetsA set of possible underlying functions (or images) {bu} is said to be convex if,ub in the set, all the linearly interpolated combinationsfor any two elements bua and bub(1 − η)bua + ηb0≤η≤1(18.5.23)are also in the set. Many deterministic constraints that one might want to impose onthe solution bu to an inverse problem in fact define convex sets, for example:• positivity• compact support (i.e., zero value outside of a certain region)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).(cf. 18.5.8). Evaluation of both the function and the gradient should of course takeadvantage of the sparsity of A, for example via the routines sprsax and sprstxin §2.7.

We will discuss the conjugate gradient technique further in §18.7, in thecontext of the (nonlinear) maximum entropy method. Some of that discussion canapply here as well.The conjugate gradient method notwithstanding, application of the unsophisticated steepest descent method (see §10.6) can sometimes produce useful results,particularly when combined with projections onto convex sets (see below).

If thesolution after k iterations is denoted bu(k), then after k + 1 iterations we have814Chapter 18.Integral Equations and Inverse Theory|Pi bu−bu| ≤ |bua − bu|for all bua in Ci(18.5.24)While this definition sounds complicated, examples are very simple: A nonexpansiveprojection onto the set of positive bu’s is “set all negative components of bu equalto zero.” A nonexpansive projection onto the set of ub(x)’s bounded by uL (x) ≤ub(x) ≤ uU (x) is “set all values less than the lower bound equal to that bound, andset all values greater than the upper bound equal to that bound.” A nonexpansiveprojection onto functions with compact support is “zero the values outside of theregion of support.”The usefulness of these definitions is the following remarkable theorem: Let Cbe the intersection of m convex sets C1 , C2 , .

. . , Cm . Then the iteration(k+1)(k)b= (P1 P2 · · · Pm )buu(18.5.25)will converge to C from all starting points, as k → ∞. Also, if C is empty (thereis no intersection), then the iteration will have no limit point. Application of thistheorem is called the method of projections onto convex sets or sometimes POCS [7].A generalization of the POCS theorem is that the Pi ’s can be replaced bya set of Ti ’s,Ti ≡ 1 + βi (Pi − 1)0 < βi < 2(18.5.26)A well-chosen set of βi ’s can accelerate the convergence to the intersection set C.Some inverse problems can be completely solved by iteration (18.5.25) alone!For example, a problem that occurs in both astronomical imaging and X-raydiffraction work is to recover an image given only the modulus of its Fouriertransform (equivalent to its power spectrum or autocorrelation) and not the phase.Here three convex sets can be utilized: the set of all images whose Fourier transformhas the specified modulus to within specified error bounds; the set of all positiveimages; and the set of all images with zero intensity outside of some specified region.In this case the POCS iteration (18.5.25) cycles among these three, imposing eachconstraint in turn; FFTs are used to get in and out of Fourier space each time theFourier constraint is imposed.The specific application of POCS to constraints alternately in the spatial andFourier domains is also known as the Gerchberg-Saxton algorithm [8].

Характеристики

Список файлов книги