Effect of object orientation angle on t2* image and reconstructed magnetic susceptibility: numerical simulations.

The magnetic field resulting from material magnetization in magnetic resonance imaging (MRI) has an object orientation effect, which produces an orientation dependence for acquired T2* images. On one hand, the orientation effect can be exploited for object anisotropy investigation (via multi-angle imaging); on the other hand, it is desirable to remove the orientation dependence using magnetic susceptibility reconstruction. In this report, we design a stick-star digital phantom to simulate multiple orientations of a stick-like object and use it to conduct various numerical simulations. Our simulations show that the object orientation effect is not propagated to the reconstructed magnetic susceptibility distribution. This suggests that accurate susceptibility reconstruction methods should be largely orientation independent.

Introduction

T2* weighted magnetic resonance imaging (T2*MRI) is designed to image an internal inhomogeneous fieldmap that is usually formed from the magnetic susceptibility property of an object (typically biological tissue in the medical imaging realm) via material magnetization in a main field.1 For biological tissue imaging (like brain imaging), the underlying source of T2*MRI is attributed to the tissue’s magnetic susceptibility property.2,3 For brain imaging with T2* MRI, the vascular structure in gray matter tissue reveals orientation dependence in the interim fieldmap and in the output T2* images.4–7 On one hand, the object orientation effect (also called the angular effect) that is inherent with the T2*MRI has been used for tissue anisotropy investigation via multi-angle imaging6–13 (by rotating the object orientation so as to change the orientation angle with respect to the main field direction). On the other hand, the angular effect is undesirable in the magnetic susceptibility reconstruction from the T2* images.4,5,14–16 The challenge of susceptibility reconstruction involves the inverse problem of T2*MRI, as described by computed inverse MRI (CIMRI).15 It has been experimentally shown16 that the magnetic susceptibility distribution in an agar cylinder can be consistently reconstructed from a T2* phase image acquired at different orientation angles. However, understanding of the orientation effect on T2*MRI data acquisition and quantitative susceptibility mapping remains elusive. In this report, we design a stick-star phantom and perform numerical simulations on T2*MRI and CIMRI. We thereby demonstrate the effect of object orientation on forward T2* image acquisition and backward susceptibility source reconstruction.

From the viewpoint of imaging system, a T2*MRI procedure produces a T2* image (output) from a susceptibility-expressed source (input) via the object magnetization that provides an inhomogeneous fieldmap (interim). The object magnetization can be expressed by a 3D convolution with an anisotropic bipolar-valued dipole kernel.15,17 In the Fourier domain, the 3D convolution corresponds to a 3D spatial filtering (with a filter corresponding to the Fourier transform of the magnetic dipole field), which expedites understanding of the element-wise data entry multiplication of 3D filtering. Due to the zero surface embedded in the 3D filter,15 the 3D filtering involves a “multiply-by-zero” data degeneracy, and subsequently the 3D inverse filtering involves a “divide-by-zero” singularity. In this report, we will elucidate the 3D filtering aspects of two CIMRI methods for magnetic susceptibility reconstruction. The filter truncation solver (inverse filtering with a truncated inverse filter) involves the regularization of “divide-by-zero” singularity, and the total variation (TV) iteration solver suffers from data degeneracy of “multiply-by-zero”.

Theory

Complex T2* image formation

The T2*MRI is designed to sense an inhomogeneous fieldmap that is largely caused by an inhomogeneous susceptibility source via a magnetization in a main field for tissue imaging. Let B0 denote the main field, χ(r) the susceptibility distribution of tissue structure, and b(r) the susceptibility-induced fieldmap (specifically, the z-component of the vector field), then the fieldmap establishment via susceptibility magnetization is expressed by17–19

where * denotes 3D spatial convolution, and h(r) represents the point magnetic dipole field.20 The asymmetry about x-, y-, z-axis indicates the anisotropy of the kernel. Exposed to the fieldmap b(r), the proton spins (via the host of water molecules) undergo Larmor precessions. The intravoxel spin precession dephasing average signal of voxel at r, called T2* voxel signal, is expressed by

where Ω(r) denotes a voxel at r = (x,y,z) with the voxel size of |Ω|, B0 the main field, T the echo time, γ the proton gyromagnetic ratio. The spatial tessellation of T2* voxel signal (via a field-gradient encoding scheme of MRI) makes a multivoxel T2* image. Upon the acquisition of the complex-valued dataset C(r;TE), we can extract its magnitude image and phase image by

where we explicitly keep the T dependence in T2* image. It is noted that A(r) represents the magnitude loss of the complex signals and P(r) the phase angle accrual. We will show that both A(r) and P(r) are dependent upon the structural orientation of tissue susceptibility distribution with respect to the main field direction.

CIMRI

The magnetic susceptibility reconstruction is implemented by reversing the T2*MRI procedure by computer, as termed CIMRI and notated by

where “recon” connotes reconstruction. For simplifying the fieldmap calculation, we assume that the T2* image is acquired under a small phase angle condition15,18 (where the phase angle is far from wrapped), which leads to a linear relationship between the fieldmap and the T2* phase image (different by a constant factor). Upon the fieldmap calculation, CIMRI is implemented by solving the ill-posed inverse problem of 3D convolution in Equation (1). In general, the methods for CIMRI implementation can be classified into three categories:15 matrix inverse, inverse filtering, and iteration. In this report, we will compare two CIMRI solvers for magnetic susceptibility source reconstruction: filter truncation and TV iteration.

Filter truncation solver

In the Fourier domain, the magnetic susceptibility reconstruction can be achieved by inverse filtering with a truncated filter15

where k=(kx,ky,kz), sgn(t) denotes a sign function (sgn(t) = 1 for t > 0 and = −1 for t > 0), and ɛ0 is a small positive constant (say 0.01), called a truncation threshold. The tilde denotes Fourier transform and “trunc” indicates truncation by thresholding. The asymmetry of h̃ (k) about kx-, ky-, and kz-axes indicates the 3D filter anisotropy (in Fourier domain). The filter truncation in Equation (5) is used to regularize the singularity associated with the data entries within the vicinity of the zero surface cone (defined by h̃(k)= 0). It is known that any data entry alteration on the inverse filter will result in energy non-conservation after an inverse filtering.15 The filter truncation in Equation (5) involves a selection of the threshold ɛ0, which delineates a truncated cone zone (thickness > 1 voxel, ɛ0-dependent) that encompasses the zero surface cone (thickness = 1 voxel).

Total variation regularization

The TV-regularized iteration method for magnetic susceptibility reconstruction15 is expressed by

where ||χ|| denotes a TV norm and λ the regularization parameter. The magnetic susceptibility reconstruction by TV iteration can be efficiently implemented by split Bregman iteration algorithms.15,21

Goodness of magnetic susceptibility reconstruction

For numerical simulation, we can predefine a source truth for χ, and then calculate its reconstruction version χ from a T2* phase image by CIMRI. In order to numerically characterize the goodness of susceptibility reconstruction in reference to the predefined truth, we suggest a measure of 3D spatial correlation by

where corr designates a operation for finding the correlation coefficient.

Numerical Simulations

In order to demonstrate the angular dependence in T2*MRI, we design a stick-star digital phantom in Figure 1 for numerical simulations. The stick-star pattern consists of five straight sticks at polar angles (with respect to the main field direction) of {0, ±27.4°, ±54.7°, ±82.2°, ±90°}, with one stick deliberately specified at an orientation angle of 54.7° (which represents the magic angle as determined by cos2θ = 1/3). The geometry of the stick-star phantom is numerically characterized by a 3D binary volume, denoted by χ(x,y,z), which assumes 1 ppm susceptibility values in the stick interior and 0 ppm values in the stick exterior. All the stick cylinders assume the same susceptibility values (1 ppm). We consider the star-stick phantom as a collection of rotated vascular vessels (or fiber axons, tendons) at five selected orientation angles, thereby we can observe the angle effect through a profile across the sticks along a half-circle scanline (as indicated in Figure 1(a1)).

Figure 1

A stick-star digital phantom, which assumes a binary volume of magnetic susceptibility distribution: 1 ppm inside stick interior with a 0 ppm background. The labels ‘x’, ‘y’, and ‘z’ indicate 3D axial directions. The sticks are placed at (y,z)-plane, with a polar angles of {0, ±27.4°, ±54.7°, ±82.2°, ±90°}; of which 54.7° is determined by cos2(θ) = 1/3.

In our numerical simulation, we add an additive Gaussian to the fieldmap calculation, that is

In implementation, we select NoiseLevel = 0.02. The fieldmap calculated by Equation (8) in object space is shown in show in Figure 2, in which the stick orientation effect is depicted by profile along a half-circle scanline, as plotted in Figure 2(b): the field value varies with the stick orientation angle. In particular, we observe noticeable field value drops at angles 54.7°, which are a result of “magic angle effect”.22 In addition, we also notice that the field value profiles across the sticks at large polar angles are non uniform and demonstrate edge enhancement, which may be explained by the spatial derivative property of the 3D convolution in Equation (8), as well as by a slight contribution of long-distance magnetization influence from other sticks.17,23

Figure 2

(a1–a3) Three principal cross sections of the 3D fieldmap b(x,y,z) resulting from the stick-star phantom magnetization (calculated by 3D convolution in Equation (8) with NoiseLevel = 0.02 and B0 = 1T). The susceptibility distribution along the half-circle scanline (indicated in (a1)) across the sticks is shown in (b), which shows that the field value is a function of stick orientation angle. The magic angle effect manifests as a field value drop at the polar angles (±54.7° with respect to z-axis).

The complex-valued T2* image is calculated from the fieldmap by Equation (2) by implementing a Bloch simulation algorithm24 (a liner approximation of intravoxel dephasing). Of the complex image, the magnitude and phase images are calculated by Equation (3), with the results shown in Figure 3. Again, we observe that the stick orientation effect manifests as an image value change with respect to the stick orientation angle: the height of waveform varies with the orientation angle. It is noted that the T2* magnitude image reveals strong non-negative edges, which are due to the following aspects: (1) the inheritance of the edge effect in fieldmap, (2) the linear inter-voxel difference approximation of the Bloch simulation algorithm,24 and (3) the inherent negativity and edge effect of intravoxel dephasing magnitude signal.18 Meanwhile, the T2* phase image conforms very well to the fieldmap (except for a scale factor under the phase-unwrapped condition). Since the T2* image formation from the fieldmap only involves spatial partition (voxelization) and intravoxel spin signal average, we assume that there is no cause for creating an orientation effect in the intravoxel dephasing signal formation stage. In other words, the angular effect only occurs in the fieldmap formation during magnetization, and it propagates to the T2* image during intravoxel dephasing signal formation. We should mention that the effect of tissue structure orientation on the MR signals has been widely used for tissue structure imaging, in particular for the structural study of neuronal fibers7,9 and for tendon and cartilage.10,13,25 Nevertheless, in this report, we are instead concerned with how the orientation effect affects the magnetic susceptibility reconstruction. A good magnetic susceptibility reconstruction is expected to stop the orientation effect propagating from the T2* image to the reconstructed source, or remove the orientation effect from source reconstruction.

Figure 3

T2* magnitude image A(x,y,z) and phase image P(x,y,z) of the stick-star phantom (calculated by Equations (2) and (3)). It is noted that the magnitude image suffers from an edge effect and non-negativity, and both the magnitude and phase show orientation dependence.

In a small phase angle regime, the MR phase image is linearly related to the fieldmap by a scale factor,15 implying that there is no information loss between the fieldmap and the MR phase image. From the fieldmap, we can reconstruct the susceptibility phantom (the phantom geometry expressed in terms of magnetic susceptibility property) by two CIMRI solvers: filter truncation and TV iteration. The reconstructed susceptibility maps are shown in Figure 4. By visual comparison with the predefined susceptibility truth in Figure 1, we can see that both the filter truncation and the TV iteration solvers can reconstruct the susceptibility phantom in which there is no noticeable residual of the orientation effect. It is noted that the filter truncation method suffers from heavy background noise15 (with a global uniform randomness), and the TV iteration method suffers from a sparse noise (“comeback” noise).26–28 It is mentioned that the stick with an orientation at the magic angle (54.7°) is reconstructed as accurately as other oriented sticks, implying that magnetic susceptibility reconstruction is free from the object orientation effect (except for a very small information loss, as is demonstrated later in this paper).

Figure 4

Reconstructed magnetic susceptibility distributions of the stick-star phantom.

Notes: Upper row: reconstruction by filter truncation solver (by Equation(5)); lower row: reconstruction by TV iteration solver (by Equation (6)). It is seen that the reconstructed susceptibility distribution is free from both the orientation effect and the edge effect.

In order to scrutinize the orientation effect, we provide in Figure 5 the profiles along the half-circle scanlines for the predefined susceptibilty truth and the reconstructed susceptibility distribtions. In Figure 5, we numerically characterize the reconstruction goodness by ρ (defined in Equation (7)), which shows that the TV iteration solver outperforms the filter truncation solver. A good susceptibility source reconstruction is expected to produce a square waveform with equal heights along the half-circle scanline, indicating the recovery of uniform susceptibility distributions over all the stick interiors. In other words, the angle effect removal is indicated by reconstruction of the square waves with equal height at different orientation angles on the scanlinie profile across sticks at different orientation angles. It is noticed that the noise in the reconstructed stick-star pattern is of globalism, without any locality or angular depedence. That is, the effect of stick orientation angle on the magnetic susceptibility reconstruction manifests as a global reconstruction noise due to the inevitable information loss associated with the ill-posed inverse problem.

Figure 5

Numerical profiles along the half-circle scanlines in Figure 4. (a1) Predefined susceptibility truth and the reconstructed susceptibility distribution by TV iteration solver; (a2) predefined truth and reconstruction by filter truncation solver. It is noted that both TV iteration and filter truncation solvers can reconstruct the susceptibility source that is free from the orientation effect (indicated by the equal height susceptibility reconstruction in square waves) except a globally distributed reconstruction noise (residual).

Discussion

Experiments have observed the object orientation effect on T2* images: both MR magnitude and phase images of an anisotropic susceptibility distribution vary with the orientation angle with respect to the main field. We attribute the orientation effect to the interim fieldmap establishment during susceptibility magnetization. The orientation effect propagates to the T2* image during intravoxel dephasing signal formation. By reversing the T2*MRI procedure by CIMRI, the angular effect is not backwardly propagated to the reconstructed susceptibility source (a concomitant achievement of 3D deconvolution).

In the Fourier domain, the 3D convolution of susceptibility magnetization is represented in an element-wise multiplication of 3D spatial filtering, which allows us to look into data manipulations rendered by the CIMRI solvers: the filter truncation solver deals with a “divide-by-zero” problem by truncating the zero surface cone, and the TV iteration solver only involves a “multiply-by-zero” data degeneracy problem (due to “multiply-by-zero” surface filtering). As a result, the filter truncation solver suffers an information change due to the data entry alteration in the truncation cone zone, and the TV iteration solver suffers an information loss due to the data degeneracy (“multiply-by-zero” surface filtering). It is understood that the zero surface (with 1-vixel thickness) embedded in the 3D filter only occupies a small fraction number of data entries in relative to the digital Fourier space (eg, an estimate of occupancy ratio ~1E-4 in a 128 × 128 × 128 digital Fourier domain), which is much less than that of the truncation cone zone (cone thickness > 1 voxel). That is, information loss due to the suppression of the data entries on the 3D filter zero surface is unavoidable in 3D deconvolution. In comparison, the TV iteration solver is implemented at a minimal inevitable information loss due to multiply-by-zero surface filtering, and the filter truncation solver suffers more additional information change (loss) due a larger truncation zone. In principle, the “multiply-by-zero” data degeneracy problem can be solved by a multi-angle imaging scheme.4,14 In practice, the multi-angle imaging strategy is limited due to the requirement of object rotation scanning. Furthermore, the information loss due to the truncation cone zone of the filter truncation solver and that due to zero surface for the TV iteration solver are very small, such that the multi-angle imaging is not rewarded. Therefore, magnetic susceptibility reconstruction can be achieved from a T2* image acquired from a single orientation angle.5,15,29 In particular, we have shown elsewhere that the TV iteration solver can provide a very satisfactory reconstruction.15,18,30

The multiply-by-zero surface filtering implies that all data entries (spatial frequencies of the object pattern) falling on the zero surface will be completely suppressed (forced to become zeros after filtering) and remain irretrievable thereafter. The zero surface cone embedded in the 3D dipole filter assumes a cone angle that is equal to the magic angle (54.7° determined by cos2θ = 1/3). Since a line in the Fourier domain corresponds to a global textural pattern (uniform randomness) in the object space domain (globalism property of Fourier transform), it is very unlikely for a practical object pattern to generate a straight line on the zero surface cone (on a kx-kz or ky-kz plane) in the Fourier domain such that the object pattern disappears completely in the T2* image. Let us put the zero surface filtering in a different way: the zero surface in the Fourier domain corresponds to a global cluttered pattern in the object space, which is completely suppressed in the fieldmap (as a result of 3D filtering) and it is never recoverable by the 3D deconvolution anyhow. That is, the exclusion of the data entries falling on the zero surface cone plays a subtraction of global random noise pattern during the fieldmap establishment. However, the information loss due to the zero surface filtering is very small for TV iteration solver, which can be numerically characterized by an occupancy ratio (estimated by~1E-4) relative to the full information over the whole Fourier domain.

The orientation effect manifests as the field value drops close to the magic angle (54.7°) due to multiplicative filtering by small filter values approaching to zero. If not degenerated to zeros, the information transformation due to non-zero multiplications can be recovered by an inverse transformation (non-zero division). In principle, an ideal solution to a well-behaved inverse problem can completely recover the source. This is, intuitively, a number s is transformed by s/p can be perfectly recovered by (s/p) p = s provided that p≠0. For our 3D ill-posed deconvolution problem, it is the singularity of zero surface that prevents the perfect recovery of the magnetic susceptibility source. The data entries falling in the regions outside the zero surface are subject to “multiply- by-nonzero” filtering and are totally recoverable when the deconvolution solution is used. For the sake of numerical stability, the singularity of inverse filtering is regularized by delineating a truncation cone zone that is much larger than the zero surface cone. In comparison, the TV iteration solver only encounters the multiply-by-zero filtering, which is a data degeneracy problem but not a singular problem, so the TV iteration suffers minimal inevitable information loss (this is negligible, as discussed previously).

In our stick-star phantom in Figure 1, we include a vessel at a polar angle of 54.7° in object space. Our simulation shows that the information loss due to “zeroing” at the magic angle in the Fourier space does not lead to excessive reconstruction error on the 54.7° stick in comparison with other oriented sticks. That is, the inevitable information loss due to “multiply-by-zero” surface filtering is applied globally to the object pattern, not preferably to the specific stick at the magic angle in object space. Therefore, we understand that the stick at the magic angle in the object space is reconstructed as well as other oriented sticks (see Figs. 4 and 5). This observation has been demonstrated by an agar cylinder phantom experiment.16

The object orientation effect in T2*MRI has been explored for structure anisotropy study.6–9,11 In biomedical imaging realm, there has been effort on exploiting the orientation effect phenomenon for tissue anisotropy imaging, as widely reported in terms of magic angle spinning and magic angle effect.7–10,12,25 On the other hand, there have also been efforts to remove the orientation effect for obtaining orientation invariant images.31 In this report, we show that magnetic susceptibility reconstruction by CIMRI can produce a reconstructed susceptibility distribution that is free from object orientation effects. Considering the stick-star phantom as a rotated blood vessel, a rotated neuron axon, or a rotated tendon, our methodology and conclusions are applicable to all kinds of biomedical tissue structure study by T2*MRI, such as blood vasculature, neural plexus, and musculature.

Conclusion

With a stick-star digital phantom representing a magnetic susceptibility structure of cylinders at different orientation angles, we quantitatively evaluated the object orientation effect on the inhomogeneous fieldmap, the MR magnitude and phase images. Based on two CIMRI methods, ie, filter truncation and TV iteration, we show that the orientation effect is not propagated to the reconstructed susceptibility source distribution in any form of orientation-dependent pattern. By looking into the insights of T2*MRI, we clarify that the orientation effect is due to magnetic susceptibility magnetization with an anisotropic bipolar-valued point dipole field kernel, and that the CIMRI-based magnetic susceptibility reconstruction renders an inverse of the anisotropic magnetization (essentially a 3D deconvolution). Based on 3D filtering formulation of the convolutional susceptibility magnetization, we explain the origin of noise in CIMRI-based magnetic susceptibility reconstructions: the noise associated with filter truncation method is due to data entry alteration in the truncation cone zone, and that with TV iteration method is due to a “multiply-by-zero” surface data degeneracy. Since the TV iteration suffers a minimal inevitable information loss due to the “multiply-by-zero” data degeneracy (inherent to T2*MRI), it renders better reconstruction performance than the filter truncation method. As a concomitant achievement, the magnetic susceptibility reconstruction can stop the orientation effect propagation from the T2* image to the reconstructed susceptibility source. In conclusion, the CIMRI-based magnetic susceptibility reconstruction can reproduce a susceptibility source that is essentially free from object orientation effects. Nevertheless, the magnetic susceptibility reconstruction suffers from a small inevitable information loss, which introduces a global noise pattern in the form of point randomness or stripe clutter, rather than in any orientation-dependent object-relevant pattern.