Suhail S. Saquib

ML parameter estimation for Markov random fields with applications to Bayesian tomography

Markov random fields (MRFs) have been widely used to model images in Bayesian frameworks for image reconstruction and restoration. Typically, these MRF models have parameters that allow the prior model to be adjusted for best performance. However, optimal estimation of these parameters(sometimes referred to as hyper parameters) is difficult in practice for two reasons: i) direct parameter estimation for MRFs is known to be mathematically and numerically challenging; ii)parameters can not be directly estimated because the true image cross section is unavailable.In this paper, we propose a computationally efficient scheme to address both these difficulties for a general class of MRF models,and we derive specific methods of parameter estimation for the MRF model known as generalized Gaussian MRF (GGMRF).The first section of the paper derives methods of direct estimation of scale and shape parameters for a general continuously valued MRF. For the GGMRF case, we show that the ML estimate of the scale parameter, sigma, has a simple closed-form solution, and we present an efficient scheme for computing the ML estimate of the shape parameter, p, by an off-line numerical computation of the dependence of the partition function on p.The second section of the paper presents a fast algorithm for computing ML parameter estimates when the true image is unavailable. To do this, we use the expectation maximization(EM) algorithm. We develop a fast simulation method to replace the E-step, and a method to improve parameter estimates when the simulations are terminated prior to convergence.Experimental results indicate that our fast algorithms substantially reduce computation and result in good scale estimates for real tomographic data sets.

Parallelizable Bayesian tomography algorithms with rapid, guaranteed convergence

Bayesian tomographic reconstruction algorithms generally require the efficient optimization of a functional of many variables. In this setting, as well as in many other optimization tasks, functional substitution (FS) has been widely applied to simplify each step of the iterative process. The function to be minimized is replaced locally by an approximation having a more easily manipulated form, e.g., quadratic, but which maintains sufficient similarity to descend the true functional while computing only the substitute. We provide two new applications of FS methods in iterative coordinate descent for Bayesian tomography. The first is a modification of our coordinate descent algorithm with one-dimensional (1-D) Newton-Raphson approximations to an alternative quadratic which allows convergence to be proven easily. In simulations, we find essentially no difference in convergence speed between the two techniques. We also present a new algorithm which exploits the FS method to allow parallel updates of arbitrary sets of pixels using computations similar to iterative coordinate descent. The theoretical potential speed up of parallel implementations is nearly linear with the number of processors if communication costs are neglected.

Model-based image reconstruction from time-resolved diffusion data

This paper addresses the issue of reconstructing the unknown field of absorption and scattering coefficients from time- resolved measurements of diffused light in a computationally efficient manner. The intended application is optical tomography, which has generated considerable interest in recent times. The inverse problem is posed in the Bayesian framework. The maximum a posteriori (MAP) estimate is used to compute the reconstruction. We use an edge-preserving generalized Gaussian Markov random field to model the unknown image. The diffusion model used for the measurements is solved forward in time using a finite-difference approach known as the alternating-directions implicit method. This method requires the inversion of a tridiagonal matrix at each time step and is therefore of O(N) complexity, where N is the dimensionality of the image. Adjoint differentiation is used to compute the sensitivity of the measurements with respect to the unknown image. The novelty of our method lies in the computation of the sensitivity since we can achieve it in O(N) time as opposed to O(N2) time required by the perturbation approach. We present results using simulated data to show that the proposed method yields superior quality reconstructions with substantial savings in computation.

A non-homogeneous MRF model for multiresolution Bayesian estimation

The popularity of Bayesian methods in image processing applications has generated great interest in image modeling. A good image model needs to be non-homogeneous to be able to adapt to the local characteristics of the different regions in an image. In the past however, such a formulation was difficult since it was not clear as to how to choose the parameters of the non-homogeneous model. But now motivated by results in maximum likelihood parameter estimation of MRF models, we formulate in this paper a non-homogeneous Markov random field (MRF) image model in the multiresolution framework. The advantage of the multiresolution framework is two fold: first, it makes it possible to estimate the parameters of the nonhomogeneous MRF at any resolution by using the image at the coarser resolution. Second, it yields multiresolution algorithms which are computationally efficient and more robust than their single resolution counterparts. Experimental results in tomographic image reconstruction and optical flow computation problems verify the superior modeling provided by the new model.

Provably convergent coordinate descent in statistical tomographic reconstruction

Statistical tomographic reconstruction algorithms generally require the efficient optimization of a functional. An algorithm known as iterative coordinate descent with Newton-Raphson updates (ICD/NR) has been shown to be much more computationally efficient than indirect optimization approaches based on the EM algorithm. However, while the ICD/NR algorithm has experimentally been shown to converge stably, no theoretical proof of convergence is known. We prove that a modified algorithm, which we call ICD functional substitution (ICD/FS), has guaranteed global convergence in addition to the computational efficiency of the ICD/NR. The ICD/FS method works by approximating the log likelihood at each pixel by an alternative quadratic functional. Experimental results show that the convergence speed of the globally convergent algorithm is nearly identical to that of ICD/NR.

Efficient ML estimation of the shape parameter for generalized Gaussian MRFs

A certain class of Markov random fields (MRF) known as generalized Gaussian MRFs (GGMRF) have been shown to yield good performance in modeling the a priori information in Bayesian image reconstruction and restoration problems. Though the ML estimate of temperature T of a GGMRF has a closed form solution, the optimal estimation of the shape parameter p is a difficult problem due to the intractable nature of the partition function. We present a tractable scheme for ML estimation of p by an off-line numerical computation of the log of the partition function. In image reconstruction or restoration problems, the image itself is not known. To address this problem, we use the EM algorithm to compute the estimates directly from the data. For efficient computation of the expectation step, we propose a fast simulation technique and a method to extrapolate the estimates when the simulations are terminated prematurely prior to convergence. Experimental results show that the proposed methods result in substantial savings in computation and superior quality images.

Tractable models and efficient algorithms for Bayesian tomography

Bayesian methods have proven to be powerful tools for computed tomographic reconstruction in realistic physical problems. However, Bayesian methods require that a number of modeling and computational problems be addressed. The paper summarizes a coherent system of statistical modeling and optimization techniques designed to facilitate efficient, unsupervised Bayesian emission and transmission tomographic reconstruction. New results are included on improved convergence behavior of these methods.

An efficient multiresolution algorithm for compensating density-dependent media blurring

The sharpness of a printed image may suffer due to the presence of material layers above and below the dye layers. These layers contribute to scattering and surface reflections that make the degradation in sharpness density-dependent. We present data that illustrate this effect, and model the phenomenon numerically. A digital non-linear sharpening filter is proposed to compensate for this density-dependent blurring. The support and shape of this filter is constrained to lie in a space spanned by a set of basis filters that can be computed efficiently. Burt and Adelsons Laplacian pyramid is used to develop an efficient scale-recursive algorithm in which the image is decomposed into the high-pass basis images in a fine-to-coarse scale sweep, and the sharpened image along with a local density image is subsequently synthesized by a coarse-to-fine scale sweep using these basis images. The local density image is employed, in combination with a scale dependent gain function, to modulate the high-pass basis images in a space-varying fashion. A robust method is proposed for the estimation of the gain functions directly from measured data. Experimental results demonstrate that the proposed algorithm successfully compensates for media-related density dependent blurring.

A real-time multiresolution algorithm for correcting distortions produced by thermal printers

As printing proceeds in a thermal printer, heat from previously printed lines of image data accumulates in the print head and alters the thermal state of the heating elements. The changing state of the heating elements manifests itself as errors in the printed image. We have modeled the heat diffusion within the thermal printer and the density response of the receiver medium to derive a computationally efficient inverse thermal printer model. In this model, the heat diffusion problem for the moving receiver is simplified by showing that a moving receiver is equivalent to a stationary receiver with higher conductivity. The thermal print head is modeled as having a finite number of discrete layers with differing time constants. The layer temperature updates can be decoupled and are time-recursive if expressed in relative rather than absolute temperatures. Decoupling allows the layers to be updated at multiple spatial and temporal resolutions. The inverse printer model then reduces to an elegant algorithm that comprises three interleaved recursions; namely, absolute temperature propagation from coarse-to-fine scale, energy propagation from fine-to-coarse scale and relative temperature update in time. Experimental results demonstrate that the proposed algorithm successfully corrects the errors introduced by thermal printers.