Bayesian estimation problems involving Gaussian distributions often result in linear estimation techniques. Nevertheless, there are no general statements as to whether the linearity of the Bayesian estimator is restricted to the Gaussian case. The two common strategies for non-Gaussian models are either finding the best linear estimator or numerically evaluating the Bayesian estimator by Monte Carlo methods. In this paper, we focus on Bayesian interpolation of non-Gaussian first-order autoregressive (AR) processes where the driving innovation can admit any symmetric infinitely divisible distribution characterized by the Lévy-Khintchine representation theorem. We redefine the Bayesian estimation problem in the Fourier domain with the help of characteristic forms. By providing analytic expressions, we show that the optimal interpolator is linear for all symmetric -stable distributions. The Bayesian interpolator can be expressed in a convolutive form where the kernel is described in terms of exponential splines. We also show that the limiting case of Lévy-type AR(1) processes, the system of which has a pole at the origin, always corresponds to a linear Bayesian interpolator made of a piecewise linear spline, irrespective of the innovation distribution. Finally, we show the two mentioned cases to be the only ones within the family for which the Bayesian interpolator is linear.

Achieving accurate interpolation is an important requirement for many signal-processing applications. While nearest-neighbor and linear interpolation methods are popular due to their native GPU support, they unfortunately result in severe undesirable artifacts. Better interpolation methods are known but lack a native GPU support. Yet, a particularly attractive one is prefiltered cubic-spline interpolation. The signal it reconstructs from discrete samples has a much higher fidelity to the original data than what is achievable with nearest-neighbor and linear interpolation. At the same time, its computational load is moderate, provided a sequence of two operations is applied: first, prefilter the samples, and only then reconstruct the signal with the help of a B-spline basis. It has already been established in the literature that the reconstruction step can be implemented efficiently on a GPU. This article focuses on an efficient GPU implementation of the prefilter, on how to apply it to multidimensional samples (e.g. RGB color images), and on its performance aspects.

We present a new class of continuously defined parametric snakes using a special kind of exponential splines as basis functions. We have enforced our bases to have the shortest possible support subject to some design constraints to maximize efficiency. While the resulting snakes are versatile enough to provide a good approximation of any closed curve in the plane, their most important feature is the fact that they admit ellipses within their span. Thus, they can perfectly generate circular and elliptical shapes. These features are appropriate to delineate cross sections of cylindrical-like conduits and to outline bloblike objects. We address the implementation details and illustrate the capabilities of our snake with synthetic and real data.

Edge-preserving smoothers need not be taxed by a severe computational cost. We present, in this paper, a lean algorithm that is inspired by the bi-exponential filter and preserves its structure-a pair of one-tap recursions. By a careful but simple local adaptation of the filter weights to the data, we are able to design an edge-preserving smoother that has a very low memory and computational footprint while requiring a trivial coding effort. We demonstrate that our filter (a bi-exponential edge-preserving smoother, or BEEPS) has formal links with the traditional bilateral filter. On a practical side, we observe that the BEEPS also produces images that are similar to those that would result from the bilateral filter, but at a much-reduced computational cost. The cost per pixel is constant and depends neither on the data nor on the filter parameters, not even on the degree of smoothing.

We propose an active contour (a.k.a. snake) that takes the shape of an ellipse. Its evolution is driven by surface terms made of two contributions: the integral of the data over an inner ellipse, counterbalanced by the integral of the data over an outer elliptical shell. We iteratively adapt the active contour to maximize the contrast between the two domains, which results in a snake that seeks elliptical bright blobs. We provide analytic expressions for the gradient of the snake with respect to its defining parameters, which allows for the use of efficient optimizers. An important contribution here is the parameterization of the ellipse which we define in such a way that all parameters have equal importance; this creates a favorable landscape for the proceedings of the optimizer. We validate our construct with synthetic data and illustrate its use on real data as well.

Interpolation is the means by which a continuously defined model is fit to discrete data samples. When the data samples are exempt of noise, it seems desirable to build the model by fitting them exactly. In medical imaging, where quality is of paramount importance, this ideal situation unfortunately does not occur. In this paper, we propose a scheme that improves on the quality by specifying a tradeoff between fidelity to the data and robustness to the noise. We resort to variational principles, which allow us to impose smoothness constraints on the model for tackling noisy data. Based on shift-, rotation-, and scale-invariant requirements on the model, we show that the L(p)-norm of an appropriate vector derivative is the most suitable choice of regularization for this purpose. In addition to Tikhonov-like quadratic regularization, this includes edge-preserving total-variation-like (TV) regularization. We give algorithms to recover the continuously defined model from noisy samples and also provide a data-driven scheme to determine the optimal amount of regularization. We validate our method with numerical examples where we demonstrate its superiority over an exact fit as well as the benefit of TV-like nonquadratic regularization over Tikhonov-like quadratic regularization.

In the field of image segmentation, most level-set-based active-contour approaches take advantage of a discrete representation of the associated implicit function. We present in this paper a different formulation where the implicit function is modeled as a continuous parametric function expressed on a B-spline basis. Starting from the active-contour energy functional, we show that this formulation allows us to compute the solution as a restriction of the variational problem on the space spanned by the B-splines. As a consequence, the minimization of the functional is directly obtained in terms of the B-spline coefficients. We also show that each step of this minimization may be expressed through a convolution operation. Because the B-spline functions are separable, this convolution may in turn be performed as a sequence of simple 1-D convolutions, which yields an efficient algorithm. As a further consequence, each step of the level-set evolution may be interpreted as a filtering operation with a B-spline kernel. Such filtering induces an intrinsic smoothing in the algorithm, which can be controlled explicitly via the degree and the scale of the chosen B-spline kernel. We illustrate the behavior of this approach on simulated as well as experimental images from various fields.

Because more output data must be created than is available from the input, magnification is an ill-posed problem. Traditional magnification relies on resampling an interpolation model at the appropriate rate; unfortunately, this simple solution is blind to the presence of the analog filter that was implicitly present when the samples of the function to be magnified were acquired. Consistent resampling has been introduced to take this into account, but it turns out that this solution is still under-constrained. In this paper, we propose regularization as a way to devise a deterministic magnification method that fully satisfies consistency constraints in the absence of noise, and at the same time that produces an output that best fulfills a wide class of criteria for regularity. Contrarily to many other methods, ours has been designed without ever leaving the continuous domain. We conduct experiments that show the benefit of our approach.

A snakuscule (a minuscule snake) is the simplest active contour that we were able to design while keeping the quintessence of traditional snakes: an energy term governed by the data, and a regularization term. Our construction is an area-based snake, as opposed to curve-based snakes. It is parameterized by just two points, thus further easing requirements on the optimizer. Despite their ultimate simplicity, snakuscules retain enough versatility to be employed for solving various problems such as cell counting and segmentation of approximately circular features. In this paper, we detail the design process of a snakuscule and illustrate its usefulness through practical examples. We claim that our didactic intentions are well served by the simplicity of snakuscules.

Traditional snakes, or active contours, are planar parametric curves. Their parameters are determined by optimizing the weighted sum of three energy terms: one depending on the data (typically on the integral of its gradient under the curve, or on its integral over the area enclosed by the curve), one monitoring the shape of the curve (typically promoting its smoothness, or regularizing ambiguous solutions), and one incorporating prior knowledge (typically favoring a given shape). We present in this paper a snake that we designed to be as simple as possible without losing too many of the characteristics of more complicated, fuller versions. It retains an area data term and requires regularization to avoid an ill-posed optimization problem. It is parameterized by just two points, thus further easing requirements on the optimizer. Despite its extreme simplicity, this active contour can efficiently solve a variety of problems such as cell counting and segmentation of approximately circular features.

We present a simple, original method to improve piecewise-linear interpolation with uniform knots: we shift the sampling knots by a fixed amount, while enforcing the interpolation property. We determine the theoretical optimal shift that maximizes the quality of our shifted linear interpolation. Surprisingly enough, this optimal value is nonzero and close to 1/5. We confirm our theoretical findings by performing several experiments: a cumulative rotation experiment and a zoom experiment. Both show a significant increase of the quality of the shifted method with respect to the standard one. We also observe that, in these results, we get a quality that is similar to that of the computationally more costly "high-quality" cubic convolution.

Every now and then, a new design of an interpolation kernel appears in the literature. While interesting results have emerged, the traditional design methodology proves laborious and is riddled with very large systems of linear equations that must be solved analytically. We propose to ease this burden by providing an explicit formula that can generate every possible piecewise-polynomial kernel given its degree, its support, its regularity, and its order of approximation. This formula contains a set of coefficients that can be chosen freely and do not interfere with the four main design parameters; it is thus easy to tune the design to achieve any additional constraints that the designer may care for.

We address the task of rendering by ray tracing the isosurface of a high-quality continuous model of volumetric discrete and regular data. Based on first principles, we identify the quadratic B-spline as the best model for our purpose. The nonnegativity of this basis function allows us to confine the potential location of the isosurface within a binary shell. We then show how to use the space-embedding property of splines to further shrink this shell to essentially a single voxel width. Not all rays traced through a given shell voxel intersect the isosurface; many may only graze it, especially when the ray-tracing vantage point is close to or within the volume to be rendered. We propose an efficient heuristic to detect those cases. We present experiments to support our claims.

We present a simple, original method to improve piecewise linear interpolation with uniform knots. We shift the sampling knots by a fixed amount, while enforcing the interpolation property. Thanks to a theoretical analysis, we determine the optimal shift that maximizes the quality of our shifted linear interpolation. Surprisingly enough, this optimal value is nonzero and it is close to 1/5. We confirm our theoretical findings by performing a cumulative rotation experiment, which shows a significant increase of the quality of the shifted method with respect to the standard one. Most interesting is the fact that we get a quality similar to that of high-quality cubic convolution at the computational cost of linear interpolation.

We address the task of rendering by ray tracing the isosurface of a high-quality continuous spline model of volumetric discrete and regular data. By expressing the spline model as a sum of nonnegative B-splines, we are able to confine the potential location of the isosurface within a thin binary shell. We then show how to use the space-embedding property of splines to further shrink this shell to essentially a single-voxel width. We also propose a new illumination model that highlights the outline of the rendered isosurface, which provides for a sensitive test of the perceived quality of the rendering. We present experiments to support our claims, along with an efficient algorithm to compute simultaneously an array of B-splines and of its derivatives.

We consider the problem of interpolating a signal using a linear combination of shifted versions of a compactly-supported basis function phi(x). We first give the expression for the cases of phi's that have minimal support for a given accuracy (also known as "approximation order"). This class of functions, which we call maximal-order-minimal-support functions (MOMS) is made of linear combinations of the B-spline of the same order and of its derivatives. We provide an explicit form of the MOMS that maximizes the approximation accuracy when the step-size is small enough. We compute the sampling gain obtained by using these optimal basis functions over the splines of the same order. We show that it is already substantial for small orders and that it further increases with the approximation order L. When L is large, this sampling gain becomes linear; more specifically, its exact asymptotic expression is 2/(pie)L. Since the optimal functions are continuous, but not differentiable, for even orders, and even only piecewise continuous for odd orders, our result implies that regularity has little to do with approximating performance. These theoretical findings are corroborated by experimental evidence that involves compounded rotations of images.

We define texture mapping as an optimization problem for which the goal is to preserve the maximum amount of information in the mapped texture. We derive a solution that is optimal in the least-squares sense and that corresponds to the pseudo-inverse of the texture-mapping transformation. In practice, a first-order approximation of the least-squares solution is used as an initial estimate for the mapped texture. This initial solution is refined by successive approximation to yield the least-squares optimal result. In essence, the proposed multi-pass method acts like an adaptive antialiasing filter.

The most essential ingredient of interpolation is its basis function. We have shown in previous papers that this basis need not be necessarily interpolating to achieve good results. On the contrary, several studies have confirmed that non-interpolating bases, such as B-splines and O-moms, perform best. This opens up a much wider choice of basis functions. We give to the designer the tools that will allow him to characterize this enlarged space of functions. In particular, he will be able to specify up-front the four most important parameters for image processing: degree, support, regularity, and order. The theorems presented then allow him to refine his design by dealing with additional coefficients that can be selected freely, without interfering with the main design parameters.

We propose a new method for the intermodal registration of images using a criterion known as mutual information. Our main contribution is an optimizer that we specifically designed for this criterion. We show that this new optimizer is well adapted to a multiresolution approach because it typically converges in fewer criterion evaluations than other optimizers. We have built a multiresolution image pyramid, along with an interpolation process, an optimizer, and the criterion itself, around the unifying concept of spline-processing. This ensures coherence in the way we model data and yields good performance. We have tested our approach in a variety of experimental conditions and report excellent results. We claim an accuracy of about a hundredth of a pixel under ideal conditions. We are also robust since the accuracy is still about a tenth of a pixel under very noisy conditions. In addition, a blind evaluation of our results compares very favorably to the work of several other researchers.

Echo-planar imaging (EPI) is a fast nuclear magnetic resonance imaging (MRI) method. Unfortunately, local magnetic field inhomogeneities induced mainly by the subject's presence cause significant geometrical distortion, predominantly along the phase-encoding direction, which must be undone to allow for meaningful further processing. So far, this aspect has been too often neglected. In this paper, we suggest a new approach using an algorithm specifically developed for the automatic registration of distorted EPI images with corresponding anatomically correct MRI images. We model the deformation field with splines, which gives us a great deal of flexibility, while comprising the affine transform as a special case. The registration criterion is least squares. Interestingly, the complexity of its evaluation does not depend on the resolution of the control grid. The spline model gives us good accuracy thanks to its high approximation order. The short support of splines leads to a fast algorithm. A multiresolution approach yields robustness and additional speedup. The algorithm was tested on real as well as synthetic data, and the results were compared with a manual method. A wavelet-based Sobolev-type random deformation generator was developed for testing purposes. A blind test indicates that the proposed automatic method is faster, more reliable, and more precise than the manual one.

Based on the theory of approximation, this paper presents a unified analysis of interpolation and resampling techniques. An important issue is the choice of adequate basis functions. We show that, contrary to the common belief, those that perform best are not interpolating. By opposition to traditional interpolation, we call their use generalized interpolation; they involve a prefiltering step when correctly applied. We explain why the approximation order inherent in any basis function is important to limit interpolation artifacts. The decomposition theorem states that any basis function endowed with approximation order can be expressed as the convolution of a B-spline of the same order with another function that has none. This motivates the use of splines and spline-based functions as a tunable way to keep artifacts in check without any significant cost penalty. We discuss implementation and performance issues, and we provide experimental evidence to support our claims.

We extend the classical interpolation method to generalized interpolation. This extension is done by replacing the interpolating function by a non-interpolating function that is applied to prefiltered data, in order to preserve the interpolation condition. We show, both theoretically and practically, that this approach performs much better than classical methods, for the same computational cost.

Due to magnetic field inhomogeneities, EPI images are geometrically distorted, predominantly along the phase-encoding direction. Currently, the distortion is either ignored or compensated manually using a warping function defined through a set of landmarks. We propose an automatic method to unwarp the geometric distortion of EPI images by registering them with corresponding undistorted anatomical MRI images. We show that cubic splines are optimal interpolating functions for both landmark interpolation and approximation. We will consequently use the same warping space in our algorithm which replaces landmarks by an image difference criterion. B-splines are used as generating functions, which leads to a fast and accurate computation. Multiresolution gives robustness and additional speedup. The algorithm performance was evaluated using both real and synthetic data and was found superior to the manual method.

We propose a new optimizer in the context of multimodal image registration. The optimized criterion is the mutual information between the images to be align. This criterion requires that their joint histogram be available. For its computation, we introduce differentiable and separable Parzen windows that satisfy the partition of unity. Along with a continuous model of the images based on splines, this allows us to derive exact and tractable expressions for the gradient and the Hessian of the criterion. Then, we develop an optimizer based on the Marquardt-Levenberg (1963) strategy. Our new optimizer is specific to mutual information, in the same sense that Marquardt-Levenberg is specific to least-squares. We show that our optimizer is particularly well-adapted to an iterative coarse-to fine approach. We validate its accuracy by comparing its performance to that of several results available in the literature.

We present an automatic subpixel registration algorithm that minimizes the mean square intensity difference between a reference and a test data set, which can be either images (two-dimensional) or volumes (three-dimensional). It uses an explicit spline representation of the images in conjunction with spline processing, and is based on a coarse-to-fine iterative strategy (pyramid approach). The minimization is performed according to a new variation (ML*) of the Marquardt-Levenberg algorithm for nonlinear least-square optimization. The geometric deformation model is a global three-dimensional (3-D) affine transformation that can be optionally restricted to rigid-body motion (rotation and translation), combined with isometric scaling. It also includes an optional adjustment of image contrast differences. We obtain excellent results for the registration of intramodality positron emission tomography (PET) and functional magnetic resonance imaging (fMRI) data. We conclude that the multiresolution refinement strategy is more robust than a comparable single-stage method, being less likely to be trapped into a false local optimum. In addition, our improved version of the Marquardt-Levenberg algorithm is faster.

We investigate aspects of multi-modal image registration based on a new criterion named mutual information (or sometimes Shannon information). This criterion is intensity-based and requires no landmarks; hence, its application can be automated without resorting to segmentation. We present a form amenable to derivation with respect to the geometric transformation parameters (affine transformation). This form involves Parzen windows; we explore the dependence of the registration accuracy on these windows and propose that they be tuned to each resolution level in a pyramid approach. We conduct experiments and show that both the window width and the number of windows is relevant. In addition, we show that it is beneficial to use a spline-based high-order interpolation scheme for applying the geometric transformation.

We present examples of a new type of wavelet basis functions that are orthogonal across shifts but not across scales. The analysis functions are piecewise linear while the synthesis functions are polynomial splines of degree n (odd). The approximation power of these representations is essentially as good as that of the corresponding Battle-Lemarie orthogonal wavelet transform, with the difference that the present wavelet synthesis filters have a much faster decay. This last property, together with the fact that these transformations are almost orthogonal, may be useful for image coding applications.

Presents a general framework for the fast, high quality implementation of geometric affine transformations of images (p=2) or volumes (p=3), including rotations and scaling. The method uses a factorization of the p/spl times/p transformation matrix into p+1 elementary matrices, each affecting one dimension of the data only. This yields a separable implementation through an appropriate sequence of 1-D affine transformations (scaling+translation). Each elementary transformation is implemented in an optimal least squares sense using a polynomial spline signal model. The authors consider various matrix factorizations and compare their method with the conventional nonseparable interpolation approach. The new method provides essentially the same quality results and at the same time offers significant speed improvement.

The use of the wavelet transform to detect differences between sequentially acquired functional magnetic resonance images (fMRIs) is explored. A statistical data model is developed that makes use of the orthogonality and regularity conditions of the wavelets to achieve a signal decomposition into uncorrelated components, enabling application of standard parametric tests of significance on wavelet coefficients directly. This overcomes the problems associated with high intervoxel correlations in the spatial domain, and achieves economy in statistical testing by limiting the search for significant signal components to a subspace where the signal power is located. Thus, a smaller p-value adjustment for multiple testing is required, resulting in a lower detection threshold for a given overall level of statistical significance. For the fMRIs investigated, a 10:1 reduction in the number of statistical tests was achieved, and about 1% of the wavelet coefficients were significant (p<0.05 per volume), which then served to resynthesize the difference images by inverse wavelet transform.

We present an automatic sub-pixel registration algorithm that minimizes the mean square difference of intensities between a reference and a test data set (volumes or images). It uses spline processing, is based on a coarse-to-fine pyramid strategy, and performs minimization according to a variation of the iterative Marquardt-Levenberg (1963) scheme. The geometric deformation model is a general affine transformation that one may optionally restrict to a rigid-body (isometric scale, rotation and translation), procrustean (rotation and translation) or translational case; it also includes an optional parameter for the linear adaptation of intensity. We present several PET and fMRI experiments and show that this algorithm provides excellent results. We conclude that the multi-resolution refinement strategy is faster and more robust than a comparable single-scale one.

This paper focuses on the design of fast algorithms for rotating images and preserving high quality. The basis for the approach is a decomposition of a rotation into a sequence of one-dimensional translations. As the accuracy of these operations is critical, we introduce a general theoretical framework that addresses their design and performance. We also investigate the issue of optimality and present an improved least-square formulation of the problem. This approach leads to a separable three-pass implementation of a rotation using one-dimensional convolutions only. We provide explicit filter formulas for several continuous signal models including spline and bandlimited representations. Finally, we present rotation experiments and compare the currently standard techniques with the various versions of our algorithm. Our results indicate that the present algorithm in its higher-order versions outperforms all standard high-accuracy methods of which we are aware, both in terms of speed and quality. Its computational complexity increases linearly with the order of accuracy. The best-quality results are obtained with the sine-based algorithm, which can be implemented using simple one-dimensional FFTs.