How specifically does an MRI machine build an image from received radio waves?

Unlike the excellent Wikipedia page on ultrasound imaging, the one on MRI only explains the principle theory behind MRI - that oscillating strong magnetic fields cause water molecules to emit radiowaves - without explaining how his is applied to build a detailed 3D image.

So, how do we get from the billions of excited hydrogen atoms spitting radio waves (presumably in all directions) to building up a 3D image... and what precisely does each 3D pixel record?

I little background for those interested - I want to be able to perform a "virtual MRI" of a computer-modelled patient. For x-ray and ultrasound I have enough understanding to do this but for MRI I don't.

• I feel as if this is more of an engineering question. – Aaron Stevens Aug 9 at 16:32
• One may look to Computed Tomography for some of the answers here. However CT is not physics per se, so it should be closed as a physics question. – Roy Simpson Aug 15 at 14:20

In MRI, an image is created by using gradient magnetic fields. By adding a gradient magnetic field the magnetic field is different at different positions in the body. The most important term in understanding the use of this is the larmor frequency. This is the frequency with which the hydrogen atoms will precess in a certain magnitude magnetic field and is proportional to the gyromagnetic ratio. A gradient magnetic field is first applied in the length of the patient (head to feet). Then, by applying a radiofrequency wave the atoms in exactly that slice of the patient for which the rf pulse has the larmor frequency are excited. This is called slice selection. From this point we know that all information must come from this one slice, so one of the three dimensions is known. The next step is to apply a gradient in lets say the left right direction. Because of the different larmor frequencies, the hydrogen atoms at a different lateral position will now precess at a different frequency, which is the same frequency at which they will emit a radiofrequency wave. So from the frequency of the recieved pulses you can know the second dimension of its source. This is called frequency encoding As for the third and last dimension a similar but slightly different technique is applied, called phase encoding. If you really want to get into it look that up, but for now you might want to start with understanding the first two dimensions.

The answer given by John is (partially) true for CT scans, but most certainly not for MRI. This is the reason that a CT scanner has a rotating head whereas an MRI has no moving parts. If you want to get into CT image building, look up filtered back projection.

• Thanks for your answer. Are you saying that when doing say the longitudinal measurement (let's say on the z-axis) at a particular frequency f, all we get from this is a collection of datapoints which lie in the plane z=d(f) - we have no information at all how they are distributed until we cross-reference the other directions? What exactly do you record... if millions of atoms in that slice are resonating surely you get an unholy mess of millions of overlapping signals? – Mr. Boy Jun 11 '14 at 9:42
• Exactly. The radiofrequency pulse you apply has such a frequency that only the atoms in a certain plane (slice) are excited. So any returning pulse can only origin from that slice. When the atoms relax, they emit a radiofrequency wave at exactly their larmor frequency. You collect these emitted waves, which indeed is a complete mess of different frequencies. With fourier analysis the second coordinate (lets say x). The last direction is done through phase encoding which, to be honest, i'd have to put some time in to read up on before i can explain it. – fhdrsdg Jun 11 '14 at 15:50
• What I don't get is how on earth you you use all that to map bits of signal to specific points. If you're not receiving on a grid/matrix like with x-ray, how can you isolate anything?! I guess I need to find the actual equations used, I thought they would be well known as it's pure physics but maybe they are someone's IP? – Mr. Boy Jun 12 '14 at 13:37
• Ok, what happens is that the phase encoding gradient (say x) is turned on, then off again. This makes that the precessing H atoms at different position x have a different phase. Next the frequency encoding gradient is turned on and the rf pulses form the H atoms are received. The amplitude of the received pulse is placed horizontally in a so called k space. By changing the strength of the frequency encoding gradient, the vertical position in the k space of the pulse is encoded. Finally, a 2D fourier transform is applied to the k space to get the image. – fhdrsdg Jun 12 '14 at 14:21
• If that didn't make any sense to you, try searching online for k space filling and check this out. As you've probably figured out by now, MRI is a much more complicated imaging technique than CT for instance. Good luck! – fhdrsdg Jun 12 '14 at 14:24

The key to understanding image generation in MRI lies in realizing that the signal sent to the antenna (or coil) by the patient's tissues includes two different types of information:

1. Information regarding the magnitude of the transverse magnetization of the tissue under the influence of a structured sequence of RF stimulation pulses. This magnitude depends on the tissue composition of every single voxel (3D pixel) of the anatomy being interrogated, and will be coded for medical interpretation by mapping it to gray-scale values to generate the final image.
2. The necessary spatial encoding that will allow tracing back each component of the signal to a specific location in physical space.

In the following diagram of an MRI of the brain in progress, two voxels (red and magenta) are highlighted, and their individual signals added up to form a resultant wave, filling in a line in k-space (Fourier space): The spatial information is induced via linear gradient fields. At any particular point in space, i.e. $$\color{red}{\vec r}$$ corresponding to the location of the red voxel, the precession frequency of the atoms of hydrogen is (rotating frame):

$$\omega =\gamma\, \vec G_z \cdot \color{red}{\vec r}$$

with $$\gamma$$ corresponding to the gyromagnetic ratio; and $$\vec G_z,$$ the 3D gradient of the magnetic field $$\vec G_z = \nabla B_z.$$

The phase of the hydrogen atoms at location $$\color{red}{\vec r},$$ is the time integral of

\begin{align} \phi(\color{red}{\vec r}, t) &= \int_0^t \omega(\color{red}{\vec r},\tau)\,\mathrm d\tau\\ &=\int_0^t \gamma\, \vec G_z(\tau) \cdot \color{red}{\vec r} \, \mathrm d\tau \\ &= \left( \gamma \int_0^t \vec G_z(\tau) \, \mathrm d\tau \right) \cdot \color{red}{\vec r} \\ &= \vec k(t) \cdot \color{red}{\vec r} \end{align}

where $$\vec k(t)=\gamma \, \int_0^t \vec G_z(\tau)\, \mathrm d\tau$$ is defined such that it parametrically defines a path through spatial frequency space.

The signal acquired in MRI is the sum of all transverse magnetization:

\begin{align} \color{purple}{\text{Signal}}(t) & \sim \int_{-\infty}^\infty M_{xy}(\vec r, t)\,\mathrm dx\mathrm dy\mathrm dz\\[2ex] &=\int_{-\infty}^\infty M_{xy}(\vec r,t)\, \mathrm d\vec r \\[2ex] & = \large \int_{-\infty}^\infty \underbrace{M_{xy}(\vec r, 0)}_{\color{blue}{\text {Image}}} \, \underbrace{\mathrm e^{-\mathrm 2\pi i \,\vec k(t) \,\cdot\, \vec r}}_{\color{orange}{\text{Phase rotation}}} \mathrm d\vec r \tag 1 \end{align}

The phase rotation depends on time and space and does not affect the magnitude of the magnetization (or signal).

Equation (1) is a Fourier transform relating $$\color{purple}{\text{Signal}}$$ as a function of time and the $$\color{blue}{\text {Image}}.$$

Each point in physical space (patient) will contribute one frequency component to the time domain signal. Conversely, each point in k-space represents the magnitude of one given spatial frequency over the whole image.

Critically, there is no FFT between signal reception ("Resultant signal" captured by the antenna (coil) in the diagram above) and the line in K-space being filled in. The signal received is already in Fourier space thanks to the action of the gradients that induce different spatial frequencies across the $$x$$- and $$y$$-axes in real space. Here is the diagrammatic depiction of the resulting phase differences induced by the application the phase and frequency encoding gradients used to obtain a line close to the center of K-space at the start of the $$180^o$$ train in an turbo (or fast) spin echo sequence (taken from this great animation by Tyler Moore on youtube): Notice the change along the $$y$$-axis (vertical waveform along the spins diagram in the left-upper quadrant) when a more peripheral line of K-space is obtained after the second refocusing pulse (right at the start of the readout): On the other hand, the spatial frequency along the $$x$$-axis is identical to the prior line of K-space at the same point in time.

A finer point to understand is the need for as many phase encoding steps as the matrix in the phase encoding direction. The explanation for these repetitive experiments, each filling in a single line of K-space, is the fact that for every single step the spins colored the same on the diagram below will only differ in phase - not in frequency - as they go through the different lobes of the frequency encoding gradients, because the phase gradient is not turned on during readout, and the spins return to processing at the same frequency, albeit with residual phase shifts. When two constituent sinusoids contributing to the signal emitted by the patient have the same frequency, and vary only in phase, they cannot be distinguished. A single sinusoidal of a different phase and amplitude will be received, losing the information contingent solely on phase differences.

Finally, the signal extracted from the patient during the application of the middle lobe of the readout gradient can be conceptualized by simplifying the process, and thinking of it as the aggregate of different discrete sine waves: the first ones captured in the temporal line correspond to the highest frequencies collected, which are replaced by progressively slower sine waves of increasing amplitude until they peak at the point of maximum coherence in the transverse magnetization of the spins after the $$180^0$$ pulse. After that the process is reverted: The component frequencies captured map to points in K space. The ADC transforms the analog signal, and the spacing between samples is referred to as the dwell time (sampling rate), whereas the reciprocal of the dwell time $$(\Delta t)$$ is the receiver bandwidth $$(\text{RBW}=1/\Delta t).$$ Each frequency can be sampled by the receiver for a shorter or longer period of time, and the longer each frequency represented in K space is sampled (lower bandwidth), the higher the signal-to-noise (SNR) will be, at the cost of a slower readout rate.

An ADC effectively averages the input signal over the dwell time, thus the effective noise on the digitized signal is proportional to $$1/\sqrt{\Delta t}$$ making it desirabel to increase $$\Delta t.$$ Signal noise [sic] is also directly proportional to the receiver bandwidth. Richard Ansorge and Martin Graves - The Physics and Mathematics of MRI.

The Nyquist theorem imposes restrictions on the number of samples needed to capture a given frequency, but since two quadrature samples are taken at each $$\Delta t$$ interval, here would be possible sample of a matrix with 5 acquisition points along the frequency direction: [Please note that on the image above, sampling is carried out only on the first part of the echo - this is theoretically possible given the symmetry of the signal, although in practice, more than 50% of the echo is sampled. Also in partial-echo techniques it is last part of the signal that is sampled to decrease the TE by shortening the readout gradient.]

Contrarily, a higher bandwidth will sample more points per msec at the expense of lower SNR, but with the upside of obtaining more closely packed echoes (shorter echo spacing, ESP), enabling longer ETL's for a given TR: $$\text{ESP}\sim \frac{\text{freq encoding steps}}{\text{receiver bandwidth}}.$$

Resources:

A K-space analysis of Small-Tip-Angle-Excitation by John Pauly et al. JMR 1989

MRI: Introduction to K-space by Eric Wong

Fast Spin Echo - Spaces Animation by Tyler Moore

• This is pretty useful information (have a bounty!), though I wonder what audience you were writing for and whether the text is really accessible to them. – Emilio Pisanty Aug 5 at 13:30
• @EmilioPisanty Thank you, Emilio! The audience I have in mind is that Venn diagram intersection between the radiologist in practice whose curiosity in the topic is not fully satiated with the equation-free intuitions usually found in MRI textbooks for physicians, on the one hand; and the basic sciences / mathematically savvy who know a thousand time more than me about real analysis, but who may be encountering MRI as a clinical tool for the first time. I'm not sure about the second subset, but I have a pretty good idea of the first. – Antoni Parellada Aug 5 at 13:34
• stellar answer! – lurscher Sep 16 at 16:10

The answer has grown and grown in my attempt to use words rather than mathematically formulae with the comment made by @Emilio Pisanty, . . . though I wonder what audience you were writing for and whether the text is really accessible to them.

I think that describing the MRI scanner is difficult for several reasons which include the fact that it is 3D object on top of which you need to include a fourth dimension, time evolution. At each step in any description there are complications which if addressed diverts one from the main thrust/overview of explaining how the final image is produced. It is for this reason that I intend to produce an overview and as appropriate provide more detail.

To answer the question, How specifically does an MRI machine build an image from received radio waves? I have omitted much of the detail as to how the signals which are used for the construction of the final image are produced.

Nuclei have a magnetic moment due to their spin. In general, for a large number of nuclei their magnetic moments are orientated randomly and so the thermal dynamic equilibrium magnetisation of a specimen $$M$$ is zero.

If many nuclei are exposed to a uniform magnetic field along the z-axis overall there is a thermal dynamic equilibrium magnetisation $$M$$ produced with the magnetisation vector $$\vec M$$ pointing along the z-axis.

For the isotope $$^1_1H$$ the gyromagnetic ratio is particularly high and as hydrogen is an abundant ingredient in the human body in the form of water, MRI scanners generally are designed to detect radio frequency waves, frequency $$\approx 100\, \rm MHz$$, the magnetic resonance signal, produced as a result of the precession of hydrogen nuclei in a magnetic field.

A magnetization vector $$\vec M$$ in a magnetic field $$\vec B$$ undergoes precession as described by the equation $$\frac{d\vec M}{dt}= \gamma \vec M \times \vec B$$ with the frequency of precession, called the Larmor frequency, $$\omega_L = \gamma \, B$$ where $$\gamma$$ is the gyromagnetic ratio and the applied magnetic field $$B$$.

When a large constant external magnetic field, $$B_0 \approx \rm tesla$$, is applied is much easier to continue the description of MRI in a frame of reference which is rotating about the $$z$$ axis at the Larmor frequency, $$\gamma B_0$$ because in this rotating frame it is as though there is no external magnetic field $$B_0$$ and no precession of a magnetisation vector unless an extra (perturbing) magnetic field added.
To change the magnetisation of the specimen under investigation a pulse of the appropriate frequency with a sinc envelope is applied to the specimen via the rf coils. This is the perturbation which produces a decaying precession of the protons in the hydrogen atoms. The magnetisation undergoes free induction decay at the Larmor frequency towards the original thermal dynamic equilibrium magnetisation. Whist undergoing the decay radio frequency waves are emitted which are picked up by the same rf coils as were responsible for producing the initial change in magnetisation of the specimen. The decay is characterised by relaxation times $$T1,\,T2$$ and $$T2^*$$ which depend on where the hydrogen atoms (protons) are sited eg water, fat etc. and measurement of these relaxation times enables the identification of the substance (tissue) in which the decay occurred. An mri scanner picks up the rf signals from the decay and can identify the substance and its location responsible for production of the rf signal.

There are techniques used, like echo formation, which are used to ensure that the rf signal produced during the decay process have as large an amplitude as possible, but these are details in the context of this answer.

A rf pulse at the Larmor frequency is applied to a specimen which results in a change in the magnetisation of the whole specimen with hydrogen nuclei precessing at the Larmor frequency ie it is a resonance. As time progresses the magnetisation changes and electromagnetic waves at the Larmor frequency are produced. The waves from all the hydrogen nuclei in the specimen are picked up by the rf coils.

Given that the Larmor frequency is $$\gamma B$$ why not apply additional magnetic field superimposed on the large static magnetic field $$B_0$$ so that each part of the specimen is in a different value magnetic field ie the Larmor frequency and hence the emitted signals from different parts of the specimen are all different. This would mean that, out of the chaos of all parts of the specimen producing a decaying precession signal of the same frequency, each part of the specimen would produce a unique frequency decaying precession signal which could be relatively easily processed.

To change the frequency of precession and hence the frequency of the emitted signal gradient magnetic field in the z-direction can be added to the static magnetic field $$B_0$$. As a start let an additional field $$G_{\rm z} z$$ which varies with position $$z$$ be added. This would mean that within a slice of the specimen at position $$z$$ all the Larmor frequencies and hence the decaying precession signal frequencies would be the same, $$\gamma G_{\rm z} z$$. By making the rf sinc excitation pulse be at this frequency only those hydrogen nuclei in that slice will be excited. On top of the gradient field which varies in the z-direction now add a gradient field $$G_{\rm x}x$$ which varies in the x-direction so now the Larmor frequency and the decaying precision frequency vary with both position coordinates $$z$$ and $$x$$, $$\gamma (G_{\rm z} z+G_{\rm x}x)$$.

Using a band pass filter each of these individual frequencies could be analysed and so the decaying precession signals from parts of the specimen for a given $$z$$ and $$x$$ could be analysed. These signals would be the sum of the results of the precession decay over the whole range of $$y$$ values. All that needs to be done now is to assign unique Larmor frequencies for all positions $$(x,y,z)$$ by applying a gradient field in the y-direction.
This is where there is a stumbling block. When using three gradient fields it is impossible to produce unique Larmor frequencies for all positions $$(x,y,z)$$ as shown below. The numbers represent the change in the magnetic field in the z-direction due to field magnetic gradients in the z and x-directions. Chose the z-slice and the x-column where the field is larger by $$+2$$ due to the x-gradient field. Now impose a y-gradient field and the problem becomes apparent. The Larmor frequencies characterised by the magnetic fields in the x-column under consideration, $$-1,0,1,2,3$$, are not unique to that column. This means that if one chose the Larmor frequency characterised by a total magnetic field difference of $$-1$$ there are four locations which produce decaying precession signal of that frequency. In early days this method was used and the ambiguity in position was overcome in a variety of ways including moving the specimen relative to the scanner and taking multiple readings.

This is where phase encoding (signal related to phase) comes in but please note it is not really so much different from the frequency encoding (signal related to frequency). In most text it is the x-axis which is labelled frequency and the y-axis phase, but it could be just as easily be labelled the other way around. For a lot of us it is the phase encoding which is the real stumbling block in the understanding of how an mri scanner works.

The process by which the phase encoded information is obtained is as follows.
A rf sinc pulse at the precessional frequency of the z-slice excites the protons and they start precessing at the Larmor frequency. The z and x gradient fields are applied and the output from the rf coils stored. The process is repeated but now with a small y gradient field $$G_{\rm y}$$ applied for a time $$\tau$$ and the output recorded. The process is then repeated with larger and larger y gradient fields with the number of times this process is repeated depended on the requirements of the final resolution of the image.
When there is no y-gradient field the Larmor frequency of all the protons is a given x-column are the same. Application of the y-gradient field for a fixed time $$\tau$$ changes the Larmor frequency and so the rate of precession in the x-column is different for each y-position in that column. There is a change of phase between the precessions at different positions $$y$$. When the y-gradient field is switched off after a time $$\tau$$ the precession frequency goes back to its value before the y-gradient field was switched on and the phase difference between at difference y-positions stay the same. The phase difference in the x-column is proportional to $$G_{\rm y} \,y\,\tau$$. If now the y-gradient field is doubled and applied for the same time $$\tau$$ the final difference in phase between y-positions is also doubled.

To measure the phase difference requires more than one reading.
Suppose that there are two sources of amplitude $$A$$ and $$B$$ and phase difference $$\phi$$ and only the amplitude sum of the two waves waves can be measured.
The phasor diagram above shows that the amplitude of the sum $$X$$ can be found from application of the cosine rule $$X^2 = A^2 +B^2 +2AB \cos \phi$$ but on measuring the amplitude all information about the phase difference is lost.
However now suppose that the phase difference is doubled to $$2\phi$$ then the resulting amplitude $$Y$$ is now given by $$Y^2 = A^2 +B^2 +2AB \cos 2\phi$$.

If only the relative magnitudes of $$A$$ and $$B$$ are needed then one cam make $$A=1$$ and there are two equations and two unknowns $$B$$ and $$\phi$$ so $$\phi$$ can be found. With more sources one needs to advance the phase more times to get the phase differences.

The effect of repeating than scanning process to get the phase information is shown in the series of images below. Within the z-slice the time evolution of the phase of the magnetisation at $$x,y$$ is given by $$\Delta \phi (t) = \gamma G_{\rm x} x\,t + \gamma G_{\rm y} y \,\tau$$, the first term being the frequency term and the second term being the phase term.

If the density function which is going to be displayed as the final image is $$\rho(x,y)$$ then as a function of time the signal received from the whole z-slice will be $$\displaystyle \iint \rho(x,y)e^{i\gamma G_{\rm x}t x + i\gamma G_{\rm y}\tau y}dx dy$$ which can be written as $$\iint \rho(x,y)e^{-i k_x x- ik_y y}dx dy$$ which resembles a 2D Fourier integral. Here $$k_{\rm x} = \gamma G_{\rm x} t$$ and $$k_{\rm y} = \gamma G_{\rm y} \tau$$. The information about the final image is stored in a 2D k-space and from that the final image needs to be reconstructed.

In one dimension to produce a square wave the following summation has to be done $$\sin kx + \frac 13 \sin 3kx + \frac 15 \sin 5kx . . . .$$ and this is illustrated in the gif image below.

The right hand diagram being a representation of the square wave in a 1D k-space showing the amplitudes (y-axis) of the k values (x-axis).

The output from an MRI scanner starts off as a representation in 2D k-space and what needs to be done is to sum the amplitudes over all the k-values to get the final image.

The process is as follows. Choose a point in k-space eg blue dot.
Measure length of k-vector from the origin. Evaluate the reciprocal of the length of the k-vector. This is the value of the "wavelength" of the plane wave which you are going to construct.
The direction of the plane wave fronts will be at right angles to the direction of the k-vector.
A section through the plane wave will be at right angles to the xy-plane with an amplitude as encoded as a value (greyscale in the image)in k-space.

Repeat for all the k-plane data points adding together the plane waves algebraically as was done to produce the 1D square wave.
The final image will appear as shown top right image below.

If the "amplitude" at a particular k-value is made artificially larger the resulting wave which has an amplitude which is larger than it should be is seen as part of the final image.

In generating this answer I have used many resources but would like to list two of particular note.

For the image reconstruction I think that this is a fabulous video.

Although not the easiest resource to use this website Questions & Answers in MRI has lots of very useful information and many very good links to other webpages.

The Tyler Moore animations I found useful after I have understood the basics of what was going on.

There is also an edX course Fundamentals of Biomedical Imaging : Magnetic Resonance Imaging (MRI) which worth looking at.