Advances in Sensitivity Encoding With Arbitrary k-Space Trajectories
Klaas P.Pruessmann,1Markus Weiger,1Peter Bo¨rnert,2and Peter Boesiger1*
New,efficient reconstruction procedures are propod for nsi-tivity encoding(SENSE)with arbitrary k-space trajectories.The prented methods combine gridding principles with so-called conjugate-gradient iteration.In this fashion,the bulk of the work of reconstruction can be performed by fast Fourier transform (FFT),reducing the complexity of data processing to the same order of magnitude as in conventional gridding reconstruction. Using the propod method,SENSE becomes practical with non-standard k-space trajectories,enabling considerable scan time reduction with respect to mere gradient encoding.This is illus-trated by imaging simulations with spiral,radial,and random k-space patterns.Simulations were also ud for investigating the convergence behavior of the propod algorithm and its depen-dence on the factor by which gradient encoding is reduced.The in vivo feasibility of non-Cartesian SENSE imaging with iterative re-construction is demonstrated by examples of brain and cardiac imaging using spiral trajectories.In brain imaging with six receiver coils,the number of spiral interleaves was reduced by factors ranging from2to6.In cardiac real-time imaging with four coils, spiral SENSE permitted reducing the scan time per image from 112ms to56ms,thus doubling the frame-rate.Magn Reson Med 46:638–651,2001.©2001Wiley-Liss,Inc.
Key words:nsitivity encoding;SENSE;k-space trajectories; spiral imaging;iterative reconstruction
Arrays of simultaneously operated receiver coils have re-cently attracted increasing attention as a means of enhanc-ing scan speed in MRI.Several authors have propod imaging schemes designed to utilize parallel signal acqui-sition with multiple coils for the purpo of reducing scan time(1–8).Exploiting knowledge of spatial coil nsitiv-ity,the techniques enable the reduction of the number of gradient-encoding steps without compromising spatial resolution or thefield of view(FOV).
Due to the specific role of coil nsitivity,array imaging with reduced k-space sampling requires particular consid-erations for image generation.As coil nsitivity is ud as a means of signal encoding,the net encoding functions are no longer harmonics of the FOV as in conventional Fourier imaging.Conquently,image reconstruction from under-sampled multiple-coil data cannot be performed by mere Fourier transform(FT)but requires specialized processing.
A general approach to the reconstruction problem was described in Ref.6,introducing the SENSitivity Encoding (SENSE)method.In this work,reconstruction formulae were derived for array imaging with arbitrary coil config-urations and k-space trajectories.However,in their general form the formulae are numerically rather demanding. An efficient implementation was previously described for the spec
ial ca of sampling k-space along a regular Cartesian grid(6).In this ca,SENSE processing may be performed in a fashion similar to earlier Cartesian ap-proaches(2,4).In afirst step,conventional Fourier recon-struction is performed for each coil element,yielding a t of aliad images with reduced FOV.In a cond step,the aliasing is undone using knowledge of the individual coil nsitivities.This latter step requires little computation becau the aliasing effect superimpos the signal of only a small number of volume elements in each pixel of the reduced FOV.In other words,the point spread function (PSF),which corresponds to Cartesian undersampling,is well behaved in that it exhibits sharply localized,equidis-tant peaks.
With general k-space trajectories,such as spirals(9,10) or radial(11,12)and stochastic(13)schemes,the situation is more complicated.Similarly to the Cartesian ca, SENSE permits reducing the density of general sampling patterns while maintaining the covered k-space area.How-ever,arbitrary k-space patterns entail highly variable PSFs. With spiral acquisition,for example,the PSF exhibits con-tinuous,ring-shaped structures rather than evenly spaced, isolated peaks(14).As a conquence,the aliasing that results from reducing k-space density is considerably more complex and prevents straightforward unfolding in the image domain.This obrvation illustrates why general SENSE reconstruction has so far been riously hampered by its numerical complexity(15).
In this work we propo a novel,highly efficient imple-mentation of non-Cartesian SENSE processing.The new approach is bad on the idea of performing reconstruction iteratively,as recently propod independently by the au-thors and by Kannengieer et al.(16,17).It is shown that so-called weak reconstruction with SNR optimization can be accomplished using the conjugate-gradient(CG)itera-tion method(18).The key improvement in reconstruction speed is then achieved by combining fast Fourier trans-form(FFT)with forward and rever gridding operations for efficient execution of the CG iteration loops.Further speed benefits are accomplished by veral measures of preconditioning and numerical optimization.
Using the propod procedure,the computation times required for non-Cartesian SENSE reconstruction are con-siderably diminished,making nsitivity-encoded imag-ing practical with arbitrary acquisition patterns.The gen-eral applicability of the method has been verified by im-aging simulations with spiral,radial,and random trajectories,as well as Cartesian sampling for comparison. The in vivo feasibility of the new approach is demon-
1Institute of Biomedical Engineering,University of Zu¨rich and Swiss Federal
Institute of Technology Zu¨rich,Switzerland.
2Philips Rearch Laboratory Hamburg,Division Technical Systems,Ham-
burg,Germany.
Grant sponsor:EUREKA;Grant number:EU1353;Grant sponsor:KTI;Grant
number:3030.2.
*Correspondence to:Dr.P.Boesiger,Institute of Biomedical Engineering,
University of Zu¨rich and ETH Zu¨rich,Gloriastras35,CH-8092Zu¨rich,Swit-
zerland.
Received22November2000;revid17January2001;accepted5February
2001.
2001ISMRM Young Investigator I.I.Rabi Award Finalist.
Magnetic Resonance in Medicine46:638–651(2001)©2001Wiley-Liss,Inc.638
strated by examples of spiral SENSE imaging of the brain and heart.
THEORY AND METHODS
Reconstruction Problem
新疆帽
The problem of nsitivity-bad image reconstruction from multiple-coil data,obtained with general k-space tra-jectories,was previously studied theoretically in Ref.
6.The concepts propod in that work were bad on considerations of voxel functions,reflecting the spatial signal respon associated with pixel values.Two strate-gies for ensuring suitable voxel functions,both character-ized by underlying voxel conditions,were suggested.The first approach strictly requires the resulting voxel func-tions to be best approximations of given , box or Dirac shapes.The cond approach us a weaker voxel condition,requiring only that the resulting voxel functions share the orthogonality relations of their ideal counterparts.
This latter approach forms the basis of the prent work. Following the notation ud in Ref.6,the weak voxel condition reads:
FEϭId[1] where E and F denote the encoding and reconstruction matrices,respectively,and Id denote
s the identity matrix. The reconstruction matrix F reprents the linear mapping of sample values ud for creating the desired image.The encoding matrix E is formed by the scalar products of the hybrid gradient and nsitivity encoding functions and the chon voxel prototypes.In the scope of this work,Dirac distributions are ud as voxel prototypes,forming a reg-ular Cartesian grid corresponding to a common NϫN image matrix.With the Dirac choice,the weak voxel con-dition may be restated more intuitively as:each voxel function must be equal to1in the center of the voxel to which it belongs,and equal to zero in all other voxels’centers.The entries of the encoding matrix then reduce to
E͑␥,͒,ϭe i krs␥͑r͒,[2] where rdenotes the position of the-th voxel,kthe-th sampling position in k-space,and s␥the complex spatial nsitivity of the␥-th coil.
Let n
C and n
K
denote the number of receiver coils ud
and the number of sampling positions in k-space,respec-
tively.Often the total number of encoding functions,n
C n K
,
exceeds N2.Then Eq.[1]is underdetermined,leaving de-grees of freedom in the choice of the reconstruction matrix, which may be exploited for optimizing SNR.As derived in Ref.6,in terms of SNR the best solution of Eq.[1]is given by
Fϭ͑E H⌿˜Ϫ1E͒Ϫ1E H⌿˜Ϫ1,[3] where the superscript H indicates the complex conjugate transpo,and⌿˜denotes the sample noi matrix,which describes the levels and the correlation of stochastic noi in signal samples.Using this expression for the reconstruc-tion matrix,image reconstruction is performed by:
ϭ͑E H⌿˜Ϫ1E͒Ϫ1E H⌿˜Ϫ1m,[4]
where m and v are vectors of length n
C
n
K
and N2,respec-tively,listing the complex sample values acquired and the complex pixel values of the reconstructed image.Stochas-tic noi in thefinal image is described by the image noi matrix:
XϭF⌿˜F Hϭ͑E H⌿˜Ϫ1E͒Ϫ1.[5]
For efficient numerical evaluation,Eq.[4]is rewritten as a linear system of equations with the unknown image vector v:
͑E H⌿˜Ϫ1E͒ϭE H⌿˜Ϫ1m.[6] Iterative Solving
Solving the system described by Eq.[6]is numerically challenging,as it is usually rather large.The size of E is
n
C
n
K
ϫN2,and the matrix in brackets has size N2ϫN2. Handling the matrices in a straightforward manner re-quires large amounts of memory and processing time.In particular,using direct methods for solving the equation entails massive computation with a number of operations on the order of N6.Calculating a128ϫ, would require more than1012operations and about1GB of memory,making image reconstruction rather tedious with prent-day workstations.
One key to more efficient processing is solving Eq.[6] iteratively.Starting with some initial vector,iterative al-gorithms yield a progression of approximate solutions, which converge to the exact solution.A variety of such techniques exist for the treatment of large linear systems. For the prent ca,the so-called conjugate-gradient(CG) method(18)is particularly suited.On the one hand,as shown later in this ction,it may be combined with FFT for very efficient calculations.On the other hand,the CG iteration does not require particular provisions for ensur-ing convergence.It converges safely given that the coeffi-cient matrix involved is positive definite,which holds for the system descri
bed by Eq.[6](e Appendix A).
The CG algorithm theoretically yields the exact solution of an N2ϫN2system after at most N2iterations.For N in the range of128,though,it is not practical to carry out the entire procedure.However,for well-behaved problems, sufficiently accurate approximations are obtained after a relatively small number of iterations(19).
Each CG iteration step consists in multiplying the matrix to be inverted with a residuum vector and veral further calculations of minor complexity.Thus the iteration speed depends crucially on how fast matrix-vector multiplica-tion can be performed.The number of iterations necessary to achieve a given accuracy is related to the so-called condition of the matrix to be inverted and the suitability of the starting vector.Fortunately,in the prent ca,the factors can be greatly influenced by veral measures of numerical optimization,as discusd in the following.
SENSE With Arbitrary k-Space Trajectories639
杨梅岭
Elimination of the Noi Matrix
合肥旅游攻略In afirst step,Eq.[6]is simplified by eliminating the sample noi matrix.Note that noi levels and corr
elation are considered exclusively for the sake of minimizing noi in thefinal image.The correct spatial assignment of the MR signal is independently ensured by the voxel con-dition in Eq.[1],which is fulfilled also if the noi statis-tics are not correctly taken into account.This can be easily verified by inrting Eq.[3]into Eq.[1].
幼儿营养食谱Hence,one option for eliminating the noi matrix is to skip noi asssment and simply replace⌿˜by identity. As a result,noi in thefinal image will not be exactly minimized.The verity of this drawback is difficult to estimate a priori without actually calculating image noi. Generally,the benefit of considering noi statistics scales with encoding redundancy and is greater the less equiva-lent the coils are with respect to load,gain,mutual cou-pling,and electronic noi.
To eliminate the noi matrix from Eq.[6]without com-promising the SNR in thefinal image,it is necessary to consider the structure of⌿˜.This matrix has one row and one column for each hybrid encoding carried out,thus its
size is n可以光脚跳绳吗
C n
K
ϫn C n K.Each diagonal entry reprents the
variance of noi in the corresponding sample value,while the off-diagonal elements reflect noi correlation between samples.
Noi correlation occurs mainly between samples taken simultaneously with different coils.On the one hand,coils with overlapping nsitive regions receive partly identical thermal radiation.On the other hand,coil coupling gives ri to mutual noi transmission.To a certain degree, noi correlation does also occur between successive sam-ples acquired with the same coil,especially when over-sampling is applied.However,this effect shall be ne-glected in the scope of this work.Furthermore,it is as-sumed that noi correlation between different coils is constant over time,as the underlying mechanisms are approximately invariant.Receiver noi is then described by a spar matrix with a simple structure:
⌿˜͑␥,͒,͑␥Ј,Ј͒ϭ⌿␥,␥Ј␦,Ј.[7] In this formulation,the Kronecker delta symbol express the assumption that noi correlation occurs only between samples taken simultaneously at the same position in k-
space.The time-independent n
C ϫn C matrix⌿can be
女人的心思
determined experimentally by statistical analysis of refer-ence noi samples taken in the abnce of MR signal.Let ␥denote the noi output of the␥-th channel.Then the entries of⌿are given by:
⌿␥,␥Јϭ␥*␥Ј,[8] where the bar indicates time averaging.Note that this formula for the determination of⌿is statistically equiva-lent to Eq.[A9]derived in Ref.6.Using averaging accord-ing to Eq.[8],the relative error in each entry of⌿decreas only as the inver square root of the number of noi samples taken.Thus,for reasonable accuracy a number of noi samples in the range of103is advisable,which can readily be obtained during the common scan preparation pha.Note that noi correlation,just like coil nsitivity, depends on the particular coil and load configuration. Therefore,it is important to collect noi samples in the actual imaging tup.
With noi statistics simplified according to Eq.[7],the noi variance matrix can be eliminated from Eq.[6]by a simple manipulation.The basic idea is to create a t of virtual receiver channels by linear combination of the original ones,such that the virtual channels exhibit unit noi levels and no mutual noi correlation.As shown in Appendix B,suitable weighting coefficients for this pur-po are given by the inver of the matrix L obtained by Cholesky decomposition(20)of⌿:
半岛都市报电子版
⌿ϭLL H.[9] Virtual sampling data with decorrelated unit noi are thus obtained from the original samples by
m␥,decorrϭ␥Ј͑LϪ1͒␥,␥Јm␥Ј,.[10] The net coil nsitivities associated with the virtual chan-nels are given likewi by
s␥decorr͑r͒ϭ␥Ј͑LϪ1͒␥,␥Јs␥Ј͑r͒,[11] leading to the modified encoding matrix
E͑␥,͒,decorrϭe i krs␥decorr͑r͒.[12] With sample values and nsitivities modified in this fash-ion,the newly combined channels can be treated exactly like physical ones.The noi variance matrix of the virtual channels is equal to identity and can thus be omitted when Eq.[6]is reformulated for image reconstruction from the modified data.Note that by the transition to virtual chan-nels the solution of Eq.[6]is not altered(e Appendix B). In particular,the optimization of SNR is prerved.
Both options—decorrelation and neglecting noi corre-lation—lead to the same simplified formula.Dropping the superscript for decorrelation,it reads in both cas:
͑E H E͒ϭE H m.[13]
FT for Matrix-Vector Multiplication
In each CG loop,the matrices E and E H,which appear on the left side of Eq.[13],need to be multiplied with some temporary vectors,say,x and y,respectively.For perform-ing the operations efficiently,advantage is taken of the fact that,despite the nsitivity contribution to hybrid encoding,E is still characterized by Fourier terms.The (␥,)-th component of the product E x,reads:
͑E x͒͑␥,͒ϭe i krs␥͑r͒x.[14] The sum on the right may be regarded as the integral of a t of weighted Dirac distributions in the spatial domain
640Pruessmann et al.
͑E x͒͑␥,͒ϭ͵e i kr͑s␥͑r͒x␦͑rϪr͒͒d r,[15]
which is the inver FT of the distribution t,sampled at k:
͑E x͒͑␥,͒ϭFTϪ1ͫs␥͑r͒x␦͑rϪr͒ͬ͑k͒.[16] Similarly,the-th component of the product of E H and y,
͑E H y͒ϭ␥,eϪi krs*␥͑r͒y␥,,[17]
can be regarded as a sum of k-space integrals weighted by complex conjugate coil nsitivity,
͑E H y͒ϭ␥s*␥͑r͒
ͩ͵eϪi kr͑y͑␥,͒␦͑kϪk͒͒d kͪ,[18]
where the term in brackets reprents a forward FT,sam-pled at r:
͑E H y͒ϭ␥s*␥͑r͒
ͩFTͫ͑y͑␥,͒␦͑kϪk͒͒ͬ͑r͒ͪ.[19]
The FTs in Eqs.[16]and[19]can esntially be calculated by FFT.However,while the voxel positions,r,form a regular Cartesian grid compatible with FFT,the k-space positions,k,form an arbitrary sampling pattern.There-fore,additional gridding is required preceding and follow-ing the forward and rever transforms,respectively.This step can be accomplished by convolution and resampling following the same principles as ud for common grid-ding reconstruction(21).Note that for this operation the experimental trajectory need not fulfill the Nyquist density criterion.Nevertheless,the Cartesian k-space grid ud for resampling must be sufficiently den in order to avoid aliasing artifact in the calculated FTs.In both directions the FFT and gridding procedure must be repeated for each coil,as the functions to be transformed are coil-dependent. By applying Eqs.[16]and[19]successively,the matrix-vector multiplications required for iterative solving of Eq.
[13]can be performed very efficiently.Indeed,the FFT approach reduces the complexity of a single iteration step drastically from O(N)ϭN4to O(N)ϭN2logN.Note also that when using FFT for evaluating E H E x,the matrix E need not be explicitly calculated or stored.
Density and Intensity Correction
A common way of reducing the number of iterations re-quired in the CG method is so-called preconditioning(19). The basic idea of preconditioning is to manipulate the system of equations such that the right side approximates the solution of the system.The right side of Eq.[13]is obtained by multiplication of E H with the vector m of sample values.According to Eq.[19],this operation is somewhat similar to conventional processing of non-Car-tesian,multiple-coil data.For each coil element an inter-mediate image is created by gridding and FT(21),followed by nsitivity-weighted summation(22).
Thus,the right side of Eq.[13]reprents a rudimentary image reconstruction.Clearly,upon reduced gradient en-coding,the resulting image will exhibit aliasing artifact due to the violation of the Nyquist criterion by each t of single-coil data.Nevertheless,this image forms afirst ap-proximation of exact reconstruction.Hence,the prent problem is a priori favorably conditioned.However,un-like standard
gridding,Eq.[19]does not involve any cor-rection for inhomogeneous k-space density.Furthermore, in the image domain this equation does not imply what is commonly referred to as intensity ,the com-pensation for inhomogeneity of overall coil nsitivity.In the following it is shown how the known correction mechanisms can be incorporated into Eq.[13]for im-proved preconditioning.
Density correction of the sample vector m is performed by the diagonal matrix
D͑␥,͒,͑␥,͒ϭ
1
淘宝网店装修d͑k͒,[20]
where d(k)denotes the relative density of the k-space sampling pattern at the position k.Unfortunately,this matrix cannot be inrted into Eq.[13]in a straightforward manner,as any manipulation must be performed equally on both sides of the equation.Instead,the following for-mulation is ud:
͑E H DE͒ϭE H D m.[21]
This modification does not reprent an exact transform. However,note that the matrix D was put into the same places where earlier the noi matrix appeared.As stated previously,any invertible matrix can replace the noi matrix without violating the voxel condition,although it affects SNR.In fact,the introduction of D is equivalent to pretending that the noi level in sample values increas linearly with k-space density,thus misleading SNR opti-mization to a certain degree.
SNR optimization is bad on encoding redundancy, i.e.,on the fact that parts of information required for spa-tial resolution may be available more than just once,be it from different coils or from neighboring positions in k-space.In such cas the best way to average redundant information depends on the levels and correlation of the concomitant noi.On the one hand,the introduction of density correction in Eq.[21]does not change the relative estimation of noi between coils,but only between k-space positions.On the other hand,noi relations be-tween cloly neighboring positions in k-space are not changed significantly unless the density function d(k) exhibits considerable local changes.Therefore,the SNR drawback of density correction is expected to be negligible with regular types of trajectories,such as spirals and radial patterns,yet may be considerable ,random sam-pling.
SENSE With Arbitrary k-Space Trajectories641
Note that,as oppod to the noi matrix,which was
eliminated for the sake of numerical optimization,the
density correction matrix D is diagonal and therefore
caus only minimal extra computation.
The incorporation of intensity correction is relatively
simple,as it can be accomplished by exact transforms.
According to Eq.[19],the image E H D m is weighted by the
sum-of-squares of coil nsitivities.Let I denote an N2ϫN2diagonal matrix given by the inver square-root of this
weighting factor:
I,ϭ
1
ͱ␥͉s␥͑r͉͒2.[22]
Left-multiplication of Eq.[21]with I2yields
͑I2E H DE͒ϭI2E H D m.[23] In this formulation,the nonaliad components of the image on the right side are intensity-corrected.However, the matrix in brackets is no longer safely positive-definite, as required for using the CG method.To restore this prop-erty,the equation is left-multiplied by I–1,and on the left side the identity I I–1is inrted:
͑IE H DEI͒͑IϪ1͒ϭIE H D m.[24] Note that by this operation the right side and the interme-diate solution in brackets were both modified in the same ,left-multiplied by I–1.Therefore,the right side still reprents a good approximation,corresponding to favorable conditioning,while the leftmost matrix in brack-ets is now positive-definite(e Appendix A).
Using Eq.[24],both density and intensity correction can be applied in CG iteration.Both correction methods are optional and work independently.For switching either mechanism off,the respective matrix is simply replaced by identity.
Starting Image
Generally,for fast convergence the starting image of the CG iteration should be a good approximation of the exact solution.However,with preconditioning,as discusd previously,the best choice is simply a zero image.This is becau the approximation resulting from thefirst loop then is a suitably scaled copy of the right side of the system.In this fashion,the iteration is esntially started with the known approximate solution,yet without the need to consider its scaling,which involves numerous
factors,including N,n
C ,n
K
,and the units chon for coil
nsitivity and k-space density.
Implementation
According to the previous considerations,image recon-struction consists of esntially three steps.First,the right side of Eq.[24]is calculated:
aϭIE H D m,[25]yielding the intermediate image a.In the cond step,an approximate solution b
approx
of
͑IE H DEI͒bϭa[26]
is determined by CG iteration.An approximate solution of
Eq.[13]is then obtained by intensity correction of b
approx
:
approxϭI b approx.[27]
Figure1schematically shows the implementation of this procedure as ud for the prent work.Its central part is a CG unit,which controls the iteration process and per-forms minor calculations required by the CG scheme. The CG algorithm ud is formally described in Appendix C.In each iter
ation step,it requires the multiplication of the matrix IE H DEI with a vector residuum resulting from the previous iteration.This multiplication is carried out by the loop which forms the major part of Fig.1.First,the current residuum is multiplied by I.The cond step is the applica-tion of E,consisting of coil-wi multiplication by coil n-sitivity,followed by individual FFT and gridding along the experimental k-space trajectory.The transform and gridding steps are jointly reprented as FT2.In k-space,density cor-rection is carried out independently for each channel.The following application of E H consists of coil-wi gridding along the Cartesian grid and FFT back to image domain (FT1),followed by individual weighting by complex conju-gate coil nsitivity and subquent summation of the single-coil results.Finally,the sum is again multiplied by I and then fed back into the CG unit.
For the initial calculation of a,advantage is taken of the fact that the matrices applied to the sample vector m in Eq.
[25]are identical to the leftmost three matrices in Eq.[26]. The sample data obtained from the receiver channels are fed into the corresponding k-space paths,before density correction.The resulting image a is ud for initializing
the CG unit,which starts the iteration process with b
approx
(0)
ϭ0.
Each iteration loop results in a refined approximation
b
approx
(i)to the exact solution of Eq.[26].As soon as the process has converged,the current approximation is out-put and multiplied by I to yield the intensity-corrected
reconstruction v
approx
.
As the exact solution is not known,the accuracy of a current approximation can only be estimated.One plausi-ble measure of accuracy is the relative violation of Eq.[26]:
␦ϭ
͉IE H DEI b approxϪa͉
͉a͉,[28]
which can be calculated with little effort as a by-product of the CG operations(Appendix C).Here,the brackets|⅐| denote the Euclidean vector length.The relative deviation from the exact solution,
⌬ϭ
͉I b approxϪ͉
͉͉,[29]
is interesting for studies of convergence behavior,yet re-quires either knowledge of the exact solution v or of a
642Pruessmann et al.
