MRI Simulation
“I made it very clear that the physics of MRI is not an easy subject. NO ONE GETS IT THE FIRST TIME! In fact, I have
stated many times that you never really understand MRI, you just get used to it…. Don’t worry. Most people get it on the fifth or sixth time through.” — A Professor’s Lecture Notes (from the Web)
Introduction
 |
Block diagram of an MRI scanner. (from the Web) |
MRI is based on the quantum mechanical property of proton spin. By applying various magnetic fields, the protons within the brain, can be stimulated into generating an RF (radio frequency) signal which can be detected by receiver circuitry.
The MRI scanning technique “divides” the brain into many small “cubes” called voxels (volume elements). The applied magnetic fields are designed such that the protons within each voxel will generate an RF signal which is unique for each voxel. Using the unique signals, image reconstruction software can recreate a 3-D image of the brain.
Additionally, the applied magnetic fields are designed so that, depending on the type of substance within each voxel, the
protons within each voxel will “modulate” the RF signal, allowing the image reconstruction software to determine the type of substance within each voxel.
The Sections in this Overview (adapted from Liang page 8):
- Spin physics: Theory
This section describes the physics (mostly classical) within a voxel before magnetic fields are applied. - Bloch Equation
The Bloch equation is the classical (non-quantum mechanical) equation which describes how a voxel will respond to an applied magnetic field. - Signal Generation (The Pulse Sequence):
This section describes the magnetic fields generated by the MRI equipment.- Static Magnetic Field (B0): This is the main (high intensity, super-conducting) magnetic field. (by convention called the B0 field). This static field is always on and is not under user control.
- The Pulse Sequence: The pulse sequence is the user controlled process of dynamically generating magnetic fields, and is the heart of the MRI technique. The pulse sequence consists of two fields: The B1 field, and the Gradient fields.
- Gradient Fields (X,Y,Z): The Gradient Fields (GradientX, GradientY, GradientZ) are applied by adding a small “gradient” to the B0 field in such a way that each voxel has a unique value of B0. (In practice, usually only one gradient is applied at a time, therefore selecting either a plane, row, or column of voxels)
- RF Magnetic Field (B1): The RF magnetic field is applied to the voxels (by convention called the B1 field), and in response causes the Voxels to generate an RF signal, which can be detected by the MRI receiving circuitry.
- Signal Acquisition:
The RF signal generated by the voxels is acquired by MRI receiver circuitry. This signal is then “converted” to a form which can be utilized by the image reconstruction software. - Image Reconstruction:
The acquired signal data is processed by a two-dimensional Fast Fourier Transform (2-D FFT), which reconstructs an image of the scanned voxels. - High-resolution MRI (Limitations):
(under construction) Signal strength, signal to noise ratios…
Open Source MRI Simulator Source Code
- MRI Files mri.cpp, mri.h, brainmri128.h
Note: Code for “educational” purposes only. Some files may need to be hack’ed depending on directories structures and filenames.
Signal Processing Flow Diagram:
From a “signal processing” viewpoint, the MRI process can be reduced to three equations (NRCIOM pp. 209-210):
- Equation 3.1: Signal Generation
- Equation 2.1: The Bloch equation
- Equation 4.1: Signal Acquisition
____________________________________
| |
| MRI transmitter (section 3) |
| Equation 3.1 |
|____________________________________|
||B0 ||B1 ||Gradient
|| || ||
____/_____/______/_______________
| |
| Subject (voxels) (section 2) |
| Bloch equation, Equation 2.1 |
|____________________________________|
||signal = dM/dt
||
___________/_______________________
| |
| MRI receiver (section 4) |
| Equation 4.1 |
|....................................| "Laboratory" Frame
| Demodulator, LowPassFilter | _______________________
| Analog to Digital converter |
|____________________________________| "Rotating" Frame
||real ||imaginary
|| ||
___________/______/_______________
| |
| 2-D FFT |
|____________________________________|
||
/
Reconstructed Image
Rotating Frame vs. Laboratory Frame
In the literature, the laboratory frame is the stationary frame of reference.
The rotating frame, is a frame of reference rotating at the larmor Frequency (equation 3.2).
It is useful for conceptual simplicity. Liang pp. 72-76. The rotating frame is essentially
equivalent to viewing signals after they have gone through the demodulator.
MRI coordinate system:
The following are somewhat arbitrary, for visualizing the coordinate system:
- Z axis: in through neck, out through top of head
- X axis: in one ear, and out the other
- Y axis: in through back of head, and out between the eyes
1. Basic spin physics
This section is largely from Feynman_1964 II, chs. 34,35, and Liang pages 57-67, and is mostly a classical (non-quantum mechanical) description of spin physics.
[Note: boldface = vector; For vectors: · = dot product; x = cross product]
MRI physics is largely concerned with the physics of the hydrogen atom,
and in particular the hydrogen nuclei, i.e. a proton.
Protons possess a quantum mechanical property known as
spin angular momentum: J. The goal of this section is to start with
the spin angular momentum of an individual proton and finish with the bulkMagnetization and equilibriumMagnetization
of a voxel. The reason to do this is that
an individual proton’s spin is not observable by the MRI equipment, however
the bulkMagnetization and equilibriumMagnetization
are parameters of the Bloch equation, and the Bloch equation is “observable”
by the MRI equipment.
The components of J: Jx, Jy, Jz are related by the formula: (See Feynman_1964 II-34-10)
1.1) J · J = Jx^2 + Jy^2 + Jz^2
Each component (for example, Jz) is quantized to certain discrete values, by the formula:
1.2) Jz = j * hbar
where hbar = Planck’s constant divided by 2 pi = 1.055e-34 joules-second.
For hydrogen the values for j are: +1/2 and -1/2.
From this the equation:
1.3) J · J = j * (j + 1) * hbar^2
can be derived. See Feynman_1964 II-34-11 [note: This equation is common in the literature, but doesn't seem to be particularly useful in MRI calculations(?) However, equation 1.2 is useful.]
Another property of the proton, useful in MRI, is the magnetic moment. It is related to
spin angular momentum by the formula:
1.4) magneticMoment = g * (qe / (2 * mp)) * J;
where: magneticMoment is in units of joules/Tesla, g = nuclear g-factor (a dimensionless quantity, For proton:2 * 2.79; neutron:-2 * 1.93; electron:2 * 1);
qe = charge (electron, 1.6e-19 coulombs); mp = mass (proton, 1.67e-27 kg); See Feynman_1964 II-34-4
In MRI, the quantity called gyromagnetic ratio is common and is defined by the equation:
1.5) gyromagneticRatio = g * (qe / (2 * mp));
For hydrogen the gyroMagneticRatio is: 2.675e8 Radians/(seconds-Tesla), 42.57MHz/Tesla,
and therefore from equations 1.4 and 1.5:
1.6) magneticMoment = gyromagneticRatio * J; Liang page 59
The bulk magnetization vector is related to the magnetic moment by the formula: Liang p.64
1.7) for (i = 0, '''bulkMagnetization''' = 0; i < NUMBEROFPROTONS; i++) {
'''bulkMagnetization''' += '''magneticMoment[i]''';
} /* iForLoop */
where, bulkMagnetization is in units of amps/meter, NUMBEROFPROTONS = number of hydrogen atoms within 1 voxel (about 6.7e19 for a 1×1x1 mm voxel)
Initially, before any external magnetic fields are applied,
the magneticMoment of each proton is pointing in a random direction,
so of the net bulkMagnetization vector is 0. (Liang p.64)
At this point, the question could be asked, why isn’t the
net bulkMagnetization vector always 0? For example, suppose
an external field is applied, and the protons align themselves with the field.
Then by equation 1.2, it would seem that half of the protons should have j = +1/2 and
half of them should have j = -1/2, and so by equation 1.6 and 1.7 the net
bulkMagnetization should be 0 (and MRI wouldn’t work).
Actually, this is pretty much what happens,
but by a phenomenon known as Zeeman splitting (see equilibriumMagnetization below), there are slightly more protons
with +1/2 than -1/2, so there is a small net nonzero bulkMagnetization
(and this is what allows MRI to work)
The equilibriumMagnetization vector (also called the net magnetization vector)
is the value that the bulkMagnetization vector will obtain when an
external static magnetic field (such as B0, See section 3a, Static Magnetic Field below)
is applied, and sufficient time is allowed for equilibrium to be established.
The equilibriumMagnetization vector is determined from the phenomenon
known as Zeeman splitting, and the Boltzmann relationship as follows: (Liang pp. 64-67)
From quantum (classical?) theory, (Feynman_1964 II-15-2,II-34-11):
1.8) potentialEnergy = -magneticMoment · B0
Since by convention in MRI, B0 is only in the ‘z’ direction, this simplifies to:
1.9) potentialEnergy = -magneticMomentz * B0 = -gyromagneticRatio * Jz * B0 = -gyromagneticRatio * j * hbar * B0
From equation 1.2, the values for j are: +1/2 (spin up) and -1/2 (spin down).
1.10) potentialEnergySpinUp = -1/2 * gyromagneticRatio * hbar * B0
1.11) potentialEnergySpinDown = 1/2 * gyromagneticRatio * hbar * B0
1.12) dEnergy = potentialEnergySpinDown – potentialEnergySpinUp = gyromagneticRatio * hbar * B0
This nonzero result for dEnergy is the Zeeman splitting phenomenon.
From the Boltzmann relationship (Feynman I-40-9, quantum statistical mechanics):
1.13) spinRatio = numberSpinUp / numberSpinDown = exp(dEnergy / (boltzmannConstant * temperature));
Where, boltzmannConstant = 1.38e-23 joules/K, temperature = 310 Kelvin (body temperature)
Since after applying B0, the protons within a voxel can only be either spin up or spin down,
then: NUMBEROFPROTONS = numberSpinUp + numberSpinDown, and:
1.14) spinDensity = numberSpinUp – numberSpinDown = NUMBEROFPROTONS * (spinRatio – 1) / (spinRatio + 1);
And from equation 1.14, equation 1.7 reduces to:
1.14.1) bulkMagnetization = (numberSpinUp * |magneticMomentz| + numberSpinDown * -|magneticMomentz|) * k;
Where, k is the unit vector in the z direction
Which simplifies to:
1.15) equilibriumMagnetization = bulkMagnetization = spinDensity * |magneticMomentz| * k;
Example: Calculating the value of equilibriumMagnetization:
From equations 1.6 and 1.2:
|magneticMomentz| = gyromagneticRatio * 1/2 * hbar = 2.675e8 * 0.5 * 1.055e-34 = 1.41e-26 joules/tesla
Let B0 = 1.5 tesla (a typical MRI value), and from equation 1.12,
dEnergy = gyromagneticRatio * hbar * B0 = 2.675e8 * 1.055e-34 * 1.5 = 4.23e-26 joules
From equation 1.13: spinRatio = exp(dEnergy / (boltzmannConstant * temperature)) =
exp(4.23e-26 / (1.38e-23 * 310)) = exp(9.9e-6) = 1.00000990
For a 1mm3 voxel, NUMBEROFPROTONS = 6.7e19, and from equation 1.14:
spinDensity = NUMBEROFPROTONS * (spinRatio – 1) / (spinRatio + 1) = 6.7e19 * 4.95e-6 = 3.32e14 protons per mm3
And therefore from equation 1.15:
equilibriumMagnetization = 3.32e14 * 1.41e-26 * k = 4.68e-12 Joules/tesla-mm3 = 4.68e-3 amps/meter
The number for spinDensity (3.3e14) vs. NUMBEROFPROTONS (6.7e19), shows that only about
1 proton in 200,000 contributes to the equilibriumMagnetization. This is unfortunate
as it greatly reduces the signal strength available for developing high-resolution MRI.
2. The Bloch Equation
The Bloch equation (NRCIOM p. 209) describes the time rate of change of the
bulkMagnetization vector for a voxel. It is this time rate of change
which is detected by the MRI receiver circuitry.
From the MRI perspective, the Bloch equation describes
everything that is needed to know about the behavior of a voxel:
|
2.1) dM/dt = gyromagneticRatio * (M x B) – (Mx * i + My * j)/t2 – (Mz * k – MzEq * k)/t1 |
Where:
- i, j, k are the unit vectors in the x, y, z directions respectively.
- M is the bulkMagnetization vector.
- Mx, My, Mz are the individual components of the bulkMagnetization vector.
- MzEq is the ‘z’ component of the equilibriumMagnetization vector.
- B is the magneticField vector (in units of Tesla). See equation 3.1 below.
- t1, t2: relaxation times, empirically derived, generally unique for each substance type. See table 3.3
In general, there is no closed-form solution to the Bloch equation, NRCIOM page 38.
The Bloch equation broken into individual components:
- 2.2) dMx/dt = gyromagneticRatio * (My * Bz – Mz * By) – Mx/t2;
- 2.3) dMy/dt = gyromagneticRatio * (Mz * Bx – Mx * Bz) – My/t2;
- 2.4) dMz/dt = gyromagneticRatio * (Mx * By – My * Bx) – (Mz – MzEq)/t1;
3. Signal Generation (The Pulse Sequence):
The MRI scanning equipment generates a magnetic field described by the equation (NRCIOM page 210):
| 3.1) B(r,t) = ((B0 + G(r,t)) * k) + (B1(t) * cos(localFrequency * t) * i) + (B1(t) * sin(localFrequency * t) * j); |
Where,
- B(r,t) is the magneticField vector generated by the MRI equipment(in units of Tesla)
- B0 is the static magnetic field
- G(r,t) is the Gradient field
- B1(t) is the RF field
- i, j, k are the unit vectors in the x, y, z directions respectively.
- localFrequency: set to the larmorFrequency (see equation 3.2) of the desired slice (figure 3.2, table 3.2)
- r: position (in 3-D space)
- t: time
3.2) larmorFrequency = gyromagneticRatio * B0. (example: if B0 = 1.5 tesla, the larmorFrequency = 63.855 MHz);
3a. Static Magnetic Field (B0)
After application of the B0 vector in the ‘z’ direction, the Bloch equation
yields (NRCIOM page 210, also see Goldman page7):
3.3) Mx = exp(-t/t2) * (Mx0 * cos(larmorFrequency * t) + My0 * sin(larmorFrequency * t))
3.4) My = exp(-t/t2) * (-Mx0 * sin(larmorFrequency * t) + My0 * cos(larmorFrequency * t))
3.5) Mz = exp(-t/t1) * Mz0 + MzEq * (1 – exp(-t/t1))
Where, Mx0, My0, Mz0 are the initial values of the M vector
Free Induction Decay (FID):
Equations 3.3 and 3.4 are exponentially decaying sinewaves oscillating at the larmorFrequency.
In the MRI literature, this is known as a free induction decay signal. (Elster page 46)
It is this free induction decay that is detected by the MRI receiver circuitry (see equation 4.1)
Note: initially at the beginning of time, Mx0, My0, Mz0 are 0. After application of the
excitation process (See “Spin Echo Step 1″ below), these values will no longer be 0,
allowing the free induction decay to occur.
From equations 3.3-3.5 it can be seen that at equilibrium (time = large),
then Mx = My = 0 and Mz = MzEq
3b. The Pulse Sequence
The pulse sequence is the user controlled process of dynamically generating magnetic
fields, and is the heart of the MRI technique.
The pulse sequence consists of two fields: The B1 (RF) field, and the Gradient fields.
Various pulse sequences can be designed, corresponding to the various
imaging techniques available on the MRI hardware.
Imaging Techniques (adapted from Morris):
- Single Point
- Line
- Planar (2-D)
- Projection Reconstruction
- Fourier method
- Spin Echo Imaging
- Gradient Echo Imaging
- Echo Planar Imaging
- Spiral Imaging
- And others…
- Volume (3-D)
The single point and line imaging techniques are designed to scan a single voxel
or a line of voxels respectively at a time. These are older(?) techniques which are
generally slower than the planar or volume techniques. Currently, the most popular
technique is the Fourier imaging method, of which there are many different variants.
This document will discuss the basic Spin Echo imaging technique.
Figure 3.1a: Pseudo-Code for the Spin Echo Process
for (loopZ = 0; loopZ < NUMBEROFVOXELSZ; loopZ++) {
for (loopY = 0; loopY < NUMBEROFVOXELSY; loopY++) { _________
SelectSlice(loopZ) ; /
SelectPhase(loopY); ||
for (loopX = 0; loopX < NUMBEROFVOXELSX; loopX++) { Figure 3.1b
SelectFrequencyAndReadout(loopX); ||
} /* xForLoop */ ____/___
} /* yForLoop */
} /* zForLoop */
Figure 3.1b: Spin Echo Pulse Sequence Timing Diagram (note7):
90deg | 180deg 90degree
| ||| |
RF ---------|||------------------|||---------------------- ... ---------|||---
pulse | ||| |
(B1) | sliceN sliceN+1
->| |<Tgz/2
Tgz->| |< Tgz->| |<
____ ____ __+gz
Gz _______| | _____________| |_____ _____ ... ____
note1 |_| note5 __-gz
note2
->| |< Tgz/2
->| |< Tgzgy
___
Gy ________________| |_________________________________ ... ____
note3 |___|
(Tgy) note6
->| |< Tgx/2 ->| |< Tgx
___ note4 ______
Gx ________________| |___________________| |______ ... ____
|<----- Te/2 ------>|
|<------------- Te ------------->|
|<------------------------- Tr -------------------------->|
| fid1 ||| spinEcho
NMR || |||||
Signal _________|||__________________________|||||||______ ... ____
Notes:
- Gz/RF: Tgz time = slice selection (z direction) (and 90 degree RF pulse), Vlaardingerbroek pp. 59-62
- Gz: Post-excitation re-phasing time, Liang page 149
- Gy: Tgy time = “phase selection” (y direction), Gy amplitude is a variable, Jin page 34
- Gx: Tgx/2: De-phasing time, Vlaardingerbroek page 66, Elster page 81
- Gz/RF: 180 degree RF pulse, Vlaardingerbroek pp. 65-67, Jin pages 29-30,34
- Gx: Tgx: Readout time and “frequency selection” (x direction), Vlaardingerbroek page 66, 67
- adapted from Vlaardingerbroek page 67, and Jin page 34
Table 3.1: Spin Echo Pulse Sequence Timing Numbers (typical)
| Symbol | Parameter | Value | Notes |
|---|---|---|---|
| Tgz | Slice Duration | 2.3-4.7 milliseconds | Vlaardingerbroek page 61 |
| Tgzgy | Slice to Phase time | 0? | |
| Tgy | Phase Duration | Tgx / 2 | Set to Tgx/2 for simplicity |
| Tgx | Readout Time | 5 milliseconds | |
| Te | Echo Time | 10-100 milliseconds | Set to maximize tissue contrast, Elster page 100, Vlaardingerbroek page 58 |
| Tr | Repetition Time | ~1 seconds | Vlaardingerbroek page 58 |
note: timing values currently ignore rise/fall/skew times.
Spin Echo Step 1: Slice Selection and the RF Pulse (from figure 3.1b, note 1)
In order for the MRI receiver circuitry to detect a signal (see section 4),
it is necessary for the voxels to generate a RF signal.
From equations 3.3-3.5, at equilibrium (time = large), the voxels are static (Mx = My = 0 and Mz = MzEq).
The voxels can be perturbed from this equilibrium state by applying a RF
magnetic field (B1) in the plane perpendicular to B0 (i.e. in the X-Y plane, a.k.a. transverse plane).
This process is known in the literature as excitation.
Elster pages 25-27, Liang pages 70-72,77
From equation 3.1, the MRI transmitter “excitation” RF magnetic field (B1) is:
3.6a) B1x = B1(t) * cos(localFrequency * t);
3.6b) B1y = B1(t) * sin(localFrequency * t);
And note that the localFrequency is set to oscillate at the larmorFrequency,
which from equation 3.2 is: B0 * gyroMagneticRatio.
The larmorFrequency is the “resonance” frequency of a voxel, and the excitation effect rapidly drops off
as frequencies move away from the larmorFrequency. This is the “resonance” part of Magnetic Resonance Imaging,
and is key to the MRI process, as it allows the user to isolate a voxel(s) by adding a
small gradient to B0 and thereby giving each voxel a unique larmorFrequency. Elster pages 25-26
After the voxels have been excited, and the RF magnetic field turned off,
the voxels will generate a Free Induction Decay (FID), oscillating at the
larmorFrequency as per equations 3.3, 3.4. This RF signal can then be detected
by the MRI receiver circuitry.
Depending on the imaging technique involved, various numbers of voxels
will be excited at a time. For example, with “single point imaging”,
only a single voxel will be excited at a time. With planar imaging techniques
(such as Spin Echo), an entire slice (plane) will be excited at a time.
By convention, the slice selection process selects a plane perpendicular to the Z
direction, i.e. in the x-y (transverse) direction. Elster page 75
The slice selection first begins by applying a small gradient Gz (in the Z direction) to the B0 magnetic field.
From equation 3.1: (B0 + G(r,t)) * k.
Now each x-y plane of voxels will have a unique larmorFrequency:
Table 3.2: Slice Selection, Gradient Gz vs. larmorFrequency
| Slice | Gradient, Gz (Tesla) | larmorFrequency (Hz) | notes | ||
|---|---|---|---|---|---|
| minimum | middle | maximum | |||
| -10 | -0.001 | -44702.6 | -42573.9 | -40445.2 | bottom of brain (neck) |
| … | … | … | … | … | … |
| -2 | -0.0002 | -10643.5 | -8514.79 | -6386.09 | |
| -1 | -0.0001 | -6386.09 | -4257.39 | -2128.7 | |
| 0 | 0 | -2128.7 | 0 | 2128.7 | middle of brain |
| +1 | 0.0001 | 2128.7 | 4257.39 | 6386.09 | |
| +2 | 0.0002 | 6386.09 | 8514.79 | 10643.5 | |
| … | … | … | … | … | … |
| +9 | 0.0009 | 36187.9 | 38316.6 | 40445.2 | top of brain |
Example assumes: Gz = 10.0 mT/meter, SLICETHICKNESS = 10mm
Note: frequency numbers are in the rotating frame (i.e. after the demodulator). For the laboratory frame add to the numbers: B0 * gyroMagneticRatio = 63.38609MHz @ B0 = 1.5Tesla
And from table 3.2, the localFrequency can then be set to 63.38609MHz + 0Hz to select slice 0, or set to 63.38609MHz + 8514.79Hz to select slice 2, etc.
To complete the slice selection process the function B1(t) from
equations 3.6a, 3.6b, needs to be determined.
A common function is sinc(t), as this function will generate
a rectangular band of frequencies in the range of the selected slice. Elster pp. 74-75, Vlaardingerbroek pp. 61-62
Other functions used are: rectangular pulse and SLR pulses. Elster page 75
In addition to the type of function for B1(t), the magnitude of B1 and the value
of Tgz (SLICEDURATION of the RF pulse) need to be determined.
The SLICEDURATION value can be found by: (from Vlaardingerbroek pp. 60-62)
An example with: Gz = 10.0 mT/meter, slice width = 10mm, sinc(t) with 2 side lobes on each side.
GRADIENTZVALUE = (10.0e-6); /* Tesla/mm */
SLICETHICKNESS = (10); /* mm */
EFFECTIVELENGTH = (2.0 * PI / ( GYROMAGNETICVALUE * SLICETHICKNESS * GRADIENTZVALUE)) = 0.000234885 seconds;
NUMBEROFSIDELOBES = (2 * 2); /* 2 lobes on each side */
NUMBEROFLOBES = (NUMBEROFSIDELOBES + 2.0);
SLICEDURATION = (NUMBEROFLOBES * EFFECTIVELENGTH) = 0.00140931 seconds;
The magnitude of B1 is determined by the desired “flip angle” (Elster page25, Liang pages79-80).
The flip angle is defined by: angle = gyroMagneticRatio * B1 * EFFECTIVELENGTH;
(Note: EFFECTIVELENGTH is equal to the length of a rectangular pulse that
would generate an equivalent flip angle to a sinc pulse applied for a length of SLICEDURATION. Vlaardingerbroek page 61)
The flip angle determines the quantity of Mz transferred to the transverse plane (Mx, My).
For a 90 degree pulse, Mz moves from its equilibrium value to zero,
and Mx, My move from 0 to their maximum values.
For a 90 degree RF flip angle, the magnitude of B1 can be found by: (Vlaardingerbroek pp. 60-62, Elster page25)
B1MAGNITUDE = (PI / (2.0 * GYROMAGNETICVALUE * EFFECTIVELENGTH)) = 2.5e-005 Tesla
(note: B1MAGNITUDE * GYROMAGNETICVALUE = 1064.35 Hz);
Figure 3.2: Pseudo-Code for Slice Selection/RF excitation:
for (loopZ = 0; loopZ < NUMBEROFVOXELSZ; loopZ++) { // # of slices
for (loopY = 0; loopY < NUMBEROFVOXELSY; loopY++) { // # of phase encoding steps
...
for (t = 0 ; t < SLICEDURATION; t += dt) { /* Slice Selection */
localFrequency = function(loopZ); /* see table 3.2 */
localOscillatorX = cos(localFrequency * t) ;
localOscillatorY = sin(localFrequency * t) ;
sincTime = (pi * NUMBEROFLOBES / SLICEDURATION) * (t - SLICEDURATION / 2.0);
sincPulse = sin(sincTime) / sincTime;
/* quadrature coil (from Jin pages 11-16, 25) */
B1x = B1Magnitude * sincPulse * localOscillatorX;
B1y = -B1Magnitude * sincPulse * localOscillatorY;
for (k = 0; k < NUMBEROFVOXELSZ; k++) {
for (j = 0; j < NUMBEROFVOXELSY; j++) {
for (i = 0; i < NUMBEROFVOXELSX; i++) {
voxel[k][j][i].BlochEquation(B1x, B1y, B0 + Gz(k - NUMBEROFVOXELSZ/2));
} /* iForLoop */
} /* jForLoop */
} /* kForLoop */
} /* tforLoop */
...
SelectPhase(loopY);
SelectFrequencyAndReadout();
} /* yforLoop */
} /* zForLoop */
Below are results of a simulation (in the laboratory frame of reference) of the slice selection process, including the sinc waveform,
B1x, Mx, Mz. Note that Mz moves from its equilibrium value to approximately zero,
and Mx moves from 0 to its maximum value (consistent with a “90 degree” pulse).
Figure 3.3: Slice Selection and RF pulse timing
Spin Echo Step 1a: Post-excitation Re-phasing Time (from figure 3.1b, note 2)
This process is somewhat difficult to visualize.
During the slice selection process different regions within the slice
may be oscillating at different frequencies. For example, from table 3.2,
slice zero may vary from -2128.7Hz to + 2128.7Hz. Over time this effect
will cause the regions to lose “phase coherence” with one another.
If not corrected this can lead to an undesirable signal loss.
Example:
from the trigonometric formula: sin(a) + sin(b) = 2 * sin(a/2 + b/2) * cos(a/2 – b/2)
Then let, a = frequency1 * t, b = frequency2 * t, and frequency2 = frequency1 + dfrequency
so that, signal = 2 * sin((f1 + f1 + df) * t/2) * cos((f1 – f1 – df) * t/2)
= 2 * sin(f1 * t + df * t/2) * cos(-df * t/2)
Therefore, whenever the quantity “df * t” equals pi, the signal disappears.
To correct this problem, assume that halfway through the 90 degree pulse (t = ~0.725msec = tgz/2, figure 3.3),
the regions were in perfect phase coherence (df * t equals 2pi), then at the end of the 90 degree pulse
(t = ~1.4msec = tgz, figure 3.3) the regions will have lost phase coherence.
At this point applying a negative gradient
for a time equal to tgz/2 will cause the regions to regain phase coherence.
The analogy is runners in a race. After beginning from the start line they run at
different speeds. If after running for say five minutes, they all turn around
and run for another five minutes, then they should all arrive back at the start line at the same moment.
(under construction) add: wave packet spreading analogy
Spin Echo Step 2: Phase Selection (and de-phasing X) (from figure 3.1b, note 3,4)
At this point a slice has been selected. Now in order to reconstruct
an image, it is necessary to be able to uniquely distinguish between
the voxels within the x-y plane.
By convention, the X direction is distinguished by frequency encoding,
during the Readout time. This is accomplished by applying a gradient Gx (in the X direction)
to B0, which then gives each voxel in the X direction a unique larmorFrequency. (see: spin echo step 4)
This still leaves the Y direction to be determined.
Conceptually, the Y direction is encoded by giving each row of voxels a unique
“frequency” in a very low frequency range (frequency < 1 Hz). Rajan page 23.
Figure 3.4, below, shows the spatial frequency distribution for a 16×16 voxel X-Y slice.
| X Y |
0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 |
| 15 | 0.437 -53217 |
0.437 -46565 |
0.437 -39913 |
0.437 -33261 |
0.437 -26609 |
0.437 -19957 |
0.437 -13304 |
0.437 -6652 |
0.437 0 |
0.437 6652 |
0.437 13304 |
0.437 19957 |
0.437 26609 |
0.437 33261 |
0.437 39913 |
0.437 46565 |
| 14 | 0.375 -53217 |
0.375 -46565 |
0.375 -39913 |
0.375 -33261 |
0.375 -26609 |
0.375 -19957 |
0.375 -13304 |
0.375 -6652 |
0.375 0 |
0.375 6652 |
0.375 13304 |
0.375 19957 |
0.375 26609 |
0.375 33261 |
0.375 39913 |
0.375 46565 |
| 13 | 0.312 -53217 |
0.312 -46565 |
0.312 -39913 |
0.312 -33261 |
0.312 -26609 |
0.312 -19957 |
0.312 -13304 |
0.312 -6652 |
0.312 0 |
0.312 6652 |
0.312 13304 |
0.312 19957 |
0.312 26609 |
0.312 33261 |
0.312 39913 |
0.312 46565 |
| 12 | 0.250 -53217 |
0.250 -46565 |
0.250 -39913 |
0.250 -33261 |
0.250 -26609 |
0.250 -19957 |
0.250 -13304 |
0.250 -6652 |
0.250 0 |
0.250 6652 |
0.250 13304 |
0.250 19957 |
0.250 26609 |
0.250 33261 |
0.250 39913 |
0.250 46565 |
| 11 | 0.187 -53217 |
0.187 -46565 |
0.187 -39913 |
0.187 -33261 |
0.187 -26609 |
0.187 -19957 |
0.187 -13304 |
0.187 -6652 |
0.187 0 |
0.187 6652 |
0.187 13304 |
0.187 19957 |
0.187 26609 |
0.187 33261 |
0.187 39913 |
0.187 46565 |
| 10 | 0.125 -53217 |
0.125 -46565 |
0.125 -39913 |
0.125 -33261 |
0.125 -26609 |
0.125 -19957 |
0.125 -13304 |
0.125 -6652 |
0.125 0 |
0.125 6652 |
0.125 13304 |
0.125 19957 |
0.125 26609 |
0.125 33261 |
0.125 39913 |
0.125 46565 |
| 9 | 0.062 -53217 |
0.062 -46565 |
0.062 -39913 |
0.062 -33261 |
0.062 -26609 |
0.062 -19957 |
0.062 -13304 |
0.062 -6652 |
0.062 0 |
0.062 6652 |
0.062 13304 |
0.062 19957 |
0.062 26609 |
0.062 33261 |
0.062 39913 |
0.062 46565 |
| 8 | 0.000 -53217 |
0.000 -46565 |
0.000 -39913 |
0.000 -33261 |
0.000 -26609 |
0.000 -19957 |
0.000 -13304 |
0.000 -6652 |
0.000 0 |
0.000 6652 |
0.000 13304 |
0.000 19957 |
0.000 26609 |
0.000 33261 |
0.000 39913 |
0.000 46565 |
| 7 | -0.062 -53217 |
-0.062 -46565 |
-0.062 -39913 |
-0.062 -33261 |
-0.062 -26609 |
-0.062 -19957 |
-0.062 -13304 |
-0.062 -6652 |
-0.062 0 |
-0.062 6652 |
-0.062 13304 |
-0.062 19957 |
-0.062 26609 |
-0.062 33261 |
-0.062 39913 |
-0.062 46565 |
| 6 | -0.125 -53217 |
-0.125 -46565 |
-0.125 -39913 |
-0.125 -33261 |
-0.125 -26609 |
-0.125 -19957 |
-0.125 -13304 |
-0.125 -6652 |
-0.125 0 |
-0.125 6652 |
-0.125 13304 |
-0.125 19957 |
-0.125 26609 |
-0.125 33261 |
-0.125 39913 |
-0.125 46565 |
| 5 | -0.187 -53217 |
-0.187 -46565 |
-0.187 -39913 |
-0.187 -33261 |
-0.187 -26609 |
-0.187 -19957 |
-0.187 -13304 |
-0.187 -6652 |
-0.187 0 |
-0.187 6652 |
-0.187 13304 |
-0.187 19957 |
-0.187 26609 |
-0.187 33261 |
-0.187 39913 |
-0.187 46565 |
| 4 | -0.250 -53217 |
-0.250 -46565 |
-0.250 -39913 |
-0.250 -33261 |
-0.250 -26609 |
-0.250 -19957 |
-0.250 -13304 |
-0.250 -6652 |
-0.250 0 |
-0.250 6652 |
-0.250 13304 |
-0.250 19957 |
-0.250 26609 |
-0.250 33261 |
-0.250 39913 |
-0.250 46565 |
| 3 | -0.312 -53217 |
-0.312 -46565 |
-0.312 -39913 |
-0.312 -33261 |
-0.312 -26609 |
-0.312 -19957 |
-0.312 -13304 |
-0.312 -6652 |
-0.312 0 |
-0.312 6652 |
-0.312 13304 |
-0.312 19957 |
-0.312 26609 |
-0.312 33261 |
-0.312 39913 |
-0.312 46565 |
| 2 | -0.375 -53217 |
-0.375 -46565 |
-0.375 -39913 |
-0.375 -33261 |
-0.375 -26609 |
-0.375 -19957 |
-0.375 -13304 |
-0.375 -6652 |
-0.375 0 |
-0.375 6652 |
-0.375 13304 |
-0.375 19957 |
-0.375 26609 |
-0.375 33261 |
-0.375 39913 |
-0.375 46565 |
| 1 | -0.437 -53217 |
-0.437 -46565 |
-0.437 -39913 |
-0.437 -33261 |
-0.437 -26609 |
-0.437 -19957 |
-0.437 -13304 |
-0.437 -6652 |
-0.437 0 |
-0.437 6652 |
-0.437 13304 |
-0.437 19957 |
-0.437 26609 |
-0.437 33261 |
-0.437 39913 |
-0.437 46565 |
| 0 | -0.500 -53217 |
-0.500 -46565 |
-0.500 -39913 |
-0.500 -33261 |
-0.500 -26609 |
-0.500 -19957 |
-0.500 -13304 |
-0.500 -6652 |
-0.500 0 |
-0.500 6652 |
-0.500 13304 |
-0.500 19957 |
-0.500 26609 |
-0.500 33261 |
-0.500 39913 |
-0.500 46565 |
For each voxel: the top number is frequencyY, the bottom number is frequencyX.
for frequencyX: Hz/Voxel: bandwidthX / NUMBEROFVOXELSX: 6652.18, see example in step 4 below
for frequencyY: Hz/Voxel: bandwidthY / NUMBEROFVOXELSY: 0.0625, (where bandwidthY = 1/Tr = 1) see below
Note: frequency numbers are in the rotating frame (i.e. after the demodulator). For the laboratory frame add to the numbers: B0 * gyroMagneticRatio = 63.38609MHz @ B0 = 1.5Tesla
The actual implementation of generating the Y frequencies using the phase encoding
process is somewhat difficult to visualize. (Elster pp. 82-89, Rajan pp. 21-24, Horowitz chapter 11)
The idea is to intentionally introduce a phase shift to the Y voxels.
Hence, the term phase encoding.
(This phase encoding process is somewhat similar to the process of re-phasing the Z spins in step 1a).
At each phase encoding step (See loopY forLoop below), the magnitude of the gradient Gy is changed,
starting from say a maximum negative gradient (at loopY=0) decreasing to 0, then increasing to finally
a maximum positive gradient (at loopY=NUMBEROFVOXELSY).
Over time as loopY changes from 0 to NUMBEROFVOXELSY, a greater or lesser phase shift is induced,
and this has the effect of “manually” generating a frequency for the voxels in the Y direction.
Additionally, within each phase encoding step, the gradient Gy varies depending on
Y position, which has the effect of giving each row of voxels a unique frequency. The complete equation for Gy is:
Gy(loopY, j – NUMBEROFVOXELSY/2) =
(-GRADIENTPERVOXELY + loopY * dGRADIENTPERVOXELY) * (j – NUMBEROFVOXELSY/2);
An example with: FOVy = 250 mm, NumberOfVoxelsY = 16, READOUTTIME = 0.000150327 (seconds)
PHASEDURATION = (READOUTTIME / 2.0)
= 7.51633e-005 (seconds); /* seconds, see step4 for READOUTTIME value */
GRADIENTYMAXIMUM = (PI * NUMBEROFVOXELSY / (FIELDOFVIEWY * PHASEDURATION * GYROMAGNETICVALUE)) = 1e-005 (Tesla/mm);
GRADIENTPERVOXELY = (GRADIENTYMAXIMUM * FIELDOFVIEWY / NUMBEROFVOXELSY)
= 0.00015625 (Tesla/Voxel, maximum);
dGRADIENTPERVOXELY = (2.0 * GRADIENTPERVOXELY / NUMBEROFDISPLAYVOXELSY)
= 1.95312e-005 (Tesla/Voxel);
The value for GRADIENTYMAXIMUM is derived from the relationship from Elster page 83:
GRADIENTYMAXIMUM * PHASEDURATION * GYROMAGNETICVALUE/(2 * PI)
= 0.5 * NUMBEROFVOXELSY / FIELDOFVIEWY
This relationship is designed to maximize spatial resolution, so that at maximum gradient
strength there is a 180 degree phase shift
between voxels. This satisfies the Nyquist sampling theorem (Stearns chapter 2) requirement that: maximum voxel frequency = 0.5 * samplingFrequency, (where, in this case, samplingFrequency = bandwidthY = 1 / Tr, Tr = repetition time, the time between phase encoding steps, see figure 3.1b)
Figure 3.5: Pseudo-Code for “Implementation” of the phase encoding process
for (loopZ = 0; loopZ < NUMBEROFVOXELSZ; loopZ++) { // # of slices
for (loopY = 0; loopY < NUMBEROFVOXELSY; loopY++) { // # of phase encoding steps
SelectSlice(loopZ) ;
...
for (t = 0 ; t < PHASEDURATION; t += dt) { // phase encoding process
for (k = 0; k < NUMBEROFVOXELSZ; k++) {
for (j = 0; j < NUMBEROFVOXELSY; j++) {
for (i = 0; i < NUMBEROFVOXELSX; i++) {
voxel[k][j][i].BlochEquation(B0 + Gy(loopY, j - NUMBEROFVOXELSY/2));
} /* iForLoop */
} /* jForLoop */
} /* kForLoop */
} /* forLoop */
...
SelectFrequencyAndReadout();
} /* yforLoop */
} /* zForLoop */
De-phasing X
While the phase encoding process is in progress, an additional step is to purposely
de-phase X. This is done by applying the same gradient Gx that will be used for
readout (see step 4), but for half the time, see figure 3.1b note 4. (This is a
process similar to the re-phasing process of Z in step 1a.)
The reason to do this de-phasing is somewhat difficult to visualize.
If this step was left out, then at the beginning of the readout period
the X spins would immediately begin to de-phase (similar to what happens to Z in step 1a).
If however, this initial de-phasing step is done, then after the 180 degree pulse in step three (which has the effect of turning the “runners” around, see step 1a),
the X spins will begin to re-phase at the start of the readout period, and reach perfect phase coherence
at the center of the readout period (at which point they will again start to de-phase).
Spin Echo Step 3: 180 degree pulse and Contrast Resolution (from figure 3.1b, note 5)
At this point since the slice has been selected and the phase encoding process finished,
why not immediately began the Readout phase?
There are two problems: the first has to do with contrast resolution, and the second
concerns the rapid FID decay due to a quantity known as T2* (T2star).
Contrast Resolution
From the Bloch equation (equation 2.1), there are 3 quantities which are available
to distinguish one voxel from another: t1, t2, and MzEq (MzEq == spin density).
For most of the tissues of interest (For example: gray matter, white matter),
the spin densities are very similar and not useful for contrast resolution.
(see table 3.3)
Table 3.3: Values for Tissue Specific Parameters
| Tissue Type | T1 (milliseconds) (at 1.5T) |
T2 (milliseconds) (at 1.5T) |
Chemical Shift (note 1) | Spin Density (relative to water) |
|---|---|---|---|---|
| Water | 4000 | 2300 | 1.0 | 1.0 |
| CSF | 3000? | 200 | 1.0 | 1.0 |
| Fat | 260 | 60 | 1.0 + 3.0e-6 (220Hz @ 1.5T) |
0.9 |
| Grey Matter | 1000 | 106 | - | 0.84 |
| White Matter | 650 | 69 | - | 0.72 |
| Muscle | - | 40 | - | - |
note1: gyromagneticRatio(of tissue) = ChemicalShift * gyromagneticRatio(of water)
T1 relaxation: spin-lattice relaxation, longitudinal relaxation
T2 relaxation: spin-spin relaxation, transverse relaxation
However, from table 3.3, for gray matter and white matter, the difference between values for T2 should provide
sufficient contrast.
From Elster page 100, Liang page 223:
3.7a) signalIntensity = spinDensity * (1 – exp(-Tr/T1)) * exp(-Te/T2)
For T2-weighting this reduces to:
3.7b) signalIntensity = exp(-Te/T2)
From the Web, for maximum contrast between 2 tissues: maximize(signalIntensity1 – signalIntensity2);
Solving for Te (ECHOTIME):
ECHOTIME = (T2GRAYMATTER * T2WHITEMATTER * log(T2GRAYMATTER / T2WHITEMATTER) / (T2GRAYMATTER – T2WHITEMATTER)) = 85 milliseconds;
Therefore, in order to obtain sufficient contrast it is necessary to wait for a time equal to ECHOTIME = 85 milliseconds.
Unfortunately, due to T2*, it is not possible to simply wait.
T2* (T2star)
Following a 90 degree pulse, the FID signal should decay with a time constant
T2 (equations 3.3, 3.4). However, due to inhomogeneities in the B0 field, spins at different locations
experience slightly different magnetic fields, causing various regions to
have slightly different larmor frequencies, consequently causing phase
de-coherence, similar to Z in step 1a. Therefore, instead of decaying at T2, the FID decays at T2star.
From Elster page 36, Liang page 111:
3.8) 1/T2star = 1/T2 + gyroMagneticRatio * dB0.
Where, dB0 = variation in B0 (in parts per million, ppm)
Example: dB0 = 10 ppm, T2 = 69 milliseconds, then T2star = 0.37 milliseconds.
Since, ECHOTIME is much larger than T2star, the FID after the 90 degree pulse will have died out before
the Readout period is reached.
Therefore, a solution is to insert a 180 degree pulse halfway to ECHOTIME,
this will cause the spins to “turn around” (see step 1a) and by the middle of the Readout period
will have regained phase coherence.
Note: creating a 180 degree pulse is identical to a 90 degree pulse except: double the value of B1MAGNITUDE in step 1,
everything else in step one remains the same.
(under construction:) add CPMG, Goldman page 11
Spin Echo Step 4: Readout, Frequency Selection (Spin Echo) (from figure 3.1b, note 6)
Encoding the X direction is much more straightforward than the Y direction.
By convention, the X direction is encoded by frequency encoding,
during the Readout time. This is accomplished by applying a gradient Gx (in the X direction)
to B0, which then gives each voxel in the X direction a unique larmorFrequency.
See Figure 3.6 Pseudo-Code for Signal Acquisition (and Readout/FrequencySelection)
From Elster page 81:
GRADIENTXVALUE * GYROMAGNETICVALUE / (2.0 * PI ) = BANDWIDTH / FIELDOFVIEWX;
An example with: Gx = 10mT/meter, FOVx = 250 mm, NumberOfVoxelsX = 16;
FIELDOFVIEWX = (250.0); /* mm */
NUMBEROFVOXELSX = (16);
GRADIENTXVALUE = (10.0e-6); /* Tesla/mm (10mT/m) */
GRADIENTPERVOXELX = (GRADIENTXVALUE * FIELDOFVIEWX / NUMBEROFVOXELSX) = 0.00015625 (Tesla/Voxel);
BANDWIDTH = (GYROMAGNETICVALUE * GRADIENTXVALUE * FIELDOFVIEWX / (2.0 * PI)) = 106435 (Hz);
SAMPLEPERIOD = (1.0 / BANDWIDTH) = 9.39542e-006 (seconds), sampleFrequency:106435 (Hz);
READOUTTIME = (SAMPLEPERIOD * NUMBEROFVOXELSX) = 0.000150327 (seconds);
function: Gx() = (i – NUMBEROFVOXELSX/2) * GRADIENTPERVOXELX; (figure 3.6)
Hz/Voxel: BANDWIDTH / NUMBEROFVOXELSX = 6652.18; (see figure 3.4)
4. Signal Acquisition
The MRI receiver hardware detects a signal described by the equation (NRCIOM page 210):
|
4.1) S(t) = SUM( -dM(r,t)/dt · C(r) dr ); (sum over all voxels) |
Where,
- S(t) is the “raw” NMR signal from the receiver coil, Liang page 95, from Faraday’s law of induction, Feynman II 17-2, Jin pp. 139-143
- M(r,t) is the bulkMagnetization vector of the individual voxels
- C(r) is the receiver coil sensitivity, which by the principle of reciprocity is equal to: B1(r) / I, where,
- B1(r) is the magneticField vector generated by the MRI RF transmitter coil (transmitter coil assumed to be the same as the receiver coil), see equation 3.1
- Note: As B1 is usually homogeneous across all voxels, C(r) can be moved outside the integral
- I = unit current in coil which produces magnetic field B1
- rfCoilSensitivity estimate = 0.000234885 (Tesla/amp), C=B1/I, assumes: B1 = 0.23mT, I = 1amp (from: transmitter output = 1kw, (1kv, 1kohm))
- r: position (in 3-D space)
- t: time
Figure 3.6 Pseudo-Code for Signal Acquisition (and Readout/FrequencySelection)
for (loopZ = 0; loopZ < NUMBEROFVOXELSZ; loopZ++) { // # of slices
for (loopY = 0; loopY < NUMBEROFVOXELSY; loopY++) { // # of phase encoding steps
localFrequency = SelectSlice(loopZ) ;
SelectPhase(loopY);
...
sampleTime = SAMPLEPERIOD;
for (loopX = 0, t = 0 ; t < READOUTTIME; t += dTime) { // Readout
magnetX = 0.0;
magnetY = 0.0;
localOscillatorX = cos(localFrequency) ;
localOscillatorY = sin(localFrequency) ;
for (k = 0; k < NUMBEROFVOXELSZ; k++) {
for (j = 0; j < NUMBEROFVOXELSY; j++) {
for (i = 0; i < NUMBEROFVOXELSX; i++) {
voxel[k][j][i].BlochEquation(B0 + Gx(i - NUMBEROFVOXELSX/2),magnetX,magnetY);
} /* iForLoop */
} /* jForLoop */
} /* kForLoop */
/* ****** Signal Acquisition section ********* */
/* note: quadrature coil */
/* signal = -dM/dt * coilSensitivity (equation 4.1) */
rawSignalX = rfCoilSensitivity * (magnetX - lastmagnetX) / dTime ;
rawSignalY = rfCoilSensitivity * (magnetY - lastmagnetY) / dTime ;
lastmagnetX = magnetX;
lastmagnetY = magnetY;
/* demodulator (Liang page 98) */
demodulatedSignalX = rawSignalX * localOscillatorX;
demodulatedSignalY = -rawSignalY * localOscillatorY;
/* LowPassFilter (see Stearns chapter 8 ) */
filteredSignalX = LowPassFilter(demodulatedSignalX);
filteredSignalY = LowPassFilter(demodulatedSignalY);
/* analog to digital converter */
/* Stearns chapter 2, Nyquist sampling thereon, page24 */
if (t >= sampleTime) {
sampleTime += SAMPLEPERIOD ;
/* output to 2-D fft */
fftDataX[loopY][loopX] = filteredSignalX; /* real part */
fftDataY[loopY][loopX] = filteredSignalY; /* imaginary part */
loopX++;
} /* if */
} /* xforLoop */
} /* yforLoop */
} /* zForLoop */
Demodulator analysis (from Web):
Assume, rawNmrSignal = a(t) * cos((larmorFrequency + signal)t + phase) ;
localOscillator = cos(larmorFrequency * t) ;
demodulatedSignal = rawNmrSignal * LocalOscillator =
a(t)/2 * [ cos((2larmorFrequency + signal)t + phase) + cos(signal * t + phase) ] ;
And after LowPassFilter:
filteredSignal = a(t)/2 * cos(signal * t + phase) ; (result: larmorFrequency removed from original rawNmrSignal)
5. Image Reconstruction
Figure 3.7 Pseudo-Code for 2-D FFT
Adapted from Web, keywords: Paul Bourke, July 1998, 2 dimensional FFT
Algorithm: Transform each column (x), then transform each row (y)
for (loopY = 0; loopY < NUMBEROFVOXELSY; loopY++) {
for (loopX = 0; loopX < NUMBEROFVOXELSX; loopX++) {
real[loopX] = fftDataX[loopY][loopX];
imaginary[loopX] = fftDataY[loopY][loopX];
} /* xForLoop */
fft( real, imaginary, NUMBEROFVOXELSX);
for (loopX = 0; loopX < NUMBEROFVOXELSX; loopX++) {
fftDataX[loopY][loopX] = real[loopX] ;
fftDataY[loopY][loopX] = imaginary[loopX] ;
} /* xForLoop */
} /* yForLoop */
for (loopX = 0; loopX < NUMBEROFVOXELSX; loopX++) {
for (loopY = 0; loopY < NUMBEROFVOXELSY; loopY++) {
real[loopY] = fftDataX[loopY][loopX];
imaginary[loopY] = fftDataY[loopY][loopX];
} /* yForLoop */
fft( real, imaginary, NUMBEROFVOXELSY);
for (loopY = 0; loopY < NUMBEROFVOXELSY; loopY++) {
fftDataX[loopY][loopX] = real[loopY] ;
fftDataY[loopY][loopX] = imaginary[loopY] ;
} /* yForLoop */
} /* xForLoop */
Note: the data will probably need to be reorganized, as the default
organization from most FFT routines place the DC (frequency equals zero)
in the upper left hand corner, whereas the MRI picture expects
the DC value in the center of the picture. (In addition, it may be necessary
to reorganize the data depending on the order of how the phase data was collected. (In the literature the order of collecting phase data is called traversing K-space. see k-space below))
Simulation Results
The following is the result of a simulation, using data from the
BrainWeb Website (also see Cocosco_1997).
The original data was from a 3-D anatomical model of the brain, with dimensions: x:181mm, y:217mm, z:181mm.
(A discrete model, in which each voxel is composed of only one type of tissue was used.)
A slice was taken at z:90mm. The data was then scaled to x:128mm, y:128mm
and saved in a file: brainmri128.h.
After the simulation/fft the 128×128 picture was scaled up to 256×256 .gif picture:
MRI Simulation Result
Output generated from the MRI Simulator…
Glossary and Units:
- Contrast Resolution: The ability to distinguish between features with “slightly” different optical densities.
- FID: Free Induction Decay
- Field of View (FOV): The size of the 2 or 3 dimensional spatial area of the image.
- Image Plane Resolution: Resolution in the plane of the slice.
- Image Resolution: Same as spatial resolution?
- Longitudinal: MRI, Z direction
- Lorentzian (line): Fourier transform of a exponentially decaying sinusoid: Signal = A * exp(-kt) * sin(2pi * frequency * t); Bloch pages 88,141
- Slice Resolution: Thickness of the imaging slice.
- Spatial Resolution: How closely two features can be placed to one another and still be identified as different features.
- Specific Absorption Rate (SAR): Regulatory limits on tissue absorption of electromagnetic fields (Watts/kg)
- Transverse: MRI, XY plane
- Voxel: volume element
Prefixes
- mega: 1.0e6
- giga: 1.0e9
- tera: 1.0e12
- peta: 1.0e15
Units, Constants, Formulas
- 1 Angstrom = 0.1 nanometers (nm)
- 1 Tesla = 10,000 Gauss = 1 Weber/meter^2
- 1 M (molar) = 1 mole/liter = 1 umole/ul
- 1 ul = 1 mm3
- 1 ml = 1,000 mm3 = 1 cm3 = 1cc
- 1 Liter = 1 million mm3
- 1 Dalton = 1 molecular weight
- volume of sphere = 4/3 * pi * r^3
- surfaceArea of sphere = 4 * pi * r^2
- Avagadro’s number = 6.022e23; 1mole = 6.022e23 particles; 1umole = 6.022e17 particles
References:
- Cocosco1997: BrainWeb: Online Interface to a 3D MRI Simulated Brain Database; C.A. Cocosco, V. Kollokian, R.K.-S. Kwan, A.C. Evans; NeuroImage, vol.5, no.4, part 2/4, S425, 1997
- Merkle1989: Large Scale Analysis of Neural Structures, csl-89-10, [p89-00173]
Books/Monographs
Neural Imaging
- Abragam – The Principles of Nuclear Magnetism, 1961 (reprint 1999)
- Elster – Questions and Answers in MRI, 2001
- Goldman – Quantum Description of High-resolution Nmr In Liquids, 1991
- Haacke – MRI: Physical Principles and Sequence Design, 1999
- Horowitz – MRI physics for radiologists, 1994
- Jin – Electromagnetic Analysis and Design in MRI, 1999
- Liang – Principles of Magnetic Resonance Imaging, 2000
- Morris – Nuclear magnetic resonance imaging, 1986
- NRCIOM – Mathematics and Physics of Emerging Biomedical Imaging, 1996
- Rajan – MRI: a conceptual overview, 1997
- Slichter – Principles of Magnetic Resonance, 1990
- Stearns – Signal Processing Algorithms, 1988
- Vlaardingerbroek – MRI: theory and practice, 1999
Miscellaneous
- Feynman1964 – The Feynman Lectures on Physics (volumes I,II,III), 1964
Web Sites
- Link to BrainWeb @ mcgill