CONSTRAINED RESTORATION AND THE RECOVERY OF DISCONTINUITIES

The linear image restoration problem is to recover an original brightness distribution X0 given the blurred and noisy observations Y = KX0 + B, where K and B represent the point spread function and measurement error, respectively. This problem is typical of ill-conditioned inverse problems that frequently arise in low-level computer vision. A conventional method to stabilize the problem is to introduce a priori constraints on X0 and design a cost functional H(X) over images X, which is a weighted average of the prior constraints (regularization term) and posterior constraints (data term); the reconstruction is then the image X, which minimizes H. A prominent weakness in this approach, especially with quadratic-type stabilizers, is the difficulty in recovering discontinuities. One seeks an estimate of X0, which not only recovers the shape of the original image over smooth patches, for example, those that are planar or quadric, but also recovers sharp transitions between these components. We therefore examine prior smoothness constraints of the form phi(D(k)X), where phi(u) = -(1 + \u\)-1, and D(k) denotes a kth order derivative k = 1,2, or 3. The important attributes of phi are its concavity on (0, infinity) and its finite asymptotic behavior (lim(u--> infinity) phi(u) < infinity). Such constraints permit the recovery of discontinuities without introducing auxilliary variables for marking the location of jumps and suspending the constraints in their vicinity. (In fact, our optimization problem is equivalent to one involving a noninteracting "line process.") In this sense, discontinuities are addressed implicitly rather than explicitly. Selecting the parameters, especially the relative weight-lambda between the prior and posterior terms (the "smoothing parameter"), is also problematical. Moreover, in our view, there is a conspicuous absence of theoretical results on model validation, even for idealized X0. By exploiting the concavity of phi and assuming that X0 is an ideal (but prototypical) pattern, we calculate-lambda by requiring that Pr{X0 is-an-element-of W} almost-equal-to 1, where W is the set of coordinate-wise minima of H. This procedure then yields-lambda (actually an upper bound) as a function of the other model parameters, such as the noise variance and blur coefficients.