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.
 A: 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.
A: 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:

*

*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.

*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 describes 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, t=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}}.$
Decomposing equation $(1)$:
$$\text{Slice signal}_{k_x,k_y}=\int_{\text{slice}}\int_{\text{slice}}\rho_{(x,y)}\rm e^{-2\pi i \color{blue}{k_x} x} \;\rm e^{-2\pi i \color{red}{K_y} y}\rm dx \rm dy\tag 2$$
where $\color{blue}{\small \text K_x = \frac{\gamma}{2\pi}G_x \times t}$ and $\color{red}{\small \text K_y = \frac{\gamma}{2\pi}G_y \times \tau_{pe}},$ in which $\tau$ is the time for which the phase encoding gradient is applied, $\rho_{(x,y)}$ are the characteristics of the tissue in the anatomy we are trying to image, and $\gamma$ is the gyromagnetic ratio. This expression clearly relates the variable time in the read-out gradient with the values of the complex sinusoids obtained by integrating at each point in time over the entire slice.
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 (inspired by this great animation by Tyler Moore on youtube) of the resulting phase differences induced by the application the phase and frequency encoding gradients used to obtain a line of K-space in a spin echo sequence ($*$ applicable to turbo spin echo - see below for explanation of turbo spin echo):
Initially, the slice selection gradient is turned on during the application of the $90^\circ$ (excitation) and $180^\circ$ (refocusing) pulses. There is no travelling in k-space yet:

In the next step, the phase-encoding gradient is turned on: the gradient strength or amplitude as well as the time during which the gradient is applied determine the distance away from the $x$-axis the line about to be read will be in k-space because of the equation $\color{red}{\small \text K_y = \frac{\gamma}{2\pi}G_y \times \tau_{pe}}$ above $(2)$ (this is schematized as a red diamond travelling up along the $y$-axis in k-space while the phase gradient is on):

A pre-phasing lobe of the frequency-direction gradient will be turned on right after the phase-encoding gradient is applied. The amplitude or intensity of the gradient will determine the bandwidth of frequencies (given a fixed FoV). The longer it is applied, the more dephasing between anatomical points along the $x$-axis (frequency). During the application of the pre-phasing lobe, the position along k-space will become negative getting ready for the acquisition of a line of from left to right (purple arrowhead pointing left):

At the beginning of the read-out gradient, the first data point collected corresponds to the maximum negative value of $\color{blue}{\small \text K_x = \frac{\gamma}{2\pi}G_x \times t}$ due to the negative prephasing (or dephasing lobe) shown above, the second data point to a slightly more positive value of $\small \text k_x$ (arrowhead pointing right):

At the $\text {TE}$ point the spins still subject the same gradient amplitude (although reverse sign to the prephasing gradient) will have had a chance to revert the initial dephasing effect caused by the initial negative frequency prephasing gradient, showning no residual dephasing in the frequency direction - at that point only the effect exerted by the initial phase gradient will be apparent. The signal returned to the antenna will have maximum intensity:

During the second half of the read-out gradient dephasing in the opposite direction will increase over time and cause loss of signal strength in a mirror-like fashion that can be harnessed to save acquisition time. In practice, more than $50\%$ of the echo is sampled, but not all of it is needed. In partial-echo techniques it is the last part of the signal that is sampled to decrease the TE by shortening the readout gradient.

After the line of k-space is completed, a negative lobe in the frequency direction will be again applied to rephase the spins:

The dephasing along the phase initially introduced will be now reverted to return to the center of k-space:

After that, the process will be repeated to obtain another line of k-space (within the same TR in TSE, or in a different TR in classic spin-echo) (second green line in k-space below). The order in which k-space will be filled-in is not random, but rather follows a strategy (for instance 'linear' in which the center of k-space ($\text k_0$ or "contrast line" is acquired in the middle of the read-out; or 'low-high', in which the contrast line is acquired first, followed by more peripheral lines):

There is a need for as many phase encoding steps as the matrix in the phase encoding direction. The explanation is that 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 (see here).
I have tried to reproduce the process of filling k-space in this applet, illustrating a single-echo SE acquisition.

In Eq $(2)$ the signal is expressed as a function of each coordinate in k-space, which in turn depends on the readout time $t$ and the phase gradient applied $\tau_{pe}.$ Going back to $\color{blue}{\small \text K_x = \frac{\gamma}{2\pi}G_x \times t}$ the strength of the gradient will remain unchanged during the acquisition of a line of k-space, and therefore will the difference in precession frequency between adjacent points in the anatomy in the frequency direction, or the actual frequency at each anatomical point. During the read-out all the frequencies from each anatomical point within the slice are sampled and digitized by the ADC (see below) to produce a point in k space.
A true-to-fact visualization of the process can be as follows: Consider three tissue cells anatomically separate in the body. The red cell is part of the liver and happens to lie at one end of the gradient (for simplicity, just the frequency gradient $\vec G_x$) with its hydrogen atoms precessing at a frequency much higher than a blue cell in the left kidney (opposite side of the belly), which spins are slowest. In between an orange cell in the pancreas will precess at an intermediate frequency. Essentially, and throughout the emission of the signal back from the patient, these three cells will act as radio stations on the FM dial, emitting at a fixed frequencies. The intensity of the signal for each one of these three cells will increase progressively to the middle of the echo without changing their individual frequency, only to decrease progressively after that point in a symmetrical fashion.
The ADC (see below) will listen to the sum of all the frequencies coming back from the slice (red, orange and blue) at the same time in each digitization sample, collected at regular intervals (dwell time) (see below). Each sample will be directly allocated to the proper position on k-space according to complex plane coordinates, and only at a later point a FFT (or wavelet transform) will generate the final image.

Given that the integral of frequency over time is the definition of phase, it will actually be the dephasing along the phase and frequency directions induced by the gradients that will produce the signal. This is more clearly apparent in the expression on page 41 of Richard Ansorge's The Physics and Mathematics of MRI, in which the signal becomes a function solely of time:
$$\begin{align}
S(t)&=\int\int\int_{\text{FoV}} \rho(\vec r)\; \mathrm e^{-\mathrm i \omega_0 - \mathrm i \theta_{\text{acc}}}\mathrm dx\mathrm dy \mathrm dz\\[2ex]
&= \mathrm e^{-\mathrm i \omega_0 t } \int\int\int_{\text{FoV}} \rho(\vec r)\; \mathrm e^{-\mathrm i 2\pi \vec k \cdot \vec r}\mathrm dx\mathrm dy \mathrm dz
\end{align}$$

The effect of a particular gradient is to cause a local additional instantaneous precession $\gamma_p  {\bf\vec G}(t) \cdot \bf \vec r,$ thus if the transverse spins were in phase at some time $t_0$ (say due to a $90^\circ$ RF pulse) than at a later time, $t,$ the accumulated phase offset $\phi_{\text{acc}}(\bf\vec r,t)$ is given by

$$\theta_{\text{acc}}({\bf \vec r}, t)= \gamma_p \Bigg[\color{orange}{\int_{t_0}^t \vec {\bf G}(\tau)\mathrm d\tau} \Bigg]\cdot \bf \vec r= 2\pi \;\color{orange}{\bf\vec k} \cdot \bf \vec r;$$
where $\rho(\vec r)$ is the proton density in different locations; and $\omega_0$ is the Larmor frequency induced by the main magnetic field $B_0.$

The amplitude and phase of the signal will be treated as real and imaginary components in the receiver. The ADC transforms the analog to digital signal:

The ADC rating depends on the sample rate or sampling frequency (how closely the samples can be taken) and the amplitude or resolution (for example 16 or 24 bits). The actual process calls for a preamplification and a demodulation of the Larmor carrier frequency $(\small 50–300 \text{MHz})$ to be able to capture the information in the MR signal contained within a small bandwidth $\small \leq 1–2 \text{ MHz}),$ determined by the maximum gradient strength and the field of view (FOV), before the ADC conversion.
At this point in the process a filter (cosine, cosine squared, Fermi,
Gaussian, Hamming, Lorentzian, Riesz, Tukey,or others) operating in k-space will be applied to remove salt-and-pepper noise (SPN) occurring during the process of acquisition and transmission.
The spacing between samples is referred to as the dwell time (time between digitized samples), whereas the reciprocal of the dwell time $\small (\text{T}_s=\Delta t)$ is the sampling bandwidth $\small (\text{RBW}=1/\Delta t),$ or sampling rate $(\text f_s).$
The sampling bandwidth and the number of sampled points (matrix size in the frequency direction) determine the sampling time: $\small \text{Total sampling time} =\text{Matrix}_{\text{fr.}} \times \frac{1}{\text{RBw}} =\text{Matrix}_{\text{fr.}} \times \text{Dwell time}.$
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 $\small 1/\sqrt{\Delta t}$ making it desirable to increase $\small \Delta t.$ Signal noise [sic] is also directly proportional to the receiver bandwidth. Richard Ansorge and Martin Graves - The Physics and Mathematics of MRI.

A higher bandwidth results in reductions in the minimum $\small \text{TR}$ and $\small \text{TE},$ potentially decreasing scan time (provided the contrast afforded by minimal $\small \text{TR}$ values are medically desirable). A higher receiver bandwidth will sample more points per msec at the expense of lower SNR (see here), but with the upside of obtaining more closely packed echoes (shorter echo spacing, ESP), allowing a longer ETL's for a given $\small\text{TR}$: $\small \text{ESP}\sim \frac{\text{freq encoding steps}}{\text{sampling bandwidth}}.$ Reducing the acquisition time by increasing the bandwidth is not typically done, although there are some exceptions, such as 3D gradient (FFE) sequences here: although a higher bandwidth will decrease the sampling time, it is going to make a difference only as long as it allows a shorter minimum $\small \text{TR},$ and that lower $\small \text{TR}$ is desirable to achieve the contrast in the image intended - otherwise the few milliseconds saved are not comparable to the $\small \text{TR}.$
From a different perspective, the receiver bandwidth is defined as $\small \text{BW} =\text f_s \times  \text{Matrix}_{\text{fr}}=\text f_s \times  \text{No. samples}.$ This definition is different from the one above (sampling rate). There is a third concept of bandwidth as the range of precession frequencies that need be included for a chosen FoV:
$$\text{rBW (range of freq)}= \frac{\gamma}{2\pi} \times\text G_x \times \text{FoV} $$

Both expressions will match: sampling at the Nyquist rate, and given a matrix or number of points to sample, $\small \text{T}_s=\text{Pixels (no. sampled points)} / (2\times \text{max freq}).$
The gradient steepness and the bandwidth will have to move in synch unless the FoV is to be changed. On the other hand, if a change of FoV is introduced by the operator, the MRI unit will first try to achieve the desired FoV by changing the gradient. Only if changing the gradient is not enough to achieve the target FoV at a given matrix due to engineering limits, will the rBW be altered. The reason for this hierarchy of tasks is that changing the rBW will also alter the SNR or introduce chemical shift.

$(*)$ In a turbo spin echo (TSE) (same as fast spin echo (FSE)) each excitation pulse (90-degree pulse) is followed by a series of equally-spaced refocusing 180-degree pulses (echo-train length or ETL) and a corresponding echo will be collected after each of these refocusing pulses. This is just a strategy to save time by extracting several echoes out of every excitation pulse instead of the single-echo scheme of a conventional SE. From The Physics and Mathematics of MRI a CSE sequence would be of the form:

The refocusing $\small 180^\circ$ pulse is meant to re-phase the spins, causing them to regain coherence and thereby to recover transverse magnetization, which is really what is being captured. The $T_2^*$ of free-induction decay (FID) (see illustration above) is thus modified into a less rapid decay. The $\small 180^\circ$-pulse is half-way between the $\small 90^\circ$ pulse and the echo.
On a FSE sequence the diagram is:

The acquisition time that would be calculated in conventional spin echo (CSE) as $\small\text{TR} \times \text{no. phase encoding steps N}_y \times \text{no. of averages (NEX)}$ indicating that to fill in k-space the number of different excitation pulses necessary is now reduced by a factor of ETL:
$$\text{Acquisition time}=\frac{ \text{TR} \times \text{matrix (Ph.)} \times \text{NEX}\times \text{no. packages}}{\text{ETL}\times \text{SENSE accel.}}\tag 3$$
The concept of package is explained below. SENSE (and compressed SENSE) are not explained.
Clinical MR exams are obtained as a set of stacks of images through the anatomy of interest at different desired planes and angulations and with different types of contrast between normal and diseased tissues and organs (e.g. T1, T2, proton density (PD)) $(***)$. Each of these stacks of images have identical parameters, and are called sequences.
Typically, the slices or images in a given sequence will be obtained in packages (interleaved) - approximately $5$ or $6$ for a $\text{T}1$-weighted sequence and $2$ or $3$ for a  $\text{T}2$ or $\text{PD}$ sequence. This takes advantage of the dead time imposed by the $\text{TR}$ or time between $\small 90^\circ \text{ pulses}.$ For example, if the $\text{TE}$ is $80$ msec, and the $\text{TR}$ $3,000$ msec, depending on the echo train there may be some $2,900$ msec of dead time, which can be used to consecutively select other slices (slice selection gradient) and proceed with the sequence of $\small 90^\circ$ excitation pulse, followed by the corresponding train of $\small 180^\circ$ pulses in as many slices as can be fit in this otherwise idle time interval. Within each of these packages or $\text{TR}$-sets each of the slices included will have a number of lines of k-space recorded equal to the echo train length. Importantly, these slices are not consecutive to avoid crosstalk. In a $2$-package sequence, the slices will be separate by one: $\{1,3,5, \dots\}$ for package $1$ and $\{2,4,6, \dots\}$ for package $2.$ For a $3$-package sequence, $2$ anatomically consecutive slices will be skipped.
The sequence of events is: With the slice selection gradient on, a $\small 90^\circ \text{ pulse}$ for the first slice in a package starts the "timer" for the repetition time $\small \text{TR}.$ This is followed by a train of $\small 180^\circ$ refocusing pulses and corresponding echo collection to fill a number of lines in k-space equal to the train length for the first image in the package. Then the same happens for the second image (anatomically separate) and so forth and so on until the last image in the package. At that point the "timer" for the $\small\text{TR}$ stops as a new $\small 90^\circ \text{ pulse}$ is emitted to selectively excite slice the first image all over again to obtain a new number of lines equal to the train length to continue filling k-space for as many lines as the matrix in the phase direction dictates. Hence $\small\text{TR}$ is the distance between two consecutive $\small 90^\circ$ pulses, and in every $\small \text{TR}$ a number of lines equal to the train length will be added to k-space for each of the images contained in the same package.

The need to selectively excite each of the slices and sample the corresponding echoes in the train limits how short the $\small \text{TR}$ can be set by the MRI unit. A separate issue limiting how short the $\small \text{TR}$ can go to save time is the SAR (Specific Absorption Rate) leading to tissue heating.

$(**)$ The SNR doesn't solely depend on the bandwidth. In fact, the SNR is controlled by
$$\small \text{SNR} \sim \text{VS} \times \sqrt{\frac{\text{No. Phase steps (Matrix)}\times\text{ No. Averages (NEX)}}{\text{rBW} } }\tag 4$$

$(***)$ The signal is generated by the transverse tipping of the hydrogen spins. The transverse magnetization decays over time as the longitudinal magnetization along $B_0$ recovers, as encapsulated in the Bloch equations:

Different normal and pathologic tissues differ in $\small \text T1, \; \text T2$ or concentration of hydrogen atoms (proton density or $\small \text{PD}).$ This is exploited to obtain different types of diagnostic contrast on the MRI images. To do so the repetition and echo times are adjusted to capture maximal differences (from here):


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
The Physics and Mathematics of MRI by Richard Ansonge and Michael Graves
MRI Physics: MRI Image Formation Parameters, SNR, and Artifacts
A: 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.
A: Just for fun, you might read on MRI Diffusion Tensor Imaging. https://www.news-medical.net/health/Diffusion-Tensor-Imaging-(DTI)-Explained.aspx Rather than looking at the density of a tissue voxel, it looks at the direction of the voxel's orientation WRT water diffusion. Water diffuses along the axon between excitation and re-emission. This is magically interpreted to display the long-distance path of tracts in the brain.

I find the images breathtakingly beautiful.
