The present application claims priority to Chinese patent application No. 201310185159.0, filed on May 17, 2013, and entitled “PARALLEL ACQUISITION IMAGE RECONSTRUCTION METHOD AND DEVICE FOR MAGNETIC RESONANCE IMAGING”, the entire disclosure of which is incorporated herein by reference.
The present disclosure generally relates to Magnetic Resonance Imaging (MRI) technology, and more particularly, to a parallel acquisition image reconstruction method and a parallel acquisition image reconstruction device for magnetic resonance imaging.
Medical MRI is an imaging technology using characteristics of magnetic sensitive nuclei (proton nuclei) of human body in a magnetic field. In MRI technology, imaging speed is an important criterion to evaluate an imaging method. So far, methods for improving MRI speed can be classified into three categories: first, improving performance of hardware; second, improving scanning technology in k-space; and third, partial k-space data imaging. In MRI technology, magnetic resonance signal space (original data space) is known as k-space, which is also known as Fourier transformation space. Signal data sampled in k-space is processed by inverse Fourier transform and magnitude calculation, in order to obtain a magnetic resonance magnitude image. Partial k-space data acquisition and reconstruction is an imaging technology which just samples a part of k-space data, so that the scanning speed can be substantially increased without hardware and other additional changes.
Data sampling and k-space filling are important factors that limit the speed of MRI data acquisition. Techniques of parallel acquisition image reconstruction for magnetic resonance imaging fill partially the k-space data with a coil reconstruction-merger method using multi-channel phased array coils, so as to obtain full k-space data from the partially sampled k-space data for image reconstruction. In this way, only part of the k-space data is actually sampled, leading to increased imaging speed.
Generalized Auto-calibrating Partially Parallel Acquisitions (GRAPPA) is a commonly used parallel acquisition image reconstruction method. A conventional GPAPPA method is shown in
Recently, an Iterative Self-consistent Parallel Imaging Reconstruction from Arbitrary k-Space (SPIRiT) method is provided. The SPIRiT method reconstructs images in an iterative manner, and has a better reconstruction result than the GRAPPA method. As shown in
However, a large number of channels are required in the parallel acquisition process to obtain a better imaging result, but reconstruction speed is reduced when the number of channel is large. Meanwhile, because some channels may have poor signal-to-noise ratios, the signal-to-noise ratio of the final image can be reduced if data of all channels are combined. In a word, the conventional method can result in low reconstruction speed and non-optimal signal-to-noise ratio.
Therefore, a parallel acquisition image reconstruction method for magnetic resonance imaging, which can improve image reconstruction speed and signal-to-noise ratio, is needed.
The present disclosure aims to solve the problems of a low image reconstruction speed and a low signal-to-noise ratio image of the conventional parallel acquisition image reconstruction method for magnetic resonance imaging.
In order to solve the problems mentioned above, a parallel acquisition image reconstruction method for magnetic resonance imaging is provided in the present disclosure. The method includes:
sampling magnetic resonance signals from a plurality of channels, and filling the magnetic resonance signals in an initial k-space, where the initial k-space includes a fully-sampled area and an under-sampled area, each data point in the fully-sampled area is sampled, and the under-sampled area includes sampled data points and un-sampled data points;
obtaining a first virtual space by performing a mathematical transformation on the initial k-space, where the first virtual space comprises a plurality of virtual channels each having a first parameter for evaluating signal-to-noise ratio, and obtaining a second virtual space by reserving virtual channels of the first virtual space which have first parameters higher than a predetermined threshold value;
calculating a first combination coefficient from the data of the initial k-space to data of the second virtual space, and a second combination coefficient from the data of the second virtual space to the data of the initial k-space;
putting the first combination coefficient and the second combination coefficient into a predetermined objective function to obtain the data of the second virtual space; and
transforming the data of the second virtual space to an image domain to obtain a reconstructed image.
In some embodiments, the mathematical transformation includes: Wavelet transformation, Curvelet transformation or Karhunen-Loeve (KL) transformation.
In some embodiments, the Karhunen-Loeve transformation is adopted to obtain the first virtual; data of the fully-sampled area serves as calibration data, and the first parameter is amplitude of an eigenvalue of the KL transformation of the calibration data of each channel.
In some embodiments, the first combination coefficient is a convolution kernel of a fitting operation from the initial k-space to the second virtual space of each channel, where the first combination coefficient is obtained according to an equation shown below:
Src*G=Dst,
where Src represents the data of the initial k-space of each channel, Dst represents the data of the second virtual space, and G represents the first combination coefficient;
the second combination coefficient is a convolution kernel of a fitting operation from the second virtual space to the initial k-space of each channel, where the second combination coefficient is obtained according to an equation shown below:
Dst*P=Src,
where P represents the second combination coefficient.
In some embodiments, the objective function is expressed as an equation shown below:
Val=∥GPx−x∥2+λ·Reg(x)
where G represents the first combination coefficient, P represents the second combination coefficient, x represents the data of the second virtual space, Reg(x) represents a cost function, λ represents a coefficient of the cost function, Val represents a target value, and the data x of the second virtual space is obtained according to the objective function under a condition that the target value Val has a minimum value.
In some embodiments, the objective function is expressed as an equation shown below:
Val=∥(GP−I)(DTa+DcTm)∥2+λ·∥ψFH(DTa+DcTm)∥1
where DTa represents sampled data of the second virtual space, DcTm represents un-sampled data of the second virtual space, Ψ represents a regularization transformation matrix, F represents a Fourier transformation, Val represents a target value, the un-sampled data DcTm of the second virtual space is obtained according to the objective function under a condition that the target value Val has a minimum value, and the data x of the second virtual space is obtained by combining the un-sampled data DcTm and the sampled data DTa of the second virtual space.
In some embodiments, the objective function is expressed as an equation shown below:
Val=∥(GP−I)(DTa+DcTm)∥2+λ·∥ψ·S·FH(DTa+DcTm)∥1
where DTa represents sampled data of the second virtual space, DcTm represents un-sampled data of the second virtual space, Ψ represents a regularization transformation matrix, S represents a coil sensitivity coefficient matrix, F represents a Fourier transformation, Val represents a target value, the un-sampled data DcTm of the second virtual space is obtained according to the objective function under a condition that the target value Val has a minimum value, and the data x of the second virtual space is obtained by combining the un-sampled data DcTm and the sampled data DTa of the second virtual space.
In some embodiments, assuming that each data point of the second virtual space is unknown, the objective function is expressed as an equation shown below:
Val=∥(GP−I)x∥2+λ·∥ψFHx∥1+β·∥Dx−a∥2
where x represents the data of the second virtual space, Ψ represents a regularization transformation matrix, F represents a Fourier transformation, D represents a sampling matrix of sampled data, β represents an adjustment coefficient, Val represents a target value, and the data x of the second virtual space is obtained according to the objective function under a condition that the target value Val has a minimum value.
In some embodiments, assuming that each data point of the second virtual space is unknown, the objective function is expressed as an equation shown below:
Val=∥(GP−I)x∥2+λ·∥ψ·S·FHx∥1+β·∥Dx−a∥2
where x represents the data of the second virtual space, Ψ represents a regularization transformation matrix, S represents a coil sensitivity coefficient matrix, F represents a Fourier transformation, D represents a sampling matrix of sampled data, β represents an adjustment coefficient, Val represents a target value, and the data x of the second virtual space is obtained according to the objective function under a condition that the target value Val has a minimum value.
In some embodiments, if the parallel acquisition image reconstruction has a time dimension, the objective function is expressed as an equation shown below:
Val=Σf=1N
where xf represents data of an fth two dimensional frame of the second virtual space in parallel acquisition image reconstruction having the time dimension.
In one embodiment, a parallel acquisition image reconstruction device for magnetic resonance imaging is provided, the device may includes:
a data sampling unit, adapted for sampling magnetic resonance signals from a plurality of channels, and filling them in an initial k-space;
a data transformation unit, adapted for performing a mathematical transformation on the initial k-space to obtain a first virtual space;
a channel selection unit, adapted for reserving data of channels of the initial k-space, which have a first parameter higher than a predetermined threshold value, to obtain a second virtual space;
a coefficient calculation unit, adapted for calculating a first combination coefficient from the initial k-space to the second virtual space, and a second combination coefficient from the second virtual space to the initial k-space;
a second virtual space calculation unit, adapted for putting the first combination coefficient and the second combination coefficient, which are obtained by the coefficient calculation unit, into a predetermined objective function to obtain data of the second virtual space; and
an image reconstruction unit, adapted for transforming the data of the second virtual space, which is obtained by the second virtual space calculation unit, into an image domain to obtain a reconstructed image.
Compared with the prior art, the present disclosure has the following advantages: obtaining the first virtual space by performing a mathematical transformation on the initial k-space; and reserving data of channels of the first virtual space, which have a higher signal-to-noise ratio, to serve as the second virtual space by calculating the first parameter of the first virtual space which is used to evaluate signal-to-noise ratio of each channel, whereby the image reconstruction speed is improved and the signal-to-noise ratio of image is improved.
In order to make a better understanding of the present disclosure for the persons skilled in the art, a conventional parallel acquisition image reconstruction method for magnetic resonance imaging is briefly described below first.
After filling sampled magnetic resonance signals in k-space, data of k-space can be expressed as Equation (1):
x=D
T
a+D
c
T
m Equation (1)
where x represents a N*1 vector; a represents a A*1 vector; m represents a M*1 vector; D represents a A*N matrix, which is a sampling matrix of sampled data; and Dc represents a M*N matrix, which is a sampling matrix of un-sampled data.
At first, the sampled data DTa of each channel is used to obtain un-sampled data DcTm by fitting. Then data x of the k-space is obtained. The data x of k-space is processed by inverse Fourier transformation and magnitude operation, so as to obtain a magnetic resonance magnitude image. A conventional method for obtaining the un-sampled data DcTm by fitting is described below.
As mentioned above, the GPAPPA method is usually adopted in the prior art. The method includes:
letting: x=GDTa Equation (2)
obtaining: DTa+DcTm=GDTa Equation (3)
obtaining after a transformation: DcTm=(G−I)DTa Equation (4)
where G represents a convolution kernel of the GRAPPA method.
Referring to
A SPIRiT method in the prior art is described as below:
letting: x=Kx Equation (5)
obtaining: (K−I)(DTa+DcTm)=0 Equation (6)
obtaining after a transformation: (K−I)DcTm=−(K−I)DTa Equation (7)
where K represents a convolution kernel of the SPIRiT method. As shown in
In order to clarify the objects, characteristics and advantages of the disclosure, the embodiments of the present disclosure will be described in detail in conjunction with the accompanying drawings.
In S01, magnetic resonance signals from a plurality of channels are sampled and the magnetic resonance signals are filled into an initial k-space, where the initial k-space includes a fully-sampled area and an under-sampled area, each data point in the fully-sampled area is sampled, and the under-sampled area includes sampled data points and un-sampled data points.
It should be noted that, in the embodiment, in order to avoid confusion, a space, in which initial magnetic resonance signals (data) obtained from all channels by directly parallel acquisition are filled, is referred to as the initial k-space. Hereinafter, a space, which is obtained by a mathematical transformation on the initial k-space, is referred to as a first virtual space, and a space, which reserves channels of the first virtual space which have a first parameter higher than a predetermined threshold value, is referred to as a second virtual space, where the first parameter is used to evaluate signal-to-noise ratio of each channel.
In one embodiment, the initial k-space x is a [NRONPENC] matrix, where NC represents a channel number, and NRO and NPE represent a dimension size of a frequency encoding direction and a dimension size of a phase encoding direction, respectively. In one embodiment, data of the fully-sampled area can be used as calibration data. Therefore, a calibration data matrix may be represented by [NrRONrPENC], where r represents a calibration data part of the initial k-space.
In S02, a first virtual space is obtained by performing a mathematical transformation on the initial k-space, where the first virtual space comprises a plurality of virtual channels each having a first parameter for evaluating signal-to-noise ratio, and a second virtual space is obtained by reserving virtual channels of the first virtual space which have first parameters higher than a predetermined threshold value.
In one embodiment, a mathematical transformation method, such as Wavelet transformation, Curvelet transformation, Karhunen-Loeve (KL) transformation and etc, may be used obtain a first virtual space. Other embodiments are possible, which is not limited herein. The KL transformation is taken as an example and described in detail below.
The data of the initial k-space is transformed by a KL transformation to obtain the first virtual space. Amplitude of an eigenvalue of the KL transformation of the calibration data of the initial k-space of each channel serves as a first parameter value.
KL transformation is performed on the calibration data part r of the initial k-space to obtain:
r′=r×K Equation (8)
where K represents a NC×NC KL transformation matrix, r′ represents data of the first virtual space after the KL transformation. After the KL transformation, the amplitude of the eigenvalue of K can be used to evaluate a signal-to-noise ratio of a channel of the first virtual space obtained after the mathematical transformation. The higher the eigenvalue, the higher the signal-to-noise ratio. Therefore, after the mathematical transformation, channels having a low signal-to-noise ratio may be eliminated, in order to improve the quality of the reconstructed image.
As shown in
In S03, a first combination coefficient from the initial data k-space to the second virtual space is calculated, and a second combination coefficient from the second virtual space to the initial k-space is calculated.
In one embodiment, a process of filling the second virtual space is shown in
Src*G=Dst Equation (9)
where Src represents the calibration data of the initial k-space, Dst represents the calibration data of the second virtual space, and G represents a convolution kernel from Src to Dst, which is also the first combination coefficient.
Similarly, a fitting operation may also be performed from Dst to Src, which can be expressed by Equation (10):
Dst*P=Src Equation (10)
where P represents a convolution kernel from Dst to Src, which is also the second combination coefficient.
In S04, the first combination coefficient and the second combination coefficient are put into a predetermined objective function to obtain the data of the second virtual space.
In one embodiment, the predetermined objective function may be expressed as Equation (11):
Val=∥GPx−x∥2+λ·Reg(x) Equation (11)
where G represents the first combination coefficient, P represents the second combination coefficient, x represents the data of the second virtual space, Reg(x) represents a cost function, λ represents a coefficient of the cost function, and Val represents a target value.
The data x of the second virtual space may be obtained according to Equation (11) under a condition that the target value Val has a minimum value.
In one embodiment, different calculation methods may be adopted in different situations. For example, if the sampled data of the second virtual space is not changed, an objective function (11-1) shown as below may be used to obtain the data x of the second virtual space:
Val=∥(GP−I)(DTa+DcTm)∥2+λ·∥ψFH(DTa+DcTm)∥1 Equation (11-1)
where DTa represents sampled data of the second virtual space, DcTm represents un-sampled data of the second virtual space, Ψ represents a regularization transformation matrix, F represents a Fourier transformation, and Val represents a target value. The un-sampled data DcTm of the second virtual space may be obtained according to Equation (11-1) under a condition that the target value Val has a minimum value, and the data x of the second virtual space may be obtained by combining the un-sampled data DcTm and the sampled data DTa of the second virtual space.
When Equation (11-1) is adopted, the regularization transformation matrix needs to be applied to each channel of the second virtual space, resulting in a large time cost. In order to improve image reconstruction speed, a coil sensitivity coefficient is introduced in the calculation. An objective function (11-2) shown as below may be used:
Val=∥(GP−I)(DTa+DcTm)∥2+λ·∥ψ·S·FH(DTa+DcTm)∥1 Equation (11-2)
where DTa represents sampled data of the second virtual space, DcTm represents un-sampled data of the second virtual space, Ψ represents a regularization transformation matrix, S represents a coil sensitivity coefficient matrix, F represents a Fourier transformation, and Val represents a target value. The un-sampled data DcTm of the second virtual space may be obtained according to Equation (11-2) under a condition that the target value Val has a minimum value, and the data x of the second virtual space may be obtained by combining the un-sampled data DcTm and the sampled data DTa of the second virtual space.
In one embodiment, a method for obtaining the coil sensitivity coefficient matrix S may include: setting the data of the under-sampled area of the initial k-space to zero (only retaining the calibration data); performing a transformation operation to obtain image domain data Image_i of the multiple channels; combining (e.g., summing, summing of squares, or others) data of the multiple channels to obtain Image—0; and calculating Si=Image_i./Image—0, where “./” represents a division method in a point by point manner, Image_i represents an image of an ith channel, Image—0 represents a total image, and Si represents a coil sensitivity coefficient of the ith channel.
The above is only an example to obtain the coil sensitivity coefficient matrix S. In some embodiments, the coil sensitivity coefficient matrix S may be obtained with other methods.
After the introduction of the coil sensitivity coefficient matrix S, the regularization transformation matrix Ψ only needs to be applied once on an image after channel combination, thereby the image reconstruction speed is improved.
In S05, transforming the data of the second virtual space is transformed to an image domain to obtain a reconstructed image.
In one embodiment, inverse Fourier transformation and magnitude operation are performed on the data of the second virtual space, in order to obtain a magnetic resonance magnitude image. In some embodiments, the step of magnitude operation is optional.
According to embodiments of the present disclosure, the data of the initial k-space is processed by mathematical transformation to serve as the first virtual space, and the channels having high signal-to-noise ratios are reserved to serve as the second virtual space which is used for image reconstruction, in which the number of channels used for image reconstruction is reduced. Therefore, the image reconstruction speed is improved, and the image quality is improved.
In addition, in the process of obtaining the data x of the second virtual space, the regularization transformation matrix Ψ only applied once on the image after channel combination by introducing the coil sensitivity coefficient matrix S. Therefore, time cost is reduced, and the image quality is further improved.
A device corresponding to the method mentioned above for parallel acquisition image reconstruction is provided.
a data sampling unit 701, adapted for sampling magnetic resonance signals from a plurality of channels, and filling them in an initial k-space;
a data transformation unit 702, adapted for performing a mathematical transformation on the initial k-space to obtain a first virtual space;
a channel selection unit 703, adapted for reserving data of channels of the initial k-space, which have a first parameter higher than a predetermined threshold value, to obtain a second virtual space;
a coefficient calculation unit 704, adapted for calculating a first combination coefficient from the initial k-space to the second virtual space, and a second combination coefficient from the second virtual space to the initial k-space;
a second virtual space calculation unit 705, adapted for putting the first combination coefficient and the second combination coefficient, which are obtained by the coefficient calculation unit 704, into a predetermined objective function to obtain data of the second virtual space; and
an image reconstruction unit 706, adapted for transforming the data of the second virtual space, which is obtained by the second virtual space calculation unit 705, into an image domain to obtain a reconstructed image.
In one embodiment, a Wavelet transformation, a Curvelet transformation, a Karhunen-Loeve (KL) transformation or other mathematical transformation methods may be used by the data transformation unit 702 to transform the data of the initial k-space to the first virtual space.
In the above embodiments, the sampled data is not changed in the process of calculating the objective function.
In some embodiments, the sampled data may be changed in the process of calculating the objective function, which means each data point of the second virtual space is assumed to be unknown. Then the data x of the second virtual space may be obtained according to an objective function (11-3) shown as below:
Val=∥(GP−I)x∥2+λ·∥ψFHx∥1+β·∥Dx−a∥2 Equation (11-3)
where x represents the data of the second virtual space, Ψ represents a regularization transformation matrix, F represents a Fourier transformation, β represents an adjustment coefficient, D represents a sampling matrix of the sampled data, and Val represents a target value. And the data x of the second virtual space may be obtained according to Equation (11-3) under a condition that the target value Val has a minimum value.
Similarly, a coil sensitivity coefficient may be introduced into Equation (11-3), in order to further improve imaging speed. Specifically, the data x of the second virtual space may be obtained according to an objective function (11-4) shown as below:
Val=∥(GP−I)x∥2+λ·∥ψ·S·FHx∥1+β·∥Dx−a∥2 Equation (11-4)
where x represents the data of the second virtual space, Ψ represents a regularization transformation matrix, S represents a coil sensitivity coefficient matrix, F represents a Fourier transformation, D represents a sampling matrix of sampled data, β represents an adjustment coefficient, and Val represents a target value. The data x of the second virtual space may be obtained according to Equation (11-4) under a condition that the target value Val has a minimum value.
The method according to embodiments of the present disclosure may be used in parallel acquisition image reconstruction in which a time dimension is introduced. For an fth two dimension (2D) frame, data of the second virtual space may be expressed as Equation (12):
x
f
=D
T
a
f
+D
c
T
m
f Equation (12)
where xf represents the fth frame data of the second virtual space.
In the process of reconstructing a 2D frame at time T, all xf is combined to form a great vector according to Equation (13):
∥(GP−I)x∥2=Σf=1N
By putting Equation (12) or Equation (13) to Equation (11) or any one of Equations from (11-1) to (11-4), an objective function for calculating data x of the second virtual space is obtained in parallel acquisition image reconstruction having a time dimension.
For example, by putting Equation (13) into Equation (11), an objective function for calculating the data xf of the second virtual space, which has a time dimension, is obtained as below:
Val=Σf=1N
Embodiments of the present disclosure are described in a progressive way. Differences between the embodiments are described in detail, while similar parts of the hereafter mentioned embodiments may refer to the aforementioned embodiments, for the sake of clarity.
Although the present disclosure has been disclosed above with reference to preferred embodiments thereof, it should be understood that the disclosure is presented by way of example only, and not limitation. Those skilled in the art can modify and vary the embodiments without departing from the spirit and scope of the present disclosure.
Number | Date | Country | Kind |
---|---|---|---|
201310185159.0 | May 2013 | CN | national |