EIGEN-VECTOR APPROACH FOR COIL SENSITIVITY MAPS ESTIMATION

Abstract
A method for estimating a coil sensitivity map for a magnetic resonance (MR) image includes providing (61) a matrix A of sliding blocks of a 2D image of coil calibration data, calculating (62) a left singular matrix V∥ from a singular value decomposition of A corresponding to τ leading singular values, calculating (63) P=V∥V∥H, calculating (64) a matrix S that is an inverse Fourier transform of a zero-padded matrix P, and solving (65) MHcr=(Sr)Hcr for cr, where cr is a vector of coil sensitivity maps for all coils at spatial location r, and
Description
TECHNICAL FIELD

This disclosure is directed to methods for estimating coil sensitivity maps (CSM) of magnetic resonance imaging (MRI) apparatuses.


DISCUSSION OF THE RELATED ART

Parallel imaging makes use of multiple receiver coils to acquire the image in parallel. It can be used to accelerate the image acquisition by exploiting the spatially varying sensitivities of the multiple receiver coils since each coil image is weighted differently by the coil sensitivity maps (CSM). The accurate estimation of coil sensitivity maps can ensure the success of approaches that make use of the sensitivity maps either in the reconstruction formulation or in coil combination. Some methods explicitly make use of the coil sensitivities in reconstruction, and the goodness of the CSM's is important to the quality of the reconstructed image. Other approaches implicitly make use of the coil sensitivities by performing an auto-calibrating coil-by-coil reconstruction, but the CSM's are needed if one wants to obtain the coil combined complex image or the phase image.


The most common way to determine the sensitivity maps is to obtain low-resolution pre-scans. However, when the object is not static, the sensitivity functions are different between pre-scan and under-sampled scans, and this could lead to reconstruction errors. To compensate for this, joint estimation approaches have been proposed, however, these approaches usually have high computation cost and are restricted to explicit reconstructions.


The eigenvector approach proposed tries to build a connection between implicit and explicit approaches, by showing that the CSM can be computed with the auto-calibrated coil-by-coil reconstruction. It has been shown that the coil sensitivities can be computed as the eigenvector of a given matrix in the image space corresponding to the eigenvalue's ‘1”s. However, to the best of the inventor's knowledge: (1) the detailed mathematical derivations for the eigenvector approach are not well understood; (2) the optimization criterion for computing the CSM is not very clear; and (3) there lacks an efficient approach.


SUMMARY

Exemplary embodiments of the invention as described herein generally include methods for coil sensitivity maps (CSM) estimation for 2-D MR images. Embodiments of the invention mathematically derive an eigenvector approach and show that the coil sensitivity maps at a given spatial location can be obtained by solving a generalized eigenvalue system. Embodiments of the invention establish relationships between the generalized eigenvalue system and two related eigenvalue systems where the associated matrices are Hermitian. Embodiments of the invention can reduce the computational and storage costs and avoid the computation of large matrices by using equivalent representations. Experimental results on simulated and real data sets show that, a method according to an embodiment of the invention can obtain CSM's with high quality and the resulting reconstructed images have fewer artifacts compared to those reconstructed with CSM's obtained by a B1-Map.


According to an aspect of the invention, there is provided a method for estimating a coil sensitivity map for a magnetic resonance (MR) image, including providing a matrix A of sliding blocks of a 2D nx×ny image of coil calibration data, where







A
=

(




a

1
,
1





a

1
,
2








a

1
,

n
b








a

2
,
1





a

2
,
2








a

2
,

n
b
























a


n

c
,



1





a


n
c

,
2








a


n
c

,

n
b






)


,




nc is a number of coils, nb is a number of sliding blocks extracted from the coil calibration data, and ai,j is a kxky×1 column vector that represents a jth sliding block of an ith coil, determining a unit eigenvector α that maximizes custom-character(Sr)Hα,MHαcustom-character, where Sr is an inverse Fourier transform of a zero-padded matrix P=VVH at spatial location r in the 2D nx×ny image, where







V





=


[





v
1

,





v
2

,





,




v
τ




]

=

(




v

1
,
1





v

1
,
2








v

1
,
τ







v

2
,
1





v

2
,
2








v

2
,
τ























v


n
c

,
1





v


n
c

,
2








v


n
c

,
τ





)






is a matrix composed of left singular vectors of A corresponding to τ leading singular values where vi,k is a kxky×1 column vector that is an i-th block in vector vk,






M


(


(



1


1





1




0


0





0




















0


0





0



)



(



0


0





0




1


1





1




















0


0





0



)













(



0


0





0




0


0





0




















1


1





1



)


)





is a nc×kxkync sparse matrix of the same size as Sr, and finding an optimizer cr that maximizes a correlation between (Sr)Hα and MHα where the optimizer cr is a column vector of coil sensitivities of a coil at spatial position r in the 2D nx×ny image.


According to another aspect of the invention, there is provided a method for estimating a coil sensitivity map for a magnetic resonance (MR) image, including providing a matrix A of sliding blocks of a 2D nx×ny image of coil calibration data, where







A
=

(




a

1
,
1





a

1
,
2








a

1
,

n
b








a

2
,
1





a

2
,
2








a

2
,

n
b
























a


n

c
,



1





a


n
c

,
2








a


n
c

,

n
b






)


,




nc is a number of coils, nb is a number of sliding blocks extracted from the coil calibration data, and ai,j is a kxky×1 column vector that represents a jth sliding block of an ith coil, calculating a left singular matrix V from a singular value decomposition of A, where A=VΣUH and







V





=


[





v
1

,





v
2

,





,




v
τ




]

=

(




v

1
,
1





v

1
,
2








v

1
,
τ







v

2
,
1





v

2
,
2








v

2
,
τ























v


n
c

,
1





v


n
c

,
2








v


n
c

,
τ





)






is a matrix of left singular vectors of A corresponding to τ leading singular values where vi,k is a kxky×1 column vector that is an i-th block in vector Vk, calculating






P
=



V




V

H


=

(


(




p

1
,
1
,
1





p

1
,
1
,
2








p

1
,
1
,


k
x



k
y









p

1
,
1
,
1





p

1
,
2
,
2








p

1
,
2
,


k
x



k
y























p

1
,

n
c

,
1





p

1
,

n
c

,
2








p

1
,

n
c

,


k
x



k
y







)













(




p


n
c

,
1
,
1








p


n
c

,
1
,


k
x



k
y









p


n
c

,
2
,
1








p


n
c

,
2
,


k
x



k
y




















p


n
c

,

n
c

,
1








p


n
c

,

n
c

,


k
x



k
y







)


)






where pi,j,t is a kxky×1 column vector, calculating Si,j,t=FH(Pt(pi,j,t)) where FH represents an inverse Fourier transform and Pt represents a zero-padding operator, and solving MHcr=(Sr)Hcr for cr, where cr is a vector of coil sensitivity maps for all coils at spatial location r,







S
r

=


(


(




s

1
,
1
,
1

r




s

1
,
1
,
2

r







s

1
,
1
,


k
x



k
y



r






s

1
,
2
,
1

r




s

1
,
2
,
2

r









s

1
,
2
,


k
x



k
y



r




















s

1
,

n
c

,
1

r




s

1
,

n
c

,
2

r







s

1
,

n
c

,


k
x



k
y



r




)













(




s


n
c

,
1
,
1

r







s


n
c

,
1
,


k
x



k

y








r






s


n
c

,
2
,
1

r







s


n
c

,
2
,


k
x



k

y








r

















s


n
c

,

n
c

,
1

r







s


n
c

,

n
c

,


k
x



k

y








r




)


)






and







M
=


(


(



1


1





1




0


0





0




















0


0





0



)



(



0


0





0




1


1





1




















0


0





0



)













(



0


0





0




0


0





0




















1


1





1



)


)

.





According to a further aspect of the invention, the method includes finding an optimizer cr that maximizes a correlation between (Sr)Hα and MHα where the optimizer cr is vector of coil sensitivity maps for all coils at spatial location r in the 2D nx×ny image.


According to a further aspect of the invention, solving MHcr=(Sr)Hcr comprises solving eigenvalue equations











M


(

S
r

)


H



k
x



k
y




β

=


λ
_


β


,




and (Sr(Sr)H)γ={tilde over (λ)}γ, wherein λ and {tilde over (λ)} denote eigenvalues, β and γ denote eigenvectors.


According to a further aspect of the invention, the method includes calculating









S
r



M
H






using






m

i
,
j



=





t
=
1



k
x



k
y









s

i
,
j
,
t



=





t
=
1



k
x



k
y










F
H



(


P
t



(

p

i
,
j
,
t


)


)



=


F
H



(




t
=
1



k
x



k
y










P
t



(

p

i
,
j
,
t


)



)





,




where mi,j is an nxny×1 column vector containing the entries of mi,jr at all spatial locations, and mi,jr is the i,j elements of the 2D MR image at spatial location r.


According to a further aspect of the invention, the method includes calculating Sr(Sr)H as Sr(Sr)H=kxkyZr(Zr)H, where







Z
r

=

(




z

1
,
1

r




z

1
,
2

r









z

1
,
τ

r






z

2
,
1

r




z

2
,
2

r









z

2
,
τ

r




























z


n
c

,
1

r




z


n
c

,
2

r









z


n
c

,
τ

r




)





is an nc×τ matrix with the j-row and k-th column being Zj,kr, Zj,kr denotes the r-th entry of Zj,k=FH (Pc(vj,k)) and Pc denotes the zero-padding operator that maps a center of vj,k to a center of zj,k.


According to a further aspect of the invention, the method includes computing the i-th row and j-th column of Zr(Zr)H as









k
=
1

τ








z

i
,
k

r




z

i
,
k

r

_






where zi,kr denotes the i-th row and k-th column in Zr.


According to a further aspect of the invention, the coil calibration data is acquired from a center square of frequency space data.


According to a further aspect of the invention, the matrix A is constructed from calibration data of size cx×cy×cz by gathering sliding blocks of kernel size kx×ky×kz, where A has [(cx−kx+1)×(cy−ky+1)×(cz−kz+1)] columns and kxkykznc rows.


According to another aspect of the invention, there is provided a non-transitory program storage device readable by a computer, tangibly embodying a program of instructions executed by the computer to perform the method steps for estimating a coil sensitivity map for a magnetic resonance (MR) image.





BRIEF DESCRIPTION OF THE DRAWINGS


FIG. 1 is a table summarizing the notation used in this disclosure, according to an embodiment of the invention.



FIG. 2 illustrates zero padding for the case kx=ky=3, according to an embodiment of the invention.



FIG. 3 illustrates an eigenvector approach where each k-space point is represented as a linear combination of the k-space points of all the coils in its neighborhood, according to an embodiment of the invention.



FIG. 4 illustrates the principle of an eigenvector approach of EQ. (13), according to an embodiment of the invention.



FIG. 5 is a table summarizing the value of mean absolute correlation (mac) under different number of kept singular vectors and different kernel sizes, according to an embodiment of the invention.



FIG. 6 is a flowchart of an eigenvector approach to estimating coil sensitivity maps, according to an embodiment of the invention.



FIGS. 7(
a)-(d) depict CSM's obtained by various approaches, according to an embodiment of the invention.



FIGS. 8(
a)-(c) are graphs of the mean absolute correlation as a function of different number of kept vectors and different kernel sizes, according to an embodiment of the invention.



FIG. 9 illustrates a reconstruction of a 2D cardiac image, according to an embodiment of the invention.



FIGS. 10(
a)-(b) illustrates a reconstruction of a 3D abdomen image, according to an embodiment of the invention.



FIG. 11 illustrates a typical Cartesian sampling pattern used in dynamic cardiac MR imaging, according to an embodiment of the invention.



FIGS. 12(
a)-(b) illustrate the influence of the amount of used calibration data on the reconstruction, according to an embodiment of the invention.



FIG. 13 compares the square and the rectangular calibration data in terms of the mac criterion, according to an embodiment of the invention.



FIG. 14 illustrates the conversion of the calibration data into a calibration matrix for 3D CSM estimation, according to an embodiment of the invention.



FIGS. 15(
a)-(h) illustrate a comparison of a 2D eigenvector approach and a 3D eigenvector approach according to an embodiment of the invention.



FIG. 16 is a block diagram of an exemplary computer system for implementing an eigenvector approach to coil sensitivity maps estimation for 2-D MR images, according to an embodiment of the invention.





DETAILED DESCRIPTION OF EXEMPLARY EMBODIMENTS

Exemplary embodiments of the invention as described herein generally include systems and methods for an eigenvector approach to coil sensitivity maps (CSM) estimation for 2-D MR images. Accordingly, while the invention is susceptible to various modifications and alternative forms, specific embodiments thereof are shown by way of example in the drawings and will herein be described in detail. It should be understood, however, that there is no intent to limit the invention to the particular forms disclosed, but on the contrary, the invention is to cover all modifications, equivalents, and alternatives falling within the spirit and scope of the invention.


As used herein, the term “image” refers to multi-dimensional data composed of discrete image elements (e.g., pixels for 2-dimensional images and voxels for 3-dimensional images). The image may be, for example, a medical image of a subject collected by computer tomography, magnetic resonance imaging, ultrasound, or any other medical imaging system known to one of skill in the art. The image may also be provided from non-medical contexts, such as, for example, remote sensing systems, electron microscopy, etc. Although an image can be thought of as a function from R3 to R or R7, the methods of the inventions are not limited to such images, and can be applied to images of any dimension, e.g., a 2-dimensional picture or a 3-dimensional volume. For a 2- or 3-dimensional image, the domain of the image is typically a 2- or 3-dimensional rectangular array, wherein each pixel or voxel can be addressed with reference to a set of 2 or 3 mutually orthogonal axes. The terms “digital” and “digitized” as used herein will refer to images or volumes, as appropriate, in a digital or digitized format acquired via a digital acquisition system or via conversion from an analog image.


Notation:

Throughout this disclosure, scalars are denoted by italic lowercase letters (e.g., nc, nx, ny, kx, ky), vectors by boldface lowercase letters (e.g., xi, ci), and matrices by italic uppercase letters (e.g., S Sr). The main notations used in this paper are summarized in Table I, shown in FIG. 1.


Let { }H denote the complex-conjugate transpose operator, x denote the complex-conjugate of a scalar x, and x be the complex-conjugate of a vector x.


Let i and j denote the indices for the coils in the interval [1, nc], r denote the spatial location in the 2D image of size nx×ny, and t denote the spatial location in the image block of size kx×ky.


For an nxny×1 column vector xi, representing an image by concatenating its columns, xir corresponds to its value at the spatial location r, where r goes from 1 to nxny. For a kxky×1 column vector vi,k, representing an image block by concatenating its columns, vi,kt denotes its value at the spatial location t, where t goes from 1 to kxky.


Let y=Pt(x) denote the zero padding operator defined as follows. The input x, a column vector of size kxky×1, represents a kx×ky image block by concatenating its columns; t, a scalar in the interval [1, kxky], determines which value in x to be put at the center of the matrix of size nx×ny; the output y is an nx×ny×1 column vector representing an nx×ny matrix, and its non-zero entries are the values of x. In particular, we let y=Pc(x) denote the zero padding operator that put the center point of x to the center location of y. For an image of size nx×ny, define the index of its center location as












n
x

+
1

2



+


(






n
y

+
1

2



-
1

)



n
x



;




and for an image block of size kx×ky, define the index of its center location as











k
x

+
1

2



+


(






k
y

+
1

2



-
1

)




k
x

.






Embodiments of the present disclosure assume that nx and ny are even; but kx and ky can be either even or odd.



FIG. 2 illustrates zero padding for the case kx=ky=3. The nine squares correspond to the input x, the square 21 is placed at the center of the nx×ny image, t in custom-character determines which value in x to be put at the center of the 2D image; the remaining entries are filled with zero except for the nine squares; and the three rectangles 22, 23, 24 denote the outputs of P1, P5, and P9, respectively.


Eigenvector Approach

One eigenvector approach is illustrated in FIG. 3, where each k-space point is represented as a linear combination of the k-space points of all the coils in its neighborhood. Let xir denote the value of xi at the spatial location r. Let wi,j be the convolution kernel for representing the k-space point of the i-th coil by the j-th coil. Mathematically, xir is estimated as:











x
i
r

=




j
=
1


n
c









w

i
,
j


*

x
j




,


r

,




(
1
)







where * denotes the convolution operator. Define






q
i,j
=F
H(PC(wi,j)).  (2)


qi,j is the result of applying an inverse Fourier transformation to the convolution kernel wi,j with zero padding, illustrated in FIG. 2. Applying the inverse Fourier transform to both sides of the equation in FIG. 3, one has











F
H



(

x
i

)


=




j
=
1


n
c









q
ij




F
H



(

x
j

)








(
3
)







where the convolution theorem has been applied. FH(xi)=mcustom-characterci can be inserted into EQ. (3) to obtain











c
i

=




j
=
1


n
c









q
ij



c
j




,





i
=
1

,
2
,





,


n
c

.





(
4
)







EQ. (4) indicates that all the entries in this equation can be decoupled with respect to the spatial location. Let cir, qi,jr and cjr denote the entries of ci, qi,j, and cj at the spatial location r, respectively. EQ. (4) can be rewritten as











c
i
r

=




j
=
1


n
c









q

i
,
j

r



c
j
r




,





i
=
1

,
2
,





,

n
c

,



r
.






(
5
)







Denote Qr as a matrix with the value at the i-row and the j-th column being qi,jr. Let cr be an nc×1 column vector composed of ci's at the spatial location r. EQ. (5) can be rewritten as the following compact matrix format:






c
r
=Q
r
c
r
,∀r  (6)



FIG. 3 illustrates the principle of CSM: the three rectangles denote the k-space data of the i-th coil, the first coil, and the last coil, respectively; xir the square 31 in the left plot, is estimated as the linear combination of the k-space points of all the coils in its neighborhood of size kx×ky (here kx=ky=3); the nine small squares 32 in X1 and 33 in Xnc correspond to the values in the convolution kernel wi,l and wi,nc.


In an eigenvector approach according to an embodiment of the invention, sliding blocks of size kx×ky) are gathered from the calibration data to form a matrix:









A
=

(




a

1
,
1





a

1
,
2








a

1
,

n
b








a

2
,
1





a

2
,
2








a

2
,

n
b
























a


n
c

,
1





a


n
c

,
2








a


n
c

,

n
b






)





(
7
)







where nb denotes the number of sliding blocks extracted from the calibration data, and ai,k, which is a kxky×1 column vector, denotes the k-th sliding block of the i-th coil.


Let






A=VΣU
H  (8)


be the Singular Value Decomposition of A. Denote









V





=


[


v
1

,

v
2

,





,

v
τ


]

=

(




v

1
,
1





v

1
,
2








v

1
,
τ







v

2
,
1





v

2
,
2








v

2
,
τ























v


n
c

,
1





v


n
c

,
2








v


n
c

,
τ





)






(
9
)







as a matrix composed of the left singular vectors of A corresponding to the leading τ singular values. Here, vi,k, a kxky×1 column vector, is the i-th block in the vector vk. It can be observed from EQS. (7) and (9) that, A and V share a similar structure in that: (1) the block ai,k and vi,k have the same size (kxky×1); and (2) the number of such blocks for each column is same.


Let P=VVH, and denote









P
=

(


(




p

1
,
1
,
1





p

1
,
1
,
2








p

1
,
1
,


k
x



k
y









p

1
,
1
,
1





p

1
,
2
,
2














p

1
,
2
,


k
x



k
y























p

1
,

n
c

,
1





p

1
,

n
c

,
2








p

1
,

n
c

,


k
x



k
y







)













(




p


n
c

,
1
,
1








p


n
c

,
1
,


k
x



k

y














p


n
c

,
2
,
1








p


n
c

,
2
,


k
x



k

y

























p


n
c

,

n
c

,
1








p


n
c

,

n
c

,


k
x



k

y












)


)





(
10
)







where pi,j,t is a kxky×1 column vector. Embodiments have used a similar block representation for P as A and V in that they all have nc “rows”. Note that the right hand side of EQ. (10) also represents PH, as P is Hermitian (i.e., P=PH).


As V and V have orthogonal column vectors,






V

V

H
A=V

V

H
ΣU
H
=[V∥0]ΣUH≈A,  (11)


where I denotes an identity matrix. When τ equals the rank of A, the “≈” becomes “=”. EQ. (11) indicates that, each column of A lies in the space spanned by orthogonal projection VVH when τ is set to the rank of A.


Next, embodiments generate EQ. (11) for all k-space blocks. Specifically, let xi,r be a kxky×1 column vector that denotes the k-space block of the i-th coil at the spatial location r. Assuming that the k-space blocks lie in the space spanned by V,











(




x

1
,
r







x

2
,
r












x


n
c

,
r





)

=


V








V




H



(




x

1
,
r







x

2
,
r












x


n
c

,
r





)




,



r
.






(
12
)







The case when EQ. (12) holds approximately will be discussed in the next subsection.


Picking a column of P defined in EQ. (10) and making use of P=PH=VVH P, EQ. (12) yields











x

i
,
r

t

=




j
=
1


n
c










(

p

i
,
j
,
t


)

H



x

i
,
r





,


r

,




(
13
)







where Xi,rt denotes the t-th element in the vector xi,r. FIG. 4 illustrates the principle of an eigenvector approach according to an embodiment of the invention of EQ. (13) with kx=ky=3: xi,rt denoted as the square 41 in the left plot, is represented by the summation of the convolutions between the k-space data of the j-th coil with kernels (pi:j:t)H for j=1, 2, . . . , nc; the plots in the first, second, and third row correspond to t=1, 5, 9, respectively; and the nine small squares 42 in X1 and 43 in Xnc, which denote (pi,j,t)H, constitute the convolution kernel. Note that, by moving the sliding block xi,r over all the spatial locations r=1, 2, . . . , nxny, the left hand side of EQ. (13) covers the k-space data of the i-th coil; and the right hand side of EQ. (13) corresponds to the summation of the convolutions between the k-space data of the j-th coil and the kernel (pi:j:t)H for j=1, . . . , nc.


For a given i and t, applying the inverse Fourier transformation to both sides of EQ. (13) yields











c
i

=




j
=
1


n
c










s

i
,
j
,
t


_



c
j




,




(
14
)







where i=1, 2, . . . , nc; t=1, 2, . . . , kxky, and






s
i,j,t
=F
H(Pt(pi,j,t))  (15)


is the result of applying the inverse Fourier transform to pi,j,t with proper zero padding, and si,j,t denotes the complex conjugate of si,j,t.


EQ. (14) indicates that the values at different spatial locations can be decoupled. Let cir, si,j,tr and cjr denote the entries of ci, si,j,t, and cj at the spatial location r, respectively. Embodiments can rewrite EQ. (14) as











c
i
r

=




j
=
1


n
c










s

i
,
j
,
t

r

_



c
j
r




,




(
16
)







where i=1, 2, . . . , nc; t=1, 2, . . . , kxky.


Let cr be a column vector composed of ci's at the spatial location r. Then,











S
r

=

(


(




s

1
,
1
,
1

r




s

1
,
1
,
2

r







s

1
,
1
,


k
x



k
y



r






s

1
,
2
,
1

r




s

1
,
2
,
2

r









s

1
,
2
,


k
x



k
y



r




















s

1
,

n
c

,
1

r




s

1
,

n
c

,
2

r







s

1
,

n
c

,


k
x



k
y



r




)






(




s


n
c

,
1
,
1

r







s


n
c

,
1
,


k
x



k
y



r






s


n
c

,
2
,
1

r







s


n
c

,
2
,


k
x



k
y



2

















s


n
c

,

n
c

,
1

r







s


n
c

,

n
c

,


k
x



k
y



r




)



)







Let




(
17
)






M
=

(


(



1


1





1




0


0





0




















0


0





0



)



(



0


0





0




1


1





1




















0


0





0



)






(



0


0





0




0


0





0




















1


1





1



)



)





(
18
)







be a sparse matrix of the same size as Sr (i.e., nc×kxkync). It can be verified that MMH=kxkyI, where I denotes an identity matrix. EQ. (16) can be rewritten as






M
H
c
r=(Sr)Hcr.  (19).


EQ. (19) indicates that cr is the eigenvector of a generalized eigenvalue equation





(Sr)Hα=λMHα,  (20)


corresponding to the eigenvalue 1, where α is an nc×1 column vector representing the generalized eigenvector and λ is a scalar representing the generalized eigenvalue.


Since Sr and M are nc×kxkync non-square matrices, it is unclear (1) whether there exist a real eigenvalue λ and (complex) eigenvector α that satisfy EQS. (19) and (20); and (2) how one can solve EQ. (19) efficiently and robustly. The above situation is further complicated by the fact that EQ. (12) usually only holds approximately as τ is usually set less than or equal to the rank of A, even for the calibration data. Therefore, the left and right sides of EQS. (13), (14), (16), and (19) may only be approximately equal


To deal with the aforementioned issues, embodiments of the invention solve the following two eigenvalue equations:














M


(

S
r

)


H



k
x



k
y




β

=


λ
_


β


,




(
21
)









(



S
r



(

S
r

)


H

)


γ

=


λ
~


γ


,




(
22
)







where λ and λ denote the eigenvalues, β and γ denote the eigenvectors. Specifically, embodiments of the invention can compute cr as the eigenvector of EQ. (21) corresponding to the largest eigenvalue. In other words, cr is the solution to the following equation










c
r

=



arg





max




α


:


||
α


||
2


=
1









(

S
r

)

H


α

,


M
H


α









(
23
)







and cr is the optimizer that maximizes the correlation between (Sr)Hα and MHα, i.e., the correlation between the left and right hand sides of EQ. (19). Embodiments of the invention may refer to this approach as MACO, which stands for maximal correlation. In addition, EQ. (23) also applies to the case in which EQ. (19) holds exactly.


It will now be shown that M(Sr)H is Hermitian:


Theorem 1: The matrix M(Sr)H is Hermitian, i.e., M(Sr)H=SrMH. Thus, the eigenvalues of EQ. (21) are all real.


The proof of Theorem 1 is provided in the Appendix, below. The eigen-decomposition of








M


(

S
r

)


H



k
x



k
y






has the following features: (1) the matrix is Hermitian, and the existence of real eigenvalues is guaranteed; (2) there exists numerically stable and efficient approaches for the eigenvalue decomposition of a Hermitian matrix; and (3) the size of M(Sr)H, nc×nc, is much smaller than Sr. In addition, experiments suggest that M(Sr)H may be positive semi-definite.


Next, the relationships among EQS (20), (21), and (22) will be defined.


Theorem 2: For any λ and α satisfying the generalized eigenvalue equation of EQ. (20), λ is real, and














M


(

S
r

)


H



k
x



k
y




α

=
λα

,




(
24
)







i.e., λ and α are also the eigenvalue and eigenvector of EQ, (21), respectively.


The proof of Theorem 2 is given in the Appendix, below. Theorem 2 shows that any λ and α that satisfy EQ. (19) can be computed by applying the eigenvalue decomposition to the square Hermitian matrix









M


(

S
r

)


H



k
x



k
y



.




Theorem 3: For any λ and α satisfying the generalized eigenvalue equation of EQ. (20), then















S
r



(

S
r

)


H



k
x



k
y




α

=


λ
2


α


,




(
25
)







i.e., λ2 and α are also the eigenvalue and eigenvector of EQ. (22), respectively. In addition, λ2 and α are the eigenvalue and eigenvector of EQ. (25), if and only if λ2 and α are the eigenvalue and eigenvector of













(



M


(

S
r

)


H



k
x



k
y



)

2


α

=


λ
2


α


,




(
26
)







The proof of Theorem 3 is given in the Appendix, below. This theorem shows that: (1) if a generalized eigenvalue equation has solutions, it can also be solved by EQ. (22); (2) the eigenvalue decomposition of









S
r



(

S
r

)


H



k
x



k
y






has the same eigenvalues and eigenvectors as








(



M


(

S
r

)


H



k
x



k
y



)

2

;




and (3) the eigenvector decomposition of









S
r



(

S
r

)


H



k
x



k
y






is also the eigenvector of









M


(

S
r

)


H



k
x



k
y



,




as








M


(

S
r

)


H



k
x



k
y






is Hermitian. If the conjecture that








M


(

S
r

)


H



k
x



k
y






is positive semidefinite is true, then









S
r



(

S
r

)


H



k
x



k
y






has the same eigenvalues as









M


(

S
r

)


H



k
x



k
y



.




A naive implementation would be to compute Sr for all the spatial locations and then apply eigenvalue decomposition to EQ. (21) to obtain cr. However, the computational cost is nc2kxky 2D Fourier transformations, leading to a computation cost of O(ne2kxkynxny log(nxny)) and a storage cost of O(nc2kxkynxny). This can be too demanding for both computation and storage. For example, when nx=ny=256, kx=ky=6, and nc=32, it costs about 10 GB for storage using single precision in Matlab and O(109) floating operations for computation. When the size kernel size ky×k3, is large, the computational and storage costs are even higher.


To deal with the above issues, a MACO approach according to embodiments of the invention need only to compute SrMH or Sr(Sr)H. A useful operation for efficiently computing SrMH is











m

i
,
j

r

=




t
=
1



k
x



k
y





s

i
,
j
,
t

r



,


r

,




(
27
)







where mi,jr is the i,j elements of the 2D MR image at spatial location r, which can be rewritten as











m

i
,
j


=




t
=
1



k
x



k
y





s

i
,
j
,
t




,




(
28
)







where mi,j is an nxny×1 column vector containing the entries of mi,jr at all spatial locations. Incorporating EQ. (15):











m

i
,
j


=





t
=
1



k
x



k
y






F
H



(


P
t



(

p

i
,
j
,
t


)


)



=


F
H

(




t
=
1



k
x



k
y






P
t



(

p

i
,
j
,
t


)



)



,




(
29
)







which shows that only one Fourier transformation is needed to obtain ni,j. As i and j goes from 1 to nc and SrMH is Hermitian,








n
c



(


n
c

+
1

)


2




Fourier transformation are needed. With mi,j, it can be shown that the i-th row and j-th column of the matrix SrMH is mi,jr. The computation cost is O(nc2nxny log(nxny)), and the memory cost is O(ne2nxny).


When τ is small, EQ (22) can be more efficiently computed than EQ. (21). Let






z
j,k
=F
H(Pc(vj,k))  (30)


and Zj,kr denotes the r-th entry of zi,k. Let Zr be an nc×τ matrix with the j-row and k-th column being zj,kr, i.e.,










Z
r

=

(




z

1
,
1

r




z

1
,
1

r









z

1
,
τ

r






z

2
,
1

r




z

2
,
2

r









z

2
,
τ

r




























z


n
c

,
1

r




z


n
c

,
2

r









z


n
c

,
τ

r




)





(
31
)







Theorem 4: Sr(Sr)H can be computed as:






S
r(Sr)H=kxkyZr(Zr)H  (32)


The proof of Theorem 4 can be found in the Appendix, below.


By transforming the computation of Sr(Sr)H to Zr(Zr)H in EQ. (32), the number of 2-D Fourier transforms can be reduced from nc2kxky ncτ. As τ is typically far less than nckxky, this can save the computation cost. In addition, the storage cost can be reduced from O(nc2kxkynxny) to O(nc2kxky), when Zr(Zr)H is implemented as follows. An implementation according to an embodiment of the invention is based on the observation that the i-th row and j-th column of Zr(Zr)H can be computed as












k
=
1

τ




z

i
,
k

r




z

i
,
k

r

_






(
33
)







where zi,kr denotes the i-th row and k-th column in Zr. This shows that, instead of pre-computing Zr, embodiments of the invention can compute one column of Zr each time.


Table II, shown in FIG. 5, summarizes the computation cost for computing Sr, M(Sr)H, and Sr(Sr)H. It can be seen that the computational and storage costs of computing M(Sr)H and Sr(Sr)H are much slower than for Sr. In addition, when







τ
<



n
c

+
1

2


,




the computation cost of Sr(Sr)H is lower than M(Sr)H.



FIG. 6 is a flowchart of an eigenvector approach to estimating coil sensitivity maps, according to an embodiment of the invention. Referring now to the figure, an eigenvector approach begins at step 61 by providing the matrix A of sliding blocks of a 2D nx×ny image of coil calibration data as defined in EQ. 7. At step 62, the left singular matrix V of EQ. (9) is calculated from the singular value decomposition of A of EQ. (8), and matrix P is calculated at step 63 from V using P=VVH. At step 64, the elements si,j,tr of matrix Sr are calculated by applying an inverse Fourier transform to a zero-padded P, and at step 65, the eigenvalue equation MHcr=(Sr)Hcr of EQ. (19) is solved using the eigenvalue systems of EQS. (21) and (22). The system of EQ. (21) can be solved using EQS. (27), (28), and (29) at step 66, and the system of EQ. (22) can be solved using EQ. (32), at step 67.


EXPERIMENTS

Simulation:


Experiments were performed on a simulated phantom. A 2D brain phantom and coil sensitivity maps of eight coil channels, were generated using publicly available codes, and the size of the brain phantom is set to nx=ny=256. The simulated eight coil sensitivity maps are shown in FIG. 7(a). With the simulated brain phantom and the simulated sensitivity maps, k-space data was generated for the eight coils. To estimate the coil sensitivity maps, the center 32×32 k-space data was used as the calibration data. The estimated CSM's of a MACO approach according to an embodiment of the invention are shown in FIG. 7(b), where the kernel size is set to kx=ky=6 and τ=55. For comparison, FIG. 7(c) depicts CSM's estimated by the B1-Map approach with spatial scaling, and FIG. 7(d) depicts CSM's obtained by dividing the low resolution coil images by the sum-square images. It can be observed that, a MACO approach according to an embodiment of the invention generates CSM's that are very close to the ground truth. Further, the estimated sensitivity maps in FIG. 7(d) in fact correspond to: (1) those estimated by the B1-Map approach if the spatial smoothing is not applied; and (2) those estimated by a MACO approach according to an embodiment of the invention when the calibration kernel size is set to the size of the calibration data.


Next, CSM's estimated by a MACO approach according to an embodiment of the invention were quantitatively compared with the ground truth and the B1-Map approach. Let c*r be the ground truth coil sensitivities for all the coils at the spatial location r, assuming that c*r has a unique length, i.e., ∥c*r∥=1. Note that, in both a MACO approach according to an embodiment of the invention and a B1-Map, if cr is the obtained eigenvector, then ηcr is also a valid eigenvector for the corresponding eigenvalue if η is a complex scalar with unit length (i.e., |η|=1). The similarity between the estimated cr and the ground truth c*r was measured using the absolute correlation defined as:





ρ(cr,*cr)=custom-characterc*r,crcustom-character|,  (34)


As ∥c*r∥=∥cr∥=1, if ρ(cr,c*r)=1, then cr=ηc*r, where η is a complex scalar with unit length. It is clear that 0≦ρ(cr,c*r)≦1, and the larger the absolute correlation is, the closer is the estimated cr to the ground truth. With the absolute correlation, the goodness of c′ r=1, 2, . . . , nxny can be measured using the mean absolute correlation (mac) defined as










m





a





c

=




r
=
1



n
x



n
y







ρ


(


c
r

,

c
*
r


)




n
x



n
y



.






(
35
)







The kernel size is set to kx=ky=6 and the values of the mean absolute correlation (mac) were explored using different numbers of kept singular vectors, with the results shown in FIGS. 8(a)-(b). In each of FIGS. 8(a), (b), and (c), plot 81 represents mac values as calculated using a MACO approach according to an embodiment of the invention, plot 82 represents mac values as calculated using a B1-Map, and plot 83 represents mac values as calculated by dividing the low resolution coil image by the sum-of-square image. Since A in EQ. (7) is a 288×729 matrix, τ can be selected in the interval [1, 288]. As can be observed from FIG. 8(a), the value of mac is low, when τ is either very small or very large. This is reasonable, because on one hand, when τ is very small, V cannot capture the full range space of the sliding blocks, and thus the left and right hand sides of EQ. (9) do not match well. On the other hand, when τ is very large, P=VVH defined in EQ. (10) tends to an identity matrix, and in this case, the left and right hand sides of EQ. (9) match very well. However, a near identity matrix for P gives too much freedom to the sliding blocks in EQ. (12), and thus P contains less information about the calibration data. In the extreme case, when P is an identity matrix, EQ. (12) always holds and it does not have any information about the calibration data. Nevertheless, FIGS. 8(a)-(b) indicate that one can find a large range of is τ's with which a MACO approach according to an embodiment of the invention yields larger values of mac than both the B1-Map approach and the approach which divides the low resolution coil image by the sum-of-square image. The optimal value of τ may be related to the noise in the calibration data (note that, for the simulated phantom, no noise is added). Experiments on dozens of data sets show that setting τ as the smallest value such that Σk=1τσkk=1288σk≧0.95 works well, where σk is the k-th largest singular vector of A in EQ. (7).


In addition, the performance of a MACO approach according to an embodiment of the invention was tested under varying kernel sizes, and the results are depicted in FIG. 8(c). As can be observed from the figure, the kernel size can change from 2×2 to 10×10 with kx=ky. The value of τ was set to the smallest value for which Σk=1τσkk=1288σk≧0.95. The results again show that a MACO approach according to an embodiment of the invention can achieve a higher mac than the B1-Map. The optimal kernel size may depend on the size of the calibration data, the size of the image, and the noise contained in the calibration data.


Reconstruction For MR reconstruction, the following formulation was used:












min
m




1
2






i
=
1


n
c









F
u



(


c
i


m

)


-
yi



2
2




+

λ




Wm


1



,




(
36
)







where Fu is an under-sampling Fourier transformation operator, ci denotes the coil sensitivity map for the i-th coil, m denotes the image to be reconstructed, y, is the observed undersampled k-space data for the i-th coil, and W is a redundant Haar wavelet transformation.


Reconstruction of Two-dimensional Cardiac Image:


Data was acquired from a healthy volunteer on a 1.5T clinical MR scanner. The imaging parameters include: field of view=313×370 mm2, matrix size=176×208, and 20 coils.



FIG. 9 illustrates the reconstruction of a 2D Cardiac Image. Images in the first row are the coil sensitivity maps estimated by a MACO approach according to an embodiment of the invention for coils 1, 10, and 20, respectively; images in the second row are the coil sensitivity maps estimated by the B1-Map for coils 1, 10, and 20, respectively; the first image in the last row illustrates incoherent Cartesian sampling of the k-space; the second image and the third image in the last row correspond to the reconstructions with the CSM's estimated a MACO approach according to an embodiment of the invention and the B1-Map, respectively.


The first two rows of FIG. 9 depicts the estimated coil sensitivity maps by a MACO approach according to an embodiment of the invention and B1-Map with the center 24 lines as calibration data. It can be seen that the CSMs estimated by a MACO approach according to an embodiment of the invention, shown in the first row, are smoother than the B1-Map, shown in the second row. The first image in the last row illustrates the incoherent Cartesian sampling pattern for generating the under-sampled k-space data. Such incoherent sampling leads to an acceleration factor of 6.3. Images obtained by solving EQ. (36) with CSM's estimated by a MACO approach according to an embodiment of the invention, the second image in the last row, are smoother than those estimated by the B1-Map, shown in the third image in the last row. The CSM estimated by the B1-Map leads to several artifacts, while the image reconstructed with the CSM's estimated by a MACO approach according to an embodiment of the invention contains fewer artifacts due to a more accurate estimation of the sensitivity maps.


Reconstruction of Three-Dimensional Abdomen Image:


Data was acquired from a healthy volunteer on a 1.5T clinical MR scanner. Imaging parameters include: field of view=380×320 mm2, matrix size=320×270, 72 partitions, and 24 coils. A Partial Fourier transform is used in the acquisition of the k-space data, and the sampling pattern in the phase encoding partition plane is shown in FIG. 10(a). Such sampling pattern leads to an acceleration factor of 9. An inverse Fourier transformation is applied to decouple the samples in the readout direction, and the CSM's are estimated and reconstructed for each decoupled slice along this direction. Specifically, the CSM's for each slice are estimated from the center calibration data of size 24×24. Images in the first and second rows of FIG. 10(b) show the coil sensitivity maps estimated by a MACO approach according to an embodiment of the invention and the B1-Map for one partition, and the reconstructed images for the same partition are shown in the last row of FIG. 10(b). Note that: (1) the CSM's obtained by a MACO approach according to an embodiment of the invention are smoother than those by obtained from a B1-Map; (2) some infolding artifacts are observed in the image reconstructed with the CSM's estimated by the B1-Map; and (3) the CSM's estimated by a MACO approach according to an embodiment of the invention have fewer artifacts.


Like the B1-Map approach, if cr is the obtained eigenvector by a MACO approach according to an embodiment of the invention, then ηcr is also a valid eigenvector for the corresponding eigenvalue if η is a complex scalar with unit length, i.e., |η|=1. This indicates that a MACO approach according to an embodiment of the invention may not be able to recover the phase of the reconstructed image. According to embodiments of the invention, coil sensitivity maps were normalized to the first coil, and thus the phase of the reconstructed image is relative to the first coil. This may not be an issue when the magnitude of the reconstructed image is the focus.


Influence of the Amount of Calibration Data on the Estimated CSM Quality

In dynamic cardiac MR imaging, the k-space data can be acquired in an interleaving way, as illustrated in the left plot of FIG. 11, which illustrates a typical Cartesian sampling pattern used in dynamic cardiac MR imaging. In the figure, solid lines represent acquired data, and dashed lines represent missing data. The average k-space data over time, shown in the right plot of FIG. 11, can be used for calibration purpose. One might think that using more the calibration data will improve the quality of the estimated CSM. However, in real applications, the data are usually accompanied by noise. More importantly, the boundary high frequency data might be more sensitive to noise than the center low frequency data, as the high frequencies are typically much weaker in strength than the low frequencies, leading to less robustness when the same noise is added. Therefore, one could conjecture that “the center square calibration data could perform better than the center rectangular calibration data for coil sensitivity maps estimation”.


Experiments were conducted using k-space data acquired on a 1.5T clinical MR scanner. The acquired k-space data was subsampled to generate the following under-sampled k-space data: 24 temporal phases each containing 18 frequency encodings, 30 coils, and image size=192×164. FIGS. 12(a)-(b) illustrate the influence of the amount of used calibration data on the reconstruction. FIG. 12(a) depicts the mask for generating the calibration k-space data and the estimated coil sensitivity maps, where white denotes that the corresponding average k-space data are used for calibration, and FIG. 12(b) depicts the reconstruction images. FIG. 12(a) illustrates four coil sensitivity maps estimated by an eigenvector approach according to an embodiment of the invention, using (1) the center rectangular 44 lines (top row), (2) the center rectangular 32 lines (2nd row), (3) the center 44×44 square block (3rd row), and (4) the center 32×32 square block (bottom row). It can be seen that, compared with the more rectangular calibration data, the square center calibration data can achieve a smoother CSM. For example, the black holes 121, 122 in coil 3 and coil 4 are observed when using center rectangular calibration data, but not when using center square calibration data. FIG. 12(b) compares the reconstruction results achieved by a reconstruction approach according to an embodiment of the invention. In FIG. 12(b), the boxed regions 124 at the upper left are zoomed out and displayed in the rest of each image. The noise observed using a CSM with rectangular calibration data disappears when using the square calibration data.


In addition, experiments were conducted on synthetic static 2D data as follows. A phantom of size 256×256 was generated, and then 8 coil sensitivity maps were synthesized from which the k-space data was generated. The CSM was estimated by using (1) the center kxky square data, and (2) the center rectangular k lines, where k ranges from 32 to 64. The mac (maximum correlation) criterion defined above was computed. The value of mac is between 0 and 1, and a larger mac value means that the estimated coil profile is closer to the simulated ground truth. FIG. 13 compares the square calibration data 131 and the rectangular calibration data 132 in terms of the mac criterion, from which can be seen a clear advantage of the square calibration data. The underlying reason might be that the center low frequency data are less sensitive to the noise than the boundary high frequency data due to the larger intensities.


Extension to 3-Dimensions

When the data is acquired using 3D Cartesian sampling, one way for estimating the coil profiles is to decouple the data along the frequency-encoding direction and then perform the estimation in a 2D manner. Then, the 3D coil profiles can be obtained by stacking the estimated 2D coil profiles. The independent estimation of the coil profiles for each slice causes discontinuity in the profile along the frequency encoding direction. This discontinuity can be propagated into the reconstruction, causing artifacts. In addition, applying a 2D eigen-vector approach ignores the correlation between the calibration data in the frequency encoding direction and downgrades the effect of rank reduction in the eigen-vector approach achieved by singular value decomposition.


A 3D eigen-vector approach according to an embodiment of the invention utilizes all the calibration data together and estimates the coil profiles from different slices simultaneously. FIG. 14 illustrates the conversion of the calibration data into a calibration matrix for 3D CSM estimation, according to an embodiment of the invention. As shown in FIG. 14, the calibration matrix A is constructed from calibration data of size cx×cy×cz by gathering sliding blocks of kernel size kx×ky×kz. The number of columns in A is [(cx−kx+1)×(cy−ky+1)×(cz−kz+1)], and the number of rows of A is kxkykznc, where nc is the number of coils. A low rank subspace of A can be computer by singular value decomposition A=VΣUH, and setting V to the left singular vectors of A corresponding to the leading τ singular values. Let Sc be the coil sensitivity at coil cε[1, 2, . . . , c], which are computed following an eigenvector approach according to an embodiment of the invention as discussed above, i.e., set to the eigenvectors corresponding to the largest eigenvalues.


A 3D CSM estimation method according to an embodiment of the invention was compared with a 2D method using 3D data. The data were acquired from a healthy volunteer on a 1.5T clinical MR scanner. The imaging parameters include a field of view of 320×260×96 mm2, and 26 coils. In the phase encoding-partition plane, a total of 2370 (out of 24960) points were acquired, and the central 24×24 block was treated as calibration data. FIGS. 15(a)-(h) illustrate a comparison of a 2D eigenvector approach and a 3D eigenvector approach according to an embodiment of the invention. FIGS. 15(a), (c), and (e) depict coil profiles of coils nos. 1, 5 and 10 using a 2D eigenvector approach, and FIGS. 15(b), (d), and (f) depict coil profiles of coils nos. 1, 5 and 10 using a 3D eigenvector approach according to an embodiment of the invention. As shown in FIGS. 15(a)-(f), the vertical direction is the frequency encoding direction. The coil profiles from the 3D method are smoother than the ones from the 2D method, especially along the frequency encoding direction. FIGS. 15(g)-(h) compare the central region of the reconstructed image at slice 72 along the transversal axis, using a 2D eigenvector approach for FIG. 15(g) and a 3D eigenvector approach for FIG. 15(h). The overall image using coil profiles from a 3D method according to an embodiment of the invention is smoother than the image from the 2D method. The artifact marked by the circle 152 in FIG. 15(h) is less obvious than the circle 151 in the image of FIG. 15(g).


Reconstruction with the coil profile from a 3D method according to an embodiment of the invention substantially eliminates artifacts caused by discontinuities along the frequency encoding direction from the 2D method. The overall image noise is lower in a 3D method according to an embodiment of the invention due to a smoother coil profile. Efficient computation of a 3D eigenvector method according to an embodiment of the invention can be achieved by utilizing methods according to embodiments of the invention disclosed above, which has a storage cost of O(kxnync2).


System Implementations

It is to be understood that embodiments of the present invention can be implemented in various forms of hardware, software, firmware, special purpose processes, or a combination thereof. In one embodiment, the present invention can be implemented in software as an application program tangible embodied on a computer readable program storage device. The application program can be uploaded to, and executed by, a machine comprising any suitable architecture.



FIG. 16 is a block diagram of an exemplary computer system for implementing an eigenvector approach to coil sensitivity maps (CSM) estimation for 2-D MR images according to an embodiment of the invention. Referring now to FIG. 16, a computer system 161 for implementing the present invention can comprise, inter alia, a central processing unit (CPU) 162, a memory 163 and an input/output (I/O) interface 164. The computer system 161 is generally coupled through the I/O interface 164 to a display 165 and various input devices 166 such as a mouse and a keyboard. The support circuits can include circuits such as cache, power supplies, clock circuits, and a communication bus. The memory 163 can include random access memory (RAM), read only memory (ROM), disk drive, tape drive, etc., or a combinations thereof. The present invention can be implemented as a routine 167 that is stored in memory 163 and executed by the CPU 162 to process the signal from the signal source 168. As such, the computer system 161 is a general purpose computer system that becomes a specific purpose computer system when executing the routine 167 of the present invention.


The computer system 161 also includes an operating system and micro instruction code. The various processes and functions described herein can either be part of the micro instruction code or part of the application program (or combination thereof) which is executed via the operating system. In addition, various other peripheral devices can be connected to the computer platform such as an additional data storage device and a printing device.


It is to be further understood that, because some of the constituent system components and method steps depicted in the accompanying figures can be implemented in software, the actual connections between the systems components (or the process steps) may differ depending upon the manner in which the present invention is programmed. Given the teachings of the present invention provided herein, one of ordinary skill in the related art will be able to contemplate these and similar implementations or configurations of the present invention.


While the present invention has been described in detail with reference to exemplary embodiments, those skilled in the art will appreciate that various modifications and substitutions can be made thereto without departing from the spirit and scope of the invention as set forth in the appended claims.


APPENDIX
Proof of Theorem 1:

The proof begins with two technical lemmas.


Lemma 1:

Assume nx and ny are even numbers. Let A, B, P and Q denote nx×ny complex matrices, with the entry at the i-th-row and the j-th column being ai,j, bi,j, pi,j and qi,j, respectively. Let P and Q be the 2-D inverse Fourier transformations of A and B:











a

m
,
n


=




u
=
1


n
x







υ
=
1


n
y





a

u
,
υ





exp


-
2


π






j


(



(

u
-
1

)



m
/

n
x



+


(

υ
-
1

)



n
/

n
y




)




/



n
x



n
y








,




(
37
)








q

m
,
n


=




u
=
1


n
x







υ
=
1


n
y





b

u
,
υ





exp


-
2


π






j


(



(

u
-
1

)



m
/

n
x



+


(

υ
-
1

)



n
/

n
y




)




/



n
x



n
y








,






where





j

=



-
1


.




If






(
38
)









a

u
,
υ


=


b



n
x

+
2
-
u

,


n
y

+
2
-
υ



_


,

u
=
2

,





,

n
x

,

υ
=
2

,





,

n
y







and




(
39
)








q

1
,
υ


=


a

u
,
1


=


b

1
,
υ


=


b

u
,
1


=
0




,

u
=
1

,





,

n
x

,

υ
=
1

,





,

n
y

,






one





has






p

m
,
n



=


q

m
,
n


_


,

m
=
1

,





,

n
x

,

n
=
1

,





,


n
y

.





(
40
)







Proof:

Inserting EQS. (39) and (40) into EQ. (37), one obtains EQ. (41):













p

m
,
n


=






u
=
2


n
x







υ
=
2


n
y






b



n
x

+
2
-
u

,


n
y

+
2
-
υ



_




exp


-
2


π






j


(



(

u
-
1

)



m
/

n
x



+


(

υ
-
1

)



n
/

n
y




)




/



n
x



n
y













=







u


=
2


n
x








υ


=
2


n
y






b


u


,

υ




_




exp


-
2


π






j


(



(


n
x

+
1
-

u



)



m
/

n
x



+


(


n
y

+
1
-

υ



)



n
/

n
y




)




/



n
x



n
y













=







u


=
2


n
x








υ


=
2


n
y






b


u


,

υ




_




exp


-
2


π






j


(



(

1
-

u



)



m
/

n
x



+


(

1
-

υ



)



n
/

n
y




)




/



n
x



n
y













=







u


=
2


n
x








υ


=
2


n
y







b


u


,

υ






exp


-
2


π






j


(



(


u


-
1

)



m
/

n
x



+


(


υ


-
1

)



n
/

n
y




)





_

/



n
x



n
y












=







u


=
1


n
x








υ


=
1


n
y







b


u


,

υ






exp


-
2


π






j


(



(


u


-
1

)



m
/

n
x



+


(


υ


-
1

)



n
/

n
y




)





_

/



n
x



n
y












=





q

m
,
n


_

.








(
41
)







Lemma 2:

Let u and v denote kxky×1 column vectors with complex entries. Let ut and vt denote the t-th elements of u and v, respectively. Define










p
=




t
=
1



k
x



k
y







H



(



t



(



u
t

_


v

)


)




,




(
42
)







q
=




t
=
1



k
x



k
y







H



(



t



(



υ
t

_


u

)


)




,




(
43
)







which are nxny×1 column vectors. Then,






p= q.  (44)


Proof:

Let










a
=




k
=
1



k
x



k
y







t



(



u
t

_


v

)




,




(
45
)






b
=




t
=
1



k
x



k
y








t



(



υ
t

_


u

)


.






(
46
)







Let A and B be nx×ny complex matrices that are the matrix version of a and b respectively. Making use of the definition of the zero padding operator, and






a
u,v
=b
n

x

+2−u,n

y

+2−v
,u=2, . . . , nx,v=2, . . . , ny,  (47)





and






a
1,v
=a
u,1
=b
1,v
=b
u,1=0,u=0,u=1, . . . , nx,v=1, . . . , ny.  (48)


Applying Lemma 2, one can verify EQ. (44).


Theorem 1 can now be proved. Denote qj,ir as the j-th row and the i-th column of SrMH. It follows from EQS. (17) and (18) that











q

j
,
i

r

=




t
=
1



k
x



k
y





s

i
,
j
,
t

r



,



r
.






(
49
)







Denote qj,i as an nx×ny column vector with the r-th entry being qj,ir. Then










q

j
,
i


=




t
=
1



k
x



k
y






s

i
,
j
,
t


.






(
50
)







To prove that SrMH is Hermitian, it suffices to verify






q
j,i= qi,j  (51)


It follows from EQS. (9) and (10) that











p

i
,
j
,
t


=




k
=
1

τ





υ

i
,
k

t

_



v

j
,
k





,




(
52
)







where vi,kt denotes the t-th element in vi,k. Therefore,













q

j
,
i


=






t
=
1



k
x



k
y







H



(



t



(

p

i
,
j
,
t


)


)









=






t
=
1



k
x



k
y







H

(



t



(




k
=
1

τ





υ

i
,
k

t

_



v

j
,
k




)










=






k
=
1

τ



(




t
=
1



k
x



k
y








H

(



t



(



υ

i
,
k

t

_



v

j
,
k



)


)

.











(
53
)







By utilizing Lemma 2, one can obtain EQ (51).


Since







M


(

S
r

)


H



k
x



k
y






is Hermitian, its eigenvalues are real. This ends the proof of this theorem.


Proof of Theorem 2:

For any pair of λ and α that satisfy EQ. (20), left multiplying both sides of EQ. (20) by M yields






M(Sr)Hα=λMMHα=λkxkyα  (54)


where MMH=kxkyI was used for establishing the last equality. Therefore, λ and α are also the eigenvalue and eigenvector of EQ. (21), respectively.


Proof of Theorem 3:

For any pair of λ and α that satisfy EQ. (20), left multiplying both sides of EQ. (20) by Sr yields






S
r(Sr)Hα=λSrMHα=λ1kxkyα  (55)


where Theorem 2 establishes the second eqaulity. Therefore, λ and α are also the eigenvalue and eigenvector of EQ. (22), respectively.


For any pair of λ and α that satisfy EQ. (25),






S
r(Sr)HMMHα=(λkxky)2α=SrMHSrMHα  (56)


where MMH=kxkyI and (Sr)HM=MHSr have been used (see Theorem 1). The last equality in EQ. (56) indicates that α is the eigenvector of SrMHSrMH corresponding to eigenvalue (λkxky)2.


Proof of Theorem 4:

The proof begins with a technical lemma, whose proof is omitted here.


Lemma 3: FH(Pt(vj,k))custom-characterFH(Pt(vj′,k′)) is invariant to the spatial location t, i.e., FH(Pt(vj,k))custom-characterFH(Pt(vj′,k′))=FH(Pt′(vj,k))custom-characterFH(Pt′(vj′,k′)).


Now, the theorem may be proved. The j-th row and j′-th column of Sr(Sr)H can be computed as













i
=
1


n
c




m

i
,
j
,

j



r


,




where




(
57
)








m

i
,
j
,

j



r

=




t
=
1



k
x



k
y






s

i
,
j
,
t

r




s

i
,

j


,
t

r






,



r
.






(
58
)







EQ. (58) can be rewritten as the following matrix format:










m

i
,
j
,

j




=




t
=
1



k
x



k
y






s

i
,
j
,
t






s

i
,

j


,
t


_

.







(
59
)







Inserting EQ. (52) into EQ. (15) yields










s

i
,
j
,
t


=



F
H



(


P
t



(

p

i
,
j
,
t


)


)


=




k
=
1

τ





v

i
,
k

t

_




F
H



(


P
t



(

v

j
,
k


)


)









(
60
)







In computing si,j,tcustom-charactersi,j′,t, a useful operation is






F
H(Pt(vj,k))custom-characterFH(Pt(vj′,k′))  (61)


which is identical to






F
H(Pc(vj,k))custom-characterFH(Pc(vj′,k′))  (62)


using Lemma 3.


Let Γj,j′r be a τ×τ matrix whose k-th row and k′-th column are the value of zj,krzj′,k′r. Denote vi,: as the i-th row of V. Then






m
i,j,j′
r
=k
x
k
y

v
i,:
Γj,j′rvi,:H.  (63)


It is interesting to note that
















i
=
1


n
c




m

i
,
j
,

j



r


=




k
x



k
y






i
=
1


n
c






v

i
,
:


_



Γ

j
,

j



r




v

i
,
:

H

_










=




k
x



k
y






i
=
1


n
c




trace


(


Γ

j
,

j



r





v

i
,
:

H



v

i
,
:



_


)










=




k
x



k
y



trace
(


Γ

j
,

j



r






i
=
1


n
c






v

i
,
:

H



v

i
,
:



_



)








=




k
x



k
y



trace
(

Γ

j
,

j



r

)









=




k
x



k
y






k
=
1

τ




z

j
,
k

r




z

j
,
k

r

_





,







(
64
)







where trace(A) denotes the summation of the diagonal entries of a matrix A. Therefore, EQ. (32) holds.

Claims
  • 1. A method for estimating a coil sensitivity map for a magnetic resonance (MR) image, comprising the steps of: providing a matrix A of sliding blocks of a 2D nx×ny image of coil calibration data, wherein
  • 2. A method for estimating a coil sensitivity map for a magnetic resonance (MR) image, comprising the steps of: providing a matrix A of sliding blocks of a 2D nx×ny image of coil calibration data, wherein
  • 3. The method of claim 2, further comprising finding an optimizer cr that maximizes a correlation between (Sr)Hα and MHα wherein said optimizer cr is vector of coil sensitivity maps for all coils at spatial location r in the 2D nx×ny image.
  • 4. The method of claim 2, wherein solving MHcr=(Sr)Hcr comprises solving eigenvalue equations
  • 5. The method of claim 4, further comprising calculating SrMH using
  • 6. The method of claim 4, further comprising calculating Sr(Sr)H as Sr(Sr)H=kxkyZr(Zr)H, wherein
  • 7. The method of claim 6, further comprising computing the i-th row and j-th column of
  • 8. The method of claim 2, wherein the coil calibration data is acquired from a center square of frequency space data.
  • 9. The method of claim 2, wherein the matrix A is constructed from calibration data of size cx×cy×cz by gathering sliding blocks of kernel size kx×ky×kz, wherein A has [(cx−kx+1)×(cy−ky+1)×(cz−kz+1)] columns and kxkykznc rows.
  • 10. A non-transitory program storage device readable by a computer, tangibly embodying a program of instructions executed by the computer to perform the method steps for estimating a coil sensitivity map for a magnetic resonance (MR) image, the method comprising the steps of: providing a matrix A of sliding blocks of a 2D nx×ny image of coil calibration data, wherein
  • 11. The computer readable program storage device of claim 10, the method further comprising finding an optimizer cr that maximizes a correlation between (Sr)Hα and MHα wherein said optimizer cr is vector of coil sensitivity maps for all coils at spatial location r in the 2D nx×ny image.
  • 12. The computer readable program storage device of claim 10, wherein solving MHcr=(Sr)Hcr comprises solving eigenvalue equations
  • 13. The computer readable program storage device of claim 12, the method further comprising calculating SrMH using
  • 14. The computer readable program storage device of claim 12, the method further comprising calculating Sr(Sr)H as Sr(Sr)H=kxkyZr(Zr)H, wherein
  • 15. The computer readable program storage device of claim 14, the method further comprising computing the i-th row and j-th column of Zr(Zr)H as
  • 16. The computer readable program storage device of claim 10, wherein the coil calibration data is acquired from a center square of frequency space data.
  • 17. The computer readable program storage device of claim 10, wherein the matrix A is constructed from calibration data of size cx×cy×cz by gathering sliding blocks of kernel size kx×ky×kz, wherein A has [(cx−kx+1)×(cy−ky+1)×(cz−kz+1)] columns and kxkykznc rows.
CROSS REFERENCE TO RELATED UNITED STATES APPLICATIONS

This application claims priority from “On the Eigen-Vector Approach for Coil Sensitivity Maps Estimation”, U.S. Provisional Application No. 61/618,002 of Liu, et al., filed Mar. 30, 2012, “Coil Sensitivity maps Estimation in Dynamic CMRI: Square Could be Better than Rectangular”, U.S. Provisional Application No. 61/724,023 of Liu, et al., filed Nov. 8, 2012, and “An Eigen-Vector Approach for Coil Sensitivity Estimation in the 3D Scenario”, U.S. Provisional Application No. 61/724,001 of Liu, et al., filed Nov. 8, 2012, the contents of all of which are herein incorporated by reference in their entireties.

Provisional Applications (3)
Number Date Country
61618002 Mar 2012 US
61724023 Nov 2012 US
61724001 Nov 2012 US