Patent Application
20120008844
This disclosure is directed to methods for reconstructing time-varying signals from limited number of measurements. More specifically, this disclosure is directed to using motion estimation and compensation techniques to reconstruct more accurate estimates of support for image sequences with substantial motion, such as cardiac magnetic resonance imaging (MRI) sequences.
Compressive Sensing or Compressive Sampling (CS) is a recently introduced framework for reconstructing sparse or compressible signals accurately from far fewer samples than suggested by the Nyquist-Shannon sampling theorem. CS theory states that signals with sparse representations can be recovered from a small number of incoherent linear measurements. The CS theory has attracted significant attention for applications such as Magnetic Resonance Imaging (MRI) where long acquisition times have been problematic. This is especially true for applications such as dynamic MRI where high spatio-temporal resolution is needed. For example, in cardiac cine MRI, it is desirable to acquire the whole cardiac volume within a single breath-hold in order to avoid artifacts due to respiratory motion. Conventional MRI techniques do not allow reconstruction of high resolution image sequences from such limited amount of data. However, since MR images often have sparse representations, CS is an ideal candidate for accelerating data acquisition in MRI.
There has been recent interest in incorporating prior information into the CS framework to handle situations with partially known support, referred to as a sparsity pattern. One extension assumes that the support of the signal in the sparse domain is partially known. For example, such prior information can be available when recursively reconstructing a time sequence of sparse signals. Under the assumption that the support of the signal changes slowly over time, the support of the signal at the previous time-point can be used as the known part of the support for the current time-point. If the known part of the support is denoted using the set T, the basic idea is to find a signal with the sparsest support outside of T while satisfying the data consistency constraint. While this “modified CS” (MCS) approach works well for image sequences with little or no motion, motion causes significant change in support between adjacent frames. When the support estimate is not accurate, the performance of MCS deteriorates.
Exemplary embodiments of the invention as described herein generally include methods and systems for applying motion estimation to a sequence of images reconstructed using MCS. Motion estimation and compensation techniques that have been successfully used in video coding applications can be used to improve support estimation, and thus, the overall performance. Using this motion information, motion compensated frames are created. Support estimates are then obtained from these motion compensated frames instead of the initially reconstructed frames and the MCS framework is used with this updated support information to produce the final reconstructed image sequence. Experimental results using phantoms as well as real MRI data sets illustrate the improved performance of a method according to an embodiment of the invention.
According to a first aspect of the invention, there is provided a method for reconstructing a digital image from a set of measurements, the method including providing a previous image frame corresponding to a previous time point in a time series of measurements of an image signal, and a current image frame corresponding to a current time point in the time series, calculating an estimated motion vector dt=[d1(n1, n2, t), d2(n1, n2, t)]T for spatial point (n1, n2) and current time point t between the previous image frame and the current image frame, calculating a motion compensated current image frame ft(n1, n2) from the previous image frame ft−1(n1, n2) as f1(n1, n2)=ft−1(n1+d1(n1, n2, t), n2+d2(n1, n2, t)), estimating a known support set of a sparse signal estimate of the motion compensated current image frame, where the support set of the sparse signal estimate comprises indices of non-zero elements of the sparse signal estimate, calculating a sparse signal corresponding to the current image frame whose support contains a smallest number of new additions to the known support set while satisfying a data consistency constraint, and correcting the motion compensated current image frame image frame from the sparse signal.
According to a further aspect of the invention, the time series of measurements of the image signal comprises a measurement vector yt and a measurement matrix At for each time point t, the sparse signal xt is calculated for each time point t of the time series of measurements from the measurement vector yt and the measurement matrix At, and the image frame ft is calculated from the sparse signal xt for each time point in the time series.
According to a further aspect of the invention, calculating the sparse signal xt from the measurement vector y and the measurement matrix At includes finding an xt that minimizes an ll norm |xt| subject to a data consistency constraint yt=Atxt.
According to a further aspect of the invention, calculating the sparse signal xt from the measurement vector y and the measurement matrix At for each time point t includes finding xt that minimizes |xt| subject to yt=Atxt for t=0, where xt is an N-dimensional vector, yt is an M-dimensional vector, and M<<N, initializing a support set of xt at t=0 as Ŝ0={iε[1,N]:({circumflex over (x)}0)i2>α}, where α is a predetermined threshold, where the support set of xt comprises indices of non-zero elements of xt, finding an {circumflex over (x)}t that minimizes |xt| over the complement of the support set subject to yt=Atxt, updating the support set Ŝt of xt from Ŝt={iε[1,N]:({circumflex over (x)}t)i2>α}, and repeating the steps of finding an {circumflex over (x)}t that minimizes |xt| and updating the support set Ŝt of xt for all t>0.
According to a further aspect of the invention, the image frame for a given time point is related to the sparse signal for that time point by a linear transformation.
According to a further aspect of the invention, estimating the known support set Ŝt of the sparse signal of the motion compensated current image frame comprises calculating Ŝt of xt from Ŝt={iε[1,N]:({circumflex over (x)}tMC)i2>α}, where t is the current time point, α is a predetermined threshold, {circumflex over (x)}tMC is an N-dimensional vector that is the sparse signal of the motion compensated current image frame, and the support set of {circumflex over (x)}tMC comprises indices of non-zero elements of {circumflex over (x)}tMC.
According to a further aspect of the invention, calculating a sparse signal corresponding to the current image frame comprises finding {circumflex over (x)}t that minimizes |xt| over the complement of the support set Ŝt while satisfying the data consistency constraint.
According to a further aspect of the invention, the data consistency constraint is
y
t
=A
t
x
t
According to a further aspect of the invention, finding the sparse signal corresponding to the current time point comprises finding {circumflex over (x)}t that minimizes |Wxt| over a complement of a support set V={iε[1,N]:(xP)i2>α} subject to yt=Atxt, where xP is a N-dimensional vector of predictions of the sparse signal estimate of the motion compensated image frame, and W is a weight matrix that is a function of the prediction vector xP.
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 executable by the computer to perform the method steps for reconstructing a digital image from a set of measurements.
a)-(d) depict an original frame of the phantom sequence and three reconstructed frames, according to an embodiment of the invention.
Exemplary embodiments of the invention as described herein generally include systems and methods for using motion estimation and compensation techniques to reconstruct more accurate estimates of support for image sequences with substantial motion. 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.
Let x denote an N-dimensional vector that represents the desired signal in the sparse domain. This signal can be linearly measured using an M×N measurement matrix A,
y=Ax, (1)
to obtain the M-dimensional measurement vector y, where M<<N. Given the measurement vector y and the measurement matrix A, the CS theory suggests that the sparse vector x can be recovered by solving the constrained ll minimization:
arg
xmin∥x∥l such that y=Ax, (2)
provided that certain conditions on sparsity and restricted isometry are satisfied. In EQ. (2), minimizing the ll norm promotes sparsity while the data constraint y=Ax is used to ensure consistency with measurements. If x is K-sparse, that is, only K out of N coefficients are non-zero, practical results suggest that x can be recovered from 3K to 5K random measurements. The CS theory has been extended to cover signals that are not strictly sparse, but are compressible.
Before describing MCS, notation should be established. Let AC denote the complement of a set A, A\ B=A∩BC denote set difference, and |A| denote the cardinality of the set A. Let (β)A denote a sub-vector that contains the elements with indices in A.
MCS makes use of partial information about the support of the sparse vector x that is available from prior knowledge. Let S denote the support set of x, i.e. the set S contains the indices of the non-zero elements of x. This support can be expressed as S=T∩Δ\Δe, where T denotes the known part of the support of x, Δ=S\T denotes the unknown part of the support, and Ae=T\ S denotes the error in the known part of the support. The goal of MCS is to find a signal β whose support contains the smallest number of new additions to the known support T while satisfying the data consistency constraint:
arg
βmin∥(β)T
Similar to the CS, minimizing the ll norm of (β)T
An MCS method can be used to reconstruct a time sequence of images. Let xt denote the sparse representation of the image frame at time instance t and yt denote the measurement vector at time instance t. Similarly, let St and At denote the support set of xt and the measurement matrix at time instance t, respectively. An algorithm for dynamic MCS can then be stated as follows, with reference to the flow chart of
A dynamic method begins at step 11, for time t=0, by solving the MCS minimization argβmin∥(β)T
As described above, a dynamic MCS framework relies on the assumption that the support of the signal changes slowly over time, and, thus, the estimate of the support set of xt obtained using the previous time point is accurate. In other words, dynamic MCS works well when |St\Ŝt−1| is small. However, in many image sequences, the objects being imaged move during data acquisition and such motion can cause significant differences in support sets of adjacent time points. If the motion between image frames can be estimated, the support estimate can be updated to compensate for such motion. According to an embodiment of the invention, motion compensation methods can be incorporated into the MCS framework.
While a framework according to an embodiment of the invention can easily accommodate more sophisticated motion models, a simple motion model can be used in this work, by only considering image-plane motion and assume that it can be well approximated by 2-dimensional (2-D) correspondence fields. Let ƒ(n1, n2, t) denote the spatio-temporal sample at time t and spatial location (n1, n2). If one assumes that there exists real-valued correspondence vectors dt=[d1(n1, n2, t), d2(n1, n2, t)]T such that
ƒ(n1,n2,t+1)=ƒ(n1+d1(n1,n2,t),n2+d2(n1,n2,t). (4)
According to an embodiment of the invention, the correspondence vector can be determined by a non-rigid registration of the two images. Such registration can be computationally complex, and according to other embodiments of the invention, it is sufficient to determine the correspondence vector using simple block matching, under the assumption of translational motion. Using ft to denote the image frame at time t in vector notation, let dt=ME(ƒt,ƒt+1) and ƒt+1=MC(ƒt, dt) denote the motion estimation and motion compensation operators, respectively. Furthermore, let Ψ denote the orthonormal sparsity transform such that xt=Ψƒt and ƒt=Ψ−1xt. Note that the images f themselves need not be sparse. Thus, f represents the (possibly non-sparse) image frame and x=Ψf is the sparse representation of that image frame. For example, an MRI image (i.e. f) is not generally sparse. However, if one takes its wavelet transform, the resulting vector of wavelet coefficients (i.e. x) is (approximately) sparse.
An algorithm according to an embodiment of the invention for MC-MCS can then be stated as follows, with reference to the flow chart of
The minimization in EQ. (3) can be rewritten as such that y=Aβ, (5)
arg
βmin∥WT
where WT
In other words, the diagonal weights corresponding to the elements within the estimated support are set to zero and outside the estimated support (i.e. in TC are set to one.
A framework according to an embodiment of the invention can be further extended to not only include the sparse support information, but also information about the magnitudes of the sparse coefficients. Let xP denote a prediction of the sparse vector x. The information about both the support and the magnitudes of the sparse coefficients contained in this prediction can be incorporated into the above framework through the use of the weights. The weighted ll minimization can be written as
argβmin∥Wβ∥1, such that y=Aβ. (7)
The weights in the above equation can be selected in one of several ways.
(1) Set the diagonal weights as
(2) Estimate the support set of xP using V={iε[1,N]:(xP)i2>α}. Then the weights can be determined using
(3) Given V, the support set of xP, the weights can be determined using
(4) Given V, the support set of xP, the weights can be determined using
It should be noted that a framework according to an embodiment of the invention can incorporate other prediction models, depending on the application. For example, predictive models for multi-view images, predictive models based on image registration, or other predictive models (single/double exponential signal decay models for T2, diffusion models, etc.) can be incorporated into a framework according to an embodiment of the invention.
While the weight formulae written above assume a noise free measurement environment, they can be easily modified to incorporate noisy measurements using by changing the data consistency constraint to ∥y−Aβ∥2<ε, where the value of ε is adjusted depending on the amount of noise in the measurements.
Experiments were conducted to evaluate the performance of an MCMCS technique according to an embodiment of the invention. Three image sequences were used in the experiments. A phantom sequence was created by incorporating motion into a cropped Shepp-Logan phantom. The dimensions of each image frame were 128×128 pixels and there were 10 image frames the sequence. A heart sequence was a cardiac cine MR image sequence with 256×256 pixels and 20 frames, i.e., phases within a cardiac cycle. A mobile sequence was a video sequence. This sequence was cropped to 128×128 pixels and 10 frames were used in the experiments. Sample frames from each sequence are displayed in
A framework according to an embodiment of the invention is also general enough to accommodate different motion estimation techniques. Two different motion estimation techniques were used to estimate motion: a simple translational block-based motion estimation with block size of 8×8 and exhaustive search using the Mean Absolute Difference (MAD) metric was used for the heart and mobile sequences. For the phantom sequence, whose complex motion cannot be well approximated by translational block motion, a more general motion estimation model was used that allowed for translation and an affine coordinate transformation for the phantom sequence.
The first results are from a set of experiments designed to investigate the accuracy of support estimation with and without motion compensation. Let b %-energy support denote the set Ut={iε[1,N]:(xt)i2>ζt} where ζ is the largest real number for which at least b % of the signal energy at time t is contained in Ut. The value of b can be selected to contain most of the signal energy, e.g. b=99. The b %-energy support is useful in evaluating supports of signals that are not in a strict sense sparse, but are compressible. Since MCS relies on the assumption that the support of the signal changes slowly over time, the change in b % energy support between consecutive time frames can be considered as a measure of expected MCS performance. In other words, if |Ut\Ut−1|/|Ut| is small, MCS can be expected to perform well. Similarly, the mismatch in b % support between an image frame and its motion compensated estimate can be used as a measure of expected MC-MCS performance. This mismatch can be defined as |Ut\UtMC|/|Ut|, where UtMC={iε[1,N]:({circumflex over (x)}tMC)i2>ζt}. The quantities |Ut\Ut−1|\|/|Ut| vs. |Ut\UtMC|/|Ut| were evaluated for the three image sequences described earlier and the results are presented in
Experiments were performed to compare a MC-MCS method according to an embodiment of the invention to MCS and conventional CS algorithms. For the first set of experiments, random Gaussian measurements were used. At t=0, M=0.5N (i.e. 50% measurements) and the set T was empty. For the remaining frames (t>0), different undersampling ratios ranging from M=0.15N to M=0.25N were tested. The NESTA toolbox, available at http://acm.caltech.edu/nesta/, was used to solve the conventional CS minimization and an MCS and an MC-MCS algorithm according to an embodiment of the invention were implemented using the NESTA toolbox. Normalized Root Squared Error (NRSE), defined as et=√{square root over (∥ft−{circumflex over (f)}t∥22/∥ft∥22)} was used to compare the performances of different algorithms. For each undersampling ratio, experiments were repeated using 5 random realizations of the measurement matrix, and the NRSE values obtained from these realizations were averaged.
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.
The computer system 91 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.
This application claims priority from “Motion-Compensated Compressed Sensing for Dynamic Imaging”, U.S. Provisional Application No. 61/399,516 of Kim, et al., filed Jul. 12, 2010, the contents of which are herein incorporated by reference in their entirety.
| Number | Date | Country | |
|---|---|---|---|
| 61399516 | Jul 2010 | US |
