CA2147556C - Method for reconstruction for mri signals - Google Patents

Abstract

A method and an apparatus are used in magnetic resonance imaging. In generating the resonance signals, the magnetic gradient fields are varied to represent a scan on a non-linear, continuously varying curve. One specific example given is a multi-leafed rose. New mathematical treatments of the resultant signal are presented, which are used to define a novel method for the reconstruction of an image from these signals based on k-plane coordinates naturally associated with the signal acquisition. The superior image quality that the new reconstruction method gives means that non-linear, continuously varying acquisition in k-space may be used so that magnetic field switching and high cost hardware is not required, while signal acquisition times are significantly reduced.

Description

METHOD FOR RECONSTRUCTION FOR MRI SIGNALS
FIELD OF INVENTION
The present invention relates to image processing and more particularly to the reconstruction of magnetic resonance images.
BACKGROUND
Magnetic resonance imaging (MRI) is an application of nuclear magnetic resonance (NMR) techniques to the production of images of objects, especially living tissue. In conducting MRI, a uniform magnetic field is applied to the object being imaged. Magnetic dipoles of nuclei in the object tend to align with the magnetic field. Resonance excitation pulses are then applied to the image region of the object to cause the magnetic dipoles to precess about the strong magnetic field.
As the dipoles precess, their magnetic lines cut the wires of a receiver coil, producing resonance signals. Gradient magnetic fields are applied transversely to the main field to encode in two-dimensional frequency space, or k-space, the resonance signals in accordance with the spatial positions of the dipoles. The resonance signals are then used to reconstruct the desired image.

Various techniques are currently used for generating the signals and reconstructing the images. These include Fourier Imaging and Spin Warp Imaging in which the image is reconstructed from parallel lines in the spatial frequency domain (the k-plane) on an image or picture plane along parallel lines, actually at discrete points in a rectangular grid using the discrete Fourier transform. In Zeugmatography the image is reconstructed from radial, angularly spaced lines in the k-plane generally using backprojection methods. In each of these cases, each line is a single signal acquisition. It is necessary to allow a significant time to elapse between signal acquisitions in order to allow spins to realign with the main magnetic field.
The Echo Planar Imaging (EPI) technique has been developed to address the problem of lengthy periods between signal acquistions. This technique uses rapidly switching gradient magnetic fields during measurement time to produce a series of straight line segments across the k-plane. This covers enough of the k-plane in one or two passes for image reconstruction via the discrete Fourier transform. The EPI
techniques used include the Fast Low-angle Excitation Echo planar Technique (FLEET), the Blipped Echo planar Single pulse Technique (BEST) and Modulus BEST (MBEST). There is an inherent difficulty in rapidly switching a magnetic field and the hardware for producing the switching gradient fields for EPI is consequently very expensive.
The present invention applies to alternative data acquisition techniques that do not suffer from the lengthy periods between signal acquisitions found in the Fourier Imaging, Spin Warp Imaging and Zeugmatography techniques, and does not require the complex and expensive switching equipment used in EPI. Higher quality images from these techniques are now possible because of the applicant's development of new mathematical techniques that will allow high quality reconstruction of an image from an MRI signal representing a continuous, non-linear curve on the k-plane. This means that gradient fields can be modulated continuously rather than discontinuously switched, which involves simpler hardware and less power consumption than is required for EPI.
SUMMARY
Thus, according to one aspect of the present invention there is provided a method of magnetic resonance imaging an image region of an object, said method comprising:
generating a main magnetic field through the image region of the object to align magnetic dipoles of nuclei in the image region of the object with the main magnetic field;

exciting precession of the magnetic dipoles of the nuclei about the main magnetic field such that each nucleus generates a component of a resonance signal;
producing magnetic field gradients across the main magnetic field so as to vary the temporal phase and frequency of the resonance signal components to encode two-dimensional spatial frequencies according to the positions of the nuclei in the image region;
varying the magnetic field gradients to generate a two-dimensional, non-linear, continuous scanning curve;
monitoring the resonance signal; and operating on the resonance signal using a novel transform based on natural coordinates in k-space and in the use of the corresponding Jacobian to weight the sampled resonance signal as part of the transform to convert the resonance signal into a representation of the spatial positions and densities of the imaged nuclei.
The scanning curve may have various forms. Curves that are suitable include multiple leaf roses, spirals, sinusoidal wave forms and others.
According to another aspect of the present invention there is provided an apparatus for magnetic resonance imaging of an image region of an object, said apparatus comprising main field generator means for generating a main magnetic field through the image region of the object for aligning magnetic dipoles of nuclei in the image region of the object with the main magnetic field;

5 excitation means for exciting precession of the magnetic dipoles of the nuclei about the main magnetic field such that each nucleus generates a component of a resonance signal;
gradient field generating means for producing magnetic field gradients across the main magnetic field so as to vary the temporal phase and frequency of the resonance signal components to encode two-dimensional spatial frequencies according to the positions of the nuclei in the image region;
gradient field varying means for varying the gradient magnetic fields with time to produce a two-dimensional, continuous, non-linear scanning curve;
receiving means for monitoring the resonance signal; and computer means for operating on the resonance signal using a novel transform based on natural coordinates in k-space and in the use of the corresponding Jacobian to weight the sampled resonance signal as part of the transform to convert the resonance signal into a representation of the spatial positions and densities of the nuclei.
In the following, there will be given a description of the concepts involved in magnetic resonance imaging using current techniques and a description of the mathematical principles on which the present invention is based.
BRIEF DESCRIPTION OF THE DRAWINGS
In the accompanying drawings, which illustrate exemplary embodiments of the present invention:
Figure 1 is a schematic representation showing the arrangement of magnets in a conventional MRI machine;
Figure 2 is a simplified block diagram of the MRI receiver system;
Figure 3 illustrates the k-plane curves associated with Spin Warp or Fourier Imaging;
Figure 4 illustrates the curves associated with Zeugmatography;
Figure 5 illustrates a curve associated with FLEET Echo Planar Imaging;
Figure 6 illustrates the gradient fields associated with BEST Echo Planar Imaging;

Figure 7 illustrates the k-plane curve associated with BEST Echo Planar Imaging;
Figure 8 illustrates the gradient fields associated with MBEST Echo Planar Imaging;
Figure 9 illustrates the k-plane curve associated with MBEST Echo Planar Imaging;
Figure 10 illustrates a k-plane curve associated with an N-Leafed Rose scan according to the present invention; and Figure 11 illustrates the gradient fields required to produce the scan of l0 Figure 10.
DETAILED DESCRIPTION
In the following is presented a summary of the basic principles of magnetic resonance imaging along with a detailed description of the mathematical analysis illustrating how the present invention may be implemented and brought into practice.

Generation of the Signal The signal that is eventually used to reconstruct spatial distribution of the X-Y components of the spin density is essentially caused by the precession of the spins as their magnetic lines cut the wires of the receiver coil. Following is a detailed description of the generation of the signal as described by Hinshaw and Lent6. A schematic representation showing the arrangement of magnets in a conventional MRI machine is shown in Figure 1, while a simplified block diagram of the receiver system is illustrated in Figure 2.
In Figure 1, the MRI machine 10 is used for imaging the head of a subject 12. The machine includes four main field coils 14 that surround the object to be imaged. Four X-gradient field coils 16 are located with two at the front and two at the back of the machine, with the coils at each end vertically spaced. The four Y-gradient field coils 18 are similarly arranged but are horizontally spaced.
The Z-gradient field coils 20 are annular and are located coaxially with the main field coils, at either end of the main field coils.
Inside the main field coils is the radio frequency receiver coil 22 which surrounds the object to be imaged. The radio frequency transmitter is similarly located.

With reference to Figure 2, the signal from the receiver coil 22 is delivered through a matching network 24 and a pre-amplifier 26 to two phase sensitive detectors 28 and 30. One detector also receives a 90° phase shifted reference signal from a phase shifting circuit 32. The detectors 28 and 30 deliver signals to respective low-pass filters 34 and 36 which pass the signals SA(t) and SB(t) respectively.
In the receiver coil (the first block of the diagram in Figure 2) the spins generate a voltage via Faraday's law that is given by h(t) _ - dfff 3M(t~h)'Bc(h) dP 1.
JJJR
dt where ~ denotes the usual scalar product on R3. The function B~ is the receiver coil sensitivity function. To simplify the discussion, it is reasonable to assume that the coil is uniformly sensitive to the nuclear magnetization over a finite volume of space. That is, it is taken Bc (p) _ (~ + ~ ~,'s Ch) 2.
where xs is the characteristic function of the set S (which is equal to 1 on S
and 0 elsewhere) over which the coil is sensitive, a and b are constants and I and ~
are unit vectors in the so called read and phase directions, which we denote as the x and y directions. Also assumed here is that only the spins in a selected slice have been excited by the previous application of a radio frequency pulse via the radio frequency transmitter. Thus the x and y directions represent Cartesian 5 coordinates on that slice.
Using Equation 2., Equation 1. reduces to v(t) - a jjfs (aMX (t, P) + bMY (t, h)) dp~ 3.
dt Note that only the components of the nuclear magnetization density in the X Y
plane induce voltage in the receiver coil. Because of this fact and because a Fourier to integral will follow later in the analysis, it has been the tradition of MRI
researchers to express the X and Y components of the magnetization density together in one complex number, M(t, p), as follows MCt~P) =MxCt~h) + iMY(t~p)~ 4.
Knowledge of the form of the solution of Bloch's equationsl'2'n is now used.
The solution at a point p consists of a vector M(t, p) precessing about the Z
axis with a frequency of - yBZ (t, p) (clockwise), where y is the gyromagnetic constant and BZ is the Z component of the applied magnetic field. The X Y component of M

will decay exponentially with time constant TZ ( p) and the Z component will increase exponentially to the equilibrium density Mo (p) with time constant T, (p) .
Considering only the X Y component of M(t, p) there is a vector of magnitude A~p~ a ~t ''~TZ~p~ at an angle to the X axis of e(p)-y~BZ(~, p) dz radians where tL
denotes the time at which measurement begins. The quantity A( p) is the initial length of the transverse magnetization vector and e(p) is the initial angular position. The applied magnetic field during measurement is given by BZ (t, p) = Bo + G(t) ~ p 5.
where Bo is the magnetic induction produced by the main field coils and G(t) ~ p = G,~ (t)X I + Gy (t)Y J + GZ (t)Z K 6.
is the magnetic induction provided by the gradient coils. Equation 6. can be rewritten in terms of the x, y and z read, phase and slice coordinates as G(t) ~ p = GX (t)x i + Gy (t) y j + GZ (t)z k 7.
where GZ is equal to zero during measurement. It will become convenient to introduce the function r K(t) = 2~ f G(i) dz 8.
tL
where y is the gyromagnetic constant of the nuclear species being imaged, so that the X Ycomponent of M(t, p) can be expressed as MCt~ p) = A(1~) a ~t-tL)TZcp) e~e(p) e-2~i(fo~t-t~)+K(t)p) 9.
where we have introduced the Larmor frequency fo = wo/(2 ~) (radians/sec) where ~o = Y Bo (cycles/sec). We are interested in recovering A(p), the magnitude of the quantity -(t-t~)Tz~P) 18~P) P Ct~ p) = A(p) a a . I o.
In light of the definition of M(t, p) given in Equation 4., we can rewrite Equation 1.
to as V (t) dt f f ~s Re[(a tb) tYI (t, p)] ~. 11.
Defining the constants c and ~ from a = c cos(~) and b = -c sin() and using them in Equation 11. yields V(t) dt f f js Re[c e'~ M(t, p)] dp. 12.

After substituting the expression for M given by Equation 9. into Equation 12.
we end up with an expression that gives the voltage in the coil as V (t~ - - d ~~~ A(P) e-''-'L'''Zc~> ~ cos[ø - r~o (t - t~ ) - 2~ K(t) ~
p+9(P)~ ~'P. 13.
Moving the differentation operator past the integral gives y(t)=c~''o~~~ A(P)e~'-'L>~TZcp~ (cos~~-~o(t-t~)-2~cK(t)~p+e(P)~~ T
0 2 (P) y G(t) _ p 14.
- sine - ~o (t - t~ ) = 2~ K(t) ~ p + 8 ( p)~ 1 +

Based on physical magnitudes, the terms 1 /(coo TZ ( p)) and (y G(t) ~ p) l ~o are known to be small and can be neglected. This reduces Equation 14. to y(t)-~'~~~sA(P)e ~'-''»rzcn~~os~,-ay (t-tL)-2~K(t)~p+B(P)~dP 15.
where we have introduced cl = - c wo and ~1 = ~ - (~/2). And we have arrived at our working expression for the voltage generated in the receiver coil by the nuclear magnetization.
The next stages of processing are in the Matching Network and the Preamplifier (see Figure 2). These affect the signal by modifying the constants cl and ~ 1 of Equation 15. to c2 and ~2 respectively. (The Preamplifier causes the magnitude of c2 to be much larger than the magnitude of cl.) For completeness, we write the signal as it emerges from the Preamplifier as v(t) - c2 JJJs A(P) a ~' 'L~'Tzcp~ cos [~2 - ~o (t - tL ) - 2~ K(t) ~ 9 (p)~
dp~ 16.
The signal now enters the first stage of processing by entering the Quadrature Phase Detector. Here the signal is split into two channels and each channel is mixed (multiplied) with a sinusoidal signal. Without loss of generality, it can be assumed that channel A is mixed with rl(t)=acos[w(t-tL)] and that channel B is mixed with r2(t) = a sin[cu (t-tL)]. Mixing then produces the signals sA(t) _ rl(t) v(t) and sB(t) = r2(t) v(t) or using Equation 16., sA (t) _ c2 f ~~s A(P)e ~' ''~'TZcn> a cos[cc~(t-t~)~cos[~z-~o(t-tr.)-2~K(t)~p+6(P)]17.
And sB (t) -c2 JJ Js A(P) a ~' '' ~'Tz cp> a sin[~ (t - t~ )~ cos~~z - ~o(t - tr. ) - 2n K(t) ~ p + 6(P)] 18.
Equations 17. and 18. can be expanded to SA(t)-~za~~~s A(P)e ~' ''~'TZcn> ~cos~(~+cooXt-tr.)Wz-e(P)+2~cK(t)~p~
19.
+2cos~(co-cooxt-tL)+~z+8(P)-2~K(t)~p~
and SB (t) - cz a ~~~ A(P) a ~' ''''Tz cn> ~ sin~(~ + ~o xt - tc ) - ~z - a (P) +
2~ K(t) ~ p 20.
+ ~sin~(ar-c~oXt-t~)+~z +9(p)-2~cK(t)~ p, ~.
Next the low pass filters will remove the part of the signal at the frequency co + coo 5 to leave SA(t)-c3JJJsA(P) a ~' ''>'TZcn~ cos~(to-~o)(t-tc)+~2 +e(P)-2~cK(t)~ p~a~2l.
and SB (t) - c3 ~~~ A(P) a ~' '' ~'TZ cps sin~(cc~ - ~o )(t - tL ) + ~z + e(P) -2~ K(t) ~ p~ ~ 22.
where c3 = (1/2) a c 2. Finally the signals on the two channels are viewed as one 10 complex valued signal via S(t) = SA(t) + i SB(t) 23.
which is expressed as ) KJJJS '4~t~ a (t t~)T2(p) eiB(P) ei(~-~o)(t-tc) e-2~rK'(t)~p dp 24.
where K = c 3 a '~z .
If the frequency of the reference signal, cu, is taken to be equal to the Larmor frequency coo of the main magnetic field Bo (to produce a resonance signal), and the definition of p given by Equation 10., the expression for the signal, Equation 24., becomes S(t) - K f j js P (t, p) e-2~ ~K(t).P dp 25.
Note that this signal exists only after the time tL and lasts for some finite time interval of length TM. We now assume that TM is small in comparison to TZ ( p) . So the X Y component of the magnetization density can be taken to be a function of space only:
P ~h) = A~p) e~e(p) . 26.
With this assumption, the expression for the signal is rewritten as S(t) - KJJJS P ~p) e-2~c~K(r).p ~ 27.
Equation 27. now defines the MRI signal from which p ( p) is to be recovered.
Usually the quantity A(p) is of interest but in some flow experiments the phase 6(p) is of interest as well.
Finally, we want to take into account that only a thin slice of spin has been excited and to integrate the information in the thickness direction. This leads to a two dimensional set A that we will call the active set and to the consideration of p as a point in R2 instead of R3. Thus Equation 27. reduces to S(t) JJA P ( j]) a 2niK(t).P dp 28.
where K has been put equal to 1 since it does not significantly affect the reconstruction process of recovering p from S. Note that it may be desired not to use some of the information at the beginning of the measurement. We will call this time the null measurement time and will denote the time at which the useful signal begins as tl where tl >_ tL. The total time during which useful signal is collected is therefore T = tM - (tl - tL) = tF - tl where tF denotes the final measurement time. With these definitions, we can rewrite Equation 8. as t, t K(t) = Y ~ G(i) di + Y ~ G(i) di 29.
2~ tL 2~ t, when t >_ tl.
Current Reconstruction Practice In the following, the current techniques used to generate MRI signals and the discrete Fourier transform methods that are used to effect the reconstruction will be reviewed.
Fourier Imaging and Spin Warp Imaging These two methods produce the same signal but differ in the kind of gradient field used during the null measurement time. As discussed above, the null measurement time is a time in which the spins are exposed to a gradient field for the purpose of measurement but no signal is recorded. It is a time before tl that is not used for spin preparation. In what follows we will put tl = -Tl2.
In Fourier Imaging, the y gradient field is held at a constant value 2 ~
ry l y for a variable amount of time ~tq~ during the null measurement time if q >_ 0, and at a constant value of -2 n by / y if q < 0. Here t9 = qtl, -Q<-q<-Q-1 30.

with tl being a constant positive value (Q is even). After the y gradient field is switched off, the x gradient field is then switched on at a constant value of -2 ~ rx / y for a total time of T/2.
In Spin Warp Imaging, the y gradient field is held at a constant value 2 ~c Gy(q)ly (which is different for each q) for a fixed total time tl during the null measurement time. Here Gy~q)= qry~ -Q~q~Q-1 31.

to with ry being a constant value. At the same time, the x gradient field is also switched on at a level of -(2 ~ rx)/(tl y).
In both cases, the x gradient field is held constant at 2 ~ hx / y when t >_ t I
Therefore, with K = Kxi + Ky j Equation 29. gives:
In the case of Fourier Imaging:

0 for -~tq ~-T St<-T
KX(t,q)= -rx (t+T) for -T St<--2 I~'Xt for - 2 <_ t <_ 2 sign(tq ) rv (t+ ~ t9 ~ +T ) for - ~ tg ~ -T <- t <_ -T
32.
Ky(t,q) _ rytq for -T < t <
In the case of Spin Warp Imaging:
rxT t+T +t, for -t~ -T <t<-T
2t, 2 2 2 Kx (t, q) _ rxt for -T <t<T

Gv(q) t+~+t, for -t,-2 <-t<-2 33.
5 Kv(t,q)=
G,, (q) t, for - T < t < T

In both cases the k-plane trajectories are as shown in Figure 3 and both cases produce the signals -2niCrXtx+~ryt~ yl S(t, q) - ~ f A p (x, y) a J dx dy 34.
for each q E { -Ql2, ..., (Ql2)-1 } and for -Tl2 <_ t <_ Tl2. Current practice prefers the Spin Warp technique because the T2 relaxation effects are minimized by virtue of the shorter null measurement time.
We see that each signal given by Equation 34. is the Fourier transform of p(x,y) along the curve (Kx(t), Ky(t)) _ (rx t, q ry tllQ) in the k-plane since p has support in A. Specifically, the curve is a line parallel to the kx axis as shown in Figure 3. If we were to discretize (sample) the signal N times then, taking all Q
signals into account, we would get a rectangular grid of points on the k-plane; just what we need in order to use an Inverse Discrete Fourier Transform (IDFT) algorithm such as the Fast Fourier Transform (FFT) to recover an approximation of p. The definition of the Discrete Fourier Transform (DFT) and the IDFT can be found, for example, in Brigham3.
Let us look a little more closely at the IDFT reconstruction process. In terms of its Fourier transform, p is given by P ~x~ Y) _ ~ ,~ P (kx, ky ) e2~r i(kXx+kyY) dk dk x y 35.
where p is the Fourier transform of p. The signal, Equation 34., can be written as 36.
S t~ q~ - P rxt~ q rytl To get a Riemann sum approximation for Equation 35. we sample the signal at t = 1 Ot for 1 E { -Nl2, ..., (Nl2)-1 } where Ot = TlN. (Notice that we are ignoring the last data point if tF = Tl2 but not if t F = Tl2 - TlN.) This corresponds to choosing k x =1 Okx, where Okx = rxOt, for each a . More explicitly, Llkx = l x T . 37.
N
For the q'h sample, ky = q Oky where t-' 38.
to Oky = ry Q .
Therefore S( ~Ot, q) = P (~Okx, qOky ) 39.
and if we only use the information contained in the signal, we get the following Riemann sum approximation for Equation 35.

_ p~~~ ~~_ 2 n i ( e~kxx + q~k,, ) ~l' ~~~
P Cx~ Y) P e~x ~ q~y ) a ~ L1I(.x !~y 4~.
q=-Q l 2 ~---N l 2 or Ql2-1 Nl2-1 2~i(PrXTx+gryt~y) P Cx~ Y) _ ~ ~ S ~~~~ ~') a N Q
x N y - 41.
q=-Ql2 2=-Nl2 If we subject the N samples of the qth signal to an IDFT (modified by multiplying by a suitable unitary number), we get the N periodic function on the integers given by Nl2-1 6n q 1 S(,e~t~C~) e2~inPlN for n E Z 42.
N ~=-Nl2 To get the approximation to the two-dimensional transform we then subject the nth sequence of a".9 values to another IDFT (again suitably modified via multiplying by a unitary number) to obtain 1 Q/2-~
6 n'm - - ~ 6 n'q a 2n imq l Q for m E Z 43.
Q q=-Ql2 or substituting Equation 42. into Equation 43:
1 Q/2-1 N/2-1 2~i~~+
~n,m ~ S ~Ot~ ~~ a . 44.
q=-Ql2 2=-Nl2 Comparing this to Equation 41. gives n m rxT ryt'6 n,m ~ P I-' T ' h t 45.
x y 1 And we have recovered an approximate reconstruction of p. Note the dependence of the pixel resolution on the parameters rXT and rY tl.
The behaviour of this approximation is well understood (see for example Brigham3). We must choose n and m such that n m EA
l0 . 46.
rx T ' ry t, As a function on Z2, an,m is N x Q periodic and if the number of points in the set n m n m EA 47 rxT ' r,,t, rxT ' rv t, is larger than NQ then abasing will occur. (Note that if there are less than NQ
points that satisfy Equation 46. then there is no harm in redefining the set A
to be large enough to include exactly NQ points of the form ( n/(rXT), m/(rytl) ). ) A
5 "rippling effect" (equivalent to a convolution of p with sinc functions) caused by the finiteness of T and Q will also occur.
Zeugmatography The word "Zeugmatography" was coined by Lauterbur7 to describe the technique first used to produce NMR images in 1973. According to Lauterbur, the 10 term comes from the Greek ~svy~a meaning "that which is used for joining", refering to the joining of the Radio Frequency (RF) irradiation technique used in the preparation phase of the experiment and the gradient fields technique used in the measurement phase. As such, Lauterbur had intended the word to apply to all NMR imaging but it has stuck only with the technique described in this section.
15 The technique that Lauterbur invented and used to generate the first published NMR image was motivated by the "projection reconstruction"
techniques then (and now) commonly used to process the signals obtained from X-Ray Computer Assisted Tomography machines ("CAT scan" machines). The process involves scanning the k-plane along radial lines. This process gives signals that are very closely related to the Radon Transform (or "projection") of p as we will see below. Once that connection is made, projection reconstruction, an approximation of the Inverse Radon Transform can be used to reconstruct the image. Let us look at the details of this process.
Zeugmatography is a multiple signal process in which Q signals are collected with the following constant gradient fields applied:
Gx (t) _ 2n rx (q) Y
GY (t) _ 2~c rY (q)~ 1 ~ q ~ Q
Y
where rx(q) and rY(q) are chosen such that r = rx2 + ryz is constant and A(q) = tari 1(rx(q) l ry(q)) satisfies A(q) = 2 ~ q l Q. Let ~p : R2\{0{ ~ R x [0, 2 ~]
be the rectangular to polar coordinate transformation with ~p (kX,ky) _ (k,8) where k and 8 are the usual radius and angle of polar coordinates. Then we can express the resulting k-plane curves as k(t,q)=rt, o<_t<_T
e(t,R')= q 2T~~ l~R'~Q. 4s.
Q
The associated k-plane trajectories are depicted in Figure 4. Thus the signals generated are s(t, q) = P rr, Q
49.
where p is expressed in polar coordinates.
To see how this kind of data are used to reconstruct p, we will take a look at the limiting case of having information about p on every line through the origin. If we could have tl = -oo and tF = oo then the signals would give p (k,B) for all k E R.
Introduce polar coordinates (r, a) on the image space where r and a represent the radial and angular parts related to the coordinates x and y in the usual way.
Thus we would recover p via the inverse Fourier transform written in polar coordinates:
P (Y, a ~ k 9 ) a 2~ ikr ~cos B cos a + sin a sin a ~ ~ k ~ dk d 8 50.

Equation 50 is then used as the basis for an approximation using the limited data actually given. But before we present that reconstruction, we will expose the relationship that the signal (the ideal, infinitely long one) has to the Radon transform of p. The Radon transform of p (or projection of p, as some call it) in the direction 8 is given by p (r,8) _ ~ p (s(8), r(9)) ds 51.
where r and s are related to x and y via the rotation r rose sing x 52.
s - sine core y Using the s(9) and r(8) coordinates, we can write the Fourier transform of p as p (k, a ) _ ~' f p (s (e ), r(e » e-2'~ ~~ dr ds 53 .
or P (k~e) _ ~ p (r,9) e-2~~~ dr s4.
so, with respect to the k variable, p(k,9) is the one-dimensional Fourier transform of the Radon transform of p in the direction 8. Thus one would recover this Radon transform by taking the inverse Fourier transform of the signal. Standard algorithms for computing the inverse Radon transform (such as filtered back projection, as discussed by Louis and Natterer9 or by Lewittg for example) could then be used for reconstruction. However, it is more efficient to use Equation 50.
directly as we will now show.
Assuming that we are given N samples of each signal as given by Equation 49., a straight forward numerical approximation of Equation 50. can be written (using Equation 49.) as 2,~ T2r Q N nT 2~lnr~r~cos~2Qq lcosa + sinC2Qq lsinaJ
PO~a) = 2 ~~n S 'q a J J 55.
QN
where samples of each signal are spaced TlN apart. In practice (in order to make use of the FFT) the results of Equation 55. are usually computed on a polar grid and then linearly interpolated onto a rectangular grid of points so that the picture can be displayed on a computer screen composed of rectangularly arranged pixels.
Single FID Methods An NMR signal produced via pulsed RF is called an FID for Free Induction Decay. When one uses the multiple FID Fourier Imaging, Spin Warp Imaging and Zeugmatography methods, a time of length 3 to 5 times Tl (the "spin-lattice"

relaxation time) must be allowed to elapse between signal acquisitions. This is to allow the spins to realign with the Z axis so that the NMR experiment can begin anew. However, if one allows the GX and Gy fields to vary in time during the measurement in such a manner that a suitable area of the k-plane is traced out by 5 K(t) , then we can hope to reconstruct the image from one signal. These methods offer considerable savings in signal collection time over the multiple FID
experiments.
The signal for a single FID experiment is given by -2~Ci~Kx(t)x+KY(t)y~
S(t) = f fA p(x, y) a dx dy . s6.
10 Traditional methods of recovering an approximation of p from signals in the form of Equation 56. involve an application of the FFT in a heuristic manner using knowledge of the trace of K(t) in the k-plane. A complicated alternative to the direct FFT approach for Echo Planar type signals is given by Fiener and Lochers.
Echo Planar Imaging 15 Echo Planar Imaging uses rapidly switching gradient fields during the measurement time, resulting in an associated k-plane curve that covers enough of the k-plane during one or two passes for image reconstruction. The method was first introduced by Mansfieldl° in 1977 and since improved by many workers.
Several variations are in present use and here we will briefly review those that are discussed by Callaghan4. It should be noted that due to the high gradients and rapid switching required for this method that the associated gradient field hardware is expensive.
The first method we'll look at is called FLEET for Fast Low angle Excitation Echo planar Technique. The low angle refers to the fact that the spins are tipped less than 90° by the RF pulse during the preparation phase.
This is done because a second RF pulse is applied before a second sweep of the k-plane (i.e.
to some preparation is done between the two measurement phases). This double preparation allows one to acquire the second signal immediately without having to wait 3 to 5 times the Tl time.
Side-stepping consideration of the preparation involved, let's take a look at the gradient fields used in FLEET. For the first sweep:
2~ rx fog (P -1)T <- t < ~T , ~ E {1,4,5,8,9,....,452}
y 4S2 4S2 57.
2~ rx for ~~ -1)T < t < ~T , ~ E {2,3,6,7,....,452 - 2,452 - l~
y 4SZ 4SZ

2~ r Gy (t) _ - '' for 0 <- t <- T 58.
Y
where Gy is held constant and Gx is switched to produce a square wave. The number of square wave cycles is S2. When Equations 57. and 58. are integrated as per Equation 29. we obtain:
rT t - (m 1)T for (m 2)T <- t < mT , m E {5,9,13,...,452-3}
4S2 4S2 4Sz and 0<-t< 4~,m=1 KX (t) _ and ~m4~~T <- t < T, m = 4S2 59.
-rx t- (m 1)T for (m 2)T <_ t < mT , m E {3,7,11,...,452-1}

Ky (t) = rvt for o _< t <- T 60.
where Kx has become a triangular wave of S2 cycles.
For the second sweep, the sign of rx is changed. This results in two k-plane curves as shown as dashed and solid in Figure 5. Comparing Figure 5 to Figure 3, we see that the gradient field strength, rx needs to be at least an order of magnitude higher for Echo Planar imaging than it does for Spin Warp imaging.
Reconstruction then proceeds by pairing the rightward kx positive branches of one sweep with the leftward kx negative branches of the other and by pairing the leftward kX positive branches of one sweep with the rightward kx negative branches of the other sweep. This results in nearly horizontal data lines in k-space and samples along these spliced lines are subjected to an FFT as if the data did come from horizontal lines. Also since only data from positive kY values are collected, it must be assumed in this method that p is real so that complex conjugates can be used to form negative ky data. The error introduced by the actual slight angle of the data lines to the kx axis is difficult to quantify and we will not attempt to do so here.
More recently, according to Callaghan4, Mansfield and coworkers have refined their technique with methods that place even more demands on the production and switching of the gradient fields. These methods are known as the BEST, for Blipped Echo planar Single-pulse Technique, and MBEST, for Modulus BEST, methods. The x and y gradient field configurations for BEST are shown in Figure 6. The resulting k-plane trajectories for BEST are shown in Figure 7.
The x and y gradient field configurations for MBEST are shown in Figure 8 where a null measurement time, the time of gradient manipulation of the spins before signal measurement is shown. The resulting k-plane trajectories for MBEST are shown in Figure 9.

Reconstruction using sampled data from BEST or MBEST
signals using the FFT is a straightforward exercise, with the quality of reconstruction being similar to the FFT reconstruction computed from Spin Warp signals.
The Natural K-Plane Coordinate Reconstruction Method (central to the present invention) The basic problem in 2 dimensional MRI is trying to produce a 2 dimensional picture from a 1-dimension signal. It is somewhat remarkable that it can be done at all! Part of the motivation for the invention of the natural k-plane coordinate reconstruction method comes from a desire to use all of the information contained in the signal as opposed to only using samples from the signal as is done in conventional MRI reconstruction. This, in turn, has led to a reconstruction method that uses the given k-plane samples as opposed to those interpolated onto a grid.
If we knew the value of p (k) for all values of k = (kx, ky ) then we could recover p exactly via the inverse Fourier transform:

2~c i(xkX + yky) P Cx~ Y) _ ~ ~ P (kx' ky ) a dkx dk y 61.
If it were known before hand that ,o had support on some set B in the k-plane then Equation 61. would reduce to ., 2~ i (xkx + yky ) P Cx~ Y) = B P Ckx ~ ky ) a dkx dky . 62.
5 Since we know that the unknown p in MRI has compact support in a picture rectangle A, we know that p cannot have support on any compact set B in the k-plane. Nevertheless, if we computed the RHS of Equation 62. we would end up with a function which we will call PB p which is a band-limited version of p.
Thus we have .. 2~ci(xkx + yky) to I'B P Cx~ Y) _ ~~B P ~kx' ky ) a dkx dky . 63.
We will let the set B be defined by a swath of width 2w in the k-plane carved out by the given curve K(~) . The curve is already parametrized by t, so we let t be one coordinate. The other coordinate is then defined to be in a transverse direction to the curve in some convenient (but arbitrary) manner. This other coordinate direction will be labeled (3 and, for convenience, we will label the points on the given curve with [3 = 0. (We therefore have (t, (3) E [tl, tF] x [-w, w].) The relationship between the (kx, kY) and the (t, (3) coordinate systems is given by a coordinate transformation:
v(~.R)=(k.(t,P)~k.(t.R)) sa.
With the coordinate transformation of Equation 64. in hand, we can rewrite Equation 63. in terms of t and (3 as ~'B P (x~Y) _ JJB P ~kx ~t~ ~)~ky (t, /3)) e2~~(xkx~r,a~+YkVc~,~» ~ J~ (t~ ~3) ~ dt dpi 65.
where .l~ = DAP is the Jacobian matrix of the transformation and ~ .I~ ~ is its determinant, the Jacobian of the transformation.
Equation 65. also allows for the introduction of a point-spread function. This can be seen by substituting the definition of p into Equation 65.:
w tF oo ao I'BP Cx~Y) - ~ f f f P (a, b) 2 2~l(kX(t,~)a + ky(t,/3)b) da db -w t, ~o-oo 66.
X e2~i(xkx(t>a)+Yky(t>~)) ~ Jy (t' dt d~ = P'~' 'YB
~P

where ~ Cab) = w tFe2~ci(kx(t~Ij)a+ky(t~~)b) ~ J~(t~l~) dt dpi 6 l _w t, The full information needed for Equation 65. is available only if one could run a continuous set of MRI experiments corresponding to the variation of the parameter (3. But we are interested in constructing an approximation for p from a single experiment when ~3 = 0 where we only know p (t,0) (using the (t, (3) coordinates) for tl <_ t <_ tF. To motivate the proposed reconstruction method let's make an assumption about the nature of p in B. Assume it has the tensor product form:
1 o PCt~ I~ = s'(t) T'~'(L~ 69.
for some assumed function W. The function S is, of course, meant to represent the signal emerging from the MRI machine. If p satisfies Equation 69. then 1'BP Cx~ Y) - ~~ Sit) W C~) 22~i(xlrx(t~~) +Ykyt~~)) ~ J~ (t~ ~) I dt d/3 .
~o.
This suggests that we could define a map, PB, W, on suitable functions p as follows:

~'B,wP ~x~Y) _ JJB p ~t'~) W ~~) e2~~~XxX~r,l~> +vkycl~l~» I J~ (t~ ~) ~ dt d~3 . 71.
If the width of the swath defined by B (the width 2w of the interval in which the values of (3 lie) is in some sense small then one would expect that the value of PB, ~.p(x,y) does not depend strongly on the choice of W. In the extreme, if we do not want to make any assumptions at all about the nature of p in the (3 direction then we could choose W((3) = 8((3). This leads to the map 1'K P (x~ Y) = I'a,s P ~x~ Y) _ J Ja P (t'~) W O) e2~t t (xkx (t.~ ) +Yky O,~
)) ~ J~ ~t~ ~) ~ d t d/3 IF
= JP~t~~) e2~~~x~~~>+vxvct» ~ J~(t~p) ~ dt . 72.
r, Equation 72. defines the reconstruction we are after, it uses coordinates that are more naturally suited to the given k-plane curve than the usual (kx, kv) rectilinear coordinates. Since S(t) = p (t,0) in an MRI experiment, the continuous version of the natural k-plane coordinate reconstruction method is to compute as the reconstruction:
tF
R~x~Y) = f S(t) e2~i(xKX(t)+yKy(t)) ~ J~(t~p) ~ dt, 73.
t, Note that while Equation 73. uses no information about the nature of p in the ~3 direction, the choice of (3 itself affects the final reconstruction through .I~ .
The natural k-plane coordinate reconstruction also leads to a point spread function:
tp op o0 pgP(x~v) f f f p(a,b)e 2ni(aKX(t)+bKY(t)) dadb e2"i(xKx(t)+yKy(t)) ~ J~(t'0) dt tt -ao--m co ao tF
= f f p(a,b) f a 2~i(Kx(t)(a-x)+Kv(t)(b-y)) I J~(t,0) ~ dt dadb t, 74.
tF
= J J p(a,b) f e2"i(Ks(t)(x-a)+xy(t)(y-b)) ~ J~(t,0) ~ dt dadb rt = P*~K(x~Y) where tF
~ (x y) _ e2~a(Kx(t)x+Ky(t)y) ~ J (t 0) ~ dt ' ,~ ~ ' 75.
t, is the point-spread function. Thus we see that if ~ has a large peak at the origin (i.e.
is close somehow to the Dirac b function at the origin) then 1'K p will give a reasonable representation of p. This will occur when a sufficient region of the k-plane is covered by the curve K .

Discretization In order to compute the reconstruction tF
RS(p) = f S (t) e2~~(p.K(t)) I J~ (t~p) ~ dt 76.
t, exactly we would need an analogue computer. In this era of digital computers, the 5 signal emerging from the MRI machine has made a trip through an analogue-to-digital (A/D) converter. This gives discrete samples S(tj), j E f 1,2, ... , N}, of the signal. As a result, instead of computing the integral of Equation 76., it is natural to compute a Riemann sum approximation of this. This leads to a new reconstruction operator RN, that represents the discrete version of the natural k-plane coordinate to reconstruction method, given by N
RNS(p) - ~ S(t j ) e2~ci(p.K(t~ )) ~ J~ (t J'O) ~ ~t 77.
j=1 where Ot = (tF - tl)lN. (We are assuming equally spaced samples here since this is what most A/D converters produce and because the theory is the same for unequally spaced samples; the formula of Equation 77. would just be more messy.) 15 The behaviour of the map RN is slightly different from that of R even though limN ~ ~ RN S = RS for continuous functions S by the Riemann Integration Theorem, when ~ J~ ~ and K are continuous.
Consider the ma.p pN,K given by N ~ \K\t ' )) e27c i( p~K(tJ )) ~ J~ !t ' ~~ ~.
N,KP Cp) ~P J tP J ~ . 78.
j=1 Equation 78. can be rearranged to give a point-spread function:
N °° °°
I'N xP (P) - ~ ~ j p (a) a 2~ri(a~K(t~ )) dll e2~i(p.K(t; )) ~ J~ It J'O) J=1 ~o -oo N -2~ci(a'~K(t~)) 0) ~ ~ C~CZ 7 - ~ ~ P Via) ~e ~ J~ ~t~
~o -~ ~=1 In other words, - P ~N~K 80.
where N
-2~ i(a~K(t~ )) to ~N,K gyp) _ ~ a ~ J~ ~tj ~~) ~ ~t . 81.
j=1 It was mentioned above that the A/D converter provided samples of the signal at times t~. Actually this is only an approximation. The A/D
converter gives the integral of the signal over a small sampling time interval. If this time interval is small and the signal is continuous (which it is in our model by the Riemann Lebesgue Lemma) then this integral, divided by the length of the time interval, provides a very good approximation to the value of the signal at the middle of the sampling time interval. However, this characteristic of the A/D
converter might affect the kind of approach one might take if one were to find a more intellegent way of using the signal "samples" than in a Riemann sum (i.e.
if one were to try and use another numerical quadrature scheme for the reconstruction map).
Taking the samples as being exactly the value of the signal at time t~, the natural k-plane coordinate reconstruction method gives as the reconstruction of p the function PN,K p where I'N,K is an integral operator with the difference kernel N
2~ i(K(t~ )~(p-a)) is ~~a~ p) - ~ a ~ J~ ~t~ ~~) ~ ~t 82.
where Ot = (tF - tl) l N. That is, with the kernel of Equation 82, I'N,K P (p) = J JA P (a ) n(a, p) da . 83.
Modelling the A/D converter more exactly by taking the j'h sample as oT
r;+-z S; = S (t) dt 84.
0 i _ez t' - 2 (instead of S~ = S(t J) as we have been doing) where, usually, Di « 0t, and using Equation 77. to reconstruct leads to a modified integral operator, ~N~, whose kernel, K, is slightly different from the difference kernel of Equation 82.
The reconstruction then is really ~N gP CP) _' JJA p Via) K Via, P) da . 86.
Since the signal is continuous, we expect that ~NKP~pN~P as oz ~ o.
l0 The Effect of an Out-of Phase Detector While, usually, the unknown density is known to be real, the detection process may introduce a phase factor ~2 into the signal. In that case the signal will not be S(t) = p (K(t)) but rather S(t) =e'~ p(K(t)) where ~2 may or may not be known. If c~2 is known then we can simply multiply the signal by a -'~z Otherwise we will reconstruct ei~ p * ~ where ~ is the point-spread function associated with the reconstruction process. In the case that ~2 is zero we would reconstruct (p * ~R ) + i(p * ~ I ) where ~ R + i~ 1= c~. When c~ 2 is unknown then it is best to compute the modulus of the reconstruction, when p is real , giving pPl= (P*~R)2 +~P*~1)2 . 87.
since I e'~z p * ~ I= I P * ~ I . In the case that ~ is real and p is positive, Equation 87.
will give exactly p * ~. Thus we can see advantages in choosing k-plane curves that lead to real point spread functions.
N-Leafed Rose Scans In terms of polar coordinates, the equation for an n-leafed rose is given by r = A cos( n 8 ) . 88.
A plot of Equation 88. showing the associated n-leafed rose k-space trajectory, is given in Figure 10. The rose pattern actually has n leaves if n is odd and 2n leaves if n is even. Polar coordinates on the k-plane are related to the Cartesian coordinates by kx = r cos(9) 89.
ky = r sin(6) . 90.
So if we take time t = 8, the substitution of Equation 88. into Equations 89.
and 90.
leads to the k-plane curve given by KX (t) = A cos( nt) cos( t) 91.
Ky (t) = A cos(nt) sin(t) , 92.
The gradient fields required to generate this curve are given by Gx (t) = dKx (t) _ -A~cos(nt) sin(t) + n sin(nt) cos(t)] 93.
2~ dt G y (t) = dd '' (t) _ -A ~cos( nt) cos( t) + n sin( nt) sin( t)]. 94.
2~
1 o A plot of these gradient field functions is given in Figure 11 with the top graph showing Gx and the bottom graph showing Gy. To get a complete rose pattern on the k-plane in one experiment, we need to take t E [tl, tF] _ [0, 2~c]. In an actual MRI experiment, this time interval would be scaled down to a time interval that is significantly shorter than the TZ relaxation time. Also, the time interval may be 15 shifted so that the curve begins at the origin of the k-plane.

Since the k-plane curve moves largely in the radial direction, much like the Zeugmatography scan, the natural coordinate to use for the (3 direction is angle. Then ~pX (t, ~3) _ cos(/3) - sin(~3) Kx (t) _ A cos(nt) cos(/3 + t) spy (t, ~3) sin(/3) cos(~(i) Ky (t) A cos(nt) sin(~3 + t) J ' 95.
The partial derivatives of ~p in that case are then a~x (t, ~3 ) _ -A n sin(nt) cos(/3 + t) - A cos(nt) sin(~3 + t) 96.
at a~'' (t, ~i ) = A cos( nt) cos( ~3 + t) 9'7, a~
( -A cos nt sin + t - 98.
a~
a~'' (t, /3) _ -A n sin(nt) sin( f3 + t) + Acos(nt) cos(~3 + t) 99.
at which leads to the Jacobian ~~ (t~ ~) I= a~x (t~ ~) a~Y (t~ ~) - a~x (t~ ~) any (t~ ~) at a~ a~ at . loo.
= I AZ n cos(nt) sin(nt)) The reconstruction map is therefore given by 2n RS(x,.y) _ ~,5,(t)e2~iAcos(nt)(xcos(t)+ysin(t)) ~A2 y~COS(nt)Sln(nt)~ C~t 1~1.

for the analogue case and by ( >) 2~ciAcos(nt~)(xcos(t~)+ysin(t~)) z ( ~) ( ~)I 102.
R S x S t a A n cos nt . sin nt . ~t j=1 where Ot = 2~1N, in the digital case.
While one embodiment of the present invention has been described in the foregoing, it is to be understood that other embodiments are possible within the scope of the invention, especially with regards to the variety of possible k-plane curves, including multiple, interleaved curves and with regard to efficient methods l0 of computing the required sums. The invention is to be considered limited solely by the scope of the appended claims.

BIBLIOGRAPHY
1. Abragam, A., The Principles of Nuclear Magnetism, Oxford Clarendon Press, 1961.
2. Bloch, F., "Nuclear Induction", Phys. Rev. 70, 460-474, 1946.
3. Brigham, O.E., The Fast Fourier Transform, Prentice-Hall, Inc., Englewood Cliffs, N.J., 1974.
4. Callaghan, P.T., Principles of Nuclear Magnetic Resonance Microscopy, Clarendon Press, Oxford, 1991.
5. Feiner, L.F. and Locher, P.R., "On NMR Spin Imaging by Magnetic Field Modulation", Appl. Phys. 22, 257-271, 1980.
6. Hinshaw, W.A. and Lent, A.H., "An Introduction to NMR Imaging: From the Bloch Equation to the Imaging Equation", Proc. IEEE 71, 338-350, 1983.
7. Lauterbur, P.C., Nature 242, 190, 1973.
8. Lewitt, R.M., "Reconstruction Algorithms: Transform Methods", Proc. IEEE
71, 390-407, 1973.
9. Louis, A.K. and Natterer, F., "Mathematical Problems of Computerized Tomography", Proc. IEEE 71, 379- 389, 1983.
10. Mansfield, P., "Mufti-planar image formation using NMR spin echos", J.
Phys C. 10, L55, 1977.

11. Slichter, C.P., Principles of Magnetic Resonance, Harper and Row, New York, 1963.

Claims (10)

1. A method of magnetic resonance imaging an image region of an object, said method comprising:
generating a main magnetic field through the image region of the object to align magnetic dipoles of nuclei in the image region of the object with the main magnetic field;
exciting precession of the magnetic dipoles of the nuclei about the magnetic field such that each nucleus generates a component of a resonance signal;
producing magnetic field gradients across the main magnetic field so as to vary the temporal phase and frequency of the resonance signal components to encode two-dimensional spatial frequencies according to the positions of the nuclei in the image region;
varying the magnetic field gradients to generate a two-dimensional, non-linear, continuous scanning curve;
monitoring the resonance signal; and operating on the resonance signal using a transform based on natural coordinates in k-space and in the use of the corresponding Jacobian to weight the sampled resonance signal as part of the transform to convert the resonance signal into a representation of the spatial positions and densities of the imaged nuclei.

2. A method according to Claim 1 comprising cyclically varying the magnetic field gradients.

3. A method according to Claim 2 comprising varying the magnetic field gradients to generate a multi-leafed rose scanning curve.

4. A method according to Claim 2 comprising varying the magnetic field gradients sinusoidally.

5. A method according to Claim 2 comprising varying the magnetic field gradients to generate a spiral scanning curve.

6. An apparatus for magnetic resonance imaging of an image region of an object, said apparatus comprising main field generator means for generating a main magnetic field through the image region of the object for aligning magnetic dipoles of nuclei in the image region of the object with the main magnetic field;
excitation means for exciting precession of the magnetic dipoles of the nuclei about the main magnetic field such that each nucleus generates a component of a resonance signal;
gradient field generating means for producing magnetic field gradients across the main magnetic field so as to vary the temporal phase and frequency of the resonance signal components to encode two-dimensional spatial frequencies according to the positions of the nuclei in the image region;
gradient field varying means for varying the gradient magnetic fields with time to produce a two-dimensional, continuous, non-linear scanning curve;
receiving means for monitoring the resonance signal; and computer means for operating on the resonance signal using a transform based on natural coordinates in k-space and in the use of the corresponding Jacobian to weight the sampled resonance signal as part of the transform to convert the resonance signal into a representation of the spatial positions and densities of the nuclei.

7. Apparatus according to Claim 6 wherein the magnetic gradient field varying means comprise means for cyclically varying the gradient magnetic fields.

8. Apparatus according to Claim 7 wherein the gradient field varying means comprise means for varying the magnetic field gradients to generate a multi-leafed rose scanning curve.

9. Apparatus according to Claim 7 wherein the gradient field varying means comprise means for varying the magnetic field gradients sinusoidally.

10. Apparatus according to Claim 7 wherein the gradient field varying means for varying the magnetic field gradients to generate a spiral scanning curve.