CROSS-REFERENCE TO RELATED APPLICATION This application claims priority to U.S. Provisional Application No.

61 /983,946, titled "Method for Field Mapping" and filed on April 24, 2014, the entire contents of which is incorporated herein by reference.

BACKGROUND

The present disclosure relates to magnetic resonance (MR) imaging.

More particularly, the present disclosure relates to methods of estimating and mapping the main magnetic field (BO), where BO inhomogeneity is introduced when a body is placed in a MR scanner.

An accurate BO off-resonance map has many applications, such as in shimming, image distortion correction, quantitative susceptibility mapping and fat/water separation. However, accurate estimation of the BO map can be hampered by the additional phase term related to the chemical shift of fat (1 ). To avoid the effect of chemical shift, conventional BO mapping is performed by phase unwrapping of the complex data where signals from fat and water are "in-phase", which is ensured by using a single-peak model for both water and fat, e.g. 3-point Dixon technique (2,3).

Recently, multi-echo gradient echo (GRE) acquisition has attracted the attention of the academic community because of its potential for performing R2 ^* (decay rate) mapping, BO mapping, quantitative fat/water separation and quantitative susceptibility mapping from a single scan. Many multi-point Dixon techniques (1 ,4,5) have been established to quantify fat-fraction (FF) through iterative multi-parameter fitting procedures by using a multi-peak fat model, where the fitting parameters include the frequency shift (f _b o) due to BO inhomogeneity, the signals of the water (p _w ) and fat (PF) components, and R ₂ ^* . The results of the iterative procedure on a voxel-by-voxel basis (6) can be further enhanced using region-growing (7), region-merging (8), region-labeling (9), fat likelihood analysis (1 0), applying a multi-resolution strategy for guiding the search process (1 1 ,12), and algebraic decomposition (1 3). The state of the art is to estimate the aforementioned multi-parameters jointly for all voxels by using the graph-cut algorithm in BO map estimation (4,14-16).

The multi-point Dixon techniques mainly aim to generate an accurate FF map instead of a realistic and accurate BO map. To improve the robustness for preventing fat/water swaps, most of these techniques apply two constraints during the BO searching process: the prior knowledge of spatial smoothness and a limited range of f _b o values. Because the effect of the phase term in the complex signal equation is periodic, an accurate FF measurement can be achieved from either "non-smooth" (9) or "over smooth" (8) BO maps. Further, if the true BO inhomogeneity is outside the range, the resulting BO map will be wrapped and might not be useful for further data processing.

SUMMARY

Various embodiments of the present disclosure present unwrapping- based BO mapping systems and methods. An initial BO map is obtained by unwrapping a pseudo-in-phase complex data that is generated from selected echoes of multi-echo gradient echo image data. Global and local phase errors that are present in the initial BO map are corrected to produce a phase-error- corrected BO map. A bias frequency shift is employed to produce a corrected fat-fraction map based on the phase-error-corrected BO map. The corrected fat-fraction map is employed to generate a final BO map.

Accordingly, in a first aspect, there is provided a method of mapping of the main magnetic field (BO) of a magnetic resonance scanner, the method comprising:

receiving multiple gradient echo image data;

processing the multiple gradient echo image data to generate a magnitude-based water mask;

generating a complex data set based on selected echoes of the multiple gradient echo image data, wherein the selected echoes are selected to have an effective echo delay time that corresponds to a fat-induced phase shift that is closest to a multiple of 2 ^* pi;

unwrapping the phase of the complex data set to obtain an initial BO map;

calculating a bias frequency shift based on a calculated phase shift of fat during acquisition of the selected echoes and the effective echo time; calculating a bias-corrected phase shift of fat that is produced during acquisition of the selected echoes;

employing the magnitude-based water mask to correct for a global phase error present in the initial BO map, thereby obtaining a globally phase- error-corrected BO map;

employing the bias-corrected phase shift of fat to correct for local phase errors in the globally phase-error-corrected BO map, thereby obtaining a phase-error corrected BO map;

calculating a corrected fat-fraction map based on the phase-error corrected BO map and the bias frequency shift; and

employing the corrected fat-fraction map to perform bias removal on the phase-error corrected BO map and generating a final BO frequency map.

In another aspect, there is provided a method of mapping of the main magnetic field (BO) of a magnetic resonance scanner, the method comprising: receiving multiple gradient echo image data;

processing the multiple gradient echo image data to generate a magnitude-based water mask;

generating a complex data set based on selected echoes of the multiple gradient echo image data, wherein the selected echoes are selected to have an effective echo delay time that corresponds to a fat-induced phase shift that is closest to a multiple of 2 ^* pi;

unwrapping the phase of the complex data set to obtain an initial BO map;

employing the magnitude-based water mask to correct for a global phase error present in the initial BO map, thereby obtaining a globally phase- error-corrected BO map;

employing a bias-corrected phase shift of fat to correct for local phase errors in the globally phase-error-corrected BO map, thereby obtaining a phase-error corrected BO map;

calculating a corrected fat-fraction map based on the phase-error corrected BO map and a bias frequency shift; and employing the corrected fat-fraction map to perform bias removal on the phase-error corrected BO map and generating a final BO frequency map.

In another aspect, there is provided a method of mapping of the main magnetic field (BO) of a magnetic resonance scanner, the method comprising: receiving multiple gradient echo image data;

generating a complex data set based on selected echoes of the multiple gradient echo image data, wherein the selected echoes are selected to have an effective echo delay time that corresponds to a fat-induced phase shift that is closest to a multiple of 2 ^* pi;

unwrapping the phase of the complex data set to obtain an initial BO map;

correcting for a global phase error and local phase errors present in the initial BO map, thereby obtaining a phase-error corrected BO map;

calculating a corrected fat-fraction map based on the phase-error corrected BO map and a calculated bias frequency shift; and

employing the corrected fat-fraction map to perform bias removal on the phase-error corrected BO map and generating a final BO frequency map.

In another aspect, there is provided a method of mapping of the main magnetic field (BO) of a magnetic resonance scanner, the method comprising: receiving multiple gradient echo image data;

generating a complex data set based on selected echoes of the multiple gradient echo image data, wherein the selected echoes are selected to have an effective echo delay time that corresponds to a fat-induced phase shift that is closest to a multiple of 2 ^* pi;

unwrapping the phase of the complex data set to obtain an initial BO map;

correcting for a global phase error and local phase errors present in the initial BO map, thereby obtaining a phase-error corrected BO map; and

performing bias removal on the phase-error corrected BO map and generating a final BO frequency map.

In another aspect, there is provided a method of mapping of the main magnetic field (BO) of a magnetic resonance scanner, the method comprising: receiving multiple gradient echo image data;

generating a complex data set based on selected echoes of the multiple gradient echo image data, such that the complex data set is a pseudo-in-phase complex data set;

unwrapping the phase of the complex data set to obtain an initial BO map;

correcting for a global phase error and local phase errors present in the initial BO map, thereby obtaining a phase-error corrected BO map; and

performing bias removal on the phase-error corrected BO map and generating a final BO frequency map.

In another aspect, there is provided a magnetic resonance imaging system comprising:

a magnetic resonance scanner; and

computing hardware operatively connected to said magnetic resonance scanner, said computing hardware comprising a processor coupled to a memory, wherein the processor, in response to executing instructions stored in the memory, is configured to:

process multiple gradient echo image data to generate a magnitude-based water mask;

generate a complex data set based on selected echoes of the multiple gradient echo image data, wherein the selected echoes are selected to have an effective echo delay time that corresponds to a fat- induced phase shift that is closest to a multiple of 2 ^* pi;

unwrap the phase of the complex data set to obtain an initial BO map;

calculate a bias frequency shift based on a calculated phase shift of fat during acquisition of the selected echoes and the effective echo time;

calculate a bias-corrected phase shift of fat that is produced during acquisition of the selected echoes;

employ the magnitude-based water mask to correct for a global phase error present in the initial BO map, thereby obtaining a globally phase-error-corrected BO map;

employ the bias-corrected phase shift of fat to correct for local phase errors in the globally phase-error-corrected BO map, thereby obtaining a phase-error corrected BO map;

calculate a corrected fat-fraction map based on the phase-error corrected BO map and the bias frequency shift; and

employ the corrected fat-fraction map to perform bias removal on the phase-error corrected BO map and generating a final BO frequency map.

A further understanding of the functional and advantageous aspects of the disclosure can be realized by reference to the following detailed description and drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

Embodiments will now be described, by way of example only, with reference to the drawings, in which:

FIG. 1 is a block diagram of an example system for performing various methods of the present disclosure.

FIG. 2A is a flowchart illustrating a method of performing BO mapping according to an example embodiment of the present disclosure.

FIG. 2B is a flowchart illustrating a method of performing BO mapping according to another example embodiment of the present disclosure.

FIG. 3 is a flowchart illustrating an example method of echo selection for the generation of a pseudo-in-phase complex data set.

FIG. 4 is a flowchart illustrating an example method of generating a magnitude-based water map.

FIG. 5 is a flowchart illustrating an example method of performing global phase error correction.

FIG. 6A is a flowchart illustrating an example method of performing local phase error corrections.

FIG. 6B is a flowchart illustrating an example method of performing additional local phase error corrections.

FIG. 6C is a flowchart illustrating another example method of performing additional local phase error corrections.

FIGS. 7A-C illustrates an example of processing a unipolar 6-echo GRE data set (IDEAL-IQ sequence) acquired using a GE 3.0 T scanner. (A) and (B) represent BO-NICE generated BO maps and R2 ^* maps for three consecutive slices (from left to right), respectively; (C) shows a comparison of FF maps between the IDEAL, IDEAL with BO-NICE, and BO-NICE for the three consecutive slices (from top to bottom), respectively. Fat/water swaps are observed in the FF map generated by using GE product software

(indicated by arrows).

FIGS. 8A-C show an example of processing a bipolar 5-echo GRE data set acquired using a Siemens 3.0 T scanner: (A), (B) and (C) are the axial, coronal and sagittal orientation, respectively. From left to right are the magnitude image at the first echo, BO map, Fat-fraction map and R2 ^* map, respectively.

FIGS. 9A-G demonstrate an example of processing the bipolar dual- echo GRE data acquired using a Philips 1 .5 T: (A) magnitude image from an individual channel; (B) phase image from the same channel as (A), where the bipolar-related phase error has been corrected; (C) initial BO map from the same channel as (A), arrow pointing to phase error; (D) error-corrected BO map from the same channel as (A); (E) the channel-combined BO map; (F) and (G) are the final channel-combined water-only and fat-only images, respectively.

DETAILED DESCRIPTION

Various embodiments and aspects of the disclosure will be described with reference to details discussed below. The following description and drawings are illustrative of the disclosure and are not to be construed as limiting the disclosure. Numerous specific details are described to provide a thorough understanding of various embodiments of the present disclosure. However, in certain instances, well-known or conventional details are not described in order to provide a concise discussion of embodiments of the present disclosure.

As used herein, the terms "comprises" and "comprising" are to be construed as being inclusive and open ended, and not exclusive. Specifically, when used in the specification and claims, the terms "comprises" and

"comprising" and variations thereof mean the specified features, steps or components are included. These terms are not to be interpreted to exclude the presence of other features, steps or components.

As used herein, the term "exemplary" means "serving as an example, instance, or illustration," and should not be construed as preferred or advantageous over other configurations disclosed herein.

As used herein, the terms "about" and "approximately" are meant to cover variations that may exist in the upper and lower limits of the ranges of values, such as variations in properties, parameters, and dimensions. Unless otherwise specified, the terms "about" and "approximately" mean plus or minus 25 percent or less.

It is to be understood that unless otherwise specified, any specified range or group is as a shorthand way of referring to each and every member of a range or group individually, as well as each and every possible sub-range or sub -group encompassed therein and similarly with respect to any subranges or sub-groups therein. Unless otherwise specified, the present disclosure relates to and explicitly incorporates each and every specific member and combination of sub-ranges or sub-groups. As used herein, the term "on the order of", when used in conjunction with a quantity or parameter, refers to a range spanning approximately one tenth to ten times the stated quantity or parameter.

Referring now to FIG. 1 , an example magnetic resonance (MR) system is provided for performing methods disclosed herein. The example system includes a magnetic resonance scanner 100 that employs a main magnet 102 to produce a main magnetic field BO, which generates a polarization in a patient 120 or the examined subject. The example system includes gradient coils 1 04 for generating magnetic field gradients. A reception coil 108 detects the MR signals from patient 120. The reception coil 1 08 can also be used as a transmission coil. Alternatively, body coil 106 may be employed to radiate and/or detect radio frequency (RF) pulses. The RF pulses are generated by an RF unit 1 1 0, and the magnetic field gradients are generated by a gradient unit 1 1 2. The manner by which MR signals are detected using the sequence of RF pulses and magnetic field gradients, and how MR images are reconstructed in general, are known to those skilled in the art.

It will be understood that the MR system can have additional units or components that are not shown for clarity, such as, but not limited to, additional control or input devices, and additional sensing devices, such as devices for cardiac and/or respiratory gating. Furthermore, the various units can be realized other than in the depicted separation of the individual units. It is possible that the different components are assembled into units or that different units are combined with one another. Various units (depicted as functional units) can be designed as hardware, software or a combination of hardware and software. Control and processing unit 200 obtains magnetic resonance images of patient 120 according to a gradient echo pulse sequence, as described in further detail below. Control and processing unit 200 is interfaced with magnetic resonance imaging scanner 100 for receiving acquired images and for controlling the acquisition of images. Control and processing unit 200 receives image data from RF unit 1 10 and processes the imaging data according to the methods described below.

FIG. 1 provides an example implementation of control and processing unit 200, which includes one or more processors 202 (for example, a

CPU/microprocessor), bus 204, memory 206, which may include random access memory (RAM) and/or read only memory (ROM), one or more optional internal storage devices 208 (e.g. a hard disk drive, compact disk drive or internal flash memory), a power supply 210, one more optional

communications interfaces 21 2, optional external storage 214, a display 216, and one or more input/output devices and/or interfaces 218 (e.g., a receiver, a transmitter, a speaker, a display, a user input device, such as a keyboard, a keypad, a mouse, a position tracked stylus, a position tracked probe, a foot switch, and/or a microphone for capturing speech commands). Control and processing unit 200 may include many more or less components than those shown.

Although only one of each component is illustrated in FIG. 1 , any number of each component can be included in control and processing unit 200. For example, a computer typically contains a number of different data storage media. Furthermore, although bus 204 is depicted as a single connection between all of the components, it will be appreciated that the bus 204 may represent one or more circuits, devices or communication channels which link two or more of the components. For example, in many computers, bus 204 often includes or is a motherboard.

In one embodiment, control and processing unit 200 may be, or include, a computer or any other hardware equivalents. Control and processing unit 200 may also be implemented as one or more physical devices that are coupled to processor 202 through one of more

communications channels or interfaces. For example, one or more

components of control and processing unit 200 can be implemented using application specific integrated circuits (ASIC). Alternatively, control and processing unit 200 can be implemented as a combination of hardware and software, where the software is loaded into the processor from the memory or over a network connection.

Control and processing unit 200 may be programmed with a set of instructions which when executed in the processor causes the system to perform one or more methods described in the disclosure. For example, as shown in FIG. 1 , control and processing unit 200 may be programmed with instructions in the form of a set of executable modules 240 (the executable instructions may be stored in memory 206), such as, but not limited to, pulse sequence generation module 242, image acquisition module 244, image reconstruction module 246, and B0 mapping module 248.

Modules 242, 244 and 246 may be implemented using algorithms known to those skilled in the art for pulse sequence generation, image acquisition, and image reconstruction, respectively. Pulse generation module 242 establishes the sequence of RF pulses and magnetic field gradients depending on the desired imaging sequence, and image acquisition module 244 that stores the MR signals detected by the coils 108 and/or 106 in raw data space. Image reconstruction module 246 processes the acquired image data to perform image reconstruction of an MR image. According to various aspects of the present disclosure, a BO mapping module is provided for processing image data and generating a corrected BO map. The corrected BO map be employed by image reconstruction module 246. While the control and processing system of FIG. 1 is shown as implementing the BO mapping algorithm in the present example implementation, it will be understood that the methods of the present disclosure need not be performed on a computing device that is directly or indirectly connected to the magnetic resonance scanner, and that the methods of the present disclosure, or portions of the methods of the present disclosure, may be performed on a separate and disconnected computing system to post-process the image data,

Some aspects of the present disclosure can be embodied, at least in part, in software. That is, some methods can be carried out in a computer system or other data processing system in response to its processor, such as a microprocessor, executing sequences of instructions contained in a memory, such as ROM, volatile RAM, non-volatile memory, cache, magnetic and optical disks, or a remote storage device. Further, the instructions can be downloaded into a computing device over a data network in a form of compiled and linked version. Alternatively, the logic to perform the processes as discussed above could be implemented in additional computer and/or machine readable media, such as discrete hardware components as large- scale integrated circuits (LSI's), application-specific integrated circuits (ASIC's), or firmware such as electrically erasable programmable read-only memory (EEPROM's).

A computer readable medium can be used to store software and data which when executed by a data processing system causes the system to perform various methods. The executable software and data can be stored in various places including for example ROM, volatile RAM, non-volatile memory and/or cache. Portions of this software and/or data can be stored in any one of these storage devices. As used herein, the phrases "computer readable material" and "computer readable storage medium" refers to all computer- readable media, except for a transitory propagating signal per se. Examples of computer-readable media include but are not limited to recordable and non- recordable type media such as volatile and non-volatile memory devices, read only memory (ROM), random access memory (RAM), flash memory devices, floppy and other removable disks, magnetic disk storage media, optical storage media (e.g., compact discs (CDs), digital versatile disks (DVDs), etc.), among others.

As noted above, various embodiments of the present disclosure provide systems and methods for the generation of a corrected BO map, and for the generation of images, and/or other measures of interest, based on a corrected BO map, and based on the processing of multi-echo gradient echo data.

Referring now to FIG. 2A, a flow chart is shown that illustrates an example method of performing BO mapping according to the present disclosure. As shown at 300, multi-echo gradient echo image data is received. For example, an MR scanning system may be employed to collect MR image data from a subject using a multi-echo pulse sequence known to those skilled in the art. Non-limiting examples of such sequences include IDEAL/IDEAL-IQ on a GE scanner, mDIXON on a Philips scanner and MEGE on a Siemens scanner.

As shown at 310, a complex data set is generated based on the multi- echo image data. The complex data set is generated from the multi-echo image data by selecting echoes such that fat signals associated with the selected echoes are approximately in phase (pseudo-in-phase), such that the unwrapping of the phase complex data set is suitable for generating an initial BO phase map {φω,ίηΐ) - In some embodiments, the complex data set is obtained by calculating a Hermitian product ("synthetic Hermitian

multiplication") of the signals from the selected echoes. An effective echo time {ATE _b o) can then be derived from the known echo times and

subsequently used to define a scaling factor ( 1/2πΑΤΕ _Μ ) for converting phase to frequency. This scaling factor may be employed, for example, for converting a BO phase map into a BO frequency map.

Various example methods of selecting appropriate echoes for generating the complex data set, henceforth referred to as S _b o, are now described. The selection of the echoes for calculating S _b0 is based on the determination of a suitable or optimal ATE _b o, which is the ATE closest to the time period required for fat to experience a phase shift of 2π (or multiples of 2π). Different echo-selection strategies are implemented in according to the example methods presented herein, in order to accommodate acquisitions with both uniformly and non-uniformly spaced echoes. It will be shown that both strategies may be based on estimating the values of the relative phase difference between fat and water for different echoes.

The phase difference can be estimated as follows. When BO off- resonance is present in raw multi-echo data, the signal (S) from a pixel or voxel (it is noted that the terms voxel and pixel are used interchangeably in the present disclosure, where a pixel pertains to a spatial location within a two-dimensional slice, and a voxel refers to a spatial location within a three- dimensional volume) containing water and fat at the f ^h echo can be written, using a multi-peak fat model, as follows (5): where pw and pF are the signals of the water and fat components,

respectively; a _m is the known relative intensity of the m ^th peak of the fat spectrum, Af _m is the corresponding relative frequency shift from water, and M is the total number of peaks of the multi-peak fat model.

Neglecting the unknown phase terms in p _w and pF, Eq. 1 can be rewritten into a magnitude term {S _mag ) and a phase term (φ), as follows:

S _mag (FF,TE _} ,R ₂ ^* ) =

[2a]

<p(FF,TE _} ,f _b0 = wrap((p _fw (FF,TE _j ) + 2nf _b0 TE _j ) , [2b] where the fat-fraction FF is defined as and where "wrap" is an operator that constrains phase values to the range of [-π π). Given this expression for the fat-fraction, the phase term (<p _w ) related to the chemical shift is then:

<p _wJ (FF, TED = angle((l - FF) + FF∑£ a _{m β} ^ί2πΔ ^ ^ΤΕ ί). [2c]

The phase difference among different echoes due to the presence of fat, A<pfw ni, can thus be calculated or estimated based on Eq. 2c, where FF = 1 . It will be understood that equation 2c may be employed as modeling one or more fat peaks. For example, in some embodiments, a 6-peak fat model (19) may be employed.

A specific example implementation for the selection of suitable echoes for the generation of the complex data set is illustrated in the flow chart shown in FIG. 3. For datasets acquired with uniformly spaced echoes, as shown at 350 and 360, ATE _b0 is always equal to a multiple of the echo spacing ATE , i.e. ATEbo = IboATE, where Ibo is an integer. Thus, optimizing ATE _b o through echo-selection equates to estimating an optimal value of l _b o- With the single- peak fat model, the optimal R value is equal to

with ATE≠ 1/ f _cs . Similarly, with the multi-peak model, the optimal Ibo value is estimated by round{2n/(Z ^N _j i ^~ ₁ ¹ |wrap(<p _wJ+1 - <p _wj ) |/(N _e - 1))), where N _e denotes the total number of echoes. The echoes for generating Sbo are then identified given the knowledge of the optimal o and the known N _e .

As shown in FIG. 3 at 370, if Ibo is less than N _e for image data based on evenly spaced echoes, two echoes with an index difference equal to Ibo are selected. On the other hand, as shown at 375, if Ibo is equal to or larger than N _e , selecting more than two echoes and designing a specific formula for Sbo may be needed.

For example, when both the optimal Ibo and N _e are equal to three, one can generate S _b o from all three echoes using the formula (( ⁵ 3 : .) ( . ⁵ i .)) ; in this case O TE _b0 is equal to ( ^2TE 3 ^{~ TE} z ^{~ TE} i) (i.e. 3Ό ΤΕ). Practically, when an operator prescribes scanning parameters, the TEs are known or have been pre-selected to the optimized values. Therefore, selecting appropriate echo times and building a formula for constructing the complex data set for calculating BO can be achieved automatically. For example, in the case of N _e = 4, l _b o = 5, one can use (S ₄ S ^* i)( S ₃ S ^* i). Typically, when a sequence with more than three TEs is implemented on a scanner or used with a preset protocol, a pair of echoes is usually already suitable for BO mapping.

As shown at 355 in FIG. 3, in order to select suitable echoes from a dataset acquired with non-uniformly spaced echoes, values o\ Δφ^ _ιηι axe calculated over all possible echo-pairs, and the paired echoes with less than a pre-selected (e.g. user-defined) threshold value relative to a multiple of 2π (for example, ±0.1 π) are identified. In one example implementation, if only a single pair is identified, the two echoes included in this pair are selected; while if multiple paired echoes are identified, the pair with the minimum TE difference and the shortest TE value among all identified pairs, is selected, because high SNR in Sbo is desired.

Referring again to 310 of FIG. 2A, once the echoes are selected, the complex data set, S _b0 is generated, as described above. It is noted that the selection of suitable echoes results in a complex data set based on fat signals that are "pseudo-in-phase".

After having generated the complex data set, the phase of the complex data set it is unwrapped, as shown at 31 5 in FIG. 2A, using a suitable phase unwrapping method. It will be understood that there are many different phase unwrapping methods that are known in the art, and which may be employed to perform the phase unwrapping step. A non-limiting example of a suitable phase unwrapping method is the PUROR unwrapping algorithm (20). The phase unwrapping step results in an initial BO map, which may be generated in the form of a phase map, frequency map, or any other suitable method map form.

The unwrapped phase term of Sbo corresponds to the phase difference between the echoes and is the sum of phase-shifts caused by different sources:

<PbO,ini ⁼ VbO + ψηη + <Pfw,ini > [3] where φ ₀ ο is caused by BO inhomogeneity; (p _U n ^rmr represents the unwrapping phase error (unknown multiple of 2π phase shifts); and is a bias term related to the chemical-shift between fat and water and is dependent on ATEbo- The challenge is to remove the two error terms {<p _un ^error and (pfw ni) - The first challenge in obtaining an accurate BO map is to remove ( un " ^m - Unwrapping errors can affect the initial BO map globally or locally. Global phase unwrapping errors arise because phase values in the unwrapped image are relative, i.e. the entire initial BO phase map may be shifted by a multiple of 2π (we define the term <p _g ^un , the subscript g indicating a global effect). Second, local phase errors {φ , the subscript / referring to a local effect) may exist in regions where larger-than-π phase gradients exist in the real BO map. It is important to note that for all phase unwrapping methods - which incorporate a constraint that the spatial phase gradient is less than π - the performance of phase unwrapping is sensitive to noise, large jumps at the interface between tissue and air, large BO inhomogeneity, and singularities.

The presented example method identifies and corrects both global and local phase errors, as shown at step 325 in FIG. 2A, by evaluating a B0- corrected phase image of the Hermitian-product between the j ^th and (j+1 ) ^th echoes (A(pjj ₊ i ^orl ). The BO-corrected phase (Δφ ^οοη ) is given by:

Αφ^ = Αφ°™ ₊₁ - φ„ _0Μί ΔΤΕ _]ι]+1 /ΔΤΕ _Μ [4a]

= (A<p°y ₊₁ - 2π χ / _Μ ΔΤ¾ ₊₁ ) + <p ^ror [4b] where ATE _{j +1} is the difference of TEs between two consecutive echoes {TE _j+1 - TE _j ); (Pun ^or includes two terms caused by <p _g ^un and φ '-

<pu ^er n ^ror = -« ⁿ + <ΡΓ)Δ¾ _{+ 1} /ΔΤΕ _Μ . [4C]

The evaluation process is based on the following observations: if no phase error is present in the initial BO map all water-only pixels have near zero phase values (Δφ ^οοιτ ~ 0) ; the fat-only pixels have phase values close to a phase term (A<p ^fat ) approximated using Eq. 2c with FF = 1 ; the other pixels have phase values between 0 and A<p ^fat . These observations lead to a first criterion for guiding phase-error correction: the values of Δφ ^οοι of all tissue- pixels should be in the range between zero and A<p ^fat . However, the first criterion is not sufficient to differentiate fat/water swaps, especially when the images or regions are predominantly occupied by water or fat. To solve this ambiguity, the present example method generates two separate water masks - a phase based water mask from a thresholded Δφ ^οοιτ and a magnitude-based water mask from a thresholded magnitude- based FF map - and imposes a second criterion that the two masks must match. The magnitude-based water mask is generated as shown in step 305 of FIG. 2A, and is subsequently employed in global phase error correction, as shown at step 325.

The magnitude-based water mask can be generated according to a number of different methods. In one example, in which three or more echoes are detected, the following method may be employed, as illustrated in the flow chat shown in FIG. 4. Theoretically, to achieve full-range FF mapping (0- 100%) only based on magnitude information a multi-peak fat model must be used. Prior to applying pixel-by-pixel magnitude fitting using Eq. 2a, it was demonstrated that fat pixels fall along a non-exponential curve, while the water curve follows the expected exponential decay.

The first step in the present example magnitude fitting method shown in FIG. 4 involves normalization of the (|pw| + IP F|) term in Eq. 2a over all echoes by using the magnitude value from the first echo; correspondingly TEi is subtracted from the TEs of all echoes. To estimate the two remaining unknown variables (FF and R2 ^* ) in Eq. 2a, a set of trial fat-fraction maps are generated by setting the fat-fraction of all voxels equal to a set of discrete values between 0 and 1 (FF _f , t denotes the trial value), as shown at 400 in FIG. 4. For each FF _t , a residue map is also calculated by varying the R2 ^* value and comparing against the measured magnitude - the R2 ^* yielding the lowest residue {RES _t ) for a given pixel is used to generate the trial R2 ^* map (/¾f ^* ) _> as shown at 405 and 410. By repeating the above procedure several times with different FF values, an R2 ^* map and a magnitude-based FF map can be obtained by calculating the weighted average of the R2 ^* and FF used in each FF _t iteration above, as shown at 415 and 420:

and

FF _mag = (∑i ₌₁ w _t FF _t )/(∑[ ₌₁ w _t ), [5b] where T is the total number of the tested FF values and the weighting factor

From the magnitude-based FF map, a threshold value (e.g. a preselected threshold, such as a user-specified threshold) can be used to generate a binary magnitude-based water mask, as shown at 425, for use in verifying the phase errors in the aforementioned global phase error correction method (i.e. in step 325 of FIG. 2A. It is noted that the R2 ^* map generated from Eq. 5a can be used as the final R2 ^* map.

In one example implementation, the range of FF values is set from 0 to 1 .0, for example, in steps of 0.1 ; the range of R2 ^* is set, for example, from 0 to 300 s ^"1 with an interval of 5 s ^"1 ; pixels with a magnitude-based FF less than 0.3 are included in the magnitude-based water-mask. It is noted that the aforementioned range of R2 ^* values is prescribed to pertain to tissue properties, which are commonly falling in this range. The interval of 5 s ^"1 is selected for considering both accuracy and computational cost: a smaller value is better but will increase computational time, and the same trade off pertains to the selection of the values for the FF mask.

In the preceding embodiment, unipolar multi-echo (total echo number >

2) GRE data are acquired and magnitude-based masks are generated through pixel-by-pixel fitting. Alternatively, the various embodiments of the present disclosure can be adapted to process a conventional dual-echo GRE data set. More particularly, because pixel-by-pixel fitting cannot be performed for dual-echo acquisitions, the magnitude-based masks can be generated based on the difference of signal intensity or texture features between fat and water based on the difference of signal intensity or texture features between fat and water.

The magnitude-based water map, which may be generated as described above, is then employed for global phase error correction in step 325 of FIG. 2A and 2B. An example method of global phase-error correction is illustrated in the flow chart shown in FIG. 5. In this present example method, both the global and local phase error corrections are performed sequentially, on a slice-by-slice basis. To determine the unknown global phase shift, a set of BO map templates is generated by shifting the initial BO phase map with a set of multiples of 2π, as shown at 500 in FIG. 5:

Vbofni = VbO.ini + ²ⁿⁿ g,t > [6] where n _{g t} is an integer in the range of ±(round(ATE _b0 /ATE)-1). Each trial B0 map is used to perform background removal using Eq. 4a. The pixels with phase values of Δφ ^οοι less than 0.1 π are included in a phase-based water mask, as shown at 505.

For each trial BO map, the matched ratio between the corresponding trial phase-based water mask and the magnitude-based water mask is calculated as the total number of pixels identified by both water masks divided by the total number of pixels included in the magnitude-based mask. Among these trials, the one with the highest match ratio is considered as the correct one, as illustrated at 510 in FIG. 5. In the case where the same match ratio is obtained from two different values of n _g,t , separated by round(ATE _b o/ATE), the one with the smaller absolute value of n _g is selected. Once the value of n _g is determined, the initial BO phase map is shifted by 2πη ₉ , resulting in a <p _g ^un - corrected BO map.

After performing global phase-error correction, local phase-error correction is performed on the globally phase-error corrected BO map.

However, prior to performing this step, a bias-corrected value for the fat phase shift {AqJ ^at ) that accrues during acquisition is first calculated, as shown at 320 in FIG. 2A.

The process of local phase-error correction may be performed as illustrated in the flow charts shown in FIGS. 6A-C. As illustrated at 600 in FIG. 6A, pre-selected (e.g. user-specified) thresholds of the phase gradient (π) are used to divide the ^ "-corrected B0 map into regions. For each region, the total number of pixels that have phase values of Δφ ^οοι (calculated as described above) in the range of [-0.25xAq> ^fat A<p ^fat ] and [A<p ^fat -0.25xAq> ^fat ] for positive and negative A<p ^fat , respectively, are determined; and an "in-range ratio" is calculated as the total number of pixels that fall in the above range divided by the total number of pixels, as shown at 61 0.

The value for Aq> ^fat is calculated as follows. Eq. 2c is employed, with FF=1 , to calculate the phase shift Αφ^ between echoes due to the chemical shift due to fat. Using the values of A<p _{fw n} i and ATE _b0 , the R value is the recalculated using the bias-corrected fat model in order to obtain the value of Αφ ^ί3ί . This can be achieved by calculating f _Uas = A(p _fW j _n (2%ATE _b o) then replacing the original chemical shift with the bias-corrected values (i.e. A† _m,b ias

= A† _m — fbias)■

The in-range ratio is calculated for each trial BO map with ni _,t - an integer in the range of ±(round(ATE _b o/ATE) -1 ). To satisfy the first criterion, the ni t that results in the highest in-range ratio is selected and used to shift the phase values of the region by 2πηι, as shown at 615. It is noted that the process of identifying rii and n _g is the same.

After having performed the local phase-error correction illustrated in

FIG. 6A, further local phase-error corrections may be optionally determined. A first optional additional local phase-error correction method is illustrated in the example flow chart shown in FIG. 6B. In step 630, "out-of-range pixels" having phase correction values outside of the aforementioned range are identified. As shown at 635, the out-of-range pixels are then divided into regions (e.g. by identifying connected out-of-range pixels). The same strategy used in the FIG. 6A may then applied to determine the value of ni, as shown at 640-650. It is noted that if two different values of ni _,t result in the same in-range number, the one with the smaller absolute value may be selected for the example methods illustrated in FIGS. 6A and 6B. FIG. 6C illustrates another optional local phase-error correction that may be applied. As shown at 670, the phase-based water mask is regenerated using the phase map that is locally corrected using the

aforementioned methods. Any pixels that are identified by the phase-based water mask, but not the magnitude-based water mask, are considered "unmatched" pixels and are grouped into regions of connected pixels, as shown at 675 and 680. If the region size is smaller than a pre-selected threshold (e.g. a threshold of at least 1 00, 250, or 500 pixels) then no further processing is performed. For each larger region, as shown at 680-695, the value of ni is determined as described in FIGS. 6A and 6B, in order to shift the mean phase value of Δφ ^οοι ίο be close to zero.

Following correction for the global and local phase errors, and having thereby obtained a phase-corrected B0 map, removal of the bias term ( , _η/ ) must be performed in order to generate an accurate B0 map (as can be seen from eqn. 3). Removal of the term can be achieved in the pulse- sequence design using the conventional Dixon-based B0 mapping

approaches which rely on true in-phase Sbo {i.e., 0). For the

conventional Dixon-based techniques, both echo times and the formulas for calculating Sbo are designed to ensure that for a single-peak fat/water chemical shift {† _cs ) 27rf _cs ATE _b o is equal to a multiple of 2π. However, for a general multi-echo acquisition, the /™ = ⁰ condition is rarely met and a bias term is introduced in the initial B0 map - primarily affecting fat-dominant pixels and resulting in a biased FF map if the standard multi-peak fat model (19) is used.

Referring again to FIG. 2A, the bias correction is applied as follows. A corrected FF map is obtained first obtained, as shown at 330. To obtain a corrected FF map from the phase-error corrected BO map, the fat model is modified by calculating† _bias = A<p _{fw n} /{2nATE _b0 ) (as described above and employed in the calculation of Δφ ^ί3ί ) and then replacing the original chemical shift with the bias-corrected values (i.e. A† _m , _b ias = A† _m - f _b i _as ) - The corrected FF map may thus be calculated, for example, using a Penrose-Moore

pseudoinverse method (5), employing the phase corrected BO map and the bias chemical shift. The final BO may be expressed as a frequency map using ATEbo.

In embodiments in which the multi-echo data includes three or more echoes, and where and R2 ^* map (or a T2 ^* map) is obtained during the calculation of the magnitude-based water map (e.g. as shown in FIG. 4), then the relaxation map may also be used in the calculation of the corrected fat- fraction map. This example embodiment is illustrated in FIG. 2B, where the relaxation map is generated in 305A and employed for the calculation of the corrected fat-fraction in 330A. It is noted that for cases in which more than two echoes are acquired, it has been widely accepted that including the relaxation map will improve accuracy of the fat-fraction map, resulting in more accurate bias removal. However, the aforementioned method may alternatively be performed without including employing the relaxation map. It is further noted that when two echoes are acquired (dual echo), the relaxation effects are typically relatively small and negligible. It will be understood that the aforementioned example methods may be implemented using image data that is collected using two echoes, or three or more echoes (the example magnitude method illustrated in FIG. 4 may be employed using three or more echoes). Furthermore, the fat model may include a single peak, or may be a more precise multi-peak model. Generally, it has been found that superior results are obtained when three or more echoes are used, and where a multi-peak fat model is employed.

In cases where image data is acquired with three or more echoes, iterative nonlinear data fitting algorithms with a multi-peak model can be adapted (10, 18). In general cases with two or more echoes, as noted above, the magnitude-based water mask can be estimated based on the difference of image intensity or texture features between water and fat.

In some example embodiments, the BO mapping methods disclosed herein may be implemented using data acquired with bipolar excitation pulses. In such implementations, before performing the BO mapping method, phase errors due to gradient delays should be removed. Specifically, zero and first-order phase corrections should be performed on an echo-by-echo basis.

As noted above, the BO mapping methods disclosed herein may be employed for data acquired using single- or multi-channel RF receive coils. When MR signals are collected by a multi-channel RF coil, complex images are reconstructed on a channel-by-channel basis. The resulting channel-by- channel complex images can be combined into a channel-combined data set by using phase sensitive approaches. When the present BO mapping methods are employed to process channel-combined data set, the

aforementioned implementation can be directly used to generate the final BO map.

When the BO-NICE method is applied to each individual channel data, the final BO map is calculated as magnitude-weighted mean over all channels.

EXAMPLES

The following examples are presented to enable those skilled in the art to understand and to practice embodiments of the present disclosure. They should not be considered as a limitation on the scope of the disclosure, but merely as being illustrative and representative thereof. Throughout these examples, the aforementioned methods of BO mapping are referred to as the "BO-NICE" method, where the acronym "NICE" denotes the Non-Iterative Correction of phase-Errors.

FIGS. 7A-C illustrate an example implementation based on imaging (unipolar 6-echoes, ΤΕ1 /ΔΤΕ = 1 .2/1 .0 ms) of a whole liver, which was processed as an example implementation of BO-NICE method illustrated in FIG. 2B. FIGS. 7A and 7B show the representative final BO maps and R2 ^* maps, respectively, of three consecutive slices. FIG. 7C shows the FF maps of the three slices obtained by using the GE IDEAL method, the IDEAL method modified according to the present BO-NICE method (in which the BO map is generated using BO-NICE, then the BO map is removed from the data set, resulting in a BO-corrected data set, and where The IDEAL technique is then applied to the BO-corrected data set), and the BO-NICE methods, from left to right, respectively. The results shown in FIG. 7C suggest the accuracy of BO maps and demonstrate that BO-NICE method significantly improves the performance of IDEAL, where IDEAL introduces fat/water swaps. Note that channel-combination was conducted before applying the BO-NICE method. FIGS. 8A-C illustrate an example (bipolar 5-echoes, ΤΕ1 /ΔΤΕ = 3.0/1 .4 ms) of a whole head image acquisition, showing three different slices, respectively. As with the unipolar multi-echo acquisition, R2 ^* map and magnitude-based FF map were generated from the magnitude only. However, an extra step has to be applied for bipolar acquisition in order to correct gradient-delay-related phase errors.

After the aforementioned step has been conducted, the BO-NICE method is applied to each individual channel data, followed by using magnitude-weighted channel combination. As demonstrated by FIGS. 8A-C, the BO-NICE implemented on a channel-by-channel basis works well for a 3D whole brain bipolar data set.

FIGS. 9A-G illustrate an example (bipolar 2-echoes, TE1 /TE2 = 2.3/4.6 ms) of a whole-heart acquisition. Among Dixon techniques, the bipolar dual- echo technique is intrinsically the fastest, making it optimal for high resolution or 3D applications. To perform dual-echo fat-suppression, phase errors due to gradient delays should be removed; and if the removal is perfect, there should be no difference of conducting fat/water separation between bipolar and unipolar acquisition. To remove the unwanted BO off-resonance effect, the convention dual-echo approach maps BO simply using unwrapping but without phase error-correction, i.e., the outcome is determined by how well unwrapping performs. As shown in FIG. 9A, magnitude-based water- and fat- masks are generated, based on the fact that signal intensity from fat surrounding the chest wall is much higher than that from the heart for T1 - weighted images. These masks are then used to correct phase errors (arrow in FIG. 9C) for generating phase error-corrected BO maps (see FIG. 9D) for each channel, followed by combining them into final BO maps (FIG. 9E) by using magnitude-weighted mean method. Note that the initial BO map shown in FIG. 9C is obtained from FIG. 9B with phase unwrapping. Water-only and fat-only images are reconstructed using the two-point Dixon approach on a channel-by-channel basis. Final channel combination of water-only and fat- only images (FIG. 5F and 5G) is performed using the sum of squares approach. This example demonstrates that extending BO-NICE to dual echo acquisition will improve the robustness of BO mapping and fat-suppression using phase unwrapping.

The specific embodiments described above have been shown by way of example, and it should be understood that these embodiments may be susceptible to various modifications and alternative forms. It should be further understood that the claims are not intended to be limited to the particular forms disclosed, but rather to cover all modifications, equivalents, and alternatives falling within the spirit and scope of this disclosure.

REFERENCES

Reeder SB, Pineda AR, Wen Z, Shimakawa A, Yu H, Brittain JH, Gold GE, Beaulieu CH, Pelc NJ. Iterative decomposition of water and fat with echo asymmetry and least-squares estimation (IDEAL):

application with fast spin-echo imaging. Magn Reson Med

2005;54(3):636-644.

Glover GH, Schneider E. Three-point Dixon technique for true water/fat decomposition with BO inhomogeneity correction. Magn Reson Med 1991 ;1 8(2):371 -383.

Szumowski J, Coshow WR, Li F, Quinn SF. Phase unwrapping in the three-point Dixon method for fat suppression MR imaging. Radiology 1994;1 92(2):555-561 .

Hernando D, Kellman P, Haldar JP, Liang ZP. Robust water/fat separation in the presence of large field inhomogeneities using a graph cut algorithm. Magn Reson Med 2010;63(1 ):79-90.

Yu H, Shimakawa A, McKenzie CA, Brodsky E, Brittain JH, Reeder SB. Multiecho water-fat separation and simultaneous R2 ^* estimation with multifrequency fat spectrum modeling. Magn Reson Med

2008;60(5):1 122-1 134.

Reeder SB, Wen Z, Yu H, Pineda AR, Gold GE, Markl M, Pelc NJ. Multicoil Dixon chemical species separation with an iterative least- squares estimation method. Magn Reson Med 2004;51 (1 ):35-45. Yu H, Reeder SB, Shimakawa A, Brittain JH, Pelc NJ. Field map estimation with a region growing scheme for iterative 3-point water-fat decomposition. Magn Reson Med 2005;54(4):1032-1039.

Lu W, Lu Y. JIGSAW: Joint Inhomogeneity estimation via Global Segment Assembly for Water-fat separation. IEEE Trans Med Imaging 201 1 ;30(7):1417-1426.

Soliman AS, Yuan J, Vigen KK, White JA, Peters TM, McKenzie CA. Max-IDEAL: A max-flow based approach for IDEAL water/fat separation. Magn Reson Med 2013.

Yu H, Reeder SB, Shimakawa A, McKenzie CA, Brittain JH. Robust multipoint water-fat separation using fat likelihood analysis. Magn Reson Med 2012;67(4):1065-1076.

Tsao J, Jiang Y. Hierarchical IDEAL: fast, robust, and multiresolution separation of multiple chemical species from multiple echo times. Magn Reson Med 2013;70(1 ):155-159.

Lu W, Hargreaves BA. Multiresolution field map estimation using golden section search for water-fat separation. Magn Reson Med 2008;60(1 ):236-244.

Jacob M, Sutton BP. Algebraic decomposition of fat and water in MRI. IEEE Trans Med Imaging 2009;28(2):173-184.

Cui C, Wu X, Newell JD, Jacob M. Fat water decomposition using globally optimal surface estimation (GOOSE) algorithm. Magn Reson Med 2014.

Berglund J, Kullberg J. Three-dimensional water/fat separation and T2 ^* estimation based on whole-image optimization-application in breathhold liver imaging at 1 .5 T. Magn Reson Med 2012;67(6):1684- 1693.

Hernando D, Haldar JP, Sutton BP, Ma J, Kellman P, Liang ZP. Joint estimation of water/fat images and field inhomogeneity map. Magn Reson Med 2008;59(3):571 -580.

Reeder SB, Cruite I, Hamilton G, Sirlin CB. Quantitative assessment of liver fat with magnetic resonance imaging and spectroscopy. J Magn Reson Imaging 201 1 ;34(4):729-749.

Yu H, Shimakawa A, Hines CD, McKenzie CA, Hamilton G, Sirlin CB, Brittain JH, Reeder SB. Combination of complex-based and

magnitude-based multiecho water-fat separation for accurate quantification of fat-fraction. Magn Reson Med 201 1 ;66(1 ):199-206. Hernando D, Liang ZP, Kellman P. Chemical shift-based water/fat separation: a comparison of signal models. Magn Reson Med

2010;64(3):81 1 -822.

Liu J, Drangova M. Intervention-based multidimensional phase unwrapping using recursive orthogonal referring. Magn Reson Med 2012;68(4):1 303-1316.