METHOD FOR OPTIMIZED BIAS AND SIGNAL INFERENCE IN MAGNETIC RESONANCE IMAGE ANALYSIS

Abstract

An approach to estimate noise, Rician signal bias and true signal in magnitude signal data obtained with magnetic resonance imaging. The method uses multiple measurements at different scan parameter settings, also referred to as weightings, and an iterative algorithm to estimate noise, expected signal and associated Rician signal bias. Measurements at all measured weighting levels contribute to the ultimate estimation of the bias-free signal decay function. Therefore, of the so processed magnetic resonance image data, weighted signals can be computed at arbitrary weighting levels and with considerably better signal-to-noise ratio than the originally obtained data at corresponding weightings. Bias-free weighted image data at desired weighting levels, maps of the decay function fit parameters, or maps of a combination of such decay function parameters can be used for rapid and highly sensitive tissue characterization.

Claims

1. A method for highly reliable non-mono-exponential fitting of diffusion-weighted images, said method comprising the steps of: a) providing a magnetic resonance imaging apparatus capable of operation at diffusion encoding levels (b-factors) within a range ranging at least between an upper and a lower limit; b) acquiring a set of magnitude signal decays from a patient using said magnetic resonance imaging apparatus at each of a selected plurality of encoding levels distributed across the range of b-factors within its capability, each said set of signal decays being representative of an image of a preselected cross-section of the patient, and each said signal decay of each said set corresponding to a pixel of its associated image; c) processing said acquired or bias-corrected signal decays on a pixel by pixel basis to obtain the best possible fit between them and a non-mono-exponential equation wherein a non-linear least-squares method is used to obtain a fit.

2. The method of claim 1, wherein said lower limit is below 500 sec/mm.sup.2.

3. The method of claim 1, wherein said lower limit is below 300 sec/mm.sup.2.

4. The method of claim 1, wherein said lower limit is below 150 sec/mm.sup.2.

5. The method of claim 1, wherein said upper limit is above 2000 sec/mm.sup.2.

6. The method of claim 1, wherein said upper limit is above 3000 sec/mm.sup.2.

7. The method of claim 1, wherein said upper limit is about 5000 sec/mm.sup.2 or higher.

8. The method of claim 1, wherein the processing of said acquired or bias-corrected signal decays on a pixel by pixel basis in step 1c) is made to obtain the best possible fit between them and a non-mono-exponential equation of the form: $S\left(b\right)=S_0\exp \left(-b.Math.\mathrm{AD}{C}_K+b^2.Math.\mathrm{AD}{C}_K^2.Math.\frac{K_1}{6}-b^3.Math.\mathrm{AD}{C}_K^3.Math.\frac{K_2^2}{27}\right)$ wherein, in this equation, S.sub.0 is a constant that depends on constant values associated with (a) the magnetic imaging apparatus, (b) the spin-spin relaxation time T2, (c) the spin lattice relaxation time T1, and (d) the spin density ρ, wherein b is the diffusion encoding level, preferably determined for each pixel location, and wherein ADC.sub.K is a parameter for the apparent kurtosis diffusion coefficient, and wherein, K.sub.1 is a parameter for the excess kurtosis, and wherein, K.sub.2 is a second parameter for the excess kurtosis, and wherein K.sub.1 and K.sub.2 are expressed as predefined functions of ADC.sub.K.

9. The method according to claim 1, further using step of claim 8 to compute images with signal values between sampled values and outside the acquired range.

10. The method according to claim 1, wherein results obtained are displayed such that a zero value is displayed as a black pixel, and the larger the value the brighter the display of the associated pixel becomes.

11. An apparatus for highly reliable non-mono-exponential fitting of diffusion-weighted images, said apparatus comprising a controller arranged to: a) receive a set of magnitude signal decays from a patient obtained from a magnetic resonance imaging apparatus at each of a selected plurality of encoding levels distributed across a range of b-factors within its capability, each said set of signal decays being representative of an image of a preselected cross-section of the patient, and each said signal decay of each said set corresponding to a pixel of its associated image; b) wherein said range of b-factors ranges at least between an upper and a lower limit, said lower limit preferably being below 500 sec/mm.sup.2, and more preferably below 300 sec/mm.sup.2, and most preferably below 150 sec/mm.sup.2, such as about 100 sec/mm.sup.2 or lower, and the upper limit preferably being above 2000 sec/mm.sup.2, and more preferably above 3000 sec/mm.sup.2, and most preferably about 5000 sec/mm.sup.2 or higher; and c) processing said acquired or bias-corrected signal decays on a pixel by pixel basis to obtain the best possible fit between them and a non-mono-exponential equation, wherein a non-linear least-squares method is used to obtain a fit.

12. The apparatus of claim 11, wherein said lower limit is below 500 sec/mm.sup.2.

13. The apparatus of claim 11, wherein said lower limit is below 300 sec/mm.sup.2.

14. The apparatus of claim 11, wherein said lower limit is below 150 sec/mm.sup.2.

15. The apparatus of claim 11, wherein said upper limit is above 2000 sec/mm.sup.2.

16. The apparatus of claim 11, wherein said upper limit is above 3000 sec/mm.sup.2.

17. The apparatus of claim 11, wherein said upper limit is about 5000 sec/mm.sup.2 or higher.

18. A method to infer numerically an estimation of a mean value of absolute noise residuals in a Rician signal distribution in magnitude signal data obtained with magnetic resonance imaging, said method comprising the steps of: a) simulating a multitude of signal levels θ≥0; b) simulating for each signal level a multitude of realizations with the addition of Gaussian distributed random noise with σ.sub.g=1; c) forming the magnitude of each signal so realized in step b); d) determining for each signal level θ the absolute value of the difference between original signal obtained in step a) and each value obtained in the previous step; e) storing the arithmetic average of these differences obtained in step d) for each signal level θ as function values of α(θ); f) determining suitable linear interpolations for α(θ).

19. The method according to claim to 18, wherein after step 11c), a step forming arithmetic averages within groups of k such magnitude signals is inserted.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0149] These and other objects and advantages of the present invention will be better understood by those skilled in the art in view of the following detailed description of a preferred embodiment of the invention rendered in conjunction with the appended drawings. In appended drawings, like reference numerals are utilized to refer to like elements throughout, and:

[0150] FIG. 1 is an illustrative graphical representation of the logarithmic decay signal intensity from pure water at 17.5 C as a function of b-factor;

[0151] FIG. 2 is an illustrative graphical representation of the logarithmic decay signal intensity in tissue as a function of b-factor, with mono-exponential and non-mono-exponential model fits;

[0152] FIG. 3 is an illustrative graphical representation of the logarithmic decay signal intensity in normal prostate and prostate tumor tissue as a function of b-factor, with respective mono-exponential fits for b-factors up to 1000 sec/mm.sup.2;

[0153] FIGS. 4A and 4B are illustrative graphical representations of Rician signal probability density distributions and their expectation values (FIG. 4A), as well as resulting bias (FIG. 4B) at low signal to noise ratios;

[0154] FIGS. 5A and 5B are illustrative graphical representations of ADC maps of prostate cancer obtained with maximum b-factors of 500 sec/mm.sup.2 (FIG. 5A) and 1400 sec/mm.sup.2 (FIG. 5B);

[0155] FIG. 6 is an illustrative graphical representation of the ADC measured with b-factors 0 and 1000 sec/mm.sup.2 at different locations along the magnet principal axis in a water phantom surrounded with ice-water and a known ADC value of 1.1 μm.sup.2/ms at OC temperature;

[0156] FIG. 7 is a high-level block diagram of an illustrative embodiment of a magnetic resonance imaging system suitable for use in the method of present invention;

[0157] FIG. 8 is an imaging pulse sequence diagram for the acquisition of one image of data from a slice using the system illustratively depicted in FIG. 7;

[0158] FIGS. 9A and 9B are illustrative graphical representations of the correction function a without averaging (FIG. 9A) and with averaging (FIG. 9B) determined with Monte Carlo simulations for signal-to-noise ratios between 0 and 6;

[0159] FIG. 10 is an illustrative graphical representation of simulated data for a b-factor range of 0 to 3000 sec/mm.sup.2 and the result of using the method described herein to infer bias-free signal decay and bias;

[0160] FIGS. 11A and 11B are illustrative graphical representations of geometrically averaged raw data of three diffusion-encoding directions obtained with the preferred imaging protocol at b-factors 1500 sec/mm.sup.2 (FIG. 11A) and 3000 sec/mm.sup.2 (FIG. 11B);

[0161] FIGS. 12A-D are illustrative comparisons of diffusion encoding direction averaged σ.sub.g maps obtained with the method described herein using the bi-exponential model (FIG. 12A), the gamma distribution related model (FIG. 12B) and the kurtosis model without (FIG. 12C) and with (FIG. 12D) spatial low pass filtering;

[0162] FIGS. 13A-D are illustrative comparisons of diffusion-weighted images at a b-factor of 1500 sec/mm.sup.2 obtained with conventional averaging (FIG. 13A) and processed with the novel method using the bi-exponential model (FIG. 13B), the gamma distribution related model (FIG. 13C) and the kurtosis model (FIG. 13D);

[0163] FIGS. 14A-D are illustrative comparisons of diffusion-weighted images at a b-factor of 3000 sec/mm.sup.2 obtained with the novel method for a single diffusion-encoding direction with bias (FIG. 14A) and without bias (FIG. 14B), as well for three orthogonal diffusion-encoding directions geometrically averaged with bias (FIG. 14C) and without bias (FIG. 14D);

[0164] FIG. 15 is an illustrative graphical representation of an array of 20 images reconstructed with the new method proposed herein with b-factors equally spaced between 150 and 3000 sec/mm.sup.2.

[0165] FIG. 16 is an illustrative graphical representation of a simulation to estimate the averaging effect of different fitting approaches.

[0166] FIG. 17 is an illustrative graphical representation of a computed signal-to-noise ratio map together with a scale of absolute signal-to-noise ratio values for a single diffusion-encoding direction and a b-factor of 0.

[0167] FIGS. 18A and 18B are illustrative graphical representations of kurtosis parameters K vs. ADC.sub.K for region of interest data measured in prostate (FIG. 18A) and pixel data measured in prostate and normal brain (FIG. 18B);

[0168] FIG. 19 is an illustrative graphical representation of kurtosis parameters K.sub.2 vs. K.sub.1 for pixel data measured in normal brain;

[0169] FIG. 20A-C are illustrative graphical representations of ADC.sub.K maps obtained with different ranges of b-factors (FIG. 20A) and diffusion-weighted images that were either measured, extrapolated with the proposed kurtosis model with assumed correlation between K and ADC.sub.K, or extrapolated with a basic mono-exponential model for b=1633 sec/mm2 (FIG. 20B) and b=2567 sec/mm2 (FIG. 20C).

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

[0170] As alluded to above, the present invention arises primarily from the observation that the tissue water diffusion signal strength vs. b-factor encoding level follows a non-mono-exponential function and that with signal-to-noise ratios prevalent on clinical imaging systems, this signal decay along an arbitrary diffusion encoding direction, but with b-factors higher than 100 to 300 sec/mm.sup.2, can adequately be described with basic non-mono-exponential model functions that use four parameters or less. Specifically, diffusion experiments in tissues at different signal to noise ratios suggest that model-related differences are exceedingly small and, therefore, for typical signal-to-noise any differences between model fit and measured signals are primarily noise related.

[0171] If this observation is true, there is significant potential for the use of such model functions to infer noise in diffusion-weighted magnetic resonance imaging experiments. This in turn, permits the determination and elimination of signal bias that is introduced by the conversion of the complex magnetic resonance signal to a magnitude signal. The method, furthermore, very efficiently can obtain bias-free diffusion-weighted images with excellent signal-to-noise ratio for arbitrary diffusion weighting levels for optimal viewing of lesions, whereas conventional acquisition protocols rely on signal averaging at selected b-factors to obtain bias-contaminated diffusion-weighted images with acceptable signal-to-noise ratio. Moreover, the method provides accurate so-called “apparent diffusion coefficients” that can be largely independent of imaging protocol design and magnetic field gradient non-linearities. The ramifications of this are significant as will be discussed in greater detail in the example below.

[0172] The various steps of the method of the invention now will be discussed. These steps include: [0173] (1) simulating a multitude of signal levels θ and simulating for each signal level a multitude of realizations with the addition of Gaussian distributed random noise with σ.sub.g=1. Determining for each signal level and each realization the absolute difference between original signal θ and the magnitude of the signals with noise added. Storing the arithmetic average of these differences for each signal level θ as function values of α(θ); [0174] (2) alternatively, simulating a multitude of signal levels θ and simulating for each signal level a multitude of realizations with the addition of Gaussian distributed random noise with σ.sub.g=1. Forming the magnitude of each signal so realized and generating arithmetic averages within groups of k such magnitude signals. Determining for each signal level θ the absolute difference between original signal and each averaged magnitude signal with incorporated noise. Storing the arithmetic average of these differences for each signal level θ as function values of α(θ); [0175] (3) determining suitable linear interpolations for α(θ); [0176] (4) the provision of an appropriate magnetic resonance imaging apparatus; [0177] (5) acquiring a plurality of sets of signals from a patient using the magnetic resonance imaging apparatus at each of a selected number of encoding levels distributed across a wide range of b-factors along at least one diffusion encoding direction; [0178] (6) setting or initializing the value of the noise σ.sub.g for each pixel to 0 or a suitable very small value; [0179] (7) processing the acquired or bias-corrected signal decays on a pixel by pixel basis to obtain the best possible fit between them and a non-mono-exponential equation of the form:

S(b)=S.sub.0(ƒ.Math.exp(−bD.sub.f)+(1−f).Math.exp(−bD.sub.s))  In this equation, S.sub.0 is a constant that depends on constant values associated with (a) the magnetic imaging apparatus, (b) the spin-spin relaxation time T2, (c) the spin lattice relaxation time T1, and (d) the spin density ρ. Also, b is the diffusion encoding level, preferably determined for each pixel location, D.sub.ƒ and D.sub.s describe the diffusion coefficients of two compartments with fast and slow diffusion, respectively, and the fraction ƒ is used to account for the relative signal contribution of each compartment. Preferably, the best fit is obtained with an advanced non-linear least-square method such as the Levenberg-Marquardt algorithm; [0180] (8) alternatively, processing the acquired or bias-corrected signal decays on a pixel by pixel basis to obtain the best possible fit between them and a non-mono-exponential equation of the form:

[00014]$S\left(b\right)=\frac{S_0}{{\left(1+\theta b\right)}^k}$  In this equation, S.sub.0 is a constant that depends on constant values associated with (a) the magnetic imaging apparatus, (b) the spin-spin relaxation time T2, (c) the spin lattice relaxation time T1, and (d) the spin density ρ. Also, b is the diffusion encoding level, preferably determined for each pixel location, and k and θ are a shape and scaling parameter, respectively, of a gamma distribution that is associated with this function. Preferably, the best fit is obtained with an advanced non-linear least-square method such as the Levenberg-Marquardt algorithm; [0181] (9) alternatively, processing the acquired or bias-corrected signal decays on a pixel by pixel basis to obtain the best possible fit between them and a non-mono-exponential equation of the form:

[00015]$S\left(b\right)=S_0\exp \left(-b.Math.\mathrm{AD}{C}_K+b^2.Math.\mathrm{AD}{C}_K^2.Math.\frac{K}{6}\right)$  In this equation, S.sub.0 is a constant that depends on constant values associated with (a) the magnetic imaging apparatus, (b) the spin-spin relaxation time T2, (c) the spin lattice relaxation time T1, and (d) the spin density ρ. Also, b is the diffusion encoding level, preferably determined for each pixel location, K a parameter termed excess kurtosis and ADC.sub.K the apparent kurtosis diffusion coefficient. Preferably, the best fit is obtained with an advanced non-linear least-square method such as the Levenberg-Marquardt algorithm; [0182] (10) alternatively, processing the acquired or bias-corrected signal decays on a pixel by pixel basis to obtain the best possible fit between them and a non-mono-exponential equation of the form:

S(b)=S.sub.0 exp(−(bD).sup.α)

In this equation, S.sub.0 is a constant that depends on constant values associated with (a) the magnetic imaging apparatus, (b) the spin-spin relaxation time T2, (c) the spin lattice relaxation time T1, and (d) the spin density ρ. Also, b is the diffusion encoding level, preferably determined for each pixel location, and D and a parameters that determine the shape of the decay curve. This model is recommended to account for perfusion effects and, therefore, a measurement at a very low or zero b-factor should be included. Preferably, the best fit is obtained with an advanced non-linear least-square method such as the Levenberg-Marquardt algorithm; [0183] (11) determining θ.sub.i=ŝ.sub.i/σ.sub.g for each fitted signal ŝ.sub.i as an estimate for the signal-to-noise ratio at each measured signal s.sub.i; [0184] (12) obtaining for each pixel a first estimation of noise σ.sub.g by performing the following computation:

[00016]${\sigma }_g=\frac{1}{N_f}\underset{1}{\overset{N}{.Math.}}\frac{{\left({\hat{s}}_i-s_i\right)}^2}{\begin{array}{c}2\frac{{\hat{s}}_i}{{\theta }_i}+2{\theta }_i({\hat{s}}_i-\\ \frac{{\theta }_i^2}{2{\hat{s}}_i^2}\left(\exp \left(-\frac{{\theta }_i^2}{4}\right)\sqrt{\frac{\pi }{2}}\frac{{\hat{s}}_i}{{\theta }_i}\left[\left({\hat{s}}_i^2+2\frac{{\hat{s}}_i^2}{{\theta }_i^2}\right)I_0\left(\frac{{\theta }_i^2}{4}\right)+{\hat{s}}_i^2{I}_1\left(\frac{{\theta }_i^2}{4}\right)\right]\right))\end{array}}$  where N equals the number of b-factors to be analyzed, N.sub.ƒ=N−N.sub.p, and N.sub.p equals the number of free fit parameters. Moreover, wherein α was determined according to step (1) or, if averaging is performed, according to step (2), and ŝ.sub.i is the best fit equation result and s.sub.i is the measured signal; [0185] (13) alternatively, obtaining for each pixel a first estimation of noise σ.sub.g by performing the following computation:

[00017]${\sigma }_g=\frac{1}{N_f}\underset{1}{\overset{N}{.Math.}}.Math.{\hat{s}}_i-s_i.Math./\alpha \left({\theta }_i\right)$ [0186] where N equals the number of b-factors to be analyzed, N.sub.ƒ=N−N.sub.p, and N.sub.p equals the number of free fit parameters. Moreover, wherein α was determined according to step (1) or, if averaging is performed, according to step (2), and ŝ.sub.i is the best fit equation result and s.sub.i is the measured signal; [0187] (14) determining for each pixel the Rician signal bias b.sub.i for each fitted signal level ŝ.sub.i according the equation (see, Koay C G and Basser P J.; Analytically Exact Correction Scheme for Signal Extraction from Noisy Magnitude MR Signals, Journal of Magn Reson, 179:317-322, 2006):

[00018]$b_i=\frac{1}{2{\sigma }_g^2}\left(\exp \left(-\frac{{\hat{s}}_i^2}{4{\sigma }_g^2}\right)\sqrt{\frac{\pi }{2}}{\sigma }_g\left[\left({\hat{s}}_i^2+2{\sigma }_g^2\right)I_0\left(\frac{{\hat{s}}_i^2}{4{\sigma }_g^2}\right)+{\hat{s}}_i^2{I}_1\left(\frac{{\hat{s}}_i^2}{4{\sigma }_g^2}\right)\right]\right)-{\hat{s}}_i$ [0188] (15) obtaining for each pixel a bias-corrected signal value s.sub.i.sup.c by subtracting the Rician signal bias b.sub.i from the measured signal level s.sub.i; [0189] (16) repeating steps (7) to (15), but without changing the fitting function, until the change of s.sub.i.sup.c falls under a predefined threshold; [0190] (17) averaging arithmetically on a pixel by pixel basis the σ.sub.g values determined along different diffusion encoding directions; [0191] (18) performing a spatial averaging of the directionally averaged σ.sub.g values; [0192] (19) dividing on a pixel by pixel basis s.sub.i measured at the lowest b-factor with σ.sub.g as processed in steps (17) and (18) and to store the resulting values as signal-to-noise ratio map; [0193] (20) repeating steps (14) and (15) with σ.sub.g as processed in steps (17) and (18) to obtain for each pixel a final set of bias-corrected signal values s.sub.i.sup.c; [0194] (21) processing bias-corrected signal decays on a pixel by pixel basis to obtain the best possible final fit with any non-mono-exponential function, preferably one with low number of free parameters such as an equation of the form:

[00019]$S\left(b\right)=S_0\exp \left(-b.Math.\mathrm{AD}{C}_K+b^2.Math.\mathrm{AD}{C}_K^2.Math.\frac{K_1}{6}-b^3.Math.\mathrm{AD}{C}_K^3.Math.\frac{K_2^2}{27}\right)$ [0195]  In this equation, S.sub.0 is a constant that depends on constant values associated with (a) the magnetic imaging apparatus, (b) the spin-spin relaxation time T2, (c) the spin lattice relaxation time T1, and (d) the spin density ρ. Also, b is the diffusion encoding level, preferably determined for each pixel location, ADC.sub.K the apparent kurtosis diffusion coefficient, and K.sub.1 and K.sub.2 excess kurtosis parameters that are expressed as functions of ADC.sub.K. Preferably, the best fit is obtained with an advanced non-linear least-square method such as the Levenberg-Marquardt algorithm; [0196] (22) displaying arithmetic averages of resulting direction-dependent parameter maps and geometric averages of the direction-dependent fitted signals for each b-factor sampled, or any other b-factor level inside or outside the sampled range of b-factors.

[0197] Some of the data presented herein was obtained with line scan diffusion imaging. This technique and an apparatus suitable for practice thereof is described in detail in U.S. Pat. No. 5,786,692 issued on Jul. 28, 1998 to Stephan E. Maier, et al. and entitled “Line Scan Diffusion Imaging”, which is hereby incorporated by reference into this specification. The suitability of this technique for diffusion imaging at high b-factors is documented in detail in U.S. Pat. No. 6,751,495 issued on Jun. 15, 2004 to Stephan E. Maier et al. and entitled “Method of Fast and Reliable Tissue Differentiation Using Diffusion-Weighted Magnetic Resonance Imaging”, which is hereby incorporated by reference into this specification.

[0198] The preferred embodiment of the invention herein described illustratively adopts the single-shot echo planar imaging technique (see, Turner R, Le Bihan D, Maier J, Vavrek R, Hedges L K, and Pekar J.; Echo-Planar Imaging of Intravoxel Incoherent Motions, Radiology, 177:407-414, 1990). This imaging method appears to be suitable for rapidly obtaining largely motion artifact free images even at very high b-factors of up to 5000 sec/mm.sup.2. An apparatus and pulse sequence suitable for the practice thereof, as will be seen from FIGS. 7 and 8, includes a magnetic resonance imaging system 10, generally having a magnet assembly, interface circuitry, and a computer 40. The magnet assembly includes a very strong magnet 13 that creates a homogeneous magnetic field within and around a sample (e.g. an inert sample or patient). X, Y, and Z magnetic field gradient coils 14, 16, and 18 also form a part of the assembly and are positioned proximate or surrounding the sample 20. The assembly further comprises one or more RF coils 22, which are positioned near or around the sample.

[0199] The interface circuitry includes a gradient waveform generator 24 that has signal outputs connected to the gradient coils 14, 16, and 18, and a control input connected to the computer and an output connected to an input of an RF power amplifier 28. The RF power amplifier has an output connected to an input of an RF switch 30. The RF power amplifier has an output connected to the RF coil 22, and has an output connected to the input of an RF detector 32.

[0200] The computer 40 includes computing hardware 42 and storage. The computing hardware can comprise special purpose hard-wired computing circuitry dedicated to magnetic resonance acquisition, an imaging and a special programmed general purpose computer, or a combination of both. The storage 46 can include various types of storage, such as disk storage and random access memory.

[0201] The storage can be used to store data and programs, including programs used to interact with the system's interface circuitry 12. The computer has a video output for providing video signals to a display 46, as well as control outputs connected respectively to control inputs of the gradient waveform generator 24 and the RF signal generator 26. The computer also has acquisition input operatively connected to an output of the RF detector 32.

[0202] In operation, referring to FIGS. 7 and 8, the imaging system 10 builds images sequentially under control of the computer 40 according to a single-shot echo-planar imaging (EPI) sequence. At the beginning of an acquisition sequence for an image, the computer 40 sends a signal to the RF signal generator 26, which responds by generating a water-selective 900 pulse 50. This pulse is amplified by the RF power amplifier 28 and provided to the RF coil 22 via the RF switch 30. As this pulse is being provided, the computer instructs the gradient waveform generator 24 to drive the Z coil 16 with a slice and frequency-selective multi-polar pulse.

[0203] Next, the gradient waveform generator 24 provides a first set of diffusion gradient pulses 54, 56, and 58 respectively to the X, Y, and Z gradient coils 14, 16, and 18. These gradient signals each include a signal rectangular pulse, which is provided in order to sensitize the magnetic resonance imaging process to diffusion. After the falling edge of the diffusion gradient signals, a 180° pulse 60 is provided to the RF coil 22, in much the same way that the 90° pulse was. At the same time, the gradient waveform generator provides a rectangular pulse 62 on the Z gradient coil. This pulse is of shorter duration than the diffusion gradient pulses. Then, the waveform generator provides a second set of diffusion gradient signals 66, 68, and 70 respectively to the X, Y, and Z gradient coils 14, 16 and 18. These second diffusion gradient signals are of the same amplitude and duration as the first diffusion gradient signals. Once the second diffusion gradient signals are turned off, the gradient waveform generator provides repeated pulses 72 on the Y coil 16 for spatial encoding along the phase direction. At the same time, the gradient waveform generator provides repeated pulses 74 on the X coil 14 for spatial encoding along the frequency direction.

[0204] As a result of this excitation sequence, a train of echoes 76 is received from the slice that was excited by the 90° and 180° pulses. The RF coils receives these echoes and provides them via RF switch 30 to the RF detector. The computer 40 receives the output of the detector, and processes it to obtain one image to be displayed on the display 46. After the echo has been received, optional crusher gradient signals 78, 80, and 82 can be applied to the gradient coils 14, 16 and 18.

[0205] The method then proceeds with the acquisition of images from a patient using the magnetic resonance imaging apparatus at each of a selected plurality of encoding directions and of encoding levels distributed across the range of b-factors within its capability. This is accomplished by varying the gradient amplitude G along an arbitrary direction in the equation for the b-factor b=γ.sup.2G.sup.2δ.sup.2(Δ−δ/3). Each of the acquired sets of signal decays is representative of an image of a selected anatomical cross-section of the patient. Further, each individual signal decay in each set corresponds to a pixel of its associated image.

[0206] In the preferred embodiment α(θ) is determined without averaged signals.

[0207] This involves a) simulating a multitude of signal levels θ; b) simulating for each signal level a multitude of realizations with the addition of Gaussian distributed random noise with σ.sub.g=1; c) determining for each signal level and each realization the absolute difference between original signal θ and the magnitude of the signals with noise added; d) storing the arithmetic average of these differences for each signal level θ as function values of α(θ) and e) determining suitable linear interpolations for α(θ). Next, data processing begins with a first setting of noise σ.sub.g for each pixel to 0 or a suitable very small value.

[0208] In the preferred embodiment the fitting model function is a bi-exponential decay function. The acquired or bias-corrected signal decays are processed on a pixel by pixel basis to obtain the best possible fit between them and a non-mono-exponential equation of the form:

S(b)=S.sub.0(ƒ.Math.exp(−bD.sub.f)+(1−f).Math.exp(−bD.sub.s))

In this equation, S.sub.0 is a constant that depends on constant values associated with (a) the magnetic imaging apparatus, (b) the spin-spin relaxation time T2, (c) the spin lattice relaxation time T1, and (d) the spin density ρ. Also, b is the diffusion encoding level, preferably determined for each pixel location, D.sub.ƒ and D.sub.s describe the diffusion coefficients of two compartments with fast and slow diffusion, respectively, and the fraction ƒ is used to account for the relative signal contribution of each compartment. Preferably, the best fit is obtained with an advanced non-linear least-square method such as the Levenberg-Marquardt algorithm.

[0209] After completion of this step, θ.sub.i=ŝ.sub.i/σ.sub.g is determined for each fitted signal s.sub.i as an estimate for the signal-to-noise ratio at each measured signal s.sub.i. Then for each pixel a new estimate of noise σ.sub.g is obtained by performing the following computation:

[00020]${\sigma }_g=\frac{1}{N_f}\underset{1}{\overset{N}{.Math.}}\frac{{\left({\hat{s}}_i-s_i\right)}^2}{\begin{array}{c}2\frac{{\hat{s}}_i}{{\theta }_i}+2{\theta }_i({\hat{s}}_i-\\ \frac{{\theta }_i^2}{2{\hat{s}}_i^2}\left(\exp \left(-\frac{{\theta }_i^2}{4}\right)\sqrt{\frac{\pi }{2}}\frac{{\hat{s}}_i}{{\theta }_i}\left[\left({\hat{s}}_i^2+2\frac{{\hat{s}}_i^2}{{\theta }_i^2}\right)I_0\left(\frac{{\theta }_i^2}{4}\right)+{\hat{s}}_i^2{I}_1\left(\frac{{\theta }_i^2}{4}\right)\right]\right))\end{array}}$ [0210] where N equals the number of b-factors to be analyzed, N.sub.ƒ=N−N.sub.p, and N.sub.p equals the number of free fit parameters. Moreover, where ŝ.sub.i is the best fit equation result and s.sub.i is the measured signal.

[0211] Alternatively, for each pixel a new estimate of noise σ.sub.g is obtained by performing the following computation:

[00021]${\sigma }_g=\frac{1}{N_f}\underset{1}{\overset{N}{.Math.}}.Math.{\hat{s}}_i-s_i.Math./\alpha \left({\theta }_i\right)$

Where N equals the number of b-factors to be analyzed, N.sub.ƒ=N−N.sub.p, and N.sub.p equals the number of free fit parameters. Moreover, where ŝ.sub.i is the best fit equation result and s.sub.i is the measured signal. Subsequently, for each pixel Rician signal bias b.sub.i is computed for each fitted signal level s.sub.i according the equation:

[00022]$b_i=\frac{1}{2{\sigma }_g^2}\left(\exp \left(-\frac{{\hat{s}}_i^2}{4{\sigma }_g^2}\right)\sqrt{\frac{\pi }{2}}{\sigma }_g\left[\left({\hat{s}}_i^2+2{\sigma }_g^2\right)I_0\left(\frac{{\hat{s}}_i^2}{4{\sigma }_g^2}\right)+{\hat{s}}_i^2{I}_1\left(\frac{{\hat{s}}_i^2}{4{\sigma }_g^2}\right)\right]\right)-{\hat{s}}_i$

Finally, for each pixel a bias-corrected signal value s.sub.i.sup.c is obtained by subtracting the Rician signal bias b.sub.i from the measured signal level s.sub.i. Fitting, signal-to-noise ratio calculation, noise estimation, bias calculation, and bias subtraction from the measured signal are repeated until the change of s.sub.i.sup.c falls under a predefined threshold.

[0212] In the preferred embodiment, σ.sub.g measured along the different diffusion encoding directions is arithmetically averaged and the resulting maps are low pass filtered. The σ.sub.g so obtained is used to perform a final bias calculation and bias subtraction from the measured signal. The bias-corrected signal value s.sub.i.sup.c is then again fitted with a non-mono-exponential function that uses as few parameters as possible. With the resulting fit parameters and the fit function, images of b-factors outside acquired range are computed. For each pixel, the signal value s.sub.i measured at the lowest b-factor is divided with the pixel specific σ.sub.g and the resulting data is stored as signal-to-noise ratio map. As last step, arithmetic averages of resulting direction-dependent parameter maps and geometric averages of the direction-dependent fitted signals for each b-factor sampled, or any other b-factor level inside or outside the sampled range of b-factors are displayed on the computer monitor 46.

[0213] The following example will further describe the invention in terms of the actual experimental process, which led to its creation.

Example

[0214] With wide b-factor range diffusion scans it has been demonstrated that in prostate and human brain water signal attenuation is better described with non-mono-exponential models than a mono-exponential model as it is generally used in clinical applications. The data presented herein is the result of processing diffusion-weighted image data obtained in patients and normal subjects. Diffusion-weighted images of the prostate were obtained with the single-shot echo planar imaging technique on two different 3 Tesla magnetic resonance imaging systems. Diffusion-weighted images of the brain were obtained with line scan diffusion imaging on a 1.5 Tesla magnetic resonance imaging system.

[0215] For the prostate scans in patients, one setup included a 3 Tesla Philips Achieva (Philips, Eindhoven, NL) scanner with the capability to deliver 40 mT/m maximum gradient strength along each axis. Conventional diagnostic imaging with an external RF coil array and two-fold parallel coil acceleration included a diffusion scan of 4 min duration, where images at b-factors 0, 500, 1000 and 1500 sec/mm.sup.2 at isotropic spatial resolution of 3×3×3 mm.sup.3 were obtained. For each of these b-factors, measurements were performed along three orthogonal diffusion encoding directions. Moreover, each measurement was repeated six times and resulting magnitude signals were averaged. The same spatial resolution and coil setup with two-fold acceleration was employed for a wide range b-factor scan that was performed to obtain data for the proposed processing. For this scan of 3 min duration, 21 evenly spaced b-factor measurements between 0 and 3000 sec/mm.sup.2 were performed along three orthogonal diffusion encoding directions.

[0216] A second setup employed for prostate scans in patients, included a 3 Tesla GE Discovery MR750 (General Electric Medical Systems, Milwaukee, Wis.) system with the same gradient capability, but with a single-channel endo-rectal coil (Medrad, Pittsburgh, Pa.). In a 3 min scan, diffusion-weighted images were obtained with an in-plane resolution of 4.4×4.4 mm.sup.2 and 5 mm slice thickness for 15 evenly spaced b-factors between 0 and 3500 sec/mm.sup.2 and three orthogonal diffusion encoding directions.

[0217] Scans in brains of normal volunteers were performed on a 1.5 Tesla Horizon Echospeed (GE Medical Systems, Milwaukee, Wis.) system with 22 mT/m maximum gradient strength along each axis. Line scan diffusion images with an in-plane resolution of 3.4×3.4 mm.sup.2 and 6 mm slice thickness were obtained with a conventional quadrature head coil without parallel acceleration. Diffusion encoding was performed for 32 evenly spaced b-factors between 5 and 5000 sec/mm.sup.2 and along 6 non-collinear and non-coplanar directions. The scan time for a single mid-ventricular axial slice was 27 min.

[0218] As outlined in the summary of the invention, in a preparatory processing step, a as a function of signal-to-noise ratio has to be determined. FIG. 9A shows the function so obtained without averaging for a signal-to-noise ratio between 0 and 6. FIG. 9B shows the same curve (k=1), but also curves for different number of averages k.

[0219] The processing algorithm was implemented on a personal computer with an Intel 4-core i7-6700 CPU 3.40 GHz chip and a Gentoo Linux operating system with 4.14.65 kernel. The final version was implemented with Python version 3.6.5 software. Non-linear least-square fitting was performed with the predefined algorithm “scipy.optimize.curve_fit”, version 0.19.1. For each non-mono-exponential function, suitable parameter start values and value bounds were defined. Typical processing time with this setup was around 30 min per slice and typical number of iterations required equaled three.

[0220] FIG. 10 shows simulated signal decay with a signal-to-noise ratio of 10 at the lowest assumed b-factor. The magnitude of simulated complex signals at 13 b-factors between 0 and 3000 s/mm.sup.2 was processed with the proposed method to estimate bias and true signal decay. As can be appreciated in the graph, true signal and estimated true signal decay match closely, supporting the notion that the proposed process is effective.

[0221] FIGS. 11A and 11B show directionally averaged signals of a prostate scan with an external coil array at b-factors 1500 and 3000 sec/mm.sup.2, respectively. The overall inferior signal to noise ratio, particularly at a b-factor of 3000 sec/mm.sup.2 can readily be appreciated. In the center of the image, a normal prostate, without any obvious high intensity cancer lesion, is clearly visible. The same data used to generate FIGS. 11A and 11B was used to evaluate the new method. To avoid any perfusion effects, which could be detrimental to finding the true tissue diffusion signal decay, only b-factor signals down to the second-lowest level, i.e., 150 s/mm.sup.2, were included. FIGS. 12A-C show resulting noise maps after averaging the noise measured along three diffusion direction directions. Such maps were obtained with the bi-exponential model function (FIG. 12A), the gamma distribution related model function (FIG. 12B), and the kurtosis model function truncated to the second order polynomial (FIG. 12C). The similarity of these maps attests to the assumption that observed residuals are primarily noise related, rather than model function related. Finally, FIG. 12D shows a low pass filtered version of FIG. 12C. The low spatial variability of these estimated noise maps are in good agreement with noise profiles that result with coil arrays and parallel acceleration, i.e, there is clear tendency for noise amplification in the central area that is farthest away from all coils (see Pruessmann K P, Weiger M, Scheidegger M B, and Boesiger P.; SENSE: Sensitivity Encoding for Fast MRI, Magn Reson Med, 42:952-62, 1999).

[0222] FIG. 13A shows the b-factor 1500 sec/mm.sup.2 image obtained with the conventional scan protocol, i.e., after magnitude averaging signals from three directions and six repetitions of each direction. FIGS. 13B-C shows images of the geometric average of signals along three diffusion encoding directions obtained with the proposed distributed b-factor sampling and the new processing method for a b-factor of 1500 sec/mm.sup.2 using the bi-exponential model function (FIG. 13B), the gamma distribution related model function (FIG. 13C), and the kurtosis model function truncated to the second order polynomial (FIG. 13D). The high detail visible on these images, particularly when compared to the image obtained with conventional imaging (FIG. 13A), underline the potential importance of the new method. The quality of these images in comparison to the image obtained with conventional averaging is particularly noteworthy in view of the substantially shorter duration of the experimental scan protocol.

[0223] FIGS. 14A-D show additional images of the same reconstruction as shown in FIGS. 13B-C, but at the highest b-factor measured, i.e., 3000 sec/mm.sup.2. Images shown in FIGS. 14A and B were reconstructed from a single superior-inferior diffusion encoding direction without removal of bias (FIG. 14A) and with removal of bias (FIG. 14B). The central region of the prostate shows clear signs of diffusion anisotropy related signal loss. Without bias removal this effect is clearly diminished. Images shown in FIGS. 14C and D were generated from the geometric average of the reconstructed signals of three diffusion encoding directions without removal of bias (FIG. 14C) and with removal of bias (FIG. 14D). At a diffusion-weighting level of 3000 sec/mm.sup.2, a cancer lesion, if present, would impress as an easy distinguishable spot of high signal intensity. High quality images for the entire range of diffusion-weightings employed in the reconstruction with the new method, i.e., from 150 to 3000 sec/mm.sup.2, sorted in ascending order from left to right and from top to bottom, are shown in FIG. 15.

[0224] To obtain the data shown in FIG. 16, a bi-exponential decay was simulated for normal tissue (ƒ=0.8, D.sub.ƒ=2.2 μm.sup.2/ms, D.sub.s=0.4 μm.sup.2/ms,) and tumor tissue (ƒ=0.6, D.sub.ƒ=2.2 μm.sup.2/ms, D.sub.s=0.2 μm.sup.2/ms). Noise was added to the decay signal in 10000 realizations and for each realization a gamma distribution related decay model, a kurtosis model, a correlated kurtosis model as proposed in this application and a bi-exponential model was fitted to the data. The standard deviation of the differences between fitted signals and the noise-free signal relative to noise was used to estimate the averaging effect. The lowest averaging effect is seen for low b-factors. The averaging effect, as might be expected, increases with decreasing number of model parameters.

[0225] Finally, FIG. 17 shows a synthesized map of signal-to-noise ratio for a single-direction measurement. With the help of the scale that shows absolute signal-to-noise ratio values, it can be concluded that for this experimental setup signal-to-noise ratios in the prostate reach levels between 20 and 30, and in some areas of the prostate peripheral zone as low as 15 or lower.

[0226] The following figures are based on data with exquisitely high signal-to-noise ratios of 100:1 and higher. These high signal-noise-ratios are the result of using large voxels and in the case of the prostate data, a small coil next to the tissue of interest. These high signal-to-noise ratios are believed to be helpful to best reflect the observed correlation between ADC.sub.K and K, K.sub.1 and K.sub.2 respectively. FIG. 18A shows a scatter plot of K vs. ADC.sub.K measured in 55 patients for regions of interest of normal prostate and prostate cancer lesions. For this data a high correlation coefficient r.sup.2=0.92 can be determined (see Langkilde F, Kobus T, Fedorov A, Dunne R, Tempany C, Mulkern R V, and Maier S E.; Evaluation of Fitting Models for Prostate Tissue Characterization Using Extended-Range b-Factor Diffusion-Weighted Imaging, Magn Reson Med. 79:2346-2358, 2018). Meanwhile, FIG. 18B shows the same data, but now together with prostate pixel data of an individual patient and brain data of one normal subject. Except for pixels that can be attributed to gray matter, the parameters K vs. ADC.sub.K appear to exhibit a strong correlation. Finally, to obtain the data shown in FIG. 19, brain diffusion signals were processed with a kurtosis fit truncated to the third order polynomial. Evidently, there is a significant correlation between K.sub.1 and K.sub.2.

[0227] FIGS. 20A-C demonstrate the application of the novel kurtosis model fitting method based on a predetermined correlation between ADC.sub.K and K. Subsets of signal decay data were generated with a maximum b-factor of 700 sec/mm.sup.2, 1633 sec/mm.sup.2, and 2567 sec/mm.sup.2 with 4, 8, and 12 b-factors, respectively. Based on the data shown in FIGS. 18A and B, a suitable second order polynomial function was determined to describe the parameter K as a function of the independent variable ADC.sub.K. FIG. 20A shows nearly identical ADC.sub.K maps generated with the three data subsets using kurtosis fits, where K follows the predetermined polynomial function of ADC.sub.K. In FIGS. 20B and C, true measured signals at b-factors 1633 sec/mm.sup.2 (FIG. 20B) and 2567 sec/mm.sup.2 (FIG. 20C) are compared to signals extrapolated from the data set with a maximum b-factor of 700 sec/mm.sup.2 using the novel kurtosis fitting method and mono-exponential fitting. The images shown confirm the well-known inadequacy of mono-exponential analysis for extrapolating and describing diffusion at higher b-factors and the high potential the proposed method has to achieve exactly this goal.

[0228] Having thus described a preferred embodiment of the invention and an example indicative of its significance in the art, numerous modifications, alterations, variations, changes and the like will occur to those skilled in the art. For example, various other non-mono-exponential functions and close representations may occur to those skilled in the art within the scope of this invention in its broadest aspects. Accordingly, it is to be understood that the foregoing specification is to be understood as being illustrative only, and the scope of the invention is intended to be limited only by the terms of the appended claims.