Principal Investigator: Prof. Michael Unser
Participants: Arash Amini, Emrah Bostan, Julien Fageot, Ulugbek S. Kamilov, Masih Nilchian, Ha Q. Nguyen, Pedram Pad, Zsuzsanna Püspöki, Martin Storath, Pouya Deghani Tafti, John Paul Ward
Period: April 1, 2011—March 31, 2016.
In recent years, the focus of the research in signal processing has shifted away from the classical linear paradigm and its intimate links with the theory of stationary Gaussian processes. Instead of considering Fourier transforms and performing quadratic optimization, researchers are presently favoring wavelet-like representations and have adopted "sparsity" as design paradigm.
Our objective in this project was to develop a unifying operator-based framework for signal processing that provides the "sparse" counterpart of the classical theory. This framework was missing until now. To that end, we have introduced a new family of sparse stochastic processes that are continuously defined and ruled by differential equations. We have proposed a principled approach to construct the corresponding wavelet-like sparsifying transforms. Within this framework, we have been able to establish a rigorous correspondence between maximum a posteriori (MAP) estimation and the variational reconstruction of signals with sparsity-promoting regularization. We have also investigated the minimum-mean-square solution and proposed algorithmic solutions for the denoising of Lévy processes.
As primary application of our framework, we have proposed a principled approach to discretize ill-posed inverse problems and to compute statistical estimates of signals. We have demonstrated the benefits of our approach with the development of new reconstruction algorithms for emerging bioimaging modalities including x-ray phase-contrast tomography, superresolution fluorescence microscopy, digital-holography microscopy, refractive-index tomography, as well as phase-contrast magnetic-resonance imaging for the measurement of flow fields.
Another contribution is the construction of an extended family of steerable wavelets, first in 2D and then in 3D, which had not been achieved before. We have taken advantage of the parametric form of our wavelets to optimize them for specific biomedical image-processing tasks such as the attenuation of noise, the detection of junctions, the analysis of textures and morphological components, both in 2D and 3D.
This document describes the work performed during the 5 years of support by the European Research Council for the project ERC-AG No. 267439.
Michael Unser and Pouya Tafti
The defining element of our formulation is some admissible differential operator L. Three fundamental ingredients are then associated to this operator: an extended family of sparse stochastic processes that is whitened by L; some corresponding regularization functionals (Gibbs energies) that embody our prior knowledge on the form of the signal; and a wavelet-like sparsifying transform. A prototypical example is the derivative operator; it goes hand-in-hand with the Lévy processes, the total-variation semi-norm, and the Haar wavelet transform.
The main achievement of FUN-SP is the development of a comprehensive theory of sparse stochastic processes that fits our premises. These processes are specified as solutions of stochastic differential equations driven by non-Gaussian (or Lévy) noise [26, 30]. The driving noise constitutes the unpredictable component of the model and is referred to as the "innovation."
We have introduced a theoretical framework for the rigorous characterization of SSP as generalized functions (a.k.a. tempered distributions). Our main tool is the characteristic functional , which is the infinite-dimensional counterpart of the characteristic function used in the classic theory of probabilities. We have shown that the probability distributions (joint or marginal) of such processes are infinite-divisible and that this property implies a heavy-tail behavior that is the statistical indicator of sparsity .
The proposed framework also includes the more traditional families of Gaussian and non-Gaussian CARMA processes which are ruled by ordinary differential equations. We have shown that the information that is contained in the sampled version of these processes is sufficient to correctly identify the underlying continuous-domain differential system [4, 27].
Our body of research is substantial enough to justify the writing of a textbook entirely devoted to the subject (Unser-Tafti, Cambridge University Press, 2014).
As SSPs are continuous-domain entities, they lend themselves naturally to the discretization of signal-recovery problems by projecting the statistical model onto a suitable reconstruction basis. This allowed us to specify several corresponding signal estimators, including the maximum a posteriori (MAP) and the minimum-mean-square-error (MMSE) solutions . In that respect, one of our main achievements has been a message-passing algorithm that returns the optimal solution (MMSE estimator) for the denoising of Lévy processes .
In general, however, non-Gaussian MMSE estimators are intractable numerically. This has motivated us to investigate alternative solutions. We have focused on MAP estimators, which are solutions of a Gibbs energy-minimization problem. A pleasing outcome of our theory is that MAP for SSP is compatible with sparsity-promoting regularization and ℓ1 minimization. Accordingly, we have invested a substantial effort in the resolution of non-quadratic energy-minimization problems. Specific contributions include methods to regularize vector fields  and to partition vector and manifold data [21, 23] subject to the total-variation constraint.
A specific example of SSP is the compound-Poisson process, which is piecewise-constant and compatible with the Potts model. We have shown how to reconstruct such signals efficiently with the help of the Vitterbi algorithm [31, 34].
In a slightly different vein, we have exploited the connection between SSP and exponential splines. Our strongest result is the proof of the optimality of such splines for the interpolation of alpha-stable AR(1) processes . We have also applied exponential splines to resize images using more conventional identification and regression techniques . Finally, we have characterized the approximation properties of nonuniform Sobolev splines , which can be seen as the multidimensional extension of the exponential splines.
The key assumption that underlies sparsity is the existence of a dictionary of functions that allows for a concise description of the signal. It has motivated us to propose general mechanisms to design wavelet-like transforms.
Since a SSP is specified as the solution of a stochastic differential equation (SDE), it is possible, at least formally, to decouple this SSP by the application of the "whitening" operator L. This uncovers the innovation w = L s, which is the driving term of the differential equation.
While ideal "whitening" is not feasible in practice, the alternative is to expand the process in an operator-like wavelet basis that acts as a smoothed version of the operator. To systematize this approach, we have introduced a general multidimensional framework for the construction of operator-like wavelet bases [7, 19]. The enabling component is a generalized B-spline that is associated to the whitening operator L and that can be interpreted as a "localized" version of its Green's function. The remarkable aspect is that our approach is generic and can be applied to any admissible linear SDE.
When the operator L commutes with dilations (as is the case for pure or fractional derivatives), the solution of the SDE is self-similar—it is a fractal process. In such a scenario, the operator-like wavelets are conventional wavelets that are shifted and rescaled replicates of a single mother wavelet. Our formalism lends itself well to the explicit derivation of the wavelet-domain statistics of such processes . An interesting finding is that these become increasingly Gaussian-like as the scale gets coarser, in agreement with the behavior observed for natural images.
For the special case of alpha-stable processes, we were able to derive the expression of the mutual information in any transformed domain. This allowed us to verify that the operator-like wavelet expansion is a good approximation of an independent-component analysis (ICA) . In particular, we have evaluated the ICA for Lévy processes and observed a remarkable transition from the discrete-cosine transform (the Karhunen-Loève transform in the Gaussian regime with alpha = 2) to the Haar wavelet transform, which becomes optimal at sparser regimes with alpha < 1.
A fruitful recipe to improve the invariance properties of wavelet-like representations is to introduce redundancy. To achieve perfect rotation invariance, we have constructed steerable wavelets by applying the multiorder Riesz transform to an isotropic tight wavelet frame. We have characterized the decay of such wavelets  and proposed a unifying 2D parameterization in terms of circular harmonics . We have extended the approach to higher dimensions by working out a connection with spherical harmonics and nodal basis functions on the sphere . We have also optimized the shape of the radial-frequency profiles of the filters to improve the spatial localization of the wavelets .
Additional contributions to the topic of wavelet design include divergence-free wavelets to denoise flow fields  and new complex wavelets whose phase information is directly related to the scale of the object (ISBI 2015).
We have developed specific wavelet-based algorithms to analyze biomedical images, which demonstrates the practical merits of our formulation.
Operator-like wavelets facilitate the transposition of classic image-analysis methods to wavelet-based ones. For instance, one can apply gradient-like wavelets to perform edge detection at multiple scales. An interesting challenge is then to attempt the recovery of the original image from the mere multi-scale contour map. We have developed an algorithmic solution that takes advantage of the tight-frame property of our transform and applies sparsity constraints to compensate for the missing information. We have perfected this representation by optimizing the wavelet profiles for optimal localization in the spatial domain [ICIP 2014]. This work has been rewarded with the Best Student Award at the IEEE International Conference on Image Processing (ICIP) in 2014.
We have applied our wavelet framework to the detection of local centers of symmetry . Our approach relies on a remarkable property of circular harmonic wavelets; namely, their ability to perfectly decorrelate isotropic random fields. (The same holds true for isotropic fractal processes, including fractional-Brownian fields and their sparse counterparts.) We have also applied our parameterization of steerable wavelets to design junction detectors . These detectors were optimized for maximal directionality.
Equipped with our operator-based framework, we could design the first instance of steerable wavelets in 3D . These were applied to problems such as the attenuation of noise in images, the analysis of textures, as well as the identification of directional structures (filaments, axons) using a morphological-component analysis.
We have also extended our framework to higher dimensions with the help of spherical harmonics .
We did release some dedicated software packages to facilitate the diffusion of our mathematical tools and algorithms.
Part of this offering is geared towards the designers of algorithms. It includes Matlab toolboxes that implement the steerable wavelets in WPs 2b, 3a, 3b.
Another instance is a general framework for the benchmarking of algorithms for single-molecule localization microscopy (see WP 4c).
To address the demands of the users of bioimage-analysis software, including biologists, we have also spent a significant effort to implement and document user-friendly Java plugins for ImageJ/Icy, which are public-domain, multiplatform pieces of software. The resulting packages for interactive analysis of images are
In accordance with our general guidelines, these software contributions are all multiplatform and freely available on the web.
A significant part of the project was to apply the theoretical tools and signal transformations of WPs 1-3 to concrete problems in bioimaging.
We have introduced a general approach to discretize ill-posed linear inverse problems. Our approach is reminiscent of the "finite-element" method for solving PDEs . The idea is to select a set of basis functions and to project both the statistical model and the imaging physics onto the reconstruction space they define. The stochastic process is then decoupled by expressing the whitening operator in this basis, too, which enables one to derive the statistical distribution of the discretized version of the signal. By using the property that the underlying statistical distributions are infinitely divisible, we have shown that the MAP estimator of the signal minimizes a sparsity-promoting functional. We have also formulated an efficient iterative-thresholding algorithm, with the twist that the optimal thresholding function is tied to the underlying statistical model.
The framework also lends itself to the discretization in wavelet bases (Unser-Tafti, Chapter 11) and results in a reconstruction algorithm that is a slight variation of FISTA. We have used cycle spinning to improve the quality of wavelet-based reconstructions, while justifying the scheme theoretically .
In addition, we have investigated reconstruction methods based on belief propagation  and approximate message passing . To assess the full potential of non-Gaussian statistical reconstructions, we have derived the optimal MMSE solution for the denoising of Lévy processes . This has allowed us to measure the performance gap with the more standard MAP estimator and to propose a simple modification of the thresholding function to achieve optimal performance . The enabling component there is consistent cycle spinning , which makes the connection with standard wavelet-based denoising.
Additional contributions to the topic include the introduction of a new regularization functional—the structure-tensor total variation —and a variational formulation for the simultaneous segmentation and reconstruction of an image based on a piecewise-constant model [39, 50].
Within our framework, we have developed an iterative method to reconstruct high-resolution phase-contrast tomograms in collaboration with Prof. Stampanoni at the Swiss Synchrotron Light Source . We could demonstrate that our second-generation methods improve image quality, as compared to the direct filtered-back-projection (FBP) family of algorithms that had been used until then. Alternatively, one is able to decrease the number of views—and, hence, the acquisition time as well as the radiation exposure—up to a factor ten without significant degradation in image quality. To further accelerate our method, we have also implemented a preconditioned variant that boils down to an iterative form of the FBP . Additionally, we have investigated the influence of the basis functions on the quality of reconstruction . The largest part of this work is described in Masih Nilchian's Ph.D. thesis, defended in May 2015, which was rewarded by the 2015 Research Prize of the Swiss Society of Biomedical Engineering.
Finally, we have developed a novel reconstruction algorithm for interior tomography that operates along lines and minimizes a generalized version of total variation . We have applied sparsity priors to improve the reconstruction of the phase in x-ray differential-phase radiography, too .
Our primary contribution to the topic of deconvolution microscopy has been the introduction of a second-order extension of total variation, which is much better suited to the smoothness and directionality properties of biological structures. By focusing on second-order derivative operators, we have come up with a regularization functional—the nuclear norm of the Hessian—that is invariant to translation, scaling, and rotation . We find that our new sparsity-promoting regularizer consistently outperforms the more standard total-variation-based solution. In particular, our approach did rank second in an international 3D-deconvolution contest (ISBI 2014); it was only surpassed by a dictionary-based method that was cleverly attuned to the data at hand.
Another related contribution is an integrated framework for the joint deconvolution-and-segmentation of 3D fluorescence micrographs .
We have investigated signal-processing methods to improve the quantitative estimation of the phase in digital-holography microscopy (DHM). In particular, we have developed a novel adaptive scheme that relies on the Riesz transform to demodulate holograms . We have also reformulated the reconstruction of in-line holograms as a nonlinear (e.g., quadratic) inverse problem and developed a corresponding iterative solver . In addition, we have implemented a novel 2D phase-unwrapping algorithm that uses sparsity constraints to improve the quality of the reconstruction .
In collaboration with the team of Prof. Psaltis, we have developed a novel DHM approach for 3D refractive tomography. Our scheme benefits from an improved, nonlinear forward model that is based on the beam-propagation method. The computational pillar of our formulation is a "learning" algorithm that applies a variant of the back-propagation algorithm of neural networks to iteratively reconstruct the image . We have applied our system to the imaging of cells and could obtain a substantial improvement over the state of the art, which is based on diffraction tomography.
The research of this WP (and to some extent WP 4a) resulted in two EPFL Ph.D. theses: "Compressed Optical Imaging," defended by Aurélien Bourquard in December 2013, and "Sparsity-Driven Statistical Inference for Inverse Problems," defended by Ulugbek S. Kamilov in March 2015.
Since the state of the art of the reconstruction of magnetic-resonance images (MRI) is more mature than in optics, we have concentrated our activity on phase-contrast imaging and the recovery of flow fields, which are still emerging topics. In particular, we have introduced a new sparsity-promoting regularizer that penalizes the singular values of the Jacobian of the flow field . We have also investigated the use of spatiotemporal regularization to improve the quality of reconstruction, which did necessitate the development of new algorithms (ISBI 2013).
All the work of FUN-SP done on MRI reconstruction is described in Emrah Bostan's thesis (see also ). There, he did also address an outstanding problem in optical imaging (not initially mentioned in the grant), namely, the recovery of the phase of a wavefront based on intensity information only. Starting from the transport-of-intensity equation, he proposed a novel variational formulation of the phase-recovery problem . His results are the best obtained so far with incoherent imaging; they are qualitatively on par with those of DHM that take advantage of coherent light. For the full details, see the EPFL Ph.D. thesis "Sparsity-Based Data Reconstruction Models for Biomedical Imaging," defended by Emrah Bostan in May 2016.