Background field removal based on local complex phase unwrapping and spherical mean value property

Introduction Quantitative Susceptibility Mapping (QSM) is valuable technique for quantification of iron or calcification content, and venous oxygen saturation (1-4). The quality of QSM relies on the accuracy in background field removal, since the original phase contains the component from the background field which is mainly induced by the air-tissue interfaces. Recently the SHARP method was proposed, which removes the background field based on the spherical mean value property of the harmonic functions (4). Phase unwrapping is generally required, which is usually a time-consuming procedure. Although phase unwrapping could be avoided by using the Laplacian of the field which could be calculated directly from the original phase (5, 6), the Laplacian based method has errors in regions with sharp phase changes, such as the edges of the veins (5, 6). Besides, the Laplacian of the field is equivalent to a spherical filter with radius of 1 pixel, may lead to reduced accuracy in the SHARP processed phase images (6, 7). In this study, we proposed a method which allows for using an optimal kernel size other than 1 pixel in SHARP without explicit phase unwrapping. This helps to reduce the processing time and improve the accuracy in background field removal using SHARP. Theories and Methods To remove the background field using SHARP, the difference (B’) between the field B and its local mean value Bm in a spherical VOI has to be calculated. This is a convolution process using a normalized spherical kernel, f: B’=B-Bm=B-B*f. Since B=Blocal+Bb, where Bb is the background field, and Bb=Bb*f. Thus, B’=Blocal-Blocal*f. The local field can be calculated from B’ through a deconvolution using the kernel (δ-f). In the above calculations, the unwrapped phase P was used, and B’=P’/(-γB0TE). For a given pixel, P’=P-P*f=∑[P(i,j,k)-Ps]/N, where Ps denotes for the pixels in the spherical VOI centered at P(i,j,k). Fundamentally, the difference between the phase of every pixel and the center pixel in the spherical VOI (P’) needs to be calculated. Assuming that |P(i,j,k)-Ps|<π, and there is only one wrap within the spherical region, this phase difference can be calculated through complex division: P(i,j,k)-Ps=arctan(exp(1i*(P-Ps))). This suggests that it is possible to calculate this phase difference P’ using the original phase Pw. Since the phase images with wraps Pw=P+2πn, so∑[Pw(i,j,k)-Pw,s]/N equals P’, if n(i,j,k)=ns(i,j,k), or P’±2πm/N when n(i,j,k)=ns(i,j,k)±1. ns(i,j,k) denotes for the number of wraps for a certain pixel inside the spherical region. We only need to determine “m”, the number of pixels where n(i,j,k)=ns(i,j,k)±1. This is done by counting the number of pixels where |Ps(i,j,k)-P(i,j,k)|>π. Consequently, an initial estimate of P’ could be calculated using the original phase images Pw through convolution and this estimate of P’ will be corrected using m. Then the local field could be recovered from P’ through a deconvolution as it is done in SHARP. The proposed algorithm was tested on a 3D brain model, which has several basal ganglia structures, veins, grey matter, white matter and sinuses, with realistic susceptibility values. The air-tissue interface induced background field was simulated by assigning 9ppm to the sinuses. The field of this 3D brain model was calculated using the fast forward calculation. Phase images were generated at TE=13ms. Further tests were carried out on one in vivo dataset, which was collected with 0.5 isotropic resolution on a 3T SIEMENS VERIO scanner using 3D gradient echo sequence, with matrix size 512x384x224, TE=15ms. The in vivo phase images were first unwrapped using Phun (8), when the traditional SHARP was used; while the proposed algorithm was directly applied on the original phase images with wraps. Same spherical kernel and regularization parameter were used in both the proposed and traditional SHARP. The kernel sizes were chosen to be 4 pixels and 6 pixels for 3D brain model and in vivo data respectively; while the regularization parameters were chosen to be 0.01 and 0.05 for the model and in vivo data respectively. Susceptibility maps were generated using the phase images processed with the proposed algorithm, for the in vivo data only, using the truncated kspace division with a truncation threshold of 0.1.