Recovery of missing information in diffusion magnetic resonance imaging data

Abstract

There is described herein a method for recovering missing information in diffusion magnetic resonance imaging (dMRI) data. The data are modeled according to the theory of moving frames and regions where frame information is missing are reconstructed by performing diffusions into the regions. Local orthogonal frames computed along the boundary of the regions are rotated into the regions. Connection parameters are estimated at each new data point obtained by a preceding rotation, for application to a subsequent rotation.

Claims

1. A method for recovering missing information in diffusion magnetic resonance imaging (dMRI) data, the method comprising: fitting, at each data point of the dMRI data, a local orthogonal frame expressed as (i, j, k) to represent three directions; determining connection parameters c.sub.ijk at each data point to represent a rotation of the local orthogonal frame in a spatial neighborhood of the data point in accordance with the dMRI data; and performing inward diffusions into a region where frame information is missing by: rotating the local orthogonal frame in a heading direction at data points starting along a boundary of the region and moving into the region; and estimating connection parameters c.sub.ijk at each new data point obtained by a preceding rotation of the local orthogonal frame, for application to a subsequent rotation of the local orthogonal frame at the new data point.

2. The method of claim 1, wherein rotating the local orthogonal frame comprises computing new values for the three directions (i, j, k) independently.

3. The method of claim 2, further comprising adjusting the computed new values to enforce orthogonality of the three directions (i, j, k).

4. The method of claim 1, wherein determining connection parameters c.sub.ijk at each data point of the dMRI data comprises computing closed-form connections in linear space.

5. The method of claim 1, wherein estimating connection parameters c.sub.ijk at each new data point comprises selecting a computation scheme from a group comprising finite differentiation, energy minimization, and closed-form connections in linear space.

6. The method of claim 5, wherein selecting the computation scheme comprises selecting as a function of local parameters of each one of the new data points.

7. The method of claim 1, wherein performing inward diffusions into a region of missing information comprises rotating the local orthogonal frame and estimating the connection parameters c.sub.ijk for a plurality of layers of the region in an iterative fashion.

8. The method of claim 7, wherein each one of the plurality of layers has a thickness corresponding to one voxel.

9. The method of claim 1, wherein fitting, at each data point of the dMRI data, a local orthogonal frame expressed as (i, j, k) to represent three directions comprises using the dMRI data to obtain a first one of the three directions and basing a second one and a third one of the three directions on a local geometry of a heart wall.

10. The method of claim 9, wherein the region is a volumetric region or a distributed region, and wherein missing information is recovered for a plurality of regions in the heart wall in order to reconstruct a three-dimensional image of fiber orientations.

11. A system for recovering missing information in diffusion magnetic resonance imaging (dMRI) data, the system comprising: a processing unit; and a non-transitory memory communicatively coupled to the processing unit and comprising computer-readable program instructions executable by the processing unit for: fitting, at each data point of the dMRI data, a local orthogonal frame expressed as (i, j, k) to represent three directions; determining connection parameters c.sub.ijk at each data point to represent a rotation of the local orthogonal frame in a spatial neighborhood of the data point in accordance with the dMRI data; and performing inward diffusions into a region where frame information is missing by: rotating the local orthogonal frame in a heading direction at data points starting along a boundary of the region and moving into the region; and estimating connection parameters c.sub.ijk at each new data point obtained by a preceding rotation of the local orthogonal frame, for application to a subsequent rotation of the local orthogonal frame at the new data point.

12. The system of claim 11, wherein rotating the local orthogonal frame comprises computing new values for the three directions (i, j, k) independently.

13. The system of claim 12, wherein the program instructions are further executable for adjusting the computed new values to enforce orthogonality of the three directions (i, j, k).

14. The system of claim 11, wherein determining connection parameters c.sub.ijk at each data point of the dMRI data comprises computing closed-form connections in linear space.

15. The system of claim 11, wherein estimating connection parameters c.sub.ijk at each new data point comprises selecting a computation scheme from a group comprising finite differentiation, energy minimization, and closed-form connections in linear space.

16. The system of claim 15, wherein selecting the computation scheme comprises selecting as a function of local parameters of each one of the new data points.

17. The system of claim 11, wherein performing inward diffusions into a region of missing information comprises rotating the local orthogonal frame and estimating the connection parameters c.sub.ijk for a plurality of layers of the region in an iterative fashion.

18. The system of claim 17, wherein each one of the plurality of layers has a thickness corresponding to one voxel.

19. The system of claim 11, wherein fitting, at each data point of the dMRI data, a local orthogonal frame expressed as (i, j, k) to represent three directions comprises using the dMRI data to obtain a first one of the three directions and basing a second one and a third one of the three directions on a local geometry of a heart wall.

20. The system of claim 19, wherein the region is a volumetric region or a distributed region, and wherein missing information is recovered for a plurality of regions in the heart wall in order to reconstruct a three-dimensional image of fiber orientations.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) Further features and advantages of the present invention will become apparent from the following detailed description, taken in combination with the appended drawings, in which:

(2) FIG. 1 is a flowchart of an exemplary computer-implemented method for recovering missing information in diffusion magnetic resonance imaging (dMRI) data;

(3) FIG. 2A illustrates the turning of frame axes at x expressed in the local basis f.sub.1, f.sub.2, f.sub.3 when x moves in the direction v;

(4) FIG. 2B illustrates frame field variation characterized by the connections c.sub.ijk for i=1;

(5) FIG. 3 illustrates a differential of the frame axis f.sub.1(x) expressed in local coordinates;

(6) FIG. 4 is a flowchart of an example embodiment of the method of inward diffusion for regions having missing frame information;

(7) FIG. 5 illustrates reconstructed volumes and errors using the proposed method for a damaged portion of a rat short axis slab and a sparse rule-based synthetic field;

(8) FIG. 6 illustrates mean angular reconstruction error .sub.1 between reconstructed fibers and the ground truth for increasingly Poisson-sparse volumes (top) using the proposed method, and reconstruction for increasingly interleaved volumes (bottom) in a long-axis cutout;

(9) FIG. 7 illustrates mean volume error .sub.1 of reconstruction for various methods;

(10) FIG. 8 illustrates an environment comprising a dMRI data reconstruction system;

(11) FIG. 9 is a block diagram of an exemplary embodiment of the dMRI data reconstruction system;

(12) FIG. 10 is a block diagram of an exemplary embodiment of a data modeling unit;

(13) FIG. 11 is a block diagram of an exemplary embodiment of a data recovery unit; and

(14) FIG. 12 is a block diagram of an exemplary hardware and software implementation of the method of FIG. 1.

(15) It will be noted that throughout the appended drawings, like features are identified by like reference numerals.

DETAILED DESCRIPTION

(16) There is illustrated in FIG. 1 a flowchart of a computer-implemented method 100 for recovering missing information in diffusion Magnetic Resonance Imaging (dMRI) data. The method 100 is based on modelling the dMRI data in accordance with the theory of moving frames. Step 102 comprises fitting, at each data point of the dMRI data, a local orthogonal frame expressed as (i, j, k) to represent three directions. Step 104 comprises determining connection parameters c.sub.ijk at each data point to represent the rotation of the local orthogonal frame in the spatial neighborhood of the data point, so as to best fit the dMRI data. Steps 102 and 104 are presented mathematically below.

(17) Let a point x=.sub.ix.sub.ie.sub.i.sup.3 be expressed in terms of e.sub.1, e.sub.2, e.sub.3, the natural basis for .sup.3. Next, we define a right-handed orthonormal frame field f.sub.1,f.sub.2,f.sub.3:.sup.3.fwdarw..sup.3. Each frame axis can be expressed by the rigid rotation f.sub.i=.sub.ja.sub.ije.sub.j where A={a.sub.ij}.sup.33 is a differentiable attitude matrix such that A.sup.1=A.sup.T. Treating f.sub.i and e.sub.i as symbols, we can write Equation (1) below. Further, since each e.sub.i is constant, the differential geometry of the frame field is completely characterized by A. Taking the exterior derivative on both sides yields Equation (2) where d denotes the exterior derivative, and C=(dA)A.sup.1={c.sub.ij}.sup.33 is the Maurer-Cartan matrix of connection forms c.sub.ij. Writing f.sub.i as symbols, then Equation (2) can be understood as df.sub.1=.sub.jc.sub.ijf.sub.j.

(18) $\begin{array}{cc}{\left[f_1,f_2,f_3\right]}^T={A\left[e_1{e}_2{e}_3\right]}^T& \left(1\right)\\ {d\left[f_1,f_2,f_3\right]}^T=\left(\mathrm{dA}\right){A^{-1}\left[f_1,f_2,f_3\right]}^T={C\left[f_1,f_2,f_3\right]}^T& \left(2\right)\\ C=\left[\begin{array}{ccc}0& c_{12}& c_{13}\\ -c_{12}& 0& c_{23}\\ -c_{13}& -c_{23}& 0\end{array}\right]& \left(3\right)\\ {\tilde{f}}_i\left(x_0+v\right)=f_i+d{f}_i.Math.v.Math.+\left({.Math.v.Math.}^2\right)f_i+\underset{ji}{.Math.}{f}_j\underset{k}{.Math.}{v}_k{c}_{\mathrm{ijk}}& \left(4\right)\end{array}$

(19) The Maurer-Cartan matrix is a skew symmetric and hence referring to Equation (3) there are at most 3 independent, non-zero one-forms: c.sub.12, c.sub.13, and c.sub.23. One-forms operate on tangent vectors through a process denoted contraction, written as dwv for a general one-form dwv=.sub.iw.sub.ide.sub.i and tangent vector v on .sup.3, which yields dwv=.sub.iw.sub.ide.sub.i.sub.jv.sub.je.sub.j=.sub.iw.sub.iv.sub.i since de.sub.ie=.sub.ij, where .sub.ij is the Kronecker delta.

(20) The space of linear models for smooth frame fields is fully parametrized by the one-forms c.sub.ij. This space can be explored by considering the motion of f.sub.i in a direction v=.sub.kv.sub.kf.sub.k using the first order terms of a Taylor series centered at x.sub.0 as given by Equation (4) where f.sub.i and df.sub.i are evaluated at x.sub.0, and c.sub.ijkc.sub.ijf.sub.k are the connection forms of the local frame. Since only 3 unique non-zero combinations of c.sub.ij are possible, there are in total 9 connections c.sub.ijk. These coefficients express the rate of turn of the frame vector f.sub.i towards f.sub.j when x moves in the direction f.sub.k. Referring to FIGS. 2A and 2B, there is illustrated the behavior of the frame field described by c.sub.ijk. For example, with f.sub.1 taken to be the local orientation of a fiber and f.sub.3 taken to be the component of the heart wall normal orthogonal to f.sub.1, then c.sub.131 measures the circumferential curvature of a fiber and c.sub.123 measures the change in its helix angle. FIG. 2A illustrates the turning of frame axes at x expressed in the local basis f.sub.1, f.sub.2, f.sub.3 when x moves in the direction v. FIG. 2B illustrates the frame field variation characterized by the connections c.sub.ijk for i=1.

(21) A first order generator for frame fields using Equation (4) requires knowledge of the underlying connection forms c.sub.ijk. Three different ways of computing the connection forms c.sub.ijk may be used. The computation schemes are referred to herein as finite differentiation, energy minimization, and closed-form connections in linear space.

(22) The finite differentiation computation scheme is a direct estimate based on finite differences. In smooth frame fields, the connection one-forms c.sub.ij can be directly obtained using Equation (2), i.e., df.sub.i.Math.f.sub.k=(.sub.j.sup.3c.sub.ijf.sub.i).Math.f.sub.k=.sub.j.sup.3c.sub.ij.sub.jk=c.sub.ik. The differentials df.sub.i can be computed by applying the exterior derivative for a function, i.e., for the k'th component f.sub.i, f.sub.ik:.sup.3.fwdarw.,

(23) $f_{\mathrm{ik}}=\underset{l}{\stackrel{3}{.Math.}}\frac{f_{\mathrm{ik}}}{x_l}d{e}_l$
as given by Equation (5) where

(24) $J_i=\left[\frac{f_{\mathrm{ip}}}{x_q}\right]{33}^{}$
is the Jacobian matrix of partial derivatives of f.sub.i. Now setting v=f.sub.k we obtain Equation (6). The Jacobian matrix J, can be approximated to first order using, e.g., finite differences on f.sub.i with a spacing of size x:

(25) $\frac{f_{\mathrm{ij}}}{x_k}\left(x\right)\frac{f_{\mathrm{ij}}\left(x+e_k\right)-f_{\mathrm{ij}}\left(x\right)}{x}.$

(26) $\begin{array}{cc}{\mathrm{df}}_i.Math.f_j.Math.v.Math.=f_j^T{\mathrm{df}}_i.Math.v.Math.=\stackrel{3}{\underset{k}{.Math.}}\stackrel{3}{\underset{l}{.Math.}}{f}_{\mathrm{jk}}\frac{f_{\mathrm{ik}}}{{\mathrm{dx}}_l}{\mathrm{de}}_l.Math.v.Math.=f_j^T{J}_{i}v& \left(5\right)\\ c_{\mathrm{ijk}}=f_j^T{J}_i{f}_k& \left(6\right)\end{array}$

(27) The energy minimization computation scheme is a regularized optimization scheme. The connection forms c.sub.ijk at x.sub.0 can be obtained as the minimizer of an extrapolation energy contained within a neighborhood as given by Equation (7) where is a regularization weight used to penalize high curvature. Denoting {tilde over (f)}.sub.i as the normalized approximation to f.sub.i at x.sub.0+v using Equation (4) then we can choose to minimize the angular error between {tilde over (f)}.sub.i and f.sub.i:

(28) $\left(x_0,\right)=\frac{1}{.Math..Math.}\underset{v}{.Math.}\underset{i}{\overset{3}{.Math.}}{\mathrm{.Math.}}_i\left(x_0+v\right),$
with .sub.i(x.sub.0+v)=ar cos(f.sub.i(x.sub.0+v).Math.{tilde over (f)}.sub.i(x.sub.0+v)).

(29) $\begin{array}{cc}{c}_{\mathrm{ijk}}^{*}\left(x_0\right)=\arg \underset{c_{\mathrm{ijk}}}{\min }\left(x_0,\right)+.Math.c_{\mathrm{ijk}}.Math.& \left(7\right)\end{array}$

(30) The closed form connections in linear space computation scheme is based on trigonometrical considerations in the first-order structure of 3D frame fields and enforces the requirement that c.sub.ijv=.sub.kc.sub.ijkv.sub.k. This requirement is not explicitly enforced in the finite differentiation and energy minimization computation schemes. The closed form connections in linear space computation scheme also provides exact c.sub.ijk measurements in manifolds that have low second-order curvatures (d.sup.2f.sub.i.fwdarw.0). Given a local basis f.sub.i and data-driven neighboring bases f.sub.i(v), the one-forms c.sub.ijv can be solved for using linear least-squares. We begin by expanding Equation (4) to yield Equation (8) and analyze this expression geometrically using FIG. 3.
f.sub.i(v)=f.sub.i+c.sub.ijvf.sub.j+c.sub.ikvf.sub.k(8)

(31) Let f{(v)} denote the projection of f.sub.i(v) in the f.sub.if.sub.j plane, i.e. f.sub.i.sup.j(v)=f.sub.i(v)(f.sub.i(v).Math.f.sub.k)f.sub.k, k(1,2,3)ij and let .sub.ij(v) denote the signed angle between f.sub.1 and f.sub.i.sup.j(v) with positive values assigned to .sub.ij(v) rotating f.sub.i towards f.sub.j, obtained as Equation (9). Using trigonometry, we obtain Equation (10). Then by expanding the contraction c.sub.ijv we obtain Equations (11) to (13) which when substituted into Equation (10) yield Equation (14) wherein we have a linear system in the nine unknowns c.sub.ijk. Now, aggregating n measurements v.sub.i(x) yields Equation (15).

(32) $\begin{array}{cc}{}_{\mathrm{ij}}\left(v\right)=\mathrm{sgn}\left(f_i^j\left(v\right).Math.f_j\right).Math.\arccos \left(.Math.\frac{f_i.Math.f_i^j\left(v\right)}{.Math.f_i^j\left(v\right).Math.}.Math.\right)& \left(9\right)\\ \tan \left({}_{\mathrm{ij}}\left(v\right)\right)=\frac{c_{\mathrm{ij}}\left(v\right)}{.Math.f_i.Math.}=c_{\mathrm{ij}}.Math.v.Math.\mathrm{since}{f}_i\mathrm{is}a\mathrm{unit}\mathrm{vector}& \left(10\right)\\ \begin{array}{c}{c}_{\mathrm{ij}}.Math.v.Math.=\left(v.Math.f_1\right)c_{\mathrm{ij}1}+\left(v.Math.f_2\right)c_{\mathrm{ij}2}+\left(v.Math.f_3\right)c_{\mathrm{ij}3}\\ =v_1{c}_{\mathrm{ij}1}+c_2{c}_{\mathrm{ij}2}+v_3{c}_{\mathrm{ij}3}\end{array}& \begin{array}{c}\left(11\right)\\ \left(12\right)\end{array}\\ \mathrm{with}{v}_{k}v.Math.f_k& \\ .Math.\left[c_{12}.Math.v.Math.c_{13}.Math.v.Math.c_{23}.Math.v.Math.\right]=\left[v_1{v}_2{v}_3\right].Math.\left[\begin{array}{ccc}{c}_{121}& c_{131}& c_{231}\\ c_{122}& c_{132}& c_{232}\\ c_{123}& c_{133}& c_{233}\end{array}\right]& \left(13\right)\\ {\left[\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right]}^T\left[\begin{array}{ccc}{c}_{121}& c_{131}& c_{231}\\ c_{122}& c_{132}& c_{232}\\ c_{123}& c_{133}& c_{233}\end{array}\right]={\left[\begin{array}{c}\tan \left({}_{12}\left(v\right)\right)\\ \tan \left({}_{13}\left(v\right)\right)\\ \tan \left({}_{23}\left(v\right)\right)\end{array}\right]}^T& \left(14\right)\\ \left[\begin{array}{ccc}{v}_{11}& v_{12}& v_{13}\\ \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}\\ v_{i1}& v_{i2}& v_{i3}\\ \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}\\ v_{n1}& v_{n2}& v_{n3}\end{array}\right]\left[\begin{array}{ccc}{c}_{121}& c_{131}& c_{231}\\ c_{122}& c_{132}& c_{232}\\ c_{123}& c_{133}& c_{233}\end{array}\right]=\left[\begin{array}{ccc}\tan \left({}_{12}\left(v_1\right)\right)& \tan \left({}_{13}\left(v_1\right)\right)& \tan \left({}_{23}\left(v_1\right)\right)\\ \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}\\ \tan \left({}_{12}\left(v_i\right)\right)& \tan \left({}_{13}\left(v_i\right)\right)& \tan \left({}_{23}\left(v_i\right)\right)\\ \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}\\ \tan \left({}_{12}\left(v_n\right)\right)& \tan \left({}_{13}\left(v_n\right)\right)& \tan \left({}_{23}\left(v_n\right)\right)\end{array}\right]& \left(15\right)\end{array}$

(33) Accordingly, we require n3 otherwise the system will be undetermined. In general, v.sub.i.Math.v.sub.j0 such that V is not a full row rank. A QR factorization or singular value decomposition (SVD) pseudoinverse process can be used to solve for C C=(V.sup.TV).sup.1V.sup.T{tilde over (C)}.

(34) Referring back to FIG. 1, once the dMRI data have been modeled in accordance with the theory of moving frames, and more specifically that local orthogonal frames are fitted at each data point, as per step 102, and connection parameters c.sub.ijk are computed at each data point, as per step 104, the method 100 moves on to step 106 of performing inward diffusions into one or more regions where frame information is missing. An example embodiment of step 106 is illustrated in FIG. 4.

(35) The region of missing information may be a volumetric region or a distributed region. In some data sets, there may be a combination of volumetric regions and distributed regions for which frame information is missing. The regions may be reconstructed sequentially or in parallel.

(36) Each region is defined by a boundary and inward diffusions are performed starting from the boundary and into the region. As per step 402, each local orthogonal frame is rotated in a heading direction for the data points at the boundary of the region. New values for the three directions (i, j, k) are computed at the new data points. At step 404, connection parameters c.sub.ijk are estimated at each new data point obtained by the rotation of step 402. The estimated parameters c.sub.ijk are then used to further rotate the local orthogonal frames into the region, as per step 406.

(37) In some embodiments, the method 106 is performed iteratively, in a layer-by-layer fashion. For example, if the region to be reconstructed were to be spherical, then the layers would be concentric spherical shells of decreasing diameters. For a non-spherical region, the layers may have a shape that corresponds to the shape of the boundary of the region. The layers may be set out evenly, such that each layer has a same thickness, or they may be set out unevenly, such that the thickness varies from one layer to another. Similarly, a given layer may have a constant thickness or a thickness that varies throughout the layer. The thickness of each layer may be set as desired, for example at one voxel, two voxels, three voxels, and the like.

(38) Alternatively, the region is divided into a plurality of sub-regions, and the inward diffusions are performed separately in each sub-region. The order of diffusion for the sub-regions may be set in accordance with available data points from the dMRI data, such that a first rotation of the local frames as per step 402 is performed from known data points. This embodiment may be more suitable for distributed regions, whereas the layer-by-layer approach may be favored for volumetric regions. Other approaches for performing the inward diffusion of regions where frame information is missing may also be used.

(39) In some embodiments, rotating the local orthogonal frame comprises computing new values for the three directions (i, j, k) independently. There may be a need to adjust the newly computed values of the three directions (i, j, k) so as to enforce orthogonality, as per step 408. As illustrated in FIG. 4, step 408 may be performed after the first rotation from the data points on the boundary of the region, and/or it may be performed after each additional rotation from a new data point further into the region.

(40) Step 106 of inward diffusions into regions where frame information is missing will now be described mathematically. We define as a region where frame field data are available, and A as a reconstruction domain. Starting with A.sup.0=A each iteration propagates information on the boundary A such that A.sup.n+1=AB.sub.r is eroded with a ball element B.sub.r of radius r=1 and c.sub.ijk is updated in and A.sup.n+1 until A.sup.n+1=. Frames f.sub.i are transported from xA in a direction v across A using a neighbor accumulation of Equation (4) as given by Equation (16) which is then normalized.

(41) $\begin{array}{cc}{\tilde{f}}_i^{n+1}{|}_{x+v}=\underset{x{A}^n}{.Math.}\left(f_i^n{{.Math.}_x+\underset{ji}{.Math.}{f}_j^n.Math.}_x\underset{k}{.Math.}{c}_{\mathrm{ijk}}\left(x+v\right).Math.f_k^n{|}_x\right)& \left(16\right)\end{array}$

(42) In some embodiments, connection forms c.sub.ijk at xA are obtained by combining the three computation schemes described above, namely finite differentiation, energy minimization, and closed-form connections in linear space. For each data point, the computation scheme may be selected based on the following heuristics in Equation (17), where k=V.sup.1.sub.FV.sub.F is the condition number of V in Equation (15) and F is the Frobenius norm where N=(2r+1).sup.3 and .sub.0=3 were determined empirically and generally offer a good tradeoff between neighborhood connectivity and well-conditioning of V.

(43) 0$\begin{array}{cc}{c}_{\mathrm{ijk}}\left(x\right)=\{\begin{array}{c}\begin{array}{c}\mathrm{Closed}\mathrm{form}\mathrm{of}\mathrm{Equation}\left(15\right)\text{:}x>N\mathrm{neighbours}\\ \mathrm{and}<{}_0\end{array}\\ \begin{array}{c}\mathrm{Energy}\mathrm{of}\mathrm{Equation}\left(7\right)\text{:}=0.0001\&\mathrm{seeding}\\ \mathrm{with}\mathrm{Equation}\left(6\right)\text{:}\mathrm{else}\end{array}\end{array}& \left(17\right)\end{array}$

(44) Note that unrealistically large c.sub.ijk values can still arise in spite of the regularization. To see this, take f.sub.1(x) and make it parallel to f.sub.2(x) in a neighboring voxel. Using Equation (4) we have f.sub.1(v)f.sub.1+c.sub.12(v)f.sub.2+c.sub.13(v)f.sub.3m such that f.sub.1(v)\\f.sub.2(v)c.sub.12v.fwdarw..

(45) In some embodiments, a hard threshold may be applied on and c.sub.ijk in Equation (17). To do this we set

(46) $c_{\mathrm{ijk}}=0\mathrm{if}>\frac{}{4}$
or if c.sub.ijk exceeds bounds obtained as follows. In the discrete case using forward differences, the frame axis differential d(f.sub.i1e.sub.1+f.sub.i2e.sub.2+f.sub.i3e.sub.3) is bounded by

(47) $.Math.\frac{f_i}{x_k}.Math.\frac{1}{2}.Math.f_{\mathrm{ij}}\left(x+e_k\right)-f_{\mathrm{ij}}\left(x\right).Math.1$
and f.sub.i=1. We thus have |c.sub.ijv=|f.sub.j.sup.T.sub.f.sub.i(x.sub.1,x.sub.2,x.sub.3)v|f.sub.j.sup.T.Math.[v.sub.1v.sub.2v.sub.3].sup.Tv.sub.1, and set c.sub.ijk=min(c.sub.ijk,v.sub.1).

(48) The diffusion process guided by Equations (16) and (17) does not enforce orthogonality of the resulting frame field. Since this is a first-order method, there may be some orthogonality drift as the process gets deeper into the region A. To see this, using Equation (8), we get
f.sub.1(v).Math.f.sub.2(v)=(f.sub.1+c.sub.12(v)f.sub.2+c.sub.13(v)f.sub.3).Math.(f.sub.3+c.sub.21(v)f.sub.1c.sub.23(v)f.sub.3)

(49) and similarly for the other axis products. Since f.sub.i is by definition orthogonal at 0, we have f.sub.i.Math.f.sub.j=.sub.ij such that

(50) $f_i\left(v\right).Math.f_2\left(v\right)=\{\begin{array}{c}1+c_{\mathrm{ij}}^2+c_{\mathrm{ik}}^2\text{:}i=j\\ c_{\mathrm{ik}}{c}_{\mathrm{jk}}\text{:}ij\end{array}\mathrm{for}ki,kj.$

(51) The extrapolated frame f.sub.i(v) will therefore never be exactly orthonormal. To enforce orthonormality we therefore fix f.sub.1(v) and find its orthogonal complement f(v) using Equations (18) and (19) where c.sub.ij is taken as c.sub.ij(v). We similarly proceed with f(v).

(52) $\begin{array}{cc}f\frac{1}{2}\left(v\right)=f-\left(f_1\left(v\right).Math.f_2\left(v\right)\right)f_1\left(v\right)& \left(18\right)\\ =\left(-c_{12}-c_{13}{c}_{23}\right)f_1+\left(1-c_{13}{c}_{23}{c}_{13}\right)f_2+\left(c_{23}-c_{13}^2{c}_{23}\right)f_3& \left(19\right)\end{array}$

(53) The method 100 may be applied to reconstruction of various internal organs, such as, but not limited to the brain, the kidney and the pancreas. In some embodiments, the method 100 is applied to cardiac fiber reconstruction, i.e. the heart. For example, the dMRI data may be used to obtain a first direction for the local orthogonal frames whereas the other two directions are based on the local geometry of a heart wall. Therefore, the method 100 may be used to reconstruct a three-dimensional image of fiber orientations.

(54) Given a partial volume of fiber orientations f.sub.1 in a mask H of the heart, f.sub.1 may be reconstructed everywhere in A=H. For cardiac fiber reconstruction, step 106 of method 100 may be guided by rule-based priors for fiber orientations based on H and estimated heart wall normal directions of which one relates to the circumferential arrangement of myofibers and the other to their helix angle turning.

(55) The circumferential component may be estimated as follows. Using a smoothing kernel G.sub., the Euclidean distance transforms G.sub.*D.sub.+ and G.sub.*D.sub. to the outer and inner walls are first computed. From the average D=(D.sub.+D.sub.) local wall normal directions are computed using {tilde over (f)}.sub.3=D. The apex .sub.0 and an upward are identified, and used to obtain heart centerline measurements .sub.t parametrized over t steps along , such that:

(56) ${}_t=\frac{\underset{x}{.Math.}w\left(x\right)\left(x\right)x}{\underset{x}{.Math.}w\left(x\right)\left(x\right)x},\mathrm{where}\left(x\right)=\mathrm{sgn}\left(.Math.\left(x-{}_0-t\hat{u}\right).Math.\hat{u}.Math.\right)$

(57) In other words, (x) is 1 in the current short axis plane and 0 elsewhere, and:

(58) $w\left(x\right)=\{\begin{array}{c}1\text{:}x\mathrm{in\_the}\mathrm{\_myocardian}\\ 0\text{:}\mathrm{else}\end{array}.$

(59) Then, a smooth heart centerline is obtained as L(t)=G.sub.*.sub.t. We can now obtain a local long-axis direction f.sub.L using

(60) $f_L=\frac{L\left(t\right)}{t}$
and finally estimate the circumferential direction f.sub.C from the cross product of f.sub.L and the local wall normal {circumflex over (f)}.sub.3 as f.sub.C=f.sub.Lf.sub.3.

(61) The helical component may be estimated as follows. A rule-based helix angle variation prior is used, from .sup.+ to .sup. from outer to inner wall. A voxel x is first parametrized over the local depth of the heart wall in the range [0, 1], where 0 indicates that the voxel is lying on the outer wall and f.sub.L on the inner wall, using:

(62) $\left(x\right)=\frac{D_{+}}{D_{+}+D_{-}}\left(x\right)={\left(1+D_{-}{D}_x^{-1}\right)}^{-1}\left(x\right)\left[0,1\right]$

(63) Then, the local helix angle at x is linearly interpolated using:
a(x)=(.sub.+y.sub.++y.sub.)(x).

(64) Finally, the helix fiber direction {circumflex over (f)}.sub.1 is obtained using a helical rotation of f.sub.c about the local transmural axis {circumflex over (f)}.sub.3 from the axis angle {circumflex over (f)}.sub.3, using Rodrigues' formula:
{circumflex over (f)}.sub.1=cos f.sub.C+sin ({circumflex over (f)}.sub.3f.sub.C)+(1cos )({circumflex over (f)}.sub.3.Math.f.sub.C){circumflex over (f)}.sub.3

(65) Each diffusion pass n+1 combines current frame field estimates f.sub.i, differentials df.sub.i and rule-based priors {circumflex over (f)}.sub.i using

(66) ${\tilde{f}}_1^{n+1}={}_1{\hat{f}}_1\left(x\right)+\left(1-\right)\underset{yA}{.Math.}\left({\tilde{f}}_1^n+c_{12}.Math.v.Math.{\tilde{f}}_2^n+c_{13}.Math.v.Math.{\tilde{f}}_3^n\right)$${\tilde{f}}_1^{n+1}={}_0\left({\tilde{f}}_3-\left({\tilde{f}}_3.Math.{\tilde{f}}_1^n\right){\tilde{f}}_1^n\right)+\left(1-{}_3{\tilde{f}}_3^n\left(x\right)\right),{\tilde{f}}_2^{n+1}={\tilde{f}}_3^{n+1}{\tilde{f}}_1^{n+1}$

(67) Here, v=xy, (x) denotes the current (diffused) boundary around x from which data are recovered, and .sub.1=0.1, .sub.3=0.7 are prior weights determined empirically. The higher the confidence in the rule-based model, the larger these coefficients should be. Each {acute over (f)}.sub.i.sup.n+1 is normalized after each diffusion pass.

(68) Damaged diffusion volumes were simulated using Poisson disk stochastic sampling, where each sample point p satisfies a minimum distance constraint to others. At p, an ellipsoid with random semi-axis lengths (range=: 1 to 10 voxels) is carved out. A prototypical synthetic in vivo mask was also obtained by regularly slicing H along its long-axis. The corruptions were applied to dMRI volumes of a healthy rat heart. Subsequently, the method 100 was compared against a standard vector interpolation scheme based on distance weighting, against a pure vector diffusion scheme using Equation (16) with c.sub.ijk=0, and against a ruled-based model. The robustness to noise was also tested by combining Poisson sampling and random angular perturbations to f.sub.1, prior to reconstruction.

(69) FIG. 5 depicts the reconstruction of synthetic cardiac fiber fields respectively generated in a short axis slab near the mid-section from which a chunk of tissue was removed, and from a synthetic transmural cardiac sample obtained using a rule-based helix angle variation (total turning 120 degrees). The reconstruction method 100 shows a significantly reduced angular error compared to other methods.

(70) FIG. 6 shows Poisson and slicing error maps using the method 100 according to one embodiment. The error increases with the amount of damage, and is higher near boundaries.

(71) Referring to FIG. 7 there are compared mean reconstruction errors for the different methods. The method 100, referred to in FIG. 7 as Frame Inpainting, outperforms the other methods in minimizing the reconstruction error in all corruption scenarios, and also shows comparable robustness to noise.

(72) Turning to FIG. 8, there is illustrated a dMRI data reconstruction system 802. The system 802 may be implemented as a downloaded software application, a firmware application, or a combination thereof. The system 802 may be accessible remotely from any one of a plurality of devices 806 over connections 804. The devices 806 may comprise any device, such as a personal computer, a tablet, a smart phone, or the like, which is configured to communicate over the connections 804. In some embodiments, the system 802 may itself be provided directly on one of the devices 806, either as a downloaded software application, a firmware application, or a combination thereof.

(73) The connections 804 may comprise wire-based technology, such as electrical wires or cables, and/or optical fibers. The connections 804 may also be wireless, such as RF, infrared, Wi-Fi, Bluetooth, and others. Connections 804 may therefore comprise a network, such as the Internet, the Public Switch Telephone Network (PSTN), a cellular network, or others known to those skilled in the art. Communication over the network may occur using any known communication protocols that enable devices within a computer network to exchange information. Examples of protocols are as follows: IP (Internet Protocol), UDP (User Datagram Protocol), TCP (Transmission Control Protocol), DHCP (Dynamic Host Configuration Protocol), HTTP (Hypertext Transfer Protocol), FTP (File Transfer Protocol), Telnet (Telnet Remote Protocol), SSH (Secure Shell Remote Protocol).

(74) One or more databases 808 may be integrated directly into the system 802 or any one of the devices 806, or may be provided separately therefrom (as illustrated). In the case of a remote access to the databases 808, access may occur via connections 804 taking the form of any type of network, as indicated above. The various databases 808 described herein may be provided as collections of data or information organized for rapid search and retrieval by a computer. The databases 808 may be structured to facilitate storage, retrieval, modification, and deletion of data in conjunction with various data-processing operations. The databases 808 may be any organization of data on a data storage medium, such as one or more servers. The databases 808 illustratively have stored therein any one of datasets of dMRI data, computed local orthogonal frame directions (i, j, k), c.sub.ijk parameters, data defining boundaries of regions to reconstruct, reconstructed images, partially reconstructed images, and any other information used for the methods described herein.

(75) FIG. 9 is an example embodiment of the system 702. A data modeling unit 902 is configured to receive input data, such as dMRI data, and perform initial modeling steps of the data in accordance with the theory of moving frames, as per steps 102 and 104 of the method 100. Once the dMRI data are modeled, the data are transmitted to a data recovery unit 904, configured to perform step 106 of method 100. In some embodiments, the boundaries of regions to be reconstructed are defined by the data modeling unit 902 and provided to the data recovery unit 904. Alternatively, the data recovery unit 904 is configured to identify the regions to be reconstructed and define the boundaries.

(76) FIG. 10 is an example embodiment of the data modeling unit 902, which illustratively comprises a fitting unit 1002 and a connection parameters computation unit 1004. The fitting unit 1002 may be configured to receive the dMRI data and compute the local orthogonal frame directions at each data point of the dMRI data. The connection parameters computation unit 1004 may be configured to compute the connection parameters c.sub.ijk in accordance with any one of the computation schemes described above, namely finite differentiation, energy minimization, and closed-form connections in linear space. Other computation schemes may also be used. In some embodiments, the closed-form connections computation scheme is used exclusively for computing the connection parameters c.sub.ijk of the known data, whereas a combination of computation schemes is used to estimate the connection parameters c.sub.ijk in the regions to be reconstructed. The connection parameters computation unit 1004 outputs a partially reconstructed image, with or without the boundaries of regions for reconstruction having been defined.

(77) FIG. 11 is an example embodiment of the data recovery unit 904. There is provided a frame rotation unit 1102 and a connection parameters estimation unit 1104. The frame rotation unit 1102 receives the partially reconstructed image, which may be represented as the computed frame directions (i, j, k) at each known data point and corresponding connection parameter c.sub.ijk. Alternatively, the information input to the frame rotation unit 1102 is composed only of data points found at the boundaries of the regions to be reconstructed. Irrelevant data points, i.e. those not needed to perform the steps of rotation 402 and estimation 404, may be omitted. The connection parameters estimation unit 1104 is configured to receive the newly computed frame values after a rotation, and estimate connection parameters c.sub.ijk for the newly computed frame values. As indicated above, the process may be performed iteratively until each region has been fully reconstructed, or until a threshold has been reached regarding recovered data. A reconstructed image may be output.

(78) Note that while illustrated separately, the connection parameters computation unit 1004 and the connection parameters estimation unit 1104 may be provided as a single unit that can perform both computation based on obtained dMRI data and estimation based on rotated local orthogonal frames. Similarly, the two units 1004 and 1104 may share certain resources. Other units may also be combined or separated as desired to suit the targeted functions.

(79) FIG. 12 shows a schematic representation of an example implementation of the method 100 as a combination of software and hardware components. A computing device 1200 is illustrated with one or more processing units (referred to as the processing unit 1202) and one or more computer-readable memories (referred to as the memory 1204) having stored thereon program instructions 1206 configured to cause the processing unit 1202 to generate one or more outputs based on one or more inputs. The inputs may comprise the dMRI data and any additional information needed to perform the method 100. The outputs may comprise partially reconstructed images, fully reconstructed images, and any other information resulting from performing any step of the method 100.

(80) The processing unit 1202 may comprise any suitable devices configured to cause a series of steps to be performed so as to implement the computer-implemented method 100 such that instructions 1206, when executed by a computing device 1200 or other programmable apparatus, may cause the functions/acts/steps specified in the methods described herein to be executed. The processing unit 1202 may comprise, for example, any type of general-purpose microprocessor or microcontroller, a digital signal processing (DSP) processor, a central processing unit (CPU), an integrated circuit, a field programmable gate array (FPGA), a reconfigurable processor, other suitably programmed or programmable logic circuits, or any combination thereof.

(81) The memory 1204 may comprise any suitable known or other machine-readable storage medium. The memory 1204 may comprise non-transitory computer readable storage medium such as, for example, but not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any suitable combination of the foregoing. The memory 1204 may include a suitable combination of any type of computer memory that is located either internally or externally to device such as, for example, random-access memory (RAM), read-only memory (ROM), compact disc read-only memory (CDROM), electro-optical memory, magneto-optical memory, erasable programmable read-only memory (EPROM), and electrically-erasable programmable read-only memory (EEPROM), Ferroelectric RAM (FRAM) or the like. Memory may comprise any storage means (e.g., devices) suitable for retrievably storing machine-readable instructions executable by processing unit.

(82) Each computer program described herein may be implemented in a high level procedural or object-oriented programming or scripting language, or a combination thereof, to communicate with a computer system. Alternatively, the programs may be implemented in assembly or machine language. The language may be a compiled or interpreted language. Computer-executable instructions may be in many forms, including program modules, executed by one or more computers or other devices. Generally, program modules include routines, programs, objects, components, data structures, etc., that perform particular tasks or implement particular abstract data types. Typically the functionality of the program modules may be combined or distributed as desired in various embodiments.

(83) While illustrated in block diagrams as groups of discrete components communicating with each other via distinct data signal connections, it will be understood by those skilled in the art that the present embodiments are provided by a combination of hardware and software components, with some components being implemented by a given function or operation of a hardware or software system, and many of the data paths illustrated being implemented by data communication within a computer application or operating system. The structure illustrated is thus provided for efficiency of teaching the present embodiment.

(84) The above description is meant to be exemplary only, and one skilled in the relevant arts will recognize that changes may be made to the embodiments described without departing from the scope of the invention disclosed. For example, the blocks and/or operations in the flowcharts and drawings described herein are for purposes of example only. There may be many variations to these blocks and/or operations without departing from the teachings of the present disclosure. For instance, the blocks may be performed in a differing order, or blocks may be added, deleted, or modified.

(85) The present disclosure may be embodied in other specific forms without departing from the subject matter of the claims. Also, one skilled in the relevant arts will appreciate that while the systems, methods and computer readable mediums disclosed and shown herein may comprise a specific number of elements/components, the systems, methods and computer readable mediums may be modified to include additional or fewer of such elements/components. The present disclosure is also intended to cover and embrace all suitable changes in technology. Modifications which fall within the scope of the present disclosure will be apparent to those skilled in the art, in light of a review of this disclosure, and such modifications are intended to fall within the appended claims.