rus
DOI: 10.22184/1993-8578.2021.14.3-4.188.195
Diffusion MRI is a magnetic resonance imaging method that estimates random molecular motion of molecules in biological tissue. The most popular its application consists in disentangling of tissue anisotropy. Nowadays high resolution imaging techniques allow to make conclusions not only about brain white matter structure but also a gray matter attracts a great attention. However, such potentially valuable MRI methods are not used in clinical environment routinely because they require a long acquisition time. In this work we process low-resolution diffusion MRI and calculate high-resolution images by means of super-resolution reconstruction. We calculate diffusion-tensor and estimate quality of its fitting comparing with alternative interpolation methods. We demonstrate benefit of super-resolution technique in terms of root-mean-square error and by analyzing fractional anisotropy and other tensor invariants.
Diffusion MRI is a magnetic resonance imaging method that estimates random molecular motion of molecules in biological tissue. The most popular its application consists in disentangling of tissue anisotropy. Nowadays high resolution imaging techniques allow to make conclusions not only about brain white matter structure but also a gray matter attracts a great attention. However, such potentially valuable MRI methods are not used in clinical environment routinely because they require a long acquisition time. In this work we process low-resolution diffusion MRI and calculate high-resolution images by means of super-resolution reconstruction. We calculate diffusion-tensor and estimate quality of its fitting comparing with alternative interpolation methods. We demonstrate benefit of super-resolution technique in terms of root-mean-square error and by analyzing fractional anisotropy and other tensor invariants.
Теги: diffusion mri diffusion-tensor fractional anisotropy magnetic resonance iomography root-mean-square error tensor invariants volt-ampere characteristic диффузионная мрт дробная анизотропия магнитно-резонансная томография среднеквадратичная ошибка тензор диффузии тензорный инвариант
Application of super-resolution method in diffusion MRI
INTRODUCTION
Investigation of diffusion properties of brain tissue are often limited by the magnetic resonance image (MRI) resolution. For example, resolution plays a key role in brain gray matter segmentation and anisotropy analysis. Such major limiting factor as a partial volume effect hinders a fine image consideration. Another important challenge in image processing is a signal-to-noise (SNR) ratio which reduces with increase of resolution and downsizing of voxel dimension. One of the possible solution of the problem is to acquire images with thick slices and high resolution in plane and perform retrospective post-processing. However, initial 3D MR data are not isotropic which might bias diffusion tensor parameters [1].
On the first view, applying interpolation technique [2] to diffusion tensor may be necessary to improve tissue parameters estimation. In general image interpolation is a very common procedure in many processing pipelines and it may have a significant influence on other processing procedures [3] such as registration, segmentation or diffusion tensor field calculation. Applying interpolation improves visual inspection of the image but, on the other hand, blur and contrast loss artifacts become pronounced.
Another method which improves image quality is a super-resolution (SR) [4, 5]. Opposite to the interpolation it concatenates a new information from several voxels providing a higher reliability of MRI which, in general, adds to the images neuroscientific or clinically useful diagnostic value.
Modern MRI scanners are equipped with sequences allowing sub-millimeter resolution. However, high-resolution diffusion weighted imaging (DWI) is still challenging because of high level of physiological and systematic noise, involuntary motion of the subject during scanning, and long acquisition time. Small lesions in gray matter is difficult or even not possible to visualize and then perform a fair analysis. Subtle changes in biological tissue during maturation or aging, Alzheimer decease, and plaques development, for example, are very tiny for MRI. This hinders to improve specificity and sensitivity adjustment of diffusion sequences’ protocol parameters.
In this circumstances SR possess a number of advantages, such as processing high SNR low-resolution images acquired in a reasonable acquisition time, independence from main magnetic field strength and a number of receivers.
In this short communication we describe how process low-resolution through plane DWIs and reconstruct a high-resolution image by means of SR method. To achieve isotropic resolution several images with axial sub-voxel shifts were created. We calculate diffusion-tensor and estimate quality of its fitting comparing with alternative interpolation (linear, cubic, spline) methods. We demonstrate benefit of SR technique in terms of root-mean-square error (sse) and by analyzing fractional anisotropy and diffusion tensor invariants in the human brain tissue in vivo.
METHODS
Data generation
The main approach to implement SR consists in combining several low resolution images in high-resolution one if it is known the sub-voxel transform of the low-resolution images (Fig.1a). In particularly, in our work we used through plane (slice encoding direction) translation of low resolution images. Modern MR scanner consoles are developed in a way to be flexible and meet all demands of the MR researcher. The position and the slice orientation of the images can be easily adjusted to achieve consecutive number of shifts in desired direction.
In vivo human MRI data were downloaded from open database db.humanconnectome.org. Images were acquired on a 3T Siemens scanner (Siemens, Erlangen, Germany) with a 32-channel RX head coil. Ninety volumes of whole-brain DWIs were acquired in 90 noncollinear and noncoplanar directions ([gx, gy, gz]: b = 1000 s/mm2) and interleaved three non-diffusion weighted images (b = 0 s/mm2) using the single-refocused Stejskal-Tanner spin-echo sequence (FoV = 210 × 180 × 139 mm3, voxel size = 1,25 × 1,25 × 5 mm3 after coarse graining, in plane resolution 168 × 144, GRAPPA2, TE/TR = 89,5/5520 ms, partial Fourier 6/8, bandwidth 1488 Hz/px) [6]. For excitation of 3 slices simultaneously multi-band radio frequency (mb-rf) pulse was used. Additionally, a set of reversed phase-encoded gradient b = 0 s/mm2 images were acquired for correction of geometrical distortions. Diffusion gradient was characterized with diffusion time Δ = 22 ms, diffusion gradient duration δ = 6 ms. In order to perform a super-resolution reconstruction images were shifted along the slice encoding direction 4 times with a step 1.25 mm (Fig.1a).
Data processing
Before applying high resolution reconstruction DWIs were processed using mrtrix [7] and FSL [8] (Fig.1b). The steps of correction included (1) noise and Gibbs ringing artifact reduction; (2) susceptibility and eddy current distortion correction was performed using topup and eddy algorithms which processed a set of reversed phase-encoded gradient b = 0 s/mm2 images; (3) for distortion free images radio-frequency bias field correction was applied.
We performed calculation of DTI of the low resolution images. After reconstruction of diffusion tensor, it was interpolated in 1D slice encoding direction with subsequent tensor estimation. We also obtained tensor invariants, i.e. mean diffusivity (MD), fractional anisotropy (FA), radial and axial diffusivities (RD, AD). Interpolation was done with a log-Euclidean transformation to satisfy the relation between the DWI signal and the tensor elements [9]. Interpolation was performed with three different methods, i.e. linear, cubic, and spline. Procedure of interpolation was implemented with matlab [10] home developed scripts.
SR method was implemented by the following procedure. Every voxel of images was axially divided on the number of sub-voxels corresponding the number of axial shifts. Then shifted images were translated to the same coordinate system (Fig.1a). Afterwards overdetermined liner system was resolved by three methods. In the method I we calculated least mean square value for every voxel. For the method II we applied Gaussian convolution kernel with maximal weights corresponding maximal signal. The method III consisted in multiplicative averaging.
Mathematically such procedure can be represented by the following system of linear equations:
x = ABy, (1)
where B is an operator of dividing of the voxel, and A is an operator of shift in order to translate images in the uniform space. We found an optimal value c by averaging of vector x:
c = T+TKx, (2)
where T+T is a pseudo-inverse matrix, K is a kernel operator. In case of method I of super resolution K = L, where L is a box-car matrix. For method II K = G, where G is a set of coefficients of a Gaussian kernel. Method III assumed multiplicative averaging:
с = (Рх)а, (3)
where P is a product operator and a is equal to the inverse number of shifts. Algorithms of SR were implemented in matlab [10] home developed scripts.
After reconstruction of SR images, diffusion tensor was calculated with subsequent eigenvalues and eigenvector estimation. We also estimated tensor invariants, i.e. MD, FA, RD, and AD.
Additionally, sse of diffusion tensor fitting was calculated for all cases. We used histograms of estimated values (sse, MD, FA, RD, AD) to compare different methods.
RESULTS
In Fig.2 we demonstrate histograms of different diffusion tensor metrics and sse of tensor fitting:
sse = log(S)T[I-X(XTX)-1XT)]log(S),
X = gxgT. (4)
In Eq.(4) g = [gx, gy, gz] is a gradient direction, S is a diffusion weighted signal, I identity matrix, and T is a transpose. Histograms were calculated for different types of tensor interpolation and SR methods. Sum of all frequencies there was normalized to 1. Histograms of sse (Fig.2a) demonstrate a strong asymmetry with maximum located in vicinity of zero. Major statistical parameters, i.e. mean, median, standard deviation, are presented in figure legend and provide an evidence of highest quality of tensor fitting in case of SR I method. Superiority of SR can be explained by the procedure of reconstruction which assumes combining of several voxels in one. Unexpectedly broader distribution of sse for cubic and spline tensor interpolation can be explained polynomial lobes near nodes’ points keeping the property of continuity and smoothness of curves. It is worth mentioning that histogram of sse for tensor fitting for original low resolution dataset is close to linear interpolation method.
In Figs.2b-e we show histograms of different tensor metrics calculated with different methods. Metrics for original dataset is also included. Obviously that they demonstrate very close distributions in general. As a consequence, FA maps of metrics (Fig.2f) are very similar with variation in small areas.
DISCUSSION AND CONCLUSIONS
Super-resolution method is efficient, when acquired voxel size is larger then small areas, and partial volume effect hinders correct characteristic of clinically significant structures, for example in gray matter. Our results demonstrate that the application of well-established super-resolution methods to diffusion weighted imaging datasets, which might do not have high spatial resolution, can potentially reveal anatomical details only seen for the images with high resolutions. This is not possible to achieve with clinical routine diffusion imaging protocol, which produces images coarser compared to the size of the gray and white matter anatomical structures. In the processing steps registration of images, acquired from different scans, plays important role if motion, physiological and systematic noise are present. Also accuracy of distortion correction may influence on the final reconstruction. ■
Acknowledgement. This work is supported by the ERC funding and BIF-IZKF UKA RWTH.
Declaration of Competing Interest. The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
INTRODUCTION
Investigation of diffusion properties of brain tissue are often limited by the magnetic resonance image (MRI) resolution. For example, resolution plays a key role in brain gray matter segmentation and anisotropy analysis. Such major limiting factor as a partial volume effect hinders a fine image consideration. Another important challenge in image processing is a signal-to-noise (SNR) ratio which reduces with increase of resolution and downsizing of voxel dimension. One of the possible solution of the problem is to acquire images with thick slices and high resolution in plane and perform retrospective post-processing. However, initial 3D MR data are not isotropic which might bias diffusion tensor parameters [1].
On the first view, applying interpolation technique [2] to diffusion tensor may be necessary to improve tissue parameters estimation. In general image interpolation is a very common procedure in many processing pipelines and it may have a significant influence on other processing procedures [3] such as registration, segmentation or diffusion tensor field calculation. Applying interpolation improves visual inspection of the image but, on the other hand, blur and contrast loss artifacts become pronounced.
Another method which improves image quality is a super-resolution (SR) [4, 5]. Opposite to the interpolation it concatenates a new information from several voxels providing a higher reliability of MRI which, in general, adds to the images neuroscientific or clinically useful diagnostic value.
Modern MRI scanners are equipped with sequences allowing sub-millimeter resolution. However, high-resolution diffusion weighted imaging (DWI) is still challenging because of high level of physiological and systematic noise, involuntary motion of the subject during scanning, and long acquisition time. Small lesions in gray matter is difficult or even not possible to visualize and then perform a fair analysis. Subtle changes in biological tissue during maturation or aging, Alzheimer decease, and plaques development, for example, are very tiny for MRI. This hinders to improve specificity and sensitivity adjustment of diffusion sequences’ protocol parameters.
In this circumstances SR possess a number of advantages, such as processing high SNR low-resolution images acquired in a reasonable acquisition time, independence from main magnetic field strength and a number of receivers.
In this short communication we describe how process low-resolution through plane DWIs and reconstruct a high-resolution image by means of SR method. To achieve isotropic resolution several images with axial sub-voxel shifts were created. We calculate diffusion-tensor and estimate quality of its fitting comparing with alternative interpolation (linear, cubic, spline) methods. We demonstrate benefit of SR technique in terms of root-mean-square error (sse) and by analyzing fractional anisotropy and diffusion tensor invariants in the human brain tissue in vivo.
METHODS
Data generation
The main approach to implement SR consists in combining several low resolution images in high-resolution one if it is known the sub-voxel transform of the low-resolution images (Fig.1a). In particularly, in our work we used through plane (slice encoding direction) translation of low resolution images. Modern MR scanner consoles are developed in a way to be flexible and meet all demands of the MR researcher. The position and the slice orientation of the images can be easily adjusted to achieve consecutive number of shifts in desired direction.
In vivo human MRI data were downloaded from open database db.humanconnectome.org. Images were acquired on a 3T Siemens scanner (Siemens, Erlangen, Germany) with a 32-channel RX head coil. Ninety volumes of whole-brain DWIs were acquired in 90 noncollinear and noncoplanar directions ([gx, gy, gz]: b = 1000 s/mm2) and interleaved three non-diffusion weighted images (b = 0 s/mm2) using the single-refocused Stejskal-Tanner spin-echo sequence (FoV = 210 × 180 × 139 mm3, voxel size = 1,25 × 1,25 × 5 mm3 after coarse graining, in plane resolution 168 × 144, GRAPPA2, TE/TR = 89,5/5520 ms, partial Fourier 6/8, bandwidth 1488 Hz/px) [6]. For excitation of 3 slices simultaneously multi-band radio frequency (mb-rf) pulse was used. Additionally, a set of reversed phase-encoded gradient b = 0 s/mm2 images were acquired for correction of geometrical distortions. Diffusion gradient was characterized with diffusion time Δ = 22 ms, diffusion gradient duration δ = 6 ms. In order to perform a super-resolution reconstruction images were shifted along the slice encoding direction 4 times with a step 1.25 mm (Fig.1a).
Data processing
Before applying high resolution reconstruction DWIs were processed using mrtrix [7] and FSL [8] (Fig.1b). The steps of correction included (1) noise and Gibbs ringing artifact reduction; (2) susceptibility and eddy current distortion correction was performed using topup and eddy algorithms which processed a set of reversed phase-encoded gradient b = 0 s/mm2 images; (3) for distortion free images radio-frequency bias field correction was applied.
We performed calculation of DTI of the low resolution images. After reconstruction of diffusion tensor, it was interpolated in 1D slice encoding direction with subsequent tensor estimation. We also obtained tensor invariants, i.e. mean diffusivity (MD), fractional anisotropy (FA), radial and axial diffusivities (RD, AD). Interpolation was done with a log-Euclidean transformation to satisfy the relation between the DWI signal and the tensor elements [9]. Interpolation was performed with three different methods, i.e. linear, cubic, and spline. Procedure of interpolation was implemented with matlab [10] home developed scripts.
SR method was implemented by the following procedure. Every voxel of images was axially divided on the number of sub-voxels corresponding the number of axial shifts. Then shifted images were translated to the same coordinate system (Fig.1a). Afterwards overdetermined liner system was resolved by three methods. In the method I we calculated least mean square value for every voxel. For the method II we applied Gaussian convolution kernel with maximal weights corresponding maximal signal. The method III consisted in multiplicative averaging.
Mathematically such procedure can be represented by the following system of linear equations:
x = ABy, (1)
where B is an operator of dividing of the voxel, and A is an operator of shift in order to translate images in the uniform space. We found an optimal value c by averaging of vector x:
c = T+TKx, (2)
where T+T is a pseudo-inverse matrix, K is a kernel operator. In case of method I of super resolution K = L, where L is a box-car matrix. For method II K = G, where G is a set of coefficients of a Gaussian kernel. Method III assumed multiplicative averaging:
с = (Рх)а, (3)
where P is a product operator and a is equal to the inverse number of shifts. Algorithms of SR were implemented in matlab [10] home developed scripts.
After reconstruction of SR images, diffusion tensor was calculated with subsequent eigenvalues and eigenvector estimation. We also estimated tensor invariants, i.e. MD, FA, RD, and AD.
Additionally, sse of diffusion tensor fitting was calculated for all cases. We used histograms of estimated values (sse, MD, FA, RD, AD) to compare different methods.
RESULTS
In Fig.2 we demonstrate histograms of different diffusion tensor metrics and sse of tensor fitting:
sse = log(S)T[I-X(XTX)-1XT)]log(S),
X = gxgT. (4)
In Eq.(4) g = [gx, gy, gz] is a gradient direction, S is a diffusion weighted signal, I identity matrix, and T is a transpose. Histograms were calculated for different types of tensor interpolation and SR methods. Sum of all frequencies there was normalized to 1. Histograms of sse (Fig.2a) demonstrate a strong asymmetry with maximum located in vicinity of zero. Major statistical parameters, i.e. mean, median, standard deviation, are presented in figure legend and provide an evidence of highest quality of tensor fitting in case of SR I method. Superiority of SR can be explained by the procedure of reconstruction which assumes combining of several voxels in one. Unexpectedly broader distribution of sse for cubic and spline tensor interpolation can be explained polynomial lobes near nodes’ points keeping the property of continuity and smoothness of curves. It is worth mentioning that histogram of sse for tensor fitting for original low resolution dataset is close to linear interpolation method.
In Figs.2b-e we show histograms of different tensor metrics calculated with different methods. Metrics for original dataset is also included. Obviously that they demonstrate very close distributions in general. As a consequence, FA maps of metrics (Fig.2f) are very similar with variation in small areas.
DISCUSSION AND CONCLUSIONS
Super-resolution method is efficient, when acquired voxel size is larger then small areas, and partial volume effect hinders correct characteristic of clinically significant structures, for example in gray matter. Our results demonstrate that the application of well-established super-resolution methods to diffusion weighted imaging datasets, which might do not have high spatial resolution, can potentially reveal anatomical details only seen for the images with high resolutions. This is not possible to achieve with clinical routine diffusion imaging protocol, which produces images coarser compared to the size of the gray and white matter anatomical structures. In the processing steps registration of images, acquired from different scans, plays important role if motion, physiological and systematic noise are present. Also accuracy of distortion correction may influence on the final reconstruction. ■
Acknowledgement. This work is supported by the ERC funding and BIF-IZKF UKA RWTH.
Declaration of Competing Interest. The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Readers feedback
