Abstract. This paper addresses the task of learning convex regularizers to guide the reconstruction of images from limited data. By imposing that the reconstruction be amplitude-equivariant, we narrow down the class of admissible functionals to those that can be expressed as a power of a seminorm. We then show that such functionals can be approximated to arbitrary precision with the help of polyhedral norms. In particular, we identify two dual parameterizations of such systems: (i) a synthesis (or atomic) form with an [Formula: see text]-penalty that involves some learnable dictionary; and (ii) an analysis form with an [Formula: see text]-penalty that involves a trainable regularization operator. After having provided geometric insights and proved that the two forms are universal, we propose an implementation that relies on a specific architecture (tight frame with a weighted [Formula: see text] penalty) that is easy to train. We illustrate its use for denoising and the reconstruction of biomedical images. We find that the proposed framework outperforms the sparsity-based methods of compressed sensing, while it offers essentially the same convergence and robustness guarantees.

Phase retrieval, a nonlinear problem prevalent in imaging applications, has been extensively studied using random models, some of which with i.i.d. sensing matrix components. While these models offer robust reconstruction guarantees, they are computationally expensive and impractical for real-world scenarios. In contrast, Fourier-based models, common in applications such as ptychography and coded diffraction imaging, are computationally more efficient but lack the theoretical guarantees of random models. Here, we introduce structured random models for phase retrieval that combine the efficiency of fast Fourier transforms with the versatility of random diagonal matrices. These models emulate i.i.d. random matrices at a fraction of the computational cost. Our approach demonstrates robust reconstructions comparable to fully random models using gradient descent and spectral methods. Furthermore, we establish that a minimum of two structured layers is necessary to achieve these structured-random properties. The proposed method is suitable for optical implementation and offers an efficient and robust alternative for phase retrieval in practical imaging applications.

State-of-the-art image reconstruction often relies on complex, abundantly parameterized deep architectures. We propose an alternative: a data-driven reconstruction method inspired by the classic Tikhonov regularization. Our approach iteratively refines intermediate reconstructions by solving a sequence of quadratic problems. These updates have two key components: (i) learned filters to extract salient image features; and (ii) an attention mechanism that locally adjusts the penalty of the filter responses. Our method matches leading plug-and-play and learned regularizer approaches in performance while offering interpretability, robustness, and convergent behavior. In effect, we bridge traditional regularization and deep learning with a principled reconstruction approach.

We consider a large class of shallow neural networks with randomly initialized parameters and rectified linear unit activation functions. We prove that these random neural networks are well-defined non-Gaussian processes. As a by-product, we demonstrate that these networks are solutions to stochastic differential equations driven by impulsive white noise (combinations of random Dirac measures). These processes are parameterized by the law of the weights and biases as well as the density of activation thresholds in each bounded region of the input domain. We prove that these processes are isotropic and wide-sense self-similar with Hurst exponent 3/2. We also derive a remarkably simple closed-form expression for their autocovariance function. Our results are fundamentally different from prior work in that we consider a non-asymptotic viewpoint: The number of neurons in each bounded region of the input domain (i.e., the width) is itself a random variable with a Poisson law with mean proportional to the density parameter. Finally, we show that, under suitable hypotheses, as the expected width tends to infinity, these processes can converge in law not only to Gaussian processes, but also to non-Gaussian processes depending on the law of the weights. Our asymptotic results provide a new take on several classical results (wide networks converge to Gaussian processes) as well as some new ones (wide networks can converge to non-Gaussian processes).

The stiffness of cells and of their nuclei is a biomarker of several pathological conditions. Current measurement methods rely on invasive physical probes that yield one or two stiffness values for the whole cell. However, the internal distribution of cells is heterogeneous. We propose a framework to estimate maps of intracellular and intranuclear stiffness inside deforming cells from fluorescent image sequences. Our scheme requires the resolution of two inverse problems. First, we use a novel optical-flow method that penalizes the nuclear norm of the Hessian to favor deformations that are continuous and piecewise linear, which we show to be compatible with elastic models. We then invert these deformations for the relative intracellular stiffness using a novel system of elliptic PDEs. Our method operates in quasi-static conditions and can still provide relative maps even in the absence of knowledge about the boundary conditions. We compare the accuracy of both methods to the state of the art on simulated data. The application of our method to real data of different cell strains allows us to distinguish different regions inside their nuclei.

We investigate the learning of a convex patch-based synthesis model for the reconstruction of images. In essence, we propose to learn a dictionary via bilevel optimization for denoising. Using implicit differentiation, we find a closed-form formula of the derivative of the minimizer of an objective with respect to the atoms of the dictionary. We also propose a novel way to handle the mean of each patch of the predictions, which improves our model when it is applied to other inverse problems. For minimizing the objective involving the learnt dictionary, we propose an early stopping criterion to further improve the performance of the model for denoising. Finally, we assess our model in a compressed sensing MRI inverse problem and show that, despite being trained on denoising only, our model yields good reconstruction performances.

Lipschitz-constrained neural networks have several advantages over unconstrained ones and can be applied to a variety of problems, making them a topic of attention in the deep learning community. Unfortunately, it has been shown both theoretically and empirically that they perform poorly when equipped with ReLU activation functions. By contrast, neural networks with learnable 1-Lipschitz linear splines are known to be more expressive. In this paper, we show that such networks correspond to global optima of a constrained functional optimization problem that consists of the training of a neural network composed of 1-Lipschitz linear layers and 1-Lipschitz freeform activation functions with second-order total-variation regularization. Further, we propose an efficient method to train these neural networks. Our numerical experiments show that our trained networks compare favorably with existing 1-Lipschitz neural architectures.

In supervised learning, the regularization path is sometimes used as a convenient theoretical proxy for the optimization path of gradient descent initialized from zero. In this paper, we study a modification of the regularization path for infinite-width 2-layer ReLU neural networks with nonzero initial distribution of the weights at different scales. By exploiting a link with unbalanced optimal-transport theory, we show that, despite the non-convexity of the 2-layer network training, this problem admits an infinite-dimensional convex counterpart. We formulate the corresponding functional-optimization problem and investigate its main properties. In particular, we show that, as the scale of the initialization ranges between $0$ and $+\infty$, the associated path interpolates continuously between the so-called kernel and rich regimes. Numerical experiments confirm that, in our setting, the scaling path and the final states of the optimization path behave similarly, even beyond these extreme points.

We formulate an optimization problem to estimate probability densities in the context of multidimensional problems that are sampled with uneven probability. It considers detector sensitivity as an heterogeneous density and takes advantage of the computational speed and flexible boundary conditions offered by splines on a grid. We choose to regularize the Hessian of the spline via the nuclear norm to promote sparsity. As a result, the method is spatially adaptive and stable against the choice of the regularization parameter, which plays the role of the bandwidth. We test our computational pipeline on standard densities and provide software. We also present a new approach to PET rebinning as an application of our framework.

Abstract. We propose to learn non-convex regularizers with a prescribed upper bound on their weak-convexity modulus. Such regularizers give rise to variational denoisers that minimize a convex energy. They rely on few parameters (less than 15,000) and offer a signal-processing interpretation as they mimic handcrafted sparsity-promoting regularizers. Through numerical experiments, we show that such denoisers outperform convex-regularization methods as well as the popular BM3D denoiser. Additionally, the learned regularizer can be deployed to solve inverse problems with iterative schemes that provably converge. For both CT and MRI reconstruction, the regularizer generalizes well and offers an excellent tradeoff between performance, number of parameters, guarantees, and interpretability when compared to other data-driven approaches.

A ridge is a function that is characterized by a one-dimensional profile (activation) and a multidimensional direction vector. Ridges appear in the theory of neural networks as functional descriptors of the effect of a neuron, with the direction vector being encoded in the linear weights. In this paper, we investigate properties of the Radon transform in relation to ridges and to the characterization of neural networks. We introduce a broad category of hyper-spherical Banach subspaces (including the relevant subspace of measures) over which the back-projection operator is invertible. We also give conditions under which the back-projection operator is extendable to the full parent space with its null space being identifiable as a Banach complement. Starting from first principles, we then characterize the sampling functionals that are in the range of the filtered Radon transform. Next, we extend the definition of ridges for any distributional profile and determine their (filtered) Radon transform in full generality. Finally, we apply our formalism to clarify and simplify some of the results and proofs on the optimality of ReLU networks that have appeared in the literature.

The usual explanation of the efficacy of wavelet-based methods hinges on the sparsity of many real-world objects in the wavelet domain. Yet, standard wavelet-shrinkage techniques for sparse reconstruction are not competitive in practice, one reason being that the lack of shift-invariance of the wavelet transform produces blocky artifacts. The standard remedy is cycle spinning, which results in a substantial reduction of these artifacts. In this letter, we propose a theoretical investigation of the sparsity of solutions to the cycle-spinning variant of wavelet-based resolutions of linear inverse problems. We derive a representer theorem that provides a complete characterization of the solution set. Our theorem indicates that the solutions are typically not sparse, where sparsity is measured with respect to the wavelet dictionary. This exposes that the role of sparsity in the success of wavelet-based solutions of linear inverse problems requires further investigation. We corroborate our theoretical results with numerical examples for the problem of image denoising.

The projection of sample measurements onto a reconstruction space represented by a basis on a regular grid is a powerful and simple approach to estimate a probability density function. In this paper, we focus on Riesz bases and propose a projection operator that, in contrast to previous works, guarantees the bona fide properties for the estimate, namely, non-negativity and total probability mass 1. Our bona fide projection is defined as a convex problem. We propose solution techniques and evaluate them. Results suggest an improved performance, specifically in circumstances prone to rippling effects.

Reservoir Computing is a class of Recurrent Neural Networks with internal weights fixed at random. Stability relates to the sensitivity of the network state to perturbations. It is an important property in Reservoir Computing as it directly impacts performance. In practice, it is desirable to stay in a stable regime, where the effect of perturbations does not explode exponentially, but also close to the chaotic frontier where reservoir dynamics are rich. Open questions remain today regarding input regularization and discontinuous activation functions. In this work, we use the recurrent kernel limit to draw new insights on stability in reservoir computing. This limit corresponds to large reservoir sizes, and it already becomes relevant for reservoirs with a few hundred neurons. We obtain a quantitative characterization of the frontier between stability and chaos, which can greatly benefit hyperparameter tuning. In a broader sense, our results contribute to understanding the complex dynamics of Recurrent Neural Networks.

MOTIVATION: Rotated template matching is an efficient and versatile algorithm to analyze microscopy images, as it automates the detection of stereotypical structures, such as organelles that can appear at any orientation. Its performance however quickly degrades in noisy image data. RESULTS: We introduce Steer'n'Detect, an ImageJ plugin implementing a recently published algorithm to detect patterns of interest at any orientation with high accuracy from a single template in 2D images. Steer'n'Detect provides a faster and more robust substitute to template matching. By adapting to the statistics of the image background, it guarantees accurate results even in the presence of noise. The plugin comes with an intuitive user interface facilitating results analysis and further post-processing. AVAILABILITY AND IMPLEMENTATION: https://github.com/Biomedical-Imaging-Group/Steer-n-Detect. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

We formulate as an inverse problem the construction of sparse parametric continuous curve models that fit a sequence of contour points. Our prior is incorporated as a regularization term that encourages rotation invariance and sparsity. We prove that an optimal solution to the inverse problem is a closed curve with spline components. We then show how to efficiently solve the task using B-splines as basis functions. We extend our problem formulation to curves made of two distinct components with complementary smoothness properties and solve it using hybrid splines. We illustrate the performance of our model on contours of different smoothness. Our experimental results show that we can faithfully reconstruct any general contour using few parameters, even in the presence of imprecisions in the measurements.

We present a new family of active surfaces for the semiautomatic segmentation of volumetric objects in 3D biomedical images. We represent our deformable model by a subdivision surface encoded by a small set of control points and generated through a geometric refinement process. The subdivision operator confers important properties to the surface such as smoothness, reproduction of desirable shapes and interpolation of the control points. We deform the subdivision surface through the minimization of suitable gradient-based and region-based energy terms that we have designed for that purpose. In addition, we provide an easy way to combine these energies with convolutional neural networks. Our active subdivision surface satisfies the property of multiresolution, which allows us to adopt a coarse-to-fine optimization strategy. This speeds up the computations and decreases its dependence on initialization compared to singleresolution active surfaces. Performance evaluations on both synthetic and real biomedical data show that our active subdivision surface is robust in the presence of noise and outperforms current state-of-the-art methods. In addition, we provide a software that gives full control over the active subdivision surface via an intuitive manipulation of the control points.

We provide a complete pipeline for the detection of patterns of interest in an image. In our approach, the patterns are assumed to be adequately modeled by a known template, and are located at unknown positions and orientations that we aim at retrieving. We propose a continuous-domain additive image model, where the analyzed image is the sum of the patterns to localize and a background with self-similar isotropic power-spectrum. We are then able to compute the optimal filter fulfilling the SNR criterion based on one single template and background pair: it strongly responds to the template while being optimally decoupled from the background model. In addition, we constrain our filter to be steerable, which allows for a fast template detection together with orientation estimation. In practice, the implementation requires to discretize a continuous-domain formulation on polar grids, which is performed using quadratic radial B-splines. We demonstrate the practical usefulness of our method on a variety of template approximation and pattern detection experiments. We show that the detection performance drastically improves when we exploit the statistics of the background via its power-spectrum decay, which we refer to as spectral-shaping. The proposed scheme outperforms state-of-the-art steerable methods by up to 50% of absolute detection performance.

Quantitative phase imaging (QPI) is an emerging label-free technique that produces images containing morphological and dynamical information without contrast agents. Unfortunately, the phase is wrapped in most imaging system. Phase unwrapping is the computational process that recovers a more informative image. It is particularly challenging with thick and complex samples such as organoids. Recent works that rely on supervised training show that deep learning is a powerful method to unwrap the phase; however, supervised approaches require large and representative datasets which are difficult to obtain for complex biological samples. Inspired by the concept of deep image priors, we propose a deep-learning-based method that does not need any training set. Our framework relies on an untrained convolutional neural network to accurately unwrap the phase while ensuring the consistency of the measurements. We experimentally demonstrate that the proposed method faithfully recovers the phase of complex samples on both real and simulated data. Our work paves the way to reliable phase imaging of thick and complex samples with QPI.

We propose a novel unsupervised deep-learning-based algorithm for dynamic magnetic resonance imaging (MRI) reconstruction. Dynamic MRI requires rapid data acquisition for the study of moving organs such as the heart. We introduce a generalized version of the deep-image-prior approach, which optimizes the weights of a reconstruction network to fit a sequence of sparsely acquired dynamic MRI measurements. Our method needs neither prior training nor additional data. In particular, for cardiac images, it does not require the marking of heartbeats or the reordering of spokes. The key ingredients of our method are threefold: 1) a fixed low-dimensional manifold that encodes the temporal variations of images; 2) a network that maps the manifold into a more expressive latent space; and 3) a convolutional neural network that generates a dynamic series of MRI images from the latent variables and that favors their consistency with the measurements in k -space. Our method outperforms the state-of-the-art methods quantitatively and qualitatively in both retrospective and real fetal cardiac datasets. To the best of our knowledge, this is the first unsupervised deep-learning-based method that can reconstruct the continuous variation of dynamic MRI sequences with high spatial resolution.

We study a variant of the interpolation problem where the continuously defined solution is regularized by minimizing the Lp-norm of its second-order derivative. For this continuous-domain problem, we propose an exact discretization scheme that restricts the search space to quadratic splines with knots on an uniform grid. This leads to a discrete finite-dimensional problem that is computationally tractable. Another benefit of our spline search space is that, when the grid is sufficiently fine, it contains functions that are arbitrarily close to the solutions of the underlying unrestricted problem. We implement an iteratively reweighted algorithm with a grid-refinement strategy that computes the solution within a prescribed accuracy. Finally, we present experimental results that illustrate characteristics, such as sparsity, of the Lp-regularized interpolants.

We propose a novel sparse dictionary learning method for planar shapes in the sense of Kendall, namely configurations of landmarks in the plane considered up to similitudes. Our shape dictionary method provides a good trade-off between algorithmic simplicity and faithfulness with respect to the nonlinear geometric structure of Kendall's shape space. Remarkably, it boils down to a classical dictionary learning formulation modified using complex weights. Existing dictionary learning methods extended to nonlinear spaces map the manifold either to a reproducing kernel Hilbert space or to a tangent space. The first approach is unnecessarily heavy in the case of Kendall's shape space and causes the geometrical understanding of shapes to be lost, while the second one induces distortions and theoretical complexity. Our approach does not suffer from these drawbacks. Instead of embedding the shape space into a linear space, we rely on the hyperplane of centered configurations, including preshapes from which shapes are defined as rotation orbits. In this linear space, the dictionary atoms are scaled and rotated using complex weights before summation. Furthermore, our formulation is more general than Kendall's original one: it applies to discretely defined configurations of landmarks as well as continuously defined interpolating curves. We implemented our algorithm by adapting the method of optimal directions combined to a Cholesky-optimized order recursive matching pursuit. An interesting feature of our shape dictionary is that it produces visually realistic atoms, while guaranteeing reconstruction accuracy. Its efficiency can mostly be attributed to a clear formulation of the framework with complex numbers. We illustrate the strong potential of our approach for the characterization of datasets of shapes up to similitudes and the analysis of patterns in deforming two-dimensional shapes.

Single-particle cryo-electron microscopy (cryo-EM) reconstructs the three-dimensional (3D) structure of biomolecules from a large set of 2D projection images with random and unknown orientations. A crucial step in the single-particle cryo-EM pipeline is 3D refinement, which resolves a highresolution 3D structure from an initial approximate volume by refining the estimation of the orientation of each projection. In this work, we propose a new approach that refines the projection angles on the continuum. We formulate the optimization problem over the density map and the orientations jointly. The density map is updated using the efficient alternating-direction method of multipliers, while the orientations are updated through a semicoordinate- wise gradient descent for which we provide an explicit derivation of the gradient. Our method eliminates the requirement for a fine discretization of the orientation space and does away with the classical but computationally expensive templatematching step. Numerical results demonstrate the feasibility and performance of our approach compared to several baselines.

The motivation for this work is to improve the performance of deep neural networks through the optimization of the individual activation functions. Since the latter results in an infinite- dimensional optimization problem, we resolve the ambiguity by searching for the sparsest and most regular solution in the sense of Lipschitz. To that end, we first introduce a bound that relates the properties of the pointwise nonlinearities to the global Lipschitz constant of the network. By using the proposed bound as regularizer, we then derive a representer theorem that shows that the optimum configuration is achievable by a deep spline network. It is a variant of a conventional deep ReLU network where each activation function is a piecewise-linear spline with adaptive knots. The practical interest is that the underlying spline activations can be expressed as linear combinations of ReLU units and optimized using ℓ1-minimization techniques.

We propose a discretization method for continuous-domain linear inverse problems with multiple-order total-variation (TV) regularization. It is based on a recent result that proves that such inverse problems have sparse polynomial-spline solutions. Our method consists in restricting the search space to splines with knots on a uniform grid, which results in a standard convex finite-dimensional problem. As basis functions for this search space, we use the B-splines matched to the regularization order, which are optimally localized. This leads to a well-conditioned, computationally feasible optimization task. Our proposed iterative multiresolution algorithm then refines the grid size until a desired level of accuracy is met and converges to sparse solutions of our inverse problem. Finally, we present experimental results that validate our approach.

We propose to optimize the activation functions of a deep neural network by adding a corresponding functional regularization to the cost function. We justify the use of a second-order total-variation criterion. This allows us to derive a general representer theorem for deep neural networks that makes a direct connection with splines and sparsity. Specifically, we show that the optimal network configuration can be achieved with activation functions that are nonuniform linear splines with adaptive knots. The bottom line is that the action of each neuron is encoded by a spline whose parameters (including the number of knots) are optimized during the training procedure. The scheme results in a computational structure that is compatible with existing deep-ReLU, parametric ReLU, APL (adaptive piecewise-linear) and MaxOut architectures. It also suggests novel optimization challenges and makes an explicit link with $\ell_1$ minimization and sparsity-promoting techniques.

The detection of landmarks or patterns is of interest for extracting features in biological images. Hence, algorithms for finding these keypoints have been extensively investigated in the literature, and their localization and detection properties are well known. In this paper, we study the complementary topic of local orientation estimation, which has not received similar attention. Simply stated, the problem that we address is the following: estimate the angle of rotation of a pattern with steerable filters centered at the same location, where the image is corrupted by colored isotropic Gaussian noise. For this problem, we propose an estimator formulated as linear combinations of circular harmonics with given radial profiles. We prove that the proposed estimator is unbiased. This property allows us to use a statistical framework based on the Cramér--Rao lower bound (CRLB) to study the limits on the accuracy of the corresponding class of estimators. We aim at evaluating the performance of detection methods based on steerable filters in terms of angular accuracy (as a lower bound), while considering the connection to maximum likelihood estimation. Beyond the general results, we analyze the asymptotic behavior of the lower bound in terms of the order of steerablility and propose an optimal subset of components that minimizes the bound. We define a mechanism for selecting optimal subspaces of the span of the detectors. These are characterized by the most relevant angular frequencies. Finally, we project our template to the span of circular harmonics with given radial profiles and experimentally show that the prediction accuracy achieves the predicted CRLB. As an extension, we also consider steerable wavelet detectors.

We study continuous-domain linear inverse problems with generalized total-variation (gTV) regularization, expressed in terms of a regularization operator L. It has recently been proved that such inverse problems have sparse spline solutions, with fewer jumps than the number of measurements. Moreover, the type of spline solely depends on L (L-splines) and is independent of the measurements. The continuous-domain inverse problem can be recast in an exact way as a finite-dimensional problem by restricting the search space to splines with knots on a uniform finite grid. However, expressing the L-spline coefficients in the dictionary basis of the Green's function of L is ill-suited for practical problems due to its infinite support. Instead, we propose to formulate the problem in the B-spline dictionary basis, which leads to better-conditioned problems. As we make the grid finer, we show that a solution of the continuous-domain problem can be approached arbitrarily closely with functions of this search space. This result motivates our proposed multiresolution algorithm, which computes sparse solutions of our inverse problem. We demonstrate that this algorithm is computationally feasible for 1D signals when L is an ordinary differential operator.

Motivation: We introduce a formulation for the general task of finding diverse shortest paths between two end-points. Our approach is not linked to a specific biological problem and can be applied to a large variety of images thanks to its generic implementation as a user-friendly ImageJ/Fiji plugin. It relies on the introduction of additional layers in a Viterbi path graph, which requires slight modifications to the standard Viterbi algorithm rules. This layered graph construction allows for the specification of various constraints imposing diversity between solutions. Results: The software allows obtaining a collection of diverse shortest paths under some user-defined constraints through a convenient and user-friendly interface. It can be used alone or be integrated into larger image analysis pipelines. Availability and implementation: http://bigwww.epfl.ch/algorithms/diversepathsj. Contact: [email protected] or [email protected]. Supplementary information: Supplementary data are available at Bioinformatics online.

Current algorithmic approaches for piecewise affine motion estimation are based on alternating motion segmentation and estimation. We propose a new method to estimate piecewise affine motion fields directly without intermediate segmentation. To this end, we reformulate the problem by imposing piecewise constancy of the parameter field, and derive a specific proximal splitting optimization scheme. A key component of our framework is an efficient one-dimensional piecewise-affine estimator for vector-valued signals. The first advantage of our approach over segmentation-based methods is its absence of initialization. The second advantage is its lower computational cost which is independent of the complexity of the motion field. In addition to these features, we demonstrate competitive accuracy with other piecewise-parametric methods on standard evaluation benchmarks. Our new regularization scheme also outperforms the more standard use of total variation and total generalized variation.

Localization microscopy, such as STORM / PALM, can reconstruct super-resolution images with a nanometer resolution through the iterative localization of fluorescence molecules. Recent studies in this area have focused mainly on the localization of densely activated molecules to improve temporal resolutions. However, higher density imaging requires an advanced algorithm that can resolve closely spaced molecules. Accordingly, sparsitydriven methods have been studied extensively. One of the major limitations of existing sparsity-driven approaches is the need for a fine sampling grid or for Taylor series approximation which may result in some degree of localization bias toward the grid. In addition, prior knowledge of the point-spread function (PSF) is required. To address these drawbacks, here we propose a true grid-free localization algorithm with adaptive PSF estimation. Specifically, based on the observation that sparsity in the spatial domain implies a low rank in the Fourier domain, the proposed method converts source localization problems into Fourier-domain signal processing problems so that a truly gridfree localization is possible. We verify the performance of the newly proposed method with several numerical simulations and a live-cell imaging experiment.

We provide a generic framework to learn shape dictionaries of landmark-based curves that are defined in the continuous domain. We first present an unbiased alignment method that involves the construction of a mean shape as well as training sets whose elements are subspaces that contain all affine transformations of the training samples. The alignment relies on orthogonal projection operators that have a closed form. We then present algorithms to learn shape dictionaries according to the structure of the data that needs to be encoded: 1) projection-based functional principal-component analysis for homogeneous data and 2) continuous-domain sparse shape encoding to learn dictionaries that contain imbalanced data, outliers, or different types of shape structures. Through parametric spline curves, we provide a detailed and exact implementation of our method. We demonstrate that it requires fewer parameters than purely discrete methods and that it is computationally more efficient and accurate. We illustrate the use of our framework for dictionary learning of structures in biomedical images as well as for shape analysis in bioimaging.

The imaging of proteins within macromolecular complexes has been limited by the low axial resolution of optical microscopes. To overcome this problem, we propose a novel computational reconstruction method that yields isotropic resolution in fluorescence imaging. The guiding principle is to reconstruct a single volume from the observations of multiple rotated particles. Our new operational framework detects particles, estimates their orientation, and reconstructs the final volume. The main challenge comes from the absence of initial template and a priori knowledge about the orientations. We formulate the estimation as a blind inverse problem, and propose a block-coordinate stochastic approach to solve the associated non-convex optimization problem. The reconstruction is performed jointly in multiple channels. We demonstrate that our method is able to reconstruct volumes with 3D isotropic resolution on simulated data. We also perform isotropic reconstructions from real experimental data of doubly labeled purified human centrioles. Our approach revealed the precise localization of the centriolar protein Cep63 around the centriole microtubule barrel. Overall, our method offers new perspectives for applications in biology that require the isotropic mapping of proteins within macromolecular assemblies.

We present a new image reconstruction method that replaces the projector in a projected gradient descent (PGD) with a convolutional neural network (CNN). Recently, CNNs trained as image-to-image regressors have been successfully used to solve inverse problems in imaging. However, unlike existing iterative image reconstruction algorithms, these CNN-based approaches usually lack a feedback mechanism to enforce that the reconstructed image is consistent with the measurements. We propose a relaxed version of PGD wherein gradient descent enforces measurement consistency, while a CNN recursively projects the solution closer to the space of desired reconstruction images. We show that this algorithm is guaranteed to converge and, under certain conditions, converges to a local minimum of a non-convex inverse problem. Finally, we propose a simple scheme to train the CNN to act like a projector. Our experiments on sparse-view computed-tomography reconstruction show an improvement over total variation-based regularization, dictionary learning, and a state-of-the-art deep learning-based direct reconstruction technique.

Learning Tomography (LT) is a nonlinear optimization algorithm for computationally imaging three-dimensional (3D) distribution of the refractive index in semi-transparent samples. Since the energy function in LT is generally non-convex, the solution it obtains is not guaranteed to be globally optimal. In this paper, we describe linear and nonlinear tomographic reconstruction methods and compare them numerically. We present a review of the LT and, in addition, we investigate the influence of the initialization and exemplify the effect of regularization on the convergence of the algorithm. In particular, we show that both are essential for high-quality imaging in strongly scattering scenarios.

We propose a new statistical dictionary learning algorithm for sparse signals that is based on an $\alpha$-stable innovation model. The parameters of the underlying model—that is, the atoms of the dictionary, the sparsity index $\alpha$ and the dispersion of the transform-domain coefficients—are recovered using a new type of probability distribution tomography. Specifically, we drive our estimator with a series of random projections of the data, which results in an efficient algorithm. Moreover, since the projections are achieved using linear combinations, we can invoke the generalized central limit theorem to justify the use of our method for sparse signals that are not necessarily $\alpha$-stable. We evaluate our algorithm by performing two types of experiments: image in-painting and image denoising. In both cases, we find that our approach is competitive with state-of-the-art dictionary learning techniques. Beyond the algorithm itself, two aspects of this study are interesting in their own right. The first is our statistical formulation of the problem, which unifies the topics of dictionary learning and independent component analysis. The second is a generalization of a classical theorem about isometries of $\ell_p$-norms that constitutes the foundation of our approach.

Existing shape models with spherical topology are typically designed either in the discrete domain using interpolating polygon meshes or in the continuous domain using smooth but non-interpolating schemes such as subdivision or NURBS. Both polygon models and subdivision methods require a large number of parameters to model smooth surfaces. NURBS need fewer parameters but have a complicated rational expression and non-uniform shifts in their formulation. We present a new method to construct deformable closed surfaces, which includes exact spheres, by combining the best of two worlds: a smooth, interpolating model with a continuously varying tangent plane and well-defined curvature at every point on the surface. Our formulation is considerably simpler than NURBS and requires fewer parameters than polygon meshes. We demonstrate the generality of our method with applications including intuitive user-interactive shape modeling, continuous surface deformation, shape morphing, reconstruction of shapes from parameterized point clouds, and fast iterative shape optimization for image segmentation. Comparisons with discrete methods and non-interpolating approaches highlight the advantages of our framework.

We present a generic method to construct orthogonal projectors for two-dimensional landmark-based parametric spline curves. We construct vector spaces that define a geometric transformation (e.g., affine, similarity, and scaling) that is applied to a reference curve. These vector spaces contain all parametric curves up to the chosen transformation. We define the vector spaces implicitly through an orthogonal projection operator and present a theorem that characterizes the projector for landmark-based spline curves, which are popular for the user-interactive analysis of biomedical images. Finally, we show how shape priors are constructed with the spline projector and provide an example of application for the segmentation of microscopy images in biology.

We present a new family of snakes that satisfy the property of multiresolution by exploiting subdivision schemes. We show in a generic way how to construct such snakes based on an admissible subdivision mask. We derive the necessary energy formulations and provide the formulas for their efficient computation. Depending on the choice of the mask, such models have the ability to reproduce trigonometric or polynomial curves. They can also be designed to be interpolating, a property that is useful in user-interactive applications. We provide explicit examples of subdivision snakes and illustrate their use for the segmentation of bioimages. We show that they are robust in the presence of noise and provide a multiresolution algorithm to enlarge their basin of attraction, which decreases their dependence on initialization compared to singleresolution snakes. We show the advantages of the proposed model in terms of computation and segmentation of structures with different sizes.

We present texture operators encoding class-specific local organizations of image directions (LOIDs) in a rotation-invariant fashion. The LOIDs are key for visual understanding, and are at the origin of the success of the popular approaches, such as local binary patterns (LBPs) and the scale-invariant feature transform (SIFT). Whereas, LBPs and SIFT yield hand-crafted image representations, we propose to learn data-specific representations of the LOIDs in a rotation-invariant fashion. The image operators are based on steerable circular harmonic wavelets (CHWs), offering a rich and yet compact initial representation for characterizing natural textures. The joint location and orientation required to encode the LOIDs is preserved by using moving frames (MFs) texture representations built from locally-steered image gradients that are invariant to rigid motions. In a second step, we use support vector machines to learn a multi-class shaping matrix for the initial CHW representation, yielding data-driven MFs called steerable wavelet machines (SWMs). The SWM forward function is composed of linear operations (i.e., convolution and weighted combinations) interleaved with non-linear steermax operations. We experimentally demonstrate the effectiveness of the proposed operators for classifying natural textures. Our scheme outperforms recent approaches on several test suites of the Outex and the CUReT databases.

In this paper, we propose a novel deep convolutional neural network (CNN)-based algorithm for solving ill-posed inverse problems. Regularized iterative algorithms have emerged as the standard approach to ill-posed inverse problems in the past few decades. These methods produce excellent results, but can be challenging to deploy in practice due to factors including the high computational cost of the forward and adjoint operators and the difficulty of hyperparameter selection. The starting point of this paper is the observation that unrolled iterative methods have the form of a CNN (filtering followed by pointwise nonlinearity) when the normal operator (H*H, where H* is the adjoint of the forward imaging operator, H) of the forward model is a convolution. Based on this observation, we propose using direct inversion followed by a CNN to solve normal-convolutional inverse problems. The direct inversion encapsulates the physical model of the system, but leads to artifacts when the problem is ill posed; the CNN combines multiresolution decomposition and residual learning in order to learn to remove these artifacts while preserving image structure. We demonstrate the performance of the proposed network in sparse-view reconstruction (down to 50 views) on parallel beam X-ray computed tomography in synthetic phantoms as well as in real experimental sinograms. The proposed network outperforms total variation-regularized iterative reconstruction for the more realistic phantoms and requires less than a second to reconstruct a 512 × 512 image on the GPU.

We propose a new, cost-efficient method for computing back projections in parallel-ray X-ray CT. Forward and back projections are the basis of almost all X-ray CT reconstruction methods, but computing these accurately is costly. In the special case of parallel-ray geometry, it turns out that reconstruction requires back projection only. One approach to accelerate the back projection is through interpolation: fit a continuous representation to samples of the desired signal, then sample it at the required locations. Instead, we propose applying a prefilter that has the effect of orthogonally projecting the underlying signal onto the space spanned by the interpolator, which can significantly improve the quality of the interpolation. We then build on this idea by using oblique projection, which simplifies the computation while giving effectively the same improvement in quality. Our experiments on analytical phantoms show that this refinement can improve the reconstruction quality for both filtered back projection and iterative reconstruction in the high-quality regime, i.e., with low noise and many measurements.

We develop a fast algorithm for segmenting 3D images from linear measurements based on the Potts model (or piecewise constant Mumford-Shah model). To that end, we first derive suitable space discretizations of the 3D Potts model, which are capable of dealing with 3D images defined on non-cubic grids. Our discretization allows us to utilize a specific splitting approach, which results in decoupled subproblems of moderate size. The crucial point in the 3D setup is that the number of independent subproblems is so large that we can reasonably exploit the parallel processing capabilities of the graphics processing units (GPUs). Our GPU implementation is up to 18 times faster than the sequential CPU version. This allows to process even large volumes in acceptable runtimes. As a further contribution, we extend the algorithm in order to deal with non-negativity constraints. We demonstrate the efficiency of our method for combined image deconvolution and segmentation on simulated data and on real 3D wide field fluorescence microscopy data.

We investigate the performance of wavelet shrinkage methods for the denoising of symmetric-α-stable (SαS) self-similar stochastic processes corrupted by additive white Gaussian noise (AWGN), where α is tied to the sparsity of the process. The wavelet transform is assumed to be orthonormal and the shrinkage function minimizes the mean-square approximation error (MMSE estimator). We derive the corresponding formula for the expected value of the averaged estimation error. We show that the predicted MMSE is a monotone function of a simple criterion that depends on the wavelet and the statistical parameters of the process. Using the calculus of variations, we then optimize this criterion to find the best performing wavelet within the extended family of Meyer wavelets, which are bandlimited. These are compared with the Daubechies wavelets, which are compactly supported in time. We find that the wavelets that are shorter in time (in particular, the Haar basis) are better suited to denoise the sparser processes (say, α1.6, the limit corresponding to the Gaussian case (fBm) with α = 2.

We propose two minimum-mean-square-error (MMSE) estimation methods for denoising non-Gaussian first-order autoregressive (AR(1)) processes. The first one is based on the message passing framework and gives the exact theoretic MMSE estimator. The second is an iterative algorithm that combines standard wavelet-based thresholding with an optimized non-linearity and cycle-spinning. This method is more computationally efficient than the former and appears to provide the same optimal denoising results in practice. We illustrate the superior performance of both methods through numerical simulations by comparing them with other well-known denoising schemes.

Biomedical image segmentation is an active field of research where deformable models have proved to be efficient. The geometric representation of such models determines their ability to approximate the shape of interest as well as the speed of convergence of related optimization algorithms. We present a new tensor-product parameterization of surfaces that offers the possibility of local refinement. The goal is to allocate additional degrees of freedom to the surface only where an increase in local detail is required. We introduce the possibility of locally increasing the number of control points by inserting basis functions at specific locations. Our approach is generic and relies on refinable functions which satisfy the refinement relation. We show that the proposed method improves brain segmentation in 3D MRI images.

MOTIVATION: SpotCaliper is a novel wavelet-based image-analysis software providing a fast automatic detection scheme for circular patterns (spots), combined with the precise estimation of their size. It is implemented as an ImageJ plugin with a friendly user interface. The user is allowed to edit the results by modifying the measurements (in a semi-automated way), extract data for further analysis. The fine tuning of the detections includes the possibility of adjusting or removing the original detections, as well as adding further spots. RESULTS: The main advantage of the software is its ability to capture the size of spots in a fast and accurate way. AVAILABILITY AND IMPLEMENTATION: http://bigwww.epfl.ch/algorithms/spotcaliper/ CONTACT: [email protected] SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

In analogy with steerable wavelets, we present a general construction of adaptable tight wavelet frames, with an emphasis on scaling operations. In particular, the derived wavelets can be “dilated” by a procedure comparable to the operation of steering steerable wavelets. The fundamental aspects of the construction are the same: an admissible collection of Fourier multipliers is used to extend a tight wavelet frame, and the “scale” of the wavelets is adapted by scaling the multipliers. As an application, the proposed wavelets can be used to improve the frequency localization. Importantly, the localized frequency bands specified by this construction can be scaled efficiently using matrix multiplication.

Our motivation is the design of efficient algorithms to process closed curves represented by basis functions or wavelets. To that end, we introduce an inner-product calculus to evaluate correlations and L2distances between such curves. In particular, we present formulas for the direct and exact evaluation of correlation matrices in the case of closed (i.e., periodic) parametric curves and periodic signals. We give simplifications for practical cases that involve B-splines. To illustrate this approach, we also propose a least-squares approximation scheme that is able to resample curves while minimizing aliasing artifacts. Another application is the exact calculation of the enclosed area.

We introduce a variational phase retrieval algorithm for the imaging of transparent objects. Our formalism is based on the transport-of-intensity equation (TIE), which relates the phase of an optical field to the variation of its intensity along the direction of propagation. TIE practically requires one to record a set of defocus images to measure the variation of intensity. We first investigate the effect of the defocus distance on the retrieved phase map. Based on our analysis, we propose a weighted phase reconstruction algorithm yielding a phase map that minimizes a convex functional. The method is nonlinear and combines different ranges of spatial frequencies - depending on the defocus value of the measurements - in a regularized fashion. The minimization task is solved iteratively via the alternating-direction method of multipliers. Our simulations outperform commonly used linear and nonlinear TIE solvers. We also illustrate and validate our method on real microscopy data of HeLa cells.

A crucial component of steerable wavelets is the radial profile of the generating function in the frequency domain. In this paper, we present an infinite-dimensional optimization scheme that helps us find the optimal profile for a given criterion over the space of tight frames. We consider two classes of criteria that measure the localization of the wavelet. The first class specifies the spatial localization of the wavelet profile, and the second that of the resulting wavelet coefficients. From these metrics and the proposed algorithm, we construct tight wavelet frames that are optimally localized and provide their analytical expression. In particular, one of the considered criterion helps us finding back the popular Simoncelli wavelet profile. Finally, the investigation of local orientation estimation, image reconstruction from detected contours in the wavelet domain, and denoising indicate that optimizing wavelet localization improves the performance of steerable wavelets, since our new wavelets outperform the traditional ones.

We propose a framework for the detection of junctions in images. Although the detection of edges and key points is a well examined and described area, the multiscale detection of junction centers, especially for odd orders, poses a challenge in pattern analysis. The goal of this paper is to build optimal junction detectors based on 2D steerable wavelets that are polar-separable in the Fourier domain. The approaches we develop are general and can be used for the detection of arbitrary symmetric and asymmetric junctions. The backbone of our construction is a multiscale pyramid with a radial wavelet function where the directional components are represented by circular harmonics and encoded in a shaping matrix. We are able to detect M -fold junctions in different scales and orientations. We provide experimental results on both simulated and real data to demonstrate the effectiveness of the algorithm.

We introduce a new model of parametric contours defined in a continuous fashion. Our curve model relies on Hermite spline interpolation and can easily generate curves with sharp discontinuities; it also grants direct access to the tangent at each location. With these two features, the Hermite snake distinguishes itself from classical spline-snake models and allows one to address certain bioimaging problems in a more efficient way. More precisely, the Hermite snake construction allows introducing sharp corners in the snake curve and designing directional energy functionals relying on local orientation information in the input image. Using the formalism of spline theory, the model is shown to meet practical requirements such as invariance to affine transformations and good approximation properties. Finally, the dependence on initial conditions and the robustness to the noise is studied on synthetic data in order to validate our Hermite snake model, and its usefulness is illustrated on real biological images acquired using brightfield, phase-contrast, differential-interference-contrast, and scanning-electron microscopy.

We present a theoretical analysis and comparison of the effect of l1versus l2regularization for the resolution of ill-posed linear inverse and/or compressed sensing problems. Our formulation covers the most general setting where the solution is specified as the minimizer of a convex cost functional. We derive a series of representer theorems that give the generic form of the solution depending on the type of regularization. We start with the analysis of the problem in finite dimensions and then extend our results to the infinite-dimensional spaces l2(Z) and l1(Z). We also consider the use of linear transformations in the form of dictionaries or regularization operators. In particular, we show that the l2solution is forced to live in a predefined subspace that is intrinsically smooth and tied to the measurement operator. The l1solution, on the other hand, is formed by adaptively selecting a subset of atoms in a dictionary that is specified by the regularization operator. Beside the proof that l1solutions are intrinsically sparse, the main outcome of our investigation is that the use of l1regularization is much more favorable for injecting prior knowledge: it results in a functional form that is independent of the system matrix, while this is not so in the l2scenario.

Shape segmentation is an active field of research in biomedical imaging. In this context, we present a new parameterization of a snake that is locally refinable. We introduce the possibility of locally increasing the approximation power of the parametric model by inserting basis functions at a specific location. This is controlled by a user-interface that permits the refinement of an initial segmentation around an anchor position selected by a user. Our approach relies on scaling functions that satisfy the refinement relation and are related to wavelets. We also derive explicit formulas for the energy functions associated to our new parameterization. We demonstrate the accuracy of our snake and its robustness under noisy conditions on phantom data. We also present segmentation results on real cell images, which are our main target. The algorithm is made freely available as a plugin for the open source platform Icy.

We propose a new parametric 3D snake with cylindrical topology. Its construction is based on interpolatory basis functions which facilitates user-interaction because the control points of the snake directly lie on the surface of the deformable cylinder. We prove that the basis functions exactly reproduce a cylinder and propose a new parametrization as a tensor-product spline surface. We provide explicit formulas for the energy function based on Green's theorem that speed up the computation of the optimization algorithm. We have implemented the proposed framework as a freely available open-source plugin for the bioimaging platform Icy. Its utility has been tested on phantom data as well as on real 3D data to segment the spinal cord and the descending aorta.

It is well documented that natural images are compressible in wavelet bases and tend to exhibit fractal properties. In this paper, we investigate statistical models that mimic these behaviors. We then use our models to make predictions on the statistics of the wavelet coefficients. Following an innovation modeling approach, we identify a general class of finite-variance self-similar sparse random processes. We first prove that spatially dilated versions of self-similar sparse processes are asymptotically Gaussian as the dilation factor increases. Based on this fundamental result, we show that the coarse-scale wavelet coefficients of these processes are also asymptotically Gaussian, provided the wavelet has enough vanishing moments. Moreover, we quantify the degree of Gaussianity by deriving the theoretical evolution of the kurtosis of the wavelet coefficients across scales. Finally, we apply our analysis to one- and two-dimensional signals, including natural images, and show that the wavelet coefficients tend to become Gaussian at coarse scales.

To address the classic interior tomography problem where projections at each view extend only to the shadow of a circular region completely interior to the subject being scanned, previously we showed that the exact recovery of two- and three-dimensional piecewise smooth images is guaranteed using a one-dimensional generalized total variation seminorm penalty which allows a much faster reconstruction. To further accelerate the algorithm up to a level for clinical use, this paper proposes a novel multiscale reconstruction method by exploiting the Bedrosian identity of the Hilbert transform. More specifically, we show that the high frequency parts of the one-dimensional signals can be quickly recovered analytically with the Hilbert transform because of the Bedrosian identity. This implies that computationally expensive iterative reconstruction need only be applied to low resolution images in the downsampled domain, which significantly reduces the computational burden. Moreover, even for incomplete trajectories such as circular cone-beam geometry, we demonstrate that the proposed multiscale interior tomography approach can be combined with a novel spectral blending method in order to mitigate cone-beam artifacts from missing frequency regions. We show the efficacy of the proposed multiscale algorithm using circular fan-beam, helical cone-beam data, and circular cone-beam geometry. With a graphics processing unit implementation, we demonstrate that the speed of the algorithm can be significantly accelerated up to the level for clinical use for various acquisition geometries.

We introduce a novel generic energy functional that we employ to solve inverse imaging problems within a variational framework. The proposed regularization family, termed as structure tensor total variation (STV), penalizes the eigenvalues of the structure tensor and is suitable for both grayscale and vector-valued images. It generalizes several existing variational penalties, including the total variation seminorm and vectorial extensions of it. Meanwhile, thanks to the structure tensor's ability to capture first-order information around a local neighborhood, the STV functionals can provide more robust measures of image variation. Further, we prove that the STV regularizers are convex while they also satisfy several invariance properties w.r.t. image transformations. These properties qualify them as ideal candidates for imaging applications. In addition, for the discrete version of the STV functionals we derive an equivalent definition that is based on the patch-based Jacobian operator, a novel linear operator which extends the Jacobian matrix. This alternative definition allow us to derive a dual problem formulation. The duality of the problem paves the way for employing robust tools from convex optimization and enables us to design an efficient and parallelizable optimization algorithm. Finally, we present extensive experiments on various inverse imaging problems, where we compare our regularizers with other competing regularization approaches. Our results are shown to be systematically superior, both quantitatively and visually.

In this paper, we propose a continuous-domain version of principal-component analysis, with the constraint that the underlying family of templates appears at arbitrary orientations. We show that the corresponding principal components are steerable. Our method can be used for designing steerable filters so that they best approximate a given collection of reference templates. We apply this framework to the detection and classification of micrometer-sized particles that are used in a microfluidic diagnostics system. This is done in two steps. First, we decompose the particles into a small number of templates and compute their steerable principal components. Then we use these principal components to automatically estimate the orientation and the class of each particle.

Motivated by the interior tomography problem, we propose a method for exact reconstruction of a region of interest of a function from its local Radon transform in any number of dimensions. Our aim is to verify the feasibility of a one-dimensional reconstruction procedure that can provide the foundation for an efficient algorithm. For a broad class of functions, including piecewise polynomials and generalized splines, we prove that an exact reconstruction is possible by minimizing a generalized total variation seminorm along lines. The main difference with previous works is that our approach is inherently one-dimensional and that it imposes less constraints on the class of admissible signals. Within this formulation, we derive unique reconstruction results using properties of the Hilbert transform, and we present numerical examples of the reconstruction.

We propose a new variational framework for the problem of reconstructing flow fields from noisy measurements. The formalism is based on regularizers penalizing the singular values of the Jacobian of the field. Specifically, we rely on the nuclear norm. Our method is invariant with respect to fundamental transformations and can be efficiently solved. We conduct numerical experiments on several phantom data and report improved performance compared to existing vectorial extensions of total variation and curl-divergence regularizations. Finally, we apply our reconstruction method to an experimentally-acquired phase-contrast MRI recording for enhancing the data visualization.

We propose an efficient construction of wavelet frames in any number of dimensions that are divergence-free. Our approach is based on applying the Leray projector, which is scale-invariant, to a standard wavelet frame. We prove that the projected wavelets retain the basic characteristics (decay rate and order of vanishing moments) of the initial wavelets. Since the Leray projector is also shift-invariant, it is defined as a Fourier multiplier, and our construction is implemented efficiently using the fast Fourier transform. In order to illustrate the practicality of the method, we present vector field denoising experiments.

We present a new trigonometric basis function that is capable of perfectly reproducing circles, spheres and ellipsoids while at the same time being interpolatory. Such basis functions have the advantage that they allow to construct shapes through a sequence of control points that lie on their contour (2-D) or surface (3-D) which facilitates user-interaction, especially in 3-D. Our piecewise exponential basis function has finite support, which enables local control for shape modification. We derive and prove all the necessary properties of the kernel to represent shapes that can be smoothly deformed and show how idealized shapes such as ellipses and spheres can be constructed.

In this letter, we aim to identify the optimal isotropic mother wavelet for a given spatial dimension based on a localization criterion. Within the framework of the calculus of variations, we specify an Euler-Lagrange equation for this problem, and we find the unique analytic solutions. In the one- and two-dimensional cases, the derived wavelets are well known.

Parametric active contours are an attractive approach for image segmentation, thanks to their computational efficiency. They are driven by application-dependent energies that reflect the prior knowledge on the object to be segmented. We propose an energy involving shape priors acting in a regularization-like manner. Thereby, the shape of the snake is orthogonally projected onto the space that spans the affine transformations of a given shape prior. The formulation of the curves is continuous, which provides computational benefits when compared with landmark-based (discrete) methods. We show that this approach improves the robustness and quality of spline-based segmentation algorithms, while its computational overhead is negligible. An interactive and ready-to-use implementation of the proposed algorithm is available and was successfully tested on real data in order to segment Drosophila flies and yeast cells in microscopic images.

Kaiser-Bessel window functions are frequently used to discretize tomographic problems because they have two desirable properties: 1) their short support leads to a low computational cost and 2) their rotational symmetry makes their imaging transform independent of the direction. In this paper, we aim at optimizing the parameters of these basis functions. We present a formalism based on the theory of approximation and point out the importance of the partition-of-unity condition. While we prove that, for compact-support functions, this condition is incompatible with isotropy, we show that minimizing the deviation from the partition of unity condition is highly beneficial. The numerical results confirm that the proposed tuning of the Kaiser-Bessel window functions yields the best performance.

Our goal is to detect and group different kinds of local symmetries in images in a scale- and rotation-invariant way. We propose an efficient wavelet-based method to determine the order of local symmetry at each location. Our algorithm relies on circular harmonic wavelets which are used to generate steerable wavelet channels corresponding to different symmetry orders. To give a measure of local symmetry, we use the F-test to examine the distribution of the energy across different channels. We provide experimental results on synthetic images, biological micrographs, and electron-microscopy images to demonstrate the performance of the algorithm.

We present a new exponential B-spline basis that enables the construction of active contours for the analysis of biomedical images. Our functions generalize the well-known polynomial Hermite B-splines and provide us with a direct control over the tangents of the parameterized contour, which is absent in traditional spline-based active contours. Our basis functions have been designed to perfectly reproduce elliptical and circular shapes. Moreover, they can approximate any closed curve up to arbitrary precision by increasing the number of anchor points. They are therefore well-suited to the segmentation of the roundish objects that are commonly encountered in the analysis of bioimages. We illustrate the performance of an active contour built using our functions on some examples of real biological data.

We present a variational reconstruction algorithm for the phase-retrieval problem by using the differential interference contrast microscopy. Principally, we rely on the transport-of-intensity equation that specifies the sought phase as the solution of a partial differential equation. Our approach is based on an iterative reconstruction algorithm involving the total variation regularisation which is efficiently solved via the alternating direction method of multipliers. We illustrate the applicability of the method via real data experiments. To the best of our knowledge, this work demonstrates the performance of such an iterative algorithm on real data for the first time.

We study the statistics of wavelet coefficients of non-Gaussian images, focusing mainly on the behaviour at coarse scales. We assume that an image can be whitened by a fractional Laplacian operator, which is consistent with an ∥ω∥-γspectral decay. In other words, we model images as sparse and self-similar stochastic processes within the framework of generalised innovation models. We show that the wavelet coefficients at coarse scales are asymptotically Gaussian even if the prior model for fine scales is sparse. We further refine our analysis by deriving the theoretical evolution of the cumulants of wavelet coefficients across scales. Especially, the evolution of the kurtosis supplies a theoretical prediction for the Gaussianity level at each scale. Finally, we provide simulations and experiments that support our theoretical predictions.

In this work, we describe our approach of combining the most effective ideas and tools developed during the past years to build a variational 3D deconvolution system that can be successfully employed in fluorescence microscopy. In particular, the main components of our deconvolution system involve proper handling of image boundaries, choice of a regularizer that is best suited to biological images, and use of an optimization algorithm that can be efficiently implemented on graphics processing units (GPUs) and fully benefit from their massive parallel computational capabilities. We show that our system leads to very competitive results and reduces the computational time by at least one order of magnitude compared to a CPU implementation. This makes the use of advanced deconvolution techniques feasible in practice and attractive computationally.

We study the issue of localization in the context of isotropic wavelet frames. We define a variance-type measure of localization and propose an algorithm based on calculus of variations to minimize this criterion under the constraint of a tight wavelet frame. Based on these calculations, we design the variance-optimal wavelet (VOW). Finally, we demonstrate the advantage of better localization in a practical image-processing task.

We present a method for the unsupervised segmentation of textured images using Potts functionals, which are a piecewise-constant variant of the Mumford and Shah functionals. We propose a minimization strategy based on the alternating direction method of multipliers and dynamic programming. The strategy allows us to process large feature spaces because the computational cost grows only linearly in the feature dimension. In particular, our algorithm has more favorable computational costs for high-dimensional data than graph cuts. Our feature vectors are based on monogenic curvelets. They incorporate multiple resolutions and directional information. The advantage over classical curvelets is that they yield smoother amplitudes due to the envelope effect of the monogenic signal.

Cycle spinning is a widely used approach for improving the performance of wavelet-based methods that solve linear inverse problems. Extensive numerical experiments have shown that it significantly improves the quality of the recovered signal without increasing the computational cost. In this letter, we provide the first theoretical convergence result for cycle spinning for solving general linear inverse problems. We prove that the sequence of reconstructed signals is guaranteed to converge to the minimizer of some global cost function that incorporates all wavelet shifts.

We introduce an adaptive continuous-domain modeling approach to texture and natural images. The continuous-domain image is assumed to be a smooth function, and we embed it in a parameterized Sobolev space. We point out a link between Sobolev spaces and stochastic auto-regressive models, and exploit it for optimally choosing Sobolev parameters from available pixel values. To this aim, we use exact continuous-to-discrete mapping of the auto-regressive model that is based on symmetric exponential splines. The mapping is computationally efficient, and we exploit it for maximizing an approximated Gaussian likelihood function.We account for non-Gaussian Lévy-type processes by deriving a more robust estimator that is based on the sample auto-correlation sequence. Both estimators use multiple initialization values for overcoming the local minima structure of the fitting criteria. Experimental image resizing results indicate that the auto-correlation criterion can cope better with non-Gaussian processes and model mismatch. Our work demonstrates the importance of the auto-correlation function in adaptive image interpolation and image modeling tasks, and we believe it is instrumental in other image processing tasks as well.

We adopt an innovation-driven framework and investigate the sparse/compressible distributions obtained by linearly measuring or expanding continuous-domain stochastic models. Starting from the first principles, we show that all such distributions are necessarily infinitely divisible. This property is satisfied by many distributions used in statistical learning, such as Gaussian, Laplace, and a wide range of fat-tailed distributions, such as student's-t and α-stable laws. However, it excludes some popular distributions used in compressed sensing, such as the Bernoulli-Gaussian distribution and distributions, that decay like exp (-O(|x|p)) for 1p<; 2. We further explore the implications of infinite divisibility on distributions and conclude that tail decay and unimodality are preserved by all linear functionals of the same continuous-domain process. We explain how these results help in distinguishing suitable variational techniques for statistically solving inverse problems like denoising.

We consider the estimation of an independent and identically distributed (i.i.d.) (possibly non-Gaussian) vector x ∈ Rnfrom measurements y ∈ Rmobtained by a general cascade model consisting of a known linear transform followed by a probabilistic componentwise (possibly nonlinear) measurement channel. A novel method, called adaptive generalized approximate message passing (adaptive GAMP) is presented. It enables the joint learning of the statistics of the prior and measurement channel along with estimation of the unknown vector x. We prove that, for large i.i.d. Gaussian transform matrices, the asymptotic componentwise behavior of the adaptive GAMP is predicted by a simple set of scalar state evolution equations. In addition, we show that the adaptive GAMP yields asymptotically consistent parameter estimates, when a certain maximum-likelihood estimation can be performed in each step. This implies that the algorithm achieves a reconstruction quality equivalent to the oracle algorithm that knows the correct parameter values. Remarkably, this result applies to essentially arbitrary parametrizations of the unknown distributions, including nonlinear and non-Gaussian ones. The adaptive GAMP methodology thus provides a systematic, general and computationally efficient method applicable to a large range of linear-nonlinear models with provable guarantees.

This paper is devoted to the characterization of an extended family of continuous-time autoregressive moving average (CARMA) processes that are solutions of stochastic differential equations driven by white Lévy innovations. These are completely specified by: 1) a set of poles and zeros that fixes their correlation structure and 2) a canonical infinitely divisible probability distribution that controls their degree of sparsity (with the Gaussian model corresponding to the least sparse scenario). The generalized CARMA processes are either stationary or nonstationary, depending on the location of the poles in the complex plane. The most basic nonstationary representatives (with a single pole at the origin) are the Lévy processes, which are the non-Gaussian counterparts of Brownian motion. We focus on the general analog-to-discrete conversion problem and introduce a novel spline-based formalism that greatly simplifies the derivation of the correlation properties and joint probability distributions of the discrete versions of these processes. We also rely on the concept of generalized increment process, which suppresses all long range dependencies, to specify an equivalent discrete-domain innovation model. A crucial ingredient is the existence of a minimally supported function associated with the whitening operator L; this B-spline, which is fundamental to our formulation, appears in most of our formulas, both at the level of the correlation and the characteristic function. We make use of these discrete-domain results to numerically generate illustrative examples of sparse signals that are consistent with the continuous-domain model.

We introduce a general distributional framework that results in a unifying description and characterization of a rich variety of continuous-time stochastic processes. The cornerstone of our approach is an innovation model that is driven by some generalized white noise process, which may be Gaussian or not (e.g., Laplace, impulsive Poisson, or alpha stable). This allows for a conceptual decoupling between the correlation properties of the process, which are imposed by the whitening operator L, and its sparsity pattern, which is determined by the type of noise excitation. The latter is fully specified by a Lévy measure. We show that the range of admissible innovation behavior varies between the purely Gaussian and super-sparse extremes. We prove that the corresponding generalized stochastic processes are well-defined mathematically provided that the (adjoint) inverse of the whitening operator satisfies some Lp bound for p ≥ 1. We present a novel operator-based method that yields an explicit characterization of all Lévy-driven processes that are solutions of constant-coefficient stochastic differential equations. When the underlying system is stable, we recover the family of stationary continuous-time autoregressive moving average processes (CARMA), including the Gaussian ones. The approach remains valid when the system is unstable and leads to the identification of potentially useful generalizations of the Lévy processes, which are sparse and non-stationary. Finally, we show that these processes admit a sparse representation in some matched wavelet domain and provide a full characterization of their transform-domain statistics.

In this paper, we reconstruct signals from underdetermined linear measurements where the componentwise gains of the measurement system are unknown a priori. The reconstruction is performed through an adaptation of the messagepassing algorithm called adaptive GAMP that enables joint gain calibration and signal estimation. To evaluate our approach, we apply it to the problem of sparse recovery and compare it against an ℓ1-based approach. We numerically show that adaptive GAMP yields excellent results even for a moderate amount of data. It approaches the performance of oracle GAMP where the gains are perfectly known asymptotically.

Sinusoidal transforms such as the DCT are known to be optimal-that is, asymptotically equivalent to the Karhunen-Loève transform (KLT)-for the representation of Gaussian stationary processes, including the classical AR(1) processes. While the KLT remains applicable for non-Gaussian signals, it loses optimality and, is outperformed by the independent-component analysis (ICA), which aims at producing the most-decoupled representation. In this paper, we consider an extension of the classical AR(1) model that is driven by symmetric-alpha-stable (SαS) noise which is either Gaussian (α = 2) or sparse (0 <; α <; 2). For the sparse (non-Gaussian) regime, we prove that an expansion in a proper wavelet basis (including the Haar transform) is much closer to the optimal orthogonal ICA solution than the classical Fourier-type representations. Our criterion for optimality, which favors independence, is the Kullback-Leibler divergence between the joint pdf of the original signal and the product of the marginals in the transformed domain. We also observe that, for very sparse AR(1) processes (α ≤ 1), the operator-like wavelet transform is indistinguishable from the ICA solution that is determined through numerical optimization.

The steerable wavelet transform is a redundant image representation with the remarkable property that its basis functions can be adaptively rotated to a desired orientation. This makes the transform well-suited to the design of wavelet-based algorithms applicable to images with a high amount of directional features. However, arbitrary modification of the wavelet-domain coefficients may violate consistency constraints because a legitimate representation must be redundant. In this paper, by honoring the redundancy of the coefficients, we demonstrate that it is possible to improve the performance of regularized least-squares problems in the steerable wavelet domain. We illustrate that our consistent method significantly improves upon the performance of conventional denoising with steerable wavelets.

We propose a semi-automatic technique to segment corpus callosum (CC) using a two-stage snake formulation: A restricted affine transform (RAT) constrained snake followed by an unconstrained snake in an iterative fashion. A statistical model is developed to capture the shape variations of CC from a training set, which restrict the unconstrained snake to lie in the shape-space of CC. The geometry of the constrained snake is optimized using a local contrast-based energy over RAT space (which allows for five degrees of freedom). On the other hand, the unconstrained snake is optimized using a unified energy (region, gradient, and curvature energy) formulation. Joint optimization resulted in increased robustness to initialization as well as fast and accurate segmentation. The technique was validated on 243 images taken from the OASIS database and performance was quantified using Jaccard's distance, sensitivity, and specificity as the metrics.

We introduce a complete parameterization of the family of two-dimensional steerable wavelets that are polar-separable in the Fourier domain under the constraint of self-reversibility. These wavelets are constructed by multiorder generalized Riesz transformation of a primary isotropic bandpass pyramid. The backbone of the transform (pyramid) is characterized by a radial frequency profile function $h(\omega)$, while the directional wavelet components at each scale are encoded by an $M \times (2N+1)$ shaping matrix ${\bf U}$, where $M$ is the number of wavelet channels and $N$ the order of the Riesz transform. We provide general conditions on $h(\omega)$ and ${\bf U}$ for the underlying wavelet system to form a tight frame of $L_2(\mathbb{R}^2)$ (with a redundancy factor $4/3M$). The proposed framework ensures that the wavelets are steerable and provides new degrees of freedom (shaping matrix ${\bf U}$) that can be exploited for designing specific wavelet systems. It encompasses many known transforms as particular cases: Simoncelli's steerable pyramid, Marr gradient and Hessian wavelets, monogenic wavelets, and $N$th-order Riesz and circular harmonic wavelets. We take advantage of the framework to construct new generalized spheroidal prolate wavelets, whose angular selectivity is maximized, as well as signal-adapted detectors based on principal component analysis. We also introduce a curvelet-like steerable wavelet system. Finally, we illustrate the advantages of some of the designs for signal denoising, feature extraction, pattern analysis, and source separation.

The Riesz transform is a natural multidimensional extension of the Hilbert transform, and it has been the object of study for many years due to its nice mathematical properties. More recently, the Riesz transform and its variants have been used to construct complex wavelets and steerable wavelet frames in higher dimensions. The flip side of this approach, however, is that the Riesz transform of a wavelet often has slow decay. One can nevertheless overcome this problem by requiring the original wavelet to have sufficient smoothness, decay, and vanishing moments. In this paper, we derive necessary conditions in terms of these three properties that guarantee the decay of the Riesz transform and its variants, and, as an application, we show how the decay of the popular Simoncelli wavelets can be improved by appropriately modifying their Fourier transforms. By applying the Riesz transform to these new wavelets, we obtain steerable frames with rapid decay.

We introduce a new approach for the implementation of minimum mean-square error (MMSE) denoising for signals with decoupled derivatives. Our method casts the problem as a penalized least-squares regression in the redundant wavelet domain. It exploits the link between the discrete gradient and Haar-wavelet shrinkage with cycle spinning. The redundancy of the representation implies that some wavelet-domain estimates are inconsistent with the underlying signal model. However, by imposing additional constraints, our method finds wavelet-domain solutions that are mutually consistent. We confirm the MMSE performance of our method through statistical estimation of Lévy processes that have sparse derivatives.

We present a novel statistically-based discretization paradigm and derive a class of maximum a posteriori (MAP) estimators for solving ill-conditioned linear inverse problems. We are guided by the theory of sparse stochastic processes, which specifies continuous-domain signals as solutions of linear stochastic differential equations. Accordingly, we show that the class of admissible priors for the discretized version of the signal is confined to the family of infinitely divisible distributions. Our estimators not only cover the well-studied methods of Tikhonov and l1-type regularizations as particular cases, but also open the door to a broader class of sparsity-promoting regularization schemes that are typically nonconvex. We provide an algorithm that handles the corresponding nonconvex problems and illustrate the use of our formalism by applying it to deconvolution, magnetic resonance imaging, and X-ray tomographic reconstruction problems. Finally, we compare the performance of estimators associated with models of increasing sparsity.

Practical image-acquisition systems are often modeled as a continuous-domain prefilter followed by an ideal sampler, where generalized samples are obtained after convolution with the impulse response of the device. In this paper, our goal is to interpolate images from a given subset of such samples. We express our solution in the continuous domain, considering consistent resampling as a data-fidelity constraint. To make the problem well posed and ensure edge-preserving solutions, we develop an efficient anisotropic regularization approach that is based on an improved version of the edge-enhancing anisotropic diffusion equation. Following variational principles, our reconstruction algorithm minimizes successive quadratic cost functionals. To ensure fast convergence, we solve the corresponding sequence of linear problems by using multigrid iterations that are specifically tailored to their sparse structure. We conduct illustrative experiments and discuss the potential of our approach both in terms of algorithmic design and reconstruction quality. In particular, we present results that use as little as 2% of the image samples.

Compressed sensing can substantially reduce the number of samples required for conventional signal acquisition at the expense of an additional reconstruction procedure. It also provides robust reconstruction when using quantized measurements, including in the one-bit setting. In this paper, our goal is to design a framework for binary compressed sensing that is adapted to images. Accordingly, we propose an acquisition and reconstruction approach that complies with the high dimensionality of image data and that provides reconstructions of satisfactory visual quality. Our forward model describes data acquisition and follows physical principles. It entails a series of random convolutions performed optically followed by sampling and binary thresholding. The binary samples that are obtained can be either measured or ignored according to predefined functions. Based on these measurements, we then express our reconstruction problem as the minimization of a compound convex cost that enforces the consistency of the solution with the available binary data under total-variation regularization. Finally, we derive an efficient reconstruction algorithm relying on convex-optimization principles. We conduct several experiments on standard images and demonstrate the practical interest of our approach.

We present a new fast active-contour model (a.k.a. snake) for image segmentation in 3D microscopy. We introduce a parametric design that relies on exponential B-spline bases and allows us to build snakes that are able to reproduce ellipsoids. We design our bases to have the shortest-possible support, subject to some constraints. Thus, computational efficiency is maximized. The proposed 3D snake can approximate blob-like objects with good accuracy and can perfectly reproduce spheres and ellipsoids, irrespective of their position and orientation. The optimization process is remarkably fast due to the use of Gauss' theorem within our energy computation scheme. Our technique yields successful segmentation results, even for challenging data where object contours are not well defined. This is due to our parametric approach that allows one to favor prior shapes. In addition, this paper provides a software that gives full control over the snakes via an intuitive manipulation of few control points.

Poisson inverse problems arise in many modern imaging applications, including biomedical and astronomical ones. The main challenge is to obtain an estimate of the underlying image from a set of measurements degraded by a linear operator and further corrupted by Poisson noise. In this paper, we propose an efficient framework for Poisson image reconstruction, under a regularization approach, which depends on matrix-valued regularization operators. In particular, the employed regularizers involve the Hessian as the regularization operator and Schatten matrix norms as the potential functions. For the solution of the problem, we propose two optimization algorithms that are specifically tailored to the Poisson nature of the noise. These algorithms are based on an augmented-Lagrangian formulation of the problem and correspond to two variants of the alternating direction method of multipliers. Further, we derive a link that relates the proximal map of an l(p) norm with the proximal map of a Schatten matrix norm of order p. This link plays a key role in the development of one of the proposed algorithms. Finally, we provide experimental results on natural and biological images for the task of Poisson image deblurring and demonstrate the practical relevance and effectiveness of the proposed framework.

We introduce a novel family of invariant, convex, and non-quadratic functionals that we employ to derive regularized solutions of ill-posed linear inverse imaging problems. The proposed regularizers involve the Schatten norms of the Hessian matrix, which are computed at every pixel of the image. They can be viewed as second-order extensions of the popular total-variation (TV) semi-norm since they satisfy the same invariance properties. Meanwhile, by taking advantage of second-order derivatives, they avoid the staircase effect, a common artifact of TV-based reconstructions, and perform well for a wide range of applications. To solve the corresponding optimization problems, we propose an algorithm that is based on a primal-dual formulation. A fundamental ingredient of this algorithm is the projection of matrices onto Schatten norm balls of arbitrary radius. This operation is performed efficiently based on a direct link we provide between vector projections onto lq norm balls and matrix projections onto Schatten norm balls. Finally, we demonstrate the effectiveness of the proposed methods through experimental results on several inverse imaging problems with real and simulated data.

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.

We investigate the problem of the optimal reconstruction of a generalized Poisson process from its noisy samples. The process is known to have a finite rate of innovation since it is generated by a random stream of Diracs with a finite average number of impulses per unit interval. We formulate the recovery problem in a Bayesian framework and explicitly derive the joint probability density function (pdf) of the sampled signal. We compare the performance of the optimal Minimum Mean Square Error (MMSE) estimator with common regularization techniques such as ℓ1and Log penalty functions. The simulation results indicate that, under certain conditions, the regularization techniques can achieve a performance close to the MMSE method.

Many recent algorithms for sparse signal recovery can be interpreted as maximum-a-posteriori (MAP) estimators relying on some specific priors. From this Bayesian perspective, state-of-the-art methods based on discrete-gradient regularizers, such as total-variation (TV) minimization, implicitly assume the signals to be sampled instances of Lévy processes with independent Laplace-distributed increments. By extending the concept to more general Lévy processes, we propose an efficient minimum-mean-squared error (MMSE) estimation method based on message-passing algorithms on factor graphs. The resulting algorithm can be used to benchmark the performance of the existing or design new algorithms for the recovery of sparse signals.

We consider the denoising of signals and images using regularized least-squares method. In particular, we propose a simple minimization algorithm for regularizers that are functions of the discrete gradient. By exploiting the connection of the discrete gradient with the Haar-wavelet transform, the n-dimensional vector minimization can be decoupled into n scalar minimizations. The proposed method can efficiently solve total-variation (TV) denoising by iteratively shrinking shifted Haar-wavelet transforms. Furthermore, the decoupling naturally lends itself to extensions beyond ℓ1regularizers.

We have recently introduced a class of non-quadratic Hessian-based regularizers as a higher-order extension of the total variation (TV) functional. These regularizers retain some of the most favorable properties of TV while they can effectively deal with the staircase effect that is commonly met in TV-based reconstructions. In this work we propose a novel gradient-based algorithm for the efficient minimization of these functionals under convex constraints. Furthermore, we validate the overall proposed regularization framework for the problem of image deblurring under additive Gaussian noise.

Motivated by the analog nature of real-world signals, we investigate continuous-time random processes. For this purpose, we consider the stochastic processes that can be whitened by linear transformations and we show that the distribution of their samples is necessarily infinitely divisible. As a consequence, such a modeling rules out the Bernoulli-Gaussian distribution since we are able to show in this paper that it is not infinitely divisible. In other words, while the Bernoulli-Gaussian distribution is among the most studied priors for modeling sparse signals, it cannot be associated with any continuous-time stochastic process. Instead, we propose to adapt the priors that correspond to the increments of compound Poisson processes, which are both sparse and infinitely divisible.

We consider the estimation of an i.i.d.\ vector $\xbf \in \R^n$ from measurements $\ybf \in \R^m$ obtained by a general cascade model consisting of a known linear transform followed by a probabilistic componentwise (possibly nonlinear) measurement channel. We present a method, called adaptive generalized approximate message passing (Adaptive GAMP), that enables joint learning of the statistics of the prior and measurement channel along with estimation of the unknown vector $\xbf$. The proposed algorithm is a generalization of a recently-developed method by Vila and Schniter that uses expectation-maximization (EM) iterations where the posteriors in the E-steps are computed via approximate message passing. The techniques can be applied to a large class of learning problems including the learning of sparse priors in compressed sensing or identification of linear-nonlinear cascade models in dynamical systems and neural spiking processes. We prove that for large i.i.d.\ Gaussian transform matrices the asymptotic componentwise behavior of the adaptive GAMP algorithm is predicted by a simple set of scalar state evolution equations. This analysis shows that the adaptive GAMP method can yield asymptotically consistent parameter estimates, which implies that the algorithm achieves a reconstruction quality equivalent to the oracle algorithm that knows the correct parameter values. The adaptive GAMP methodology thus provides a systematic, general and computationally efficient method applicable to a large range of complex linear-nonlinear models with provable guarantees.

We introduce a new method for adaptive one-bit quantization of linear measurements and propose an algorithm for the recovery of signals based on generalized approximate message passing (GAMP). Our method exploits the prior statistical information on the signal for estimating the minimum-mean-squared error solution from one-bit measurements. Our approach allows the one-bit quantizer to use thresholds on the real line. Given the previous measurements, each new threshold is selected so as to partition the consistent region along its centroid computed by GAMP. We demonstrate that the proposed adaptive-quantization scheme with GAMP reconstruction greatly improves the performance of signal and image recovery from one-bit measurements.

We introduce a new wavelet-based method for the implementation of Total-Variation-type denoising. The data term is least-squares, while the regularization term is gradient-based. The particularity of our method is to exploit a link between the discrete gradient and wavelet shrinkage with cycle spinning, which we express by using redundant wavelets. The redundancy of the representation gives us the freedom to enforce additional constraints (e.g., normalization) on the solution to the denoising problem. We perform optimization in an augmented-Lagrangian framework, which decouples the difficultn-dimensional constrained-optimization problem into a sequence ofneasier scalar unconstrained problems that we solve efficiently via traditional wavelet shrinkage. Our method can handle arbitrary gradient-based regularizers. In particular, it can be made to adhere to the popular principle of least total variation. It can also be used as a maximum a posteriori estimator for a variety of priors. We illustrate the performance of our method for image denoising and for the statistical estimation of sparse stochastic processes.

We consider the problem of sampling continuous-time auto-regressive processes on a uniform grid. We investigate whether a given sampled process originates from a single continuous-time model, and address this uniqueness problem by introducing an alternative description of poles in the complex plane. We then utilize Kronecker's approximation theorem and prove that the set of non-unique continuous-time AR(2) models has Lebesgue measure zero in this plane. This is a key aspect in current estimation algorithms that use sampled data, as it allows one to remove the sampling rate constraint that is imposed currently.

We introduce a systematic and practical design for steerable wavelet frames in 3D. Our steerable wavelets are obtained by applying a 3D version of the generalized Riesz transform to a primary isotropic wavelet frame. The novel transform is self-reversible (tight frame) and its elementary constituents (Riesz wavelets) can be efficiently rotated in any 3D direction by forming appropriate linear combinations. Moreover, the basis functions at a given location can be linearly combined to design custom (and adaptive) steerable wavelets. The features of the proposed method are illustrated with the processing and analysis of 3D biomedical data. In particular, we show how those wavelets can be used to characterize directional patterns and to detect edges by means of a 3D monogenic analysis. We also propose a new inverse-problem formalism along with an optimization algorithm for reconstructing 3D images from a sparse set of wavelet-domain edges. The scheme results in high-quality image reconstructions which demonstrate the feature-reduction ability of the steerable wavelets as well as their potential for solving inverse problems.

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.

We present nonquadratic Hessian-based regularization methods that can be effectively used for image restoration problems in a variational framework. Motivated by the great success of the total-variation (TV) functional, we extend it to also include second-order differential operators. Specifically, we derive second-order regularizers that involve matrix norms of the Hessian operator. The definition of these functionals is based on an alternative interpretation of TV that relies on mixed norms of directional derivatives. We show that the resulting regularizers retain some of the most favorable properties of TV, i.e., convexity, homogeneity, rotation, and translation invariance, while dealing effectively with the staircase effect. We further develop an efficient minimization scheme for the corresponding objective functions. The proposed algorithm is of the iteratively reweighted least-square type and results from a majorization-minimization approach. It relies on a problem-specific preconditioned conjugate gradient method, which makes the overall minimization scheme very attractive since it can be applied effectively to large images in a reasonable computational time. We validate the overall proposed regularization framework through deblurring experiments under additive Gaussian noise on standard and biomedical images.

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.

B-splines are attractive basis functions for the continuous-domain representation of biomedical images and volumes. In this paper, we prove that the extended family of box splines are closed under the Radon transform and derive explicit formulae for their transforms. Our results are general; they cover all known brands of compactly-supported box splines (tensor-product B-splines, separable or not) in any number of dimensions. The proposed box spline approach extends to non-Cartesian lattices used for discretizing the image space. In particular, we prove that the 2-D Radon transform of an N-direction box spline is generally a (nonuniform) polynomial spline of degree N-1. The proposed framework allows for a proper discretization of a variety of tomographic reconstruction problems in a box spline basis. It is of relevance for imaging modalities such as X-ray computed tomography and cryo-electron microscopy. We provide experimental results that demonstrate the practical advantages of the box spline formulation for improving the quality and efficiency of tomographic reconstruction algorithms.

The quantitative validation of reconstruction algorithms requires reliable data. Rasterized simulations are popular but they are tainted by an aliasing component that impacts the assessment of the performance of reconstruction. We introduce analytical simulation tools that are suited to parallel magnetic resonance imaging and allow one to build realistic phantoms. The proposed phantoms are composed of ellipses and regions with piecewise-polynomial boundaries, including spline contours, Bézier contours, and polygons. In addition, they take the channel sensitivity into account, for which we investigate two possible models. Our analytical formulations provide well-defined data in both the spatial and k-space domains. Our main contribution is the closed-form determination of the Fourier transforms that are involved. Experiments validate the proposed implementation. In a typical parallel magnetic resonance imaging reconstruction experiment, we quantify the bias in the overly optimistic results obtained with rasterized simulations-the inverse-crime situation. We provide a package that implements the different simulations and provide tools to guide the design of realistic phantoms.

Most bioimaging modalities rely on indirect measurements of the quantity under investigation. The image is obtained as the result of an optimization problem involving a physical model of the measurement system. Due to the ill-posedness of the above problem, the impact of the noise on the reconstructed images must be con trolled. The recent emphasis in biomedical image reconstruction is on regularization schemes that favor sparse solutions, which renders the optimization problem non smooth. In this work, we show how step-size adaptation can be used to speed up the most recent multi-step algorithms (e.g. FISTA) employed in sparse image recovery. We present experiments in MRI and Fluorescence Molecular Tomography with specifically tailored step-adaptation strategies. Our results demonstrate the possibility of an order-of-magnitude speed enhancement over state-of-the-art algorithms.

We suggest a continuous-domain stochastic modeling of images that is invariant to spatial resolution. Specifically, we are proposing an estimator that is calibrated with respect to the sampling step, and that can potentially handle aliased data. Motivated by Markov random fields, we assume a continuous-domain ARMA model and suggest an algorithm for estimating the continuous-domain parameters from the sampled data. The continuous-domain parameters we estimate provide features that can further be used for image classification, segmentation and interpolation, regardless of sampling interval values and of aliasing effects that may appear in the digital image. Experimental results indicate that the proposed approach is preferable over a discrete-domain ARMA modeling.

In this paper, we propose a novel second-order regularizer based on the maximum response of the second-order directional derivative, assuming that the image under consideration belongs to the class of piecewise-linear signals. Compared to total-variation regularization that preserves edges but transforms piecewise-smooth regions into piecewise-constant regions, the proposed model is able to restore piecewise-linear regions and finer details. Deconvolution experiments demonstrate the performance of our approach in terms of the quality of reconstruction.

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.

We propose a new framework to extract the activity-related component in the BOLD functional magnetic resonance imaging (fMRI) signal. As opposed to traditional fMRI signal analysis techniques, we do not impose any prior knowledge of the event timing. Instead, our basic assumption is that the activation pattern is a sequence of short and sparsely distributed stimuli, as is the case in slow event-related fMRI. We introduce new wavelet bases, termed “activelets”, which sparsify the activity-related BOLD signal. These wavelets mimic the behavior of the differential operator underlying the hemodynamic system. To recover the sparse representation, we deploy a sparse-solution search algorithm. The feasibility of the method is evaluated using both synthetic and experimental fMRI data. The importance of the activelet basis and the non-linear sparse recovery algorithm is demonstrated by comparison against classical B-spline wavelets and linear regularization, respectively.

It is well known that spatial averaging can be realized (in space or frequency domain) using algorithms whose complexity does not scale with the size or shape of the filter. These fast algorithms are generally referred to as constant-time or O(1) algorithms in the image-processing literature. Along with the spatial filter, the edge-preserving bilateral filter involves an additional range kernel. This is used to restrict the averaging to those neighborhood pixels whose intensity are similar or close to that of the pixel of interest. The range kernel operates by acting on the pixel intensities. This makes the averaging process nonlinear and computationally intensive, particularly when the spatial filter is large. In this paper, we show how the O(1) averaging algorithms can be leveraged for realizing the bilateral filter in constant time, by using trigonometric range kernels. This is done by generalizing the idea presented by Porikli, i.e., using polynomial kernels. The class of trigonometric kernels turns out to be sufficiently rich, allowing for the approximation of the standard Gaussian bilateral filter. The attractive feature of our approach is that, for a fixed number of terms, the quality of approximation achieved using trigonometric kernels is much superior to that obtained by Porikli using polynomials.

We propose a general methodology (PURE-LET) to design and optimize a wide class of transform-domain thresholding algorithms for denoising images corrupted by mixed Poisson-Gaussian noise. We express the denoising process as a linear expansion of thresholds (LET) that we optimize by relying on a purely data-adaptive unbiased estimate of the mean-squared error (MSE), derived in a non-Bayesian framework (PURE: Poisson-Gaussian unbiased risk estimate). We provide a practical approximation of this theoretical MSE estimate for the tractable optimization of arbitrary transform-domain thresholding. We then propose a pointwise estimator for undecimated filterbank transforms, which consists of subband-adaptive thresholding functions with signal-dependent thresholds that are globally optimized in the image domain. We finally demonstrate the potential of the proposed approach through extensive comparisons with state-of-the-art techniques that are specifically tailored to the estimation of Poisson intensities. We also present denoising results obtained on real images of low-count fluorescence microscopy.

In this paper, we give a general characterization of regularization functionals for vector field reconstruction, based on the requirement that the said functionals satisfy certain geometric invariance properties with respect to transformations of the coordinate system. In preparation for our general result, we also address some commonalities of invariant regularization in scalar and vector settings, and give a complete account of invariant regularization for scalar fields, before focusing on their main points of difference, which lead to a distinct class of regularization operators in the vector case. Finally, as an illustration of potential, we formulate and compare quadratic (L(2)) and total-variation-type (L(1)) regularized denoising of vector fields in the proposed framework.

We present a functional framework for the design of tight steerable wavelet frames in any number of dimensions. The 2-D version of the method can be viewed as a generalization of Simoncelli's steerable pyramid that gives access to a larger palette of steerable wavelets via a suitable parametrization. The backbone of our construction is a primal isotropic wavelet frame that provides the multiresolution decomposition of the signal. The steerable wavelets are obtained by applying a one-to-many mapping (Nth-order generalized Riesz transform) to the primal ones. The shaping of the steerable wavelets is controlled by an M×M unitary matrix (where M is the number of wavelet channels) that can be selected arbitrarily; this allows for a much wider range of solutions than the traditional equiangular configuration (steerable pyramid). We give a complete functional description of these generalized wavelet transforms and derive their steering equations. We describe some concrete examples of transforms, including some built around a Mallat-type multiresolution analysis of L(2)(R(d)), and provide a fast Fourier transform-based decomposition algorithm. We also propose a principal-component-based method for signal-adapted wavelet design. Finally, we present some illustrative examples together with a comparison of the denoising performance of various brands of steerable transforms. The results are in favor of an optimized wavelet design (equalized principal component analysis), which consistently performs best.

In this paper we propose a method based on (2, 1)-mixed-norm penalization for incorporating a structural prior in FDOT image reconstruction. The effect of (2, 1)-mixed-norm penalization is twofold: first, a sparsifying effect which isolates few anatomical regions where the fluorescent probe has accumulated, and second, a regularization effect inside the selected anatomical regions. After formulating the reconstruction in a variational framework, we analyze the resulting optimization problem and derive a practical numerical method tailored to (2, 1)-mixed-norm regularization. The proposed method includes as particular cases other sparsity promoting regularization methods such as l(1)-norm penalization and total variation penalization. Results on synthetic and experimental data are presented.

In this work, we exploit the fact that wavelets can represent magnetic resonance images well, with relatively few coefficients. We use this property to improve magnetic resonance imaging (MRI) reconstructions from undersampled data with arbitrary k-space trajectories. Reconstruction is posed as an optimization problem that could be solved with the iterative shrinkage/thresholding algorithm (ISTA) which, unfortunately, converges slowly. To make the approach more practical, we propose a variant that combines recent improvements in convex optimization and that can be tuned to a given specific k-space trajectory. We present a mathematical analysis that explains the performance of the algorithms. Using simulated and in vivo data, we show that our nonlinear method is fast, as it accelerates ISTA by almost two orders of magnitude. We also show that it remains competitive with TV regularization in terms of image quality.

Many practical applications require the reconstruction of images from irregularly sampled data. The spline formalism offers an attractive framework for solving this problem; the currently available methods, however, are hard to deploy for large-scale interpolation problems in dimensions greater than two (3-D, 3-D+time) because of an exponential increase of their computational cost (curse of dimensionality). Here, we revisit the standard regularized least-squares formulation of the interpolation problem, and propose to perform the reconstruction in a uniform tensor-product B-spline basis as an alternative to the classical solution involving radial basis functions. Our analysis reveals that the underlying multilinear system of equations admits a tensor decomposition with an extreme sparsity of its one dimensional components. We exploit this property for implementing a parallel, memory-efficient system solver. We show that the computational complexity of the proposed algorithm is essentially linear in the number of measurements and that its dependency on the number of dimensions is significantly less than that of the original sparse matrix-based implementation. The net benefit is a substantial reduction in memory requirement and operation count when compared to standard matrix-based algorithms, so that even 4-D problems with millions of samples become computationally feasible on desktop PCs in reasonable time. After validating the proposed algorithm in 3-D and 4-D, we apply it to a concrete imaging problem: the reconstruction of medical ultrasound images (3-D+time) from a large set of irregularly sampled measurements, acquired by a fast rotating ultrasound transducer.

Modern deconvolution algorithms are often specified as minimization problems involving a non-quadratic regularization functional. When the latter is a wavelet-domain ℓ1-norm that favors sparse solutions, the problem can be solved by a simple iterative shrinkage/thresholding algorithm (ISTA). This approach provides state-of-the-art results in 2-D, but is harder to deploy in 3-D because of its slow convergence. In this paper, we propose an acceleration scheme that turns wavelet-regularized deconvolution into a competitive solution for 3-D fluorescence microscopy. A significant speed-up is achieved though a synergistic combination of subband-adapted thresholds and sequential TwIST updates. We provide a theoretical justification of the procedure together with an experimental evaluation, including the application to real 3-D fluorescence data.

We propose a novel algorithm for denoising Poisson-corrupted images, that performs a signal-adaptive thresholding of the undecimated Haar wavelet coefficients. A Poisson's unbiased MSE estimate is devised and adapted to arbitrary transform-domain pointwise processing. This prior-free quadratic measure of quality is then used to globally optimize a linearly parameterized subband-adaptive thresholding, which accounts for the signal-dependent noise variance. We demonstrate the qualitative and computational competitiveness of the resulting denoising algorithm through comprehensive comparisons with some state-of-the-art multiscale techniques specifically designed for Poisson intensity estimation. We also show promising denoising results obtained on low-count fluorescence microscopy images.

In this note we introduce a vector generalization of fractional Brownian motion. Our definition takes into account directional properties of vector fields-such as divergence, rotational behaviour, and interactions with coordinate transformations-that have no counterpart in the scalar setting. Apart from the Hurst exponent which dictates the scale-dependent structure of the field, additional parameters of the new model control the balance between solenoidal and irrotational behaviour. This level of versatility makes these random fields potentially interesting candidates for the stochastic modelling of physical phenomena in various fields of application such as fluid dynamics, field theory, and medical image processing.

We propose an efficient orthonormal wavelet-domain video denoising algorithm based on an appropriate integration of motion compensation into an adapted version of our recently devised Stein's unbiased risk estimator-linear expansion of thresholds (SURE-LET) approach. To take full advantage of the strong spatio-temporal correlations of neighboring frames, a global motion compensation followed by a selective block-matching is first applied to adjacent frames, which increases their temporal correlations without distorting the interframe noise statistics. Then, a multiframe interscale wavelet thresholding is performed to denoise the current central frame. The simulations we made on standard grayscale video sequences for various noise levels demonstrate the efficiency of the proposed solution in reducing additive white Gaussian noise. Obtained at a lighter computational load, our results are even competitive with most state-of-the-art redundant wavelet-based techniques. By using a cycle-spinning strategy, our algorithm is in fact able to outperform these methods.

The efficient realization of linear space-variant (non-convolution) filters is a challenging computational problem in image processing. In this paper, we demonstrate that it is possible to filter an image with a Gaussian-like elliptic window of varying size, elongation and orientation using a fixed number of computations per pixel. The associated algorithm, which is based upon a family of smooth compactly supported piecewise polynomials, the radially-uniform box splines, is realized using preintegration and local finite-differences. The radially-uniform box splines are constructed through the repeated convolution of a fixed number of box distributions, which have been suitably scaled and distributed radially in an uniform fashion. The attractive features of these box splines are their asymptotic behavior, their simple covariance structure, and their quasi-separability. They converge to Gaussians with the increase of their order, and are used to approximate anisotropic Gaussians of varying covariance simply by controlling the scales of the constituent box distributions. Based upon the second feature, we develop a technique for continuously controlling the size, elongation and orientation of these Gaussian-like functions. Finally, the quasi-separable structure, along with a certain scaling property of box distributions, is used to efficiently realize the associated space-variant elliptical filtering, which requires O(1) computations per pixel irrespective of the shape and size of the filter.

Our main goal in this paper is to set the foundations of a general continuous-domain framework for designing steerable, reversible signal transformations (a.k.a. frames) in multiple dimensions ( d >or= 2). To that end, we introduce a self-reversible, Nth-order extension of the Riesz transform. We prove that this generalized transform has the following remarkable properties: shift-invariance, scale-invariance, inner-product preservation, and steerability. The pleasing consequence is that the transform maps any primary wavelet frame (or basis) of [Formula: see text] into another "steerable" wavelet frame, while preserving the frame bounds. The concept provides a functional counterpart to Simoncelli's steerable pyramid whose construction was primarily based on filterbank design. The proposed mechanism allows for the specification of wavelets with any order of steerability in any number of dimensions; it also yields a perfect reconstruction filterbank algorithm. We illustrate the method with the design of a novel family of multidimensional Riesz-Laplace wavelets that essentially behave like the N th-order partial derivatives of an isotropic Gaussian kernel.

Reconstruction algorithms for fluorescence tomography have to address two crucial issues: 1) the ill-posedness of the reconstruction problem, 2) the large scale of numerical problems arising from imaging of 3-D samples. Our contribution is the design and implementation of a reconstruction algorithm that incorporates general Lp regularization (p ¿ 1). The originality of this work lies in the application of general Lp constraints to fluorescence tomography, combined with an efficient matrix-free strategy that enables the algorithm to deal with large reconstruction problems at reduced memory and computational costs. In the experimental part, we specialize the application of the algorithm to the case of sparsity promoting constraints (L (1)). We validate the adequacy of L (1) regularization for the investigation of phenomena that are well described by a sparse model, using data acquired during phantom experiments.

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.

We provide an amplitude-phase representation of the dual-tree complex wavelet transform by extending the fixed quadrature relationship of the dual-tree wavelets to arbitrary phase-shifts using the fractional Hilbert transform (fHT). The fHT is a generalization of the Hilbert transform that extends the quadrature phase-shift action of the latter to arbitrary phase-shifts a real shift parameter controls this phase-shift action. Next, based on the proposed representation and the observation that the fHT operator maps well-localized B-spline wavelets (that resemble Gaussian-windowed sinusoids) into B-spline wavelets of the same order but different shift, we relate the corresponding dual-tree scheme to the paradigm of multiresolution windowed Fourier analysis.

We provide an extension of the L2-spline pyramid (Unser et al., 1993) using polyharmonic splines. We analytically prove that the corresponding error pyramid behaves exactly as a multi-scale Laplace operator. We use the multiresolution properties of polyharmonic splines to derive an efficient, non-separable filterbank implementation. Finally, we illustrate the potentials of our pyramid by performing an estimation of the parameters of multivariate fractal processes.

We introduce an Nth-order extension of the Riesz transform in d dimensions. We prove that this generalized transform has the following remarkable properties: shift-invariance, scale-invariance, inner-product preservation, and steerability. The pleasing consequence is that the transform maps any primary wavelet frame (or basis) of L2(¿d) into another ¿steerable¿ wavelet frame, while preserving the frame bounds. The concept provides a rigorous functional counterpart to Simoncelli's steerable pyramid whose construction was entirely based on digital filter design. The proposed mechanism allows for the specification of wavelets with any order of steerability in any number of dimensions; it also yields a perfect reconstruction filterbank algorithm. We illustrate the method using a Mexican-hat-like polyharmonic spline wavelet transform as our primary frame.

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.

In this contribution, we study the notion of affine invariance (specifically, invariance to the shifting, scaling, and rotation of the coordinate system) as a starting point for the development of mathematical tools and approaches useful in the characterization and analysis of multivariate fractional Brownian motion (fBm) fields. In particular, using a rigorous and powerful distribution theoretic formulation, we extend previous results of Blu and Unser (2006) to the multivariate case, showing that polyharmonic splines and fBm processes can be seen as the (deterministic vs stochastic) solutions to an identical fractional partial differential equation that involves a fractional Laplacian operator. We then show that wavelets derived from polyharmonic splines have a behavior similar to the fractional Laplacian, which also turns out to be the whitening operator for fBm fields. This fact allows us to study the probabilistic properties of the wavelet transform coefficients of fBm-like processes, leading for instance to ways of estimating the Hurst exponent of a multiparameter process from its wavelet transform coefficients. We provide theoretical and experimental verification of these results. To complement the toolbox available for multiresolution processing of stochastic fractals, we also introduce an extended family of multidimensional multiresolution spaces for a large class of (separable and nonseparable) lattices of arbitrary dimensionality.

The monogenic signal is the natural 2-D counterpart of the 1-D analytic signal. We propose to transpose the concept to the wavelet domain by considering a complexified version of the Riesz transform which has the remarkable property of mapping a real-valued (primary) wavelet basis of L(2) (R(2)) into a complex one. The Riesz operator is also steerable in the sense that it give access to the Hilbert transform of the signal along any orientation. Having set those foundations, we specify a primary polyharmonic spline wavelet basis of L(2) (R(2)) that involves a single Mexican-hat-like mother wavelet (Laplacian of a B-spline). The important point is that our primary wavelets are quasi-isotropic: they behave like multiscale versions of the fractional Laplace operator from which they are derived, which ensures steerability. We propose to pair these real-valued basis functions with their complex Riesz counterparts to specify a multiresolution monogenic signal analysis. This yields a representation where each wavelet index is associated with a local orientation, an amplitude and a phase. We give a corresponding wavelet-domain method for estimating the underlying instantaneous frequency. We also provide a mechanism for improving the shift and rotation-invariance of the wavelet decomposition and show how to implement the transform efficiently using perfect-reconstruction filterbanks. We illustrate the specific feature-extraction capabilities of the representation and present novel examples of wavelet-domain processing; in particular, a robust, tensor-based analysis of directional image patterns, the demodulation of interferograms, and the reconstruction of digital holograms.

We present a multilevel extension of the popular "thresholded Landweber" algorithm for wavelet-regularized image restoration that yields an order of magnitude speed improvement over the standard fixed-scale implementation. The method is generic and targeted towards large-scale linear inverse problems, such as 3-D deconvolution microscopy. The algorithm is derived within the framework of bound optimization. The key idea is to successively update the coefficients in the various wavelet channels using fixed, subband-adapted iteration parameters (step sizes and threshold levels). The optimization problem is solved efficiently via a proper chaining of basic iteration modules. The higher level description of the algorithm is similar to that of a multigrid solver for PDEs, but there is one fundamental difference: the latter iterates though a sequence of multiresolution versions of the original problem, while, in our case, we cycle through the wavelet subspaces corresponding to the difference between successive approximations. This strategy is motivated by the special structure of the problem and the preconditioning properties of the wavelet representation. We establish that the solution of the restoration problem corresponds to a fixed point of our multilevel optimizer. We also provide experimental evidence that the improvement in convergence rate is essentially determined by the (unconstrained) linear part of the algorithm, irrespective of the type of wavelet. Finally, we illustrate the technique with some image deconvolution examples, including some real 3-D fluorescence microscopy data.

We propose a novel method of constructing exact Hilbert transform (HT) pairs of wavelet bases using fractional B- splines and state necessary and sufficient conditions for generating such wavelet pairs. In particular, we demonstrate how HT pairs of biorthogonal wavelet bases of L2(K) can be constructed using well-localized scaling functions with identical Riesz bounds. Finally, we illustrate this concept by constructing a family of analytic Gabor-like wavelets that exhibit near optimal time-frequency localization.

We consider the problem of optimizing the parameters of an arbitrary denoising algorithm by minimizing Stein's unbiased risk estimate (SURE) which provides a means of assessing the true mean-squared-error (MSE) purely from the measured data assuming that it is corrupted by Gaussian noise. To accomplish this, we propose a novel Monte-Carlo technique based on a black-box approach which enables the user to compute SURE for an arbitrary denoising algorithm with some specific parameter setting. Our method only requires the response of the denoising algorithm to additional input noise and does not ask for any information about the functional form of the corresponding denoising operator. This, therefore, permits SURE-based optimization of a wide variety of denoising algorithms (global-iterative, pointwise, etc). We present experimental results to justify our claims.

Recently, we proposed a noniterative cepstral technique for exact signal recovery in frequency-domain optical-coherence tomography. In this paper, we address the influence of measurement noise on the performance of the method. We derive analytical expressions for the bias and variance of the tomogram under a small noise approximation, and show that our technique yields unbiased and consistent estimators, which have a variance that is proportional to that of the noise and inversely proportional to the data size. We present simulation results to confirm the theoretical derivations. We also derive approximate Cramer-Rao bounds (CRBs) on the achievable accuracy of reconstruction.

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.

We demonstrate that it is possible to filter an image with an elliptic window of varying size, elongation and orientation with a fixed computational cost per pixel. Our method involves the application of a suitable global pre-integrator followed by a pointwise-adaptive localization mesh. We present the basic theory for the ID case using a B-spline formalism and then appropriately extend it to 2D using radially-uniform box splines. The size and ellipticity of these radially-uniform box splines is adaptively controlled. Moreover, they converge to Gaussians as the order increases. Finally, we present a fast and practical directional filtering algorithm that has the capability of adapting to the local image features.

We introduce a new semi-orthogonal complex wavelet basis of L2(R2). The basis functions are associated to the complex gradient-Laplace operator, which plays a central role in image processing. We define analytically a single-generator wavelet that is shifted on the coset positions of the subsampling matrix. Next, we propose the "wavelet Marr pyramid" for an extension of the new basis that achieves near shift-invariance and steerability (using a Gaussian-like smoothing kernel), for a mild redundancy factor only. This new wavelet pyramid decomposition closely mimicks the basic operations of Marx's framework for early vision. The pyramid is implemented by a fast filterbank algorithm using the FFT.

We propose a recursive data-driven risk-estimation method for non-linear iterative deconvolution. Our two main contributions are 1) a solution-domain risk-estimation approach that is applicable to non-linear restoration algorithms for ill- conditioned inverse problems; and 2) a risk estimate for a state-of-the-art iterative procedure, the thresholded Landweber iteration, which enforces a wavelet-domain sparsity constraint. Our method can be used to estimate the SNR improvement at every step of the algorithm; e.g., for stopping the iteration after the highest value is reached. It can also be applied to estimate the optimal threshold level for a given number of iterations.

A fast computational method is given for the Fourier transform of the polyharmonic B-spline autocorrelation sequence in d dimensions. The approximation error is exponentially decaying with the number of terms taken into account. The algorithm improves speed upon a simple truncated-sum approach. Moreover, it is virtually independent of the spline's order. The autocorrelation filter directly serves for various tasks related to polyharmonic splines, such as interpolation, orthonormalization, and wavelet basis design.

The problem of sampling signals that are not admissible within the classical Shannon framework has received much attention in the recent past. Typically, these signals have a parametric representation with a finite number of degrees of freedom per time unit. It was shown that, by choosing suitable sampling kernels, the parameters can be computed by employing high-resolution spectral estimation techniques. In this letter, we propose a simple acquisition and reconstruction method within the framework of multichannel sampling. In the proposed approach, an infinite stream of nonuniformly-spaced Dirac impulses can be sampled and accurately reconstructed provided that there is at most one Dirac impulse per sampling period. The reconstruction algorithm has a low computational complexity, and the parameters are computed on the fly. The processing delay is minimal just the sampling period. We propose sampling circuits using inexpensive passive devices such as resistors and capacitors. We also show how the approach can be extended to sample piecewise-constant signals with a minimal change in the system configuration. We provide some simulation results to confirm the theoretical findings.

Due to the limited depth of field of brightfield microscopes, it is usually impossible to image thick specimens entirely in focus. By optically sectioning the specimen, the in-focus information at the specimen's surface can be acquired over a range of images. Commonly based on a high-pass criterion, extended-depth-of-field methods aim at combining the in-focus information from these images into a single image of the texture on the specimen's surface. The topography provided by such methods is usually limited to a map of selected in-focus pixel positions and is inherently discretized along the axial direction, which limits its use for quantitative evaluation. In this paper, we propose a method that jointly estimates the texture and topography of a specimen from a series of brightfield optical sections; it is based on an image formation model that is described by the convolution of a thick specimen model with the microscope's point spread function. The problem is stated as a least-squares minimization where the texture and topography are updated alternately. This method also acts as a deconvolution when the in-focus PSF has a blurring effect, or when the true in-focus position falls in between two optical sections. Comparisons to state-of-the-art algorithms and experimental results demonstrate the potential of the proposed approach.

We consider the problem of optimizing the parameters of a given denoising algorithm for restoration of a signal corrupted by white Gaussian noise. To achieve this, we propose to minimize Stein's unbiased risk estimate (SURE) which provides a means of assessing the true mean-squared error (MSE) purely from the measured data without need for any knowledge about the noise-free signal. Specifically, we present a novel Monte-Carlo technique which enables the user to calculate SURE for an arbitrary denoising algorithm characterized by some specific parameter setting. Our method is a black-box approach which solely uses the response of the denoising operator to additional input noise and does not ask for any information about its functional form. This, therefore, permits the use of SURE for optimization of a wide variety of denoising algorithms. We justify our claims by presenting experimental results for SURE-based optimization of a series of popular image-denoising algorithms such as total-variation denoising, wavelet soft-thresholding, and Wiener filtering/smoothing splines. In the process, we also compare the performance of these methods. We demonstrate numerically that SURE computed using the new approach accurately predicts the true MSE for all the considered algorithms. We also show that SURE uncovers the optimal values of the parameters in all cases.

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.

A distinction is usually made between wavelet bases and wavelet frames. The former are associated with a one-to-one representation of signals, which is somewhat constrained but most efficient computationally. The latter are over-complete, but they offer advantages in terms of flexibility (shape of the basis functions) and shift-invariance. In this paper, we propose a framework for improved wavelet analysis based on an appropriate pairing of a wavelet basis with a mildly redundant version of itself (frame). The processing is accomplished in four steps: 1) redundant wavelet analysis, 2) wavelet-domain processing, 3) projection of the results onto the wavelet basis, and 4) reconstruction of the signal from its nonredundant wavelet expansion. The wavelet analysis is pyramid-like and is obtained by simple modification of Mallat's filterbank algorithm (e.g., suppression of the down-sampling in the wavelet channels only). The key component of the method is the subband regression filter (Step 3) which computes a wavelet expansion that is maximally consistent in the least squares sense with the redundant wavelet analysis. We demonstrate that this approach significantly improves the performance of soft-threshold wavelet denoising with a moderate increase in computational cost. We also show that the analysis filters in the proposed framework can be adjusted for improved feature detection; in particular, a new quincunx Mexican-hat-like wavelet transform that is fully reversible and essentially behaves the (gamma/2)th Laplacian of a Gaussian.

Our aim in this paper is to tighten the link between wavelets, some classical image-processing operators, and David Marr's theory of early vision. The cornerstone of our approach is a new complex wavelet basis that behaves like a smoothed version of the Gradient-Laplace operator. Starting from first principles, we show that a single-generator wavelet can be defined analytically and that it yields a semi-orthogonal complex basis of L2 (R2), irrespective of the dilation matrix used. We also provide an efficient FFT-based filterbank implementation. We then propose a slightly redundant version of the transform that is nearly translation-invariant and that is optimized for better steerability (Gaussian-like smoothing kernel).We call it the Marr-like wavelet pyramid because it essentially replicates the processing steps in Marr's theory of early vision.We use it to derive a primal wavelet sketch which is a compact description of the image by a multiscale, subsampled edge map. Finally, we provide an efficient iterative algorithm for the reconstruction of an image from its primal wavelet sketch.

We present a fast variational deconvolution algorithm that minimizes a quadratic data term subject to a regularization on the l(1)-norm of the wavelet coefficients of the solution. Previously available methods have essentially consisted in alternating between a Landweber iteration and a wavelet-domain soft-thresholding operation. While having the advantage of simplicity, they are known to converge slowly. By expressing the cost functional in a Shannon wavelet basis, we are able to decompose the problem into a series of subband-dependent minimizations. In particular, this allows for larger (subband-dependent) step sizes and threshold levels than the previous method. This improves the convergence properties of the algorithm significantly. We demonstrate a speed-up of one order of magnitude in practical situations. This makes wavelet-regularized deconvolution more widely accessible, even for applications with a strong limitation on computational complexity. We present promising results in 3-D deconvolution microscopy, where the size of typical data sets does not permit more than a few tens of iterations.

Tomographic reconstruction from positron emission tomography (PET) data is an ill-posed problem that requires regularization. An attractive approach is to impose an l(1) -regularization constraint, which favors sparse solutions in the wavelet domain. This can be achieved quite efficiently thanks to the iterative algorithm developed by Daubechies et al., 2004. In this paper, we apply this technique and extend it for the reconstruction of dynamic (spatio-temporal) PET data. Moreover, instead of using classical wavelets in the temporal dimension, we introduce exponential-spline wavelets (E-spline wavelets) that are specially tailored to model time activity curves (TACs) in PET. We show that the exponential-spline wavelets naturally arise from the compartmental description of the dynamics of the tracer distribution. We address the issue of the selection of the "optimal" E-spline parameters (poles and zeros) and we investigate their effect on reconstruction quality. We demonstrate the usefulness of spatio-temporal regularization and the superior performance of E-spline wavelets over conventional Battle-LemariE wavelets in a series of experiments: the 1-D fitting of TACs, and the tomographic reconstruction of both simulated and clinical data. We find that the E-spline wavelets outperform the conventional wavelets in terms of the reconstructed signal-to-noise ratio (SNR) and the sparsity of the wavelet coefficients. Based on our simulations, we conclude that replacing the conventional wavelets with E-spline wavelets leads to equal reconstruction quality for a 40% reduction of detected coincidences, meaning an improved image quality for the same number of counts or equivalently a reduced exposure to the patient for the same image quality.

Frequency domain optical coherence tomography (FDOCT) is a new technique that is well-suited for fast imaging of biological specimens, as well as non-biological objects. The measurements are in the frequency domain, and the objective is to retrieve an artifact-free spatial domain description of the specimen. In this paper, we develop a new technique for model-based retrieval of spatial domain data from the frequency domain data. We use a piecewise-constant model for the refractive index profile that is suitable for multi-layered specimens. We show that the estimation of the layered structure parameters can be mapped into a harmonic retrieval problem, which enables us to use high-resolution spectrum estimation techniques. The new technique that we propose is efficient and requires few measurements. We also analyze the effect of additive measurement noise on the algorithm performance. The experimental results show that the technique gives highly accurate parameter estimates. For example, at 25 dB signal-to-noise ratio, the mean square error in the position estimate is about 0.01 % of the actual value.

This paper introduces a new approach to orthonormal wavelet image denoising. Instead of postulating a statistical model for the wavelet coefficients, we directly parametrize the denoising process as a sum of elementary nonlinear processes with unknown weights. We then minimize an estimate of the mean square error between the clean image and the denoised one. The key point is that we have at our disposal a very accurate, statistically unbiased, MSE estimate--Stein's unbiased risk estimate--that depends on the noisy image alone, not on the clean one. Like the MSE, this estimate is quadratic in the unknown weights, and its minimization amounts to solving a linear system of equations. The existence of this a priori estimate makes it unnecessary to devise a specific statistical model for the wavelet coefficients. Instead, and contrary to the custom in the literature, these coefficients are not considered random anymore. We describe an interscale orthonormal wavelet thresholding algorithm based on this new approach and show its near-optimal performance--both regarding quality and CPU requirement--by comparing it with the results of three state-of-the-art nonredundant denoising algorithms on a large set of test images. An interesting fallout of this study is the development of a new, group-delay-based, parent-child prediction in a wavelet dyadic tree.

We present a new computational method for reconstructing a vector velocity field from scattered, pulsed-wave ultrasound Doppler data. The main difficulty is that the Doppler measurements are incomplete, for they do only capture the velocity component along the beam direction. We thus propose to combine measurements from different beam directions. However, this is not yet sufficient to make the problem well posed because 1) the angle between the directions is typically small and 2) the data is noisy and nonuniformly sampled. We propose to solve this reconstruction problem in the continuous domain using regularization. The reconstruction is formulated as the minimizer of a cost that is a weighted sum of two terms: 1) the sum of squared difference between the Doppler data and the projected velocities 2) a quadratic regularization functional that imposes some smoothness on the velocity field. We express our solution for this minimization problem in a B-spline basis, obtaining a sparse system of equations that can be solved efficiently. Using synthetic phantom data, we demonstrate the significance of tuning the regularization according to the a priori knowledge about the physical property of the motion. Next, we validate our method using real phantom data for which the ground truth is known. We then present reconstruction results obtained from clinical data that originate from 1) blood flow in carotid bifurcation and 2) cardiac wall motion.

Magnetic resonance spectroscopy imaging (MRSI) is an attractive tool for medical imaging. However, its practical use is often limited by the intrinsic low spatial resolution and long acquisition time. Spectral localization by imaging (SLIM) has been proposed as a non-Fourier reconstruction algorithm that incorporates spatial a priori information about spectroscopically uniform compartments. Unfortunately, the influence of the magnetic field inhomogeneity--in particular, the susceptibility effects at tissues' boundaries--undermines the validity of the compartmental model. Therefore, we propose BSLIM as an extension of SLIM with field inhomogeneity compensation. A B0-field inhomogeneity map, which can be acquired rapidly and at high resolution, is used by the new algorithm as additional a priori information. We show that the proposed method is distinct from the generalized SLIM (GSLIM) framework. Experimental results of a two-compartment phantom demonstrate the feasibility of the method and the importance of inhomogeneity compensation.

We consider the problem of estimating a fractional Brownian motion known only from its noisy samples at the integers. We show that the optimal estimator can be expressed using a digital Wiener-like filter followed by a simple time-variant correction accounting for nonstationarity. Moreover, we prove that this estimate lives in a symmetric fractional spline space and give a practical implementation for optimal upsampling of noisy fBm samples by integer factors

The Matern class is a parametric family of autocorrelation functions that is commonly used in geostatistics. We argue that a generalized, anisotropic version of this model is suitable for capturing the correlation structure of a variety of natural images. We specify the optimal space for the MMSE reconstruction of stochastic Matern signals from their uniformly-sampled noisy measurements (generalized sampling problem). We prove that the optimal reconstruction space is generated by the multi-integer shifts of a Matern function which form a Riesz basis. Based on this representation, we propose a practical filter-based reconstruction method that relies on the prior identification of the Matern parameters from the measured data. We present experimental results to justify the use of the Matern model and to demonstrate the performance of our signal-adapted reconstruction technique

Many experimental paradigms in biology aim at studying the response to coordinated stimuli. In dynamic imaging experiments, the observed data is often not straightforward to interpret and not directly measurable in a quantitative fashion. Consequently, the data is typically preprocessed in an ad hoc fashion and the results subjected to a statistical inference at the level of a population. We propose a new framework for analyzing time-lapse images that exploits some a priori knowledge on the type of temporal response and takes advantage of the spatial correlation of the data. This is achieved by processing the data in the wavelet domain and expressing the time course of each wavelet coefficient by a linear model. We end up with a statistical map in the spatial domain for the contrast of interest (i.e., the stimulus response). The feasibility of the method is demonstrated by an example of intrinsic microscopy imaging of mice's brains during coordinated sensory stimulation

Recently, we have proposed a new framework for detecting brain activity from fMRI data, which is based on the spatial discrete wavelet transform. The standard wavelet-based approach performs a statistical test in the wavelet domain, and therefore fails to provide a rigorous statistical interpretation in the spatial domain. The new framework provides an "integrated" approach: the data is processed in the wavelet domain (by thresholding wavelet coefficients), and a suitable statistical testing procedure is applied afterwards in the spatial domain. This method is based on conservative assumptions only and has a strong type-I error control by construction. At the same time, it has a sensitivity comparable to that of SPM. Here, we discuss the extension of our algorithm to the redundant discrete wavelet transform, which provides a shift-invariant detection scheme. The key features of our technique are illustrated with experimental results. An implementation of our framework is available as a toolbox (WSPM) for the SPM2 software

Three-directional box-splines are particularly well-suited to interpolate and approximate hexagonally sampled data. In this paper, we propose a computationally efficient end-to-end reconstruction process. First, we introduce a prefiltering step that is based on a quasi-interpolation scheme using low-complexity finite-impulse-response (FIR) filters. Second, we derive a closed analytical expression for three-directional box-splines of any order that leads to a fast evaluation of the spline surface. All operations act locally on the data, and thus are well adapted to applications dealing with large images. To demonstrate the feasibility of our method, we implemented the complete procedure and we present experimental results.

We propose here a new pointwise wavelet thresholding function that incorporates inter-scale dependencies. This non-linear function depends on a set of four linear parameters per sub-band which are set by minimizing Stein's unbiased MSE estimate (SURE). Our approach assumes additive Gaussian white noise. In order for the inter-scale dependencies to be faithfully taken into account, we also develop a rigorous feature alignment processing, that is adapted to arbitrary wavelet filters (e.g. non-symmetric filters). Finally, we demonstrate the efficiency of our denoising approach in simulations over a wide range of noise levels for a representative set of standard images.

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.

Starting from the power spectral density of Mateacutern stochastic processes, we introduce a new family of splines that is defined in terms of the whitening operator of such processes. We show that these Mateacutern splines admit a stable representation in a B-spline-like basis. We specify the Mateacutern B-splines (causal and symmetric) and identify their key properties; in particular, we prove that these generate a Riesz basis and that they can be written as a product of an exponential with a fractional polynomial B-spline. We also indicate how these new functions bridge the gap between the fractional polynomial splines and the cardinal exponential ones. We then show that these splines provide the optimal reconstruction space for the minimum mean-squared error estimation of Mateacutern signals from their noisy samples. We also propose a digital Wiener-filter-like algorithm for the efficient determination of the optimal B-spline coefficients

We introduce a three-dimensional (3-D) parametric active contour algorithm for the shape estimation of DNA molecules from stereo cryo-electron micrographs. We estimate the shape by matching the projections of a 3-D global shape model with the micrographs; we choose the global model as a 3-D filament with a B-spline skeleton and a specified radial profile. The active contour algorithm iteratively updates the B-spline coefficients, which requires us to evaluate the projections and match them with the micrographs at every iteration. Since the evaluation of the projections of the global model is computationally expensive, we propose a fast algorithm based on locally approximating it by elongated blob-like templates. We introduce the concept of projection-steerability and derive a projection-steerable elongated template. Since the two-dimensional projections of such a blob at any 3-D orientation can be expressed as a linear combination of a few basis functions, matching the projections of such a 3-D template involves evaluating a weighted sum of inner products between the basis functions and the micrographs. The weights are simple functions of the 3-D orientation and the inner-products are evaluated efficiently by separable filtering. We choose an internal energy term that penalizes the average curvature magnitude. Since the exact length of the DNA molecule is known a priori, we introduce a constraint energy term that forces the curve to have this specified length. The sum of these energies along with the image energy derived from the matching process is minimized using the conjugate gradients algorithm. We validate the algorithm using real, as well as simulated, data and show that it performs well.

Motivated by the fractal-like behavior of natural images, we develop a smoothing technique that uses a regularization functional which is a fractional iterate of the Laplacian. This type of functional was initially introduced by Duchon for the approximation of nonuniformily sampled, multidimensional data. He proved that the general solution is a smoothing spline that is represented by a linear combination of radial basis functions (RBFs). Unfortunately, this is tedious to implement for images because of the poor conditioning of RBFs and their lack of decay. Here, we present a much more efficient method for the special case of a uniform grid. The key idea is to express Duchon's solution in a fractional polyharmonic B-spline basis that spans the same space as the RBFs. This allows us to derive an algorithm where the smoothing is performed by filtering in the Fourier domain. Next, we prove that the above smoothing spline can be optimally tuned to provide the MMSE estimation of a fractional Brownian field corrupted by white noise. This is a strong result that not only yields the best linear filter (Wiener solution), but also the optimal interpolation space, which is not bandlimited. It also suggests a way of using the noisy data to identify the optimal parameters (order of the spline and smoothing strength), which yields a fully automatic smoothing procedure. We evaluate the performance of our algorithm by comparing it against an oracle Wiener filter, which requires the knowledge of the true noiseless power spectrum of the signal. We find that our approach performs almost as well as the oracle solution over a wide range of conditions.

We build a multiresolution analysis based on shift-invariant exponential B-spline spaces. We construct the basis functions for these spaces and for their orthogonal complements. This yields a new family of wavelet-like basis functions of L/sub 2/, with some remarkable properties. The wavelets, which are characterized by a set of poles and zeros, have an explicit analytical form (exponential spline). They are nonstationary is the sense that they are scale-dependent and that they are not necessarily the dilates of one another. They behave like multi-scale versions of some underlying differential operator L; in particular, they are orthogonal to the exponentials that are in the null space of L. The corresponding wavelet transforms are implemented efficiently using an adaptation of Mallat's (1998) filterbank algorithm.

We present a generalization of the Daubechies wavelet family. The context is that of a non-stationary multiresolution analysis - i.e., a sequence of embedded approximation spaces generated by scaling functions that are not necessarily dilates of one another. The constraints that we impose on these scaling functions are: (1) orthogonality with respect to translation; (2) reproduction of a given set of exponential polynomials; and (3) minimal support. These design requirements lead to the construction of a general family of compactly-supported, orthonormal wavelet-like bases of L/sub 2/. If the exponential parameters are all zero, then one recovers Daubechies wavelets, which are orthogonal to the polynomials of degree (N-1) where N is the order (vanishing-moment property). A fast filterbank implementation of the generalized wavelet transform follows naturally; it is similar to Mallat's algorithm, except that the filters are now scale-dependent. The new transforms offer increased flexibility and are tunable to the spectral characteristics of a wide class of signals.

We present a framework for feature detection in 3-D using steerable filters. These filters can be designed to optimally respond to a particular type of feature by maximizing several Canny-like criteria. The detection process involves the analytical computation of the orientation and corresponding response of the template. A post-processing step consisting of the suppression of non-maximal values followed by thresholding to eliminate insignificant features concludes the detection procedure. We illustrate the approach with the design of feature templates for the detection of surfaces and curves, and demonstrate their efficiency with practical applications.

We investigate the use of quasi-interpolating approximation schemes, to construct an estimate of an unknown function from its given discrete samples. We show theoretically and with practical experiments that such methods perform better than classical interpolation, for the same computation cost.

Shannon's sampling theory and its variants provide effective solutions to the problem of reconstructing a signal from its samples in some "shift-invariant" space, which may or may not be bandlimited. In this paper, we present some further justification for this type of representation, while addressing the issue of the specification of the best reconstruction space. We consider a realistic setting where a multidimensional signal is prefiltered prior to sampling and the samples corrupted by additive noise. We consider two formulations of the reconstruction problem. In the first deterministic approach, we determine the continuous-space function that minimizes a variational, Tikhonov-like criterion that includes a discrete data term and a suitable continuous-space regularization functional. In the second formulation, we seek the minimum mean square error (MMSE) estimation of the signal assuming that the input signal is a realization of a stationary random process. Interestingly, both approaches yield a solution included in some optimal shift-invariant space that is generally not bandlimited. The solutions can be made equivalent by choosing a regularization operator that corresponds to the whitening filter of the process. We present some practical examples that demonstrate the optimality of the approach.

The dilation matrix associated with the three-dimensional (3-D) face-centered cubic (FCC) sublattice is often considered to be the natural 3-D extension of the two-dimensional (2-D) quincunx dilation matrix. However, we demonstrate that both dilation matrices are of different nature: while the 2-D quincunx matrix is a similarity transform, the 3-D FCC matrix is not. More generally, we show that is impossible to obtain a dilation matrix that is a similarity transform and performs downsampling of the Cartesian lattice by a factor of two in more than two dimensions. Furthermore, we observe that the popular 3-D FCC subsampling scheme alternates between three different lattices: Cartesian, FCC, and quincunx. The latter one provides a less isotropic sampling density, a property that should be taken into account to properly orient 3-D data before processing using such a subsampling matrix.

We propose a novel method for image reconstruction from nonuniform samples with no constraints on their locations. We adopt a variational approach where the reconstruction is formulated as the minimizer of a cost that is a weighted sum of two terms: (1) the sum of squared errors at the specified points and (2) a quadratic functional that penalizes the lack of smoothness. We search for a solution that is a uniform spline and show how it can be determined by solving a large, sparse system of linear equations. We interpret the solution of our approach as an approximation of the analytical solution that involves radial basis functions and demonstrate the computational advantages of our approach. Using the two-scale relation for B-splines, we derive an algebraic relation that links together the linear systems of equations specifying reconstructions at different levels of resolution. We use this relation to develop a fast multigrid algorithm. We demonstrate the effectiveness of our approach on some image reconstruction examples.

We present a new family of two-dimensional and three-dimensional orthogonal wavelets which uses quincunx sampling. The orthogonal refinement filters have a simple analytical expression in the Fourier domain as a function of the order lamda, which may be noninteger. We can also prove that they yield wavelet bases of L2(R2) for any lambda > 0. The wavelets are fractional in the sense that the approximation error at a given scale a decays like O(a(lamda)); they also essentially behave like fractional derivative operators. To make our construction practical, we propose a fast Fourier transform-based implementation that turns out to be surprisingly fast. In fact, our method is almost as efficient as the standard Mallat algorithm for separable wavelets.

This paper deals with fast image and video segmentation using active contours. Region-based active contours using level sets are powerful techniques for video segmentation, but they suffer from large computational cost. A parametric active contour method based on B-Spline interpolation has been proposed in to highly reduce the computational cost, but this method is sensitive to noise. Here, we choose to relax the rigid interpolation constraint in order to robustify our method in the presence of noise: by using smoothing splines, we trade a tunable amount of interpolation error for a smoother spline curve. We show by experiments on natural sequences that this new flexibility yields segmentation results of higher quality at no additional computational cost. Hence, real-time processing for moving objects segmentation is preserved.

We present a new, robust, computational procedure for tracking fluorescent markers in time-lapse microscopy. The algorithm is optimized for finding the time-trajectory of single particles in very noisy dynamic (two- or three-dimensional) image sequences. It proceeds in three steps. First, the images are aligned to compensate for the movement of the biological structure under investigation. Second, the particle's signature is enhanced by applying a Mexican hat filter, which we show to be the optimal detector of a Gaussian-like spot in 1/omega2 noise. Finally, the optimal trajectory of the particle is extracted by applying a dynamic programming optimization procedure. We have used this software, which is implemented as a Java plug-in for the public-domain ImageJ software, to track the movement of chromosomal loci within nuclei of budding yeast cells. Besides reducing trajectory analysis time by several 100-fold, we achieve high reproducibility and accuracy of tracking. The application of the method to yeast chromatin dynamics reveals different classes of constraints on mobility of telomeres, reflecting differences in nuclear envelope association. The generic nature of the software allows application to a variety of similar biological imaging tasks that require the extraction and quantitation of a moving particle's trajectory.

The quantitative assessment of cardiac motion is a fundamental concept to evaluate ventricular malfunction. We present a new optical-flow-based method for estimating heart motion from two-dimensional echocardiographic sequences. To account for typical heart motions, such as contraction/expansion and shear, we analyze the images locally by using a local-affine model for the velocity in space and a linear model in time. The regional motion parameters are estimated in the least-squares sense inside a sliding spatiotemporal B-spline window. Robustness and spatial adaptability is achieved by estimating the model parameters at multiple scales within a coarse-to-fine multiresoluion framework. We use a wavelet-like algorithm for computing B-spline-weighted inner products and moments at dyadic scales to increase computational efficiency. In order to characterize myocardial contractility and to simplify the detection of myocardial dysfunction, the radial component of the velocity with respect to a reference point is color coded and visualized inside a time-varying region of interest. The algorithm was first validated on synthetic data sets that simulate a beating heart with a speckle-like appearance of echocardiograms. The ability to estimate motion from real ultrasound sequences was demonstrated by a rotating phantom experiment. The method was also applied to a set of in vivo echocardiograms from an animal study. Motion estimation results were in good agreement with the expert echocardiographic reading.

In this paper, we use polyharmonic B-splines to build multidimensional wavelet bases. These functions are nonseparable, multidimensional basis functions that are localized versions of radial basis functions. We show that Rabut's elementary polyharmonic B-splines do not converge to a Gaussian as the order parameter increases, as opposed to their separable B-spline counterparts. Therefore, we introduce a more isotropic localization operator that guarantees this convergence, resulting into the isotropic polyharmonic B-splines. Next, we focus on the two-dimensional quincunx subsampling scheme. This configuration is of particular interest for image processing because it yields a finer scale progression than the standard dyadic approach. However, up until now, the design of appropriate filters for the quincunx scheme has mainly been done using the McClellan transform. In our approach, we start from the scaling functions, which are the polyharmonic B-splines and, as such, explicitly known, and we derive a family of polyharmonic spline wavelets corresponding to different flavors of the semi-orthogonal wavelet transform; e.g., orthonormal, B-spline, and dual. The filters are automatically specified by the scaling relations satisfied by these functions. We prove that the isotropic polyharmonic B-spline wavelet converges to a combination of four Gabor atoms, which are well separated in the frequency domain. We also show that these wavelets are nearly isotropic and that they behave as an iterated Laplacian operator at low frequencies. We describe an efficient fast Fourier transform-based implementation of the discrete wavelet transform based on polyharmonic B-splines.

We propose a new spatio-temporal elastic registration algorithm for motion reconstruction from a series of images. The specific application is to estimate displacement fields from two-dimensional ultrasound sequences of the heart. The basic idea is to find a spatio-temporal deformation field that effectively compensates for the motion by minimizing a difference with respect to a reference frame. The key feature of our method is the use of a semi-local spatio-temporal parametric model for the deformation using splines, and the reformulation of the registration task as a global optimization problem. The scale of the spline model controls the smoothness of the displacement field. Our algorithm uses a multiresolution optimization strategy to obtain a higher speed and robustness. We evaluated the accuracy of our algorithm using a synthetic sequence generated with an ultrasound simulation package, together with a realistic cardiac motion model. We compared our new global multiframe approach with a previous method based on pairwise registration of consecutive frames to demonstrate the benefits of introducing temporal consistency. Finally, we applied the algorithm to the regional analysis of the left ventricle. Displacement and strain parameters were evaluated showing significant differences between the normal and pathological segments, thereby illustrating the clinical applicability of our method.

We present a new result characterized by an exact integral expression for the approximation error between a probability density and an integer shift invariant estimate obtained from its samples. Unlike the Parzen window estimate, this estimate avoids recomputing the complete probability density for each new sample - only a few coefficients are required, making it practical for real-time applications. We also show how to obtain the exact asymptotic behavior of the approximation error when the number of samples increases and provide the trade-off between the number of samples and the sampling step size.

Motivated by the fractal-like behavior of natural images, we propose a new smoothing technique that uses a regularization functional which is a fractional iterate of the Laplacian. This type of functional has previously been introduced by Duchon in the context of radial basis functions (RBFs) for the approximation of non-uniform data. Here, we introduce a new solution to Duchon's smoothing problem in multiple dimensions using non-separable fractional polyharmonic B-splines. The smoothing is performed in the Fourier domain by filtering, thereby making the algorithm fast enough for most multi-dimensional real-time applications.

We propose the use of polyharmonic B-splines to build non-separable two-dimensional wavelet bases. The central idea is to base our design on isotropic-polyharmonic B-splines, a new type of polyharmonic B-splines that converge to a Gaussian as the order increases. We opt for the quincunx subsampling scheme which allows us to characterize the wavelet spaces with a single wavelet, the isotropic-polyharmonic B-spline wavelet. Interestingly, this wavelet converges to a combination of four Gabor atoms, which are well separated in the frequency domain. We also briefly discuss our Fourier-based implementation and present some experimental results.

We propose a general approach for the design of 2D feature detectors from a class of steerable functions based on the optimization of a Canny-like criterion. In contrast with previous computational designs, our approach is truly 2D and provides filters that have closed-form expressions. It also yields operators that have a better orientation selectivity than the classical gradient or Hessian-based detectors. We illustrate the method with the design of operators for edge and ridge detection. We present some experimental results that demonstrate the performance improvement of these new feature detectors. We propose computationally efficient local optimization algorithms for the estimation of feature orientation. We also introduce the notion of shape-adaptable feature detection and use it for the detection of image corners.

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.

Parametric active contour models are one of the preferred approaches for image segmentation because of their computational efficiency and simplicity. However, they have a few drawbacks which limit their performance. In this paper, we identify some of these problems and propose efficient solutions to get around them. The widely-used gradient magnitude-based energy is parameter dependent; its use will negatively affect the parametrization of the curve and, consequently, its stiffness. Hence, we introduce a new edge-based energy that is independent of the parameterization. It is also more robust since it takes into account the gradient direction as well. We express this energy term as a surface integral, thus unifying it naturally with the region-based schemes. The unified framework enables the user to tune the image energy to the application at hand. We show that parametric snakes can guarantee low curvature curves, but only if they are described in the curvilinear abscissa. Since normal curve evolution do not ensure constant arc-length, we propose a new internal energy term that will force this configuration. The curve evolution can sometimes give rise to closed loops in the contour, which will adversely interfere with the optimization algorithm. We propose a curve evolution scheme that prevents this condition.

We introduce local weighted geometric moments that are computed from an image within a sliding window at multiple scales. When the window function satisfies a two-scale relation, we prove that lower order moments can be computed efficiently at dyadic scales by using a multiresolution wavelet-like algorithm. We show that B-splines are well-suited window functions because, in addition to being refinable, they are positive, symmetric, separable, and very nearly isotropic (Gaussian shape). We present three applications of these multiscale local moments. The first is a feature-extraction method for detecting and characterizing elongated structures in images. The second is a noise-reduction method which can be viewed as a multiscale extension of Savitzky-Golay filtering. The third is a multiscale optical-flow algorithm that uses a local affine model for the motion field, extending the Lucas-Kanade optical-flow method. The results obtained in all cases are promising.

This paper proposes a new family of bivariate, nonseparable splines, called hex-splines, especially designed for hexagonal lattices. The starting point of the construction is the indicator function of the Voronoi cell, which is used to define in a natural way the first-order hex-spline. Higher order hex-splines are obtained by successive convolutions. A mathematical analysis of this new bivariate spline family is presented. In particular, we derive a closed form for a hex-spline of arbitrary order. We also discuss important properties, such as their Fourier transform and the fact they form a Riesz basis. We also highlight the approximation order. For conventional rectangular lattices, hex-splines revert to classical separable tensor-product B-splines. Finally, some prototypical applications and experimental results demonstrate the usefulness of hex-splines for handling hexagonally sampled data.

We describe a new family of scaling functions, the (/spl alpha/, /spl tau/)-fractional splines, which generate valid multiresolution analyses. These functions are characterized by two real parameters: /spl alpha/, which controls the width of the scaling functions; and /spl tau/, which specifies their position with respect to the grid (shift parameter). This new family is complete in the sense that it is closed under convolutions and correlations. We give the explicit time and Fourier domain expressions of these fractional splines. We prove that the family is closed under generalized fractional differentiations, and, in particular, under the Hilbert transformation. We also show that the associated wavelets are able to whiten 1/f/sup /spl lambda//-type noise, by an adequate tuning of the spline parameters. A fast (and exact) FFT-based implementation of the fractional spline wavelet transform is already available. We show that fractional integration operators can be expressed as the composition of an analysis and a synthesis iterated filterbank.

This paper presents an interpolation method based on shifted versions of two piecewise linear generators, which provides approximation order 2 like usual piecewise-linear interpolation; i.e., this method is able to represent the constant and the ramp exactly. Our interpolation is characterized by two real parameters: /spl tau/, the location of the generators, and /spl alpha/, related to their dissymmetry. By varying these parameters, we show that it is possible to optimize the quality of the approximation, independently of the function to interpolate. We recover the optimal value of shifted-linear interpolation (/spl tau/ = 0.21 and /spl alpha/ = 1) which requires IIR prefiltering, but we also find a new configuration (/spl tau/ = 0.21 and /spl alpha/ = 0.58) which reaches almost the same quality, while requiring FIR filtering only. This new solution is able to greatly reduce the amount of Gibbs oscillations generated in the shifted-linear interpolation scheme. We validate our finding by computing the PSNR of the difference between multi-rotated images and their original version.

Hex-splines are a novel family of bivariate splines which are well suited to handle hexagonally sampled data. Similar to classical 1D B-splines, the spline coefficients need to be computed by a prefilter. Unfortunately, the elegant implementation of this prefilter by causal and anti-causal recursive filtering is not applicable for the (non-separable) hex-splines. Therefore, in this paper we introduce a novel approach from the viewpoint of approximation theory. We propose three different recursive filters and optimize their parameters such that a desired order of approximation is obtained. The results for third and fourth order hex-splines are discussed. Although the proposed solutions provide only quasi-interpolation, they tend to be very close to the interpolation prefilter.

We present a new approach for the design of optimal steerable 2-D templates for feature detection. As opposed to classical schemes where the optimal 1-D template is derived and extended to 2-D, we directly obtain the 2-D template. We choose the template from a class of steerable functions based on the analytic optimization of a Canny-like criterion. Our approach gives more orientation selective templates that have simple closed form expression. We illustrate the method with the design of operators for edge and ridge detection and demonstrate their performance improvement in practical applications.

This paper deals with fast image and video segmentation using active contours. Region based active contours using level-sets are powerful techniques for video segmentation hut they suffer from large computational cost. A parametric active contour method based on B-Spline interpolation has been proposed in F. Precioso (2002) to highly reduce the computational cost but this method is sensitive to noise. Here, we choose to relax the rigid interpolation constraint in order to robustify our method in the presence of noise: by using smoothing splines, we trade a tunable amount of interpolation error for a smoother spline curve. We show by experiments on natural sequences that this new flexibility yields segmentation results of higher quality at no additional computational cost. Hence real time processing for moving objects segmentation is preserved.

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 present an algorithm for fast elastic multidimensional intensity-based image registration with a parametric model of the deformation. It is fully automatic in its default mode of operation. In the case of hard real-world problems, it is capable of accepting expert hints in the form of soft landmark constraints. Much fewer landmarks are needed and the results are far superior compared to pure landmark registration. Particular attention has been paid to the factors influencing the speed of this algorithm. The B-spline deformation model is shown to be computationally more efficient than other alternatives. The algorithm has been successfully used for several two-dimensional (2-D) and three-dimensional (3-D) registration tasks in the medical domain, involving MRI, SPECT, CT, and ultrasound image modalities. We also present experiments in a controlled environment, permitting an exact evaluation of the registration accuracy. Test deformations are generated automatically using a random hierarchical fractional wavelet-based generator.

We propose a construction of new wavelet-like bases that are well suited for the reconstruction and processing of optically generated Fresnel holograms recorded on CCD-arrays. The starting point is a wavelet basis of L2 to which we apply a unitary Fresnel transform. The transformed basis functions are shift-invariant on a level-by-level basis but their multiresolution properties are governed by the special form that the dilation operator takes in the Fresnel domain. We derive a Heisenberg-like uncertainty relation that relates the localization of Fresnelets with that of their associated wavelet basis. According to this criterion, the optimal functions for digital hologram processing turn out to be Gabor functions, bringing together two separate aspects of the holography inventor's work. We give the explicit expression of orthogonal and semi-orthogonal Fresnelet bases corresponding to polynomial spline wavelets. This special choice of Fresnelets is motivated by their near-optimal localization properties and their approximation characteristics. We then present an efficient multiresolution Fresnel transform algorithm, the Fresnelet transform. This algorithm allows for the reconstruction (backpropagation) of complex scalar waves at several user-defined, wavelength-independent resolutions. Furthermore, when reconstructing numerical holograms, the subband decomposition of the Fresnelet transform naturally separates the image to reconstruct from the unwanted zero-order and twin image terms. This greatly facilitates their suppression. We show results of experiments carried out on both synthetic (simulated) data sets as well as on digitally acquired holograms.

We establish a link between classical osculatory interpolation and modern convolution-based interpolation and use it to show that two well-known cubic convolution schemes are formally equivalent to two osculatory interpolation schemes proposed in the actuarial literature about a century ago. We also discuss computational differences and give examples of other cubic interpolation schemes not previously studied in signal and image processing.

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.

The LeGall 5/3 and Daubechies 9/7 filters have risen to special prominence because they were selected for inclusion in the JPEG2000 standard. We determine their key mathematical features: Riesz bounds, order of approximation, and regularity (Hölder and Sobolev). We give approximation theoretic quantities such as the asymptotic constant for the L2 error and the angle between the analysis and synthesis spaces which characterizes the loss of performance with respect to an orthogonal projection. We also derive new asymptotic error formulae that exhibit bound constants that are proportional to the magnitude of the first nonvanishing moment of the wavelet. The Daubechies 9/7 stands out because it is very close to orthonormal, but this turns out to be slightly detrimental to its asymptotic performance when compared to other wavelets with four vanishing moments.

I. Introduction Wavelets are the result of collective efforts that recognized common threads between ideas and concepts that had been independently developed and investigated by distinct research communities. They provide a unifying framework for decomposing images, volumes, and time-series data into their elementary constituents across scale. Although a relatively recent construct, wavelets have become a tool of choice for engineers, physicists, and mathematicians, leading to efficient solutions in time and space frequency analysis problems, as well as a multitude of other applications. One of the consequences is that wavelet methods of analysis and representation are presently having a significant impact on the science of medical imaging and the diagnosis of disease and screening protocols. Because of a powerful underlying mathematical theory, they offer exciting opportunities for the design of new multiresolution image processing algorithms, and novel acquisition methods such as wavelet-encoded magnetic resonance imaging (MRI). This special issue of the IEEE Transactions on Medical Imaging focuses on these recent developments and highlights progress that has been accomplished in the areas related to medical imaging. II. Sizing the Wave: Some Facts and Figures Rarely has a mathematical concept generated so much response and enthusiasm within and between the engineering and mathematical research communities at large. To give a rough idea of the phenomenon, we provide a brief chronology. While wavelets have been traced all the way back to Alfred Haar in 1910 [1], for many, the starting point of their modern history coincides with two publications in the late 1980s by S. Mallat [2] and I. Daubechies [3]. These groundbreaking papers established a solid mathematical footing which would both shape and define the field. In a nutshell, S. Mallat identified the important concept of multiresolution analysis which is the corner stone of modern wavelet theory, while I. Daubechies constructed the first orthogonal wavelet bases that were compactly supported. These two contributions count among the most cited papers in the scientific literature (over 1500 SCI citations each). From that point on, the number of contributions relating to wavelet-applications and theory has increased steadily on the order of 9000 journal papers published to date. This trend is likely to continue as suggested by the strong response to the call for papers for this special issue (over 30 submissions). Wavelets have become so popular that distinct communities continue to have conferences and scientific journals entirely devoted to them. There are already historical anecdotes and folklore associated with them; an entertaining account of which can be found in the book of B. Burke Hubbard [4]. Readers who want to dive deeper into the subject have the daunting task of choosing among over 200 books written on wavelets. Our only advice in this regard is: in case of doubt, stick with the classics. Given the size of the phenomenon, it is no surprise that wavelets have had an impact on a number of disciplines, medical imaging being no exception. A first record of activity in this particular area is the workshop on wavelets in medicine and biology that took place at the annual IEEE-EMBS meeting, Baltimore, MD, 1992. The first journal paper describing a wavelet application in medical imaging—noise reduction in MRI by soft-thresholding in the wavelet domain—also appeared in 1992 [5]. Note that this work, which is often overlooked, provided the earliest description of a wavelet denoising method that has become extremely popular through the impulsion of Donoho et al. A large palette of wavelet applications in medical imaging is provided in [6]. Two complementary review articles are also available; the first gives a complete account of the activity taking place from the beginning to 1996 [7], while the second covers the more recent papers until 2000 [8]. So far, the primary applications of wavelets in medical imaging have been the following: Compression of medical images. CT reconstruction; local tomography. Wavelet denoising (MRI, ultrasound). Wavelet-based feature extraction; texture and statistical descriptors Medical image enhancement (e.g. fluoroscopy, mammography). Analysis of functional images of the brain [positron emission tomography (PET), functional MRI (fMRI)]. Wavelet-encoded MRI: Most of these topics are still active areas of research, as illustrated by the papers that are published in this special issue. III. Scanning Through the Issue Functional imaging is an area were methods of wavelet processing hold great promise. This particular line of research was initiated by U. Ruttimann, a creative researcher and good friend, who sadly passed away shortly before the publication of his paper in this very journal [9]. Another first rate statistician who was also active in this area at an early stage is J. Raz. By a sad coincidence, he also suffered a sudden death about a week before he was to present his latest results on wavelet analysis of fMRI [10]. Despite the tragic loss of these two pioneers, research in this area is alive and well as exemplified by the first three papers of this issue. Turkheimer et al. [11] consider the problem of the analysis of dynamic PET data; in particular, they advocate the use of a linear (James-Stein) wavelet estimator as an alternative to the more classical wavelet shrinkage or threshold detectors. Hossein-Zadeh et al. propose a wavelet-technique for the detection of activation in fMRI data [12]. Their contribution is twofold: first, the use of a redundant wavelet transform for better translation invariance, and second, a nonparametric detection method based on a randomization procedure. F. Meyer also considers fMRI time series but applies wavelets differently, within the context of a generalized linear model, to detrend the data [13]; that is, to get rid of signal drifts and disturbances that are not related to the stimulus. Two other areas where wavelets have achieved great success is signal denoising (typically, by simple thresholding in the wavelet domain) and tomographic reconstruction, mainly, because the Radon operator is well localized in a wavelet basis. Pizurica et al. [14] propose a novel wavelet denoising method that uses a local statistical model for improved signal estimation and noise suppression. Willett and Nowak [15] introduce a new multiscale image model, using piecewise planar basis functions, and apply their method to the reconstruction of photon-limited data (with Poisson noise). They develop penalized maximum likelihood methods for image denoising, deconvolution, and tomographic reconstruction. Kalifa et al. present a direct method for the efficient reconstruction of PET and single photom emission computed tomography data [16]. Their approach includes a nonlinear noise reduction step that is implemented by thresholding in the transformed domain; the key here is to select a transform (wavelet packet) that is optimized for the problem and data at hand (sparse representation of the signal and near diagonalization of the Radon operator). Bonnet et al. [17] also develop a direct approach for the reconstruction of cone-beam data which is known to be challenging. In essence, their approach is a wavelet adaptation of the Feldkamp algorithm. The special issue also features two contributions relating to ultrasound imaging. Michailovich and Adam consider the problem of the estimation of the spectrum of a ultrasound pulse [18]. Specifically, they develop a modified (outlier-resistant) wavelet estimator that they apply to the log-spectrum of the radio-freqency sequence. Lee et al. present a pattern recognition system that uses an M-band wavelet filterbank to extract fractal and texture features from ultrasonic images of the liver [19]. They report promising classification results, differientiating normal liver, cirrhosis, and hepatoma using a hierarchical classifier. The recent development of commercial digital mammography imaging systems not only provides a significant improvement in image quality for traditional screening, but translates into a wealth of information for the analysis and detection of mammographic features by computer. The papers by Lemaur et al. [20]. and Heinlein et al. [21] focus on the goals of early detection and visual enhancement of microcalcifications, respectively. In the former, the regularity of a wavelet basis is used to identify microcalcification in clusters. The identification of microcalcifaction in clusters as opposed to individual occurrences is of clinical significance as clusters may suggest the likelihood of malignancy. In the later paper, a discretization of the continuous wavelet transform is developed which allows a filterbank to be adapted for the enhancement of mammographic features. This implementation allows for the reconstruction of modified wavelet coefficients at arbitrary scales and orientations without the introduction of artifacts or loss of completeness. The integration of such an interactive enhancement tool into digital mammographic screening systems will be of great importance as the wealth of dynamic range (contrast) provided by digital detectors become generally available to radiologist through the introduction of lower cost softcopy display systems. The paper of Davatzikos et al. [22] offers another illustration of the versatility of wavelets. Their proposal is to represent the contours of a shape in a wavelet bases and to use this sparse representation to derive active shape models. Their results are promising and significant in terms of providing an automated solution to problems in volume quantification. The amount of data generated by modern imaging devices is often very large, and ever increasing. Thus, an important problem is to find efficient ways of compressing and encoding this information

We propose a novel multiresolution-multigrid based signal reconstruction method from arbitrarily spaced samples. The signal is reconstructed on a uniform grid using B-splines basis functions. The computation of spline weights is formulated as a variational problem. Specifically, we minimize a cost that is a weighted sum of two terms: (i) the sum of squared errors at the specified points; (ii) a quadratic functional that penalizes the lack of smoothness. The problem is equivalent to solving a very large system of linear equations, with the dimension equal to the number of grid points. We develop a computationally efficient multiresolution-multigrid scheme for solving the system. We demonstrate the method with image reconstruction from contour points.

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 define multiscale moments that are estimated locally by analyzing an image through a sliding window at multiple scales. When the analysis window satisfies a two-scale relation, we prove that these moments can be computed very efficiently using a multiresolution wavelet-like algorithm. We also show that B-spline windows are best suited for this kind of analysis because, in addition to being refinable, they are positive, symmetric and very nearly isotropic. We present two applications of our method. The first is a feature extraction method for detecting strands of DNA in noisy cryoelectron-micrographs. The second is an extension of the Lucas-Kanade optical flow algorithm (see Lucas, B. and Kanade, T., Proc. DARPA IU Workshop, p.121-30, 1981) that assumes a local affine model of the motion field. The results obtained in both cases are very promising.

We propose to design the reduction operator of an image pyramid so as to minimize the approximation error in the l(p)-sense (not restricted to the usual p=2), where p can take noninteger values. The underlying image model is specified using shift-invariant basis functions, such as B-splines. The solution is well-defined and determined by an iterative optimization algorithm based on digital filtering. Its convergence is accelerated by the use of first and second order derivatives. For p close to 1, we show that the ringing is reduced and that the histogram of the detail image is sparse as compared with the standard case, where p=2.

We present an explicit formula for B-spline convolution kernels; these are defined as the convolution of several B-splines of variable widths h(i) and degrees n(i). We apply our results to derive spline-convolution-based algorithms for two closely related problems: the computation of the Radon transform and of its inverse. First, we present an efficient discrete implementation of the Radon transform that is optimal in the least-squares sense. We then consider the reverse problem and introduce a new spline-convolution version of the filtered back-projection algorithm for tomographic reconstruction. In both cases, our explicit kernel formula allows for the use of high-degree splines; these offer better approximation performance than the conventional lower-degree formulations (e.g., piecewise constant or piecewise linear models). We present multiple experiments to validate our approach and to find the parameters that give the best tradeoff between image quality and computational complexity. In particular, we find that it can be computationally more efficient to increase the approximation degree than to increase the sampling rate.

Guest Editorial Technological advances in biomedical imaging are providing unprecedented opportunities for improving our understanding of biological processes and revealing the anatomical and functional organization of biological systems, from macro- to nano-scales. Research in this area is progressing at an extraordinary rate and is becoming more and more interdisciplinary. The IEEE International Symposium on Biomedical Imaging (ISBI) was initiated to bring together researchers from the medical and biological imaging communities and to provide an effective forum for multidisciplinary interactions. The first meeting was held at the Ritz-Carlton Hotel, Washington, DC, July 7-10, 2002. It was organized jointly by the IEEE Signal Processing Society (SPS) and the IEEE Engineering in Medicine and Biology Society (EMBS) and co-sponsored by National Institutes of Health's (NIH) National Institute of the Biomedical Imaging and Bioengineering (NIBIB), representing the first IEEE-NIH collaborative effort on a major imaging conference. The symposium was focused on the engineering aspects of biomedical imaging while promoting an integrative approach through all scales of observation. It successfully brought together a large group of biomedical imaging researchers and practitioners (535 participants) with different backgrounds to share their knowledge and to address the latest challenges in data acquisition, image reconstruction, image processing, analysis, and visualization. A highlight of the symposium was the inspired opening address by Dr. Elias Zerhouni, the newly appointed director of the NIH. He stressed the importance of imaging for the biomedical sciences and expressed a strong interest in the conference. Dr. Zerhouni is very much in favor of collaborations between engineers and biomedical scientists, having experienced them first hand as a radiologist. He holds several patents related to imaging and is widely known for having introduced magnetic resonance tagging as a diagnostic tool for assessing heart function. Another noted NIH speaker was Dr. Roderic Pettigrew, who will become NIBIB's first permanent director in September 2002. He opened the first session by expressing his strong support for ISBI; he also explained the mission of NIBIB in which imaging will play a major role. The scientific program of the conference consisted of 3 plenary talks, 10 special sessions, 19 oral sessions, and 6 poster sessions over three days. The plenary talks, one on each day, covered topics from molecular imaging (by Michael Phelps, University of California, Los Angeles (UCLA), School of Medicine), functional imaging (by Alan Koretsky, NIH), to medical image analysis (by Michael Brady, University of Oxford, U.K.). The first two talks reviewed and discussed the latest developments in imaging modalities and applications, while the latter illustrated the power and potential of image processing and analysis algorithms. The special sessions, consisting of invited papers and sponsored by NIBIB, were a unique feature of the conference. They covered a wide range of topics: "Model-Based Image Segmentation and Analysis" (organized by Christos Davatzikos), "Microarray Image Processing and Analysis" (organized by Bin Yu and Dan Bartell), "Optical Coherence Tomography" (organized by Stephen Boppart and René Salathé), "Micro Imaging" (organized by Erik Ritman and Françoise Peyrin), "Brain Connectivity and Functional Assessment" (organized by Peter Basser and John George), "Nonrigid Registration" (organized by Benoit Dawant), "Fast Acquisition and Sampling in MRI" (organized by Yoram Bresler), "Electron Microscopy" (organized by Benes Trus), "In Vivo Cellular and Molecular Imaging" (organized by John Hoffman and Gary Kelloff), and "Federal Funding Opportunities for Biomedical Imaging Research" (organized by Richard Swaja, NIH/NIBIB). In addition to the invited program, a total of 355 papers were submitted to the symposium, of which 73 were accepted for oral presentation and 142 for poster presentation. The acceptance rate for contributed papers was about 60%. The papers of the symposium are published in IEEE conference proceedings; they will be included in the IEEEXplore database which is searchable through the WEB. There are many indications that ISBI fulfills a strong need of the biomedical imaging community. Many participants expressed enthusiastic support for the conference. It was felt that the time was ripe for IEEE to have its flagship conference on biomedical imaging. In charge of planning the next meeting are Christian Roux (representing EMBS) and Richard Leahy (representing SPS). It is our belief and hope that the next ISBI, to be held in Washington, DC, once again, will be even more successful in serving our community. Efforts will be made to increase the participation and the level of cross fertilization between the medical and biological imaging communities. ISBI may also play an important role in graduate education and training, as it is often impossible for a single institution to develop a comprehensive biomedical imaging curriculum covering all the modalities and applications. It is, therefore, desirable that future ISBI's include a one-day pre-conference event with tutorials and short courses covering not only imaging fundamentals, but also emerging technologies and applications. In closing, we may say that ISBI was a success, partly due to timing, location, and the fact that the conference responded to a need of the engineering community. However, nothing could have happened without the hard work, creativity, and involvement of the volunteers. We would like to take this opportunity to thank the entire organizing committee, the program committee, the special session chairs, and the external reviewers for their large commitment of time and effort. All administrative matters were handled by the staff of both SPS and EMBS under the leadership of their executive directors, Mercy Kowalczyk and Laura Wolf. The technical program chairs were Jeffrey Fessler (for SPS) and Michael Vannier (for EMBS), with additional support provided by the special sessions chairs Erik Meijering and Jean-Louis Coatrieux and the NIBIB representative on the committee, Richard Swaja, who played a crucial role in involving the NIH. We do not know how to thank them enough for having helped us to turn what initially looked like a nice idea into a reality.

Wavelet-based methods have become most popular for the compression of two-dimensional medical images and sequences. The standard implementations consider data sizes that are powers of two. There is also a large body of literature treating issues such as the choice of the "optimal" wavelets and the performance comparison of competing algorithms. With the advent of telemedicine, there is a strong incentive to extend these techniques to higher dimensional data such as dynamic three-dimensional (3-D) echocardiography [four-dimensional (4-D) datasets]. One of the practical difficulties is that the size of this data is often not a multiple of a power of two, which can lead to increased computational complexity and impaired compression power. Our contribution in this paper is to present a genuine 4-D extension of the well-known zerotree algorithm for arbitrarily sized data. The key component of our method is a one-dimensional wavelet algorithm that can handle arbitrarily sized input signals. The method uses a pair of symmetric/antisymmetric wavelets (10/6) together with some appropriate midpoint symmetry boundary conditions that reduce border artifacts. The zerotree structure is also adapted so that it can accommodate noneven data splitting. We have applied our method to the compression of real 3-D dynamic sequences from clinical cardiac ultrasound examinations. Our new algorithm compares very favorably with other more ad hoc adaptations (image extension and tiling) of the standard powers-of-two methods, in terms of both compression performance and computational cost. It is vastly superior to slice-by-slice wavelet encoding. This was seen not only in numerical image quality parameters but also in expert ratings, where significant improvement using the new approach could be documented. Our validation experiments show that one can safely compress 4-D data sets at ratios of 128:1 without compromising the diagnostic value of the images. We also display some more extreme compression results at ratios of 2000:1 where some key diagnostically relevant key features are preserved.

We present a new family of 2D orthogonal wavelets which use quincunx sampling. The orthogonal refinement filters have a simple analytical expression in the Fourier domain as a function of the order /spl alpha/, which may be non-integer. The wavelets have good isotropy properties. We can also prove that they yield wavelet bases of L/sub 2/(R/sup 2/) for any /spl alpha/>0. The wavelets are fractional in the sense that the approximation error at a given scale /spl alpha/ decays like O(a/sup /spl alpha//); they also essentially behave like fractional derivative operators. To make our construction practical, we propose an FFT-based implementation that turns out to be surprisingly fast. In fact, our method is almost as efficient as the standard Mallat algorithm for separable wavelets.

We have designed a series of computer sessions build around ImageJ (a public-domain software for image analysis), as a practical complement to a two-semester course in image processing. The students are challenged with simple practical imaging problems as they acquire hands-on practice by experimenting with image-processing operators. In the process, they also learn how to program standard image-processing algorithms in Java. This is made possible thanks to a programmer-friendly environment and a software interface that greatly facilitates the developments of plug-ins for ImageJ. Since our students have generally not acquired programming skills yet (they typically do not even know Java), we use a learning-by-example teaching strategy, with good success.

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 present a method for the exact computation of the moments of a region bounded by a curve represented by a scaling function or wavelet basis. Using Green's theorem, we show that the computation of the area moments is equivalent to applying a suitable multidimensional filter on the coefficients of the curve and thereafter computing a scalar product. The multidimensional filter coefficients are precomputed exactly as the solution of a two-scale relation. To demonstrate the performance improvement of the new method, we compare it with existing methods such as pixel-based approaches and approximation of the region by a polygon. We also propose an alternate scheme when the scaling function is sinc(x).

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 present an optimal spline-based algorithm for the enlargement or reduction of digital images with arbitrary (noninteger) scaling factors. This projection-based approach can be realized thanks to a new finite difference method that allows the computation of inner products with analysis functions that are B-splines of any degree n. A noteworthy property of the algorithm is that the computational complexity per pixel does not depend on the scaling factor a. For a given choice of basis functions, the results of our method are consistently better than those of the standard interpolation procedure; the present scheme achieves a reduction of artifacts such as aliasing and blocking and a significant improvement of the signal-to-noise ratio. The method can be generalized to include other classes of piecewise polynomial functions, expressed as linear combinations of B-splines and their derivatives.

We define a new wavelet transform that is based on a previously defined family of scaling functions: the fractional B-splines. The interest of this family is that they interpolate between the integer degrees of polynomial B-splines and that they allow a fractional order of approximation. The orthogonal fractional spline wavelets essentially behave as fractional differentiators. This property seems promising for the analysis of 1/f/sup /spl alpha// noise that can be whitened by an appropriate choice of the degree of the spline transform. We present a practical FFT-based algorithm for the implementation of these fractional wavelet transforms, and give some examples of processing.

We present an explicit formula for spline kernels; these are defined as the convolution of several B-splines of variable widths h/sub i/ and degrees n/sub i/. The spline kernels are useful for continuous signal processing algorithms that involve B-spline inner-products or the convolution of several spline basis functions. We apply our results to the derivation of spline-based algorithms for two classes of problems. The first is the resizing of images with arbitrary scaling factors. The second is the computation of the Radon transform and of its inverse. In particular, we present a new spline-based version of the filtered back-projection algorithm for tomographic reconstruction. In both cases, our explicit kernel formula allows for the use of high-degree splines; these offer better approximation performance than the conventional lower-order formulations (e.g., piecewise constant or piece-wise linear models).

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.

We present an algorithm for nonlinear multidimensional registration. The correspondence function is represented in a spline space. We also use cubic splines to interpolate the volumetric data to be registered. We use a multiresolution strategy, which gives us robustness and additional speedup. We analyze the computational complexity of the algorithm and its performance using different multidimensional optimization methods. Finally, we present an application of the algorithm for registering technetium ethylene cysteine diethylester (ECD) SPECT images for realignment of corresponding Xe inhalation SPECT images.

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 present an exact algorithm for the computation of the moments of a region bounded by a curve represented in a scaling function or wavelet basis. Using Green's theorem, we show that the computation of the area moments is equivalent to applying a suitable multidimensional filter on the coefficients of the curve and thereafter computing a scalar product. We compare this algorithm with existing methods such as pixel-based approached and approximation of the region by a polygon.

This paper presents an account of the current state of sampling, 50 years after Shannon's formulation of the sampling theorem. The emphasis is on regular sampling, where the grid is uniform. This topic has benefitted from a strong research revival during the past few years, thanks in part to the mathematical connections that were made with wavelet theory. To introduce the reader to the modern, Hilbert-space formulation, we reinterpret Shannon's sampling procedure as an orthogonal projection onto the subspace of band-limited functions. We then extend the standard sampling paradigm for a presentation of functions in the more general class of "shift-in-variant" function spaces, including splines and wavelets. Practically, this allows for simpler-and possibly more realistic-interpolation models, which can be used in conjunction with a much wider class of (anti-aliasing) prefilters that are not necessarily ideal low-pass. We summarize and discuss the results available for the determination of the approximation error and of the sampling rate when the input of the system is essentially arbitrary; e.g., nonbandlimited. We also review variations of sampling that can be understood from the same unifying perspective. These include wavelets, multiwavelets, Papoulis generalized sampling, finite elements, and frames. Irregular sampling and radial basis functions are briefly mentioned.

We present a novel formulation for B-spline snakes that can be used as a tool for fast and intuitive contour outlining. We start with a theoretical argument in favor of splines in the traditional formulation by showing that the optimal, curvature-constrained snake is a cubic spline, irrespective of the form of the external energy field. Unfortunately, such regularized snakes suffer from slow convergence speed because of a large number of control points, as well as from difficulties in determining the weight factors associated to the internal energies of the curve. We therefore propose an alternative formulation in which the intrinsic scale of the spline model is adjusted a priori; this leads to a reduction of the number of parameters to be optimized and eliminates the need for internal energies (i.e., the regularization term). In other words, we are now controlling the elasticity of the spline implicitly and rather intuitively by varying the spacing between the spline knots. The theory is embedded into a multiresolution formulation demonstrating improved stability in noisy image environments. Validation results are presented, comparing the traditional snake using internal energies and the proposed approach without internal energies, showing the similar performance of the latter. Several biomedical examples of applications are included to illustrate the versatility of the method.

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 present an optimal spline-based algorithm for the enlargement or reduction of digital images with arbitrary scaling factors. This projection-based approach is realizable thanks to a new finite difference method that allows the computation of inner products with analysis functions that are B-splines of any degree n. For a given choice of basis functions, the results of our method are consistently better that those of the standard interpolation procedure; the present scheme achieves a reduction of artifacts such as aliasing and blocking and a significant improvement of the signal-to-noise ratio.

Quadtree-like pyramids have the advantage of re-suiting in a multiresolution representation where each pyramid node has four unambiguous parents. Such a centered topology guarantees a clearly defined up-projection of labels. This concept has been successfully and extensively used in applications of contour detection, object recognition and segmentation. Unfortunately, the quadtree-like type of pyramid has poor approximation powers because of the employed piecewise-constant image model. This paper deals with the construction of improved centered image pyramids in terms of general approximation functions. The advantages of the centered topology such a symmetry, consistent boundary conditions and accurate up-projection of labels are combined with a more faithful image representation at coarser pyramid levels. We start by introducing a general framework for the design of least squares pyramids using the standard filtering and decimation tools. We give the most general explicit formulas for the computation of the filter coefficients by any (well behaving) approximation function in both the continuous (L(2)) and the discrete (l(2)) norm. We then define centered pyramids and provide the filter coefficients for odd spline approximation functions. Finally, we compare the centered pyramid to the ordinary one and highlight some applications.

We introduce a method for multichannel restoration of images in which there is severely limited knowledge about the undegraded signal, and possibly the noise. We assume that we know the channel degradations and that there will be a significant noise reduction in a postprocessing stage in which multiple realizations are combined. This post-restoration noise reduction is often performed when working with micrographs of biological macromolecules. The restoration filters are designed to enforce a projection constraint upon the entire system. This projection constraint results in a system that provides an oblique projection of the input signal into the subspace defined by the reconstruction device in a direction orthogonal to a space defined by the channel degradations and the restoration filters. The approach achieves noise reduction without distorting the signal by exploiting the redundancy of the measurements.

We present a novel formulation for B-spline snakes that can be used as a tool for fast and intuitive contour outlining. The theory is implemented in a platform independent JAVA interface, which allows real time computation of the snake curve. In this paper, our main focus is on two points. First, we propose a novel B-spline snake formulation, where the intrinsic scale of the spline model is adjusted a priori. It leads to a reduction of the number of parameters to be optimized and eliminates the need for internal energies. The approach solves the two main drawbacks of traditional snakes (slow convergence speed and difficult to adjust weighting factors for internal energy terms). Second, we comment on our experience using JAVA for the implementation of the snake and for the design of the graphical user interface. Our technique provides a very intuitive, user-friendly, and platform independent tool for contour outlining, generally applicable to a vast range of images. Several biomedical examples of applications are included to illustrate the versatility of the 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.

The standard interpolation approach to image resizing is to fit the original picture with a continuous model and resample the function at the desired rate. However, one can obtain more accurate results if one applies a filter prior to sampling, a fact well known from sampling theory. The optimal solution corresponds to an orthogonal projection onto the underlying continuous signal space. Unfortunately, the optimal projection prefilter is difficult to implement when sine or high order spline functions are used. We propose to resize the image using an oblique rather than an orthogonal projection operator in order to make use of faster, simpler, and more general algorithms. We show that we can achieve almost the same result as with the orthogonal projection provided that we use the same approximation space. The main advantage is that it becomes perfectly feasible to use higher order models (e.g. splines of degree n=or>3). We develop the theoretical background and present a simple and practical implementation procedure using B-splines. Our experiments show that the proposed algorithm consistently outperforms the standard interpolation methods and that it provides essentially the same performance as the optimal procedure (least squares solution) with considerably fewer computations. The method works for arbitrary scaling factors and is applicable to both image enlargement and reduction.

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.

The use of the wavelet transform is explored for the detection of differences between brain functional magnetic resonance images (fMRI's) acquired under two different experimental conditions. The method benefits from the fact that a smooth and spatially localized signal can be represented by a small set of localized wavelet coefficients, while the power of white noise is uniformly spread throughout the wavelet space. Hence, a statistical procedure is developed that uses the imposed decomposition orthogonality to locate wavelet-space partitions with large signal-to-noise ratio (SNR), and subsequently restricts the testing for significant wavelet coefficients to these partitions. This results in a higher SNR and a smaller number of statistical tests, yielding a lower detection threshold compared to spatial-domain testing and, thus, a higher detection sensitivity without increasing type I errors. The multiresolution approach of the wavelet method is particularly suited to applications where the signal bandwidth and/or the characteristics of an imaging modality cannot be well specified. The proposed method was applied to compare two different fMRI acquisition modalities. Differences of the respective useful signal bandwidths could be clearly demonstrated; the estimated signal, due to the smoothness of the wavelet representation, yielded more compact regions of neuroactivity than standard spatial-domain testing.

We investigate the problem of the reconstruction of a continuous-time function f(x)/spl isin//spl Hscr/ from the responses of m linear shift-invariant systems sampled at 1/m the reconstruction rate, extending Papoulis' (1977) generalized sampling theory in two important respects. First, we allow for arbitrary (non-bandlimited) input signals (typ. /spl Hscr/=L/sub 2/). Second, we use a more general specification of the reconstruction subspace V(/spl phi/), so that the output of the system can take the form of a bandlimited function, a spline, or a wavelet expansion. The system that we describe yields an approximation f/spl isin/V(/spl phi/) that is consistent with the input f(x) in the sense that it produces exactly the same measurements. We show that this solution can be computed by multivariate filtering. We also characterize the stability of the system (condition number). Finally, we illustrate the theory by presenting a new example of interlaced sampling using splines.

Our goal in this paper is to set a theoretical basis for the comparison of re-sampling and interpolation methods. We consider the general problem of the approximation of an arbitrary continuously-defined function f(x)-not necessarily bandlimited-when we vary the sampling step T. We present an accurate L/sup 2/ computation of the induced approximation error as a function of T for a general class of linear approximation operators including interpolation and other kinds of projectors. This new quantitative result provides exact expressions for the asymptotic development of the error as T/spl rarr/0, and also sharp (asymptotically exact) upper bounds.

Describes an image processing system for the automatic assessment of ejection fraction (EF) from noisy SPECT left ventricular myocardial perfusion images. A segmentation scheme consisting of an elliptical coordinate transformation, matched filtering and dynamic contour tracking detects the endo- and epicardial boundaries. Computation of EF is performed based on the epi- rather than the endocardial boundary, which proves to be more robust for images with low signal-to-noise ratios. The computation incorporates anatomical constraints of constant myocardial mass and smooth cardiac variation. The algorithm was tested on different image modalities and shows good linear agreement with EFs obtained from conventional approaches based on planar gated blood pool imaging (PET: y=8.7+1.07x, r=0.84, Technetium-99m MIBI SPECT: y=5.0+0.80x, r=0.90, Thallium-201 SPECT: y=23.5+0.82x, r=0.77, RMS error=10.0). The proposed scheme may be an alternative for EF computation without the need for additional image acquisitions.

We propose to re-size images using an oblique projection operator instead of the orthogonal one in order to obtain faster, simpler, and more general algorithms. The main advantage is that it becomes perfectly feasible to use higher order models (e.g., splines of degree n/spl ges/3). We develop the theoretical background and present a simple and practical implementation procedure that uses B-splines. Experiments show that the proposed algorithm consistently outperforms the standard interpolation method and that it essentially provides the same performance as the optimal procedure (least squares solution) with considerably less computations.

We introduce a fast simple method for computing the real continuous wavelet transform (CWT). The approach achieves O(N) complexity per scale and the filter coefficients can be analytically obtained by a simple integration. Our method is to use P wavelets per octave and to approximate them with their oblique projection onto a space defined by a compact scaling function. The wavelet templates are expanded to larger sizes (octaves) using the two-scale relation and zero padded filtering. Error bounds are presented to justify the use of an oblique projection over an orthogonal one.

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.

The order of a wavelet transform is typically given by the number of vanishing moments of the analysis wavelet. The Strang-Fix (1971) conditions imply that the error for an orthogonal wavelet approximation at scale a=2/sup -i/ globally decays as a/sup L/, where L is the order of the transform. This is why, for a given number of scales, higher order wavelet transforms usually result in better signal approximations. We show that this result carries over for the general biorthogonal case and that the rate of decay of the error is determined by the order properties of the synthesis scaling function alone. We also derive asymptotic error formulas and show that biorthogonal wavelet transforms are equivalent to their corresponding orthogonal projector as the scale goes to zero. These results strengthen Sweldens (see Numer. Math., vol.68, no.3, p.377-401, 1994) earlier analysis and confirm that the approximation power of biorthogonal and (semi-)orthogonal wavelet expansions is essentially the same. Finally, we compare the asymptotic performance of various wavelet transforms and briefly discuss the advantages of splines. We also indicate how the smoothness of the basis functions is beneficial in reducing the approximation error.

We present an overview of the various uses of the wavelet transform (WT) in medicine and biology. We start by describing the wavelet properties that are the most important for biomedical applications. In particular we provide an interpretation of the the continuous wavelet transform (CWT) as a prewhitening multiscale matched filter. We also briefly indicate the analogy between the WT and some of the the biological processing that occurs in the early components of the auditory and visual system. We then review the uses of the WT for the analysis of 1-D physiological signals obtained by phonocardiography, electrocardiography (ECG), mid electroencephalography (EEG), including evoked response potentials. Next, we provide a survey of wavelet developments in medical imaging. These include biomedical image processing algorithms (e.g., noise reduction, image enhancement, and detection of microcalcifications in mammograms), image reconstruction and acquisition schemes (tomography, and magnetic resonance imaging (MRI)), and multiresolution methods for the registration and statistical analysis of functional images of the brain (positron emission tomography (PET) and functional MRI (fMRI)). In each case, we provide the reader with same general background information and a brief explanation of how the methods work.

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.

We introduce a general framework for the efficient computation of the continuous wavelet transform (CWT). The method allows arbitrary sampling along the scale axis, and achieves O(N) complexity per scale where N is the length of the signal. Our approach makes use of a compactly supported scaling function to approximate the analyzing wavelet. We derive error bounds on the wavelet approximation and show how to obtain any desired level of accuracy through the use of higher order representations. Finally, we present examples of implementation for different wavelets using polynomial spline approximations.

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.

We consider a general weighted least squares approximation problem with a membrane spline regularization term. The key parameters in this formulation are the weighting factors which provide the possibility of a spatial adaptation. We prove that the corresponding space-varying variational problem is well posed, and propose a novel multigrid computational solution. This multiresolution relaxation scheme uses three image pyramids (input data, weights, and current solution) and allows for a very efficient computation with an effective O(N) complexity, where N is the number of pixels. This general multigrid solver can be useful for a variety of image processing tasks. In particular, we propose new multigrid solutions for noise reduction in images (adaptive smoothing spline), interpolation/reconstruction of missing image data, and image segmentation using an adaptive extension of the K-means clustering algorithm.

This paper describes a new approach to the characterization of texture properties at multiple scales using the wavelet transform. The analysis uses an overcomplete wavelet decomposition, which yields a description that is translation invariant. It is shown that this representation constitutes a tight frame of l(2) and that it has a fast iterative algorithm. A texture is characterized by a set of channel variances estimated at the output of the corresponding filter bank. Classification experiments with l(2) Brodatz textures indicate that the discrete wavelet frame (DWF) approach is superior to a standard (critically sampled) wavelet transform feature extraction. These results also suggest that this approach should perform better than most traditional single resolution techniques (co-occurrences, local linear transform, and the like). A detailed comparison of the classification performance of various orthogonal and biorthogonal wavelet transforms is also provided. Finally, the DWF feature extraction technique is incorporated into a simple multicomponent texture segmentation algorithm, and some illustrative examples are presented.

The purpose of this paper is to derive optimal spline algorithms for the enlargement or reduction of digital images by arbitrary (noninteger) scaling factors. In our formulation, the original and rescaled signals are each represented by an interpolating polynomial spline of degree n with step size one and Delta, respectively. The change of scale is achieved by determining the spline with step size Delta that provides the closest approximation of the original signal in the L(2)-norm. We show that this approximation can be computed in three steps: (i) a digital prefilter that provides the B-spline coefficients of the input signal, (ii) a resampling using an expansion formula with a modified sampling kernel that depends explicitly on Delta, and (iii) a digital postfilter that maps the result back into the signal domain. We provide explicit formulas for n=0, 1, and 3 and propose solutions for the efficient implementation of these algorithms. We consider image processing examples and show that the present method compares favorably with standard interpolation techniques. Finally, we discuss some properties of this approach and its connection with the classical technique of bandlimiting a signal, which provides the asymptotic limit of our algorithm as the order of the spline tends to infinity.

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.

We present a general framework for the design of discrete geometrical transformation operators, including rotations and scaling. The first step is to fit the discrete input image with a continuous model that provides an exact interpolation at the pixel locations. The corresponding image model is selected within a certain subspace V(/spl phi/)/spl sub/L/sub 2/(R/sup P/) that is generated from the integer translates of a generating function /spl phi/; particular examples of this construction include polynomial spline and bandlimited signal representations. Next, the geometrical transformation is applied to the fitted model, and the result is re-projected onto the representation space. This procedure yields a solution that is optimal in the least squares sense. We show that this method can be implemented exactly using a combination of digital filters and a re-sampling step that uses a modified sampling kernel. We then derive explicit implementation formulas for the piecewise constant and cubic spline image models. Finally, we consider image processing examples and show that the present method compares very favorably with a standard interpolation that uses the same model.>

Fast algorithms for the evaluation of running windowed Fourier and continuous wavelet transforms are presented. The analysis functions approximate complex-modulated Gaussians as closely as desired and may be optimally localized in time and frequency. The Gabor filtering is performed indirectly by convolving a premodulated signal with a Gaussian-like window and demodulating the output. The window functions are either B-splines dilated by an integer factor m or quasi-Gaussians of arbitrary size generated from the n-fold convolution of a symmetrical exponential. Both approaches result in a recursive implementation with a complexity independent of the window size (O(N)).>

A new approach to the characterization of texture properties at multiple scales using an overcomplete wavelet transform is described. It is shown that this representation constitutes a tight frame of l/sub 2/, and that it has a fast iterative algorithm. A texture is characterized by a set of channel variances estimated at the output of the corresponding filter-bank. Classification experiments with 12 Brodatz textures indicate that the discrete wavelet frame (DWF) approach is superior to a standard (critically sampled) wavelet transform feature extraction. This result also suggests that this approach should perform better than most traditional single resolution techniques (co-occurrences, local linear transform, etc. . .). A detailed comparison of the classification performance of various orthogonal and biorthogonal wavelet transforms is provided. The DWF feature extraction technique is incorporated into a simple multiple-component texture segmentation algorithm. Some examples are presented.>

The authors are concerned with the derivation of general methods for the L/sub 2/ approximation of signals by polynomial splines. The main result is that the expansion coefficients of the approximation are obtained by linear filtering and sampling. The authors apply those results to construct a L/sub 2/ polynomial spline pyramid that is a parametric multiresolution representation of a signal. This hierarchical data structure is generated by repeated application of a REDUCE function (prefilter and down-sampler). A complementary EXPAND function (up-sampler and post-filter) allows a finer resolution mapping of any coarser level of the pyramid. Four equivalent representations of this pyramid are considered, and the corresponding REDUCE and EXPAND filters are determined explicitly for polynomial splines of any order n (odd). Some image processing examples are presented. It is demonstrated that the performance of the Laplacian pyramid can be improved significantly by using a modified EXPAND function associated with the dual representation of a cubic spline pyramid.>

The authors describe a novel digital filtering algorithms for the processing and representation of signals using polynomial splines. The classical polynomial spline interpolation problem is considered. It is found that it can be solved efficiently by recursive digital filtering. This result also yields a simple procedure for signal differentiation. Filters that efficiently solve the problem of smoothing spline approximations are derived. This technique is a regularized version of spline interpolation and is therefore less sensitive to noise. It is applied to the design of a robust edge detection algorithm with an adjustable scale parameter. A filtering/sampling algorithm for least squares spline approximation is described. This data reduction technique is applied to the generation of a cubic spline image pyramid that is found to compare favorably with the Gauss/Laplace pyramid.>

The least-squares polynomial spline approximation of a signal g(t) in L/sub 2/(R) is obtained by projecting g(t) on S/sup n/(R) (the space of polynomial splines of order n). It is shown that this process can be linked to the classical problem of cardinal spline interpolation by first convolving g(t) with a B-spline of order n. More specifically, the coefficients of the B-spline interpolation of order 2n+1 of the sampled filtered sequence are identical to the coefficients of the least-squares approximation of g(t) of order n. It is shown that this approximation can be obtained from a succession of three basic operations: prefiltering, sampling, and postfiltering, which confirms the parallel with the classical sampling/reconstruction procedure for bandlimited signals. The frequency responses of these filters are determined for three equivalent spline representations using alternative sets of shift-invariant basis functions of S/sup n/(R): the standard expansion in terms of B-spline coefficients, a representation in terms of sampled signal values, and a representation using orthogonal basis functions.>

A family of nonorthogonal polynomial spline wavelet transforms is considered. These transforms are fully reversible and can be implemented efficiently. The corresponding wavelet functions have a compact support. It is proven that these B-spline wavelets converge to Gabor functions (modulated Gaussian) pointwise and in all L/sub p/-norms with 1>

Least squares approximation problems that are regularized with specified highpass stabilizing kernels are discussed. For each problem, there is a family of discrete regularization filters (R-filters) which allow an efficient determination of the solutions. These operators are stable symmetric lowpass filters with an adjustable scale factor. Two decomposition theorems for the z-transform of such systems are presented. One facilitates the determination of their impulse response, while the other allows an efficient implementation through successive causal and anticausal recursive filtering. A case of special interest is the design of R-filters for the first- and second-order difference operators. These results are extended for two-dimensional signals and, for illustration purposes, are applied to the problem of edge detection. This leads to a very efficient implementation (8 multiplies+10 adds per pixel) of the optimal Canny edge detector based on the use of a separable second-order R-filter.>

Efficient algorithms for the continuous representation of a discrete signal in terms of B-splines (direct B-spline transform) and for interpolative signal reconstruction (indirect B-spline transform) with an expansion factor m are described. Expressions for the z-transforms of the sampled B-spline functions are determined and a convolution property of these kernels is established. It is shown that both the direct and indirect spline transforms involve linear operators that are space invariant and are implemented efficiently by linear filtering. Fast computational algorithms based on the recursive implementations of these filters are proposed. A B-spline interpolator can also be characterized in terms of its transfer function and its global impulse response (cardinal spline of order n). The case of the cubic spline is treated in greater detail. The present approach is compared with previous methods that are reexamined from a critical point of view. It is concluded that B-spline interpolation correctly applied does not result in a loss of image resolution and that this type of interpolation can be performed in a very efficient manner.>

An unsupervised texture segmentation system using texture features obtained from a combination of spatial filters and nonlinear operators is described. Local texture features are evaluated in parallel by a succession of four basic operations: (1) a convolution for local structure detection (local linear transform); (2) a first nonlinearity of the form f(x)= mod x mod /sup alpha /; (3) an iterative smoothing operator; and (4) a second nonlinearity g(x). The Karhunen-Loeve transform is used to reduce the dimensionality of the resulting feature vector, and segmentation is achieved by thresholding or clustering in feature space. The combination of nonlinearities f(x)= mod x mod /sup alpha / (in particular, alpha =2) and g(x)=log x maximizes texture discrimination, and results in a description with variances approximately constant for all feature components and texture regions. This latter property improves the performance of both feature reduction and clustering algorithms significantly.>

An approach is described for unsupervised segmentation of textured images. Local texture properties are extracted using local linear transforms that have been optimized for maximal texture discrimination. Local statistics (texture energy measures) are estimated at the output of an equivalent filter bank by means of a nonlinear transformation (absolute value) followed by an iterative Gaussian smoothing algorithm. This procedure generates a multiresolution sequence of feature planes with a half-octave scale progression. A feature reduction technique is then applied to the data and is determined by simultaneously diagonalizing scatter matrices evaluated at two different spatial resolutions. This approach provides a good approximation of R.A. Fisher's (1950) multiple linear discriminants and has the advantage of requiring no a priori knowledge. This feature reduction methods appears to be an improvement on the commonly used Karhunen-Loeve transform and allows efficient texture segmentation based on simple thresholding.< >

The authors present a linear feature-reduction technique for multicomponent or textured image segmentation. The transformation matrix is computed by simultaneously diagonalizing scatter matrices evaluated at two different spatial resolutions. Under reasonable conditions, this transform closely approximates the generalized Fisher linear discriminants which are optimal for region separability. Experimental examples suggest that this technique is superior to the Karhunen-Loeve transform for texture segmentation.< >

The Karhunen-Loeve transform (KLT) is applied to the analysis of dynamic sequences of thermograms describing the temporal evolution of body surface temperature following the application of an external thermal stimulus. The KLT may be evaluated either along the spatial or temporal dimensions of the data; the duality of both representations is emphasized. An example is presented to illustrate that the KLT allows an efficient data reduction and facilities tumor detection by highlighting physiologically important abnormalities in the time behavior of thermal patterns.>

The commenter points out that the recursive structure of the running Fourier transform has been investigated by a number of authors not quoted by Amin in his work (see ibid., vol.75, no.11, p.1537-1538, 1987). He also holds that the main idea behind the generalization presented in the third section stems from the properties of a yet more general class of features that can be computed using the same recursive structure. Additionally, he discusses Amin's treatment of the general form of a weighted running Fourier coefficient and points out drawbacks of expressing the weighting function as a sum of 2 M+1 geometric series and treating each term separately. The author replies to the various comments.>

An extended family of Discrete Fourier transforms is introduced. These transforms, which may be implemented by using FFTs, allow the computation of pseudo-cyclic convolutions by multiplication in the transform domain. The choice of a suitable transform (DFT1/4) or the combined use of two complementary transforms allows a fast and efficient computation of aperiodic convolutions of waveforms of duration N by using N-point transforms that require no zero padding. Finally, all members of this family are shown to be equivalent asymptotically to the Karhunen-Loève transform of an arbitrary wide sense stationary process.

The sum and difference of two random variables with same variances are decorrelated and define the principal axes of their associated joint probability function. Therefore, sum and difference histograms are introduced as an alternative to the usual co-occurrence matrices used for texture analysis. Two maximum likelihood texture classifiers are presented depending on the type of object used for texture characterization (sum and difference histograms or some associated global measures). Experimental results indicate that sum and difference histograms used conjointly are nearly as powerful as cooccurrence matrices for texture discrimination. The advantage of the proposed texture analysis method over the conventional spatial gray level dependence method is the decrease in computation time and memory storage.