In cryo-electron microscopy (cryo-EM), a microscope generates a top view of a sample of randomly oriented copies of a molecule. The problem of single particle reconstruction (SPR) from cryo-EM is to use the resulting set of noisy two-dimensional projection images taken at unknown directions to reconstruct the three-dimensional (3D) structure of the molecule. In some situations, the molecule under examination exhibits structural variability, which poses a fundamental challenge in SPR. The heterogeneity problem is the task of mapping the space of conformational states of a molecule. It has been previously suggested that the leading eigenvectors of the covariance matrix of the 3D molecules can be used to solve the heterogeneity problem. Estimating the covariance matrix is challenging, since only projections of the molecules are observed, but not the molecules themselves. In this paper, we formulate a general problem of covariance estimation from noisy projections of samples. This problem has intimate connections with matrix completion problems and high-dimensional principal component analysis. We propose an estimator and prove its consistency. When there are finitely many heterogeneity classes, the spectrum of the estimated covariance matrix reveals the number of classes. The estimator can be found as the solution to a certain linear system. In the cryo-EM case, the linear operator to be inverted, which we term the projection covariance transform, is an important object in covariance estimation for tomographic problems involving structural variation. Inverting it involves applying a filter akin to the ramp filter in tomography. We design a basis in which this linear operator is sparse and thus can be tractably inverted despite its large size. We demonstrate via numerical experiments on synthetic datasets the robustness of our algorithm to high levels of noise.

Segmentation is a fundamental task for extracting semantically meaningful regions from an image. The goal of segmentation algorithms is to accurately assign object labels to each image location. However, image-noise, shortcomings of algorithms, and image ambiguities cause uncertainty in label assignment. Estimating the uncertainty in label assignment is important in multiple application domains, such as segmenting tumors from medical images for radiation treatment planning. One way to estimate these uncertainties is through the computation of posteriors of Bayesian models, which is computationally prohibitive for many practical applications. On the other hand, most computationally efficient methods fail to estimate label uncertainty. We therefore propose in this paper the Active Mean Fields (AMF) approach, a technique based on Bayesian modeling that uses a mean-field approximation to efficiently compute a segmentation and its corresponding uncertainty. Based on a variational formulation, the resulting convex model combines any label-likelihood measure with a prior on the length of the segmentation boundary. A specific implementation of that model is the Chan-Vese segmentation model (CV), in which the binary segmentation task is defined by a Gaussian likelihood and a prior regularizing the length of the segmentation boundary. Furthermore, the Euler-Lagrange equations derived from the AMF model are equivalent to those of the popular Rudin-Osher-Fatemi (ROF) model for image denoising. Solutions to the AMF model can thus be implemented by directly utilizing highly-efficient ROF solvers on log-likelihood ratio fields. We qualitatively assess the approach on synthetic data as well as on real natural and medical images. For a quantitative evaluation, we apply our approach to the icgbench dataset.

As modern scientific image datasets typically consist of a large number of images of high resolution, devising methods for their accurate and efficient processing is a central research task. In this paper, we consider the problem of obtaining the steerable principal components of a dataset, a procedure termed "steerable PCA" (steerable principal component analysis). The output of the procedure is the set of orthonormal basis functions which best approximate the images in the dataset and all of their planar rotations. To derive such basis functions, we first expand the images in an appropriate basis, for which the steerable PCA reduces to the eigen-decomposition of a block-diagonal matrix. If we assume that the images are well localized in space and frequency, then such an appropriate basis is the prolate spheroidal wave functions (PSWFs). We derive a fast method for computing the PSWFs expansion coefficients from the images' equally spaced samples, via a specialized quadrature integration scheme, and show that the number of required quadrature nodes is similar to the number of pixels in each image. We then establish that our PSWF-based steerable PCA is both faster and more accurate then existing methods, and more importantly, provides us with rigorous error bounds on the entire procedure.

Accurate estimation of the initial pressure distribution in photoacoustic computed tomography (PACT) depends on knowledge of the sound speed distribution. However, the sound speed distribution is typically unknown. Further, the initial pressure and sound speed distributions cannot both, in general, be stably recovered from PACT measurements alone. In this work, a joint reconstruction (JR) method for the initial pressure distribution and a low-dimensional parameterized model of the sound speed distribution is proposed. By employing information about the structure of the sound speed distribution, both the initial pressure and sound speed can be accurately recovered. The JR problem is solved by use of a proximal optimization method that allows constraints and non-smooth regularization functions for the initial pressure distribution. The gradients of the cost function with respect to the initial pressure and sound speed distributions are calculated by use of an adjoint state method that has the same per-iteration computational cost as calculating the gradient with respect to the initial pressure distribution alone. This approach is evaluated through 2D computer-simulation studies for a small animal imaging model and by application to experimental in vivo measurements of a mouse.

Photoacoustic computed tomography (PACT) is an emerging imaging modality that exploits optical contrast and ultrasonic detection principles to form images of the photoacoustically induced initial pressure distribution within tissue. The PACT reconstruction problem corresponds to an inverse source problem in which the initial pressure distribution is recovered from measurements of the radiated wavefield. A major challenge in transcranial PACT brain imaging is compensation for aberrations in the measured data due to the presence of the skull. Ultrasonic waves undergo absorption, scattering and longitudinal-to-shear wave mode conversion as they propagate through the skull. To properly account for these effects, a wave-equation-based inversion method should be employed that can model the heterogeneous elastic properties of the skull. In this work, a forward model based on a finite-difference time-domain discretization of the three-dimensional elastic wave equation is established and a procedure for computing the corresponding adjoint of the forward operator is presented. Massively parallel implementations of these operators employing multiple graphics processing units (GPUs) are also developed. The developed numerical framework is validated and investigated in computer19 simulation and experimental phantom studies whose designs are motivated by transcranial PACT applications.

Single-particle cryo-electron microscopy (cryo-EM) has recently joined X-ray crystallography and NMR spectroscopy as a high-resolution structural method to resolve biological macromolecules. In a cryo-EM experiment, the microscope produces images called micrographs. Projections of the molecule of interest are embedded in the micrographs at unknown locations, and under unknown viewing directions. Standard imaging techniques first locate these projections (detection) and then reconstruct the 3-D structure from them. Unfortunately, high noise levels hinder detection. When reliable detection is rendered impossible, the standard techniques fail. This is a problem, especially for small molecules. In this paper, we pursue a radically different approach: we contend that the structure could, in principle, be reconstructed directly from the micrographs, without intermediate detection. The aim is to bring small molecules within reach for cryo-EM. To this end, we design an autocorrelation analysis technique that allows one to go directly from the micrographs to the sought structures. This involves only one pass over the micrographs, allowing online, streaming processing for large experiments. We show numerical results and discuss challenges that lay ahead to turn this proof-of-concept into a complementary approach to state-of-the-art algorithms.

We propose regularization schemes for deformable registration and efficient algorithms for their numerical approximation. We treat image registration as a variational optimal control problem. The deformation map is parametrized by its velocity. Tikhonov regularization ensures well-posedness. Our scheme augments standard smoothness regularization operators based on - and -seminorms with a constraint on the divergence of the velocity field, which resembles variational formulations for Stokes incompressible flows. In our formulation, we invert for a stationary velocity field and a mass source map. This allows us to explicitly control the compressibility of the deformation map and by that the determinant of the deformation gradient. We also introduce a new regularization scheme that allows us to control shear. We use a globalized, preconditioned, matrix-free, reduced space (Gauss-)Newton-Krylov scheme for numerical optimization. We exploit variable elimination techniques to reduce the number of unknowns of our system; we only iterate on the reduced space of the velocity field. Our current implementation is limited to the two-dimensional case. The numerical experiments demonstrate that we can control the determinant of the deformation gradient without compromising registration quality. This additional control allows us to avoid oversmoothing of the deformation map. We also demonstrate that we can promote or penalize shear whilst controlling the determinant of the deformation gradient.

Diffeomorphic registration using optimal control on the diffeomorphism group and on shape spaces has become widely used since the development of the large deformation diffeomorphic metric mapping (LDDMM) algorithm. More recently, a series of algorithms involving sub-Riemannian constraints have been introduced in which the velocity fields that control the shapes in the LDDMM framework are constrained in accordance with a specific deformation model. Here, we extend this setting by considering, for the first time, inequality constraints in order to estimate surface deformations that only allow for atrophy, introducing for this purpose an algorithm that uses the augmented Lagrangian method. We prove the existence of solutions of the associated optimal control problem and the consistency of our approximation scheme. These developments are illustrated by numerical experiments on simulated and real data.

Following the seminal work of Nesterov, accelerated optimization methods have been used to powerfully boost the performance of first-order, gradient based parameter estimation in scenarios where second-order optimization strategies are either inapplicable or impractical. Not only does accelerated gradient descent converge considerably faster than traditional gradient descent, but it also performs a more robust local search of the parameter space by initially overshooting and then oscillating back as it settles into a final configuration, thereby selecting only local minimizers with a basis of attraction large enough to contain the initial overshoot. This behavior has made accelerated and stochastic gradient search methods particularly popular within the machine learning community. In their recent PNAS 2016 paper, , Wibisono, Wilson, and Jordan demonstrate how a broad class of accelerated schemes can be cast in a variational framework formulated around the Bregman divergence, leading to continuum limit ODEs. We show how their formulation may be further extended to infinite dimensional manifolds (starting here with the geometric space of curves and surfaces) by substituting the Bregman divergence with inner products on the tangent space and explicitly introducing a distributed mass model which evolves in conjunction with the object of interest during the optimization process. The coevolving mass model, which is introduced purely for the sake of endowing the optimization with helpful dynamics, also links the resulting class of accelerated PDE based optimization schemes to fluid dynamical formulations of optimal mass transport.

In the shape analysis approach to computer vision problems, one treats shapes as points in an infinite-dimensional Riemannian manifold, thereby facilitating algorithms for statistical calculations such as geodesic distance between shapes and averaging of a collection of shapes. The performance of these algorithms depends heavily on the choice of the Riemannian metric. In the setting of plane curve shapes, attention has largely been focused on a two-parameter family of first order Sobolev metrics, referred to as elastic metrics. They are particularly useful due to the existence of simplifying coordinate transformations for particular parameter values, such as the well-known square-root velocity transform. In this paper, we extend the transformations appearing in the existing literature to a family of isometries, which take any elastic metric to the flat metric. We also extend the transforms to treat piecewise linear curves and demonstrate the existence of optimal matchings over the diffeomorphism group in this setting. We conclude the paper with multiple examples of shape geodesics for open and closed curves. We also show the benefits of our approach in a simple classification experiment.

Several imaging algorithms including patch-based image denoising, image time series recovery, and convolutional neural networks can be thought of as methods that exploit the manifold structure of signals. While the empirical performance of these algorithms is impressive, the understanding of recovery of the signals and functions that live on manifold is less understood. In this paper, we focus on the recovery of signals that live on a union of surfaces. In particular, we consider signals living on a union of smooth band-limited surfaces in high dimensions. We show that an exponential mapping transforms the data to a union of low-dimensional subspaces. Using this relation, we introduce a sampling theoretical framework for the recovery of smooth surfaces from few samples and the learning of functions living on smooth surfaces. The low-rank property of the features is used to determine the number of measurements needed to recover the surface. Moreover, the low-rank property of the features also provides an efficient approach, which resembles a neural network, for the local representation of multidimensional functions on the surface. The direct representation of such a function in high dimensions often suffers from the curse of dimensionality; the large number of parameters would translate to the need for extensive training data. The low-rank property of the features can significantly reduce the number of parameters, which makes the computational structure attractive for learning and inference from limited labeled training data.

Photoacoustic tomography (PAT) is a non-invasive imaging modality that requires recovering the initial data of the wave equation from certain measurements of the solution outside the object. In the standard PAT measurement setup, the used data consist of time-dependent signals measured on an observation surface. In contrast, the measured data from the recently invented full-field detection technique provide the solution of the wave equation on a spatial domain at a single instant in time. While reconstruction using classical PAT data has been extensively studied, not much is known for the full field PAT problem. In this paper, we build mathematical foundations of the latter problem for variable sound speed and settle its uniqueness and stability. Moreover, we introduce an exact inversion method using time-reversal and study its convergence. Our results demonstrate the suitability of both the full field approach and the proposed time-reversal technique for high resolution photoacoustic imaging.

We target the problem of estimating the center of mass of objects in noisy two-dimensional images. We assume that the noise dominates the image, and thus many standard approaches are vulnerable to estimation errors, e.g., the direct computation of the center of mass and the geometric median which is a robust alternative to the center of mass. In this paper, we define a novel surrogate function to the center of mass. We present a mathematical and numerical analysis of our method and show that it outperforms existing methods for estimating the center of mass of an object in various realistic scenarios. As a case study, we apply our centering method to data from single-particle cryo-electron microscopy (cryo-EM), where the goal is to reconstruct the three-dimensional structure of macromolecules. We show how to apply our approach for a better translational alignment of molecule images picked from experimental data. In this way, we facilitate the succeeding steps of reconstruction and streamline the entire cryo-EM pipeline, saving computational time and supporting resolution enhancement.

Steerable filter pairs that are near quadrature have many image processing applications. This paper proposes a new methodology for designing such filters. The key idea is to design steerable filters by minimizing a departure-from-quadrature function. These minimizing filter pairs are almost exactly in quadrature. The polar part of the filters is nonnegative, monotonic, and highly focused around an axis, and asymptotically the filters achieve exact quadrature. These results are established by exploiting a relation between the filters and generalized Hilbert matrices. These near-quadrature filters closely approximate three dimensional Gabor filters. We experimentally verify the asymptotic mathematical results and further demonstrate the use of these filter pairs by efficient calculation of local Fourier shell correlation of cryogenic electron microscopy.

Tensor-valued and matrix-valued measurements of different physical properties are increasingly available in material sciences and medical imaging applications. The eigenvalues and eigenvectors of such multivariate data provide novel and unique information, but at the cost of requiring a more complex statistical analysis. In this work we derive the distributions of eigenvalues and eigenvectors in the special but important case of symmetric random matrices, , observed with isotropic matrix-variate Gaussian noise. The properties of these distributions depend strongly on the symmetries of the mean tensor/matrix, . When has repeated eigenvalues, the eigenvalues of are not asymptotically Gaussian, and repulsion is observed between the eigenvalues corresponding to the same eigenspaces. We apply these results to diffusion tensor imaging (DTI), with = 3, addressing an important problem of detecting the symmetries of the diffusion tensor, and seeking an experimental design that could potentially yield an isotropic Gaussian distribution. In the 3-dimensional case, when the mean tensor is spherically symmetric and the noise is Gaussian and isotropic, the asymptotic distribution of the first three eigenvalue central moment statistics is simple and can be used to test for isotropy. In order to apply such tests, we use quadrature rules of order ≥ 4 with constant weights on the unit sphere to design a DTI-experiment with the property that isotropy of the underlying true tensor implies isotropy of the Fisher information. We also explain the potential implications of the methods using simulated DTI data with a Rician noise model.

In multiatlas segmentation, one typically registers several atlases to the novel image, and their respective segmented label images are transformed and fused to form the final segmentation. In this work, we provide a new dynamical system perspective for multiatlas segmentation, inspired by the following fact: The transformation that aligns the current atlas to the novel image can be not only computed by direct registration but also inferred from the transformation that aligns the previous atlas to the image together with the transformation between the two atlases. This process is similar to the global positioning system on a vehicle, which gets position by inquiring from the satellite and by employing the previous location and velocity-neither answer in isolation being perfect. To solve this problem, a dynamical system scheme is crucial to combine the two pieces of information; for example, a Kalman filtering scheme is used. Accordingly, in this work, a Kalman multiatlas segmentation is proposed to stabilize the global/affine registration step. The contributions of this work are twofold. First, it provides a new dynamical systematic perspective for standard independent multiatlas registrations, and it is solved by Kalman filtering. Second, with very little extra computation, it can be combined with most existing multiatlas segmentation schemes for better registration/segmentation accuracy.

The extraction of clusters from a dataset which includes multiple clusters and a significant background component is a non-trivial task of practical importance. In image analysis this manifests for example in anomaly detection and target detection. The traditional spectral clustering algorithm, which relies on the leading eigenvectors to detect clusters, fails in such cases. In this paper we propose the which sums the squared values of the first normalized eigenvectors, where can be significantly larger than . We prove that this quantity can be used to separate clusters from the background in unbalanced settings, including extreme cases such as outlier detection. The performance of the algorithm is not sensitive to the choice of , and we demonstrate its application on synthetic and real-world remote sensing and neuroimaging datasets.

We introduce a method to recover a continuous domain representation of a piecewise constant two-dimensional image from few low-pass Fourier samples. Assuming the edge set of the image is localized to the zero set of a trigonometric polynomial, we show the Fourier coefficients of the partial derivatives of the image satisfy a linear annihilation relation. We present necessary and sufficient conditions for unique recovery of the image from finite low-pass Fourier samples using the annihilation relation. We also propose a practical two-stage recovery algorithm which is robust to model-mismatch and noise. In the first stage we estimate a continuous domain representation of the edge set of the image. In the second stage we perform an extrapolation in Fourier domain by a least squares two-dimensional linear prediction, which recovers the exact Fourier coefficients of the underlying image. We demonstrate our algorithm on the super-resolution recovery of MRI phantoms and real MRI data from low-pass Fourier samples, which shows benefits over standard approaches for single-image super-resolution MRI.

In cryo-electron microscopy, the three-dimensional (3D) electric potentials of an ensemble of molecules are projected along arbitrary viewing directions to yield noisy two-dimensional images. The volume maps representing these potentials typically exhibit a great deal of structural variability, which is described by their 3D covariance matrix. Typically, this covariance matrix is approximately low rank and can be used to cluster the volumes or estimate the intrinsic geometry of the conformation space. We formulate the estimation of this covariance matrix as a linear inverse problem, yielding a consistent least-squares estimator. For images of size -by- pixels, we propose an algorithm for calculating this covariance estimator with computational complexity , where the condition number is empirically in the range 10-200. Its efficiency relies on the observation that the normal equations are equivalent to a deconvolution problem in six dimensions. This is then solved by the conjugate gradient method with an appropriate circulant preconditioner. The result is the first computationally efficient algorithm for consistent estimation of the 3D covariance from noisy projections. It also compares favorably in runtime with respect to previously proposed nonconsistent estimators. Motivated by the recent success of eigenvalue shrinkage procedures for high-dimensional covariance matrix estimation, we incorporate a shrinkage procedure that improves accuracy at lower signal-to-noise ratios. We evaluate our methods on simulated datasets and achieve classification results comparable to state-of-the-art methods in shorter running time. We also present results on clustering volumes in an experimental dataset, illustrating the power of the proposed algorithm for practical determination of structural variability.

Multi-energy computed tomography (CT) is an emerging medical image modality with a number of potential applications in diagnosis and therapy. However, high system cost and technical barriers obstruct its step into routine clinical practice. In this study, we propose a framework to realize multi-energy cone beam CT (ME-CBCT) on the CBCT system that is widely available and has been routinely used for radiotherapy image guidance. In our method, a kVp switching technique is realized, which acquires x-ray projections with kVp levels cycling through a number of values. For this kVp-switching based ME-CBCT acquisition, x-ray projections of each energy channel are only a subset of all the acquired projections. This leads to an undersampling issue, posing challenges to the reconstruction problem. We propose a spatial spectral non-local means (ssNLM) method to reconstruct ME-CBCT, which employs image correlations along both spatial and spectral directions to suppress noisy and streak artifacts. To address the intensity scale difference at different energy channels, a histogram matching method is incorporated. Our method is different from conventionally used NLM methods in that spectral dimension is included, which helps to effectively remove streak artifacts appearing at different directions in images with different energy channels. Convergence analysis of our algorithm is provided. A comprehensive set of simulation and real experimental studies demonstrate feasibility of our ME-CBCT scheme and the capability of achieving superior image quality compared to conventional filtered backprojection-type (FBP) and NLM reconstruction methods.

Intensity-based image registration is critical for neuroimaging tasks, such as 3D reconstruction, times-series alignment, and common coordinate mapping. The gradient-based optimization methods commonly used to solve this problem require a careful selection of step-length. This limitation imposes substantial time and computational costs. Here we propose a gradient-independent rigid-motion registration algorithm based on the majorization-minimization (MM) principle. Each iteration of our intensity-based MM algorithm reduces to a simple point-set rigid registration problem with a closed form solution that avoids the step-length issue altogether. The details of the algorithm are presented, and an error bound for its more practical truncated form is derived. The performance of the MM algorithm is shown to be more effective than gradient descent on simulated images and Nissl stained coronal slices of mouse brain. We also compare and contrast the similarities and differences between the MM algorithm and another gradient-free registration algorithm called the block-matching method. Finally, extensions of this algorithm to more complex problems are discussed.