Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Magn Reson Imaging. Author manuscript; available in PMC 2010 April 1.
Published in final edited form as:
PMCID: PMC2714529

Removing Background Phase Variations in Susceptibility Weighted Imaging Using a Fast, Forward-Field Calculation

Jaladhar Neelavalli, Ph.D.,1,* Yu-Chung N. Cheng, Ph.D.,1,2 Jing Jiang, M.S.,3 and E. Mark Haacke, Ph.D.1,2,3



To estimate magnetic field variations induced from air-tissue interface geometry and remove their effects from susceptibility weighted imaging (SWI) data.

Materials and Methods

A Fourier-transform-based field estimation method is used to calculate the field deviation arising from air-tissue interface geometry. This is accomplished by, first, manually drawing or automatically detecting the sinuses, the mastoid cavity and the head geometry. The difference in susceptibility, Δχ, between brain tissue and air-spaces is then calculated using a residual-phase minimization approach. SWI data are corrected by subtracting the predicted phase from the original phase image. Resultant phase images are then used to perform the SWI post-processing.


Significant improvement in the post-processed SWI data is demonstrated, most notably in the frontal and midbrain regions and to a lesser extent at the boundary of the brain. Specifically, there is much less dropout of signal after phase correction near air-tissue interfaces making it possible to see vessels and structures that were often incorrectly removed by the conventional SWI post-processing.


The Fourier-transform-based field estimation method is a powerful 3D background phase removal method for improving SW images, providing clearer images of the fore-brain and the mid-brain regions.

Keywords: Susceptibility Weighted Imaging, air tissue interface, phase aliasing, susceptibility distribution, Fourier Transform, field calculation


Susceptibility weighted imaging (SWI) has become a powerful clinical tool for revealing the presence of microhemorrhages, iron and calcium and, therefore, for studying conditions like aging and diseases such as multiple sclerosis, stroke, trauma and tumors (18). However, SWI still suffers from problems caused by high spatial frequency phase changes and phase aliasing in the background magnetic field caused by the presence of air-tissue interfaces, particularly in the mid-brain and the forebrain regions near the paranasal sinuses. Although the presence of low spatial frequency fields can be removed using a high pass (HP) filter approach (9, 10), problems still exist in terms of rapid unwanted field variations particularly near the mastoid cavity, frontal, ethmoid and sphenoid sinuses and to a lesser degree the maxillary sinuses. There have been some attempts to remove these problems by a combination of phase unwrapping the images first and then either high pass filtering the data or doing a polynomial fit to pull out the background fields (11). The goal in all of these phase machinations is to leave behind, preferably unaltered, the local phase information (arising from local susceptibility differences) from structures such as the veins, iron laden tissue and calcifications. However, each of these methods has its own drawbacks. While a simple HP filtering is able to remove low spatial frequency phase variations, usually a very strong HP filter is needed to remove the rapid phase aliasing near the air-tissue interface, which results in a concomitant loss of important local phase information (10). On the other hand, removal of aliasing by a fitting approach requires local polynomial fits throughout the brain on a slice by slice basis, often with different regions in a given slice fitted with different order polynomials.

Since the physical basis of such phase aliasing effects is clearly known, which is the static field deviation due to the air-tissue interface, the ideal solution would be to estimate the field effects within the brain from these geometries and then to subtract out the resulting phase from the SWI phase images prior to any further phase processing (i.e. unwrapping/homodyne filtering (12) etc.). The process for this would be to take the geometry of the object of interest, such as the brain, and separate the tissue from the air, and then estimate the geometrically induced local fields. We determine the geometry directly from the MR data. After subtracting these effects, the remnant phase would only be from background static field effects arising from eddy currents and the local tissue susceptibilities. In this paper, we apply a new Fourier-transform-based field estimation approach (1315) to calculate an estimate of the geometry-induced field changes and use this to remove the spurious phase arising from the air-tissue geometry. Although such a process may not be perfect, we hypothesize that even an approximation is good enough to remove almost all the harmful aliasing. The subsequent use of a mild (or small) high pass filtering should provide ‘cleaner’ phase images and consequently better SWI processed images without the loss of important phase information either near to or far from the air-tissue interfaces.

Brief Theory of the Forward-Field Calculation

When an object with a certain magnetic susceptibility distribution, χ(r), is placed in an external magnetic field B0, the resulting magnetic field, B(r) at any position r is given by:


where M(r) is the induced magnetization distribution of the object. In the case of magnetic resonance experiments, the external applied magnetic field B0 is many orders of magnitude larger in one direction, assumed to be the z direction, than in the other two orthogonal directions. Taking this into account, and observing that the second term in the above equation (referred to as Bdz(r) in reference (14)) is a convolution (16, 17), this term can be Fourier transformed into a simple product relation in the spatial frequency domain (generally referred to as k-space) to give:


Here kx, ky, kz are the coordinates in k-space and Mz(k) is the Fourier transform of the magnetization distribution of the object along the z direction. This spatially varying field Bdz(r), which is the Fourier transformation of Bdz(k), gives rise to the spatially varying phase seen in a gradient echo MR experiment. Since the susceptibilities of biological tissues are typically much less than unity, the magnetization of the object is Mz(r) = χ(r)·B00 and hence Mz(k) = FT(χ(r))·( B00) where FT stands for Fourier transformation. Thus the phase due to a spatially varying field Bdz(r), which itself arises from the presence of the susceptibility distribution χ(r), can be calculated from:


Equation 3 offers a simple, fast and powerful means for calculating the field deviation due to the presence of a known χ(r) in an otherwise uniform field B0. It is important to note here that Eq. 2 and Eq 3 are derived under the assumption that B0 is the predominant magnetic field and that the Bdz field does not affect the magnetization of the object. This is referred as the Born approximation. This approximation and the linearity of the Fourier theory makes it feasible to simulate the magnetic field of any complicated structure, as the linear summation of the fields from several sub-structures, each of which has a spatially constant susceptibility χi (index i refers to the i-th sub-structure), using:


where gi is Fourier transformation of the geometry of the i-th sub-structure.


Phantom imaging

A series of phantom experiments were run to: a) validate that the FT based field estimation method can be used for removing effects of geometry-induced field variations; and b) to determine the sensitivity of the removal process to the estimated effective-Δχ values when two different echo time magnitude data sets are used. We constructed a cylindrical gel phantom containing an air filled sphere inside, in order to approximate the geometry of an air filled sinus embedded in tissue. An 8% gelatin (Kind & Knox Gelatine, Inc.) solution was prepared and poured into a polypropylene cylindrical container, filling 1/3 of its volume (180mm height, 120mm diameter, wall thickness 1.3mm). The gel was then allowed to partially settle. A hollow ping-pong ball, made of cellulose, with diameter 40mm and wall thickness of roughly 1mm was placed at the center of the container on the partially settled gel. The gel was then allowed to solidify so that the ball was held tightly. The remaining volume of the container was subsequently filled with a lukewarm 8% gelatin solution and the whole phantom was allowed to sit for 24 hours at room temperature and then was placed in a refrigerator.

The phantom was imaged in a 1.5T Siemens Sonata system with the following parameters: (i) T1-weighted 3D gradient echo data collected at two echo times (TE), 6.58ms and 7.58ms with TR = 15ms, FA = 6o, read bandwidth (BW) = 610 Hz/pixel, field-of-view (FOV) = 200×200mm, matrix size = 256×256×192 and slice thickness (TH) = 0.78mm giving an isotropic voxel size of 0.78mm. (ii) SWI data were collected at TE = 40ms, TR = 44ms, FA = 10o, BW = 610 Hz/pixel, FOV = 200×200mm, matrix size = 256×256×192 and TH = 0.78mm. The high bandwidth used in the T1 weighted sequence, not only helps in reducing image distortion, but also allows the use of a shorter TR, reducing the total imaging time. Phase images with an effective TE=1ms were obtained from the acquired T1-weighted data by complex dividing the data at TE=7.58ms with the data obtained at TE=6.58ms. This was performed to obtain a phase image with no or minimal phase aliasing.

Since the phantom images consisted of air or gel, a complex thresholding method (CTM) which uses the noise characteristics of both magnitude and phase images to distinguish noise pixels from signal was applied to obtain the geometry of the object. Briefly, based on the SNR of the image, the CTM algorithm uses a receiver operating characteristic curve to determine the appropriate choice of magnitude and phase thresholds along with neighborhood connectivity criteria to determine what pixels represent noise in the image (18). Two geometries, Ggre and Gswi, were obtained from the T1-weighted sequence with TE 6.58ms and from SWI images (at TE 40ms), respectively. The extracted geometry was placed in a 256×256×256 matrix and a value of unity was assigned to voxels containing gel signal and a value of 0 was assigned to the noise pixels.

Estimation of appropriate Δχ and implementation of modified phase processing

For estimating the two effective Δχair_gel values for Ggre and Gswi respectively, the following procedure was employed. The given geometry is used to estimate phase at 1ms effective echo time (using Eq. 4) with Δχair_gel values varying from -2ppm to -16ppm (SI units) in steps of 1ppm. Here, Δχ is equivalent to χ of water or gel as we assume susceptibility of air to be 0. Fifteen sets of residual phase maps, corresponding to each Δχair_gel value are generated, by subtracting the simulated phase values from the TE=1ms phase dataset. In each of the 15 residual phase sets, the standard deviation of the residual phase, σresidual_phase, was evaluated in a volume of interest (VOI- regions of interest chosen in multiple consecutive slices) containing roughly 8000 voxels, chosen in close proximity to the sphere. In principle, the subtraction should lead to zero residual phase when one of the assumed Δχ values equals the actual effective Δχair_gel. However, due to the presence of noise, effect of eddy currents, the constant phase shift due to spectrometer central frequency adjustment, the coarse step size in assumed Δχ values and the inherent error of the calculated phase (14, 15, 19), the standard deviation may not be zero. Nevertheless, at the Δχ nearest to the actual effective Δχair_gel value, the σresidual_phase will be minimal. To obtain the effective Δχair_gel value, σresidual_phase was plotted as a function of Δχ and the value corresponding to the minimal σresidual_phase was chosen as the respective effective Δχair_gel value.

Using these Δχair_gel values and their corresponding Ggre and Gswi geometries, phase at an echo time of 40ms was calculated and then complex divided into SWI phase data to remove the rapid phase aliasing in the vicinity of the sphere and at the edge of the phantom. The resultant phase images, referred to as the Geometry Dependent Artifact Corrected (GDAC) phase images, were further high pass filtered using a Hanning filter of size(9, 20) 64×64 (or 32×32, when indicated so) and were qualitatively compared with the corresponding HP filtered (size 64×64) original SWI phase images. We have employed the procedure outlined in reference (9) for high pass filtering. Both the HP-GDAC and HP-original SWI phase images were subsequently used to create phase masks which were multiplied 4 times with the magnitude to generate the processed susceptibility weighted images. The procedure outlined in reference (1) was used for generation of the phase mask. The results were compared qualitatively for the amount of signal recovered using GDAC phase as compared to using the original phase images.


The accuracy of the field estimation method using Eq. 3 is dependent on the ratio of the size of the object and the matrix size (i.e. object-size/FOV) in which it is placed during the computation of its spatial Fourier transform (14, 15, 19). To evaluate the influence of the object-size/FOV on the ability of the σresidual_phase minimization method to arrive at the appropriate Δχ value, the following simulations were performed. The simulations also help to study the behavior of the σresidual_phase minimization method under ideal conditions of no noise, eddy current effects, spectrometer shift or RF phase effects. A spherical shell with a susceptibility of −9.5ppm, inner radius Ri and an outer radius Re, placed in vacuum was considered. There is no susceptibility inside the inner sphere. This choice was so made because an analytic solution of the magnetic field due to a spherical shell magnetization is known, and on the first order the hollow inner sphere may be approximated to the air filled sinus within the brain tissue and also to the hollow sphere in the phantom study. For numerical simulations an Ri of 10 voxels and an Re of 55 voxels were chosen, with the shell placed at the center of a matrix of 128×128×128 voxels. These dimensions are close to the typical object-size/FOV ratios of the sinus and head, which are 0.15 and 0.85, respectively.

The analytic field due to the spherical shell was calculated using the following relation:

Bz(r)={    0whenr<Ri13(χiχe)B0(3cos2θ1)(Rir)3whenRi<r<Re    13[(χiχe)(Rir)3+χe(Rer)3]B0(3cos2θ1)whenRe<r

where, χi is the susceptibility of the inner sphere (the hollow region), χe is the susceptibility of the shell region, and θ is the angle between the main magnetic field B0 and the position vector of the point under consideration. We chose χi=0 and χe=−9.5 ppm and the analytic field due to the spherical shell was evaluated using Eq.5 from which the consequent phase at echo time 1ms and field strength 1.0T was then calculated.

To simulate the σresidual_phase minimization process for estimating actual effective Δχ for the given geometry, phase due to the shell geometry was calculated, using Eq. 4, with Δχ varied from −4ppm to −16ppm in steps of −1ppm. Between −9 to −11 ppm the step size was reduced to 0.1ppm to improve the accuracy of the determination of Δχ. The calculated phase values, for each value of Δχ were subtracted from the analytic phase (calculated using Eq. 5) through complex division to obtain 31 sets of residual phase maps. Standard deviation of the residual phase, σresidual_phase, from roughly 7000 (and 1000) voxels close to the inner sphere in each phase set, was evaluated and plotted as a function of Δχ. The value corresponding to the minimum σresidual_phase was the estimate of the desired Δχ value.

Human Imaging

Three healthy volunteers were imaged in a 1.5T Siemens Sonata magnet with the sequence parameters given in Table 1. Slightly different parameters were used in different volunteers due to magnet time constraints. The T1 weighted complex data obtained at two echo times in each volunteer were used to create phase images with an effective TE of 1ms for volunteers 1 and 2 and an effective TE of 2ms for volunteer 3 through complex division. The magnitude images from the T1-weighted data with the shorter echo time were used to obtain the geometries of the head and relevant air-spaces in each volunteer.

Table 1
List of sequences and imaging parameters for 3 volunteers.

The 3D T1-weighted gradient echo magnitude images acquired at the shorter echo time, 5ms for volunteers 1 & 2 and 7ms for volunteer 3, were used to extract geometries of maxillary, frontal, sphenoid, ethmoid, maxillary sinuses, mastoid cavity (a signal void due to the presence of bone and air cells) and the head. Similar to the phantom case, the CTM algorithm (18) was used to distinguish noise pixels from signal. The voxels identified as tissue voxels were assigned a value of unity and the noise to a value of zero. Although the algorithm was successful in identifying the sinuses andcavity regions automatically, some manual intervention was required around the skull and near the sinuses to avoid sharp edges in the boundaries. The geometry of the whole head down to the jaw region was included (15).

To estimate the susceptibility differences (Δχ) between the brain tissue and different sinuses, Δχsinus-tissue = χtissuesinus, and between brain tissue and the mastoid cavity, Δχmastoid-tissue = χtissuemastoid, the extracted geometries were used to simulate phase at 1ms, for volunteers 1 and 2, and at 2ms for volunteer 3, using Eq. 4 with varying Δχsinus-tissue and Δχmastoid-tissue from −2ppm to −16ppm. A zero susceptibility was assigned to all sinuses (i.e., χsinus = 0) and the regions outside the head. With brain tissue predominantly containing water, its χ value with respect to air, is generally expected to be in the neighborhood of −9ppm. However, the primary reason for sampling such an broad range of Δχ values for either Δχsinus-tissue and Δχmastoid-tissue is that, it is well known that the sinuses (which include frontal, ethmoid, sphenoid, maxillary, and nasal cavity) predominantly contain air but the mastoid region is known to have both bone and air. Consequently, it is reasonable to anticipate that the mastoid cavity might have a different Δχ value with respect to brain tissue compared to Δχ between sinuses and brain tissue. The simulated phase values were complex divided into the corresponding 1ms or 2ms effective TE phase data to generate 15 sets of residual phase maps for each volunteer. The phase in the immediate vicinity of an air-space (sinus or the mastoid cavity) is predominantly influenced by the susceptibility difference between that and the brain tissue. Under this assumption, for example, in the case of quantifying Δχmastoid-tissue, only voxels around the mastoid cavity were used for the analysis. Hence the σresidual_phase from two independent volumes of interest (VOI), one chosen in close proximity to the sinuses and the other chosen in to the neighborhood of the mastoid cavity, were evaluated and plotted as a function of the corresponding Δχsinus-tissue and Δχmastoid-tissue, respectively. On an average each of the two VOIs contained roughly 6000 voxels. The susceptibility values corresponding to the minimum σresidual_phase were the desired values of Δχsinus-tissue, and Δχmastoid-tissue. The susceptibility of the mastoid cavity χmastoid is Δχsinus-tissue – Δχmastoid-tissue.

Modified SWI Phase processing

For each volunteer, the phase at TE=40ms was calculated using Eq. 4 and the estimated susceptibilities. This estimated phase was then subtracted from the original SWI phase images through complex division to give the GDAC-phase images. (In this paper, phase subtraction is always done through complex division to avoid effects from phase aliasing. Hence phase subtraction and complex division are interchangeably used unless specified otherwise). Both the GDAC phase and the original SWI phase images were subsequently used for generating the corresponding post-processed susceptibility weighted (pSWI) magnitude images: GDAC-pSWI and the conventional-pSWI respectively. Their results were qualitatively compared for the amount of the recovered brain signal.

It has been shown (12) that unwrapping the original phase images and subsequent high pass filtering can provide better HP filtered phase images and consequently produce better pSWI magnitude images. Similar to the original phase images, GDAC-phase images can also be used as input for any of the existing phase unwrapping algorithms. Since most of the rapid aliasing near air-tissues interface is already removed, using GDAC-phase should provide an improved result in such processing as well. To validate this idea, a simple two-dimensional least-squares phase unwrapping algorithm (21) was applied on both the original phase and the GDAC-phase images. The results were subsequently HP filtered to generate the corresponding pSWI images: Unwrap-HP-original-pSWI and Unwrap-HP-GDAC-pSWI, respectively. The results were also qualitatively compared. A flow chart delineating the steps involved in the new processing is shown in Figure 1.

Figure 1
Flow chart of steps involved in the new SWI processing. pSWI refers to post processed susceptibility weighted magnitude images. The resultant pSWI images obtained using HP filtered original phase (conventional-pSWI) and those obtained using HP filtered ...


Results of Simulations

A volume of interest of roughly 7000 voxels around the simulated sphere was used to mimic roughly how many voxels would be available in human studies near the sinuses or cavities. The results of the σresidual_phase minimization, within this volume of interest, for finding the susceptibility difference between the shell and surrounding region, are shown in Fig. 2. The Δχ value corresponding to the minimum σresidual_phase value is at −9.8 ppm. The precision of this result is limited by the step size of Δχ used in the minimization process, which was 0.1ppm. This simulation does not consider the effect of noise. However, since noise is additive in nature, it will only influence the minimum attainable σresidual_phase value and does not affect the Δχ at which this minima occurs. Figure 2 also shows the result from a subset of the former VOI, containing only 1024 voxels to point out the reduced chi squared error achieved from a smaller VOI. The result however, is the same.

Figure 2
Simulation Results: A plot of σresidual_phase as a function of Δχ values (in ppm). Results of analysis from two volumes of interest are plotted. The inlet shows a zoomed version of the plot between −8.9 ppm and −11.1ppm. ...

Phantom Results

The point of examining the geometry at two different echo times is to understand how much signal might be lost at the edges of the sinuses and tissues and how this affects the phase and susceptibility estimations and eventually the spurious phase removal process. Figure 3 shows transverse and coronal views of the phantom images at TE = 6.58ms. Fig. 3e and 3f show the images of the extracted geometry from two different magnitude images (TE 6.58ms and 40ms respectively) and their difference images (3g and 3h). There is a loss of signal at the edge of the sphere in the SWI data which translates to an increase in volume of the sphere when extracted from the SWI magnitude images.

Figure 3
Magnitude and phase images of the transverse (a, b) and coronal (c, d) views of the gel phantom at TE 6.58ms. Corresponding transverse view of the extracted geometries from the TE 6.58ms data set (Ggre) (e) and from the TE 40ms SWI magnitude images (G ...

Figure 4a shows the plots of σresidual_phase as a function of the Δχ value. The Δχair_gel corresponding to the minimum σresidual_phase for Ggre is at -9ppm, whereas for Gswi it is at −6ppm. This difference clearly shows the importance of the definition of the object geometry which will subsequently affect the effective Δχ value to be used and the ensuing phase calculations. Figures 4c and 4d show the two GDAC phase images obtained by subtracting the predicted phase using Δχair_gel=− 9ppm with Ggre and Δχair_gel=−6ppm with Gswi respectively from the original SWI phase images. Although the result from Δχair_gel = −6ppm still contains some remnant aliasing, the result from Δχ = −9ppm has removed almost all the aliasing in the vicinity of the sphere arising due to air-gel interface. The corresponding HP filtered phase images (f, g and h) demonstrate the dramatic improvement in the GDAC phase images in terms of the remnant phase aliasing. Moreover, Fig. 4i shows that even a mild filter can be employed to remove most of the remaining rapid phase aliasing. Figure 5 illustrates the effect of using the GDAC phase as opposed to the original phase in generating pSWI images.

Figure 4
a) σresidual_phase vs Δχ (in units of ppm) plots for Ggre and Gswi showing that the minimum errors of Δχair_gel, are at −6ppm for Gswi and −9ppm for Ggre; b) an original SWI phase image (TE 40ms); ...
Figure 5
Comparison of the minimum intensity projections over 4 slices of; a) original SWI magnitude; b) post-processed SWI magnitude (pSWI) using HP filtered original phase images with a filter size of 64×64; c) post-processed SWI magnitude using the ...

The strong aliasing from the original SWI phase images causes terrible phase artifacts in the HP filtered images which lead to a concomitant artifact in the pSWI images (see Fig. 5b). However, the pSWI images generated from the GDAC phase (using Δχair_gel=−9ppm with Ggre) is devoid of such artifacts (Fig 5c). Nevertheless, errors in the definition of the object at the boundary at different locations (typically single voxel errors) cause erroneous local dipolar field and hence erroneous additional local phase at these locations. When the predicted phase is subtracted from SWI original phase images, this extra phase gets introduced in the resultant GDAC phase images. This eventually leads to erroneous “signal loss artifact” on the boundary of the object when these HP filtered GDAC phase images are used to generate pSWI images (see edges of the phantom in Fig 5c). Conversely, there are no such erroneous local dipolar phase in the original SWI phase images or in the HP filtered original phase images. Combining the information in the HP filtered original phase images and the HP filtered GDAC phase images, by choosing minimum absolute phase between the two, one can generate a better HP filtered phase image and consequently an improved susceptibility weighted magnitude image, as shown in Figure 5d.

Human Imaging Results

Figure 6 shows the plots of σresidual_phase versus Δχ values for each of the three volunteers. The Δχsinus-tissue value corresponding to the minimum σresidual_phase is at −11ppm for volunteer 1, −13ppm for volunteer 2 and −7ppm for volunteer 3. The Δχmastoid-tissue values are −7ppm for volunteer 1, −6ppm for volunteer 2 and −5ppm in volunteer 3. Figure 6d illustrates the regions where the volumes of interest were chosen. The results show that Δχsinus-tissue is considerably different from Δχmastoid-tissue value in all three volunteers. Furthermore these values also change from volunteer to volunteer. We enlarged the FOV of volunteer 3 by a factor of 1.25 to verify if the difference between Δχsinus-tissue and Δχmastoid-tissue was from an error in field calculation due large object-size/FOV ratio. But the Δχ values were unaffected indicating that this was not a result of an error in field calculation.

Figure 6
σresidual_phase vs Δχ (in ppm) for volunteers 1 (a), 2 (b) and 3 (c). Values for both Δχsinus-tissue (>) and Δχmastoid-tissue (*) are plotted. The value of Δχ corresponding ...

The quantified Δχsinus-tissue and Δχmastoid-tissue are used to predict the phase values at 40ms echo time. Figure 7b shows the phase calculated at 40ms in volunteer 2, using the estimated Δχ values. Almost all the geometry induced background phase and its associated aliasing are removed in the GDAC phase images, generated by complex dividing SWI phase images by the calculated phase images. Figure 7e shows the improvement in the HP filtered GDAC phase image when compared to the HP filtered original phase. This translates to a significant improvement in the corresponding pSWI images (Fig. 7h) and makes it possible to properly process the frontal and midbrain regions and visualize the veins in those regions better (see Fig. 8). Figure 8 shows the improvement in the visualization of vessels in the frontal region of the brain, which were heretofore obscured by the conventional SWI processing. Despite a possible error in estimating Δχ values, the results after phase subtraction are devoid of most of the rapid aliasing near the air-tissue interface. Further high pass filtering removes any remnant low spatial frequency aliasing, giving much smoother phase images and consequently better pSWI images. Figure 9 shows results from all three volunteers. The images are taken from 3 different levels of the head along the cranio-caudal axis to show that considerable improvement in the image quality is obtained in all volunteers and in multiple regions of the brain using HP filtered GDAC phase images. The result of using phase-unwrapping first and then HP filtering for generating pSWI images is shown in Figure 10. There is improvement in the image quality when the GDAC phase images are used (Fig. 10b), compared to when the original phase images are used (Fig. 10a). The histogram plots, Figure 10d, of the regions shown in images a and b show the number of voxels recovered in the Unwrap-HP-GDAC-pSWI image compared to the Unwrap-HP-original-pSWI image.

Figure 7
a) Original SWI phase image; b) predicted phase at 40ms from the head geometry using the quantified Δχ values; c) result of subtracting (b) from (a) through complex division; d) result of HP filtering (a); e) result of HP filtering (c); ...
Figure 8
Conventional pSWI images(a and c) and corresponding GDAC-pSWI images (b and d). The venous vessels, otherwise obscured by the conventional SWI processing, are now clearly visible using the modified SWI post-processing. The arrow heads point to the regions ...
Figure 9
Conventional pSWI images (a,d and g); Corresponding GDAC-pSWI images (b,e and f); Difference images (c, f and i) of (b - a), (e - d) and (h - g) respectively. Images a through c are from volunteer 2; d through fare from volunteer 3 and g through i are ...
Figure 10
a) Unwrap-HP-original-pSWI; b) Unwrap-HP-GDAC-pSWI; c) difference image (b - c); d) histogram plots of pixels taken from the ROIs shown in a and b. The histogram curve of Unwrap-HP-original-pSWI shows a bi-modal distribution signifying the erroneous signal ...


Susceptibility weighted imaging has found its major impact in being able to present information about local changes in susceptibility. This can come from deoxyhemoglobin when imaging the veins, from hemosiderin when imaging microbleeds or from measuring local iron content in multiple sclerosis (2, 5, 22). The ability to remove the unwanted background fields, especially the deleterious effects of geometry dependent field effects from air-tissue interfaces, has relied on the use of a high pass filter approach that assumes that all background phase errors are low spatial frequency in nature. Although the results to date have been reasonable in some parts of the brain, there have been problems in the midbrain and near the sinuses. Removing rapid phase aliasing near sinuses using the geometry information can facilitate the use of a milder filter leading to better preserved contrast and more accurate local phase information. By removing this geometry dependent phase or reducing its effects, this method opens the possibility for SWI to become a viable technique in other areas in the body as well, such as: in the spine, in the breast, in the knee and perhaps even in the heart.

Previous works using the Fourier based method (15) have estimated the static field deviation due to the air-tissue interface in brain by obtaining the geometry of the head from computed tomography (CT) images and assuming a susceptibility of −9.2 ppm between the brain tissue and air and a susceptibility of −11.4 ppm between the bone and air. However, it is not always possible to obtain CT images of the subject being imaged, nor is the MR image registration with a fixed CT template always feasible. In this paper, we used the T1 weighted MR magnitude images to obtain the head geometry. However, one must be aware that the shape and the size of the extracted structures will be a function of echo time and bandwidth of the MR sequence used. For example the frontal sinus is about 1.5 cm along one dimension but can easily appear to be 1.7 cm on a gradient echo image due to the T2* signal loss. This overestimate of the sinus size can change its estimated susceptibility with respect to the brain tissue by almost 40%, i.e. a smaller susceptibility value (smaller than its actual physical value) is to be used with this increased sinus size to obtain the same phase value outside the sinus (i.e., in the brain tissue). Thus an error in the extracted geometry size in-turn affects the effective susceptibility value to be used. Furthermore, the internal structure of the extracted geometries although not included (either due to T2* loss or the presence of bone) do physically influence their effective Δχ value. Because of these factors, Δχ values were evaluated for each volunteer. There can also be small “nibbles” or localized signal losses at the tissue-sinus interfaces. In practice, we smooth over these to have a homogeneous tissue or sinus region. Also, in principle, the two short echo time datasets collected to obtain the unaliased phase maps can be cut down to one data set and its phase obtained by unwrapping(23) the image.

As pointed out previously (14, 15, 19), the accuracy of the FT based field estimation method very much depends on the ratio of the object size to the FOV. In our simulations, although the outer diameter of the shell geometry is almost 85% of the field of view, the result of susceptibility quantification is within 5% of the actual value. The primary reason is because the voxels chosen for analysis are close to the inner sphere, such that we mainly quantify the susceptibility difference between the inner sphere and the shell. Since the diameter of the inner sphere is roughly 15% of the total field of view, the difference between the calculated susceptibility and the actual susceptibility is expected to be small. Nevertheless, the calculated susceptibility, from Fig. 2, is found to be at −9.8ppm instead of −9.5ppm. This difference of about 3% could be partly due to discretization and partly due to the possible small field error in the voxels chosen.

The results from three human volunteers indicate that Δχsinus-tissue differs from Δχmastoid-tissue. It is to be noted that the values of Δχ estimated are specific to the voxel aspect ratio of the geometry, which is in this case 1:1:2, and are not to be taken as absolute Δχ values The variability in Δχmastoid-tissue and Δχsinus-tissue between volunteers is not surprising as these values can depend on the actual internal structures of the nasal cavity and the mastoid cavity, as septal cartilage distribution and distribution of air cells and bone are not included in the geometry defined by MRI. Variations in magnetic field shimming prior to imaging these volunteers might also be one of the factors influencing the inter-subject variability of the Δχ values. The Fourier based method does not include the additional field effects due to shimming and eddy currents. The different acquisition parameters used for different volunteers could also in principal, influence the inter-subject variability of Δχ values, but any such effects are expected to be very small.

An important advantage of this method is that it can be used as a precursor to existing methods of phase unwrapping or filtering to further improve these methods (see, for example, Figure 10). The time limiting steps in generating GDAC phase images are obtaining the geometry of the air-spaces and estimation of Δχ values which need considerable manual intervention. Currently, these steps take about 2 to 3 hours. Once these steps are done, generation of GDAC phase images for a dataset of size 256×256×128 takes about 3 minutes of computation in Matlab on a system with 64 bit AMD Turion processor with 2GB RAM memory. As these steps become automated, the total processing time should drop dramatically.

In conclusion, we have demonstrated that the Fourier transform based field estimation method can be successfully employed for removing the unwanted phase from air-tissue interfaces caused by the structure of the brain and air-spaces within, for susceptibility weighted imaging. This makes it possible for better evaluation of veins and quantification of iron in the midbrain and the frontal brain and can be of particular value in trauma cases where microhemorrhages and venous damage can occur in these areas. Overall, it is one more step to removing artifacts in susceptibility weighted imaging to make it more robust for clinical applications.


Grant Support: This work was supported in part by grants from the State of Michigan Grant #085P5200251, National Institutes of Health Grants #AG20948 and #2R01 HL062983-04A2.


1. Reichenbach JR, Venkatesan R, Schillinger DJ, Kido DK, Haacke EM. Small vessels in the human brain: MR venography with deoxyhemoglobin as an intrinsic contrast agent. Radiology. 1997;204:272–277. [PubMed]
2. Sehgal V, Delproposto Z, Haddar D, et al. Susceptibility-weighted imaging to visualize blood products and improve tumor contrast in the study of brain masses. J Magn Reson Imaging. 2006;24:41–51. [PubMed]
3. Sehgal V, Delproposto Z, Haacke EM, et al. Clinical applications of neuroimaging with susceptibility-weighted imaging. J Magn Reson Imaging. 2005;22:439–450. [PubMed]
4. Wycliffe ND, Choe J, Holshouser B, Oyoyo UE, Haacke EM, DK K. Reliability in detection of hemorrhage in acute stroke by a new three-dimensional gradient recalled echo susceptibility-weighted imaging technique compared to computed tomography: a retrospective study. J Magn Reson Imaging. 2004;20:372–377. [PubMed]
5. Haacke EM, Makki MI, Selvan M, et al. Susceptibility Weighted Imaging Reveals Unique Information in Multiple-Sclerosis Lesions Using High-Field MRI. Proceedings of 15th annual ISMRM; Berlin. 2007. p. #2302.
6. Raz N, Rodrigue KM, Haacke EM. Brain aging and its modifiers: insights from in vivo neuromorphometry and susceptibility weighted imaging. Annals of the New York Academy of Sciences. 2007;1097:84–93. [PMC free article] [PubMed]
7. Tong KA, Ashwal S, Obenaus A, Nickerson JP, Kido DK, Haacke EM. Susceptibility-Weighted MR Imaging: A Review of Clinical Applications in Children. Am J Neuroradiol. 2008;29:9–17. [PubMed]
8. Haacke EM, Cheng NY, House MJ, et al. Imaging iron stores in the brain using magnetic resonance imaging. Magn Reson Imaging. 2005;23:1–25. [PubMed]
9. Wang Y, Yu Y, Li D, Bae KT, Brown JJ, Lin W, Haacke EM. Artery and vein separation using susceptibility-dependent phase in contrast-enhanced MRA. J Magn Reson Imaging. 2000;12:661–670. [PubMed]
10. Haacke EM, Ayaz M, Khan A, et al. Establishing a baseline phase behavior in magnetic resonance imaging to determine normal vs. abnormal iron content in the brain. J Magn Reson Imaging. 2007;26:256–264. [PubMed]
11. Yao B, van Gelderen P, de Zwart JA, Duyn JH. Brain Iron in MR imaging: R2* and Phase Shift at Different Field Strengths. International Society of Magnetic Resonance in Medicine; Proceedings of 15th annual ISMRM; Berlin. 2007. p. #2165.
12. Rauscher A, Barth M, Herrmann KH, Witoszynskyj S, Deistung A, Reichenbach JR. Improved elimination of phase effects from background field inhomogeneities for susceptibility weighted imaging at high magnetic field strengths. Magn Reson Imaging. 2008 PMID: 18524525. [PubMed]
13. Salomir R, Senneville BD, Moonen CT. A Fast Calculation Method for Magnetic Field In homogeneity due to an Arbitrary Distribution of Bulk Susceptibility. Concepts Magn Reson Part B (Magn Reson Engineering) 2003;19B:26–34.
14. Marques JP, Bowtell R. Application of Fourier-Based Method for Rapid Calculation of Field Inhomogeniety Due to Spatial Variation of Magnetic Susceptibility. Concepts in Magnetic Resonance, Part B (Magnetic Resonance Engineering) 2005;25B:65–78.
15. Koch KM, Papademetris X, Rothman DL, de Graaf RA. Rapid calculations of susceptibility-induced magnetostatic field perturbations for in vivo magnetic resonance. Phys Med Biol. 2006;51:6381–6402. [PubMed]
16. Deville G, Bernier M, Delrieux J. NMR multiple echoes observed in solid 3He. Physical Review B. 1979;19:5666–5688.
17. Enss T, Ahn S, Warren WS. Visualizing the dipolar field in solution NMR and MR imaging: three-dimensional structure simulations. Chemical Physics Letters. 1999;305:101–108.
18. Pandian DS, Ciulla C, Haacke EM, Jiang J, Ayaz M. Complex threshold method for identifying pixels that contain predominantly noise in magnetic resonance images. J Magn Reson Imaging. 2008;28:727–735. [PubMed]
19. Neelavalli J, Cheng YCN, Haacke EM. Improved Fourier Based Method for Calculating Field Inhomogeniety from Known Susceptibility Distribution. Proceedings of the 15th Annual Meeting of ISMRM; 2007. p. #1016.
20. Haacke EM, Xu Y, Cheng Y-CN, Reichenbach JR. Susceptibility weighted imaging (SWI) Magn Reson Med. 2004;52:612–618. [PubMed]
21. Ghiglia DC, Romero LA. Robust 2-Dimensional Weighted and Unweighted Phase Unwrapping That Uses Fast Transforms and Iterative Methods. Journal of the Optical Society of America a-Optics Image Science and Vision. 1994;11:107–117.
22. Haacke EM, DelProposto Z, Chaturvedi S, Sehgal V, Tenzer M, Neelavalli J, Kido D. Imaging Cerebral Amyloid Angiopathy with Susceptibility Weighted Imaging: A Case Report. Am J Neuroradiol. 2007;28:316–317. [PubMed]
23. Rauscher A, Barth M, Reichenbach JR, Stollberger R, Moser E. Automated unwrapping of MR Phase Images Applied to BOLD MR-Venography at 3 Tesla. J Magn Reson Imaging. 2003;18:175–180. [PubMed]