Showing papers in "Siam Journal on Imaging Sciences in 2008"

A New Alternating Minimization Algorithm for Total Variation Image Reconstruction

Abstract: We propose, analyze, and test an alternating minimization algorithm for recovering images from blurry and noisy observations with total variation (TV) regularization. This algorithm arises from a new half-quadratic model applicable to not only the anisotropic but also the isotropic forms of TV discretizations. The per-iteration computational complexity of the algorithm is three fast Fourier transforms. We establish strong convergence properties for the algorithm including finite convergence for some variables and relatively fast exponential (or $q$-linear in optimization terminology) convergence for the others. Furthermore, we propose a continuation scheme to accelerate the practical convergence of the algorithm. Extensive numerical results show that our algorithm performs favorably in comparison to several state-of-the-art algorithms. In particular, it runs orders of magnitude faster than the lagged diffusivity algorithm for TV-based deblurring. Some extensions of our algorithm are also discussed.

Bregman Iterative Algorithms for $\ell_1$-Minimization with Applications to Compressed Sensing

Abstract: We propose simple and extremely efficient methods for solving the basis pursuit problem $\min\{\|u\|_1 : Au = f, u\in\mathbb{R}^n\},$ which is used in compressed sensing. Our methods are based on Bregman iterative regularization, and they give a very accurate solution after solving only a very small number of instances of the unconstrained problem $\min_{u\in\mathbb{R}^n} \mu\|u\|_1+\frac{1}{2}\|Au-f^k\|_2^2$ for given matrix $A$ and vector $f^k$. We show analytically that this iterative approach yields exact solutions in a finite number of steps and present numerical results that demonstrate that as few as two to six iterations are sufficient in most cases. Our approach is especially useful for many compressed sensing applications where matrix-vector operations involving $A$ and $A^\top$ can be computed by fast transforms. Utilizing a fast fixed-point continuation solver that is based solely on such operations for solving the above unconstrained subproblem, we were able to quickly solve huge instances of compressed sensing problems on a standard PC.

A Nonlinear Inverse Scale Space Method for a Convex Multiplicative Noise Model

Abstract: We are motivated by a recently developed nonlinear inverse scale space method for image denoising [M. Burger, G. Gilboa, S. Osher, and J. Xu, Commun. Math. Sci., 4 (2006), pp. 179-212; M. Burger, S. Osher, J. Xu, and G. Gilboa, in Variational, Geometric, and Level Set Methods in Computer Vision, Lecture Notes in Comput. Sci. 3752, Springer, Berlin, 2005, pp. 25-36], whereby noise can be removed with minimal degradation. The additive noise model has been studied extensively, using the Rudin-Osher-Fatemi model [L. I. Rudin, S. Osher, and E. Fatemi, Phys. D, 60 (1992), pp. 259-268], an iterative regularization method [S. Osher, M. Burger, D. Goldfarb, J. Xu, and W. Yin, Multiscale Model. Simul., 4 (2005), pp. 460-489], and the inverse scale space flow [M. Burger, G. Gilboa, S. Osher, and J. Xu, Commun. Math. Sci., 4 (2006), pp. 179-212; M. Burger, S. Osher, J. Xu, and G. Gilboa, in Variational, Geometric, and Level Set Methods in Computer Vision, Lecture Notes in Comput. Sci. 3752, Springer, Berlin, 2005, pp. 25-36]. However, the multiplicative noise model has not yet been studied thoroughly. Earlier total variation models for the multiplicative noise cannot easily be extended to the inverse scale space, due to the lack of global convexity. In this paper, we review existing multiplicative models and present a new total variation framework for the multiplicative noise model, which is globally strictly convex. We extend this convex model to the nonlinear inverse scale space flow and its corresponding relaxed inverse scale space flow. We demonstrate the convergence of the flow for the multiplicative noise model, as well as its regularization effect and its relation to the Bregman distance. We investigate the properties of the flow and study the dependence on flow parameters. The numerical results show an excellent denoising effect and significant improvement over earlier multiplicative models.

Sparse and Redundant Modeling of Image Content Using an Image-Signature-Dictionary

Abstract: Modeling signals by sparse and redundant representations has been drawing considerable attention in recent years. Coupled with the ability to train the dictionary using signal examples, these techniques have been shown to lead to state-of-the-art results in a series of recent applications. In this paper we propose a novel structure of such a model for representing image content. The new dictionary is itself a small image, such that every patch in it (in varying location and size) is a possible atom in the representation. We refer to this as the image-signature-dictionary (ISD) and show how it can be trained from image examples. This structure extends the well-known image and video epitomes, as introduced by Jojic, Frey, and Kannan [in Proceedings of the IEEE International Conference on Computer Vision, 2003, pp. 34-41] and Cheung, Frey, and Jojic [in Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2005, pp. 42-49], by replacing a probabilistic averaging of patches with their sparse representations. The ISD enjoys several important features, such as shift and scale flexibilities, and smaller memory and computational requirements, compared to the classical dictionary approach. As a demonstration of these benefits, we present high-quality image denoising results based on this new model.

Efficient Reconstruction of Piecewise Constant Images Using Nonsmooth Nonconvex Minimization

Abstract: We consider the restoration of piecewise constant images where the number of the regions and their values are not fixed in advance, with a good difference of piecewise constant values between neighboring regions, from noisy data obtained at the output of a linear operator (e.g., a blurring kernel or a Radon transform). Thus we also address the generic problem of unsupervised segmentation in the context of linear inverse problems. The segmentation and the restoration tasks are solved jointly by minimizing an objective function (an energy) composed of a quadratic data-fidelity term and a nonsmooth nonconvex regularization term. The pertinence of such an energy is ensured by the analytical properties of its minimizers. However, its practical interest used to be limited by the difficulty of the computational stage which requires a nonsmooth nonconvex minimization. Indeed, the existing methods are unsatisfactory since they (implicitly or explicitly) involve a smooth approximation of the regularization term and often get stuck in shallow local minima. The goal of this paper is to design a method that efficiently handles the nonsmooth nonconvex minimization. More precisely, we propose a continuation method where one tracks the minimizers along a sequence of approximate nonsmooth energies $\{J_\eps\}$, the first of which being strictly convex and the last one the original energy to minimize. Knowing the importance of the nonsmoothness of the regularization term for the segmentation task, each $J_\eps$ is nonsmooth and is expressed as the sum of an $\ell_1$ regularization term and a smooth nonconvex function. Furthermore, the local minimization of each $J_{\eps}$ is reformulated as the minimization of a smooth function subject to a set of linear constraints. The latter problem is solved by the modified primal-dual interior point method, which guarantees the descent direction at each step. Experimental results are presented and show the effectiveness and the efficiency of the proposed method. Comparison with simulated annealing methods further shows the advantage of our method.

Variational Assimilation of Fluid Motion from Image Sequence

Abstract: In this paper, a variational technique derived from optimal control theory is used in order to realize a dynamically consistent motion estimation of a whole fluid image sequence The estimation is conducted through an iterative process involving a forward integration of a given dynamical model followed by a backward integration of an adjoint evolution law By combining physical conservation laws and image observations, a physically grounded temporal consistency is imposed, and the quality of the motion estimation is significantly improved The method is validated on two synthetic image sequences provided by numerical simulation of fluid flows and on real world meteorological examples

Direct Elastic Imaging of a Small Inclusion

Abstract: In this paper we consider the problem of locating a small three-dimensional elastic inclusion, using arrays of elastic source transmitters and receivers. This procedure yields the multistatic response matrix that is characteristic of the elastic inclusion. We show how the eigenvalue structure of this matrix can be employed within the framework of a noniterative method of MUSIC (multiple signal classification) type in order to retrieve the elastic inclusion. We illustrate our reconstruction procedure with a variety of computational examples from synthetic, noiseless, and severely noisy data.

Total Variation Regularization for Image Denoising, III. Examples.

Abstract: Let $\mathcal{F}(\mathbb{R}^{2})$ be the family of bounded nonnegative Lebesgue measurable functions on $\mathbb{R}^{2}$. Suppose $s\in\mathcal{F}(\mathbb{R}^{2})$ and $\gamma:\mathbb{R}\rightarrow[0,\infty)$ is zero at zero, positive away from zero, and convex. For $f\in\mathcal{F}(\Omega)$ let $F(f)=\int_\Omega\gamma(f(x)-s(x))\,d\mathcal{L}^{2}x$; here $\mathcal{L}^{2}$ is Lebesgue measure on $\mathbb{R}^2$. In the denoising literature $F$ would be called a fidelity in that it measures how much $f$ differs from $s$ which could be a noisy grayscale image. Suppose $0<\epsilon<\infty$ and let ${\bf m}^{loc}_{\epsilon}(F)$ be the set of those $f\in\mathcal{F}(\mathbb{R}^{2})$ such that ${\bf TV}(f)<\infty$ and $\epsilon{\bf TV}(f)+F(f)\leq\epsilon{\bf TV}(g)+F(g)$ for $g\in{\bf k}(f)$; here ${\bf TV}(f)$ is the total variation of $f$ and ${\bf k}(f)$ is the set of $g\in\mathcal{F}(\mathbb{R}^{2})$ such that $g=f$ off some compact subset of $\mathbb{R}^{2}$. A member of ${\bf m}^{loc}_{\epsilon}(F)$ is called a total variation regularization of $s$ (with smoothing parameter $\epsilon$). Rudin, Osher, and Fatemi in [Phys. D, 60 (1992), pp. 259-268] and Chan and Esedog¯lu in [SIAM J. Appl. Math., 65 (2005), pp. 1817-1837] have studied total variation regularizations of $s$ where $\gamma(y)=y^2/2$ and $\gamma(y)=y$, $y\in\mathbb{R}$, respectively. It turns out that the character of a member of ${\bf m}^{loc}_{\epsilon}(F)$ changes quite a bit as $\gamma$ changes. In [SIAM J. Imaging Sci., 1 (2008), pp. 400-417] the family ${\bf m}^{loc}_{\epsilon}(F)$ was described in complete detail when $s$ is the indicator function of a compact convex subset of $\mathbb{R}^{2}$. Our main purpose in this paper is to describe, in complete detail, ${\bf m}^{loc}_{\epsilon}(F)$ when $s$ is the indicator function of either $S=([0,1]\times[0,-1])\cup([-1,0]\times[0,1])$ or $S=\{x\in\mathbb{R}^{2}:|x-{\bf c}_+|\leq 1\}\cup\{x\in\mathbb{R}^{2}:|x-{\bf c}_-|\leq 1\}$, where, for some $l\in[1,\infty)$, ${\bf c}_\pm=(\pm l,0)$. We believe these examples reveal a great deal about the nature of total variation regularizations. For example, it is said that total variation denoising preserves edges; while this is certainly true in many cases and in comparison with other denoising methods, the examples given in sections 2squares and 2circles provide evidence to the contrary. In addition, one can test computational schemes for total variation regularization against these examples. We will also establish what we believe to be a number of interesting properties of ${\bf m}^{loc}_{\epsilon}(F)$.

Combining Image Reconstruction and Image Analysis with an Application to Two-Dimensional Tomography

Abstract: The tasks of image reconstruction from measured data and the analysis of the resulting images are more or less strictly separated. One group of scientists computes by applying reconstruction algorithms to the images; the other then operates on these images to enhance the analysis. First attempts at combining image reconstruction and image analysis, in a nonsystematic way, are known as Lambda tomography or Tikhonov-Phillips methods with $\ell_1$-norms or with level-set methods. The aim of this paper is to provide a general tool to combine these two steps; i.e., even in the reconstruction step the future image analysis step is taken into account, leading to a new reconstruction kernel. Here we concentrate on linear methods. As a practical example we consider the image reconstruction problem in computerized tomography followed by an edge detection. We calculate a new reconstruction kernel and present results from simulations.

Edge Illumination and Imaging of Extended Reflectors

Abstract: We use the singular value decomposition of the array response matrix, frequency by frequency, to image selectively the edges of extended reflectors in a homogeneous medium. We show with numerical simulations in an ultrasound regime, and analytically in the Fraunhofer diffraction regime, that information about the edges is contained in the singular vectors for singular values that are intermediate between the large ones and zero. These transition singular vectors beamform selectively from the array onto the edges of the reflector cross-section facing the array, so that these edges are enhanced in imaging with travel-time migration. Moreover, the illumination with the transition singular vectors is done from the sources at the edges of the array. The theoretical analysis in the Fraunhofer regime shows that the singular values transition to zero at the index $N^\star(\om) = |{\cal A}||{\cal B}|/(\lambda L)^2 $. Here $|{\cal A}|$ and $|{\cal B}|$ are the areas of the array and the reflector cross-section, respectively, $\omega$ is the frequency, $\lambda$ is the wavelength, and $L$ is the range. Since $(\lambda L)^2/|{\cal A}|$ is the area of the focal spot size at range $L$, we see that $N^\star(\om)$ is the number of focal spots contained in the reflector cross-section. The ultrasound simulations are in an extended Fraunhofer regime. The simulation results are, however, qualitatively similar to those obtained theoretically in the Fraunhofer regime. The numerical simulations indicate, in addition, that the subspaces spanned by the transition singular vectors are robust with respect to additive noise when the array has a large number of elements.

A Parameterization Method for the Computation of Invariant Tori and Their Whiskers in Quasi-Periodic Maps: Explorations and Mechanisms for the Breakdown of Hyperbolicity

Abstract: In two previous papers [J. Differential Equations, 228 (2006), pp. 530–579; Discrete Contin. Dyn. Syst. Ser. B, 6 (2006), pp. 1261–1300] we have developed fast algorithms for the computations of invariant tori in quasi‐periodic systems and developed theorems that assess their accuracy. In this paper, we study the results of implementing these algorithms and study their performance in actual implementations. More importantly, we note that, due to the speed of the algorithms and the theoretical developments about their reliability, we can compute with confidence invariant objects close to the breakdown of their hyperbolicity properties. This allows us to identify a mechanism of loss of hyperbolicity and measure some of its quantitative regularities. We find that some systems lose hyperbolicity because the stable and unstable bundles approach each other but the Lyapunov multipliers remain away from 1. We find empirically that, close to the breakdown, the distances between the invariant bundles and the Lyapun...

Cascadic Multiresolution Methods for Image Deblurring

Abstract: This paper investigates the use of cascadic multiresolution methods for image deblurring. Iterations with a conjugate gradient-type method are carried out on each level, and terminated by a stopping rule based on the discrepancy principle. Prolongation is carried out by nonlinear edge-preserving operators, which are defined via PDEs associated with Perona-Malik or total variation-type models. Computed examples demonstrate the effectiveness of the methods proposed.

An Efficient Method for Studying Weak Resonant Double Hopf Bifurcation in Nonlinear Systems with Delayed Feedbacks

Abstract: An efficient method, called the perturbation‐incremental scheme (PIS), is proposed to study, both qualitatively and quantitatively, the delay‐induced weak or high‐order resonant double Hopf bifurcation and the dynamics arising from the bifurcation of nonlinear systems with delayed feedback. The scheme is described in two steps, namely, the perturbation and the incremental steps, when the time delay and the feedback gain are taken as the bifurcation parameters. As for applications, the method is employed to investigate the delay‐induced weak resonant double Hopf bifurcation and dynamics in the van der Pol–Duffing and the Stuart–Landau systems with delayed feedback. For bifurcation parameters close to a double Hopf point, all solutions arising from the resonant bifurcation are classified qualitatively and expressed approximately in a closed form by the perturbation step of the PIS. Although the analytical expression may not be accurate enough for bifurcation parameters away from the double Hopf point, it is...

Inpainting by Flexible Haar-Wavelet Shrinkage

Abstract: We present novel wavelet-based inpainting algorithms. Applying ideas from anisotropic regularization and diffusion, our models can better handle degraded pixels at edges. We interpret our algorithms within the framework of forward-backward splitting methods in convex analysis and prove that the conditions for ensuring their convergence are fulfilled. Numerical examples illustrate the good performance of our algorithms.

Convex Source Support and Its Application to Electric Impedance Tomography

Abstract: The aim in electric impedance tomography is to recover the conductivity inside a physical body from boundary measurements of current and voltage. In many situations of practical importance, the investigated object has known background conductivity but is contaminated by inhomogeneities. In this work, we try to extract all possible information about the support of such inclusions inside a two-dimensional object from only one pair of measurements of impedance tomography. Our noniterative and computationally cheap method is based on the concept of the convex source support, which stems from earlier works of Kusiak, Sylvester, and the authors. The functionality of our algorithm is demonstrated by various numerical experiments.

Area Distances of Convex Plane Curves and Improper Affine Spheres

Abstract: The area distance to a convex plane curve is an important concept in computer vision. In this paper we describe a strong link between area distances and improper affine spheres. Based on this link, we propose an extremely fast algorithm to compute the inner area distance. Moreover, the concepts of the theory of affine spheres lead to a new definition of an area distance on the outer part of a convex plane curve. On the other hand, area distances provide a good geometrical understanding of improper affine spheres.

Topology Preserving Linear Filtering Applied to Medical Imaging

Abstract: One of the central problems of medical imaging is the three-dimensional (3D) visualization of body parts. The 3D volume can be viewed in slices, but the extraction of a part requires a segmentation process. Inasmuch as body parts are distinguishable by their various densities, a widely accepted method for extracting an organ is to extract isodensity surfaces by a simple threshold. Unfortunately, the density of organs, arteries, etc. varies spatially due to morphology, and no unique threshold allows one to extract the organs boundaries. The snake or active contour methods have attempted to capture these boundaries as smooth and overall contrasted surfaces. The snake method suffers, however, from severe drawbacks. The contour has to be initialized near the boundary. In addition, many body parts have too complex a topology. In this paper we focus on another idea, which is to preprocess the image before thresholding. The preprocessing aims at the homogeneity of the different parts while preserving small features. Starting from a recent seminal work by Grady and Funka-Lea [in Computer Vision and Mathematical Methods in Medical and Biomedical Image Analysis: ECCV 2004 Workshops CVAMIA and MMBIA, Prague, Czech Republic, May 2004, Revised Selected Papers, Springer, Berlin, 2004, pp. 230-245], several linear heat equations on images will be compared. They stem from nonlinear partial differential equations or from their associated nonlinear filters. By linearizing these processes one obtains more accurate topology preserving methods. These linear filters will be tested comparatively to visualize challenging angiography images of arteries. A salient fact of the method will emerge. By a concentration phenomenon, peaks in the image histogram become much more concentrated under the linear heat equations, thus permitting us to fix the thresholds defining the surfaces without supervision. Automatic extraction can be performed in this way for angiography images taken at a one-year or longer delay.

Geodesic Evolution Laws—A Level-Set Approach

Abstract: Motion of curves governed by geometric evolution laws, such as mean curvature flow and surface diffusion, is the basis for many algorithms in image processing. If the images to be processed are defined on nonplanar surfaces, the geometric evolution laws have to be restricted to the surface and turn into geodesic evolution laws. In this paper we describe efficient algorithms for geodesic mean curvature flow and geodesic surface diffusion within a level-set approach. Thereby we compare approaches with an explicit representation of the surface by a triangulated surface mesh and an implicit surface representation as the zero-level surface of a level-set function. As an application we present the numerical treatment of the classical model of Rudin, Osher, and Fatemi to denoise images on surfaces.

Topological Entropy of Braids on the Torus

Abstract: A fast method is presented for computing the topological entropy of braids on the torus. This work is motivated by the need to analyze large braids when studying two-dimensional flows via the braiding of a large number of particle trajectories. Our approach is a generalization of Moussafir's technique for braids on the sphere. Previous methods for computing topological entropies include the Bestvina--Handel train-track algorithm and matrix representations of the braid group. However, the Bestvina--Handel algorithm quickly becomes computationally intractable for large braid words, and matrix methods give only lower bounds, which are often poor for large braids. Our method is computationally fast and appears to give exponential convergence towards the exact entropy. As an illustration we apply our approach to the braiding of both periodic and aperiodic trajectories in the sine flow. The efficiency of the method allows us to explore how much extra information about flow entropy is encoded in the braid as the number of trajectories becomes large.

A Model for p53 Dynamics Triggered by DNA Damage

Abstract: Several recent experiments on DNA‐damage‐induced signaling networks in mammalian cells have shown interesting dynamics in p53 protein expression during the repair cycle. Pulses of p53 are produced, whose frequency and amplitude are fairly independent of the amount of damage, but the probability of a cell exhibiting this pulsatile behavior increases with damage. This phenomenon has been described as a “digital oscillator.” We present here a simple model oscillator comprising two species, p53 and Mdm2, which is activated by the Atm kinase. The Atm kinase exhibits bistable switch‐like behavior. The network dynamics essentially consists of the core p53 oscillator, which is turned ON/OFF by the Atm switch, which is in turn activated by DNA damage. The complex dynamics are thus explained by the modular nature of the network and are fairly independent of the biological details. A stochastic model of the network dynamics reveals that the pulsatile behavior is robust to intrinsic noise of the protein components an...

Patterns and Features of Families of Traveling Waves in Large-Scale Neuronal Networks

Abstract: We study traveling wave solutions of a system of integro‐differential equations which describe the activity of large‐scale networks of excitatory neurons on spatially extended domains. The independent variables are the activity level u of a population of excitatory neurons, which have long range connections, and a recovery variable v. We have found a critical value of the parameter β ($\beta_* > 0$) that appears in the equation for v, at which the eigenvalues of the linearization of the system around the rest state $(u,v)=(0,0)$ change from real to complex. In contrast to previous studies which analyzed properties of traveling waves when the eigenvalues are real, we examine the range $\beta > \beta_*$, where the eigenvalues are complex. In this case we show that there is a range of parameters over which families of wave fronts and 1‐pulse and more general N‐pulse waves can coexist as stable solutions. In two space dimensions our numerical experiments show how single‐ring and multiring waves form in respon...

Segmentation under Occlusions Using Selective Shape Prior

Abstract: In this work, we address the problem of segmenting multiple objects, under possible occlusions, in a level set framework. A variational energy that incorporates a piecewise constant representation of the image in terms of the object regions and the object spatial order is proposed. To resolve occluded boundaries, prior knowledge of the shape of objects is also introduced within the segmentation energy. By minimizing the above energy, we solve the segmentation with depth problem, i.e., estimating the object boundaries, the object intensities, and the spatial order. The segmentation with depth problem was originally dealt with by the Nitzberg-Mumford-Shiota (NMS) variational formulation, which proposes segmentation energies for each spatial order. We discuss the relationships and show the computational advantages of our formulation over the NMS model, mainly due to our treatment of spatial order estimation within a single energy. A novelty here is that the spatial order information available in the image model is used to dynamically impose prior shape constraints only to occluded boundaries. Also presented are experiments on synthetic and real images that have promising results.

A Study of 1PI Algorithms for a General Class of Curves

Abstract: We extend an efficient cone beam transform inversion formula, developed earlier by the authors for smooth curves with positive torsion, to a more general class of helix-like curves. These curves are allowed to have negative torsion, and they can be nonsmooth at isolated points. The notions of turns and PI segments are extended. The new class is defined by several geometric conditions which impose a tradeoff between the length of critical PI lines (which reflect how severely the positivity of torsion is violated) and the extent to which the curve bends between neighboring turns. The main property of curves from this class is that critical PI lines are allowed to be arbitrarily close to the set $U$ where reconstruction is possible, but are not allowed to intersect it. Some of the conditions that define the class turn out to be common for many known trajectories, so we investigate separately the properties of the Crofton symbol of PI segments of curves that satisfy these conditions. The results of the investigation are then used to develop an efficient filtered backprojection algorithm. Numerical experiments conducted with a clock phantom demonstrate good image quality.

Some Canonical Bifurcations in the Swift–Hohenberg Equation

Abstract: We study the nature and stability of stationary solutions $u(x)$ of the fourth order Swift–Hohenberg equation on a bounded domain $(0,L)$ with boundary conditions $u=0$ and $u''=0$ at $x=0$ and $x=L$. It is well known that as L increases, the set of stationary solutions becomes increasingly complex. Numerical studies have exhibited two interesting types of structures in the bifurcation diagram for $(L,u)$. In this paper we demonstrate through a center manifold analysis how these structures arise naturally near certain bifurcation points, and that there are no others. We also analyze their stability properties.

Coherent Structures Generated by Inhomogeneities in Oscillatory Media

Abstract: We investigate the effect of spatially localized inhomogeneities on a spatially homogeneous oscillation in a reaction‐diffusion system. In dimension up to two, we find sources and contact defects, that is, the inhomogeneity may either send out phase waves or act as a weak sink. We show that small inhomogeneities cannot act as sources in more than two space dimensions. We also derive asymptotics for wavenumbers and group velocities in the far field. The results are established rigorously for radially symmetric inhomogeneities in reaction‐diffusion systems, and for arbitrary inhomogeneities in a modulation equation approximation.