Search

Geometric deep learning is a relatively nascent field that has attracted significant attention in the past few years. This is partly due to the availability of data acquired from non-euclidean domains or features extracted from euclidean-space data that reside on smooth manifolds. For instance, pose data commonly encountered in computer vision reside in Lie groups, while covariance matrices that are ubiquitous in many fields and diffusion tensors encountered in medical imaging domain reside on the manifold of symmetric positive definite matrices. Much of this data is naturally represented as a grid of manifold-valued data. In this paper we present a novel theoretical framework for developing deep neural networks to cope with these grids of manifold-valued data inputs. We also present a novel architecture to realize this theory and call it the ManifoldNet. Analogous to vector spaces where convolutions are equivalent to computing weighted sums, manifold-valued data 'convolutions' can be defined using the weighted Fréchet Mean ([Formula: see text]). (This requires endowing the manifold with a Riemannian structure if it did not already come with one.) The hidden layers of ManifoldNet compute [Formula: see text]s of their inputs, where the weights are to be learnt. This means the data remain manifold-valued as they propagate through the hidden layers. To reduce computational complexity, we present a provably convergent recursive algorithm for computing the [Formula: see text]. Further, we prove that on non-constant sectional curvature manifolds, each [Formula: see text] layer is a contraction mapping and provide constructive evidence for its non-collapsibility when stacked in layers. This captures the two fundamental properties of deep network layers. Analogous to the equivariance of convolution in euclidean space to translations, we prove that the [Formula: see text] is equivariant to the action of the group of isometries admitted by the Riemannian manifold on which the data reside. To showcase the performance of ManifoldNet, we present several experiments using both computer vision and medical imaging data sets.

Diffusion-weighted magnetic resonance imaging (dMRI) is a non-invasive technique to probe the complex micro-architecture of the tissue being imaged. The diffusional properties of the tissue at the imaged resolution are well captured by the ensemble average propagator (EAP), which is a probability density function characterizing the probability of water molecule diffusion. Many properties in the form of imaging 'stains' can then be computed from the EAP that can serve as bio-markers for a variety of diseases. This motivates the development of methods for the accurate estimation of the EAPs from dMRI, which is an actively researched area in dMRI analysis. To this end, in the recent past, dictionary learning (DL) techniques have been applied by many researchers for accurate reconstruction of the EAP fields from dMRI scans of the central nervous system (CNS). However, most of the DL-based methods did not exploit the geometry of the space of the EAPs, which are probability density functions. By exploiting the geometry of the space of probability density functions, it is possible to reconstruct EAPs that satisfy the mathematical properties of a density function and hence yield better accuracy in the EAP field reconstruction. Using a square root density parameterization, the EAPs can be mapped to a unit Hilbert sphere, which is a smooth manifold with well known geometry that we will exploit in our formulation of the DL problem. Thus, in this paper, we present a general formulation of the DL problem for data residing on smooth manifolds and in particular the manifold of EAPs, along with a numerical solution using an alternating minimization method. We then showcase the properties and the performance of our algorithm on the reconstruction of the EAP field in a patch-wise manner from the dMRI data. Through several synthetic, phantom and real data examples, we demonstrate that our non-linear DL-based approach produces accurate and spatially smooth estimates of the EAP field from dMRI in comparison to the state-of-the-art EAP reconstruction method called the MAPL method, as well as the linear DL-based EAP reconstruction approaches. To further demonstrate the accuracy and utility of our approach, we compute an entropic anisotropy measure (HA), that is a function of the well known Rényi entropy, from the EAP fields of control and injured rat spinal cords respectively. We demonstrate its utility as an imaging 'stain' via a quantitative comparison of HA maps computed from EAP fields estimated using our method and competing methods. The quantitative comparison is achieved using a two sample t-test and the results of significance are displayed for a visualization of regions of the spinal cord affected most by the injury.

In this paper, we present a novel generalization of the Volterra Series, which can be viewed as a higher-order convolution, to manifold-valued functions. A special case of the manifold-valued Volterra Series (MVVS) gives us a natural extension of the ordinary convolution to manifold-valued functions that we call, the manifold-valued convolution (MVC). We prove that these generalizations preserve the equivariance properties of the Euclidean Volterra Series and the traditional convolution operator. We present novel deep network architectures using the MVVS and the MVC operations which are then validated via two experiments. These include, (i) movement disorder classification from diffusion magnetic resonance images (dMRI), and (ii) fiber orientation distribution function (fODF) reconstruction from compressed sensed dMRIs. In both the experiments, MVVS and MVC networks outperform the state-of-the-art.

Finding the Riemannian barycenter (center of mass) or the Frechet mean (FM) of manifold-valued data sets is a commonly encountered problem in a variety of fields of science and engineering, including, but not limited to, medical image computing, machine learning, and computer vision. For example, it is encountered in tasks such as atlas construction, clustering, principal geodesic analysis, and so on. Traditionally, algorithms for computing the FM of the manifold-valued data require that the entire data pool be available a priori and not incrementally. Thus, when encountered with new data, the FM needs to be recomputed over the entire pool, which can be computationally and storage inefficient. A computational and storage efficient alternative is to consider a recursive algorithm for computing the FM, which simply updates the previously computed FM when presented with a new data sample. In this chapter we present such an alternative called the inductive/incremental Frechet mean estimator (iFME) for data residing on two well-known Riemannian manifolds, namely the hypersphere S ( d ) and the special orthogonal group SO ( d ) . We prove the asymptotic convergence of iFME to the true FM of the underlying distribution from which the data samples were drawn. Further we present several experiments demonstrating the performance iFME on synthetic and real data sets.

Accurate reconstruction of the ensemble average propagators (EAPs) from undersampled diffusion MRI (dMRI) measurements is a well-motivated, actively researched problem in the field of dMRI acquisition and analysis. A number of approaches based on compressed sensing (CS) principles have been developed for this problem, achieving a considerable acceleration in the acquisition by leveraging sparse representations of the signal. Most recent methods in literature apply undersampling techniques in the (k, q)-space for the recovery of EAP in the joint (x, r)-space. Yet, the majority of these methods follow a pipeline of first reconstructing the diffusion images in the (x, q)-space and subsequently estimating the EAPs through a 3D Fourier transform. In this work, we present a novel approach to achieve the direct reconstruction of P(x, r) from partial (k, q)-space measurements, with geometric constraints involving the parallelism of level-sets of diffusion images from proximal q-space points. By directly reconstructing P(x, r)) from (k, q)-space data, we exploit the incoherence between the 6D sensing and reconstruction domains to the fullest, which is consistent with the CS-theory. Further, our approach aims to utilize the inherent structural similarity (parallelism) of the level-sets in the diffusion images corresponding to proximally-located q-space points in a CS framework to achieve further reduction in sample complexity that could facilitate faster acquisition in dMRI. We compare the proposed method to a state-of-the-art CS based EAP reconstruction method (from joint (k, q)-space) on simulated, phantom and real dMRI data demonstrating the benefits of exploiting the structural similarity in the q-space.

Principal Component Analysis (PCA) is a widely popular dimensionality reduction technique for vector-valued inputs. In the past decade, a nonlinear generalization of PCA, called the Principal Geodesic Analysis (PGA) was developed to tackle data that lie on a smooth manifold. PGA suffers from the same problem as PCA in that, in both the methods, each Principal Component (PC) is a linear combination of the original variables. This makes it very difficult to interpret the PCs especially in high dimensions. This lead to the introduction of sparse PCA (SPCA) in the vector-space input case. In this paper, we present a novel generalization of SPCA, called sparse exact PGA (SEPGA) that can cope with manifold-valued input data and respect the intrinsic geometry of the underlying manifold. Sparsity has the advantage of not only easy interpretability but also computational efficiency. We achieve this by formulating the PGA problem as a minimization of the projection error in conjunction with sparsity constraints enforced on the principal vectors post isomorphic mapping to Rm, where m is the dimension of the manifold on which the data reside. Further, for constant curvature smooth manifolds, we use analytic formulae for the projection error leading to an efficient solution to the SEPGA problem. We present extensive experimental results demonstrating the performance of SEPGA in achieving very good sparse principal components without sacrificing the accuracy of reconstruction. This makes SEPGA accurate and efficient in representing manifold-valued data.

Analyzing data representing multifarious trajectories is central to the many fields in Science and Engineering; for example, trajectories representing a tennis serve, a gymnast’s parallel bar routine, progression/remission of disease and so on. We present a novel geometric algorithm for performing statistical analysis of trajectories with distinct number of samples representing longitudinal (or temporal) data. A key feature of our proposal is that unlike existing schemes, our model is deployable in regimes where each participant provides a different number of acquisitions (trajectories have different number of sample points or temporal span). To achieve this, we develop a novel method involving the parallel transport of the tangent vectors along each given trajectory to the starting point of the respective trajectories and then use the span of the matrix whose columns consist of these vectors, to construct a linear subspace in R(m). We then map these linear subspaces (possibly of distinct dimensions) of R(m) on to a single high dimensional hypersphere. This enables computing group statistics over trajectories by instead performing statistics on the hypersphere (equipped with a simpler geometry). Given a point on the hypersphere representing a trajectory, we also provide a “reverse mapping” algorithm to uniquely (under certain assumptions) reconstruct the subspace that corresponds to this point. Finally, by using existing algorithms for recursive Fréchet mean and exact principal geodesic analysis on the hypersphere, we present several experiments on synthetic and real (vision and medical) data sets showing how group testing on such diversely sampled longitudinal data is possible by analyzing the reconstructed data in the subspace spanned by the first few principal components.

Longitudinal data analysis is an often encountered problem in Medical Image Analysis. A differential geometric treatment of such problems has been reported in literature in the recent past. However, most of these methods require that the trajectory characterizing the evolution of features over time lie on a geodesic emanating from the initial time point on the manifold containing the trajectory. This is a stringent and not necessarily a meaningful requirement. Further, most of these methods impose the restriction that the number of samples on a trajectory be the same across the members of a group of trajectories. At times, this restriction is hard to meet from a practical view point. In this paper, we present a novel formulation of the trajectory analysis problem that overcomes the aforementioned limitations. We represent the trajectories by embedding them in a product Riemannian manifold and endowing it with a Riemannian metric, thereby facilitating the statistical analysis. Finally, we present real data (from MR brain scans of dementia patients) examples depicting the performance of our algorithms.

Principal Component Analysis (PCA) is a fundamental method for estimating a linear subspace approximation to high-dimensional data. Many algorithms exist in literature to achieve a statistically robust version of PCA called RPCA. In this paper, we present a geometric framework for computing the principal linear subspaces in both situations that amounts to computing the intrinsic average on the space of all subspaces (the Grassmann manifold). Points on this manifold are defined as the subspaces spanned by K-tuples of observations. We show that the intrinsic Grassmann average of these subspaces coincide with the principal components of the observations when they are drawn from a Gaussian distribution. Similar results are also shown to hold for the RPCA. Further, we propose an efficient online algorithm to do subspace averaging which is of linear complexity in terms of number of samples and has a linear convergence rate. When the data has outliers, our proposed online robust subspace averaging algorithm shows significant performance (accuracy and computation time) gain over a recently published RPCA methods with publicly accessible code. We have demonstrated competitive performance of our proposed online subspace algorithm method on one synthetic and two real data sets. Experimental results depicting stability of our proposed method are also presented. Furthermore, on two real outlier corrupted datasets, we present comparison experiments showing lower reconstruction error using our online RPCA algorithm. In terms of reconstruction error and time required, both our algorithms outperform the competition.

In this paper, we present novel algorithms to compute robust statistics from manifold-valued data. Specifically, we present algorithms for estimating the robust Frechet Mean (FM) and performing a robust exact-principal geodesic analysis (ePGA) for data lying on known Riemannian manifolds. We formulate the minimization problems involved in both these problems using the minimum distance estimator called the L\(_2\)E. This leads to a nonlinear optimization which is solved efficiently using a Riemannian accelerated gradient descent technique. We present competitive performance results of our algorithms applied to synthetic data with outliers, the corpus callosum shapes extracted from OASIS MRI database, and diffusion MRI scans from movement disorder patients respectively.

Manifold-valued datasets are widely encountered in many computer vision tasks. A non-linear analog of the PCA, called the Principal Geodesic Analysis (PGA) suited for data lying on Riemannian manifolds was reported in literature a decade ago. Since the objective function in PGA is highly non-linear and hard to solve efficiently in general, researchers have proposed a linear approximation. Though this linear approximation is easy to compute, it lacks accuracy especially when the data exhibits a large variance. Recently, an alternative called exact PGA was proposed which tries to solve the optimization without any linearization. For general Riemannian manifolds, though it gives better accuracy than the original (linearized) PGA, for data that exhibit large variance, the optimization is not computationally efficient. In this paper, we propose an efficient exact PGA for constant curvature Riemannian manifolds (CCM-EPGA). CCM-EPGA differs significantly from existing PGA algorithms in two aspects, (i) the distance between a given manifold-valued data point and the principal submanifold is computed analytically and thus no optimization is required as in existing methods. (ii) Unlike the existing PGA algorithms, the descent into codimension-1 submanifolds does not require any optimization but is accomplished through the use of the Rimeannian inverse Exponential map and the parallel transport operations. We present theoretical and experimental results for constant curvature Riemannian manifolds depicting favorable performance of CCM-EPGA compared to existing PGA algorithms. We also present data reconstruction from principal components and directions which has not been presented in literature in this setting., Accepted in CVPR (IEEE Conference on Computer Vision and Pattern Recognition) 2016

Brain atlas construction has attracted significant attention lately in the neuroimaging community due to its application to the characterization of neuroanatomical shape abnormalities associated with various neurodegenerative diseases or neuropsychiatric disorders. Existing shape atlas construction techniques usually focus on the analysis of a single anatomical structure in which the important inter-structural information is lost. This paper proposes a novel technique for constructing a neuroanatomical shape complex atlas based on an information geometry framework. A shape complex is a collection of neighboring shapes - for example, the thalamus, amygdala and the hippocampus circuit - which may exhibit changes in shape across multiple structures during the progression of a disease. In this paper, we represent the boundaries of the entire shape complex using the zero level set of a distance transform function S(x). We then re-derive the relationship between the stationary state wave function ψ(x) of the Schrödinger equation [formula in text] and the eikonal equation [formula in text] satisfied by any distance function. This leads to a one-to-one map (up to scale) between ψ(x) and S(x) via an explicit relationship. We further exploit this relationship by mapping ψ(x) to a unit hypersphere whose Riemannian structure is fully known, thus effectively turn ψ(x) into the square-root of a probability density function. This allows us to make comparisons - using elegant, closed-form analytic expressions - between shape complexes represented as square-root densities. A shape complex atlas is constructed by computing the Karcher mean ψ¯(x) in the space of square-root densities and then inversely mapping it back to the space of distance transforms in order to realize the atlas shape. We demonstrate the shape complex atlas computation technique via a set of experiments on a population of brain MRI scans including controls and epilepsy patients with either right anterior medial temporal or left anterior medial temporal lobectomies.

Divergence measures provide a means to measure the pairwise dissimilarity between "objects," e.g., vectors and probability density functions (pdfs). Kullback-Leibler (KL) divergence and the square loss (SL) function are two examples of commonly used dissimilarity measures which along with others belong to the family of Bregman divergences (BD). In this paper, we present a novel divergence dubbed the Total Bregman divergence (TBD), which is intrinsically robust to outliers, a very desirable property in many applications. Further, we derive the TBD center, called the t-center (using the l(1)-norm), for a population of positive definite matrices in closed form and show that it is invariant to transformation from the special linear group. This t-center, which is also robust to outliers, is then used in tensor interpolation as well as in an active contour based piecewise constant segmentation of a diffusion tensor magnetic resonance image (DT-MRI). Additionally, we derive the piecewise smooth active contour model for segmentation of DT-MRI using the TBD and present several comparative results on real data.

Groupwise registration of a set of shapes represented by unlabeled point sets is a challenging problem since, usually, this involves solving for point correspondence in a nonrigid motion setting. In this paper, we propose a novel and robust algorithm that is capable of simultaneously computing the mean shape, represented by a probability density function, from multiple unlabeled point sets (represented by finite-mixture models), and registering them nonrigidly to this emerging mean shape. This algorithm avoids the correspondence problem by minimizing the Jensen-Shannon (JS) divergence between the point sets represented as finite mixtures of Gaussian densities. We motivate the use of the JS divergence by pointing out its close relationship to hypothesis testing. Essentially, minimizing the JS divergence is asymptotically equivalent to maximizing the likelihood ratio formed from a probability density of the pooled point sets and the product of the probability densities of the individual point sets. We derive the analytic gradient of the cost function, namely, the JS-divergence, in order to efficiently achieve the optimal solution. The cost function is fully symmetric, with no bias toward any of the given shapes to be registered and whose mean is being sought. A by-product of the registration process is a probabilistic atlas, which is defined as the convex combination of the probability densities of the input point sets being aligned. Our algorithm can be especially useful for creating atlases of various shapes present in images and for simultaneously (rigidly or nonrigidly) registering 3D range data sets (in vision and graphics applications), without having to establish any correspondence. We present experimental results on nonrigidly registering 2D and 3D real and synthetic data (point sets).

Diffusion magnetic resonance imaging (MRI) is a relatively new imaging modality which is capable of measuring the diffusion of water molecules in biological systems noninvasively. The measurements from diffusion MRI provide unique clues for extracting orientation information of brain white matter fibers and can be potentially used to infer the brain connectivity in vivo using tractography techniques. Diffusion tensor imaging (DTI), currently the most widely used technique, fails to extract multiple fiber orientations in regions with complex microstructure. In order to overcome this limitation of DTI, a variety of reconstruction algorithms have been introduced in the recent past. One of the key ingredients in several model-based approaches is deconvolution operation which is presented in a unified deconvolution framework in this paper. Additionally, some important computational issues in solving the deconvolution problem that are not addressed adequately in previous studies are described in detail here. Further, we investigate several deconvolution schemes towards achieving stable, sparse, and accurate solutions. Experimental results on both simulations and real data are presented. The comparisons empirically suggest that nonnegative least squares method is the technique of choice for the multifiber reconstruction problem in the presence of intravoxel orientational heterogeneity.

We present an MPM-MAP application for the efficient estimation of piecewise parametric models for motion segmentation. This algorithm permits the simultaneous estimation of: the number of models, the parameters for each model and the regions where each model is applicable. It is based on Bayesian estimation theory, and is theoretically justified by the use of a specific cost function whose expected value decreases at every iteration and by a new model for the posterior marginal distributions which is amenable to the use of fast computational methods. We compare the performance of this method with the most similar segmentation algorithm, the well known Expectation-Maximization algorithm. We present a comparison of the performance of both algorithms using synthetic and real image sequences.

Diffusion tensor imaging can provide the fundamental information required for viewing structural connectivity. However, robust and accurate acquisition and processing algorithms are needed to accurately map the nerve connectivity. In this paper, we present a novel algorithm for extracting and visualizing the fiber tracts in the CNS, specifically in the brain. The automatic fiber tract mapping problem will be solved in two phases, namely a data smoothing phase and a fiber tract mapping phase. In the former, smoothing of the diffusion-weighted data (prior to tensor calculation) is achieved via a weighted TV-norm minimization, which strives to smooth while retaining all relevant detail. For the fiber tract mapping, a smooth 3D vector field indicating the dominant anisotropic direction at each spatial location is computed from the smoothed data. Neuronal fibers are then traced by calculating the integral curves of this vector field. Results are expressed using three modes of visualization: (1) Line integral convolution produces an oriented texture which shows fiber pathways in a planar slice of the data. (2) A streamtube map is generated to present a 3D view of fiber tracts. Additional information, such as degree of anisotropy, can be encoded in the tube radius, or by using color. (3) A particle system form of visualization is also presented. This mode of display allows for interactive exploration of fiber connectivity with no additional preprocessing.

Compressed Sensing (CS) for the acceleration of MR scans has been widely investigated in the past decade. Lately, considerable progress has been made in achieving similar speed ups in acquiring multi-shell high angular resolution diffusion imaging (MS-HARDI) scans. Existing approaches in this context were primarily concerned with sparse reconstruction of the diffusion MR signal \(S({\mathbf {q}})\) in the \({\mathbf {q}}\)-space. More recently, methods have been developed to apply the compressed sensing framework to the 6-dimensional joint \(({\mathbf {k}},{\mathbf {q}})\)-space, thereby exploiting the redundancy in this 6D space. To guarantee accurate reconstruction from partial MS-HARDI data, the key ingredients of compressed sensing that need to be brought together are: (1) the function to be reconstructed needs to have a sparse representation, and (2) the data for reconstruction ought to be acquired in the dual domain (i.e., incoherent sensing) and (3) the reconstruction process involves a (convex) optimization.

A novel adaptation of the unscented Kalman filter (UKF) was recently introduced in literature for simultaneous multi-tensor estimation and fiber tractography from diffusion MRI. This technique has the advantage over other tractography methods in terms of computational efficiency, due to the fact that the UKF simultaneously estimates the diffusion tensors and propagates the most consistent direction to track along. This UKF and its variants reported later in literature however are not intrinsic to the space of diffusion tensors. Lack of this key property can possibly lead to inaccuracies in the multi-tensor estimation as well as in the tractography. In this paper, we propose a novel intrinsic unscented Kalman filter (IUKF) in the space of diffusion tensors which are symmetric positive definite matrices, that can be used for simultaneous recursive estimation of multi-tensors and propagation of directional information for use in fiber tractography from diffusion weighted MR data. In addition to being more accurate, IUKF retains all the advantages of UKF mentioned above. We demonstrate the accuracy and effectiveness of the proposed method via experiments publicly available phantom data from the fiber cup-challenge (MICCAI 2009) and diffusion weighted MR scans acquired from human brains and rat spinal cords.

In this paper, we cast the problem of point cloud matching as a shape matching problem by transforming each of the given point clouds into a shape representation called the Schrodinger distance transform (SDT) representation. This is achieved by solving a static Schrodinger equation instead of the corresponding static Hamilton-Jacobi equation in this setting. The SDT representation is an analytic expression and following the theoretical physics literature, can be normalized to have unit 2 norm -- making it a square-root density, which is identified with a point on a unit Hilbert sphere, whose intrinsic geometry is fully known. The Fisher-Rao metric, a natural metric for the space of densities leads to analytic expressions for the geodesic distance between points on this sphere. In this paper, we use the well known Riemannian framework never before used for point cloud matching, and present a novel matching algorithm. We pose point set matching under rigid and non-rigid transformations in this framework and solve for the transformations using standard nonlinear optimization techniques. Finally, to evaluate the performance of our algorithm -- dubbed SDTM -- we present several synthetic and real data examples along with extensive comparisons to state-of-the-art techniques. The experiments show that our algorithm outperforms state-of the-art point set registration algorithms on many quantitative metrics.

Tractography refers to the process of tracing out the nerve fiber bundles from diffusion Magnetic Resonance Images (dMRI) data acquired either in vivo or ex-vivo. Tractography is a mature research topic within the field of diffusion MRI analysis, nevertheless, several new methods are being proposed on a regular basis thereby justifying the need, as the problem is not fully solved. Tractography is usually applied to the model (used to represent the diffusion MR signal or a derived quantity) reconstructed from the acquired data. Separating shape and orientation of these models was previously shown to approximately preserve diffusion anisotropy (a useful bio-marker) in the ubiquitous problem of interpolation. However, no further intrinsic geometric properties of this framework were exploited to date in literature. In this paper, we propose a new intrinsic recursive filter on the product manifold of shape and orientation. The recursive filter, dubbed IUKFPro, is a generalization of the unscented Kalman filter (UKF) to this product manifold. The salient contributions of this work are: (1) A new intrinsic UKF for the product manifold of shape and orientation. (2) Derivation of the Riemannian geometry of the product manifold. (3) IUKFPro is tested on synthetic and real data sets from various tractography challenge competitions. From the experimental results, it is evident that IUKFPro performs better than several competing schemes in literature with regards to some of the error measures used in the competitions and is competitive with respect to others.

Recursive subdivision on an initial control mesh generates a visually pleasing smooth surface in the limit. Nevertheless, users must carefully specify the initial mesh and/or painstakingly manipulate the control vertices at different levels of subdivision hierarchy to satisfy a diverse set of functional requirements and aesthetic criteria in the limit shape. This modeling drawback results from the lack of direct manipulation tools for the limit geometric shape. To improve the efficiency of interactive geometric modeling and engineering design, in this paper we integrate novel physics-based modeling techniques with powerful geometric subdivision principles, and develop a unified finite element method (FEM)-based methodology for arbitrary subdivision schemes. Strongly inspired by the recent research on Dynamic Non-Uniform Rational B-Splines (D-NURBS), we formulate and develop a dynamic framework that permits users to directly manipulate the limit surface obtained from any subdivision procedure via simulated “force” tools. The most significant contribution of our unified approach is the formulation of the limit surface of an arbitrary subdivision scheme as being composed of a single type of novel finite element. The specific geometric and dynamic features of our subdivision-based finite elements depend on the subdivision scheme used. We present our novel FEM for the modified butterfly and Catmull–Clark subdivision schemes, and generalize our dynamic framework to be applicable to other subdivision schemes. Our FEM-based approach significantly advances the state-of-the-art in physics-based geometric modeling since it provides a universal physics-based framework for any subdivision scheme. In addition, we systematically devise a mechanism that allows users to directly (not via control meshes) deform any subdivision surface; finally, we represent the limit surface of any subdivision scheme using a collection of subdivision-based novel finite elements. Our experiments demonstrate that the new unified FEM-based framework not only promises a greater potential for subdivision techniques in solid modeling, finite element analysis, and engineering design, but that it will further foster the applicability of subdivision geometry in a wide range of visual computing applications such as visualization, virtual reality, computer graphics, computer vision, robotics, and medical imaging as well.

We introduce a geometric shape modeling scheme which allows for representation of global and local shape characteristics of an object. Geometric models are well-suited for representing global shapes without local detail, but we propose a scheme which represents global shapes with local detail and permits model shaping as well as topological changes via physics-based control. The scheme represents shapes by pedal curves and surfaces, i.e. the loci of the foot of perpendiculars to the tangents of a fixed curve/surface from a fixed point called the pedal point. By varying the location of the pedal point, one can synthesize a large class of shapes which exhibit both local and global deformations. We introduce physics-based control for shaping these geometric models by letting the pedal point vary and use a snake to represent the position of this varying point. The model, a "snake pedal", allows for interactive manipulation via forces applied to the snake. We develop a fast numerical iterative algorithm for shape recovery from image data using this scheme. The algorithm involves the Levenberg-Marquardt (LM) method in the outer loop for solving the global parameters and the alternating direction implicit (ADI) method in the inner loop for solving the local parameters of the model. The combination of the global and local scheme leads to an efficient numerical solution to the model fitting problem. We demonstrate the applicability of this modeling scheme via examples of shape synthesis and shape estimation from real image data.

The authors develop integrated techniques that unify physics based modeling with geometric subdivision methodology and present a scheme for dynamic manipulation of the smooth limit surface generated by the (modified) butterfly scheme using physics based "force" tools. This procedure based surface model obtained through butterfly subdivision does not have a closed form analytic formulation (unlike other well known spline based models), and hence poses challenging problems to incorporate mass and damping distributions, internal deformation energy, forces, and other physical quantities required to develop a physics based model. Our primary contributions to computer graphics and geometric modeling include: (1) a new hierarchical formulation for locally parameterizing the butterfly subdivision surface over its initial control polyhedron, (2) formulation of dynamic butterfly subdivision surface as a set of novel finite elements, and (3) approximation of this new type of finite elements by a collection of existing finite elements subject to implicit geometric constraints. Our new physics based model can be sculpted directly by applying synthesized forces and its equilibrium is characterized by the minimum of a deformation energy subject to the imposed constraints. We demonstrate that this novel dynamic framework not only provides a direct and natural means of manipulating geometric shapes, but also facilitates hierarchical shape and nonrigid motion estimation from large range and volumetric data sets using very few degrees of freedom (control vertices that define the initial polyhedron).

Visual reconstruction refers to extracting stable descriptions from visual data ( [1] ; A. Blake and A. Zisserman, Visual Reconstruction, MIT Press, Cambridge, MA, 1987). Visual reconstruction problems are commonly formulated in an optimization framework and normally require the optimization of nonconvex functions especially when discontinuity preserving image/shape recovery is the goal. Example problems include, image restoration, surface reconstruction, shape from shading etc. Most existing deterministic methods fail to reach the global optimum and lack the generality to incorporate reasonably complex interactions between boolean valued line process variables used for representing the presence (absence) of discontinuities. Stochastic methods for solving such problems e.g. the simulated annealing algorithm or variants thereof do achieve a global optimum but are plagued by slow convergence rates. In this paper, we present a new hybrid search algorithm as an efficient solution for achieving the global optimum of the nonconvex function derived from a Markov random field formulation which allows for incorporation of complex interactions between the line process variables to better constraint the line processes. In the hybrid search, for the stochastic part, we develop an informed genetic algorithm (GA) while employing an incomplete Cholesky preconditioned conjugate gradient algorithm ( [23] ; S.H. Lai and B.C. Vemuri, Robust and efficient algorithms for optical flow computation, in: Proceedings of the International Symposium on Computer Vision, Coral Gables, FL, 1995, pp. 455–460) for the deterministic part. Our informed GA consists of a reproduction operator and an informed mutation operator. The informed mutation operator exploits specific domain knowledge in the search and is accomplished by the Gibbs sampler. Our hybrid search algorithm is highly parallelizable and leads to a globally optimal solution. The performance of our algorithm is demonstrated via experimental results on the sparse data surface reconstruction and the image restoration problem.

In this paper, we present efficient algorithms for collision detection of arbitrarily shaped rigid moving objects in a variety of interactive as well as non-interactive environments. The algorithms primarily consist of two stages. The first stage involves finding candidate objects for possible collisions. The second stage involves detecting exact (within a prespecified tolerance) collision between these candidates. The primary data structure used in the algorithms is an octree. In the first stage, we build an octree for the enclosure containing the objects, which is used to detect possible collisions. Assuming spatial/temporal coherence i.e., that the particles move slowly or that the time sampling is fast enough, the average time complexity of this stage can be shown to be O(n) (excluding the time complexity for a one time octree construction), where n is the number of particles. In the second stage, we build a surface-octree for each object. If the objects are convex and assuming coherence, the expected time complexity to detect precise (within a prespecified tolerance) collision for each pair is a constant (excluding the time complexity for a one time surface-octree construction). Therefore, the overall expected time complexity for convex object collision detection is linear with respect to n. For the concave objects, complexity analysis is nontrivial to perform and instead we provide a very practical (almost linear time) algorithm. We apply our algorithms to particle flow simulations by simulating flow density conditions often arising in granular flows.

The capacitance matrix method has been widely used as an efficient numerical tool for solving the boundary value problems on irregular regions. Initially, this method was based on the Sherman--Morrison--Woodbury formula, an expression for the inverse of the matrix (A+UVT) with ${\bf A} \in \Re^{n \times n}$ and ${\bf U, V} \in \Re^{n \times p}$. Extensions of this method reported in literature have made restrictive assumptions on the matrices A and (A+UVT). In this paper, we present several theorems which are generalizations of the capacitance matrix theorem in \cite{Buzbee} and are suited for very general matrices A and (A +UVT). A generalized capacitance matrix algorithm is developed from these theorems and holds the promise of being applicable to more general problems; in addition, it gives ample freedom to choose the matrix A for developing very efficient numerical algorithms. We demonstrate the usefulness of our generalized capacitance matrix algorithm by applying it to efficiently solve large sparse linear systems in several computer vision problems, namely, the surface reconstruction and interpolation, and the shape from orientation problem.

In this paper, we present two very efficient and accurate algorithms for computing optical flow. The first is a modified gradient-based regularization method, and the other is an SSD-based regularization method. For the gradient-based method, to amend the errors in the discrete image flow equation caused by numerical differentiation as well as temporal and spatial aliasing in the brightness function, we propose to selectively combine the image flow constraint and a contour-based flow constraint into the data constraint by using a reliability measure. Each data constraint is appropriately normalized to obtain an approximate minimum distance (of the data point to the linear flow equation) constraint instead of the conventional linear flow constraint. These modifications lead to robust and accurate optical flow estimation. We propose an incomplete Cholesky preconditioned conjugate gradient algorithm to solve the resulting large and sparse linear system efficiently. Our SSD-based regularization method uses a normalized SSD measure (based on a similar reasoning as in the gradient-based scheme) as the data constraint in a regularization framework. The nonlinear conjugate gradient algorithm in conjunction with an incomplete Cholesky preconditioning is developed to solve the resulting nonlinear minimization problem. Experimental results on synthetic and real image sequences for these two algorithms are given to demonstrate their performance in comparison with competing methods reported in literature.

Recursive subdivision schemes have been extensively used in computer graphics, computer-aided geometric design, and scientific visualization for modeling smooth surfaces of arbitrary topology. Recursive subdivision generates a visually pleasing smooth surface in the limit from an initial user-specified polygonal mesh through the repeated application of a fixed set of subdivision rules. We present a new dynamic surface model based on the Catmull-Clark subdivision scheme, a popular technique for modeling complicated objects of arbitrary genus. Our new dynamic surface model inherits the attractive properties of the Catmull-Clark subdivision scheme, as well as those of the physics-based models. This new model provides a direct and intuitive means of manipulating geometric shapes, and an efficient hierarchical approach for recovering complex shapes from large range and volume data sets using very few degrees of freedom (control vertices). We provide an analytic formulation and introduce the "physical" quantities required to develop the dynamic subdivision surface model which can be interactively deformed by applying synthesized forces. The governing dynamic differential equation is derived using Lagrangian mechanics and the finite element method. Our experiments demonstrate that this new dynamic model has a promising future in computer graphics, geometric shape design, and scientific visualization.

Surface smoothing applied to range-elevation data acquired using a variety of sources has been a very active area of research in computational vision over the past decade. Generalized splines have emerged as the single most popular approximation tool to this end. In this paper we present a new and fast algorithm for solving the surface smoothing problem using a membrane, a thin-plate, or a thin-plate—membrane spline for data containing discontinuities. Our approach involves imbedding the surface smoothing problem specified on an irregular domain (in the sense of discontinuities and boundaries) in a rectangular region using the capacitance-matrix method based on the Sherman-Morrison-Woodbury formula of matrix analysis. This formula is used in converting the problem of solving the original linear system resulting from a finite-element discretization of the variational formulation of the surface smoothing problem to solving a Lyapunov matrix equation or a cascade of two Lyapunov matrix equations. The reduced problem can then be solved very efficiently using the ADI method and the bi-conjugate-gradient technique. Our solution requires the generation of a dense capacitance matrix, for which we propose a practical and efficient method. We demonstrate the efficiency of our algorithm via experiments on sparse-data surface smoothing, making performance comparisons with the conjugate-gradient and hierarchical-basis preconditioned conjugate-gradient algorithms.

In this paper, we present new and fast numerical algorithms for shape recovery from brain MRI using multiresolution hybrid shape models. In this modeling framework, shapes are represented by a core rigid shape characterized by a superquadric function and a superimposed displacement function which is characterized by a membrane spline discretized using the finite-element method. Fitting the model to brain MRI data is cast as an energy minimization problem which is solved numerically. We present three new computational methods for model fitting to data. These methods involve novel mathematical derivations that lead to efficient numerical solutions of the model fitting problem. The first method involves using the nonlinear conjugate gradient technique with a diagonal Hessian preconditioner. The second method involves the nonlinear conjugate gradient in the outer loop for solving global parameters of the model and a preconditioned conjugate gradient scheme for solving the local parameters of the model. The third method involves the nonlinear conjugate gradient in the outer loop for solving the global parameters and a combination of the Schur complement formula and the alternating direction-implicit method for solving the local parameters of the model. We demonstrate the efficiency of our model fitting methods via experiments on several MR brain scans.

Several problems in early vision have been formulated in the past in a regularization framework. These problems, when discretized, lead to large sparse linear systems. In this paper, we present a novel physically based adaptive preconditioning technique which can be used in conjunction with a conjugate gradient algorithm to dramatically improve the speed of convergence for solving the aforementioned linear systems. A preconditioner, based on the membrane spline, or the thin plate spline, or a convex combination of the two, is termed a physically based preconditioner for obvious reasons. The adaptation of the preconditioner to an early vision problem is achieved via the explicit use of the spectral characteristics of the regularization filter in conjunction with the data. This spectral function is used to modulate the frequency characteristics of a chosen wavelet basis, and these modulated values are then used in the construction of our preconditioner. We present the preconditioner construction for three different early vision problems namely, the surface reconstruction, the shape from shading, and the optical flow computation problems. Performance of the preconditioning scheme is demonstrated via experiments on synthetic and real data sets.

It is well known that the spatial frequency spectrum of membrane and thin plate splines exhibit self-affine characteristics and, hence, behave as fractals. This behavior was exploited in generating the constrained fractal surfaces, which were generated by using a Gibbs sampler algorithm in the work of Szeliski and Terzopoulos (1989). The algorithm involves locally perturbing a constrained spline surface with white noise until the spline surface reaches an equilibrium state. We introduce a fast generalized Gibbs sampler that combines two novel techniques, namely, a preconditioning technique in a wavelet basis for constraining the splines and a perturbation scheme in which, unlike the traditional Gibbs sampler, all sites (surface nodes) that do not share a common neighbor are updated simultaneously. In addition, we demonstrate the capability to generate arbitrary order fractal surfaces without resorting to blending techniques. Using this fast Gibbs sampler algorithm, we demonstrate the synthesis of realistic terrain models from sparse elevation data.

Symmetric positive-definite (SPD) matrices are ubiquitous in Computer Vision, Machine Learning and Medical Image Analysis. Finding the center/average of a population of such matrices is a common theme in many algorithms such as clustering, segmentation, principal geodesic analysis, etc. The center of a population of such matrices can be defined using a variety of distance/divergence measures as the minimizer of the sum of squared distances/divergences from the unknown center to the members of the population. It is well known that the computation of the Karcher mean for the space of SPD matrices which is a negatively-curved Riemannian manifold is computationally expensive. Recently, the LogDet divergence-based center was shown to be a computationally attractive alternative. However, the LogDet-based mean of more than two matrices can not be computed in closed form, which makes it computationally less attractive for large populations. In this paper we present a novel recursive estimator for center based on the Stein distance – which is the square root of the LogDet divergence – that is significantly faster than the batch mode computation of this center. The key theoretical contribution is a closed-form solution for the weighted Stein center of two SPD matrices, which is used in the recursive computation of the Stein center for a population of SPD matrices. Additionally, we show experimental evidence of the convergence of our recursive Stein center estimator to the batch mode Stein center. We present applications of our recursive estimator to K-means clustering and image indexing depicting significant time gains over corresponding algorithms that use the batch mode computations. For the latter application, we develop novel hashing functions using the Stein distance and apply it to publicly available data sets, and experimental results have shown favorable comparisons to other competing methods.

In this paper, we address the problem of video tracking using covariance descriptors constructed from simple features extracted from the given image sequence. Theoretically, this can be posed as a tracking problem in the space of ($n \times n$) symmetric positive definite (SPD) matrices denoted by $P_n$. A novel probabilistic dynamic model in $P_n$ based on Riemannian geometry and probability theory is presented in conjunction with a geometric (intrinsic) recursive filter for tracking a time sequence of SPD matrix measurements in a Bayesian framework. This newly developed filtering method can be used for the covariance descriptor updating problem in covariance tracking, leading to new and efficient video tracking algorithms. To show the accuracy and efficiency of our tracker in comparison to the state-of-the-art, we present synthetic experiments on $P_n$ and several real data experiments for tracking in video sequences.

Shape modeling is an important constituent of computer vision as well as computer graphics research. Shape models aid the tasks of object representation and recognition. This paper presents a new approach to shape modeling which retains some of the attractive features of existing methods, and overcomes some of their limitations. Our technique can be applied to model arbitrarily complex shapes, which include shapes with significant protrusions, and to situations where no a priori assumption about the object's topology is made. A single instance of our model, when presented with an image having more than one object of interest, has the ability to split freely to represent each object. This method is based on the ideas developed by Osher and Sethian to model propagating solid/liquid interfaces with curvature-dependent speeds. The interface (front) is a closed, nonintersecting, hypersurface flowing along its gradient field with constant speed or a speed that depends on the curvature. It is moved by solving a “Hamilton-Jacobi” type equation written for a function in which the interface is a particular level set. A speed term synthesized from the image is used to stop the interface in the vicinity of object boundaries. The resulting equation of motion is solved by employing entropy-satisfying upwind finite difference schemes. We also introduce a new algorithm for rapid advancement of the front using what we call a narrow-band update scheme. The efficacy of the scheme is demonstrated with numerical experiments on low contrast medical images.

International audience; Estimating diffusion tensors is an essential step in many applications -- such as diffusion tensor image (DTI) registration, segmentation and fiber tractography. Most of the methods proposed in the literature for this task are not simultaneously statistically robust and feature preserving techniques. In this paper, we propose a novel and robust variational framework for simultaneous smoothing and estimation of diffusion tensors from diffusion MRI. Our variational principle makes use of a recently introduced total Kullback-Leibler (tKL) divergence for DTI regularization. tKL is a statistically robust dissimilarity measure for diffusion tensors, and regularization by using tKL ensures the symmetric positive definiteness of tensors automatically. Further, the regularization is weighted by a non-local factor adapted from the conventional non-local means filters. Finally, for the data fidelity, we use the nonlinear least-squares term derived from the Stejskal-Tanner model. We present experimental results depicting the positive performance of our method in comparison to competing methods on synthetic and real data examples.

In this paper, we present a novel dictionary learning framework for data lying on the manifold of square root densities and apply it to the reconstruction of diffusion propagator (DP) fields given a multi-shell diffusion MRI data set. Unlike most of the existing dictionary learning algorithms which rely on the assumption that the data points are vectors in some Euclidean space, our dictionary learning algorithm is designed to incorporate the intrinsic geometric structure of manifolds and performs better than traditional dictionary learning approaches when applied to data lying on the manifold of square root densities. Non-negativity as well as smoothness across the whole field of the reconstructed DPs is guaranteed in our approach. We demonstrate the advantage of our approach by comparing it with an existing dictionary based reconstruction method on synthetic and real multi-shell MRI data.

Shape modeling is an important constituent of computer vision as well as computer graphics research. Shape models aid the tasks of object representation and recognition. This paper presents a new approach to shape modeling which retains some of the attractive features of existing methods and overcomes some of their limitations. The authors' techniques can be applied to model arbitrarily complex shapes, which include shapes with significant protrusions, and to situations where no a priori assumption about the object's topology is made. A single instance of the authors' model, when presented with an image having more than one object of interest, has the ability to split freely to represent each object. This method is based on the ideas developed by Osher and Sethian (1988) to model propagating solid/liquid interfaces with curvature-dependent speeds. The interface (front) is a closed, nonintersecting, hypersurface flowing along its gradient field with constant speed or a speed that depends on the curvature. It is moved by solving a "Hamilton-Jacobi" type equation written for a function in which the interface is a particular level set. A speed term synthesized from the image is used to stop the interface in the vicinity of object boundaries. The resulting equation of motion is solved by employing entropy-satisfying upwind finite difference schemes. The authors present a variety of ways of computing the evolving front, including narrow bands, reinitializations, and different stopping criteria. The efficacy of the scheme is demonstrated with numerical experiments on some synthesized images and some low contrast medical images. >

The unscented Kalman filter (UKF) was recently introduced in literature for simultaneous multi-tensor estimation and tractography. This UKF however was not intrinsic to the space of diffusion tensors. Lack of this key property leads to inaccuracies in the multi-tensor estimation as well as in tractography. In this paper, we propose an novel intrinsic unscented Kalman filter (IUKF) in the space of symmetric positive definite matrices, which can be used for simultaneous recursive estimation of multi-tensors and tractography from diffusion weighted MR data. In addition to being more accurate, IUKF retains all the advantages of UKF for instance, multi-tensor estimation is only performed in the places where it is needed for tractography, which would be much more efficient than the two stage process involved in methods that do tracking post diffusion tensor estimation. The accuracy and effectiveness of the proposed method is demonstrated via real data experiments.

Computation of the mean of a collection of symmetric positive definite (SPD) matrices is a fundamental ingredient of many algorithms in diffusion tensor image (DTI) processing. For instance, in DTI segmentation, clustering, etc. In this paper, we present novel recursive algorithms for computing the mean of a set of diffusion tensors using several distance/divergence measures commonly used in DTI segmentation and clustering such as the Riemannian distance and symmetrized Kullback-Leibler divergence. To the best of our knowledge, to date, there are no recursive algorithms for computing the mean using these measures in literature. Recursive algorithms lead to a gain in computation time of several orders in magnitude over existing non-recursive algorithms. The key contributions of this paper are: (i) we present novel theoretical results on a recursive estimator for Karcher expectation in the space of SPD matrices, which in effect is a proof of the law of large numbers (with some restrictions) for the manifold of SPD matrices. (ii) We also present a recursive version of the symmetrized KL-divergence for computing the mean of a collection of SPD matrices. (iii) We present comparative timing results for computing the mean of a group of SPD matrices (diffusion tensors) depicting the gains in compute time using the proposed recursive algorithms over existing non-recursive counter parts. Finally, we also show results on gains in compute times obtained by applying these recursive algorithms to the task of DTI segmentation.

3D shape modeling has received enormous attention in computer graphics and computer vision over the past decade. Several shape modeling techniques have been proposed in literature, some are local (distributed parameter) while others are global (lumped parameter) in terms of the parameters required to describe the shape. Hybrid models that combine both ends of this parameter spectrum have been in vogue only recently. However, they do not allow a smooth transition between the two extremes of this parameter spectrum. We introduce a new shape-modeling scheme that can transform smoothly from local to global models or vice versa. The modeling scheme utilizes a hybrid primitive called the deformable superquadric constructed in an orthonormal wavelet basis . The multiresolution wavelet basis provides the power to continuously transform from local to global shape deformations and thereby allow for a continuum of shape models—from those with local to those with global shape descriptive power—to be created. The multiresolution wavelet basis allows us to generate fractal surfaces of arbitrary order that can be useful in describing natural detail. We embed these multiresolution shape models in a probabilistic framework and use them for recovery of anatomical structures in the human brain from MRI data. A salient feature of our modeling scheme is that it can naturally allow for the incorporation of prior statistics of a rich variety of shapes. This stems from the fact that, unlike other modeling schemes, in our modeling, we require relatively few parameters to describe a large class of shapes.

This paper proposes a novel method for computing linear basis images from tensor-valued image data. As a generalization of the nonnegative matrix factorization, the proposed method aims to approximate a collection of diffusion tensor images using nonnegative linear combinations of basis tensor images. An efficient iterative optimization algorithm is proposed to solve this factorization problem. We present two applications: the DTI segmentation problem and a novel approach to discover informative and common parts in a collection of diffusion tensor images. The proposed method has been validated using both synthetic and real data, and experimental results have shown that it offers a competitive alternative to current state-of-the-arts in terms of accuracy and efficiency.

In this paper we present a new approach for the non-rigid registration of multi-modality images. Our approach is based on an information theoretic measure called the cumulative residual entropy (CRE), which is a measure of entropy defined using cumulative distributions. Cross-CRE between two images to be registered is defined and maximized over the space of smooth and unknown non-rigid transformations. For efficient and robust computation of the non-rigid deformations, a tri-cubic B-spline based representation of the deformation function is used. The key strengths of combining CCRE with the tri-cubic B-spline representation in addressing the non-rigid registration problem are that, not only do we achieve the robustness due to the nature of the CCRE measure, we also achieve computational efficiency in estimating the non-rigid registration. The salient features of our algorithm are: (i) it accommodates images to be registered of varying contrast+brightness, (ii) faster convergence speed compared to other information theory-based measures used for non-rigid registration in literature, (iii) analytic computation of the gradient of CCRE with respect to the non-rigid registration parameters to achieve efficient and accurate registration, (iv) it is well suited for situations where the source and the target images have field of views with large non-overlapping regions. We demonstrate these strengths via experiments on synthesized and real image data.

This paper presents a novel and robust technique for group-wise registration of point sets with unknown correspondence. We begin by defining a Havrda-Charvát (HC) entropy valid for cumulative distribution functions (CDFs) which we dub the HC Cumulative Residual Entropy (HC-CRE). Based on this definition, we propose a new measure called the CDF-HC divergence which is used to quantify the dis-similarity between CDFs estimated from each point-set in the given population of point sets. This CDF-HC divergence generalizes the CDF based Jensen-Shannon (CDF-JS) divergence introduced earlier in the literature, but is much simpler in implementation and computationally more efficient.A closed-form formula for the analytic gradient of the cost function with respect to the non-rigid registration parameters has been derived, which is conducive for efficient quasi-Newton optimization. Our CDF-HC algorithm is especially useful for unbiased point-set atlas construction and can do so without the need to establish correspondences. Mathematical analysis and experimental results indicate that this CDF-HC registration algorithm outperforms the previous group-wise point-set registration algorithms in terms of efficiency, accuracy and robustness.

This paper exploits the power of optimal sampling lattices in tomography based reconstruction of the diffusion propagator in diffusion weighted magnetic resonance imaging (DW-MRI). Optimal sampling leads to increased accuracy of the tomographic reconstruction approach introduced by Pickalov and Basser [1]. Alternatively, the optimal sampling geometry allows for further reducing the number of samples while maintaining the accuracy of reconstruction of the diffusion propagator. The optimality of the proposed sampling geometry comes from the information theoretic advantages of sphere packing lattices in sampling multidimensional signals. These advantages are in addition to those accrued from the use of the tomographic principle used here for reconstruction. We present comparative results of reconstructions of the diffusion propagator using the Cartesian and the optimal sampling geometry for synthetic and real data sets.

In this paper, we present a novel variational formulation for restoring high angular resolution diffusion imaging (HARDI) data. The restoration formulation involves smoothing signal measurements over the spherical domain and across the 3D image lattice. The regularization across the lattice is achieved using a total variation (TV) norm based scheme, while the finite element method (FEM) was employed to smooth the data on the sphere at each lattice point using first and second order smoothness constraints. Examples are presented to show the performance of the HARDI data restoration scheme and its effect on fiber direction computation on synthetic data, as well as on real data sets collected from excised rat brain and spinal cord.