Abstract:
In a method and apparatus for determination of coefficients of a diffusion tensor to describe a diffusion process within a subject, spatially-resolved, variously diffusion-coded diffusion data are generated from volume elements of the subject under consideration of control data for various diffusion codings. Intermediate diffusion data are generated for each volume element from diffusion data measured with a first diffusion coding. Updated intermediate diffusion data are then determined for each volume element from diffusion data generated preceding the intermediate diffusion data and current diffusion data measured with a different diffusion coding. Components of a diffusion tensor are determined from the last determined, updated intermediate diffusion data after multiple updates with various diffusion codings.

Description:
BACKGROUND OF THE INVENTION  
       [0001]     1. Field of the Invention  
         [0002]     The present invention concerns a method for determination of coefficients of a diffusion tensor by means of magnetic resonance for description of a diffusion process within a subject, in which spatially-resolved, variously diffusion-coded diffusion data are generated from volume elements of the subject dependent on control data for various diffusion codings in chronological order.  
         [0003]     The invention likewise concerns a device for implementation of the method.  
         [0004]     2. Description of the Prior Art  
         [0005]     A method and a device of the aforementioned type are known from U.S. Pat. No. 5,539,310. The method specified therein is based on a magnetic resonance measurement sequence proposed by Stejskal and Tanner in 1965. The method described in U.S. Pat. No. 5,539,310 expands the measurement sequence proposed by Stejskal and Tanner such that coefficients or elements of a diffusion tensor are determined and graphically represented for each voxel in an examination region. A clear representation with a diffusion ellipsoid, the main axis of which represents the direction of the strongest diffusion process, ensues. The expansion in the individual directions stands for a numerical value of the diffusion process in the corresponding directions. In medical applications, the relative mobility of water molecules in endogenic tissue can thus be measured with magnetic resonance diffusion imaging. Since the diffusion in the tissue can depend on its structure (such as, for example, the fiber direction), medically-relevant conclusions can be derived from the diffusion tensor.  
         [0006]     The diffusion data are determined in the magnetic resonance measurement from the size and direction of the diffusion gradient fields used for diffusion coding. Strong gradient pulses are thereby used that are oriented symmetrically to a radio-frequency, 180° refocusing pulse. The first gradient pulse before the 180° refocusing pulse generates a phase shift for all spins; the second gradient pulse inverts this phase shift. Given stationary molecules (protons in medical imaging), the phase shift therewith cancels again. However, for molecules that, due to Brownian motion, are located at a different location during the effect of the second gradient pulse than during the effect of the first gradient pulse, the phase shift is not completely compensated. A rest phase displacement remains that leads to a signal attenuation. The diffusion coding can be controlled by the size and direction of such gradient pulses.  
         [0007]     In medical diffusion tensor measurement by means of magnetic resonance techniques and the subsequent graphical representation (DTI=Diffusion Tensor Imaging), large quantities of measurement data accumulate, from which the sought six tensor parameters or tensor coefficients (thus the independent elements or components of a symmetrical 3×3 tensor matrix) are calculated per voxel. The tensor parameters are subsequently used to calculate relevant parameter cards for the diagnostics. For example, the isotropic portion of the diffusion tensor or the anisotropic portion of the diffusion tensor is displayed in corresponding parameter cards (Average Apparent Diffusion Coefficient Map or ADCav Map, or Fractional Anisotropy Map or FA Map). Due to the large amount of measurement data, the calculation of the diffusion tensor per voxel requires the determination of the unknown parameters by means of a compensation. Methods known from multi-variant linear regression are used, for example methods that form a pseudo-inverse or implement a singular value decomposition. These methods are, however, very storage-space intensive and computationally complex because all data for the evaluation are retained in the known methods and the sought parameters are only determined by means of compensation methods after conclusion of the measurement. Rules for accounting for the intensity values are also applied as alternatives for specific, rigid sets of diffusion codings. However, such methods are very inflexible.  
       SUMMARY OF THE INVENTION  
       [0008]     An object of the present invention is to provide a method and a device in which a flexible evaluation in DTI is possible and the required storage space capacity is reduced.  
         [0009]     The object is achieved in accordance with the invention by a method wherein intermediate diffusion data are thereby generated for each volume element from diffusion data measured with a first diffusion coding, updated intermediate diffusion data are determined for each volume element from diffusion data generated preceding the intermediate diffusion data and current diffusion data measured with a different diffusion coding, and components of a diffusion tensor are determined from the last determined, updated intermediate diffusion data after multiple updates with various diffusion codings.  
         [0010]     The incoming measurement values can be processed immediately and used for updating the intermediate diffusion data with this incremental method. The tensor components can be determined from the last updated intermediate data with small processing steps that are not very calculation-intensive. It is accordingly not necessary (as in the prior art) for all measurement values to be initially determined and stored in memory before the evaluation can begin. The necessary storage requirement is significantly reduced in the inventive method. Parts of the calculations can already be implemented in parallel with the data acquisition; the time required after the measurement for the calculation of the results is likewise significantly reduced. The calculations still remaining after conclusion of the measurement no longer directly use to the measurement data. The calculation is deconstructed into parts which require an access to the measurement data and parts that can be started with intermediate results.  
         [0011]     In one embodiment the intermediate diffusion data and the updated intermediate diffusion data are determined by means of a normal smoothing method using the Moore-Penrose pseudo-inverse. For each partial measurement, thus for each measurement with a specific diffusion coding, this method requires only information for and from this partial measurement. It is therefore in principle possible to determine diffusion components even after an unforeseen termination of the measurement. The precision of the results then naturally depends on the quality of the measurement data.  
         [0012]     In a further embodiment the intermediate diffusion data and updated intermediate diffusion data are determined by means of a singular value decomposition method. Here information is in fact required regarding all diffusion codings, such as the number of the measurements times the six diffusion components to be determined; in spite of this, an incremental processing of the measurement data is also possible here with a correspondingly-reduced memory and computation requirement.  
         [0013]     In both of these embodiments, standard methods for matrix calculation from linear algebra can be used that are available in highly-optimized standard libraries or that can be embodied in hardware extensions of modern processors.  
         [0014]     In a further embodiment, a diffusion-coding effect of imaging gradient fields is additionally taken into account to determine the components of the diffusion tensor. An improvement of the quality of the diffusion tensor components so determined is connected with this since all gradient fields influencing the diffusion coding are taken into account as well.  
         [0015]     In another embodiment, the diffusion data are generated in chronological order such that diffusion data that strongly determine the components of the diffusion tensor are generated first and those that weakly determine the components of the diffusion tensor are generated afterwards. A first approximation for the components of the diffusion tensor thus can already be calculated after a short measurement time. Such a measurement order is additionally tolerant of an early measurement termination.  
         [0016]     In a further embodiment, a signal is output when the components of the diffusion tensor can be determined with a predetermined precision. In another embodiment the performance of the components is determined and is output as a quality measure.  
         [0017]     In a further embodiment, a condition number is determined from the control data as a measure for the determinability of the components. 
     
    
     DESCRIPTION OF THE DRAWINGS  
       [0018]      FIG. 1  schematically shows a diagnostic magnetic resonance apparatus with the basic components for determination of a diffusion tensor and,  
         [0019]      FIG. 2  is a block diagram showing the basic steps of a method for determination of components of a diffusion tensor in accordance with the invention. 
     
    
     DESCRIPTION OF THE PREFERRED EMBODIMENTS  
       [0020]     The diffusion behavior of living tissue gives important information regarding its microstructure. Such diffusion processes can be described by a diffusion tensor the components of which are determined with the aid of magnetic resonance technique. Diffusion-coded or diffusion-weighted magnetic resonance signals can be generated with diagnostic magnetic resonance apparatuses known today.  
         [0021]     The diffusion-weighted magnetic resonance signal can thus be described according to the following formula. 
 
 S=S   0   e   −bD    (1) 
    S 0  is thereby the initial signal strength without diffusion-weighted gradient pulses     b is thereby the diffusion weighting factor     D is the diffusion coefficient, also designated as ADC (apparent diffusion coefficient).    
 
         [0025]     Only tissue with isotropic diffusion properties can be described in the diffusion coefficients D according to equation (1). In order to also be able to detect the case of an anisotropic diffusion, the equation 1 is expanded as follows: 
 
 S=S   0   e   −bg     Dg      (2) 
 
 wherein D is the symmetrical 3×3 diffusion tensor.  
               D   =     (           D   xx           D   xy           D   xz                           D   yy           D   yz                                       D   zz           )       ⁢     
     ⁢   and           (   3   )               g   =       g   diff     /          g   diff                    (   4   )             
 
 describe the normalized diffusion-coded gradient vector. D xx , D xy , D xz , D yy , etc. are the individual components or elements of the diffusion tensor. These are also sometimes designated as tensor parameters. 
 
         [0026]     The effect of the gradients used can be described by the introduction of a diffusion matrix B. This is developed as follows:  
                   S   =       S   0     ⁢     ⅇ       -     bg   T       ⁢   Dg                     =       S   0     ⁢     ⅇ       -   b     ⁢     〈       gg   T     ,   D     〉                       =       S   0     ⁢     ⅇ     -     〈     bG   ,   D     〉                       =       S   0     ⁢     ⅇ     -     〈     B   ,   D     〉                       =       S   0     ⁢     exp   ⁡     [     -       ∑     k   =   1     3     ⁢       ∑     l   =   1     3     ⁢       B   kl     ⁢     D   kl             ]                       (   5   )             
 
         [0027]     Instead of the index of the coordinate directions x, y and z, a number index 1, 2 and 3 has been analogously employed here for reasons of the simpler representation in a formula. With the B matrix it is possible to take all diffusion imaging gradients as well as their cross terms into account in the measurement sequence.  
         [0028]     At least six measurements with various diffusion coding directions are necessary in order to determine the diffusion tensor. The elements of the diffusion tensor and the T2-weighted signals S 0  for each voxel from the measurement region or the examination region can be determined with an additional measurement without the effect of diffusion gradients (baseline image).  
         [0029]     In the event that more than seven measurements are implemented (this is desirable for reasons of the then-better signal-to-noise ratio), the coefficients are determined such that the errors then resulting are minimal. The following correlations thereby result. 
 
 Measurements or Observations:  
                         Number   ⁢           ⁢   of   ⁢           ⁢   the   ⁢           ⁢   measurements   ⁢     :           n                                 (   6   )                     Measurement   ⁢           ⁢   vector   ⁢     :               1   T     =       (       l   1     ⁢   ⋯   ⁢           ⁢     l   n       )     =     (     ln   ⁢           ⁢     S   1     ⁢   ⋯   ⁢           ⁢   ln   ⁢           ⁢     S   n       )                     (   7   )               Measurement   ⁢           ⁢   equation   ⁢     :     ⁢           ⁢               ln   ⁢           ⁢     S   i       +     v   i       =       ⁢       ln   ⁢           ⁢     S   0       -     〈       B   i     ,   D     〉                   =       ⁢       ln   ⁢           ⁢     S   0       -       B     xx   i       ⁢     D   xx       -       B   yy     ⁢     D   yy       -                     ⁢         B     zz   i       ⁢     D   zz       -     2   ⁢     B     xy   i       ⁢     D   xy       -                     ⁢       2   ⁢     B             ⁢     xz             ⁢   i           ⁢     D             ⁢   xz         -     2   ⁢     B     yz   i       ⁢     D   yz                                               (   8   )                   
 
 Unknowns: 
 
Number of the unknowns: u=7   (9) 
 
 Vector of the Unknowns: 
 
 x   T =( x   1    . . . x   u )=( D   xx   D   xy   D   xz   D   yy   D   yz   D   zz ln S   0 )   (10) 
 
Redundancy of the Equation System:  r=n−u    (11) 
 
 Solution Models: 
 
Model function:  I+v=Ax    (12) 
 
 Target Function for the Improvement v: 
 
v T v→min   (13) 
 
 Solution 
 
 x =( A   T   A ) −1   A   T   I=N   −1   n    (14) 
 
 Design Matrix:  
               A   ij     =     (       δ   ⁢           ⁢     l   i         δ   ⁢           ⁢     x   j         )             (   15   )             
 
 with the partial derivatives  
                       δ   ⁢           ⁢     l   i         δ   ⁢           ⁢     x   1         =         δ   ⁢           ⁢     l   i         D   xx       =     -     B   xx                     δ   ⁢           ⁢     l   i         δ   ⁢           ⁢     x   2         =         δ   ⁢           ⁢     l   i         D   xy       =       -   2     ⁢     B   xy                     δ   ⁢           ⁢     l   i         δ   ⁢           ⁢     x   3         =         δ   ⁢           ⁢     l   i         D   xz       =       -   2     ⁢     B   xz                         δ   ⁢           ⁢     l   i         δ   ⁢           ⁢     x   4         =         δ   ⁢           ⁢     l   i         D   yy       =     -     B   yy                     δ   ⁢           ⁢     l   i         δ   ⁢           ⁢     x   5         =         δ   ⁢           ⁢     l   i         D   yz       =       -   2     ⁢     B   yz                     δ   ⁢           ⁢     l   i         δ   ⁢           ⁢     x   6         =         δ   ⁢           ⁢     l   i         D   zz       =     -     B       z   .     ⁢     z   .                             δ   ⁢           ⁢     l   i         x   7       =         δ   ⁢           ⁢     l   i         δ   ⁡     (     ln   ⁢           ⁢     S   0       )         =   1                                                                                                               (   16   )                                                                                     
 
 The diffusion tensor is normally specified in a patient coordinate system. For continuing analyses it can be advantageous to implement a coordinate transformation in a coordinate system whose main axes correspond to the directions of the eigenvector of the tensor. For this purpose it is necessary to implement an eigenvalue-eigenvector decomposition of the diffusion matrix D.  
                   D   =       (           e   1           e   2           e   3           )     ⁢     (           λ   1         0       0           0         λ   2         0           0       0       λ3         )     ⁢       (           e   1           e   2           e   3           )     T                   =         λ   1     ⁢     e   1     ⁢     e   1   T       +       λ   2     ⁢     e   2     ⁢     e   2   T       +       λ   3     ⁢     e   3     ⁢     e   3   T                       (   17   )             
 
 with 
        e i : uniform eigenvector in the directions 1, 2, 3     e 1 : main diffusion directions     λ i : diffusity in the direction of e i , (λ 1 ≧λ 2 ≧λ 3 )        
 
         [0033]     The inverse normal equation matrix N −1  contains information about the precision of the unknown tensor elements and the estimated S 0  values  
                 Σ   xx     =       1     σ   0   2       ⁢     N     -   1           ⁢     
     ⁢   and           (   18   )                 σ   0   2     =         v   T     ⁢   v       n   -   u               (   19   )             
 
         [0034]     The matrix Σ xx  is entirely determined by the gradient pulses for diffusion coding. It supplies information about the form of the uncertainties or the spatial conformation (shape) of the uncertainties.  
         [0035]     The scalar σ 0   2  changes from pixel to pixel and supplies information about the performance of the adaptation of the tensor model with the measurement data.  
         [0036]     A condition number of the normal equation matrix can be used as a measure for whether the equation (14) leads to a stable solution with the present data. The calculation of the condition number merely requires the control data for the diffusion coding. A matrix is singular (thus not solvable) when its condition number is infinite. It is poorly conditioned when its condition number is large. This is synonymous with the situation that the reciprocity value of the condition number approximates the floating point precision of the calculation method used. A higher value for the condition number therewith indicates an assertion about a poor planning quality of the measurement method. Known methods can be used to determine the condition number.  
         [0037]     A first realization of the incremental determination of the matrix N and of the vector n is derived in the following. First introduce a vector a i ; i=1 . . . n  
               a   i   T     =     (             δ   ⁢           ⁢     l   i         δ   ⁢           ⁢     x   1                 δ   ⁢           ⁢     l   i         δ   ⁢           ⁢     x   1             ⋯           δ   ⁢           ⁢     l   i         δ   ⁢           ⁢     x   1               )             (   20   )             
 
 is defined using the partial derivations-of the equations (16). 
 
         [0038]     After each measurement of the examination area the matrix n i  can be updated because the correlation 
 
 N   i   =N   i-1   +a   i   a   i   T    (21) 
 
 exists. The vector n can additionally be updated for each voxel according to the correlation 
 
 n   i ( x,y,z ,)= n   i-1 ( x,y,z )+ln S   i ( x,y,z ) a   i    (22) 
 
         [0039]     The required updates after each measured image or volume (i.e. after each diffusion coding) are composed of only some few simple operations involving the vector a i  and the matrix N i , including calculating seven multiplications and additions per voxel. The reason is in that the vector a i  remains unchanged per measured image or volume. These calculation operations can also be implemented during the measurement (inline).  
         [0040]     After all images/volumes are measured with the individual diffusion codings and the measurement values are further processed corresponding to the equations (21) and (22) specified in the preceding, the matrix N i  is identical to the matrix N and the vector n i  is identical to the vector n in the equation (14). After inversion of the matrix N, the unknown tensor components can be determined for each voxel corresponding to equation (14).  
         [0041]     A value for the precision of the current intermediate diffusion data can likewise be determined as follows. The variance of the unit of weight can be calculated per pixel, as is indicated in the preceding in equation (19). In the event that the components of the vector v are not required, an efficient incremental solution for this can be specified. Starting from the equations (12) and (19),  
                       v   T     ⁢   v     =         (     Ax   -   1     )     T     ⁢     (     Ax   -   1     )                   =           x   T     ⁢     A   T     ⁢   Ax       ︸   I       -         2   ⁢     x   T         ︸   II       ⁢         A   T     ⁢   1       ︸   III         +         1   T     ⁢   1       ︸   IV                       (   23   )             
 
 can be derived. The discussion of the equation (23) supplies four components: 
        I: x T A T Ax This component contains the normal equation matrix N=A T A and the solution vector x. Both the matrix and the solution vector have already been calculated in connection with the incremental solution that is specified in equation 20-22.        
 
         [0043]     II: 2x T  This component contains the solution vector x which, as already specified above, was already determined as part of the incremental solution.  
         [0044]     III: A T I This component has likewise already been calculated as above.  
         [0045]     IV: I T I Only this component must be additionally calculated and stored.  
         [0046]     As stated in the preceding, the solution for {circumflex over (σ)} 0   2  can be calculated incrementally. The single additional data that must still be determined is a sum of the quadratic values of In S i (x, y, z).  
         [0047]     In a further embodiment, a singular-value decomposition (SVD) is used for incremental calculation of the elements of the diffusion tensor. The singular-value decomposition also finds a minimal error quadratic solution of an over-determined equation system. In principle, each M·N matrix A whose number of the rows M is greater than or equal to the number of the columns N can be specified as a product of an M·N column-orthogonal matrix U, an M·N diagonal matrix W with positive or zero elements (the singular values) and the transposes of an N·N orthogonal matrix V.  
               (   A   )     =       (   U   )     ·     (           w   1                                                                           w   1                                                                         ⋰                                                                       ⋰                                                                         w   N           )     ·     (     V   T     )               (   24   )             
 
         [0048]     The smallest error quadratic solution for equation (12) is specified as follows with equation (25), whereby the matrix A +  represents what is known as the Moor-Penrose or pseudo-inverse matrix:  
               (                                         x                                         )     =         (                                         V                                         )     ·     (                                                                                                                                                                   diag   ⁡     (     1   /     w   j       )                                                                                                                                                                     )     ·     (                                                                                                                                                                   U   T                                                                                                                                                                   )     ·     (                                                         1                                                         )       =       A   +     ·   1               (   25   )             
 
         [0049]     As has already been described in the preceding, the elements of the matrix A can be determined in equation 12 from the control data for the diffusion measurement, thus from the measurement protocol. These data are independent of the measurement values. The pseudo-inverse matrix A +  can be determined before the actual data acquisition as a result of thus. The solution vector x can now be calculated incrementally and pixel-by-pixel from the measurement data as follows, 
 
 x   i ( x,y,z, )= x   i−1 ( x,y,z )+l n S   i ( x,y,z )· a   i   +   (26) 
 
 wherein a i   +  designates the i-th column of the pseudo-inverse matrix A + . After the measurement has been completed with the last diffusion coding, the solution vector simultaneously exists. 
 
         [0050]     A diagnostic magnetic resonance apparatus  1  with the basic components to determine a diffusion tensor will now be described according to the preceding explanation of the bases of two embodiments and using the schematic representation in  FIG. 1 . Suitable diagnostic magnetic resonance apparatuses are, for example, commercially available from Siemens Medical Solutions under the name Magneton Symphony, Magneton Avanto, Magneton Trio, Magneton Espree. The diagnostic magnetic resonance apparatus  1  has a basic field magnet  2  that generates a homogeneous magnetic field of high field strength. To excite and to acquire (receive) magnetic resonance signals, a radio-frequency system  4  is provided that is operated in a sequence determined by the strength of the basic magnetic field as well as the gyromagnetic ratio of the atomic nucleus used for imaging. Predominantly hydrogen nuclei (protons) that are present in large number in living tissue are used in medical imaging for magnetic resonance imaging and also for determination of the diffusion tensor. With a gradient field system  6 , the nuclear spins are spatially-coded by switched magnetic gradient fields during the excitation, after the excitation and during reception (acquisition). A controller  8  (that is realized to the greatest possible extent via a computer architecture with corresponding software) controls the entire measurement workflow from the generation of the magnetic resonance examination signal in the examination region to the acquisition of the magnetic resonance signals radiated from the examination region again.  
         [0051]     As mentioned in the preceding, the measurement sequences can be fashioned such that diffusion-weighted or diffusion-coded magnetic resonance signals are generated. The magnetic resonance apparatus is correspondingly fashioned in terms of its controller  8 . Control data  10  for the gradient fields used for diffusion weighting are predetermined corresponding to the direction and strength of the diffusion coding. Magnetic resonance signals acquired from the examination region with the corresponding diffusion coding are here supplied as diffusion data  12  to an evaluation unit that calculates a diffusion tensor from a complete set of measured diffusion data, as is described in detail in the following.  
         [0052]     After a start  20  of the measurement method, the control data are initially predetermined for a specific number of diffusion codings  22 . In the method it is checked in parallel whether a condition number  24  of the control data used for the diffusion coding is sufficient in order to obtain a stable solution for the components of the diffusion tensor to be determined. The control data represent specifications for the strength, duration and the temporal interval of the gradient pulses used for diffusion coding, which are then converted by the magnetic resonance apparatus into a corresponding diffusion-weighted pulse sequence. The spatial coding of the diffusion signals ensues in a typical manner with gradient fields in the excitation, in the acquisition and in-between. The spatially-coded measurement data  12  determined for each diffusion coding are prepared regarding diffusion data in order to update intermediate diffusion data  28  and thus to generate current intermediate diffusion data  30  that take into account all previously-implemented diffusion codings. A value for the precision of the current intermediate diffusion data is also currently calculated as was explained in the preceding. Given the attainment of sufficient precision  31 , the matrix of the intermediate diffusion data is inverted in order to then determine the components of the diffusion tensor  14 .  
         [0053]     The components of the diffusion tensor are then displayed in a known representation (for example by means of diffusion ellipsoids) on an image display device corresponding to the location in the examination region.  
         [0054]     Although modifications and changes may be suggested by those skilled in the art, it is the intention of the inventors to embody within the patent-warranted hereon all changes and modifications as reasonably and properly come within the scope of their contribution to the art.