Directional Filter, Local Frequency Estimate and Algebraic Inversion of Differential Equation of Psoas Major Magnetic Resonance Elastography ()
1. Introduction
For millennia, palpation has been in clinical practice for disease diagnosis by medical doctors and health workers [1]. While palpation can easily diagnose superficial lesions and tangible diseases [2], it is impossible for touch to examine deep-lying organs, unless the presenting symptoms are noticeable to other human sensations; otherwise, a biopsy must be performed. Biopsy, being an invasive procedure, fully depends on the operator’s expertise and has become an unnecessary burden to patients. With the evolution of digital diagnostic imaging modalities in radiography, fluoroscopy, computed tomography (CT), ultrasound and magnetic resonance imaging (MRI), the lesions that were once impossible to access by manual palpation have become detectable [3]. Even tiny nodules can be distinguished by cross-sectional imaging [4].
Due to limited contrast differences in the above-mentioned modalities [5], however, the final diagnosis is often uncertain. For example, the differential diagnosis of benign versus malignant lesions may still be incomprehensible. Since biomechanical parameters of tissues can be used as biomarkers for identification and classification of diseased tissues, several parameters such as elasticity, viscosity, anisotropy, density, etc. have been identified as possessing clinical significance. Yet, elasticity has garnered keen interest among medical professionals for its ability to easily differentiate healthy and diseased tissues. Given this circumstance, elastography (a term first coined by Ophir et al. [6]) has come into play with ultrasound and MRI technologies. Together, these techniques, having the benefit of being completely non-invasive, can provide quantitative stiffness values. The stiffness values range over several orders of magnitude, allowing for excellent contrast between a lesion and its surrounding normal tissues [5]. Above all, these techniques are able to quantify the sensory palpation, making it extendible virtually to deep tissues.
Magnetic resonance elastography (MRE) [7] uses harmonic mechanical excitation to measure the mechanical properties of tissues, such as shear modulus (or stiffness). The techniques for MRE application can vary for different tissues but three basic principles remain constant: 1) production of a shear wave with vibration frequency; 2) imaging of waves inside the body using a special MRI technique; and 3) processing of data and generation of an image for determination of a tissue’s stiffness. External driver devices are typically used for generating shear waves, that include electromechanical, piezoelectric-stack and pressure-activated driver systems [8]. The waves are reflected in the wavelength. The speed of these waves depends on the medium. The waves travel faster (longer wavelength) in hard tissues [9]. The wave motion produced by the driver system is measured by an MR imaging technique called phase contrast MRI [10]. Muthupillai et al. introduced the method to encode the propagating shear waves into the phase of the MR images with the help of motion-encoding gradient (MEG) pairs [7]. An MR image containing the information about propagating wave in its phase is dealt with mathematical wave inversion algorithms based on equations of motion that allow to calculate the mechanical properties, i.e. shear modulus [11]. The images of the mechanical properties of the tissues in MRE are generally called as elastograms [12]. Recently, MRE has been applied extensively to quantitatively assess the biomechanical properties of tissues in vivo; these have included brain [13] [14] [15], liver [13] [16] [17], breast [5] [13] [18], prostate [19] [20], ovary [21] [22] and skeletal muscles [5] [13] [23] [24] [25] [26]. In addition to MRE, there are other MR imaging techniques like Diffusion-weighted imaging (DWI) [27] [28] [29], Perfusion MR imaging by using exogenous contrast agents and Arterial Spin Labelling (ASL) [30] [31] [32], MR spectroscopy (MRS) [33] and hybrid technologies [34], that have revolutionized MR imaging with their own distinct role in disease diagnosis, the detail description of which are enclosed in the references.
Various physiological, pathological and genetic conditions affect the mechanical properties of skeletal muscles. Thus, MRE of skeletal muscles estimates mechanical properties of physiologic and pathologic states in quantitative terms [5] [13] [23] [24] [25] [26]. However, according to our knowledge, there is no MRE-based research relating to genetic disorders of muscles. MEG direction determines the direction of propagation of a shear wave [5] [35], so that the selection of MEG plays a crucial role in skeletal MRE. Earlier research studies have shown that shear wave displacements are induced primarily perpendicular to the muscle fibres [13] [36].
We assume that the estimation of stiffness of the psoas major (PM) muscle could assist in diagnosis, management and treatment plans for low back pain (LBP). MRE could provide stiffness of PM muscle, thereby unravelling the pathophysiology of LBP. The directional filter could separate wave field components and could remove wave interference so that accurate stiffness values could be obtained from wave inversion algorithms [37]. Hence, we aimed to evaluate the use of directional filter in wave field maps and in determining mean stiffness values of PM muscles by two-wave inversion algorithms, namely the local frequency estimate (LFE) and algebraic inversion of differential equation (AIDE). The uniqueness of the present study is, this is the first study that evaluates the stiffness of PM muscles. The directional filter, LFE and AIDE are explained in brief in the following section.
2. Materials and Methods
2.1. MREWave Phantom Data
MREWave phantom datasets were obtained for use as reference, from the freely available MRE/Wave software for LFE wave inversion. Figure 1 shows a sample wave image from an MREWave agar gel phantom with four cylindrical inclusions (ranging from 5 mm to 25 mm in diameter) perpendicular to the slice, acquired with 100 Hz vibrational frequency. The reference value for background gel is 2.9 kPa and the inclusions are 6.4 kPa in the phantom.
2.2. Experimental Setup
All PM-MRE experiments were performed on a clinical magnetic resonance scanner (Achieva 3.0 T; Philips Healthcare, Best, The Netherlands) by using a six-element torso phased-array coil with the volunteer in the supine position. A self-built waveform generation system (LabView, USB-6221; National Instruments, Austin, TX, United States) was used to generate a vibration waveform. This system can generate sinusoidal waveforms with arbitrary frequencies and phases. In order to synchronize the vibration with repetition time, a transistor-transistor logic signal from the MRI system (radiofrequency pulse power amplifier) was used as a trigger to start the vibration [22]. Four phase offset images at π/2, π, 3π/2 and 2π were acquired at 50 Hz frequency. Thus, for a single set of MRE images, the imaging sequence was played four times. In this MRE system, the vibration offset was controlled by the waveform generator system that provided continuous (steady) vibrations throughout the whole acquisition (each imaging cycle). A power amplifier (XTi 1000; Crown International, Los Angeles, CA, United States) and a pneumatic pressure generator (Subwoofer TIT320C-4 12”; Dayton Audio, Dayton, OH, United States) was used to supply vibrations through a hose to a vibration pad. The vibration pad was designed by using a three-dimensional printer (3DTouch; 3D System, Rock Hill, SC, United States) that conformed the lower back region of the volunteer.
2.3. Jelly Phantom
A fruit jelly phantom MRE was performed using SENSE Flex S coils at 50 Hz frequency. The two jelly packets were inserted in a cylindrical structure and were vibrated to acquire the phase images. The jelly phantom data were used to test
Figure 1. Wave Images of MREWave phantom data (one of 8 phase offsets). A: Original bandpass filtered wave image; B: Top-down DF wave image; C: Bottom-up DF wave image of MRE/Wave agar gel phantom. DF: Directional filtration.
the directional filter. The cylindrical structure and the vibration pad were self-designed and produced using the 3D printer.
2.4. Volunteer Studies
All volunteer studies were carried out after obtaining ethical approval from the local Institutional Review Board and obtainment of consent to participate from all volunteers. Seventeen healthy male volunteers (17 men; mean age: 26 ± 1.77 years, age range: 20 - 45 years) with no clinical history of skeletal muscle diseases, lumbar trauma, or recent/present LBP were included in this study. During the PM data acquisition time, the volunteers were instructed to take a deep inspiration and hold their breath. The breath hold was mandatory to avoid motion artefacts and obtain good image quality. The knees were flexed by placing a cushion beneath the legs; the hands were in normal anatomical position and a head cushion was placed under the head for comfort. We selected low vibration frequency (50 Hz) because the penetrating power of this vibration frequency is greater than that of a high vibration frequency (i.e. 100 Hz). The PM muscle lies deep inside, on either side of the lumbar vertebrae; it originates from the transverse processes of the twelfth thoracic vertebra (T12)-fifth lumbar vertebra (L5) with the lateral aspects of intervertebral discs between them and inserts in the lesser trochanter of the femur. The lumbar vertebrae are vibrated by efficiently transmitted vibrations to the PM muscle attached on either side. Adequate reach of vibration frequency is essential for better propagation of wave ripples from the lumbar vertebral body perpendicular to the direction of the muscle fibre of PM muscle.
The vibration pad was placed over the dorsal side at the intervertebral disc level of the third-fourth lumbar vertebra (L3-L4) of the PM muscle, and was fixed using Velcro straps and an auxiliary positioning foam. In addition, a silicon sheet was placed between the vibration pad and the skin surface on the PM muscle. The silicon sheet acted as a substitute vibration membrane and served to prevent air from leaking through the interstices between the vibration pad and the skin on the PM muscle.
2.5. MRE
A localizer scan was acquired in axial, sagittal and coronal planes. First, the MRE sequence was performed without vibration. The acquired data was used as a reference for determining the muscle shape and for positioning a slice. The MRE sequence with vibration was then performed four times with phase offset at π, π/2, 3π/2 and 2π at the level of L3-L4. The MRE acquisitions were conducted with a gradient-echo type multi-echo magnetic resonance sequence [35]. In this sequence, multiple symmetrical gradient-echoes can be acquired by a symmetrical bipolar read-out gradient lobe. When synchronized with vibration, those gradient lobes allow for achievement of a MEG-like effect. The synchronization can be achieved by adjusting the period of bipolar read-out gradient lobes and the gap between the first and next echo (δTE) of the sequence. When δTE is synchronized with the vibration frequency, a maximum MEG-like effect can be attained. Moreover, the later generated echoes have greater MEG-like effect (1st echo < 2nd echo < 3rd echo, etc.). The second echo data was selected because the MEG-like effect of the 1st echo data was not sufficient for adequate wave propagation and those later than 3rd echo data were affected by transverse relation and phase wrapping due to magnetic susceptibility effect.
In this study, a continuous (steady state) acoustic vibration at 50 Hz was supplied, synchronized with 10 ms δTE. The magnetic resonance phase images were conducted with the following parameters: 512 × 512 acquisition matrices; 4 number of averages; 2 sensitivity encoding; 20˚ flip angle; 300 - 330 mm field of view; 10 mm slice thickness; 2.3 ms 1st transient elastography; 12.3 2nd transient elastography; 40 ms repetition time; anteroposterior read-out direction (MEG-direction); and 4 phase offsets. Prior studies have shown that shear wave displacements are generated primarily along the direction of muscle fibres [38]. However, in this study, the direction of shear wave propagated orthogonal to the direction of the PM muscle fibre and the imaging plane was set perpendicular to the muscle fibre. Thus, the read-out gradient direction was set perpendicular to the long axis of the PM muscle.
2.6. Image Processing and Analysis
All phase images were processed by the Laplacian-based-estimate phase unwrapping method [39] [40]. The phase unwrapped images were then filtered by a directional filter self-customized in MATLAB 2018b (MathWorks, Inc., Natick, MA, United States). First, the directional filter was tested in MREWave phantom and jelly phantom images for benchmark evaluation. The wave images were Fourier-transformed in temporal direction, the mask of Figure 2 was applied in Fourier domain and again inverse Fourier transformed to obtain the directional-filtered wave images. The elastograms of the MREWave phantom data were calculated by both LFE and AIDE. Before stiffness calculation by AIDE, Savitzky-Golay filter, 3 × 3 gaussian filter and 3 × 3 median filter were used for
Figure 2. Directional Filter. A: Right Directional Filter (Left-to-Right) B: Left Directional Filter (Right-to-Left). DF: Directional filtration.
suppression of image noise. Subsequently, the directional filter was applied in PM images. The directional filter [41] [42] was designed in the spatiotemporal domain and was applied to the direction of propagation of the wave. One set of images was obtained by applying both directional filter (DF) and gaussian bandpass filter (GBF). Another set of images was obtained by applying GBF only, without applying DF. In right DF, the direction of wave filtration was from left-to-right side (L → R), whereas in left DF, the direction of wave filtration was from right-to-left side (R → L). Thus, two sets of images, viz. directional filtered (DF), and non-directional filtered (non-DF) were obtained and compared after conversion into wave displacement images (wave images) (MRE/Wave; Mayo Clinic, Rochester, MN, United States). The region of interest was manually drawn along the inner surface of the PM muscle to determine stiffness value. The stiffness calculation was performed in right PM from right DF wave images and in left PM from left DF wave images. In other words, the main wave direction in PM and the direction of DF must be same for accurate stiffness reconstruction.
2.7. Directional Filter (DF)
A directional filter can select the wave fields propagating in a specific direction. A unidirectional filter [37] is an image filter in the Fourier domain (frequency plane) by using cos2 dependence in the half-domain about a selected direction and zero in the opposite half-domain as shown in Figure 2. A directional filter is similar to the lognormal filters used in LFE [43]. The filtered data was extracted from the first positive temporal frequency plane (kt = +1). Another temporal frequency plane (kt = −1), which is conjugate symmetric to (kt = +1), was discarded.
2.8. Wave Inversion Algorithms
The elastogram (shear stiffness map) was calculated by two-wave inversion algorithms, namely local frequency estimate (LFE) and the algebraic inversion of differential equation (AIDE). Both DF and non-DF wave images were converted into elastogram using the LFE and AIDE respectively. The region of interest drawn manually along the inner surface of the PM muscle was used to measure the stiffness value.
2.8.1. LFE
The wave inversion algorithm was performed using MRE/Wave (developed by Richard L. Ehman at the Mayo Clinic). In this LFE, longitudinal waves are removed by curl processing and DF is applied to remove wave interferences. LFE combines local estimates of instantaneous frequency over several scales [13] [43], derived from filters that are considered to be oriented lognormal quadrature wavelets (a product of radial and directional components). The shear stiffness calculated by LFE is a result of solving the Helmholtz equation assuming no attenuation; that is,
, (where, µ is the shear modulus, fmech is the vibration frequency and, fsp is the spatial frequency), taking ρ = 1.0 for all soft tissues and under the assumption of local homogeneity, incompressibility. The LFE algorithm is insensitive to noise but its main drawback is limited resolution; a correct estimate cannot be reached at the boundaries and the correct estimate is reached at only half a wavelength into a given region.
2.8.2. AIDE
Another method of wave inversion algorithm used in this study is AIDE [44]. Assuming incompressible AIDE, the shear modulus from a single polarization of motion can be estimated as
, where,
is complex shear modulus, ω (2πf) is angular frequency, f is vibration frequency,
is displacement field and
is the spatial Laplacian of the displacement field. Assuming ρ = 1.0 for all soft tissues, AIDE can calculate complex shear modulus, where its real part gives storage shear modulus and imaginary part is the loss shear modulus. Unlike LFE, AIDE does not depend on planar shear wave propagation but simply on the presence of motion. Notably, the AIDE result is not affected by wave interference, and we expect low amplitude areas to have low signal-to-noise ratio. This is also valid for LFE [13].
2.9. Statistical Analysis
The mean stiffness values of PM muscle from LFE and AIDE with/without using directional filter were compared by using One-way Analysis of Variance (ANOVA) test. P-value at 5% level of significance was considered. All data were analysed using the SPSS software, version 24 (IBM Corp., Armonk, NY, United States).
3. Results
Figure 1 shows the MRE/Wave phantom wave images that are bandpass filtered, top-down directional filtered, and bottom-up directional filtered. The shear waves in the bandpass filtered images showed wave interference from inclusions and the bottom wall, whereas the wave interference was removed in top-down DF images. The bottom-up DF image showed waves traveling in the direction opposite to top-down DF images.
Figure 3 shows the MRE/Wave phantom elastograms by both LFE and AIDE. There was clear improvement in the quality of elastogram in DF elastograms. The stiffness values of the phantom data sets are shown in Table 1 (the biggest inclusion was only calculated for ease). There was clear improvement in the standard deviation of inclusion by LFE and AIDE; however, there was little change in the standard deviation of background agar gel by LFE and no change by AIDE.
Figure 4 shows the wave images for jelly phantom without DF, right DF and left DF. The wave propagation was right-to-left direction for right jelly and left-to-right direction for left jelly. However, the directional filter applied to the right direction changed the direction of wave propagation to left-right direction and vice versa.
Table 1. Stiffness values of MREWave phantom with/without DF.
DF: Directional filtration; LFE: Local frequency estimate.
Figure 3. Elastograms of MREWave phantom data. A, B: Elastogram by LFE with and without DF; C, D: Elastogram by AIDE with and without DF. AIDE: Algebraic inversion of differential equation; DF: Directional filtration; LFE: Local frequency estimate.
Figure 4. T2-Weighted Image and Wave images of fruit jelly phantom (one of four phase offsets). Wave image of psoas major with and without DF. DF: Directional filtration; Lt: Left; Rt: Right.
Figure 5 shows the non-DF filtered wave images of PM, where the wave direction was radiating outwards in bilateral sides from the lumbar spine, perpendicular to the direction of muscle fibre. The right DF wave images showed clear wave propagation in right PM from left-to-right in the direction of the filter applied. The left DF wave images showed clear wave propagation in left PM from right-to-left in the direction of the filter applied.
The range of stiffness values calculated by LFE of non-DF images was 1.39 ± 0.25 kPa and 1.33 ± 0.22 kPa for right and left PM respectively, whereas for right DF images, it was 1.26 ± 0.20 kPa for right PM and for left DF images, it was 1.18 ± 0.28 kPa for left PM.
The range of stiffness values calculated by AIDE of non-DF images was 0.78 ± 0.10 kPa and 0.78 ± 0.13 kPa for right and left PM respectively, whereas for right DF images, it was 0.73 ± 0.12 kPa for right PM and for left DF images, it was 0.74 ± 0.11 kPa for left PM.
There was no statistically significant difference in mean values of stiffness with/without applying directional filter as shown in Table 2, whereas there was a statistically significant difference in mean values of stiffness between LFE and AIDE for both with/without DF (p < 0.05). The stiffness fusion maps are shown in Figure 6.
Table 2. Stiffness values by LFE and AIDE with/without using DF.
AIDE: Algebraic inversion of differential equation; DF: Directional filtration; LFE: Local frequency estimate; Lt: Left; Rt: Right.
Figure 5. Wave image of PM muscle (one of four phase offsets). Wave image of PM with and without DF. DF: Directional filtration; PM: Psoas major.
Figure 6. Magnitude-Elastogram fusion images for LFE and AIDE with DF and without DF. AIDE: Algebraic inversion of differential equation; DF: Directional filtration; LFE: Local frequency estimate, R: Right, L: Left.
4. Discussion
The purpose of this study was to evaluate the use of directional filter in wave image processing and to analyse the mean stiffness values of PM muscle by using LFE and AIDE in MRE. The wave propagation in PM muscle was bi-directional. However, uni-directional filtering could only visualize wave-propagation in single direction, i.e. either right or left PM. Directional filter successfully ensured clear wave propagation in PM muscles. In DF images, wave interferences were removed, thus we assume that the region of low motion amplitude (low signal areas) from the vertebral body and spine and the cancellation between the muscle interfaces were eliminated.
In addition, we presume that the quality of stiffness map was also improved due to application of directional filter, because the wave inversion algorithms (LFE and AIDE) are sensitive to low displacement amplitude and the directional filter removed low displacement amplitude (low signal-to-noise ratio areas). However, the stiffness values with/without using directional filter were not statistically significant. Thus, the directional filter did not affect the stiffness output. However, the stiffness output from LFE and AIDE were statistically significant in both cases of with/without applying directional filter, which implies that the stiffness value is dependent upon the type of usage of wave inversion method, and the directional filter has no significant effect on the elasticity values. However, directional filter ensured the accurate reconstruction of elastogram. Both LFE and AIDE are promising and widely adapted wave inversion methods in MR elastography, researchers could decide on their own which one to use. Since, both are inverting Helmholtz equation, we recommend using both methods according to own requirements and preferences.
The stiffness values could also be used as a bio-imaging marker for low back pain (LBP). Since LBP is almost non-specific, MR elastography technique of psoas major using directional filter and wave inversion algorithms (LFE and AIDE) could provide radiate important insights for the clinical implications for LBP. Numano et al. [45] have described the clinical application of MR elastography in psoas major in LBP in more details. We presume the usage of directional filter, LFE and AIDE would make PM-MRE technique more comprehensive.
There are a few limitations of the present study that must come under consideration. First, attenuation and loss modulus were not considered since loss modulus cannot be measured from LFE and it assumes zero attenuation [13]. Though AIDE could calculate loss modulus and an attenuation map [37], it was impossible to compare LFE and AIDE based on these parameters. That being said, Manduca et al. [37] have shown that the visual clarity of an attenuation map is enhanced by using DF. Second, we performed MRE experiments at 50 Hz frequency only with read-out gradient AP direction. Numano et al. [35] [45] had stated that the 50 Hz frequency is an optimal frequency of vibration for the PM muscle and that the anteroposterior read-out gradient direction is best for wave detection in the PM muscle. Again, we used 2D directional filter since our image data is 2D. However, 3D and 4D directional filter could remove wave reflection artifacts [46]. Furthermore, our sample size was relatively small and only one operator performed MRE experiments. Hence, repeatability and reproducibility of the experiment were not computed using multiple experimenters and evaluators. Despite these constraints, we evaluated the use of DF in PM-MRE. Future work could be performed to evaluate the repeatability and reproducibility of the data obtained using the present technique. Again, the latest wave inversion methods like multifrequency algorithms and convolutional neural networks could be used to determine the shear modulus (or stiffness).
5. Conclusion
We concluded that using directional filter, wave propagation is better visualized in the PM muscle that yielded smooth wave fields. Directional filter is useful for wave image processing. Both LFE and AIDE wave inversion could be used in psoas major MR elastography.
Acknowledgements
The authors acknowledge the Tokyo Metropolitan Government for providing grant assistance. The authors also thank all the volunteers who took part in the research, for their assistance in facilitating this study.
Supplementary Material
The movie files of the wave images of the fruit jelly phantom and psoas major muscle can be accessed at: https://github.com/surendra116083/suren.