the raw data for magnetic resonance imaging ( mri ) 
 is measured in 
 k - space , the 
 domain of spatial frequencies , where non - cartesian sampling schemes like spiral 
 or radial scans have received much attention .
in contrast to the use of the 
 computationally efficient fast fourier transform ( fft ) for the reconstruction 
 from cartesian grids , the more general sampling trajectories ask for the 
 so - called nonequispaced ffts . on the other hand ,
, iterative image 
 reconstruction in combination with the nonequispaced fft has been applied to 
 data on spiral 
 k - space 
 trajectories   and 
 in the presence of field inhomogeneities .
efficient reconstruction procedures for sensitivity 
 encoding with arbitrary 
 k - space 
 trajectories were proposed in .
its authors present methods that combine the gridding 
 principles with the conjugate gradient scheme , but mention the long computation 
 times due to their nonoptimised preliminary software . in this paper , we use a similar method , compare the 
 gridding approach 
 and an approach based on an implicit discretisation , where we focus for the 
 reader 's convenience on a simplified signal equation .
it turns out that the 
 iterative solution of the latter approach resembles the gridding method in its 
 first iteration .
both reconstruction problems are easily solved by our mature 
 software package 
 for the nonequispaced fast fourier transform ( nfft ) [ 8 , 9 ] . for readers not familiar 
 with the nfft
the 
 publicly available implementation easily allows for the efficient treatment of 
 large 3d data sets by the use of new sparse matrix storage techniques .
moreover , we compare 
 different density compensation weights arising from the discretisation of the 
 underlying integrals .
we present numerical results , based on the shepp - logan 
 phantom as well as on data acquired by an mr scanner .
the outline of this paper is as follows : section 2 
 gives a brief introduction to the theory of fourier transform image reconstruction 
 and unifies the two considered approaches in mri .
we suggest the conjugate 
 gradient method for the reconstruction problem and show that the solution is 
 efficiently computed by the iterative use of the nonequispaced fft . 

subsequently , section 3 presents the used 
 k - space 
 trajectories , considers sampling density compensation , and introduces the 
 tested simulated and acquired data sets .
section 4 shows our numerical tests 
 emphasising the reconstruction quality with respect to the number of iterations 
 and different density compensation weights .
moreover , we give detailed 
 information on the computational requirements of the suggested scheme . finally , 
 a short discussion of our results is contained in section 5 and an introduction 
 to the nonequispaced fft might be found in the appendix .
given a trajectory 
 k = k(t ) , the relation between the mr signal 
 s during the 
 readout and the object 
 p can be modeled 
 by the simplified signal equation 
 ( 1)s(k)=3p(r)e2irkdr . in the following ,
let the available samples in 
 k - space be 
 contained in the shifted unit cube , that is , 
 k  [ 1/2,1/2 ) , and the field of view be restricted to 
 n  [ n1/2,n1/2 )  [ n2/2,n2/2 )  [ n3/2,n3/2 ) , 
 where n = ( n1 , n2 , n3 )  2. then , the discretisation of integral ( 1 ) on equispaced 
 points leads to 
 ( 2)s(k)s(k):=rin3p(r)e2irk , 
 where in : = { n1/2 ,  , n1/2  1 }  { n2/2 ,  , n2/2  1 }  { n3/2 ,  ,
n3/2  1}. 
 thus , the unknown object 
 p is given implicitly by ( 2 ) .
the 
 discretisation of this integral leads to 
 ( 3)p(r)p(r):=j=0m1s(kj)e2irkjwj , where 
 wj are weights , 
 which compensate for local variations of the sampling density . here
the important difference between ( 2 ) and ( 3 ) is that the former 
 is discretised in the image domain with pixels on a uniform grid and hence with 
 unit weighting coefficients and the latter is an integral discretised in 
 k - space with 
 nonuniform samples and specific weights .
we reformulate problem ( 2 ) and ( 3 ) in matrix vector 
 notation and denote the vector of the given measurements by 
 s : = ( s(kj ) ) j=0 ,  , m1   , the reshaped vector of the unknown object by p : = ( p(r))rin   , 
 the density compensation matrix by 
 w : = diag(wj)j=0 ,  , m1 , and the nonequispaced fourier matrix 
 by 
 ( 4)a:=(e2irkj)j=0,  ,m1 ;  rin3,whereas 
 a denotes its 
 adjoint ( conjugate transpose ) . the gridding approximation ( 3 ) is easily computed 
 by one matrix vector multiplication 
 ( 5)p=aws.the adjoint nfft takes 
 (|in| log |in| + m ) floating point 
 operations for this task .
slightly more involved , the reconstruction problem 
 ( 2 ) is solved 
 by the method of least squares and hence consists in solving the weighted 
 normal equation of first kind 
 ( 6)awap = awsfor the unknown vector 
 p. in contrast to , we include density compensation weights also for the 
 implicit discretisation since this is more natural with respect to the 
  continuous residual  in 
 k - space and has 
 been proven to be better conditioned in . from the mathematical point of view , ( 6 ) is solved most 
 efficiently by the conjugate gradients ( cg ) method ( cf . [ 12 , page 288 ] ) .
we prefer to 
 solve ( 6 ) by 
 a factorised variant of conjugated gradients , where the two multiplications 
 with the ( adjoint ) nonequispaced fourier matrix per iteration are computed by 
 the nfft .
this scheme is denoted by cgnr , n for  normal equation ,  and r for 
  residual minimisation  ( cf . ) .
note that the cg method applied directly to 
 ( 6 ) as 
 suggested in 
 generates the same sequence of approximations in exact arithmetic , but the cgnr 
 approach is considered to be more stable with respect to round - off errors ( cf . 

, we suggest algorithm 1 . remarkably , this algorithm resembles an optimised 
 gridding solution after one iteration .
more formally , let the weighted residual 
 norm 
 ||r||w : = rwr be given . then
, 
 the first iteration obeys 
 ( 7)p1=argminp=psapw , since for the first residual 
 r1 = s  ap1 holds the 
 perpendicular condition 
 ( 8)r1wap=(aws)ppp(ap)wap(ap)wap=0.in other words , the gridding 
 solution 
 p 
 in ( 5 ) is scaled such that 
 its residual is minimised and we do not consider a gridding approach separately 
 anymore .
moreover , note that the proposed scheme should be 
 stopped as soon as the current residual 
 ||s
 apl||w drops below the 
 discretisation error or the level of noise in the measurements 
 s. for our first evaluation , we terminate the suggested 
 algorithm after a fixed number of iterations .
we are concerned with the reconstruction from data 
 acquired by an mr scanner as well as the reconstruction quality and the usage 
 of time and memory resources for simulated data .
algorithm 1 has been tested 
 with mr measurements of a physical phantom by the philips achieva 1.5 t device . 
 here , the sampling scheme consists of 
 36 equidistant 
 radial trajectories with 
 m = 7 557 120 points in 
 total , whereas the reconstructed image contains 
 256  256  36 = 2 359 296 voxels .
we use two 3d - shepp - logan phantoms of sizes 
 256  256  36 and 128  128  128 ( 2 097 152 voxels ) as shown in figure 1 .
comparison is done with respect to the number of 
 iterations , sampling schemes , and density compensation weights .
the simulated 
 mr data is computed by a 3d - nfft using the following 
 k - space 
 trajectories .
note that all but the first trajectory are of a special 2d 
  1d type , that 
 is , they consist of a stack of 
 36 equidistant 
 planes where the 
 m1 points within 
 each are distributed accordingly . in these cases , the first two coordinates of 
 a 
 k - space point 
 k = ( k1 , k2 , k3 ) form the 
 2d - point k=(k1,k2)[1/2,1/2)2 . 

( i ) ( 3d - radial ) the only 3d - k - space 
 trajectory , which is not of 2d 
  1d type , is 
 given by 
 ( 9)kp , q , r = r+12r(cospsinq , sinpsinq , cosq) , 
 where p = 2p / p , 
 p = 0 ,  , p  1 , q = (2q + 1)/2q , q = 0 ,  ,
furthermore , we restrict this set to the unit cube 
 ( 1/2 , 1/2 ) .
choosing 
 p = q = r = 160 in our 
 experiments yields a total number of 
 m = 3 398 033 points in 
 k - space.(ii ) ( radial ) popular also within computer tomography 
 is the 2d - radial trajectory 
 ( 10)kp , r=(1)r(rr0.5)(cos  pp , sin  pp)with 
 p = 0 , 
, p  1 , 
 r = 0 ,  , r  1 .
we set 
 p = 410 and 
 r = 512 yielding a 
 total number of 
 m = 36m1 = 36pr = 7 557 120 points in 
 k - space.(iii ) ( spiral ) this 2d - k - space 
 trajectory is given by one archimedean spiral , that is , 
 ( 11)kj = j2m1(cosj , sinj) , 
 where j=(8/5)j , and 
 j = 0, 
we have chosen 
 m1 = 65 536 yielding a 
 total number of 
 m = 36 m1 = 2 359 296 points in 
 k - space . 
 
( 3d - radial ) the only 3d - k - space 
 trajectory , which is not of 2d 
  1d type , is 
 given by 
 ( 9)kp , q , r = r+12r(cospsinq , sinpsinq , cosq) , 
 where p = 2p / p , 
 p = 0 ,  , p  1 , q = (2q + 1)/2q , q = 0 ,  ,
furthermore , we restrict this set to the unit cube 
 ( 1/2 , 1/2 ) . choosing 
 p = q = r = 160 in our 
 experiments yields a total number of 
 m = 3 398 033 points in 
 k - space .  
( radial ) popular also within computer tomography 
 is the 2d - radial trajectory 
 ( 10)kp , r=(1)r(rr0.5)(cos  pp , sin  pp)with 
 p = 0 ,  , p  1 , 
 r = 0 , 
we set 
 p = 410 and 
 r = 512 yielding a 
 total number of 
 m = 36m1 = 36pr = 7 557 120 points in 
 k - space .  
( spiral ) this 2d - k - space 
 trajectory is given by one archimedean spiral , that is , 
 ( 11)kj = j2m1(cosj , sinj) , 
 where j=(8/5)j , and 
 j = 0, 
we have chosen 
 m1 = 65 536 yielding a 
 total number of 
 m = 36 m1 = 2 359 296 points in 
 k - space .
note that each proposed sampling scheme obeys a nyquist rate near the origin of 
 the 
 k - space . 

nevertheless , all sampling sets violate the nyquist criteria in their periphery 
 since there exist boxes larger than the reciprocal field of view 
 1/n1  1/n2  1/n3 containing no 
 single sampling point 
 kj . in the following
, different weights 
 w in ( 6 ) are used to take 
 into account the local sampling density in 
 k - space .
we 
 propose formulations 
  with no 
 weights , that is , 
 w being the 
 identity , with analytic 
 weights in the case of the 3d - radial trajectory , given by 
 ( 12)wp , q , r = r+1rsin(2q+1)2q , 
 with 
 approximate weights computed in a very cheap way by counting the number of 
 sample points in the 
 256  256 box tesselation 
 of the 
 k - space 
 [ 1/2 , 1/2 ] , andwith weights 
 obtained as the area of the voronoi cell 
 ( 13)j={k[12,12)2:kkj2minl=0,  ,m1kkl2}around each sample point 
 kj ; see also [ 15 , 16 ] . 
   with no 
 weights ,
that is , 
 w being the 
 identity , with analytic 
 weights in the case of the 3d - radial trajectory , given by 
 ( 12)wp , q , r = r+1rsin(2q+1)2q , 
  with 
 approximate weights computed in a very cheap way by counting the number of 
 sample points in the 
 256  256 box tesselation 
 of the 
 k - space 
 [ 1/2 , 1/2 ] , and with weights 
 obtained as the area of the voronoi cell 
 ( 13)j={k[12,12)2:kkj2minl=0,  ,m1kkl2}around each sample point 
 kj ; see also [ 15 , 16 ] .
it was pointed out , for example , in , that the use of voronoi 
 weights as a measure of the local sampling density is a very reliable and 
 general technique . for the gridding approach ( 3 ) , a detailed discussion of sampling density 
 compensation for sampling sets violating the nyquist criteria
all tests were 
 implemented in matlab&c and tested on an amd athlon xp 2700 
 + , 2 gb memory , suse - linux , kernel 2.4.204gb - athlon , 
 fftw3.0.1 , and nfft2 ( kaiser - bessel window function with cut - off parameter 
 m = 6 and 
 oversampling factor 
   1.5 ) . besides a couple of representative examples ,
our web page on this project 
 http://www.tu-chemnitz.de/potts/projects/mri collects a 
 more thorough set of tests as pointed out in the subsequent examples . in 
 particular , this includes animated graphics showing the progress during the 
 iterations or slicing through the 3d data set .
first of all , we apply algorithm 1 to mr measurements 
 taken by the philips achieva 1.5 t device .
figure 3 shows the result after one 
 iteration , where we used voronoi weights for sampling density compensation .  

we compare the reconstruction quality with respect to the three different 
 k - space 
 trajectories and sampling density compensation weights . for the 2d  1d type of the 
 trajectories , we reconstruct each slice separately and compute the overall 
 reconstruction by 
 256   regular 1d - ffts 
 of length 
 n3 = 36 
 within the 
 third component .
table 1 shows the normalised root - mean - square 
 error 
 ( 14)rms(p,p):=pp2p2,where 
 p 
 is our 
 reconstruction and 
 p denotes the 
 original image of size 
 128  128  128 for the 
 3d - radial trajectory and of size 
 256  256  36 for all others .
additionally , we computed 
 30 iterations for the 3d - radial example which yields a rms of 0.3309 and 0.2104 for no weights 
 and analytic weights , respectively . a first hint on the stability of the 
 reconstruction process with respect to noisy data is given in the last line of 
 table 1 .
we show the rms for the reconstruction from simulated spiral data 
 which is perturbed by 
 20 percent 
 standard normally distributed noise .
the main purpose of the 3d - radial example is to demonstrate 
 that the nfft - based reconstruction is straightforward .
figure 4 shows the 64th slice and 
 the profile of the 64th row of this 
 slice after one , five , and ten iterations for analytic weights . the same 
 experiment using no weights
furthermore , we show the 18th slice and 
 the profile of the 128th row of this 
 slice after one and ten iterations for different weights and the 2d 
  1d 
 k - space 
 trajectory radial in figure 5 . the same experiment using the spiral trajectory 
 and an interleaved spiral trajectory
is presented on our web - page . in the last test , we measured cpu time and used memory of our algorithms with 
 respect to the size of the reconstruction problem 
 n , that is , the size of the reconstructed image is 
 n  n  36 
 and the total number of spiral 
 k - space samples is 
 36n .
we compare algorithm1exploiting the 
 2d  1d type of the 
 trajectory with algorithm 1 using the full 3d nfft .
the optimisation of the computing time is done by 
 different nfft flags which affect the last step of the nfft . in this last step , 

we compute by 
 (m ) floating point 
 operations the matrix vector product with the sparse matrix 
 b as defined in 
 the appendix , ( a.2 ) .
we propose different methods for the 
 precomputation and storage of this matrix that basically trade main memory for 
 computation speed .
we were able to store all nonzero entries of the matrix 
 b together with 
 their row and column index by the flag pre - full - psi in all 2d 
  1d tests ; see 
  for a similar 
 approach .
a lossless compressed form , flag pre - psi , is used in the 3d - tests up 
 to 
 n = 256 .
finally , we use a lossy compressed version which is 
 independent of the actual trajectory , uses a lookup table , and linear 
 interpolation for the entries of the sparse matrix ( pre - lin - psi ) in the 3d - tests 
 for 
 n = 512 , see the manual to [ 7 , 10 ] for details .
table 2 shows 
 our results for the shepp - logan phantom , the spiral trajectory , and voronoi 
 weights .
clearly , the method based on the 2d 
  1d model is 
 faster , however our algorithms can handle arbitrary scattered data in 
 d dimensions by 
 storing a heavily compressed form of the sparse matrix 
 b.
magnetic resonance signals are measured in 
 k - space , where 
 magnetic field gradients determine the specific trajectories that form the 
 sampling points .
we have shown that the reconstruction from non - cartesian grids 
 is easily computed by means of the nonequispaced fft where at least an 
 approximate compensation of the non - uniform sampling density proves necessary 
 for highly nonuniform sampling sets .
the gridding approach , based on an explicit 
 discretisation , and the  inverse model 
, based on an implicit discretisation , are unified and 
 solved by algorithm 1 . in particular , the gridding solution is 
 optimised in one iteration to minimise the residual , see ( 7 ) , and further refined in subsequent iterations .
we see 
 that already the gridding method may lead to very good results in table 1 and 
 in figure 5 .
however , using no density compensation weights gives reasonable 
 results only if the trajectory covers the 
 k - space 
 uniformly as for the spiral trajectory , see table 1  weights : none  and our 
 web - page .
gridding without density compensation gives non satisfactory results 
 for 3d - radial and radial trajectories as shown in the leftmost images of 
 figures 4 and 5 . including approximate weights
improves the reconstruction 
 quality , the best gridding results are obtained using the more expensive 
 voronoi weights .
however , substantial improvement of the reconstruction quality 
 is achieved during a small number of iterations for all used trajectories and 
 weights .
furthermore , only the iterative method in conjunction with a 
 reasonable approximation of weights gives an acceptable reconstruction from the 
 3d - radial trajectory ( cf . figure 4 and our web page ) .
finally , we have compared 
 the computation time and memory usage of the proposed algorithm for the 2d 
  1d model and a 
 full 3d model . in summary , the implicit discretisation and its 
 iterative solution generalises the gridding approach in a natural way and 
 refines the image quality . particularly , a poor gridding solution from a highly 
 nonuniform trajectory with approximated density compensation weights is 
 improved substantially by a few iterations of our scheme .
we have shown that 
 very efficient methods are available , which can be generalised to 
 reconstruction methods for sensitivity encoding as proposed in   and to reconstruction methods in presence of field 
 inhomogeneities [ 3 , 19 ] .
note , however , that the computation times reported in 
 table 2 depend on the stopping criterion .
directions for our future research 
 include the development of reliable stopping criteria for the iterative 
 reconstruction and proven convergence rates when the nyquist criteria is 
 violated .
 in magnetic resonance imaging ( mri ) , methods that use a non - cartesian grid in 
 k - space are becoming increasingly important . 
 in this paper , we use a recently proposed implicit discretisation scheme which 
 generalises the standard approach based on gridding . 
 while the latter succeeds for sufficiently uniform sampling sets and accurate 
 estimated density compensation weights , the implicit method further improves 
 the reconstruction quality when the sampling scheme or the weights are less regular . 
 
 both approaches can be solved efficiently with the nonequispaced fft . due to several new techniques for the storage of an involved sparse matrix , 
 
 our examples include also the reconstruction of a large 3d data set . 
 we present four case studies and report on efficient implementation of the related algorithms . 
1. INTRODUCTION
2. THEORY
3. METHODS
4. NUMERICAL RESULTS
5. DISCUSSION