A Novel Approach to Spectral MRI
Tiffany A. Fetzner

Introduction

           Magnetic resonance imaging (MRI) is a well-established non-invasive, diagnostic medical imaging technique based on the nuclear magnetic resonance (NMR) phenomenon. Every tissue in the body has a specific chemical makeup and thus the amount of hydrogen NMR signal differs from tissue to tissue. MRI allows the anatomy inside the body to be seen in either tomographic images taken along at any angle through the body, or three-dimensional volume images. The NMR information present in each pixel of one of these images is both temporal and spectral in nature. Conventional MRI relies on differences in a weighted average of this spectral and temporal information from tissue to tissue to make a diagnosis. If extracted, the spectral and temporal information can be used to facilitate a diagnosis.

           Because the identification of pathology and disease is a main goal of MRI, it is important to examine new ways to facilitate this process. This study examines the feasibility of extracting the NMR spectral information encoded in a magnetic resonance image in a way which is not known to have been attempted or researched before. The approach is based on the variable bandwidth (VBW) imaging technique, which spatially separates the spectral information in an image. In conventional VBW imaging, this separation is considered an undesirable artifact. Since this separation is equivalent to data projections in a spatial-spectral (S-S) domain, an inverse Radon transform (IRT) can be applied to recover spatial images of the different spectral components. The feasibility of this technique was tested with synthetic images.

Background

             The most common magnetic resonance (MR) image is one obtained from a spin-echo imaging sequence. This sequence records images as a function of two instrumental parameters called the echo time (TE) and the repetition time (TR). The hydrogen NMR signal from a volume element (voxel) in the body is represented as a pixel intensity in the resultant MR image. The signal from a single voxel in an imaged tissue is given by Equation 1, where T1 and T2 are two unique temporal properties of a hydrogen type. 

S(TR,TE)=Siri(1-e-(TR / T1i ))(e-TE / T2i ))        Equation ( 1)

The quantity r is the density of spin type i in a voxel. The summation is over the various spectral components giving contributions to the signal from the voxel. Equation 1 implies that spectral and temporal information is summed together in an image.  Different tissues exhibit different signal intensities, but there is no known way to determine from a single image how much of the hydrogen signal is coming from each type of tissue. The T1 and T2 temporal information can be extracted from a series of images recorded with different TR or TE times (1).  The spectral information is typically extracted by chemical shift imaging techniques (2).

            In a spin-echo image, The spectral information about an imaged object is most visibly contained in an image artifact referred to as the chemical shift artifact (CSA).  The research to validate that the CSA of an image contains spectral information about the image was performed by Kaldoudi (3). Chemical shift is the NMR spectroscopists way of describing a variation in the resonance frequency of a nuclear spin due to the chemical environment around the nucleus. Chemical shift is often denoted in the literature as d. The CSA appears in a MR image as pixels that are misrepresented. For example, the pixels representing fat or water are typically offset from each other by two pixels. The chemical shift artifact increases as the bandwidth (BW) of the imager increases.

            Owing to the digital filtering techniques used in modern digitizer circuitry, the BW of the imager may be set in software. The BW of an imager is equivalent to the quadrature sampling frequency after filtering and decimation. Possible bandwidths are given by Equation 2 where fs is the maximum quadrature sampling rate and n an integer. 

BW = fs / n                Equation ( 2 )

For state-of-the-art imagers, fs is 65,536 HZ, making possible BW values of 65,536,  32,768,  21,845,  16,384, ..., 257 Hz.
 
In an MR image, the Field of View  (FOV)  is related to the BW, by Equation 3 where g is the gyromagnetic ratio of the hydrogen nucleus (42.58 MHz / Tesla) and G is the magnetic field gradient. 

FOV = BW / (gG)                    Equation ( 3 )

            To maintain a constant FOV while changing the BW, G must be scaled with BW. Changing the BW on an imager while keeping the FOV fixed is accomplished with variable bandwidth (VBW) imaging.

          Spatial-spectral imaging (S-SI), and VBW imaging, were developed to help eliminate imaging artifacts  produced in MRI scans. The use of  VBW imaging, produced a better signal to noise ratio, but it was learned that as BW decreased the CSA increased in a non-linear fashion (3). In spatial-spectral imaging it was determined that hydrogen could be selected both spatially and spectrally with a single pulse (2); however, fat and water peaks were not easily distinguishable. Both techniques have been used to increase the image resolution of objects, but neither technique by itself can extract the spectral information desired.

            Another useful method is backprojection. Backprojection is an iterative reconstruction technique used in many medical imaging fields (4). One scientist who has done a lot of research on backprojection imaging in spatial-spectral domains is Dr. Gareth R. Eaton. His main interest is in electron paramagnetic resonance (EPR) imaging, not MRI, however, his research indicates that filtered backprojection can be useful in identifying objects in a spatial-spectral domain. Eaton’s technique (5), sometimes referred to as the projection slice algorithm (PSA) in the literature (5), is similar to, if not the same as, in some cases, the inverse of a Radon transform (6). It was the inverse Radon transform (IRT) that was used in this research.

            The inverse Radon transform, (IRT) is named after Johann Radon, in honor of his 1917 paper entitled, "On the Determination of Functions From Their Integrals Along Certain Manifolds"(6).  Filtered backprojection  is used as another name for the application of the IRT, with a Ramlak filter using a cubic interpolation of -0.5 throughout this research. This research extends spatial-spectral imaging, to spatial-spatial-spectral imaging (S-S-SI) by applying filtered backprojection to VBW imaging. Dr. Eaton alluded to a probable link between the spatial-spectral technique and the method of backprojection, in his book on EPR Imaging but concluded that the connection has not yet been investigated (5), (7).

            The three preliminary questions that framed this research were as follows. 1.) In a VBW image, would the characteristic fat and water peaks be resolvable? 2.) Would this new technique be able to resolve the distinctive CH2 and CH3 peaks of fat itself? 3.) Could metabolites actually be resolved using this technique? It should be possible to obtain spectral tissue information, satisfying some of the above three questions, by performing an inverse Radon transformation, on a set of variable bandwidth, magnetic resonance images.

Theory

            The design of this research was based around the development and implementation of a spatial-spatial-spectral (S-S-S) algorithm that used the inverse Radon transform and VBW images and extracted chemical shift information. The Figures below demonstrate pictorially how the algorithm was designed.
 
            The various steps of the ideal procedure are described verbally and pictorially next. This description assumes an object with an NMR signal located in the xy-plane. The NMR signal at any xy-location has chemical shift spectral information which can be described as a function of d. Therefore, x,y, and d define a three-dimensional data set for the spins.

 Step 1.  A set of 256 VBW MR images, as depicted in Figure 1, is read into the computer. Each VBW image has dimensions x and  P(q,y,d), where P(q,y,d) represents a projection through the yd-plane through the object at an angle q. Each image represents a different q.  

Figure 1. Original experimental image data: x vs Pq y d is a projection through the d, y space at angles, q, where q is the angle of the y axis towards +d. Pq y d vs qi.
 

Step 2. A constant x column of data is extracted from each input image. The extracted columns form a q P(q,y,d)-plane of data.

Figure 2. Extracted image column from original data.
 
 
Step 3. The extracted  q P(q,y,d )-plane of data is then used as input for the IRT, creating a yd-image. 

Figure 3. Inverse Radon transform (IRT) procedure.
 
 
 Step 4. The above procedure is carried out for each x location, forming a series of yd-images. 

Figure 4. A series of image slices, the output of the IRT, as a function of spatial coordinate y, and spectral coordinate d.
 
 
Step 5. The images of Step 4 now become an xyd-image. 

Figure 5. Remapped image coordinates along a different plane as a function of spatial coordinates x and y.
 
 
Step 6. Chemical shift information may be extracted. 

Figure. 6 An example of the extracted chemical shift information from the transformed data.
 
            The algorithm design, as modeled by the above Figures, was the conceptual design used for this research.  Not all of the conceptual steps outlined above were implemented. Only the backprojection step was implemented, as no insight would be gained by applying the procedure to 256 x locations. The written S-S-S program algorithm reads in a VBW MR image, or a test image, and computes the RT of the image, to simulate the column extraction step in Figure 2. The RT is then convolved with a Ramlak filter and the result is put through the IRT, with cubic interpolation. The remapped image of Figure 4, is then obtained. The algorithm models of Figures 1, 5 and 6, have not been obtained exactly as is outlined in the diagrams above, due to time constraints. Most of the research done on this technique concerns the results that were obtained from the series of single image preliminary tests performed. The algorithms were implemented in IDL; the programming software from Research Systems Inc., (Boulder, Colorado). The programs and IDL code used in this research can be found in the Appendix.

            Because this work is not known to have been done before, how effective the IRT would be, given the discrete set of possible bandwidths available, became a primary research question. The IRT is lossy, so some information is lost whenever it is used. It is also known from Radon's work (6), that, the utility of the IRT depends on the number of angles that are capable of being projected about an image. Since the inverse Radon transform consists of both backprojecting and filtering, the Riemann Sum procedure of IDL was modified, to accommodate the necessary filtering procedures. A Ramlak filter was used during the convolution process, and the IRT was calculated using cubic interpolation with a value of -0.5 (8).  By itself the Ramlak convolution filter is not a part of IDL. The formula for the Ramlak filter can be found in many image processing, texts. The numbers for the Ramlak filter that were utilized for this research came from Easton (9).
 
           Eaton found that backprojecting  with  6% missing data produced usable results (5). Both the Radon, and its inverse, ideally require that the number of projections used, span the field of view (10). This implies that for a 256 by 256 by 16 bit MR image obtained in a raw data  image format, at least 256 projections need to be obtained per image. A further implication for applying this technique to a full data set, is that the data set itself, needs to be square. An input set of VBW images must consist of at least 256 images. The means that the necessary total amount of input and output data, per set of images, is around 70 Megabytes.

Methods

            Exact implementation of the theory will require projections or bandwidths which are not attainable on commercial MR imagers. In order to answer how effective the IRT would be, given the bandwidth limitations of the GE Signa MR imager, two different algorithmic implementations were examined. The first involved repeating a possible projection angle for an unobtainable BW. This algorithm implementation is referred to in Appendix A as the "projection repetition" IDL code. The second implementation involved simply leaving out any projection whose BW was unobtainable. This form of implementation is referred to in Appendix B as the "missing projections" IDL code. Regardless of the implementation method used, the theoretical ideal case with all necessary projection angles was compared to the experimentally discrete set of bandwidths available.

            The theoretical set of ideal bandwidths has a number of elements equal to the number of projections required to span the FOV divided by the sine of the projection angle being used. In this study, the number of projections was arbitrarily set at 256. The number of projections can be changed by simply modifying the "target_degrees" variable. This variable just represents the number of arbitrary projections to be calculated over a 180 degree span, because the IRT gives duplicate information with just a sign change after 180 degrees.  The theoretical ideal bandwidths (BWT) were calculated for 256 q radian increments over a 180 degree span. Each projection covered,  0.0122718... radians or (0.703125 degree) steps. An example of this is 256 projections divided by the sin of 0.0122718... radians, using Equation 4 with q in radians, is 20,861 Hz. 

BWT = 256 / sinq   Equation ( 4 )

            All of the other variables were hard coded based on whatever the "target_degrees" variable is set at. A semicolon character indicates a line which is ignored by IDL when a program is executed. By removing the semicolon, inactive lines of commands again become active, allowing portions of the programs to be tested easily.  Some sections, in each implementation, are preceded by semicolons because, execution of those commands was only used in determining some of the values for the hard coded variables. These sections, however, reveal how the hard coded variables were obtained, and were intentionally left within the different implementations, to assist anyone who may use them, as part of future research..

           The IDL program used to construct the 10 Color Phantom is found in Appendix C. The 10 Color phantom consists of 10 different colors that were chosen at random to represent ten different tissues in the body. The phantom is similar to the Shepp-Logan phantom (11), and it was meant to be similar, but it should not be confused with the Shepp-Logan phantom.  Because the 10 Color phantom used is a non-standard object, the locations of every frequency spike comprising the phantom is listed in the IDL program.
 
            The IDL program code used to read-in and write-out raw data images can be found in Appendices  D &  E. In these IDL programs, the .skw file extension indicates a raw data file. All of the algorithms, were written to accept the data format already being used in the MRI Lab at the Rochester Institute of Technology.  It is the RIT MRI Lab format that the volunteer's brain MRI data was obtained in, so all test cases were written in this, .skw or raw, data format. However, a raw data file can also be saved as a .tiff image, by using programs like Adobe  Photoshop (12). This technique can read in .skw, or .tiff images, in order to accommodate manipulation of a raw data file. But, to prevent a loss of data from conversion, it is recommended that the raw data format be used. The image format can be selected, by adding or deleting a few semicolons, prior to implementation.

            As with all IDL code, compilation was performed prior to execution. Compilation was done using IDL's compile icon (8).  Execution was done, by typing the name of the procedure code, on the IDL command line. The method of  "projection repetition", was executed by typing Exp_Rie, while the "missing projections" method was executed by typing Exp_Rie_Missing. Both implementations were done after ensuring, that the dependent programs of Appendices  D &  E, had already been actively compiled.

            Upon execution, of one of the implementations, IDL prompts the user for the file name of the test image. The test image is normally read in as a raw data file, but it is possible to read and write image files as either, .skw or .tiff, file formats. IDL then prompts the user to specify a file name for each generated output image. The IDL log informs the user of which output image, I2 through I8, is currently being saved during execution.. Since the images are not generated in numerical order, the IDL log helps ensure that the images are given appropriate file names. The IDL log also provides image statistics to the user.

              The S-S-S algorithm developed for this research had to first be tested on a number of test cases. The input to the algorithm for each of the four main test cases, involved single images instead of a series of images. The purpose of the test images was to determine what d image resolution could be expected, in the final output images, after using the IRT. The test cases included a 3-object test case consisting of five different rectangles, an MR image of a healthy 43 year old male volunteer, a ten color phantom object that was made up for this research, and a single centered pixel to measure a point spread function. The output of the algorithm, for each test case, was a series of seven additional images and is explained further in the outline below.


Outline of Preliminary Test process:

            For the preliminary tests, with only one x-dimension input image, instead of 256; the Radon transform (RT) was calculated first. The result was then convolved with a Ramlak filter and then the image was reconstructed using the IRT with a cubic interpolation of -0.5. The reconstructed image after the IRT was then compared to the original test image. A difference image was then obtained. This process was repeated twice once using only experimentally possible bandwidths, and once using the theoretical ideal, or the result that would be obtainable if any BW was obtainable. The theoretical difference image, and the experimental difference image were then compared to each other, and an overall difference image was obtained.

3-object image I1
Figure 7:  The 3-Object  Test Image

For the 3-Object synthetic test image was a y,d image containing three spatially rectangular objects  at specific y-values composed of five different spectral chemical shift (d) components. It was used to determine how feasible image reconstruction through filtered-backprojection constrained by the discrete set of available bandwidths on a conventional magnetic resonance imager. This same feasibility goal existed for each of the test images examined.

 Step 1.  First, the y,d image (I1), was put through a theoretically continuous Radon transform (RT) capable of imaging at any specified bandwidth, for every q from 0 to 180 degrees. The resulting image of the theoretical RT, or (I2), is a function of (y,q,d) horizontally and q vertically.

Step 2.  Second, this theoretical RT image, (I2), was transformed back into a y,d image by using a filtered inverse Radon transform (IRT). The filter used with the IRT was a Ramlak filter.

Step 3. This new output image, (I3), was then subtracted from (I1) the original image. Each i,j difference element was then squared to remove negative values as defined by Equation 5. This image was defined as (I4). 

            I4ij = (I1ij-I3ij)2     Equation ( 5 )

Step 4. Now the process was repeated, but this time a discrete RT was used. Only the bandwidths, or q angles from 0 to 180 degrees, obtainable by the MR imager were permitted. The resulting image (I5) was then transformed via the IRT and the discrete output image (I6) was obtained.

Step 5. Images I7 and I8, as defined by Equations 6 and 7, were calculated.  These images were used to assess the image degradation that would accompany any practical application of this technique. 

            I7ij = (I1ij-I6ij)2     Equation ( 6 )

            I8ij = (I4ij-I7ij)2     Equation ( 7 )

            The aforementioned process was repeated using both the "projection repetition" implementation, and the "missing projections" implementation. This is what led to the images labeled in the results section as I2 (theoretical RT), I3(theoretical IRT reconstruction), I4(the total square difference image of the theoretical reconstruction compared to original ), I5(experimental RT), I6(experimental IRT reconstruction), I7(The total square difference image of the experimental reconstruction compared to original), and I8(The total square difference image, comparring difference image I4, to difference image I7). The sum of all the pixels in images 4, 7, and 8; Sij I4ij,  Sij I7ij, and  Sij I8ij respectively; were calculated for comparison purposes and reported in  Table1.
 

Results

            The most important results obtained from the preliminary tests can be summarized as follows. Based on the resolvability of the 3-object test Image, this technique is capable of resolving at least 25 Hz, for both the missing and repeated angle techniques. This means that the CH3 and CH2 peaks within fat itself can be resolved since the usual separation is about 0.4 ppm or 25 Hz. In addition, the characteristic fat and water peaks are also chemically resolvable, because they occur at a separation of 220 Hz apart. What has still yet to be determined is how well this technique can separate metabolites? Further research is needed in this area.

            Another important result obtained from these preliminary tests is the discovery that the most effective implementation method, using repeated projections or using missing projections depends on the geometric complexity of the object. In all test cases except for the brain image test case, the missing projection method was more effective. In the case of the brain image however the projection repetition method was more effective. The extent of this geometric dependency is unknown at this time. Further investigation should be done here.
 

Results Part 1 Tabulated Preliminary Test Data

Table 1 Data collected from preliminary test cases* 
Method Used
Test Image 
Sij I4ij/N
Sij I7ij/N
Sij I8ij/N
projection repetition 
3-object
346.219 
4762.070
 4448.502
projection repetition
10ColorPhantom
 252.951
3382.492
 3161.377
projection repetition
Brain_0rig
681.021
743.930
61.659
projection repetition
Psf 
 837.686
837.827
 0.18700
 
 
 
 
 
missing projections
Missing_3object
 346.219 
13014.252
 2547.256
missing projections
Missing_10C
 252.951
 2506.821
2157.211
missing projections
Missing_Brain_0rig
681.021
1228.557
199.448
missing projections
Missing_Psf 
 837.686
837.752
 0.09474
            * N= a scale factor of 65,536

Results Part 2 (Output Images for the Preliminary Test Cases)

 The following images are the resulting output that was obtained, from the Preliminary Test Cases performed, and whose statistics were described  in Table 1. 

3-object Images

Test image 3Object with repetition A Test image 3object with repetition B
Figure 8 Rep. 3-objectA                             Figure 9 Rep. 3-objectB

Test Image 3Object with missing projections A Test Image 3Object with Missing projections B
 Figure 10 Missing 3-objectA                        Figure 11 Missing 3-objectB

 Brain Images

Test Image Joe's Brain with repetition A Test Image Joe's Brain with Repetition B
Figure 12 Rep.BrainA                            Figure 13 Rep.BrainB

Test Image Joe's Brain with missing projections A Test Image Joe's Brain with missing projections B
  Figure 14 Missing BrainA                        Figure 15 Missing BrainB

 Point Spread Function Images

Test image Centered PSF with repetition A Test Image Centered PSF with repetition B
Figure 16 Rep. PSF A                            Figure 17 Rep. PSF B

Test Image Centered PSF with missing projections A Test Image Centered PSF with missing projections B
  Figure 18 Missing PSF A                       Figure 19 Missing PSF B

10 Color Phantom Images

Test Image 10 Color Phantom with missing projections A Test Image 10Color Phantom with missing projections B
Figure 20 Missing 10C A                          Figure 21 Missing 10C B

Test Image 10Color phantom with repetition A Test Image 10Color phantom with repetition B
 Figure 22 Rep. 10C A                              Figure 23 Rep. 10C B

Discussion of preliminary tests.

            The test images show the effectiveness of the IRT image reconstruction. The 3-object test image consists of five rectangles identical in area, but different in grayscale value. Each rectangle was resolved, even with using a discrete set of bandwidths. Each rectangle was comprised of about 22 pixels vertically by 60 pixels horizontally. There are imaging artifacts present in the experimental reconstruction, but the artifacts do not severely hamper the resolvability of the different rectangles. This test image was designed specifically to determine if different components, with some aligned in one direction, could be separated, and distinguished in the reconstruction. This was also the test image used to determine the overall known resolution of this technique. The technique was found to resolve at least 25 Hz., but the narrowest BW resolution is still unknown, at this time.

            The resolution of the 10 color phantom, indicates that only seven of the ten 'color areas' were resolvable. The only reason color was used, in the ten color phantom was to compare the resolvability of ten different 'color' intensities. This image was used as another resolution test. 

Ten Color phantom
Figure 24 The 10Color Phantom

There are two 'green' colors in the center of the test image, one outlines the main green area which is filled with a lighter but intense green hue; both 'greens' and the 'magenta', color were not resolvable. The reason for using the two 'greens', and the 'magenta' in the middle was to test both the resolution of the technique, and to observe any errors that resulted when the test image was observed in grayscale. The 'greens' and 'magenta' correspond to a horizontal oval in the lower half of the test image. The overall overall shape was distinguishable, but its three comprising areas were not. This part of the object has errors in its reconstruction. Another error occurred with the 'tan' ring that offsets the phantom from its 'white' background. The inner color beyond the ring is a 'pale yellow', but there was no distinguishing the 'tan' ring from the 'yellow' fill color. Beyond these errors the object was reconstructed, reasonably well. The irregularity in the geometrical shape of the 'gray' area came through well, the 'blue' dot, in the center of the 'black' circle was resolved correctly. The object is distinguishable from the background, and the 'light blue' circle was apparent. In comparing the two methods of reconstruction the 'missing projections' implementation had more image resolvability, but the errors mentioned previously occurred in both implementations.

            The knowledge that this technique does produce some resolution errors indicates that there are limits to the effectiveness of this technique, but those specific boundaries are not yet known. In an effort to determine the point spread function of a single point a third test image was used. In this case both the projection repetition, and the missing projection method came close to matching the theoretical ideal. The standard deviation of the ideal was 837.686, the 'projection repetition' method had a result of 837.827, while the 'missing projections method' had a result of 837.752 a difference of only 0.075. However, the missing projections method was the more effective of the two implementations.

            The indications of these three different test image cases, tends to suggest that the use of missing projections in a data set gives the best result. But, when the test image used was an MR image of a volunteer's brain the opposite was shown to be true.  The entry in Table 1 is in red to point out this contradiction to what appears to be a natural trend. The ideal standard deviation across the MRI image was calculated to be 681.012, the repetition method gave a result of 743.960, while the missing projection method gave a result of 1,228.557.  The features in the I6 image of the brain with repetition, are also more distinctive in appearance, than the I6 image where missing projections were used.
 
            The standard deviations obtained, in the results section overall were not surprising. A standard deviation across an image that is reconstructed by backprojection is on the order of 107 when no scaling factor is used (13). The data from the above results indicates that a change in the effectiveness of the algorithm exists depending on the geometry of the object being imaged. The more complicated the image the more effective repeating a previous projection instead of leaving the projection out is. For simple images such as the PSF of  a single pixel at the center of an image, allowing the missing projection angles without compensation gives a better result. For geometrically complex images, the repetition method is capable of  providing a more accurate result. A general statement of this trend is that repetition adds error to simple images, but in more geometrically complex images the error introduced is less due to cancellation effects. This is a result which could be explored in more detail, by future investigators.

            In all cases and implementations a Ramlak filter was used with the IRT in this research. It may be possible to enhance this technique further by using other filters with the same process. This is another area which future investigators can extend this research. Bilinear interpolation instead of cubic interpolation, is another option; but with a Ramlak filter cubic interpolation, gave better image resolution, but it is slower computationally.

Conclusions

            This thesis represents a proof-of-concept study to obtain NMR spectral information from variable bandwidth magnetic resonance images.  The technique itself is viable at this point for separating fat from water, and resolving the two characteristic peaks within fat itself. The resolution of this technique is at least 25 Hz. Answering whether or not this technique can successfully resolve metabolites is still unknown at this point.

            The implemented algorithms have a few limitations. The algorithm only accepts images in the Rochester Institute of  Technology MRI Lab 256 x 256 16 bit integer format. The algorithm will either leave out any angle that the imager cannot obtain, or it will repeat the projection from the last possible angle, depending on the implementation used. The algorithm can easily be modified to accept as many projections as specified by the user, with the constraint that the number of projections for  y,d,and q must be equal  to obtain an accurate output result. The spectral resolution of this technique tends to depend on the geometry of the object being imaged. 

Table of Contents