Image Processing in Magnetic Resonance Imaging
Image Processing in Magnetic Resonance Imaging |
|
Todd B. Parrish, PhD has reported no financial interest, arrangement or affiliation with a commercial organization that may have a direct or indirect influence in the subject matter of this presentation.
Introduction
In the field of MRI there is a great deal of image processing that goes into the formation of an image. For example, the demodulation of the raw signal that is detected by the coil, the filtering that takes places to eliminate aliasing, and finally the Fourier transform that creates the image. These issues are complex and are utilized every time the MR system acquires data. The image and signal processing theories behind the formation of the MR image will not be discussed in this look at Image Processing of MR images. Instead, the goal is to develop an understanding of how basic MR data is manipulated to form the complex images that represent the physiologic properties being detected.
Image processing basics:
An image is comprised of 3 dimensional (3D) objects called voxels, or volume elements. (A voxel is a three-dimensional object, whereas a pixel is the two-dimensional representation on a screen.) In MRI, the field of view and number of readout and phase encode samples determines the size in the in-plane directions. The slice thickness governs the final dimension. Typical clinical MR image dimensions are 256x512 voxels and a field of view of 22cm that results in voxel dimensions of 0.86mm x 0.43mm x 4mm of tissue. This anisotropic voxel is usually acquired in 2D MR sequences such as FLAIR or T2-weighted sequences. When viewed "head-on" the image has high quality resolution (<1mm, Figure1).
 |
 |
| Figure 1 | Figure 2 |
However, if the stack of images is viewed from another direction, the resolution of the image is compromised by the large slice thickness because the resolution viewed is now 1mmx4mm (Figure 2). With the advent of fast imaging hardware and sequences, it has become feasible to collect 3D data in a reasonable amount of time. In this case the typical voxel dimensions are 1mm in all directions. Therefore, the data can be reformatted without degradation to the resolution of the image (Figure 3).
Figure
3
Since the data can be reformatted there is no need to acquire separate scans for different anatomic orientations such as axial and sagittal. The appropriate views could be created from the 3D dataset. Having the 3D data allows one to better localize lesions since the boundaries and size can be viewed in three different planes. Furthermore, the higher resolution in the slice direction reduces partial volumes effects.
The exploitation of the 3D nature of the data allows for localization of brain regions and surgical devices. Combining the MR data with a fiducial marker such as a frame, dots, or anatomic landmarks makes it possible to coregister the MR data with the skull in the operating room or gamma knife system. The coregistration process transforms the markers in the 3D MR data to the actual marker in the operating room. The surgeon identifies the markers and records their location. A spatial transformation warps (associates) the 3D MR data to the surgical coordinates. The registration works well for the static structures such as the skull, which allows for accurate placement of the craniotomy. There is a shift that occurs once the skull is opened because of pressure changes. There are methods available that allow for registration of the brain surface with the MR data but this requires some sort of imaging system within the operating room. A similar coregistration process would be used to align the MR data with the real-time image obtained in the operating room.
Another area of active image processing research is the field of image segmentation. The goal of segmentation is to identify different tissue components within the image data. The user can use a single image (T1-weighted 3D volume), or a multi-spectral analysis, which combines multiple types of images (T1-weighted, T2-weighted and proton density) to improve the classification of the tissue. The end result is a set of images of the different tissue types: gray matter, white matter, and CSF (Figure 4).
Figure 4
Using visualization software, it is possible to highlight each tissue with a different color and attributes to make it easier to visualize. The segmented data could also be used to monitor volume changes of specific tissues simply by converting the number of voxels of a particular tissue type to a volume based on the voxel dimensions. In MR data it is difficult to create the segmented image because the contrast within the image is dependent on the sequence used (i.e. different tissue types can have the same intensity). In contrast, CT has specific Houndsfield units, which correspond to different tissue types. Therefore a simple threshold segmentation algorithm, tissue type A = image intensities between X and Y, can be used. In MR data, this is difficult to do because image contrast changes with the sequence used, RF characteristics of the coil, and different tissues may have the same signal intensity. Therefore, more sophisticated algorithms or more information such as the multi-spectral approach need to be employed. Preprocessing of the data such as scalp and skull removal (in the image, not the person) simplifies the problem as well. Several different algorithms exist that remove the scalp, muscle, bone and eyes so only brain tissue, tumor if present, and CSF remain. After completing the scalping pre-processing, the segmentation problem becomes easier. Another preprocessing step includes RF inhomogeneity correction. The issue is that the same tissue type may vary in intensity from slice to slice due to flip angle variation over the volume of the coil. The data can be used to create a bias field that estimates how to correct the 3D volume.
There are a number of different types of segmentation algorithms that have been developed for use with MR brain data. A subset of the methods used for segmentation includes thresholding, region growing, multi-spectral analysis, and template matching. Region growing requires the introduction of "seeds" to start the analysis. This can be provided by the user or done randomly by the computer. These methods require good continuity (high resolution) and good image intensity stability in the images to allow the regions to merge into a single map of a tissue type. Multi-spectral analysis takes advantage of the different contrasts provided by MRI. By creating two- or three-dimensional plots of intensity from each scan, it is possible to isolate different tissues (Figure 5). Each tissue type is identified as a cluster and labeled. More advanced algorithms can assign outliers to different classes based on criteria set by the user. Template matching algorithms warp (transform) the 3D data to a known reference template, which has already been segmented. Each voxel is then assigned to a tissue type and then warped back to the original brain space. This method works well in the normal brain or whatever population was used to generate the template. However, errors can occur in a brain with lesions, atrophy, or when abnormalities exist. These variations from the template might cause a distortion in the transformation to the template, which in turn may lead to incorrectly identified tissue, or volume errors.
Figure
5
3D data manipulation
As mentioned previously, it is possible to reslice the 3D data to any orientation desired. The image processing involved in this method is not complex but it allows one to view anatomic data in a unique way. The algorithm applies a user-supplied plane to reslice the 3D dataset (Figure 6). An advanced feature of these algorithms is that the user can specify a line to reconstruct along. This type of cut is often disorientating because it takes a curved structure and straightens it out (Figure 7). A typical application of the reslicing technique is to resample the 3D data in a different orientation to provide two or more different views of the anatomy. Reslicing allows one to track 3D structures and understand their relationship to the rest of the brain. This is further enhanced because the operator can "surf" through the brain in real-time.
 |
 |
| Figure 6 | Figure 7 |
Another method of viewing 3D data is to create surface reconstructions. The surface is created by projecting rays through the data at the orientation specified by the user (Figure 8). These rays create a projection image of the 3D object much like x-rays make on film.
Figure
8
However, the information encoded in the ray can be the distance traveled after hitting the surface, the surface voxel intensity, the sum along the ray (simulates x-ray), or the intensity after it has penetrated the surface a certain distance. The result of these different ray definitions can create completely different images from the same dataset (Figure 9).
Figure 9
Surface reconstructions are useful once the non-brain tissues have been removed and do not obscure the brain. These views are often utilized in functional brain mapping to show the relationship between the activations and the anatomy (Figure 10).
Figure 10
A post-processing ray projection method that is used every day in clinical practice is called the maximum intensity projection, MIP. Similar to surface reconstructions, MIPs project the brightest signal along the ray onto an image (see figure 8). These projections are rotated slightly and recalculated. The final result is a movie of demonstrates the 3D vasculature rotating in space . When applied to the appropriate dataset, an MR angiogram (MRA) is created (Figure 11). In the case of MRA, any static tissue that is bright relative to the vessel would also project onto the image and cause an artifact. The tissues that usually cause this are fat, the orbits, or an enhancing lesion.
Figure
11
This can be addressed in the acquisition by using a magnetization transfer pulse or fat saturation pulses or it can be addressed in the post-processing phase by using a targeted MIP. The user would specify a region of interest to create the MIP images. Therefore, bright static tissues such as the scalp or fast flowing veins could be avoided to provide a better angiogram (Figure 12).
Figure
12
A close cousin to the MIP is the minimum intensity projection, mIP. As the name implies, the mIP would project the smallest signal found along the projection ray. In this application, the background is set to a large value so the projection image is not contaminated. The mIP application has use in black blood angiographic techniques. One method that has great potential is susceptibility weighted imaging developed by Dr. Haacke to visualize the venous system with incredible detail (Figure 13).
Figure 13
Functional MRI (fMRI)
Functional MRI is a field of Neuroimaging that detects small changes in blood oxygenation in cortical regions that are synchronized with the performance of a specific mental task (paradigm). During the paradigm presentation, the scanner collects many slices (~24-36 slices) using echo planar imaging in a short amount of time, ~2 seconds. The result is a time series of volume data. This creates a 4 dimensional (4D) data set. Since the magnitude of the activation induced blood oxygenation change is quite small, ~1%, several pre-processing steps improve the detection. After the pre-processing has been applied, a statistical test or model can be used to determine the significance of the changes that occur. Finally, the statistically significant regions, "blobs", are coregistered to a high-resolution dataset for display and interpretation. These steps will be addressed individually below.
Pre-processing
In all fMRI experiments there is a dummy period where the MR signal is allowed to reach a steady-state. Typically the first 10-12 seconds (5-6 dummy scans if TR=2s) of scanning is discarded. It is important to make sure that the appropriate amount of time is allotted in the paradigm design for the signal to reach steady state. Otherwise, there would be a timing error and the resultant activations would not reflect activation due to the particular portion of the paradigm.
Once the appropriate data has been retained, one must perform motion correction of the data. The goal of motion correction is to make sure that the volumes of data acquired are spatially aligned so that when a voxel is projected over time it represents the same brain tissue over time (Figure 14).
Figure
14
Typically, the user identifies a reference volume at a specific time point (default is the first image in the series). This volume of data is compared to the next in the time series. An iterative process of shifting and rotating the new volume to match the reference volume is performed. The goal is to minimize the difference between the two datasets. Once this has been accomplished the resultant shift and rotation is applied to that volume and is used as the starting point for the next volume. This continues until all the time points have been corrected. A plot of the shifts and rotations is informative about how well the subject behaved and to determine if there is any motion correlated with the stimulus presentation (Figure 15). Typical motion correction algorithms can correct for sub-voxel movements (small magnitude) but has difficulty in sharp, large shifts. This is due to the MR signal dependence on the RF pulsing history. Therefore, slow drifts in the data are acceptable but large (0.5mm or larger) volume-to-volume jumps are not. In the event that these jumps exist, one has to throw out the data or try to salvage portions that do not exhibit motion.
Figure 15
A controversial area of pre-processing for fMRI has been the issue of spatial and temporal smoothing. The benefit of smoothing is an increase in the apparent signal to noise. There is a dramatic effect with a small amount of smoothing, however, the downside of smoothing is that the localization of the activation is compromised. The range of spatial smoothing is anywhere from 3mm up to 12mm depending on the goal of the experiment. Obviously, with the larger smoothing kernel the goal of the research is to investigate larger network issues within the brain across a population. Larger smoothing kernels are also used when data are normalized to a standard template to improve the transformation process by removing individual differences. The issue of temporal smoothing is somewhat less of an issue as the hemodynamic response is sluggish compared to the TR used in the fMRI experiment. The effect of smoothing in the time dimension is to reduce the physiologic fluctuations that occur and help reduced volume-to-volume signal changes that are a result of the imaging hardware. Temporal smoothing is also used to dominate the autocorrelations that already exist in the MR data.
Analysis
The statistical analysis of functional data has always been a hot-bed of controversy. It starts with the fact that the data are not independent samples as many statistics assume. There is an auto-correlation issue with the formation of MR data especially if one smoothes the data (i.e. mix signals in time and space from different voxels). Furthermore, there is the issue of correcting for multiple comparisons. In a typical dataset, there are 131072 (64x64x32) voxels, which can be reduced based on some intensity threshold but the number is still quite large. What does a confidence level mean with such a large number of comparisons? For example a confidence level of 95% (a=0.05) means that you are willing to accept a 5% error. In the above example that would mean 6554 voxels would be an acceptable number of false positives. How do you adjust the statistic to account for all the comparisons being made? The Bonferroni correction is a very conservative approach, which works well with a small number of comparisons, however it can be used as a guide. The correction considers each comparison separately but with a fraction of the desired confidence level. In the example above, the desired significance level is 0.05 and the number of comparisons is 131072. The Bonferroni correction would impose a significance level of 0.000000381 (0.05/131072) on each comparison to ensure that the overall significance level for the entire dataset was maintained at 0.05. In practice this is much too conservative. Each functional imaging software package has developed methods to address this issue. For example, SPM (http://www.fil.ion.ucl.ac.uk/spm/) invokes random Gaussian field theory to determine levels of significance and the thresholds that should be applied. Still other packages have utilized Monte Carlo simulations to develop a model of the noise of random data to guide the significance level and the thresholds that are applied to the data. Still other packages (Medx, Brain Voyager) provide the user with an array of methods that can used depending on the operator's choice. The bottom line is that the theory used to determine the level of significance applied to the data has to be known by the user in order to interpret the data responsibly.
Another factor that goes into the threshold calculation used on the statistical maps is the amount of connectivity imposed. There should be more significance put on 5 voxels that are suprathreshold and connected than 5 individual islands of activation, where the latter has more potential to be noise. Now this is a tricky situation especially when there is a large amount of smoothing employed. You can imagine that if you take a single voxel with a large amount of activation and smooth it to cover a large area, the surrounding voxels will benefit. This could boost these voxels above the threshold and connect them so that they meet the cluster threshold. The areas of activation could be enlarged. Again it is important that the investigator understand how the data are preprocessed and what impact that will have on the generation of the statistical maps.
Finally a statistic with an appropriate threshold needs to be applied to the data. In some functional imaging software packages, a model of the paradigm is constructed and used to fit the data. For example, it is known when the active and rest periods occur. This can be combined with a hemodynamic response to simulate the expected response from the vascular-based signal. Accounting for motion correlated signal changes, physiologic induced variations, global signal changes over time and even difficulty of the task presented can further enhance the model. The more information supplied to the model the better the estimate of the activation induced signal changes. However, if the model is designed with poor information, it is possible to increase the error, resulting in poor activation detection or, worse, activations that are not related to the paradigm which may lead to misinterpretation. The opposite end of the spectrum in functional MR data processing is to not assume too much about the signal response. This is advantageous because it can apply to any type of experiment, but the residual signal has a large amount of variance. This increased variance may reduce the detection of the activation induced signal change or it may have a structure that can be used to improve the model. Nevertheless, there is a large spectrum of functional imaging analysis packages that have been developed by individual labs to accomplish the level of modeling or statistics that are required to attack a specific issue. It is important to understand the assumptions and basis of the software package that is being used to generate the statistical maps.
Once the statistical map has been created, it needs to be displayed. The first step in this process is to make sure that the high resolution anatomic data (2D or 3D) is coregistered with the functional data. The coregistration step ensures that the underlying anatomy is the source of the activation induced signal change. Next the functional data needs to be interpolated to the resolution of the anatomic data. Typical dimensions of the fMRI data are 64x64 and the anatomic data are 256x256. Bilinear interpolation can be used to increase the matrix size of the fMRI data. This operation takes into account the neighboring data values when creating the data to fill-in between the original data. Sinc interpolation is a better method but requires more computational effort since it incorporates more of the data when determining a voxel value. A larger problem is what to do in the slice dimension. Usually, the fMRI data has a thickness of 3-4mm whereas the anatomic data is typically 1mm thick. Simple duplication of the statistics over 4mm of the brain may not represent the true activation location or continuity. The use of trilinear or sinc interpolation is preferred to ensure that the continuity of the functional activation is preserved.
With the functional significance maps overlaid on the anatomic data, how does one demonstrate activation? This is accomplished by manipulating the look up table (LUT) for the image. Typical MR images are displayed in grayscale with black=minimum value (1) and white=maximum value (256). The user controls the window width (contrast) of pixel values displayed and the center level (brightness) of the window to create the desired image contrast. In the case of fMRI data, the LUT is modified such that the lower portion (1-192) corresponds to the typical grayscale LUT for the anatomic image. The upper portion (193-256) corresponds to the area where the functional information is stored (Figure 16).
Figure 16
In most cases this means that a color LUT is used for the upper band. Typically the hot metal colormap is used such that regions that just meet significance are red/orange, higher levels of significance are yellow and the most significant are white. The cool metal colormap is used for negative activations. The 3D dataset with the modified LUT can be viewed with any 3D package to allow you to slice away the brain to expose the underlying activations or render the brain to visualize the activation patterns on the surface of the brain (Figure 17).
Figure 17
The use of diffusion information for the detection of stroke has been a mainstay of Neuroradiology. Diffusion is the random movement of hydrogen ions in water molecules. Using magnetic gradient fields, it is possible to encode micron scale movements in the large (2 x 2 x 4mm) imaging voxel. However, the diffusion information used to develop a diagnosis is just a small part of the full story that diffusion provides. There is also directional information available. A short, simplified review will be given but for a more complete review of diffusion imaging and post-processing see the references listed. In most clinical applications of diffusion imaging only the magnitude of the diffusion coefficient is displayed, since a region of brain that has suffered a stroke has a slowing of diffusion. In this case a "trace" or isotropic Apparent Diffusion Coefficient (ADC) is calculated and displayed to remove any directional information about the diffusion (Figure 18).
 Figure 18 |
 Figure 19 |
Fluid-like materials such as CSF will have a high ADC value (bright), normal brain tissue will be less but larger than regions where water is restricted, such a stroke (dark). Furthermore, the diffusion may be restricted due to structural confinement such as the case in white matter fibers. In the white matter the hydrogen nuclei have certain directions of preferred travel, such as along the direction of the axonal fibers (Figure 19). Moving in this direction is easier since going perpendicular to the fiber would cause a collision and restrict movement. It is this anisotropic diffusion information that diffusion tensor imaging (DTI) is able to quantify and display. In DTI, many gradient directions are probed to gain a better understanding of the anisotropic diffusion. The minimum number of directions is 6 while some investigators use thousands. It is a matter of imaging time and the goal of the experiment. With the multi-directional diffusion data, it is possible to calculate the ADC value for each principal axis of an ellipsoid for each voxel. If the voxel contained CSF and had isotropic diffusion, the ellipsoid would have equal length axes and describe a sphere (see Figure 19). However, in the case of white matter fibers the result will be an ellipsoid oriented along the fiber direction. Attempts have been made to reduce the directional anisotropy information into a single value The most popular has been described as the fractional anisotropy (FA). The FA is a value that describes the normalized standard deviation of the axes of the ellipsoid and has values between 0 and 1. A value close to zero is purely isotropic diffusion and a value of one is purely along a single direction. With DTI, it is also possible to determine the direction of the main axis of the ellipsoid and encode that in the image as well. This is typically done using color. An example of this can be seen in Figure 20, where red is left to right oriented fibers, green is anterior to posterior oriented and blue represent fibers oriented through the plane of the image (superior to inferior).
Figure 20
This type of image provides new information and can show the disruption, compression or rerouting of white matter fibers. Finally, the directional information can be used to create images of white matter tracks. By selecting a region of interest as a starting point, a visualization algorithm can use the specified region of interest to start the tracking algorithm to determine which voxels are connected. These algorithms use the directional information of the Diffusion tensor as well as the FA value. An example of DTI tractography is shown in Figure 21. In this example a single fiber has been launched from the corticospinal tract. The 3D fiber tracks are shown projected onto a background b=0 image, which is T2-weighted. The tractography demonstrates the ascending fibers as well as fibers that cross over to the contralateral hemisphere. The seed region is shown in yellow.
Figure 21 - - 3D view of DTI tracts projected onto a coronal b=0 image. Tractography image from software provided by MGH.
Final comment
MRI is ripe with examples of image processing, from the complex acquisition and image formation of the MR signal to the 4 dimensional (space and time) functional imaging. It is important to have a solid understanding of the methods used to develop these elegant images but it is critical to understand the pitfalls and limitation as well.
Diffusion References
Copyright © 2004 American Society of Neuroradiology, www.asnr.org