Non-parametric diffusion tensor distribution magnetic resonance imaging

Abstract

This method uses magnetic resonance (MR) data to estimate a diffusion tensor distribution (DTD). The method includes inverting a Fredholm integral of the first kind, specifically performing an nD Inverse Laplace Transform subject to several constraints. First, a positive definiteness constraint is used, which zeros out a subset of diffusion tensor components that do not lie on a manifold of symmetric positive definite matrices. Second, marginal distribution constraints are used, which further partition the manifold of symmetric positive definite matrices into even smaller domains to improve DTD estimates in each MR voxel.

Description

[0001]Leydig 774370 HHS E-068-2023-0-US-01 NON-PARAMETRIC DIFFUSION TENSOR DISTRIBUTION MAGNETIC RESONANCE IMAGING CROSS-REFERENCE TO RELATED APPLICATIONS This application claims priority to U.S. Provisional Patent Application No.63 / 698,918, filed on September 25, 2024, which is incorporated herein by reference in its entirety. STATEMENT OF GOVERNMENT SUPPORT This invention was made with Government support under project number Z01#: 1 ZIAHD008971-04 by the National Institutes of Health, Eunice Kennedy Shriver National Institute of Child Health and Human Development (NICHD). The United States Government has certain rights in the invention. FIELD This disclosure pertains to magnetic resonance (MR) based determination of a joint distribution of diffusion tensors and diffusion tensor distribution magnetic resonance imaging (DTD-MRI) for characterizing diffusion properties of a subject system. BACKGROUND Accurately measuring and mapping nervous system tissue microstructure noninvasively has been a long sought-after goal in neuroscience. Clinically, having a clear understanding of patient’s nervous tissue microstructure provides powerful information in assessing many neuropathologies, such as cancer and stroke, which are associated with changes in nervous tissue microstructure. Brain tissue microstructure plays a critical role in ensuring normal function and subtle changes in it can serve as an early biomarker for emerging pathology. However, brain structure (e.g., nervous tissue) is challenging to study and characterize in vivo because it is heterogeneous over a wide range of spatial length scales from nanometers to centimeters. See, e.g., Lichtman Jeff W., Denk Winfried, “The big and the small: Challenges of imaging the brain’s circuits,” Science 2011, 334(6056):618–623 (the entire contents of which are hereby incorporated by reference herein). It has been recognized recently that whole brain imaging at an intermediate “mesoscopic” scale, consisting of functional neuronal units, may provide valuable new insights into normal and abnormal brain function. See, e.g., Bohland Jason W. et al, “A Proposal for a Coordinated Effort for the Determination of Brainwide Neuroanatomical Connectivity in Model Organisms at a Mesoscopic Scale,” PLoS Computational Biology 2009, 5(3) (the entire contents of which are hereby incorporated by reference herein). For instance, super-resolution microscopy has allowed Leydig 774370 HHS E-068-2023-0-US-01 live cell imaging of extracellular spaces, axon morphology, sub-cellular structures, etc., which has been a critical step towards achieving this objective. See, e.g., Tønnesen Jan et al., “Super- Resolution Imaging of the Extracellular Space in Living Brain Tissue,” Cell 2018, 172(5):1108– 1121; Ch´ereau Ronan et al., “Super resolution imaging reveals activity-dependent plasticity of axon morphology linked to changes in action potential conduction velocity,” Proceedings of the National Academy of Sciences of the United States of America 2017, 114(6):1401–1406; and Sigal Yaron M. et al., “Visualizing and discovering cellular structures with super-resolution microscopy,” Science 2018, 361(6405):880–887 (the entire contents of each of which are hereby incorporated by reference herein). However, this and other optical microscopy methods have a limited field of view (FOV), are subject to light loss due to scattering and absorption, and visualization often requires exogenous dyes, contrast agents, or genetic markers, and hence are limited to organotypic and other cell cultures, organoids, or small biological model systems (e.g., ranging from C. elegans to rodents). On the other hand, magnetic resonance imaging (MRI) is uniquely suited for ex vivo and in vivo whole mammalian (e.g., human) brain clinical imaging given its non-invasiveness and large FOV. The sensitivity of MRI, however, is poor in comparison to many microscopic techniques, owing to its limited nominal spatial resolution of approximately 8 mm3voxels in clinical MRI scanners, which averages over a myriad of cells and processes. Diffusion MRI methods can overcome these sensitivity and resolution limits by using water molecules, abundant in tissue, as a probe of tissue microstructure. The advent of strong, rapidly switching magnetic field gradients (i.e., pulsed field gradients (PFG)) has the potential to obtain whole brain scans in vivo at a mesoscopic length scale, in a clinically feasible time frame. See, e.g., Foo, Thomas K. et al., “Highly efficient head-only magnetic field insert gradient coil for achieving simultaneous high gradient amplitude and slew rate at 3.0T (MAGNUS) for brain microstructure imaging,” Magnetic Resonance in Medicine 2020, 83(6):2356–2369; and Huang, Susie Y. et al., “Connectome 2.0: Developing the next generation ultra-high gradient strength human MRI scanner for bridging studies of the micro-, meso- and macro-connectome,” NeuroImage 2021, 243:118530 (the entire contents of each of which are hereby incorporated by reference herein). However, there remains a need to develop accompanying models to relate the measured macroscopic diffusion weighted (DW) signals to the underlying mesoscopic or microscopic tissue structure. Over the past few decades, several models have been proposed for this purpose. See, e.g., Novikov, Dmitry S. et al., “On modeling,” Magnetic Resonance in Medicine 2018, 79(6):3172– 3193 (the entire contents of which is hereby incorporated by reference herein). A general framework that has the potential to describe whole brain mesoscopic tissue microstructure is the Leydig 774370 HHS E-068-2023-0-US-01 diffusion tensor distribution (DTD) model. See, e.g., Jian, Bing et al., “A novel tensor distribution model for the diffusion-weighted MR signal,” NeuroImage 2007, 37(1):164–176; Magdoom, Kulam Najmudeen et al., “A new framework for MR diffusion tensor distribution,” Scientific Reports.2021;11(1):2766 (“Magdoom et al.”); and Westin, Carl-Fredrik et al, “Q-space trajectory imaging for multidimensional diffusion MRI of the human brain,” NeuroImage 2016, 135:345–362 (“Westin et al.”) (the entire contents of each of which are hereby incorporated by reference herein). The DTD model assumes tissues are composed of an ensemble of water compartments (or pools) each exhibiting anisotropic Gaussian diffusion, whose diffusion tensor ellipsoids may vary in size, shape, and orientation, as described by a probability density function (pdf) of diffusion tensors within a voxel, p(D) (i.e., the DTD). These sub-compartments, which can be viewed as microvoxels, are akin to the myriad of microscopic compartments observed in electron microscopy (EM) images of brain tissue, where water is assumed to be undergoing free or hindered diffusion. This assumption at the mesoscopic length and time scales is reasonable given the finite water permeability across lipid plasma membranes via passive diffusion across the lipid bilayer and the presence of water channels, such as aquaporins, and facilitated diffusion, which translocates water via the action of ion pumps. See, e.g., MacAulay, N. et al., “Water transport in the brain: Role of cotransporters,” Neuroscience 2004, 129(4):1029–1042 (the entire contents of which are hereby incorporated by reference herein). It should be noted, however, that the DTD model is not limited to biological systems, but by natural extension, the DTD model can be adapted to characterize diffusion transport in a system where a solvent or solute moves between compartments in the system at least in part by diffusion. In addition, owing to the action of the glymphatic system in the brain, cerebrospinal fluid (e.g., water) is continuously cycled or circulated rather than being restricted within closed compartments to prevent build-up of toxins in the brain parenchyma. See, e.g., Xie Lulu et al. “Sleep Drives Metabolite Clearance from the Adult Brain,” Science 2013, 342(6156):373–377 (the entire contents of which are hereby incorporated by reference). Experimental evidence affirms Gaussian diffusion in brain parenchyma at mesoscopic length and time scales given that apparent diffusion coefficient (ADC) reaches a constant asymptotic value at diffusion times longer than 10 ms, rather than continuously decreasing with diffusion time—which is characteristic of water being confined (i.e., restricted) in compartments. Multiple-PFG signals, which are sensitive to mesoscopic heterogeneity in live human brain is independent of the diffusion time (clinically feasible) further strengthening the argument for mesoscopic Gaussian diffusion. See, e.g., Magdoom, et al. "Water Diffusion in the Live Human Brain is Gaussian at the Mesoscale." bioRxiv (2024): 2024-04 (the entire contents of which are Leydig 774370 HHS E-068-2023-0-US-01 hereby incorporated by reference). Thus, imaging the intravoxel DTD non-parametrically has the potential to reveal various mesoscopic structural domains or water pools present within a voxel. See e.g., Clark, Chris A. et al., “Diffusion time dependence of the apparent diffusion tensor in healthy human brain and white matter disease,” Magnetic Resonance in Medicine 2001;45(6):1126–1129; and van Gelderen, Peter et al., “Water diffusion and acute stroke,” Magnetic Resonance in Medicine 1994, 31(2):154–163 (the entire contents of each of which are hereby incorporated by reference herein). However, inversion of the diffusion weighted MRI (DWI) signal data to obtain the non-parametric DTD is challenging for several reasons, including: 1) the inversion requires taking an inverse multidimensional Laplace transform (ILT), which is a well-known, ill-posed, ill-conditioned problem, and 2) the number of unknowns in the high dimensional space of diffusion tensors (i.e., 6D corresponding to the six independent components of the second order tensor) is several orders of magnitude larger than the number of experimental DWI data that can be acquired practically in a clinical time period. Several studies have attempted to estimate features of the DTD, each having its attendant limitations. For example, Westin et al. used a cumulant expansion to express the MR signal in terms of a Gaussian DTD whose 2nd-order mean, and 4th-order (isotropic) covariance tensors are estimated from the MR signal profile. Like diffusion kurtosis imaging (DKI), this method can only be used over a narrow range of b-values. Topgaard estimated a reduced 4D DTD in a phantom assuming cylindrical symmetry of the underlying microdiffusion tensors. To describe heterogeneity and microscopic anisotropy in the cortex, Avram et al. estimated 3D DTDs assuming the eigenvectors of mean diffusion tensors were coincident with the cortical (i.e., developmental) reference frame. Magdoom et al. estimated the constrained normal DTD because it results in the maximum entropy distribution and applicable for all b-values. See Westin et al. (cited above); Topgaard Daniel, “Diffusion tensor distribution imaging,” NMR in Biomedicine.2019, 32(5):e4066; Avram Alexandru V. et al., “COnstrained Reference frame diffusion TEnsor Correlation Spectroscopic (CORTECS) MRI: A practical framework for high-resolution diffusion tensor distribution imaging,” Frontiers in Neuroscience.2022.16:2140; and Magdoom et al. (cited above) (the entire contents of each of which are incorporated by reference herein). Recently, Song et al., reported a non-parametric estimation of the full DTD using only single pulsed field gradient (sPFG) DWI data. See Song, Yiqiao et al., “Measurement of Full Diffusion Tensor Distribution Using High-Gradient Diffusion MRI and Applications in Diffuse Gliomas,” Frontiers in Physics 2022, 0:196 (“Song et al.”) (the entire contents of which are hereby incorporated by reference herein). The approach of Song et al., however, is potentially subject to the ill-posedness of the 6D Leydig 774370 HHS E-068-2023-0-US-01 ILT. A promising algorithm for performing a non-parametric multi-dimensional ILTs was first introduced by Benjamini and Basser, where the 1D marginal probability density functions of the dependent variables of interest are estimated separately, and then used as constraints for the subsequent reconstruction of their nD joint density function, a technique known as marginal distribution constrained optimization (MADCO). See Benjamini Dan and Basser Peter J., “Use of marginal distributions constrained optimization (MADCO) for accelerated 2D MRI relaxometry and diffusometry,” Journal of Magnetic Resonance 2016, 271:40–45; see also U.S. Patent No. 11,415,652 B2 (the entire contents of each of which are hereby incorporated by reference herein). This approach resulted in vastly reduced data acquisition requirements, and increased precision and accuracy of the estimated joint distributions. Due to a lesser number of unknowns and their higher signal-to-noise ratio (SNR), the estimates of the marginal densities are more precise and accurate. This approach originally implemented for 2D diffusometry-relaxometry applications, and subsequently for joint 3D distribution of mean diffusivity, T1and T2relaxation times, is increasingly more powerful as dimensionality of the distribution increases given the non-linear increase in the volume of the parameter space of unknowns with the number of dimensions. However, extending this approach from estimating a joint distribution of elements of a vector, such as mean diffusivity and relaxation times, to a joint distribution of elements of a second-order tensor, such as in the DTD, is not straightforward. In this latter case the off-diagonal components of the tensor are inherently, statistically coupled to the diagonal ones, which makes estimating their individual marginal densities or distributions directly from DWI data more challenging. The methods and apparatuses disclosed herein address limitations of these prior approaches. SUMMARY Disclosed herein are methods and apparatuses that can be used to measure and map the 6D DTD in neural tissue or other specimens (including non-biological systems experiencing diffusion transport of aqueous or polymeric components between cells or other “containers”). Aspects of the present disclosure provide for robust, non-parametric estimation of the 6D DTD by systematically estimating the 6D DTD’s lower dimensional projections and using these lower dimensional projections as further constraints to build the full 6D joint DTD. The approach of the present disclosure significantly simplifies the 6D DTD estimation problem by recognizing that the diffusion tensor lies in manifold of positive semidefinite tensors, and by recognizing that the DTD within a voxel is sparse. Using lower dimensional projections of the 6D DTD, which can be estimated in a computationally efficient manner, along with the positive definiteness constraint, regions of high probability density within the parameter space can be Leydig 774370 HHS E-068-2023-0-US-01 identified, which significantly reduces the complexity of the estimation problem. Aspects of the present disclosure further provide methods and apparatuses that efficiently estimate—from a time and computational load standpoint—the lower dimensional projections using a set of single, double, and triple-pulsed field gradient (PFG) MRI measurements (or alternatively single, double, and triple diffusion encoded MRI sequences) and acquire diffusion weighted images (DWI) sensitive to arbitrary lower dimensional projections of the 6D DTD. For example, one implementation in accordance with the present disclosure provides a magnetic resonance imaging (MRI) method for estimating a diffusion tensor distribution (DTD), which is a function of diffusion tensor components that characterizes diffusion properties of a subject tissue. The MRI method includes inverting a Fredholm integral of the first kind, specifically performing an nD Inverse Laplace Transform, subject to a plurality of constraints. The constraints include: a symmetric positive semi definiteness constraint, which zeros out a subset of the discrete diffusion tensor components that do not lie on a manifold of symmetric positive semidefinite matrices; and marginal distribution constraint(s), which partition(s) the manifold of symmetric positive definite matrices into select regions where values of previously estimated marginal distributions of a joint distribution of the DTD are above a predefined threshold. The DTD that minimizes the difference between the model and the captured diffusion-weighted signals is determined as the optimally estimated DTD. Besides producing a mean diffusion tensor as in DTI, many novel imaging biomarkers can be readily calculated from a DTD implementation according to aspects of the present disclosure. Such imaging biomarkers include, inter alia, 1D marginal distributions obtained from the DTD for the mean diffusivity (MD), and microscopic fractional anisotropy (µFA) and a 3D microscopic orientation distribution function (µODF), which embodies the complex geometry of white matter pathways and other fibrous tissues that may occupy the same voxel. It is possible to probe properties of water diffusion in gray and white matter and cerebrospinal fluid (CSF) that may occupy a single voxel by studying their respective components or modes separately. Aspects of the present disclosure can be applied to improve research and diagnosis in the fields of neurology, neuroradiology, neurosurgery, oncology, brain parcellation, psychology, fMRI integration, developmental biology, etc. Just like DTI, someone skilled in the art can use DTD in clinical diagnosis, following normal and abnormal developmental trajectories, assessing diseases, disorders, trauma, and aging processes, but with higher specificity and sensitivity as compared to DTI. It can also be used in conjunction with therapeutic applications such as neuronavigation, neurosurgical planning, radiation planning, etc. While the present disclosure focuses on applications for measuring DTDs in the brain, this formalism can equally be well suited to studying heterogeneous, Leydig 774370 HHS E-068-2023-0-US-01 anisotropic tissues in the body, such as kidney medulla, cardiac, skeletal, and smooth muscle, ligaments, tendon, cartilage, peripheral nerves, spinal cord and intravertebral disc, inter alia. BRIEF DESCRIPTION OF THE DRAWINGS Aspects of the present disclosure will be described in even greater detail below based on the exemplary figures. The present disclosure is not limited to the exemplary embodiments. All features described and / or illustrated herein can be used alone or combined in different combinations in embodiments of the invention. The features and advantages of various embodiments of the present disclosure will become apparent by reading the following detailed description with reference to the attached drawings, which illustrate the following: FIGS.1 and 2 graphically illustrates a reduction of the solution space by constraints imposed according to aspects of the present disclosure; FIG.3 illustrates an example of b-tensor ellipses / ellipsoids used for estimating marginal distributions; FIG.4 illustrates an exemplary magnetic resonance imaging (MRI) system; FIG.5 illustrates exemplary PFG experiments for capturing diffusion weighted MRI signals; FIG.6 illustrates and exemplary method for capturing diffusion weighted MRI signals and estimating the DTD; FIGS.7A-B illustrate synthetic MR phantoms and the effect of signal-to-noise ratio (SNR) in the reconstruction of the synthetic phantoms; FIG.8 illustrates a voxel-wise DTD reconstruction in polyvinyl pyrrolidone (PVP) liquid phantom; FIG.9 illustrates another DTD reconstruction obtained in a representative slice of a PVP polymer solution phantom; and FIG.9 illustrates a representative computing and control environment for measuring and processing MRI signals and estimating the DTD. DETAILED DESCRIPTION An aspect of the present disclosure provides a model that estimates the full 6D DTD from MRI signals in a time and computationally efficient manner. As explained in detail below, this DTD estimator model provides robust, non-parametric estimation of the 6D DTD by systematically estimating the 6D DTD’s lower dimensional projections or marginal distributions and using these lower dimensional projections as constraints to reconstruct the full 6D joint DTD. A diffusion weighted MR signal, characterized by an ensemble of diffusion tensors Leydig 774370 HHS E-068-2023-0-US-01 distributed according to ^^(^^), may be expressed as: ^^(^^) = ^^0 ∫ ^^−^^:^^^^(^^)^^^^ (1),where ^^0 is the signal without diffusion weighting (i.e., for S(B=0)), ^^, ^^ are the second-order,symmetric b-tensor and the diffusion tensor, respectively, and “:” is the tensor dot product. See, e.g., Basser, Peter J. et al, “Estimation of the Effective Self-Diffusion Tensor from the NMR Spin Echo,” Journal of Magnetic Resonance, Series B 1994, 103(3):247–254; and Jian, Bing et al., “A novel tensor distribution model for the diffusion-weighted MR signal,” NeuroImage 2007, 37(1):164–176) (the entire contents of each of which are hereby incorporated by reference herein). It should be noted that the use of the Gaussian anisotropic diffusion kernel, ^^−^^:^^, explicitly assumes the subdomains within the voxel exhibit Gaussian diffusion. Equation (1) is related to the Fredholm integral of the first kind, which is the 6D Laplace Transform of p(D). The above signal model for a discrete set of diffusion tensors, ^^^^, and experimental b- tensors, ^^^^, is given by: ^^(^^^^) ≈ ^^0^^^^^^ ^^(^^^^) (2),where ^^^^^^ = exp(−^^^^ ∶ ^^^^) ^^^^^^ is the ^^ × ^^ anisotropic Gaussian diffusion kernel matrix,^^^^^^is the product of discretization step size of the independent components of the diffusion tensor, ^^^^, and ^^(^^^^)is an ^^ dimensional vector consisting of the probability density of the discretized components of the DTD. That is, the measured diffusion weighted signal is modeled (i.e., approximated) by ^^0^^^^^^^^(^^^^). Because direct inversion of the above signal model is ill-posed, an insight embedded in the DTD estimator model implemented herein is to use a set of physically, mathematically, and statistically-motivated constraints to reduce the degrees of freedom and the space of admissible solutions—thereby reducing the number of unknowns to estimate. First, the positive, semi-definiteness constraint is applied on the diffusion tensors by zeroing the DTD for those discrete diffusion tensor components in rank ℝ6, which do not lie on themanifold of second-order symmetric positive semidefinite matrices, ℳ+(i.e., ^^(^^) = 0, ∀ ^^ ∉ℳ+). In a study conducted, this has resulted in an approximately 64% reduction in the number of unknowns for the range of diffusivities used (numerically estimated). Second, a hierarchy of lower dimensional marginal distributions are used to constrain the 6D joint distribution estimation by partitioning ℳ+into regions (domains) with high probability density, while zeroing others, thereby vastly limiting the volume of the search space (e.g., limiting the remaining 36% of the solution space in ℳ+). The resulting optimization problem is given by: Leydig 774370 HHS E-068-2023-0-US-01 where ^^ is the Lagrange multiplier used to regularize the solution whose value is computed using the S-curve method (see, e.g., Fordham, E.J. et al., “Imaging Multiexponential Relaxation in the (y, LogeT1) Plane, with Application to Clay Filtration in Rock Cores,” Journal of Magnetic Resonance, Series A 1995113(2):139–150 (the entire contents of which are hereby incorporated by reference herein)), Δ^^^^,^^is the discretization step size used for computing the ithmarginal densityintroduced to normalize the constraints, is the ^^ × ^^ matrix integral operator, which maps agiven higher dimensional DTD to the marginal density of interest, ^^ is a tolerance variable vector introduced to relax the constraints due to the uncertainty in the estimation of marginal densities resulting from noise in the measurement, M is the number of marginals, and ^^′is the marginal density estimated independently as discussed below. The regions of interest in ℳ+were selected from the marginal densities to encompass 95% confidence intervals, which significantly reduces the rank of the Ψ and Φ matrices, making the problem computationally tractable. What is provided herein, therefore, is an estimation method of the 6D DTD that is smooth and minimizes the error between the measured and the corresponding signal generated by the DTD model for that b-tensor, subject to the constraints that a) the diffusion tensors are positive definite, b) their DTDs are non- negative, and c) the computed marginal densities from the estimated DTD agree with their corresponding measured densities. To further elucidate the impact of the constraints used herein in reducing the number of unknowns to estimate, FIGS.1 and 2 graphically illustrate a gradual pruning of the solution space by the above-discussed constraints for two cases of 3D DTDs. FIG.1 graphically illustrates the gradual pruning of the solution space by successive application of positive semi-definite and marginal density constraints illustrated for the estimation of p (Dxx, Dyy, Dxy). The region of interest (ROI) in the solution space (i.e., ℝ3) is shown shaded. The ground truth DTD is assumed to be three probability density spheres centered at Dxx = 0.2, 0.6, 3 µm2 / ms, Dyy= 1.4, 0.6, 3 µm2 / ms, Dxy= 0.2, 0.0, 0.0 µm2 / ms to represent white matter, gray matter, and CSF respectively subject to the ℳ+constraint. The vast 3D Cartesian grid of solution space is first reduced by the application of the ℳ+constraint as shown in FIG.1. A further reduction in the solution space is obtained by the application of 1D marginal density constraints for Dxx and Dyy. FIG.2 illustrates successive pruning of the space of admissible solutions by the marginal Leydig 774370 HHS E-068-2023-0-US-01 density constraints for the estimation of p(Dxx, Dyy, Dzz). The ROI in the solution space (i.e., ℝ3) is shown shaded. The ground truth DTD is assumed to be three spheres centered at Dxx= 0.2, 0.6, 3 µm2 / ms, Dyy = 0.2, 0.6, 3 µm2 / ms, Dzz = 1.4, 0.6, 3 µm2 / ms. to represent white matter, gray matter, and CSF respectively. Since all the diagonal elements of the diffusion tensor are physically required to be positive, the application of the ℳ+constraint does not result in a volume reduction of the 3D Cartesian grid of unknowns. However, the application of 1D and 2D marginal density constraints result in a significant reduction as shown in the FIG.2. The intersection of 1D marginal density constraints for Dxx, Dyy, Dzz shown in the second-row results in a smaller sculpted region in ℝ3. This is further refined by applying the 2D marginal constraints which results in a region very similar to the ground truth region as shown in the last row. Using the above framework, the full 6D DTD estimation is performed hierarchically as sixteen lower-dimensional sub-problems with many fewer unknowns. Namely, the full diffusion tensor distribution can be determined using three 1D, three 2D, four 3D, three 4D and three 5D marginal densities of the diffusion tensor (FIG.3). The diffusion tensor may be expressed as: ] Here, the diagonal elements represent the apparent diffusion coefficient along the three orthogonal directions, while the off-diagonal elements represent the correlation between the diffusivities in different directions. Furthermore, the diffusion tensor is known to be symmetric based on physical considerations. The 1D marginal densities of the three diagonal components of the diffusion tensor can be obtained from inverting the single-PFG measurements using ℓ2-regularized non-negative least squares (NNLS) approach. See, e.g., Fordham, E.J. et al., “Imaging Multiexponential Relaxation in the (y, LogeT1) Plane, with Application to Clay Filtration in Rock Cores,” Journal of Magnetic Resonance, Series A 1995113(2):139–150 (the entire contents of which are hereby incorporated by reference herein). Here, the gradients can be incremented individually along the three (x, y, z) axes. This allows for the resulting b-values along a given axis to be linearly spaced over a given interval (although other spacings can be used). As discussed below, the 1D marginal densities can then be used as constraints to successively reconstruct the 2D, 3D, 4D, and 5D marginal densities involving observable combinations of diffusion tensor components (e.g., p (Dxx, Dyy), p (Dxx, Dyy, Dxy), etc.) leading to the full 6D DTD. The lower dimensional projections of the 6D DTD can be estimated from DW signals that selectively encode the diffusion tensor components of interest and can be inverted to obtain the marginal densities. Leydig 774370 HHS E-068-2023-0-US-01 While single-PFG experiments are straightforward to implement and are signal-to-noise ratio (SNR) efficient, single-PFG experiments are not sufficient to measure correlations between and among different diffusion tensor elements because the resulting b-tensors do not fully span the ℳ+space (e.g., a diagonal b-tensor required to estimate the joint distribution of pairs of diagonal diffusion tensor elements (e.g., p(Dxx, Dyy)) cannot be realized with single-PFG data as it will inherently encode off-diagonal components). Single-PFG experiments are also not amenable for estimating the covariance tensor. Thus, the single-PFG experiments are used only to estimate the 1D marginal densities of the diagonal components. Double-PFG experiments are used to estimate the 2D / 3D marginal densities involving pairs of diagonal components, and triple-PFG experiments are used to estimate the 3D / 4D / 5D distributions involving all the diagonal components of the diffusion tensor. That is, due to the marginal densities of the off-diagonal tensor components not being directly observable, a set of double-PFG experiments is run to estimate the 2D joint density of pairs of diagonal components of the diffusion tensor, and then 3D joint density of all the coupledcomponents of the diffusion tensor (i.e., ^^(^^^^^^ , ^^^^^^, ^^^^^^), ^^(^^^^^^, ^^^^^^ , ^^^^^^), ^^(^^^^^^, ^^^^^^ , ^^^^^^))from diffusion weighted signals, which selectively encode the diffusion tensor components of interest. A set of triple-PFG experiments is used to estimate the 3D joint density of diagonalcomponents (i.e., ^^(^^^^^^ , ^^^^^^, ^^^^^^)) along with 4D (i.e., ^^(^^^^^^, ^^^^^^, ^^^^^^ , ^^^^^^), ^^(^^^^^^, ^^^^^^, ^^^^^^ , ^^^^^^, ^^^^^^), ^^(^^^^^^, ^^^^^^, ^^^^^^ , ^^^^^^ , ^^^^^^)) joint densities of all observablecomponents of the diffusion tensor. Selective encoding of the diffusion tensor components of the double- and triple-PFG experiment is achieved by zeroing the extraneous b-tensor components which are sampled randomly on ℳ+. The selectively encoded MR signal is inverted using an ℓ2-regularized non-negative least squares (NNLS) algorithm to obtain the marginal density of interest with the previously measured marginal densities used as constraints. Given the 16 marginal densities, the 6D joint density is estimated by solving Equations (3) and (4), above. The convex optimization problem for estimating the marginal and joint densities may be performed using any suitable convex optimization solver, such as CVXPY with MOSEK solver. See, e.g., Diamond Steven et al., “CVXPY: A Python-Embedded Modeling Language for Convex Optimization,” Journal of machine learning research 2016, 17:1–5; Agrawal, Akshay et al., “A rewriting system for convex optimization problems,” Journal of Control and Decision 2018, Leydig 774370 HHS E-068-2023-0-US-01 5(1):42–60; and MOSEK ApS . MOSEK Fusion API for Python.2023 (discussing CVXPY and MOSEK, respectively) (the entire contents of each of which are hereby incorporated by reference herein). A triple-PFG experiment may also be performed to sample the full 6D distribution. For the triple-PFG experiment, a set of 3×3 rank-3 matrices are chosen randomly on ℳ+to observe all the diffusion tensor components. The b-tensors for all these measurements are realized using the interfused-PFG (iPFG), which can efficiently generate b-tensors of all three ranks in a single spin- echo. The gradient strengths required to achieve a specified b-tensor are numerically calculated, such as in Magdoom 2023 et al. See Magdoom, KN et al., “A Novel Framework for In Vivo Diffusion Tensor Distribution MRI of the Human Brain, NeuroImage.2023 (the entire contents of which are hereby incorporated by reference herein). The 6D DTD space can be discretized into a Cartesian grid with the diagonal components of the diffusion tensor ranging from, e.g., 0 to 5 ^^^^2 / ^^^^ —to accommodate the range of apparent diffusion coefficients (ADCs) expected in live tissue (e.g., live human brain tissue, including to accommodate for the possibility of there being pseudo-diffusion in cerebrospinal fluid (CSF) leading to an ADC greater than for free water at 37° C, i.e., 3 ^^^^2 / ^^^^). The range of the off-diagonal components can be from, e.g.,−5 to + 5 ^^^^2 / ^^^^ based on the Sylvester’s criterion, or−2.5 to + 2.5 ^^^^2 / ^^^^ because anisotropically diffusing water pools are not expected outside thisrange in live normal human brain tissue. According to aspects of the present disclosure, the spectral resolution for all the components can be set relatively high, e.g., 0.2^^^^2 / ^^^^, while still being able to efficiently estimate the 6D DTD. Without the technology of the present disclosure, having such a high spectral reolustion, which results in approximately 2.3 billion unknowns needed to characterize the full 6D distrubtion, makes the calculation of the full 6D distribution computationaly intractable Such a high spectral resolution ordinarily would result in approximately 2.3 billion unknowns needed to characterize the full distribution (without any constraints), effectively making this problem computationally intractable. However, according the improved MRI method of the present disclosure, a total of three sets of rank-1 (n = 11 for estimating the 1D marginal densities), six sets of rank-2 (n = 21, 31 for estimating the 2D and 3D marginal densities respectively) and eight sets of rank-3 b-tensors (n = 31, 41, 51, 61 for estimating the 3D, 4D, 5D, and 6D marginal densities respectively) with b-values ranging from 0 to 2 ms / µm2can be used to estimate the full 6D DTD. Thus, in this implementation, a total of 11 × 3 + 21 × 3 + 31 × 4 + 41 × 3 + 51 × 3 + 61 = 557 b-tensors were used to reconstruct the full 6D DTD in each voxel (with initially 2.3 billion unknowns). Leydig 774370 HHS E-068-2023-0-US-01 For further elucidation of the hierarchical estimation used in embodiments of the present disclosure, FIG.3 illustrates an example of b-tensor ellipses / ellipsoids used for estimating the various marginal distributions. Additionally, in FIG 3, the hierarchical experimental design used to estimate the non-parametric 6D DTD is depicted by plots of sampled b-tensor ellipsoids for each experiment. In the top row, the rank 1 b-tensors 301 are shown within their associated planes, i.e., x- plane 302, y-plane 303, and z-plane 304. The rank 1 b-tensors 301 are displayed as sticks used to estimate the marginal density of diagonal components of the diffusion sensors—i.e., ^^(^^^^^^), In subsequent rows, these nD marginal distributions are used to reconstruct the appropriate (n+1) D marginal distributions shown using distinctly depicted arrows, ultimately leading to an estimate of the full 6D DTD. The marginal distribution estimated are indicated above each set of corresponding b-tensors. The shading in the enclosing box for each set of b-tensors indicates the plane of the b-tensor while the hatching in the plane indicates the rotation of the b-tensor along its normal vector. For example, the second row depicts the rank 2 b-tensors 311 as ellipses within their associated volume corresponding to a plane of the b-tensor, which is shaded. Thus, there are b- tensors provided as associated with the x-plane 312, y-plane 313, and z-plane 314. The rank 2 b- tensors for each of the associated planes, together with the estimated 1D marginal distributions ofthe diagonal components, are used to estimate the 2D marginal densities ^^(^^^^^^, ^^^^^^), ^^(^^^^^^ , ^^^^^^),^^(^^^^^^, As indicated by the arrows, the 2D marginal densities can then be used to estimate 3Djoint density of all the coupled components of the diffusion tensor (i.e., ^^(^^^^^^, ^^^^^^, ^^(^^^^^^, ^^^^^^ , ^^^^^^), ^^(^^^^^^ , ^^^^^^, ^^^^^^)), which are depicted in the third row. The hatching 315indicates the associated rotation of the b-tensor along its normal vector. The third row also depicts the 3D joint density of diagonal components (i.e., ^^^^^^, ^^^^^^)). This joint density contains rank 3 b-tensors 320, which are depicted as ellipsoidswithin their associated volume. As indicated by the shading of the three associated planes, and the arrows, this joint density is obtained by running a triple-PFG experiment, in addition to considering the 2D marginal densities. Continuing down, the fourth, fifth, and sixth rows respectively illustrate b-tensors for eachof the 4D, 5D, and 6D marginal densities (i.e., ^^(^^^^^^, ^^^^^^, ^^^^^^ , ^^^^^^), Leydig 774370 HHS E-068-2023-0-US-01 Accordingly, aspects of the present disclosure provide improved MRI methods and apparati that are capable of efficiently estimating the full 6D DTD from a vastly limited set of MRI experiments. The improvements provided herein enable the use of less computational power and time to estimate a full 6D DTD, with greater precision and accuracy, while also saving signal acquisition time. The improvements further advance the medical imaging field by providing a robust tool for examining a family of quantitative imaging biomarkers for use in researching and diagnosing diseases of the nervous system, inter alia. Exemplary MR Measurement Apparatus FIG.4 illustrates an exemplary MRI apparatus 1000 configured to obtaining MR measurements for use by the DTD estimator implemented according to aspects of the present disclosure. The apparatus 1000 includes a controller / interface 1002 that can be configured to apply selected magnetic fields, such as constant or pulsed field gradients, to a subject 1030 or other specimen. An axial magnet controller 1004 is in communication with an axial magnet 1006 that is generally configured to produce a substantially constant magnetic field B0. A gradient controller 1008 is configured to apply a constant or time-varying magnetic field gradient in one or more selected directions or in a set of directions using magnet coils 1010-1012 to produce respective magnetic field gradient vector components Gx, Gy, Gz or combinations thereof to produce the gradients G1, G2, G3 that are associated with b-matrices or tensors. A radiofrequency (RF) generator 1014 is configured to deliver one or more RF pulses to a specimen using a transmitter coil 1015. An RF receiver 1016 is in communication with a receiver coil 1018 and is configured to detect or measure net magnetization of spins. Typically, the RF receiver includes an amplifier 1016A and an analog-to-digital convertor 1016B that detect and digitized received signals to obtain the signals S(B). Slice selection gradients can be applied with the same hardware used to apply the diffusion gradients. The gradient controller 1008 can be configured to produce pulses or other gradient fields along one or more axes as needed for a particular b-matrix. By selection of such gradients and other applied pulses, signals associated with a plurality of b-matrices are acquired. For imaging, specimens are divided into volume elements (voxels) and MR signals for a plurality of gradient directions are acquired as discussed above, but signals can be acquired for one or only a few specimen voxels. In typical examples, signals are obtained for some or all voxels of interest. A computer 1024 or other processing system such as a personal computer, a workstation, a personal digital assistant, laptop computer, smart phone, or a networked computer can be provided for acquisition, control and / or analysis of specimen data. The computer 1024 generally includes a Leydig 774370 HHS E-068-2023-0-US-01 hard disk, a removable storage medium such as a floppy disk or CD-ROM, and other memory such as random access memory (RAM). Data can also be transmitted to and from a network using cloud- based processors and storage. Data could be uploaded to the Cloud or stored elsewhere. Computer- executable instructions for, e.g., data acquisition or control, b-matrix computations, such as random selected of directions and random selected of b-magnitudes as well as determining the associated b- matrices, and tensors and distribution estimations can be provided on a storage medium, such as a memory 1026, 1027, or delivered to the computer 1024 via a local area network, the Internet, or other network. Signal acquisition, instrument control, and signal analysis can be performed locally or with distributed processing. For example, signal acquisition and signal analysis can be performed at different locations. Signal evaluation can be performed remotely from signal acquisition by communicating stored data to a remote processor. In general, control and data acquisition with an MRI apparatus can be provided with a local processor, or via instruction and data transmission via a network. Exemplary Experimental Design FIG.5 illustrates an example PFG pulse sequence operation 100 that can be executed by an MRI apparatus (such as the exemplary MRI apparatus 1000) to generate suitable diffusion tensor weighted signals for use in estimating the full 6D DTD. As explained above, single-, double-, and triple-PFG experiments may be used to generate sufficient diffusion-weighted signals for use by a DTD estimator. FIG.5, in particular, illustrates examples of each of: (1) a single-PFG experiment—for estimating the marginal densities of the diagonal components of the diffusion tensor; (2) a double- PFG experiment—for estimating the 2D / 3D marginal densities of diffusion tensor components; and (3) a triple-PFG experiment—for sampling the 3D / 4D / 5D marginal densities and full 6D distribution. As a person of ordinary skill in the art would recognize, however, aspects of the present disclosure are not limited to the described PFG experiments. Aspects of the present disclosure may use any experimental approach for estimating the marginal density of diagonal components of the diffusion tensor, for estimating the nD marginal density, and the full 6D distribution. For example, the multi-pulsed field gradient (mPFG) experiments disclosed in International Patent Application Publication No. WO 2022 / 020102 A1 (the entire contents of which are hereby incorporated by reference herein) may be used. In magnetic resonance imaging (MRI), the marginal density of diagonal components of thediffusion tensor—i.e., ^^(^^^^^^), ^^(^^^^^^), ^^(^^^^^^) −may be measured using rank 1 b-tensors. Rank 1b-tensors can be realized from single pulse field gradient (PFG) experiments, which measure the diffusion properties of water molecules in a biological tissue. Leydig 774370 HHS E-068-2023-0-US-01 Generally, a strong magnetic field is applied to the tissue which aligns the spins of the water molecules. During a single-PFG spin echo experiment, a first RF signal is applied, and then, a gradient magnetic field is applied, which causes the precession frequencies of the spins to vary with position. A second RF signal is then applied, followed by a second gradient magnetic field. During the experiment, gradients are typically applied along three orthogonal directions, and the resulting signal is measured to obtain information about the diffusion properties of the tissue. FIG.5 illustrates an example sequence for a single-PFG spin echo experiment for measuring a diffusion weighted signal. In FIG.5, the radiofrequency (RF) pulse sequence 110 illustrates RF signals applied to the observed subject (e.g., by the RF generator and RF coil of an MRI apparatus), and the diffusion gradient pulse sequence 120 illustrates the single-pulse diffusion gradient sequence applied to the observed subject (e.g., by the gradient controller and gradient coil(s) of an MRI apparatus). In the spin echo single-PFG experiment, first, a 90° (e.g., x-axis) RF pulse 111 is applied to the subject to excite the spins of the water molecules. A gradient pulse 121 is then applied (e.g., along one orthogonal direction, such as along the x-axis), with an amplitude ^^1for a predefined duration δ . The application of the gradient pulse 121 causes a phase shift in the spins of the water molecules, which leads to an attenuation of the signal from the subject. Next, a 180° (e.g., y-axis) RF pulse 112, which inverts the phase of the spins. After applying the 180° RF pulse 112, a second gradient pulse 122 is then applied (e.g., along the same orthogonal direction) with an amplitude ^^1for a predefined duration δ . This allows the spins to rephase and generate an echo (e.g. after another time τ). The amplitude of the echo can be measured as a function of the b-value to determine the diffusion coefficient of the sample. The second gradient pulse 122 is applied a preset time after applying the first gradient pulse 121. This time is called the diffusion time Δ. The 180° RF pulse 112 is applied a predefined time τ after the 90° RF pulse 111, but within the diffusion time Δ. The attenuation in the signal is measured using an MRI receiver, and this measurement may be repeated for multiple gradient directions, strengths, duration, separation time between the gradient pulses, etc. The resulting data from the measured signals are then fitted to a mathematical model that relates the signal attenuation to the diffusion properties of the tissue. Specifically, the signal attenuation is modeled as: ^^(^^) = ^^0^^(−^^ ^^)(5), where ^^(^^)is the signal intensity for a given value of the b-value (b), ^^0is the signal intensity in the absence of diffusion weighting, and ^^ is the diffusion coefficient. The b-tensor is related to the b-value as follows: Leydig 774370 HHS E-068-2023-0-US-01 b= γ2 G2 δ2 (Δ − δ / 3) (6),where γ is the gyromagnetic ratio of water, G is the strength of the gradient, δ is the duration of the gradient pulse, and Δ is the time between the two gradient pulses (diffusion time). By measuring the signal attenuation at multiple b-values and gradient directions, it is possible to estimate the probability distribution of the marginal density along each axis (x, y, z). The estimation may be done by any means in the art, such as using a Bayesian framework to compute a posterior distribution, and then sampling from the posterior distribution using, e.g., Markov chain Monte Carlo (MCMC) methods to estimate the marginal densities. Here, ^^(^^^^^^), of the 6D DTD, which represent the probability density functions of thediffusion coefficients along the x, y, and z directions, respectively. These are each derived from rank 1 b-tensors. These functions provide information about the distribution of diffusion coefficients in the tissue, and can be used to infer the microstructure and organization of biological tissues. FIG.5 also illustrates an example double-PFG experiment used to estimate the 2D / 3D marginal densities, which also captures the off-diagonal components of the diffusion tensor. As explained above, off-diagonal components of the diffusion tensor are not directly observable by MRI experiments. However, a double-PFG experiment can capture the 2D / 3D marginal densities of the diffusion tensor components with rank 2 b-tensors. A rank 2 b-tensor is a matrix that describes the diffusion-weighting of the magnetic resonance signal as a function of the applied diffusion gradients, and it contains information about both the diagonal and off-diagonal elements of the diffusion tensor. The diagonal elements of the rank 2 b-tensor are proportional to the squared diffusion coefficients along the principal axes of diffusion (i.e., ^^^^^^, ^^^^^^, and ^^^^^^), while the off-diagonal elements are proportional to the cross-correlation between the diffusion in different directions (i.e., ^^^^^^, ^^^^^^, and ^^^^^^). By fitting a model to the signal attenuation data acquired at multiple diffusion gradient strengths, it is possible to estimate both the diagonal and off- diagonal elements of the diffusion tensor. Therefore, by applying the rank 2 b-tensor from a double pulsed-field gradient experiment, it is possible to obtain information about the off-diagonal elements of the DTD. FIG.5 illustrates an example sequence for a double-PFG experiment for measuring a non- diffusion weighted signal. Specifically, the PFG experiment shown is an interfused-PFG (iPFG) experiment. In FIG.4, the radiofrequency (RF) pulse sequence 110 illustrates RF signals applied to the observed subject (e.g., by the RF generator and RF coil of an MRI apparatus), and the diffusion gradient pulse sequence 130 illustrates the double-pulse diffusion gradient sequence applied to the observed subject (e.g., by the gradient controller and gradient coil(s) of an MRI apparatus). Leydig 774370 HHS E-068-2023-0-US-01 In the double-PFG experiment, a 90˚ (π / 2; x-axis) RF pulse 112 is applied, and followed by a set of diffusion weighting pulse pairs. The set of diffusion weighting pulse pairs includes a first pulse pair 131 and a second pulse pair 132 separated by a diffusion time, Δ. An 180˚ (π; y-axis) RF pulse 112 is applied during the diffusion time Δ between the first pulse pair 131 and the second pulse pair 132. Specifically, the 180˚ RF pulse 112 is applied after a time τ from the application of the 90˚ RF pulse 111. The first pulse pair 131 and the second pulse pair 132 are typically the same or substantially the same. For example, an initial pulse 134 of the first pulse pair 131 is typically the same or substantially the same as an initial pulse 136 of the second pulse pair 132, and a final pulse 135 of the first pulse pair 131 is the same or substantially the same as a final pulse 137 of the second pulse pair 132. In the example, each of the initial pulses 134, 136 have a magnitude G1, each of the final pulses 135, 137 have a magnitude G2, and all gradient pulses have the pulse duration δ. The initial and final pulses, however, can have various pulse shapes, magnitudes, durations, rise times, fall times, or other modulations. An initial pulse of a pulse pair need not have the same pulse characteristics as a final pulse of the pulse pair. While it is convenient that the first pulse pair and the second pulse pair as the same or substantially the same, the pulse pairs can be different but balanced or matched to provide diffusion sensitization and to compensate concomitant gradient fields. A final pulse in a pulse pair can be applied with an interfusing time τint delay after application of an initial pulse, although such delays are typically preferred to be much less than the diffusion time Δ. First and second pulse pairs are generally selected so that, absent diffusion, phase changes in spins produced by the first pulse pair are reversed or undone by the second pulse pair. Pulse orders can also be interchanged. Differing pulse orders are generally associated with different b-matrices. Typically, each pulse in a pulse pair can be associated with a gradient axis, with the initial and final pulses associated with different gradient axes. The diffusion weighting matrix (the b-matrix or b-tensor) for the set of diffusion gradient pulse pairs, is assuming rectangular gradient pulses: wherein is the b-value from diffusion gradient interactions between the ^^^^ℎand ^^^^ℎaxis, ^^(^^), ^^(^^) 12are the ^^^^ℎ-component of the first and second diffusion gradient pairs respectively in the pulse sequence, i, j are integers from 1, 2, or 3, and ^^ is a nuclear gyromagnetic ratio. The desired b-matrix of rank-2, ^^′, with b-value equal to ^^ is synthesized from two non-parallel unit vectors,^̂^1, ^̂^2, analogous to the two gradient vectors in the traditional dPFG experiment as shown below: Leydig 774370 HHS E-068-2023-0-US-01 wherein b is a diffusion weighting magnitude associated with the diffusing weighting b-matrix. The unit vectors can be randomly oriented uniformly over a sphere with b-value magnitude uniformly distributed over a desired range (via compressed sensing) to obtain the b-matrices (i.e., b-tensor) for DTD estimation. The physical diffusion gradient strengths in this pulse sequence required togenerate the desired b-matrices for a fixed ^^ and Δ are obtained by solving the equation ^^^^^^ = ^^^′^^^using a non-linear least square fitting routine, e.g., in MATLAB (Mathworks, Natick, MA), or otherwise solved. From the above-described double-PFG experiment diffusion weighted signals are obtained, from which the 2D joint density of pairs of diagonal components and 3D joint density of all thecoupled components of the diffusion tensor—i.e., ^^(^^^^^^, ^^^^^^, ^^^^^^), ^^(^^^^^^, ^^^^^^ , ^^^^^^),and ^^(^^^^^^, ^^^^^^ , ^^^^^^)—can be estimated. For example, a set (e.g., 50) of rank-2 b-tensors can beused to selectively encode the diffusion tensor components of interest by zeroing the extraneous b- tensor components, which effectively reduces them to a 2 × 2 rank-2 matrix that is sampled randomly on ℳ+. The selectively encoded rank 2 b-tensors can be inverted, e.g., using an ℓ2- regularized non-negative least squares approach, to obtain the 2D / 3D joint marginal density, with the measured densities (e.g., as previously determined) used as a constraint. FIG.5 also illustrates an example triple-PFG experiment used to sample the 3D / 4D / 5D and full 6D distribution. In the example shown, the triple-PFG experiment is an extension of the double-iPFG experiment discussed above to include a third gradient pulse such that it can generaterank 3 b-tensors for estimating the 3D, i.e., ^^(^^^^^^, ^^^^^^, ^^^^^^), 4D, i.e., ^^(^^^^^^, ^^^^^^, ^^^^^^ , ^^^^^^),^^(^^^^^^, ^^^^^^, ^^^^^^ , ^^^^^^), and ^^(^^^^^^, ^^^^^^, ^^^^^^ , ^^^^^^), 5D, i.e., ^^(^^^^^^, ^^^^^^ , ^^^^^^ , ^^^^^^, ^^^^^^),^^(^^^^^^, ^^^^^^, ^^^^^^ , ^^^^^^, ^^^^^^), and ^^(^^^^^^, ^^^^^^, ^^^^^^ , ^^^^^^ , ^^^^^^) and full 6D distribution, i.e.,^^(^^^^^^, ^^^^^^, ^^^^^^ , ^^^^^^,^^^^^^ , ^^^^^^). The rank-3 b-tensor is a 3×3 matrix which encodes the diffusionproperties in a small voxel (3D volume element) of the tissue. The diagonal element of the matrix encodes the diffusivity of water in the voxel along three orthogonal directions, while the off- diagonal elements encode the correlation between the diffusion properties in these directions. Overall, a rank 3 b-tensor obtained from a triple-PFG MRI experiment provides more detailed information about the diffusion properties of water in a tissue than a rank 2 b-tensor, as it provides information about higher order diffusion properties. In FIG.5, the radiofrequency (RF) pulse sequence 110 illustrates RF signals applied to the observed subject (e.g., by the RF generator and RF coil of an MRI apparatus), and the diffusion gradient pulse sequence 140 illustrates the triple-pulse diffusion gradient sequence applied to the Leydig 774370 HHS E-068-2023-0-US-01 observed subject (e.g., by the gradient controller and gradient coil(s) of an MRI apparatus). In the triple-PFG experiment, a 90˚ (π / 2; x-axis) RF pulse 112 is applied, and followed by a set of diffusion weighting pulses. The set of diffusion weighting pulse includes a first pulse triple 141 and a second pulse triple 142 separated by a diffusion time Δ. An 180˚ (π; y-axis) RF pulse 112 is applied during the diffusion time Δ between the first pulse triple 141 and the second pulse triple 142. Specifically, the 180˚ RF pulse 112 is applied after a time τ from the application of the 90˚ RF pulse 111. An echo occurs after a second time τ from the application of the 180˚ RF pulse 112. The first pulse triple 141 and the second pulse triple 142 are typically the same or substantially the same. For example, an initial pulse 143 of the first pulse triple 141 is typically the same or substantially the same as an initial pulse 146 of the second pulse pair 142, and similarly for the second pulses 144, 147, and final pulses 145, 148 of the pulse triples 141, 142. In the example, each of the initial pulses 143, 146 have a magnitude G1, each of the second pulses 144, 147 have a magnitude G2, each of the final pulses 145, 148 have a magnitude G3, and all gradient pulses have the pulse duration δ. Each of the pulses in the pulse triples are separated by an interfusing time τint. The initial, second, and final pulses, however, can have various pulse shapes, magnitudes, durations, rise times, fall times, or other modulations. An initial pulse of a pulse triple need not have the same pulse characteristics as a second or final pulse of the pulse triple. While it is convenient that the first pulse triple and the second pulse triple as the same or substantially the same, the pulse triples can be different but balanced or matched to provide diffusion sensitization and to compensate concomitant gradient fields. First and second pulse triples are generally selected so that, absent diffusion, changes in spins produced by the first pulse triples are reversed or undone by the second pulse triple. Pulse orders can also be interchanged. Differing pulse orders are generally associated with different b-matrices. Typically, each pulse in a pulse triple can be associated with a gradient axis, with the initial, second, and final pulses associated with different gradient axes. From the above-described triple-PFG experiment, a rank 3 b-tensor can be realized—e.g., by extending the above-described model to account for a third gradient pulse. Overall, sufficient triple-PFG experiments are run to provide a set (e.g., 100) or rank 3 b-tensors over a range of b- values randomly sampled on ℳ+. This set of rank 3 b-tensors provide samples of the full 6D DTD. This set of rank 3 b-tensors provide samples of the full 6D DTD. Exemplary MRI Method FIG.6 illustrates an exemplary magnetic resonance imaging (MRI) method for estimating a diffusion tensor distribution (DTD) 400. Leydig 774370 HHS E-068-2023-0-US-01 First, MR signal data is acquired by or otherwise received from an MRI apparatus operating on a subject tissue. In the method of FIG.6, an MRI apparatus executes single-PFG experiments to realize rank 1 b-tensors (S401). During the single-PFG experiments, gradients are incremented individually along three axis of the MRI apparatus such that resulting b-values along an axis are linearly spaced over an interval. The MRI apparatus further executes double-PFG experiments to realize rank 2 b-tensors (S402). The rank 2 b-tensors may be selectively encoded by zeroing extraneous b-tensor components. Finally, the MRI apparatus executes triple-PFG experiments to capture a set of rank 3 b-tensors (S403). The rank 3 b-tensors may also be selectively encoded by zeroing extraneous b-tensor components. 161D / 2D / 3D / 4D / 5D marginal distributions of the DTD are estimated from the measured b- tensors. Three 1D marginal distributions of the three diagonal components of the DTD are calculated from the rank 1 b-tensors (S404). The three 1D marginal distributions of the three diagonal components may be calculated by inverting the rank 1 b-tensors using an ℓ2-regularized non-negative least squares approach. Three 2D and three 3D marginal distributions of the DTD may be calculated from the rank- 2 b-tensors (S405 and S406). In particular, the 2D and 3D marginal distributions may be calculated by, in each case, inverting the signal acquired with rank-2 b-tensors using an ℓ2-regularized non- negative least squares approach and using previously calculated marginal distributions as constraints. More specifically, the 2D marginal distributions are calculated using the three previously calculated 1D marginal distributions of three diagonal components as constraints (S405), and then those three calculated 2D marginal distributions are used as constraints for calculated the 3D marginal distributions (S406) of the coupled components of the diffusion tensor. One 3D, three 4D, and three 5D marginal distributions of the DTD may be calculated from the rank-3 b-tensors (S407, S408, and S409). As above, in each case the 3D, 4D, and 5D marginal distributions are calculated by inverting the signal acquired with rank-3 b-tensors using an ℓ2- regularized non-negative least squares approach and using previously calculated marginal distributions as constraints. In particular, the 3D marginal distribution of the diagonal components of the diffusion tensor is calculated by inverting the signal acquired with rank-3 b-tensors using an ℓ2-regularized non-negative least squares approach and using the previously calculated 3D marginal distributions of the coupled components as constraints (S407). The three 4D marginal distributions of the diffusion tensor are calculated by inverting the signal acquired with rank-3 b- tensors using an ℓ2-regularized non-negative least squares approach and using all of the previously calculated 3D marginal distributions as constraints (S408). Similarly, the three 5D marginal Leydig 774370 HHS E-068-2023-0-US-01 distributions of the diffusion tensor are calculated by inverting the signal acquired with rank-3 b- tensors using an ℓ2-regularized non-negative least squares approach and using the previously calculated 4D marginal distributions as constraints (S409). Next, the full 6D DTD is estimated by optimizing a mathematical model of the diffusion- weighted signal captured by the MRI apparatus (S410). The diffusion-weighted signal may be based on the rank-3 b-tensors captured in the triple-PFG experiments. The model may be a mathematical model comprising the target DTD estimation. The optimization includes minimizing a difference between the model and the captured diffusion weighted signal. Further, the optimization is subject to several constraints. First, a positivity constraint is imposed on the probability density to ensure the probability density is physical. Second, a positive definiteness constraint is applied. The positive definiteness constraint zeros out a subset of the discrete diffusion tensor components that do not lie on the manifold of symmetric positive definite matrices. Third, a marginal distribution constraint is applied. The marginal distribution constraint partitions the manifold of symmetric positive definite matrices into select regions where values of previously estimated marginal distributions of a joint distribution of the DTD are above a predefined threshold. The previously estimated marginal distributions include the 16 marginal distributions estimated in the previous operations. The optimization operation may be executed by performing a convex joint optimization algorithm based on equations 3 and 4 defined above. The optimization operation outputs the estimated DTD (S411). The DTD may then be visualized, e.g., on a user interface (S412). The visualization may be done by any method, including that described below. As a person of ordinary skill in the art would apprehend, the above ordering of the operations of the method does not indicate a temporal order of each operation. Operations may be performed in other orders or concurrently without departing from the scope of the present disclosure. Further, not all operations are required for every method implemented according to aspects of the present disclosure. For example, MR signal data may have been previously obtained, and so the method may begin, for example, with the DTD estimation operations. Exemplary Visualization The 6D DTD can be visualized using various invariant measures and glyphs. For example, the size and shape distribution of the intravoxel diffusion tensors can be described by the isosurface contours of its three orthogonal invariants, namely the norm, fractional anisotropy or anisotropy and mode of anisotropy. See, e.g., Ennis, D and Kindlmann, G., “Orthogonal tensor invariants and the analysis of diffusion tensor magnetic resonance images, Magnetic Resonance in Medicine 2006, Leydig 774370 HHS E-068-2023-0-US-01 55(1):136–146 (the entire contents of which are hereby incorporated by reference herein). Various microscopic measures, such as microscopic fractional anisotropy (μFA) and microscopic orientation distribution function (μODF), may also be obtained by computing the ensemble average of these quantities directly from the non-parametric DTD. See Magdoom et al. (cited above). The higher-order sample moments, such as the second-order mean tensor, fourth-order covariance tensor, etc., are also calculated numerically. These measures are computed for each of the mesoscopic compartments found within a voxel to accurately quantify the heterogeneity present within each compartment. The DTD measured in each voxel is partitioned into various symmetry classes whose volume fraction are measured and mapped. The symmetry classes are defined based on the sorted eigenvalues of individual microscopic diffusion tensors as follows, 1) stick (λ1 > 0, λ2 ≈ λ3 ≈ 0), 2) planar (λ1 ≈ λ2, λ3 = 0), 3) cylindrical (λ1 > λ2 ≈ λ3 > 0 or λ1 ≈ λ2 > λ3 > 0), 4) spherical (λ1 ≈ λ2 ≈ λ3), and 5) general anisotropic (λ1> λ2> λ3). The DTD is also partitioned using non-parametric clustering approaches such as k-means on the 6D-DTD space as well as the orientation invariant 3D space of norm, FA, and anisotropy mode of the diffusion tensor. The above-mentioned scalar quantities and glyphs can be computed for each distinct probability mode or distinct water pool within a voxel or in the whole brain to more accurately quantify mesoscale heterogeneity. Exemplary MRI Measurements and image pre-processing Exemplary MRI data were acquired using a 64-channel RF coil on a 3T scanner (Prisma, Siemens Healthineers) capable of gradient strengths up to 80 mT / m per channel and a 200 T / m / s slew rate with an echo planar imaging (EPI) readout. The measurement was performed on healthy volunteers (N = 3) who provided informed consent in accordance with a research protocol approved by the Institutional Review Board (IRB) of the Intramural Research Program of the National Institute of Neurological Disorders and Stroke (NINDS). The estimation pipeline is validated using a macroscopically and microscopically isotropic 40% polyvinylpyrrolidone (PVP) liquid phantom. DWIs were acquired with a field of view (FOV) = 210 mm x 210 mm x 120 mm, 2 mm isotropic spatial resolution at 1965 Hz / pixel bandwidth, and GRAPPA acceleration factor = 4. The b-tensors are realized using interfused-PFG (iPFG) MR pulse sequence that can generate b-tensors of all ranks within a single spin-echo. The thermal noise in the images was filtered using a Marchenko-Pastur principal component analysis (PCA) algorithm implemented in DIPY software. See, e.g., Veraart, Jelle et al., “Denoising of diffusion MRI using random matrix theory,” NeuroImage 2016, 142:394–406 (discussing the Marchenko-Pastur PCA); and Garyfallidis, Eleftherios et al., “Dipy, a library for the analysis of Leydig 774370 HHS E-068-2023-0-US-01 diffusion MRI data,” Frontiers in Neuroinformatics 2014, 8(FEB):8 (discussing the DIPY software) (the entire contents of each of which are hereby incorporated by reference herein. The effect of eddy-current induced image translation, dilation, and shear were reduced by registering the individual DWI volumes with T2weighted structural image using a 3D affine transform implemented in FSL software prior to DTD estimation. See, e.g., Jenkinson Mark and Smith Stephen, “A global optimisation method for robust affine registration of brain images,” Medical Image Analysis 2001, 5(2):143–156; and Jenkinson, Mark et al, “Improved Optimization for the Robust and Accurate Linear Registration and Motion Correction of Brain Images,” NeuroImage 2002, 17(2):825–841 (discussing the 3D affine transform); see also Smith, Stephen M et al., “Advances in functional and structural MR image analysis and implementation as FSL,” NeuroImage 2004;23(SUPPL.1):S208-S219 (discussing FSL) (the entire contents of each of which are hereby incorporated by reference herein). Exemplary Validation of Methodology The efficacy of the reconstruction approach has been investigated using a set of synthetic DTD phantoms. First, a macroscopically and microscopically isotropic phantom akin to the point spread function (PSF) for linear systems is simulated for signal-to-noise ratios (SNRs) ranging from 50– 1000. This was followed by an investigation of the ability of the disclosed methodology to reliably measure an angle of crossing fibers using a two fiber motif whose crossing angle is varied from 45º–90 º. Finally, the methodology was evaluated using a motif consisting of equal proportions of gray matter, white matter, and CSF. Gray matter is represented by an emulsion with spherical diffusion tensors whose MD is Gamma distributed with a mean equal to that of brain tissue (i.e., 0.6 ^m2 / ms). White matter is represented as two fiber populations whose ellipsoid orientations cross at 90owith MD = 0.6 ^m2 / ms and D║= 7.5 D┴such as in coherent white matter fibers (e.g., corpus callosum). The CSF is represented using an emulsion of spherical diffusion tensors whose mean diffusivity is distributed according to a truncated normal distribution (i.e., only positive samples) with a mean close to that of free water at body temperature (i.e., 3 ^m2 / ms). The DW signal from each of these synthetic phantoms was generated for the various b-tensors tuilized with the disclosed methodology using Equation 2. Gaussian noise was added to the real and imaginary channels to study the effect of SNR on the reconstruction; SNR was set to 200 (low) and 1000 (high) for the non-diffusion weighted image. The effect of SNR on reconstructed DTD was investigated using a macroscopically and microscopically isotropic synthetic digital phantom (with MD = 0.6 ^m2 / ms) akin to a point spread Leydig 774370 HHS E-068-2023-0-US-01 function (PSF) for linear systems as shown in FIG.7A. FIG.7A includes the ground truth results along with those obtained for SNR ranging from 50–1000. The simulated motif is showng using micro-diffusion tensor ellipoids in the upper-right corder of FIG.7A. The reconstructed 6D-DTD for each SNR is visualized using the size-shape distribution depicted using iso-surfaces or iso- contours in the 3D parameter space of ellipsoid size, anisotropy, mode of anisotropy of the difuusion tensor, and the micro-orientation of the distribution function (^ODF) computed from the DTD. The first three moments of the mean diffusivity distribution (which are DTD derived scalar measures that characterize size) along with the anisotropic heterogeneity (i.e., FA and ^FA) are plotted as a function of SNR. Ground truth values (i.e., SNR = ∞) are also potted. The ground truth DTD has a spherical ODF and has a single peak in size-shape distribution contours centered at zero anisotropy and anisotropy mode with size = 1 ^m2 / ms (n.b. size is the norm and not the trace of the diffusion tensor as in MD). At low SNR, there was a large spread in the anisotropy and anisotropy mode which gradually reduced towards the ground truth with increasing SNR. The ^ODF resembled a cube at low SNR and approached towards the sphere with increasing SNR. The average mean diffusivity was almost constant and close to the ground truth for all the SNRs tested with errorrs in standard deviation and skewness decreasing with increasing SNR. The ^FA was close to 0.7 at the lowest SNR and continuously dropped towards zero with increasing SNR while the FA remained close to zero. In the investigation, it was shown that the reconstruction approaches ground truth at around SNR = 750. The ability of the reconstruction approach to capture DTD motifs expected in brain tissue is investigated in FIG.7B with SNR = 200 and 1000. FIG.7B shows a set of synthetic phantoms depticting motifs of micro-diffusion tensor ellipsoids and the reconstruction as a function of the SNR. The 6D-DTD is visualized using the size-shape distribution depicted using iso-surfaces or iso-contours in the 3D parameter space of ellipsoid size, anisotropy, mode of anisotropy of the diffusion tensors, and the micro-orientation distribution function (μODF) computed from the DTD. In FIG.7B, the first three motifs simulate two white matter fibers whose angle of crossing varies from 45o–75o. The last motif consists of a mixture of gray matter, white matter, and CSF in equal proportions within a voxel. At SNR of 1000, the angle of crossing was accurately reconstructed for all the cases with size-shape distribution approximating the ground truth. The lowering of SNR led to broadening of the distribution and subsequent blurring of the two fiber populations. For example, the 45º crossing angle was no longer captured and the anisotropy mode indicated the presence of orthotropic micro diffusion tensors as opposed to only linear anisotropic tensors. The most complex microstructural motif simulated in this study is the three compartment Leydig 774370 HHS E-068-2023-0-US-01 voxel consisting of gray, white matter and CSF in a single voxel. It has three peaks in the size- shape distribution as shown in the figure with the ^ODF showing the dominant crossing fiber orientations. The high SNR inversion was able to distinguish all the three compartments in the size- shape distribution plot. The gray and white matter were inseparable in low SNR inversion but CSF could be still be distinguished from tissue. The ^ODF was however accurately captured for both the SNR cases albeit with some broadening at lower SNR. The results demonstrate the orientation discrimination ability of the disclosed methodology, for example the favorable ability to resolve fibers crossing at an acute angle of 60 º even at an SNR = 200. Results are next reported from experimental data acquired with the PVP phantom for a discrete number of voxels which are shown in FIG.8. A center voxel and four voxels at the edges of the phantom where artifacts such as from gradient eddy currents are worse were examined. The results from these voxels were largely identical. The size-shape distribution showed a single peak centered at the expected 1 ^m2 / ms size and spread in the anisotropy mode. The ^ODF has a cuboid appearance when computed from the DWI data with a nominal SNR of 200. As the image is rebinned, making the voxels more coarse, that the ^ODF more closely approximates a sphere. The scalar measures and glyphs derived from the reconstructed DTD for a representative slice of PVP is shown in Figure 9. These include the non-diffusion weighted image, fractional anisotropy and microscopic fractional anisotropy, and moments of the mean diffusivity distribution. The results of Riemannian clustering of DTDs in all voxels within the slice are also included along with the map of its proportions in each voxel (i.e., volume fraction). The FA was very close to zero but the micro-FA was approximately around 0.5 throughout the slice at the voxel resolution of the experiment. The ^FA approaches zero as the SNR is increased by rebinning the image at successively lower resolutions. The MD was approximately 0.6^m2 / ms with very low standard deviation and skewness. Most of the voxels in the slice were clustered into one bin with anisotropy mode close to zero and very small FA with size approximately equal to 1 ^m2 / ms. The second cluster had the anisotropy mode shifted downward indicating planar diffusion tensors. Discussion of Exemplary Study As explained in detail above, the present disclosure provides an MRI-based methodology that can efficiently measure and map mesoscopic water pools that reside in live brain tissue. These water pools are statistically described using a non-parametric probability density function (PDF) of diffusion tensors (i.e., DTD), which is measured in each imaging voxel. Given that the DTD is sparse, the estimation of this six-dimensional probability density is performed using a hierarchy of Leydig 774370 HHS E-068-2023-0-US-01 marginal densities, which progressively shrink the space of admissible solutions. This framework has been thoroughly vetted by the inventors using a set of realistic, synthetic digital phantoms, and using a polymer liquid phantom, whose ground truth is known. The disclosed approach has also been vetted on live human brain tissue, where it holds the promise of revealing mesoscopic brain architecture and potentially detecting subtle changes in the brain which are often invisible in traditional MRI scans given the higher dimensionality of DTD. A key premise underlying the validity of the approach of the present disclosure is that over the range of experimental DWI parameters used (i.e., pulse gradient widths, diffusion times, diffusion gradient orientations and magnitudes), that the MR signal can be modeled as arising from a superposition of contributions from non-exchanging water pools each describable using an anisotropic Gaussian net displacement distribution so that Equation (1) holds. Specifically, the use of the anisotropic Gaussian diffusion kernel e−B:Din the integrand of Equation (1) is assumed to be valid. While several DTD MRI methods have been proposed to date, the underlying assumption developed by the present inventors has previously not been examined, i.e., stated, tested, and validated. At least over a clinical range of DWI parameters routinely used, and in particular, those used in the exemplary study, the present inventors have found that this assumption does hold, justifying the use of Equation (1), and the pipeline disclosed herein. Briefly, the combination of tissue homeostasis and tissue water permeability facilitated by the plasma membranes, aquaporins, etc., makes it unlikely for cell and tissue water to be purely restricted, which would violate the use of the Gaussian kernel for DTD estimation. The large surface-to-volume ratio (S / V) of the mesoscopic voxels relative to the diffusion time used prevents water from exchanging among them at these time scales. It is important to note that noise affects different DTD-derived parameters differently from their DTI counterparts. The macroscopically and microscopically isotropic PVP solution provides an ideal MRI phantom with which to assess the effects of SNR on these quantities. For example while it can be seen that with the FA ≈ 0 within the PVP spherical phantom, the μFA ≈ 0.5. In principle, the μFA should also be zero but noise creeps into the estimate of the microtensors and becomes amplified when calculating their respective FAs, so the ensemble-averaged FA is artificially over-estimated. This problem can be remedied by increasing the SNR of the DWIs, in this case by increasing voxel size through rebinning the image. A similar observation can be made for the μODF. The DTI ODF for the PVP phantom should be spherical reflecting macroscopic isotropy of the medium, as well as the μODF, which reflects its microscopic isotropy. However, it can be seen that the μODF is cuboid rather than spherical at low SNR. As the images are successively rebinned with lower and lower voxel resolution, the μODF becomes more and more Leydig 774370 HHS E-068-2023-0-US-01 spherical. These findings suggest that one must exercise caution when interpreting μFA and μODF values in DTD MRI data even when DTI-derived quantities, like MD and FA are accurate. Interestingly, to date, there have been no published studies assessing the accuracy and precision of DTD-derived parameters despite widespread attempts to disseminate DTD-MRI clinically. This is important because our the present inventors have revealed that DTD-derived parameters are affected by background noise to a greater extent than DTI-derived parameters and more than previously understood. The advantage of multidimensional NMR is the specificity it provides by separating features in nD space which are often indistinguishable in 1D NMR. This however comes at a price of non-linear increase in the number of unknowns to be estimated, relative to the finite amount of data that can practically be acquired in a clinical setting, i.e., the so-called “curse of dimensionality”. Given physical quantities are often sparse in the nD space, the framework of the present disclosure provides enormous computational economy and improvement in accuracy by resolving peaks at very high spectral resolution within regions of interest (i.e., adaptive gridding). As additional constraints from a hierarchy of marginal distributions are applied, both the number of degrees of freedom and the volume of the solution space are markedly reduced so that admissible solutions become confined to ever smaller domains as shown in FIGS.1 and 2 for the 3D case thus overcoming the curse of dimensionality. For instance, the number of unknowns needed to solve the most complex 4- compartment motif in the example simulation was approximately 15 parts-per- million (ppm). It should be noted that implementations of the present disclosure are not limited by the form of Equation (1). The number of variables can be expanded by augmenting the kernel, ^^^^^^in Equation 2 to include target distributions ^^(^^) (additional or alternate to the DTD). For example, the implementation can implemented to further characterize the relaxation processes, such as those ^^^^^^^^^^1,^^^^^^2,^^, where TR, TE and TI are the repetition, echo and inversion times. Then, the number of marginal distributions (i.e., ^^(^^1), ^^(^^2), ^^(^^),^^(^^, ^^1), and ^^(^^, ^^2))needed to compute increases with the number of dependent variables. Theincrease in the number of constraints, however, is not an obstacle. Although more computations are required to simultaneously satisfy them, these added constraints further restrict the domain of admissible solutions, again helping to overcome the “curse of dimensionality”. It should also be noted that performing the ILT is fundamentally an ill-conditioned problem and requires a priori assumptions to reduce the effects of noise in the data. In the exemplary study, Leydig 774370 HHS E-068-2023-0-US-01 it has been assumed that the solution is smooth, as indicated by the ℓ2 norm-regularization term in Equation (3), which numerically broadens the peaks in the 6D space, akin to the line broadening effect in NMR spectroscopy. This affects the resolving power of the method as observed in the inversions of simulated data with low SNR, and introduces SNR dependent errors in the inversion, such as the cube like µODF observed in the PVP phantom data, which was caused by the slight broadening in the distribution of off-diagonal diffusion tensor elements. This effect however depends on the SNR and data size (i.e., experimental design), which can be alleviated with the new generation of high-performance MRI scanners that allow smaller echo times (TE) and faster scan repetition rate (TR). As should be recognized by those skilled in the art, the experimental design of the present disclosure may be further optimized to improve accuracy of estimates of p(D), to minimize the number of DWI acquisitions and total scan time, and to increase accuracy (i.e., reduce statistical bias in estimated quantities), precision, and robustness. For example, the experimental design may be systematically pruned to minimize the number of DWIs required, while still achieving precise and accurate DTD reconstructions to enable radiological migration. This pruning may be implemented by defining the optimal number of rank-1, rank-2, and rank-3 b-tensors needed to achieve these ends. For example, rank-3 b-tensors could be excluded while still maintaining sufficient precision and accuracy for the estimated DTD, given the demands of the implementation. Additionally, in an implementation, the DTD may be restricted to be cylindrically symmetric, a priori, as is done in Topgaard et al., resulting in 4D rather than 6D DTDs. Another possible changed implementation could be to relax some of the constraints of the disclosed framework, using only single-PFG DWI MRI data, i.e., employing only rank-1 b-tensors, which would only allow for estimation of three individual marginal distributions, p(Dxx), p(Dyy), and p(Dzz), but no other individual or joint marginal distributions. Given adequate SNR, the benefit of estimating p(D) empirically, particularly in voxels containing multiple tissue types, (e.g., white and gray matter, CSF, and / or crossing white matter pathways) having distinct water pools, is that it is expected to find multimodal DTDs whose distinct peaks are centered about different mean diffusion tensors (i.e., points) within the 6D tensor space whose statistical features we could probe independently. One way to glean and summarize sub-voxel microstructural information in this case is to fit parametric distributions to these distinct probability clouds in 6D tensor space. This decomposition into parametric distributions in a principled way to summarize features of intravoxel microstructural motifs, would also help eliminate contributions from extraneous or artifactual components, such as arising from pools of Leydig 774370 HHS E-068-2023-0-US-01 ”free water”, as in so that parenchymal water diffusion can be separated from dispersion or pseudo- diffusion effects associated with free CSF compartments, such as in the ventricles. Exemplary Data Acquisition and Control Apparatus FIG.10 and the following discussion are intended to provide a brief, general description of an exemplary computing / data acquisition environment in which the disclosed technology may be implemented. Although not required, the disclosed technology is described in the general context of computer executable instructions, such as program modules, being executed by a personal computer (PC), a mobile computing device, tablet computer, or other computational and / or control device. Generally, program modules include routines, programs, objects, components, data structures, etc., that perform particular tasks or implement particular abstract data types. Moreover, the disclosed technology may be implemented with other computer system configurations, including, multiprocessor systems, network PCs, minicomputers, mainframe computers, and the like. The disclosed technology may also be practiced in distributed computing environments where tasks are performed by remote processing devices that are linked through a communications network. In a distributed computing environment, program modules may be located in both local and remote memory storage devices. An exemplary system for implementing the disclosed technology may include a general- purpose computing device in the form of an exemplary conventional PC 1100, including one or more processing units 1102, a system memory 1104, and a system bus 1106 that couples various system components including the system memory 1104 to the one or more processing units 1102. The system bus 1106 may be any of several types of bus structures including a memory bus or memory controller, a peripheral bus, and a local bus using any of a variety of bus architectures. The exemplary system memory 1104 includes read only memory (ROM) 1108 and random-access memory (RAM) 1110. A basic input / output system (BIOS) 1112, containing the basic routines that help with the transfer of information between elements within the PC 1100, is stored in ROM 1108. The exemplary PC 1100 further includes one or more storage devices 11110 such as a hard- disk drive for reading from and writing to a hard disk, a magnetic disk drive for reading from or writing to a removable magnetic disk, an optical disk drive for reading from or writing to a removable optical disk (such as a CD-ROM or other optical media), and a solid-state drive. Such storage devices can be connected to the system bus 1106 by a hard-disk drive interface, a magnetic disk drive interface, an optical drive interface, or a solid-state drive interface, respectively. The drives and their associated computer readable media provide nonvolatile storage of computer- readable instructions, data structures, program modules, and other data for the PC 1100. Other Leydig 774370 HHS E-068-2023-0-US-01 types of computer-readable media which can store data that is accessible by a PC, such as magnetic cassettes, flash memory cards, digital video disks, CDs, DVDs, RAMs, ROMs, and the like, may also be used in the exemplary operating environment. The PC 1100 may also utilize remote, networked (e.g. Cloud) storage. A number of program modules may be stored in the storage devices 11110 including an operating system, one or more application programs, other program modules, and program data. A user may enter commands and information into the PC 1100 through one or more input devices 1140 such as a keyboard and a pointing device such as a mouse. Other input devices may include a digital camera, microphone, joystick, game pad, satellite dish, scanner, or the like. These and other input devices are often connected to the one or more processing units 1102 through a serial port interface that is coupled to the system bus 1106, but may be connected by other interfaces such as a parallel port, game port, or universal serial bus (USB). A monitor 1146 or other type of display device is also connected to the system bus 1106 via an interface, such as a video adapter. Other peripheral output devices, such as speakers and printers may be included. The PC 1100 may operate in a networked environment using logical connections to one or more remote computers, such as a remote computer 1160. In some examples, one or more network or communication connections 1150 are included. The remote computer 1160 may be another PC, a server, a router, a network PC, or a peer device or other common network node, and typically includes many or all of the elements described above relative to the PC 1100, although only a memory storage device 1161 has been illustrated in FIG.9. The personal computer 1100 and / or the remote computer 1160 can be connected to a logical a local area network (LAN) and a wide area network (WAN). Such networking environments are commonplace in offices, enterprise-wide computer networks, intranets, and the Internet. When used in a LAN networking environment, the PC 1100 is connected to the LAN through a network interface. When used in a WAN networking environment, the PC 1100 typically includes a modem or other means for establishing communications over the WAN, such as the Internet. In a networked environment, program modules depicted relative to the personal computer 1100, or portions thereof, may be stored in the remote memory storage device or other locations on the LAN or WAN. The network connections shown are exemplary, and other means of establishing a communications link between the computers may be used (for example through the Cloud). The memory 1104 generally includes computer-executable instructions for selecting gradient fields, averaging acquired signals, calculation of b-matrices and mean diffusion tensor and covariance. For example, memory portion 1162 can store computer-executable instructions for estimating the DTD, and computer-executable instructions for b-matrix computation and random b- Leydig 774370 HHS E-068-2023-0-US-01 magnitude and direction selection can be stored at 1173B or previously computed b-matrix specification can be stored in memory portion 1173A. Computer-executable instructions for processing acquired signals (for example, determining mean diffusion tensor and covariance) can be stored in a memory portions 1161, 1171. Computer-executable instructions for data acquisition and control are stored in a memory portion 1170. Acquired and processed data (e.g., images based on mean diffusion tensor images) can be displayed using computer-executable instructions stored at memory portion 1171. As noted above, data acquisition, processing, and instrument control can be provided at an MRI system 1174, or distributed at one or more processing devices using a LAN or WAN. While aspects of the present disclosure have been illustrated and described in detail in the drawings and foregoing description, such illustration and description are to be considered illustrative or exemplary and not restrictive. It will be understood that changes and modifications may be made by those of ordinary skill within the scope of the following claims. In particular, the present invention covers further embodiments with any combination of features from different embodiments described above and below. Additionally, statements made herein characterizing the invention refer to an embodiment of the invention and not necessarily all embodiments. The terms used in the claims should be construed to have the broadest reasonable interpretation consistent with the foregoing description. For example, the use of the article “a” or “the” in introducing an element should not be interpreted as being exclusive of a plurality of elements. Likewise, the recitation of “or” should be interpreted as being inclusive, such that the recitation of “A or B” is not exclusive of “A and B,” unless it is clear from the context or the foregoing description that only one of A and B is intended. Further, the recitation of “at least one of A, B and C” should be interpreted as one or more of a group of elements consisting of A, B and C, and should not be interpreted as requiring at least one of each of the listed elements A, B and C, regardless of whether A, B and C are related as categories or otherwise. Moreover, the recitation of “A, B and / or C” or “at least one of A, B or C” should be interpreted as including any singular entity from the listed elements, e.g., A, any subset from the listed elements, e.g., A and B, or the entire list of elements A, B and C.

Claims

Leydig 774370 HHS E-068-2023-0-US-01 We claim:

1. A method for estimating a diffusion tensor distribution (DTD), which is a function of diffusion tensor components that characterizes diffusion properties of a subject, the method comprising: performing an nD inverse Laplace transform of equation 1, subject to a plurality of constraints, to determine the estimated DTD, where equation 1 is given as: ^^(^^) = ^^0 ∫ ^^−^^:^^^^(^^) ^^^^, wherein ^^0is an MR signal without diffusion weighting, ^^ is a second-order, symmetric b-tensor and ^^ is the a second-order, symmetric diffusion tensor, ^^(^^)is the estimated DTD, and : is a tensor dot product operation, wherein the plurality of constraints comprises mathematical, physical, and / or statistical constraints.

2. The method of claim 1, wherein the plurality of constraints comprises: a symmetric positive semi definiteness constraint, which zeros out a subset of the discrete diffusion tensor components that do not lie on a manifold of symmetric positive semidefinite matrices; a marginal distribution constraint, which partitions the manifold of symmetric positive definite matrices into select regions where values of previously estimated marginal distributions of the DTD are above a predefined threshold; and a physical constraint on possible values of components of the diffusion tensor.

3. The method of claim 2, wherein the performing the nD inverse Laplace transform is further subject to an uncertainty constraint, which requires that difference between marginal distributions of components of the estimated DTD and the previously estimated marginal distributions to be below a predetermined uncertainty threshold.

4. The method of claim 2, wherein performing the nD inverse Laplace transform comprises solving a convex joint optimization problem given by:Leydig 774370 HHS E-068-2023-0-US-01 wherein, ^^(^^^^) is the captured diffusion-weighted signal, ^^^^represents a discrete set of b-tensors, ^^0 is the MR signal without diffusion weighting, ^^^^^^ = exp〖(−^^^^ ∶ ^^^^〗) and is an^^ × ^^ matrix kernel, ^^^^ represents the discrete set of diffusion tensors, and ^^(^^^^) is an ^^ × 1vector consisting of the discretized components of the DTD, ^^^^^^ is the ^^ × ^^ matrix operator,which maps the DTD to a marginal density of interest, κ is the tolerance parameter, and ^^′is a previously estimated marginal distribution, ℳ+is the manifold of symmetric positive definitematrices, min‖^^ − ^^‖22 is an optimization operator minimizing the ℓ2-norm, and ‖^^−^^‖2is an operator for calculating a sum of squares error, and wherein ^^0^^^^^^^^(^^^^)is a model predicting the captured diffusion-weighted signal ^^(^^^^).

5. The method of claim 4, wherein the convex joint optimization problem is performed using CVXPY with a MOSEK solver or other suitable convex optimization solvers.

6. The method of claim 2, the method further comprising: calculating, from a plurality of b-tensors, sixteen marginal distributions as the previously estimated marginal distributions, the plurality of b-tensors being obtained from MR measurements on the subject tissue.

7. The method of claim 6, wherein the sixteen marginal distributions comprise three 1D marginal distributions of three diagonal components of a diffusion tensor and three 2D marginal distributions of pairs of diagonal components, four 3D marginal distributions, three 4D marginal distributions, and three 5D marginal distributions involving combinations of components of the diffusion tensor.

8. The method of claim 7, wherein the plurality of b-tensors comprises a plurality of rank 1 b-tensors obtained from single-PFG experiments conducted as the MR measurements on the subject tissue, and wherein calculating the sixteen marginal distributions from the b-tensors comprises calculating the three 1D marginal distributions of the three diagonal components by inverting the rank 1 b-tensors using an ℓ2-regularized non-negative least squares approach.

9. The method of claim 8, the method further comprising executing, by the MR apparatus on the subject tissue, the single-PFG experiments to obtain the rank 1 b-tensors, whereinLeydig 774370 HHS E-068-2023-0-US-01 during the single-PFG experiments, gradients are incremented individually along three axis of the MRI apparatus such that resulting b-values of along an axis are linearly spaced over an interval.

10. The method of claim 7, wherein the plurality of b-tensors comprises a plurality of rank 2 b-tensors obtained from a double-PFG experiments conducted as the MR measurements on the subject tissue, wherein calculating the sixteen marginal distributions from the b-tensors comprises: calculating, three 2D marginal distributions by inverting the rank-2 b-tensors using an ℓ2-regularized non-negative least squares approach using the three 1D marginal distributions of three diagonal components as a constraint; and calculating three 3D marginal distributions by inverting the rank 2 b-tensors using an ℓ2-regularized non-negative least squares approach using the three 1D marginal distributions of three diagonal components and 2D marginal distributions as a constraint.

11. The method of claim 10, the method further comprising executing, by the MR apparatus on the subject tissue, the double-PFG experiments to obtain the rank 2 b-tensors, wherein the rank 2 b-tensors are selectively encoded by zeroing extraneous b-tensor components.

12. The method of claim 1, the method further comprising capturing with the MR apparatus the diffusion-weighted signal by conducting triple-PFG experiments to capture a set of rank 3 b-tensors chosen randomly on the manifold of symmetric positive definite matrices and selectively encoding diffusion tensor elements by zeroing extraneous b-tensor components.

13. The method of claim 1, the method further comprises visualizing the DTD in a plot using glyphs and by: determining a size and shape distribution of intravoxel diffusion tensors based on isosurface contours of a norm, fractional anisotropy or anisotropy, and mode of anisotropy of the DTD; determining an orientation distribution function of the intravoxel diffusion tensors based on the micro-ODF; computing for each of a plurality of macroscopic compartments found within each voxel: a plurality of microscopic measures, which comprise at least one of a microscopic fractional anisotropy or a microscopic orientation distribution function, by computing an ensemble average of a respective microscopic measure directly from the DTD; andLeydig 774370 HHS E-068-2023-0-US-01 higher-order sampling moments, which comprise at least one of a second-order mean tensor or a fourth-order covariance tensor, numerically from the DTD; and partitioning the DTD measured in each of the voxels into various symmetry classes, the symmetry classes having a volume faction within each voxel that is measured and mapped, and the symmetry classes being defined based on sorted eigenvalues of individual microscopic diffusion tensors as follows: stick for (λ1> 0, λ2≈ λ3≈ 0), planar for (λ1 ≈ λ2, λ3 = 0); cylindrical for (λ1 > λ2 ≈ λ3 > 0 or λ1 ≈ λ2 > λ3 > 0); spherical for (λ1≈ λ2≈ λ3), and general anisotropic for (λ1 > λ2 > λ3); and partitioning the DTD measured in each of the voxels and the whole sample using non- parametric k-means clustering approach.

14. A non-transitory computer readable medium, which when executed by one or more processors, cause execution of a method for estimating a diffusion tensor distribution (DTD), which comprises a plurality of discrete diffusion tensor components and characterizes diffusion properties of a subject tissue, the method comprising: performing an nD inverse Laplace transform of equation 1, subject to a plurality of constraints, to determine the estimated DTD, where equation 1 is given as: ^^(^^) = ^^0 ∫ ^^−^^:^^^^(^^) ^^^^, wherein ^^0is an MR signal without diffusion weighting, ^^ is a second-order, symmetric b-tensor and ^^ is the diffusion tensor, ^^(^^) is the estimated DTD, and : is a tensor dot product operation, wherein the plurality of constraints comprises mathematical, physical, and / or statistical constraints.

15. A magnetic resonance (MR) apparatus, the MR apparatus comprising a processing system, comprising one or more processors, configured to execute a MR method for estimating a diffusion tensor distribution (DTD), which comprises a plurality of discrete diffusion tensor components and characterizes diffusion properties of a subject tissue, the MR method comprising: performing an nD inverse Laplace transform of equation 1, subject to a plurality of constraints, to determine the estimated DTD, where equation 1 is given as: ^^(^^) = ^^0 ∫ ^^−^^:^^^^(^^)^^^^,Leydig 774370 HHS E-068-2023-0-US-01 wherein ^^0is an MR signal without diffusion weighting, ^^ is a second-order, symmetric b-tensor and ^^ is the diffusion tensor, ^^(^^) is the estimated DTD, and : is a tensor dot product function, wherein the plurality of constraints comprises comprises mathematical, physical, and / or statistical constraints.

16. The MR apparatus of claim 15, the MR apparatus further comprising: a radiofrequency (RF) generator coupled to an RF coil configured to controllably apply RF pulses to the subject tissue; a gradient controller coupled to a plurality of gradient coils configured to apply gradient pulses to the subject tissue; an axial magnet controller coupled to an axial magnet configured to apply a magnetic field on the subject tissue; an RF receiver coupled to a receiver coil configured to detect or measure net magnetization of spins of water molecules of the subject tissue; and a controller in communication with the processing system, the RF generator, the gradient controller, the axial magnet controller, and the RF receiver, wherein the controller is configured to collectively control the RF generator, the gradient controller, the axial magnet controller, and the RF receiver to obtain a data set by conducting single-PFG experiments; double-PFG experiments; and triple-PFG experiments.

17. A method for estimating a variable distribution characterizing properties of a subject tissue from a plurality of magnetic resonance (MR) signals, the method comprising: performing an nD inverse Laplace transform of the following equation, subject to constraints, to determine the estimated variable distribution, where the equation is given as: ^^(^^) = ^^0 ∫ ^^−^^.^^^^(^^) ^^^^, wherein ^^0is an MR signal without contrast weighting, ^^ is a weighting vector and ^^ is a variable vector, ^^(^^)is the estimated variable distribution, and . is a dot product function, wherein the constraints comprise a marginal distribution constraint, constraining the estimated variable distribution to comprise values of previously estimated marginal distributions of the variable distribution that are above a predefined threshold, and a positivity constraint which ensures ^^(^^) is positive.

18. The MR method of claim 17, wherein the variable distribution comprises at least one of a diffusion tensor distribution (DTD) or a relaxation decay distribution, andLeydig 774370 HHS E-068-2023-0-US-01 wherein previously estimated marginal distributions comprise at least one of 1D, 2D, 3D, 4D, or 5D marginal DTDs, a T1 relaxation time distribution, a T2 relaxation time distribution or a T2*relaxation time distribution.

Images