Magnetic Resonance Imaging - Research Paper Example

Magnetic Resonance Imaging Table of Contents 1.Introduction 1 1.1 Tomographic Imaging 2 2.Problem Definition 2 3.Methodology 3 4.Objectives 3 5.Components of the MRI Scanner 3 5.1 The Magnets 4 5.2 The RF System (Radio Frequency System) 5 5.3 The Gradients 6 6.MRI Scanning Process 7 6.1 Image acquisition 7 6.1.1 Tissue magnetization 8 6.2 Image Reconstruction 13 6.3 Fourier Transform 14 7.Factors affecting Image quality in MRI 15 8.MRI Artifacts 16 9.Parallel Imaging 17 9.1 SENSE 17 9.2 GRAPPA 21 10.Conclusion 25 11.Appendices 26 11.1 Appendix 1: SENSE Code 26 11.2Appendix 2: GRAPPA CODE 30 12.References 41 1. Introduction Magnetic Resonance Imaging has grown to be the primary technique throughout the body in the practice of performing diagnoses. It is fast replacing the once much hyped Computed Tomography (CT). This may be due to the perceived advantages of MRI being non-intrusive, using non-ionizing radiation with a high resolution of soft tissues. MRI also has the capability of providing images in both two and three dimensions. MRI is able to provide images containing morphological (surface information) and functional information. The obtained image is based on multiple parameters, any of which can modify tissue contrast. MRI can generate thin-section images when the body is exposed to an electronic field. It creates a strong magnetic field which magnetizes the hydrogen atoms contained in the body tissues of human beings. The MRI creates a steady state of magnetism in the body within the human body. It then stimulates the body with radio waves to change the induced magnetism of the body. MRI then stops the radio waves and registers the electromagnetic transmission of the body. Finally, the MRI concentrates the transmitted signal then uses it to reconstruct the images of the internal organs of the body by computerized axial tomography. The use of the MRI technology has not come without challenges. The machine produces a lot of noise once switched on and patients are encouraged to put on earmuffs to minimize the discomfort the noise can cause. Furthermore, the machine takes significant amount of time to produce a clear and easy to understand image. Research is currently underway to produce machines which can produce clear images and yet minimizing the time taken to produce such images. This factor is giving physicists and neuroscientists sleepless nights. The other area that requires research is the accurate scanning of the brain. Being a target in motion, it necessitates for the use of a more refined technique to sample and hence give a better understanding of what goes on in the brain. The MRI machine also requires the patient to remain still during the entire scanning period if a clear image is to be produced. Some patients may be nervous and unable to remain still, hence they may be injected with medicine which can help them relax. All these are simple rules that are to be taken into account in order for a scan to be carried out successfully. More importantly, efforts are being made to minimize the scan time and ensuring not to compromise on image quality. 1.1 Tomographic Imaging This is a computer based teaching package used to provide clear understanding the principles of MRI from both microscopic, macroscopic, and imaging system. Magnetic resonance started as a tomographic imaging technique for producing NMR images of a slice through the human body. The magnetic resonance image is composed of several pictures called pixels, with an intensity proportional to the NMR signal intensity of the contents of the corresponding volume element of the object being imaged. Magnetic Resonance Imaging employs the principle of absorption and emission of energy in the radio frequency range of the electromagnetic spectrum. The misleading perception that imaging using radio waves was impossible made the development of non-ionizing imaging equipment. This myth was however broken by the invention of MRI which produces based on spatial variations in the phase and frequency of the radio frequency energy being absorbed and emitted by the imaged object. 2. Problem Definition MRI equipment are faced with a number of challenges that hinder their effectiveness. Among these are the challenge of scan time and the time it takes for a scanner to produce an image. This becomes a major hindrance to offering health care services since it takes longer to scan a single patient, making it hard to scan many patients as it would be if the machines would have been faster enough. This may make patients die as they wait for their disorders to be scanned. The other challenges are the noise produced by the scanner as it is in use, poor clarity when scanning moving objects, e.g. the brain, with the need of patients required to remain still during the entire scanning process. All these are challenges that prompted me to embark on the research to find ways to solve these problems. 3. Methodology Magnetic Resonance Imaging is a high research field and I had to get a full understanding of the latest developments in the field by reading widely. I consulted printed as well as online books, journals, magazines and newspapers. I also visited MRI research websites to ensure that I came up with a comprehensive and reliable report. 4. Objectives The aim of this project is to propose the improvements that can be done on the MRI machine to improve the speed while at the same time preserving the image quality, by use of two algorithms, i.e. SENSE and GRAPPA. This research project also aims at putting forth the calculations that are involved in these two image reconstruction techniques, and also to show the differences between GRAPPA, non-linear GRAPPA and SENSE techniques. 5. Components of the MRI Scanner The main components of MRI have undergone several improvements over the years to try and make them more efficient. Since the main issues with the MRI equipment are the time it takes to produce an image and image quality, most changes have been effected on the components of the equipment to achieve these. The diagram below shows the main instrumentation of an MRI scanner: Figure 1: The main components of an MRI scanner Although the appearance of scanners will differ depending with the manufacturer, a typical MRI scanner is generally made of the following parts: The Magnets RF System The Gradients 5.1 The Magnets This is the main component of the MRI machine. Some low field magnets are permanent or resistive, but for MRI machines above 1.0 Tesla, the magnet is usually superconductive. It is wound from an alloy, for instance Nb-Ti, which has zero electrical resistance below a critical temperature. The magnet is enclosed and cooled by a cryogen containing liquid helium, which require monthly maintenance. Poor superconductive windings will make the scanner lose over 5g per year, which causes quench; a situation where the magnet suddenly loses its superconductivity and then begin to heat up. Vents are usually attached to ensure that the magnet loses any extra heat safely. The other type of magnet found in the MRI equipment is permanent magnet constructed from ferromagnetic material. It is larger than the superconductive magnet and can function even in the absence of electricity. It has the advantage that it provides for flexibility in the MRI design. Its size and weight however makes it less popular in addition to that its ability to generate a stable magnetic field is questionable. 5.2 The RF System (Radio Frequency System) The RF system is used for transmitting radiofrequency radiation that causes the hydrogen atoms to emit signals. The RF system then receives the signals that are emitted and amplifies the signals so that it can be at a state that can be manipulated by the computer. The RF coil is the most essential component of the RF system. The coils are created to produce an oscillating magnetic field, which then induces an atom at a designated area to absorb RF radiation which in turn emits a signal. The coils are made in such a way that they are capable of receiving the emitted signal and transmit it to the computer. RF coils may be either saddle RF coils or solenoid RF coils. The coil is aligned as much as possible with the patient to increase efficiency and minimize signal to noise ratio by focusing only on the point of interest. The transmitters and the receivers in the RF system are highly sensitive to signal. 5.3 The Gradients Figure 2: Gradient Magnets showing x, y and z coils The functional role of the gradient coils is to produce linear changes in magnetic field in each of the x, y and z axes. Through combination of gradients in pairs, it is possible to produce oblique images. Gradient specifications are stated in terms of a slew rate which is equal to the maximum achievable amplitude divided by the rise in time. Typical modern slew rates are in the range of 150 T/m-s. the gradient coils are shielded in a similar manner to the main windings to reduce eddy currents induced in the cryogen which would otherwise degrade the image quality. 6. MRI Scanning Process The process from the moment the machine takes an image of the patient to the moment the image is displayed on the computer screen entail complex processing and manipulation summarized into four major steps. These include: Image acquisition K-Space domain Fourier transform Final image 6.1 Image acquisition This is the first step in Magnetic resonance imaging and forms the basis of successful image formation by MRI machine. The patient is made to lie on the table and the machine is switched on. A strong magnetic field is then applied on the patient. There are three methods to generate this field; fixed magnets, resistive magnets and super-conducting magnets. Fixed and resistive magnets are restricted to field strengths below 0.4T and thus cannot generate the higher magnetic strengths required for clear imaging. The magnetic field required must be strong and uniformly aligned in space and stable in time. The magnetic field thus applied causes hydrogen atoms in the body to be excited and thus align with the external magnetic field. A pulse of low energy radio waves is then sent to the body. The aligned hydrogen atoms then re-transmit the waves which are then received by the receiver of the RF system and transmitted to the computer for processing. The ability of MRI to distinguish between different types of tissues is based on the fact that different tissues , both normal and pathologic, become magnetized to different levels or will change their levels of magnetization (relax) at different rates. The time required to acquire images is determined by the duration of the imaging cycle or cycle repetition time (TR), and the number of cycles. The more the number of cycles, the better the image quality. 6.1.1 Tissue magnetization The placement of a tissue in the magnetic field makes the tissue reach its maximum magnetization within a few seconds and it can remain in that state of maximum magnetization unless there is a change in the magnetic field or if there are pulses of RF energy applied. The tissue to be scanned is cycled through changes in its magnetization during the imaging cycle. The direction of tissue magnetization is specified in reference to the direction of the applied magnetic field. The tissue can be subjected to either longitudinal magnetization (the direction parallel to the direction of the field), or transverse magnetization (tissue magnetized at 90⁰ with respect to the direction of the magnetic field) The direction of tissue magnetization can be changed (flipped) by applying a pulse of RF energy and this is done throughout the imaging process. The magnetization angle flip is determined by the duration and strength of the RF pulse. 90⁰ and 180⁰ flip angles are the most common. A 90⁰ pulse applied to longitudinal magnetization flips the tissue to the transverse plane. This reduces the longitudinal magnetization to zero, a condition called saturation. It also produces transverse magnetization which is unstable or excited condition. A 90⁰ pulse applied to longitudinal magnetization produces both saturation of the longitudinal magnetization and a condition of saturation. The quality of the image obtained will depend on the following times of the scanning period: Longitudinal magnetization relaxation time (T1) Transverse magnetization relaxation time (T2) 6.1.1.1 Longitudinal Magnetization and Relaxation If magnetization is temporarily redirected by an RF pulse, the tissue will return to its original longitudinal position over a period of time. If the longitudinal magnetization is considered individually, it regrows after it has been saturated. This regrowth of longitudinal magnetization is the relaxation process, which occurs after saturation. All this is determined by the characteristics of the material and the strength of the magnetic field. Longitudinal magnetization grows exponentially. Figure 3:Spin-echo Pulse Sequence: Single echo T1-weighted The brightness of the image formed in MR imaging is dependent on the level of magnetization. The brightness of the tissue changes during saturation to a dark image which recovers brightness during the relaxation period. Figure 4: Longitudinal relaxation (T1 relaxation) showing image brightness The relaxation time is hard to predict because of its exponential nature. Conventionally, the time required to reach a magnetization of 63% is taken as the longitudinal relaxation time, T1. Factors that affect either the proton resonance frequency or the frequency associated with the molecular motions will most likely affect the values of T1. Small molecules will tend to have a relatively longer T1 values than larger molecules. 6.1.1.2 Transverse Magnetization and Relaxation This is produced by applying a pulse of RF energy to the magnetization tissue. This is achieved by a 90⁰ pulse which converts longitudinal magnetization into transverse magnetization. Figure 5: Spin-echo pulse sequence: dual echo T2-weighted (Hesselink J, 2010) This is an unstable or excited condition that quickly decays after termination of the termination of the excitation pulse. The decay of transverse magnetization is a relaxation process characterized by specific relaxation times, T2. Different tissues have different T2 values that are used to tell tissues apart and adjust the image contrast. Figure 6: Transverse magnetization decay (Sprawl Educational Foundation, 2011) Transverse magnetization is used during tissue formation so as to develop the image contrast based on the differences in T2 values, and to generate the RF signals emitted by the tissue. Transverse magnetization is a spinning magnet condition that generates an RF signal. Table 1: T1 and T2 for various tissues T2 (msec) T1 (0.5T) (msec) T1 (1.5T) (msec) Adipose (fat) 80 210 260 Liver 42 350 500 Muscle 45 550 870 White Matter 90 500 780 Gray Matter 100 650 920 CSF 100 1800 2400 Figure 7: Case A; contrast-enhanced T1-weighted brain MRI shows multiple lesions, with the largest one in the left occipital lobe. Case B; T1-weighted MRI demonstrates a sub acute hematoma in the metastatic lesion (Yan et al., 2010) Figure 8: T1-weighted image (A) and T2-weighted image (B) of the spine Figure 9: Sagittal T1-weighted image (A), Sagittal T2-weighted image (B). Arrow indicates lesion found in the T2-weighted study 6.2 Image Reconstruction 6.2.1 K Space During the image acquisition process, the signals are collected, digitized and stored in the computer memory. This is a configuration process known as k space, and it is divided into lines of data that are filled one at a time. The k space must be completely filled before the image reconstruction process can be terminated. The size of the k space (number of lines) is determined by the necessity for good image quality. Figure 10: Imaging process illustrating the k space (Sprawl Educational Foundation, 2011) The image reconstruction process is usually fast compared to the acquisition process and do not require any adjustments by the machine operator. The relationship between k space data and the image data is Fourier transformation. In 2-D Fourier transform imaging, a line data correspond to the digitized MR signal at a particular encoding level. 6.3 Fourier Transform This is an operation which transforms functions from time to domain. After the k space is completely filled up, it undergoes Fourier transform. There are several reconstruction methods, but the most common for clinical applications is the 2-D Fourier transform. It is generally a mathematical procedure that sorts composite signals into individual frequency and phase components. Since each voxel in a row emits a unique signal frequency, and each voxel in a column a unique phase, the Fourier transformation is able to determine the location of each signal component and direct it to its respective pixel. Figure 11: The concept of signal encoding and image reconstruction 7. Factors affecting Image quality in MRI There are several factors that the MRI operator should be aware of when and before performing a scan as the can affect the quality of the image obtained. Field of view An image consists of field of view that related to the region of interest to be covered. Most MRI equipment have field of view ranging between 10 to 50cm. this mean for instance, if the entire spinal cord is to be scanned in the sagittal plane, its upper and lower parts need complementary series of pulse sequence. The other uncovered parts would transmit signals which will be detected as noise by the scanner and this will lead to the formation of a distorted image. Slice Thickness Giving each slice a thickness needs the excitation of a band of nuclei by an excitation pulse. Increasing the thickness of slices increases the SNR. As the thickness of a slice increases, the resolution decreases and SNR also increases. A thin slice produces an image with high resolution. Number of Excitations Each signal contributing to the formation of an image in MRI can be received once or collected several times using repeated excitations. The average signal values can then be used to generate the image. This increases precision of the image. Image Acquisition Time The longer the acquisition time, the more chances that the patient moves, and the more chances of generation of artifacts that prevent interpretation of the images. The optimum scanning time should thus be maintained so as to produce appreciable image quality free from artifacts. 8. MRI Artifacts Imaging in MRI has come with its own score of challenges. These are as discussed below: Aliasing This usually appears when the diameter of scanned area is greater than the dimensions of the scanned area. Chemical shift artifacts These appear at the interfaces of water and fat since the processional frequency of protons is slightly different in the substances. The system displays them as dark regions of signal void on one side of the water containing tissue and a region of bright signal at the other end of the water fat interface. Truncation Artifact These are bright and dark lines that are seen parallel and adjacent to boarders of abrupt intensity. These artifacts are commonly seen in phase encoding direction. 9. Parallel Imaging Parallel acquisition is a technique that uses the combination of signals from several coil elements in a phased array to reconstruct an image with the aim of improving the signal-to-noise ratio, and to speed up image acquisition and reduce scan time. Parallel reconstruction techniques are instrumental in achieving the main aim of reducing scan time while at the same time, maintaining the image quality. Parallel reconstruction algorithms can be divided into two main groups namely: SENSE (Sensitivity Encoding) – this combines the images from the individual coils to reconstruct the required image (i.e. involves reconstructing in the image domain, after Fourier transform) GRAPPA (Generalized Auto-calibrating Partially Parallel Acquisition) – undertake the image reconstruction in the frequency domain, before Fourier transform (image reconstruction from frequency signals of the individual coils). 9.1 SENSE Stands for SENSitivty Encoding. Each individual coil is covered with a systematically mapped spatial zone. A sensitivity profile is then established for each coil element. Spatial encoding in the direction of the phase encoding gradient is then under-sampled for the purpose of saving time. The intermediary images of each coil element therefore forms fold over artifact. The image is then calculated by deducing pixel value from the intermediary images of every coil element. This makes it possible to measure coil element sensitivity profiles by acquisition before the sequence of imaging as a tree dimensional acquisition in low resolution over the overall field of view or by self-calibration. Calibration enables the reconstruction of the final image, which reduces the effects of noise which would otherwise affect the image quality. SENSE makes use of spatially placed coils. These coils are placed in different positions relative to the object to be scanned. Due to this, the coils receive varied signal strength depend on how far it is from the object to be scanned. The idea employed in SENSE is the application of knowledge of sensitivities of the coil elements to calculate the aliased signal component at each point, then allocate the signals to their actual locations in the unfolded image. Increasing the gap between adjacent k space lines by an acceleration factor, R, makes the signals from R locations to overlap in the image. Reduction of field of view (FOV) can be stated mathematically by relating R-fold FOV to fold aliased image, NA as indicated below: Ij(y) = Where; NA = total number of signals present at location ‘y’, owing to aliasing Ij(y) = image signal J = 0,1,…….Nc – 1 Nc= number of elements in the coil array C = sensitivity matrix Ij(y) are known factors since they are acquired aliased images (one for each coil). are the aliased magnetization values and they are to be calculated. The equation above can be rewritten into a matrix for simplicity. The reconstruction problem in SENSE can be solved by computing a set of linear equations defined below: Where: I = the aliased signals obtained from aliased images C = the encoding matrix (also called Sensitivity Matrix) M = image to be recovered (given by; C-1I) The use of SENSE as a type of parallel imaging technique increases the noise at the expense of aliasing and reduction of scan time. Thus it would be important to increase the SNR so as to improve the image quality. For R = 1, SNR is maximum. The SNR is given by the equation below: SNRSENSE= SNRNORMAL/g Where; = the expected SNR loss that results from reduced scan time by factor R g = geometric factor that represent noise magnification that occurs when aliasing is unwrapped. It’s given by the equation: g = CTψ1C)-1]ii[(CTψ-1C)]ii The performance of the above algorithm can be evaluated using two quantification parameters. These include; the signal to noise ratio (SNR) and the artifact power. Signal to noise ratio (SNR) Before reconstruction, the user is given the option of selecting the region of interest for signal (ROS) and the region of interest for noise (RON), which is usually taken to be the background. The following formula s thus used to calculate the SNR; SNR (dB) = 20 log10 Artefact Power (AP) This concept is derived from ‘Square Difference Error’. A reference image is provided. The AP in the reconstructed image is evaluated on the basis of the reference image. AP = If Ireference= Ireconstructed, then the AP will be zero and this indicate that there is no artefact in the reconstructed image i.e. the reconstructed image and the reference image are identical. SENSE involves acquisition of the full k space and their sum of squares reconstruction used as the reference image. The specified number of k space lines is skipped to produce aliased images, depending on the acceleration factor. The inverse Fourier transform of sub-sampled k space give aliased images. To obtain the sensitivity map, the central lines of k space of the full field of view data are truncated so that a low resolution image of the coils is obtained. These values are then normalized by dividing by sum of square image to give the sensitivity map. Noise also plays a major role on the obtained images as it interferes with clarity of these images. Figure 12: Output from the SENSE showing the sensitivity map, sum of squares image, the reference image and noise-affected image 9.2 GRAPPA This stands for Generalized Auto-calibrating Partially Parallel Acquisition. In this technique, under-sampling occurs over the entire k space. After the sampling process, there occurs missing intermediary k space lines which are calculated from the recorded signals. This is achieved by combining the k space by weighting each coil signal. The above statement can be summarized in the following equation: Sj(ky + rky,kx) = * S1(ky + tRky,kx + hkx), j = 1,…,L,r tR Where; Sj(ky + r ky,kx) = the unacquired k space signal at the target coil S1(ky + tRky,kx + hkx) = the acquired undersampled signal Wi,j(l,t,h) = linear combination coefficients Conventionally in GRAPPA, sampling of the k space of each coil is done at the Nyquist rate in order to obtain auto-calibration signal (ACS) data. The outer k space is under-sampled by some outer reduction factor (ORF). GRAPPA method can be divided into two sub-types name conventional GRAPPA and Non-Linear GRAPPA (NL-GRAPPA). Conventional GRAPPA is the normal one as explained above. To transform the above formula to a nonlinear algorithm, the nonlinear mapping is applied over the acquired k space data. This transformation requires the solving of the following equation; b = Φ(A)x Where; Φ(A) = [Φ(a1), Φ(a2),….,Φ(aM)]T is M * NKmatrix, superscript T is the matrix transpose, M represent the number of acquired ACS data, NK being the dimension in reproducing the kernel Hilbert space (RKHS), usually much higher than M. a’s are row vectors of the matrix A. A kernel is defined as a continuous, symmetric and positive definite function. The variables on the x axis can be solved using the linear algorithm: x = ((Φ(A))H(Φ(A)))-1(Φ(A))Hb The estimated coefficients above are then used to reconstruct the missing data in the outer k space, just like the conventional GRAPPA. Figure 13: Axial brain image showing the reconstructed image and the reference mage. Figure 13 above illustrate the nonlinear GRAPPA in comparison to the conventional GRAPPA. It is notable that the nonlinear method does not produce a clear image in this image due to aliasing. The nonlinear method uses the normal computational algorithm to find the coefficients. It is possible to use other linear computational algorithms like reweighted least-squares and total least-squares methods. It is a normal trend of using polynomials to approximate smooth unknown functions. A polynomial kernel is thus used for Φ mapping. The following polynomial is used for this task; κ(ai,aj) = (γaiTaj + r)d Where: γ and r = scalars d = the degree of the polynomial The first and second order components are both incorporated into a single equation to construct a hypothesis space as shown; Φ(a) = [1,√2a1, √2a2,…,√2aK,ai1, aj1, ai2,aj2…]T, where i and j are chosen randomly from 1, 2,…K, such that the size of Φ(a) is equal to the desired dimension NK. The final nonlinear GRAPPA method is thereby formulated as: Where: P(n) and Q(n)are chosen randomly from Sl(ky + bR∆ky, kx = h∆kx) with l = 1,…,L; b = B1,…,B2; and h = H1,…,H2, and N is the number of randomly selected second-order terms. Some part of the nonlinear GRAPPA equation resembles the conventional GRAPPA equation. This is so as to provide the relationship between the missing and the acquired signal, in the absence of noise approximations. Figure 14:MR images acquired using Eight-channel head array It is clear from figure 14 and 15 that nonlinear method achieves a superior image quality than the conventional GRAPPA. The nonlinear method effectively removes noise experienced in the conventional method. Nonlinear GRAPPA also preserves the resolution of the axial image without blurring. Slightly blurred images are obtained when using the sagittal angle, probably due to the trade-off between noise suppression and preservation of resolution. (a) Coil 1 Coil 2 Coil 3 Coil 4 (b) Figure 15: In-vivo MR images acquired using (a) four-channel spine array and (b) eight channel head array Images reconstructed using NL-GRAPPA appear to have more noise than aliasing artifact, probably because NL-GRAPPA trades noise reduction for the removal of aliasing. 10. Conclusion The use of MRI has revolutionized the world of medicine as the methods that were earlier used to view the internal human organs were accompanied by a score of side effect, some more lethal than the disease being scanned. The use of MRI for scanning is totally harmless and the images produced are easier to understand. Although the images would occasionally be distorted due to some factors, the development of SENSE and GRAPPA has enhanced the quality of output image. The development and use of GRAPPA and SENSE has made it easier to improve the image quality while at the same time reducing the scan time. This is especially important in that the patient does not spend more time in the scanner thus reducing the chances of the patient moving. The movement of the patient is the major cause of image distortion. The reduction of scan time also means that the hospital will be able to scan more patients peer day thus increasing income. The advantages of image reconstruction by use of GRAPPA and SENSE can never be understated. 11. Appendices 11.1 Appendix 1: SENSE Code There are three steps: k-space data read-in and simulated subsampling, sensitivity estimation and SENSE reconstruction. function [recon, sucessful_flag, gmap] = sense(channelsensitivity, reduced_kspace_data, FOV_reduction_factor) %Input- % channelsensitivity: size (points_Freq Enc, lines_Phase Enc, % Number_channels), the spatial coil sensitivity function. It is % required that lines_Phase Enc must = Number_Phase Enc_Subsampled * FOV_reduction_factor; % To enforce this condition, the channel sensitivity function may % need to be interpolated, i.e., by zero padding the k-space % calibration lines before inv-FT. % reduced_kspace_data: size (Number_Freq Enc, Number_Phase Enc_Subsampled, % Number_channels); As a convention, the center k-space line is in the % middle. % FOV_reduction_factor: It must be an integer, (FOV of sensitivity function)/(FOV of % reduced_kspace_data). If the data acquisition parameters results in % fractual reduction factor, then channelsensitivity functions must % be cropped and interpolated to make the reduction factor an integer % before calling this function. %%Output- % The sense reconstruction: size (Number_Freq Enc,Number_Phase Enc) % sucessful_flag: 1 if reconstructin is successful and 0 if there is no % error in the process. % gmap: the g-factor map of the sense recon, size (Number_Freq Enc, Number_Phase Enc) %-------------------------------------------------------- % The k-space subsampling is along the second dimension.; The leading dimension is % the frequency encoding; The second dimension is phase encoding; The last % dimension is the channel index; %----------------------------------------------------- sucessful_flag = 0; % will return 1 if reconstructin is successful and there is no error in the process. [Nfe,Npe_seg,Ncoils] = size(reduced_kspace_data); Npe = size(channelsensitivity,2); Rr = (Npe/Npe_seg); if FOV_reduction_factor - round(FOV_reduction_factor)~=0; display('FOV_reduction_factor must be integer'); return; end; if Rr ~= FOV_reduction_factor; display('Sensitivity function and k-space data does not agree with the FOV reduction factor!'); display('You may need to interpolate the sensitivity functions. Type ''help sense'' to get more information.'); return; end; % generate the overlapped_img: (Nfe,Npe/R,Ncoils) % index as [-N/2:(N/2-1)]. Centered, unshifted overlapped_img = zeros(Nfe,Npe_seg,Ncoils); for coil=1:Ncoils overlapped_img(:,:,coil)=ifftshift(ifft2(ifftshift(reduced_kspace_data(:,:,coil)))); end % shift the overlapped image % In the overlapped iamges above, the aliasing occurs as if (Npe/2+1) is % the center of the unwrapped image, and the overlapped_img corresponds to % an reduced FOV with the same center. Therefore, for R=2, the reduced % FOV corresponds to (Npe/4 + 1): (Npe/2). But for R=3, the reduced FOV % corresponds to (Npe/3+1);(Npe/3*2), exactly the same as (1:Npe/3). For % the convenience of forming the SENSE equation, shift the first and sencod % halves of the overlapped image if an even reduction factor is used. % For odd reduction factors, this shift is not necessary. if mod(Rr,2)==0 %even R tmp_img = overlapped_img; overlapped_img(:,1:(Npe_seg/2),:) = tmp_img(:,(Npe_seg/2+1):Npe_seg,:); overlapped_img(:,(Npe_seg/2+1):Npe_seg,:) = tmp_img(:,1:(Npe_seg/2),:); end %unwrap recon=zeros(Nfe,Npe); for pe=1:Npe_seg for i=1:Nfe if((pe+Npe_seg*(Rr-1))>Npe) delt=(pe+Npe_seg*(Rr-1))-Npe; else delt=0; end s1=squeeze(channelsensitivity(i,pe:Npe_seg:(pe-delt+Npe_seg*(Rr-1)),:)); s=transpose(s1); I_mat=squeeze(overlapped_img(i,pe,:)); recon(i,pe:Npe_seg:(pe-delt+Npe_seg*(Rr-1))) = transpose((inv(s'*s)*s'*I_mat)); tmp = pinv(s'*s).*(s'*s); tmp = sqrt(abs(diag(tmp))); [l trash]=size(tmp); if (l ~= length(pe:Npe_seg:(pe+Npe_seg*(Rr-1)))) tmp(l+1:length(pe:Npe_seg:(pe+Npe_seg*(Rr-1))),:)=0; end gmap(i,pe:Npe_seg:(pe+Npe_seg*(Rr-1)))= tmp'; end end sucessful_flag = 1; return 11.2 Appendix 2: GRAPPA CODE function [full_fourier_data, rec_img, coef0] = grappa_function(reduced_fourier_data, ORF, pe_loc, acs_data, acs_line_loc, num_block, num_column, times_comp) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Input parameters % reduced_fourier_data : undersampled k-space data % ORF : outer reduction factor % pe_loc : undersampled phase-encoding lines' location % acs_data : auto-calibration signal data (middle region of k-space) % acs_line_loc : auto-calibration signal lines' location % num_block : number of blocks % num_column : number of columns % times_comp : times of the number of the first-order terms (the number of % the second-order terms = time_comp X the number of the first-order terms) % Output parameters % full_fourier_data : reconstructed k-space (with ACS replacement) % rec_img : reconstructed image % coef0 : coefficients for reconstruction % time_comp parameter % As the parameter time_comp increases, relevant second-order terms are % added for reconstruction. % When time_comp = 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Get dimensions and initialization [d1_reduced,d2,num_coil] = size(reduced_fourier_data); d1 = d1_reduced*ORF; if ORF==3 d1=d1_reduced*ORF-2; end if ORF==5 d1=d1_reduced*ORF-4; end if ORF==6 d1=d1_reduced*ORF-2; end %Decide which lines are possible lines to fit all_acquired_line_loc = unique(sort([pe_loc, acs_line_loc])); combined_fourier_data = zeros(d1,d2,num_coil); combined_fourier_data(pe_loc,:,:) = reduced_fourier_data; combined_fourier_data(acs_line_loc,:,:) = acs_data; ind_first = find(all_acquired_line_loc == acs_line_loc(1)); ind_last = find(all_acquired_line_loc == acs_line_loc(end)); %Form the structure that indicates where lines are fitted line_group = cell(num_block,ORF-1); for s = ind_first:ind_last for mode = 1:num_block for offset = 1:ORF-1 tentative_line_ind = [all_acquired_line_loc(s)-offset-(mode-1)*ORF : ORF : all_acquired_line_loc(s)-offset+(num_block-1)*ORF-(mode-1)*ORF]; valid_flag = 1; for t = 1:num_block if isempty(find(all_acquired_line_loc == tentative_line_ind(t))) valid_flag = 0; break; end end if valid_flag == 1 line_group{mode,offset} = unique([line_group{mode,offset}; [all_acquired_line_loc(s),tentative_line_ind] ], 'rows'); end end end end %Solve for the weighting coefficients fit_coef = zeros(num_coil,ORF-1,num_block,(times_comp+1)*num_block*num_coil*num_column); for jj = 1:num_coil for offset = 1:ORF-1 for mode = 1:num_block fit_mat = zeros(num_block*num_coil*num_column, d2*size(line_group{mode,offset},1) ); target_vec = zeros(1,d2*size(line_group{mode,offset},1) ); for nn = 1:size(line_group{mode,offset},1) temp_data = combined_fourier_data( line_group{mode,offset}(nn,[2:end]), :,:); temp_data = permute(temp_data,[1 3 2]); temp_data = reshape(temp_data, [num_block*num_coil,d2]); fit_mat((num_block*num_coil*floor(num_column/2)+1):(num_block*num_coil*ceil(num_column/2)), [1+(nn-1)*d2 : nn*d2]) = temp_data; target_vec([1+(nn-1)*d2 : nn*d2]) = combined_fourier_data( line_group{mode,offset}(nn,1), :,jj); end transfer_matrix=fit_mat((num_block*num_coil*floor(num_column/2)+1):num_block*num_coil*(floor(num_column/2)+1),:); column_label = [[floor(num_column/2):-1:1],[1:floor(num_column/2)]]; for column_idx = 1:num_column-1 if column_idx = times_comp break; end end if idx_comp >= times_comp break; end end if idx_comp >= times_comp break; end end fit_coef(jj,offset,mode,:) = inv((new_fit_mat.')'*new_fit_mat.')*(new_fit_mat.')'*target_vec.'; end end end clear temp_data; %Generate the missing lines using superpositions candidate_fourier_data = zeros(d1,d2,num_coil,num_block); for mode = 1:num_block candidate_fourier_data(:,:,:,mode) = combined_fourier_data; end for ss = 1:d1 if isempty(find(pe_loc == ss)) offset = mod(ss-1,ORF); for mode = 1:num_block tentative_line_ind = [ORF*floor((ss-1)/ORF)+1-(mode-1)*ORF : ORF : ORF*floor((ss-1)/ORF)+1+(num_block-1)*ORF-(mode-1)*ORF]; if max(tentative_line_ind) = 1 temp_data = combined_fourier_data(tentative_line_ind,:,:); temp_data = permute(temp_data,[1 3 2]); fit_mat=zeros(num_block*num_coil*num_column,d2); temp_data=reshape(temp_data,[num_block*num_coil,d2]); fit_mat((num_block*num_coil*floor(num_column/2)+1):num_block*num_coil*(floor(num_column/2)+1),:) = temp_data; column_label = [[floor(num_column/2):-1:1],[1:floor(num_column/2)]]; for column_idx = 1:num_column-1 if column_idx = times_comp break; end end if idx_comp >= times_comp break; end end if idx_comp >= times_comp break; end end for jj = 1:num_coil candidate_fourier_data(ss,:,jj,mode) = (squeeze(fit_coef(jj,offset,mode,:))).'*new_fit_mat; end else candidate_fourier_data(ss,:,:,mode) = 0; end end end end %Use ACS lines to obtain the goodness-of-fit coefficients gof_coef = zeros(num_coil,ORF-1,num_block); for jj = 1:num_coil for offset = 1:ORF-1 fit_mat =[]; target_vec = []; for ss = 1:length(acs_line_loc) if mod(acs_line_loc(ss)-1,ORF) == offset valid_flag = 1; for mode = 1:num_block if isempty(find(line_group{mode,offset}(:,1) == acs_line_loc(ss))) valid_flag = 0; break; end end if valid_flag == 1 temp_mat = []; for mode = 1:num_block temp_mat = [temp_mat; candidate_fourier_data(acs_line_loc(ss),:,jj,mode)]; end fit_mat = [fit_mat,temp_mat]; target_vec = [target_vec,combined_fourier_data(acs_line_loc(ss),:,jj)]; end end end gof_coef(jj,offset,:) = target_vec/fit_mat; end end %Combine the data from different modes using goodness-of-fit full_fourier_data = combined_fourier_data; for ss = 1:d1 if isempty(find(all_acquired_line_loc == ss)) offset = mod(ss-1,ORF); for jj = 1:num_coil for mode = 1:num_block full_fourier_data(ss,:,jj) = full_fourier_data(ss,:,jj)+gof_coef(jj,offset,mode)*candidate_fourier_data(ss,:,jj,mode); end end end end full_fourier_data(acs_line_loc,:,:) = acs_data; %Image reconstruction using IFFT2 and sum-of-squares if nargout > 1 coil_img = fftshift(fftshift(ifft2(ifftshift(ifftshift(full_fourier_data,1),2)),1),2); rec_img = sqrt(sum(abs(coil_img).^2,3)); end coef0=fit_coef; 12. References 1. B. HornRobot Vision, MIT Press, 1986, Chaps 6, 7. 2. Korin, Felmlee. Adaptive technique for three-dimensional MR imaging of moving structures. Radiology. 1990 Oct;177(1):217-21. 3. NessAiver. All you really need to know about MRI physics. 1997. 4. Hoa, D. (2014), Benefits of Parallel Imaging, June 8, 2014. Retrieved from: http://www.imaios.com/en/e-Courses/e-MRI/Parallel-imaging/benefits-parallel-techniques 5. Kastler. Comprendre l'IRM. 2006. 6. Hood, Ho. Chemical shift: the artifact and clinical tool revisited. Radiographics. 1999 Mar-Apr; 19(2):357-71. 7. Bydder, Rahal. The magic angle effect: a source of artifact, determinant of image contrast, and technique for imaging. J Magn Reson Imaging. 2007 Feb;25(2):290-300. 8. R. Gonzales, R. WoodsDigital Image Processing, Addison-Wesley Publishing Company, 1992, pp 81 - 125. 9. Certaines and Cathelineau. Safety aspects and quality assessment in MRI and MRS: a challenge for health care systems in Europe. J Magn Reson Imaging. 2001 Apr;13(4):632-8. 10. Ihalainen, Sipila. MRI quality control: six imagers studied using eleven unified image quality parameters. European radiology. 2004 Oct;14(10):1859-65 A. JainFundamentals of Digital Image Processing, Prentice-Hall, 1989, pp 15 - 20. 11. MarionAn Introduction to Image Processing, Chapman and Hall, 1991, Chap. 9. Read More