The field of the invention is imaging systems and methods, including medical imaging systems and methods. More particularly, the field of the invention is systems and methods for constrained image reconstruction that allow desirable levels of signal-to-noise ratio to be achieved without sacrificing spatial resolution.
In medical imaging, as well as other imaging technologies, signal-to-noise ratio (“SNR”) is utilized as a quantitative measure of image quality. Generally, SNR is defined as the ratio between the mean intensity value and the root-mean-square (“RMS”) noise, σ, in an image. The term “net signal” refers to the difference between an average signal value over the image, and background values, whereas the term RMS noise refers to the standard deviation of the noise value in the image. As SNR decreases in a medical image, it becomes increasingly more difficult to differentiate between anatomical features and other clinical findings of importance to the clinician. Thus, it is generally desirable to provide a relatively high SNR in medical imaging applications.
Spatial resolution is typically balanced with respect to the achievable SNR in an image. For example, the traditional relationship between spatial resolution, Δx, in a two-dimensional image and the achievable SNR in that image is given by:
Thus, as spatial resolution becomes finer, the noise variance, σ2, present in the image increases significantly. The SNR characteristics of an image are also governed by the manner in which image data is initially acquired. In many instances, it may be beneficial to alter data acquisition parameters, such as decreasing radiation dose or scan time, to the betterment of the subject of the examination, but to the detriment of SNR in the resultant images.
When parameters of an x-ray imaging study, such as tube current and tube current time product, “mAs,” are varied in order to decrease the radiation dose imparted to the subject, the SNR of the resultant images suffers. For example, decreasing tube current produces a related decrease in radiation dose; however, the noise present in the resultant images is increased, thereby affecting SNR in accordance with the following relationship:
where μ is the measured linear attenuation coefficient and σ is the RMS noise. Thus, if mAs is reduced by half, SNR will decrease by a factor of √{square root over (½)}, which corresponds to about a thirty percent decrease in SNR. Thus, while decreasing mAs during an x-ray imaging study provides a beneficial decrease in radiation dose imparted to the subject being imaged, the resultant images suffer from increased noise and, therefore, decreased SNR. Such images have limited clinical value.
Depending on the technique used, many magnetic resonance imaging (“MRI”) scans currently require many minutes to acquire the necessary data used to produce medical images. The reduction of this scan time is an important consideration, since reduced scan time increases patient throughout, improves patient comfort, and improves image quality by reducing motion artifacts. Many different strategies have been developed to shorten the scan time, including so-called parallel MRI (“pMRI”) techniques.
While the use of parallel MRI acts to decrease the amount of time required to image a subject without increasing gradient switching rates or RF power, parallel MRI methods are plagued with losses in signal-to-noise ratio (“SNR”). In general, the SNR of an image reconstructed using parallel MRI methods is decreased in accordance with the following relationship:
where g is the so-called geometry factor, or “g-factor,” and R is the acceleration factor, which describes the degree of undersampling employed and is related to, and generally limited by, the number of receiver coils in the array. Thus, parallel MRI methods suffer from a reduction in achievable SNR, offsetting the benefits provided by decreased scan time requirements.
It would therefore desirable to provide a method for image reconstruction or image processing in which an image having high SNR and high spatial resolution can be produced.
The present invention overcomes the aforementioned drawbacks by providing a method for producing an image having a high signal-to-noise ratio (SNR) by imparting the high SNR characteristics of a low spatial or temporal resolution prior image to a target image. An image to be enhanced is provided, the provided image including a previously reconstructed image or an image presently reconstructed from acquired image data. A prior image is produced from the provided image, for example, by filtering the provided image such that noise from the provided image is substantially suppressed in the prior image. Synthesized image data is produced by performing a forward projection of the provided image. A sparsified image is produced by subtracting the prior image and the provided image. A target image having a higher SNR than the provided image is reconstructed using the sparsified image, the provided image, and the synthesized image data. The provided image may be, for example, a medical image produced by at least one of an x-ray imaging system, an x-ray computed tomography system, a C-arm x-ray imaging system, an x-ray tomosynthesis imaging system, an x-ray projection imaging system, a magnetic resonance imaging system, a positron emission tomography system, a single photon emission computed tomography system, an ultrasound imaging system, and an optical imaging system.
The foregoing and other aspects and advantages of the invention will appear from the following description. In the description, reference is made to the accompanying drawings which form a part hereof, and in which there is shown by way of illustration a preferred embodiment of the invention. Such embodiment does not necessarily represent the full scope of the invention, however, and reference is made therefore to the claims and herein for interpreting the scope of the invention.
Generally speaking, the method of reconstructing an image from a set of data includes a series of numerical steps to estimate a target image, I, from the measured data samples, Y. More specifically, the image reconstruction should fulfill the following consistency condition:
AI=Y (4);
where A is a system matrix. In general, the system matrix, A, can be viewed as a forward projection operator that relates the target image, I, to the acquired data samples, Y. When dealing with computed tomography (“CT”) imaging, the system matrix can include a reprojection operation, while in magnetic resonance imaging (“MRI”), it can include a Fourier transform operation. The consistency condition of Eqn. (4), put in other words, states that when an image is faithfully reconstructed, the forward operation should substantially mimic the actual data acquisition procedure in order to generate a correct estimate of the measured projection data.
A method for reconstructing a quality desired, or target, image is provided herein. In general, a so-called “prior image” is employed to constrain an iterative image reconstruction method, in which the principles of compressed sensing (“CS”) are utilized. For example, in addition to the sparsifying transforms commonly used in CS, an image is further sparsified by subtracting the prior image from the target image. In this manner, the image reconstruction method is referred to as prior image constrained compressed sensing, or “PICCS.” Using PICCS, an image can be accurately reconstructed using a substantially fewer number of samples than required by CS methods. Additionally, it is a discovery that the noise characteristics from the prior image are imparted to the target image that is reconstructed.
Given a prior image, IP, and a target image to be reconstructed, I, the PICCS image reconstruction method may be implemented by minimizing the following objective function:
α∥Ψ1(I−IP)∥1+(1−α)∥Ψ2I∥1 (5);
where Ψ1 and Ψ2 are sparsifying transforms, ∥ . . . ∥1 is an L1-norm operation, and α is a weighting parameter that is utilized to control the relative weight of the two terms in the objective function of Eqn. (5). The data consistency condition presented in Eqn. (4) is imposed during the minimization of the objective function presented in Eqn. (5). It is noted that the L1-norm operation has the form:
for an N-dimensional vector, x. More generally, a deviation from the true L1-norm is possible while still maintaining adequate image quality in the target image. For example, the objective function of Eqn. (5) can be generalized as:
α∥Ψ1(I−IP)∥pp+(1−α)∥Ψ2I∥pp, (7);
where ∥ . . . ∥p is an Lp-norm operation having the form:
As noted above, preferably, p=1.0; however, in the alternative, different values of p are possible. It should be appreciated by those skilled in the art that the further the value of p deviates from p=1.0, generally, the more degradation will be evident in the reconstructed target image.
The sparsifying transforms, Ψ1 and Ψ2, in Eqn. (5) are, in general, different; however, in the alternative, Ψ1 and Ψ2 may be the same sparsifying transform. Exemplary sparsifying transforms include a wavelet transform, a first order finite difference, a second order finite difference, and a discrete gradient transform, such as, for example, a discrete gradient transform, ∇m,n, having the following form:
∇m,nI(m,n)=√{square root over ((I(m+1,n)−I(m,n))2+(I(m,n+1)−I(m,n))2)}{square root over ((I(m+1,n)−I(m,n))2+(I(m,n+1)−I(m,n))2)}{square root over ((I(m+1,n)−I(m,n))2+(I(m,n+1)−I(m,n))2)}{square root over ((I(m+1,n)−I(m,n))2+(I(m,n+1)−I(m,n))2)} (9);
where the indices m and n indicate the location of a pixel in an image, I. The image specified as ∇m,nI(m,n) is commonly referred to as the “gradient image.”
Both of the terms in the objective function of Eqn. (5) are important. As a result of their importance, the selection of the weighting parameter, α, is utilized to control the overall image reconstruction process. Therefore, the selection of the weighting parameter, α, will depend on the choice of the prior image, IP, and also the clinical application at hand. For example, the second term in the objective function of Eqn. (5), (1−α)∥Ψ2I∥1, mitigates streaking, or other, artifacts that are potentially inherited from the prior image, IP. By way of further example, selecting a weighting parameter of α≈0.3-0.7 is generally sufficient for cardiac imaging applications.
To better incorporate the consistency condition of Eqn. (4) into the overall image reconstruction, the method of Lagrange multipliers may be utilized. In such a manner, the consistency condition is employed to add a further constraint on the minimization of the objective function set forth in Eqn. (5). A new objective function is thus produced, which has the form:
α∥Ψ1(I−IP)∥1+(1−α)∥Ψ2I∥1+λ∥X∥22 (10);
where λ is the Lagrange multiplier, X is a difference matrix defined below in Eqn. (12), and ∥ . . . ∥22 is a squared L2-norm operation, which, for an N-dimensional vector, x, has the form:
The difference matrix in Eqn. (10) accounts for the consistency condition of Eqn. (4), and has the following form:
X=AI−Y (12).
It is noted, again, that the L1-norm operations in Eqn. (10) can similarly be replaced with the more general Lp-norm operations, as discussed above. The Lagrange multiplier, λ, is determined empirically for the particular imaging system employed when practicing the present invention. For example, the Lagrange multiplier, λ, is determined by a pre-determined tradeoff between the target data consistency requirement and the similarity to the prior image, Ip. When a large Lagrange multiplier, λ, is selected, the reconstructed image has lower noise variance; however, this may be achieved as a loss of the high spatial resolution characteristic of the prior image. Similarly, when a smaller Lagrange multiplier, λ, is used, the high spatial resolution characteristic of the prior image is well preserved, but the noise variance can be high in the target image. Such a situation affects the contrast-to-noise ratio achievable by the imaging system utilized. As will be described below in detail, however, a large Lagrange multiplier, λ, need not be employed to provide an increase, or maintenance of an existing desirable level of, signal-to-noise ratio (“SNR”). Instead, a high SNR prior image, Ip, may be produced and the SNR of this prior image is imparted to the target image, I.
The objective function presented in Eqn. (10) can further be altered in order to account for noise characteristics of the imaging system. In such a manner, the following objective function is minimized:
α∥Ψ1(I−IP)∥1+(1−α)∥Ψ2I∥1+λ(XTDX) (13);
where XT is the transpose of the difference matrix, X, and D is a system noise matrix, which is a diagonal matrix having the following matrix elements:
where σn2 an is the noise variance, and is a parameter indicative of noise in the imaging system employed when practicing the present invention. For example, in an x-ray imaging system, the noise parameter, σn2, is the noise variance associated with the nth x-ray detector. Alternatively, in an MR imaging system, the noise parameter, σn2, is estimated noise variance in the nth receiver coil. Again, the L1-norm operations in Eqn. (13) may similarly be replaced with the more general Lp-norm operations.
In PICCS image reconstruction methods, the prior image, IP, plays several roles. First, it serves as a seed image in the iterative reconstruction, which accelerates the overall image reconstruction method. Second, the prior image, IP, is employed to further sparsify the target image, I, and, thus, serves as yet another sparsifying transform. In addition, as will be described below in detail, the SNR of the prior image, Ip, is imparted to the target image, I, thereby allowing the reconstruction of a target image, I, having both a high SNR and high spatial resolution.
A brief discussion of possible prior images, IP, is provided below with respect to different imaging modalities; however, it should be appreciated by those skilled in the art that prior images, IP, other than those expressly described herein can be employed depending on the clinical application. As referred to herein, a prior image, IP, is an image of the subject that includes a priori information indicative of the target image to be reconstructed. Generally, the prior image, IP, is formed from a previously performed imaging study. For example, and as will be described in detail below, the prior image, IP, may be formed by filtering a preexisting image of the subject under examination. Exemplary filtering may include low pass filtering, resulting in a prior image, IP, having low noise and low spatial resolution characteristics. As will be described below, the low noise of the prior image, IP, is imparted to the target image of the subject, while the spatial resolution of the target image is determined not by the prior image, IP, but by characteristics of the image data utilized to reconstruct the target image.
With reference now to
The initialization of the second term in the objective function of Eqn. (5), (1−α)∥Ψ2I∥1, begins at step 114 where the estimate of the target image, I, is sparsified through the application of the sparsifying transform, Ψ2. Subsequently, the L1-norm of this sparsified target image estimate is calculated at step 116. When the discrete gradient transform, ∇m,n, is selected as the sparsifying transform, Ψ2, steps 114 and 116 can be viewed as calculating the total variation, TV, of the target image estimate, which has the form:
TV(I)=∥∇I∥1=Σ|∇I∇ (15).
After the L1-norm of the sparsified target image estimate is calculated, the result is weighted by (1−α), as indicated at step 118. The objective function of Eqn. (5) is subsequently produced in step 120 by adding the first and second terms together. This objective function is then minimized, as indicated at step 122, using, for example, a nonlinear conjugate gradient method. The minimization process proceeds until a stopping criterion is satisfied. The stopping criterion includes, for example, comparing the current estimate of the target image with the estimate of the target image from the previous iteration. Such a stopping criterion has the following form:
where, Iij(k+1) is the value of the (k+1)th estimate of the target image at the pixel location (i, j), and Iij(k) is the value of the kth estimate of the target image at the pixel location (i, j).
With reference now to
X=AI−Y (17).
Thus, the difference matrix is determined by applying the system matrix, A, to the estimate of the target image, I, and subsequently subtracting the acquired image data, Y, that corresponds to the target image. The square of the L2-norm of the difference matrix, X, is calculated next at step 210. After the square of the L2-norm of the difference matrix, X, has been produced, the Lagrange multiplier, λ, is determined and employed to weight the difference matrix, X, as indicated at step 212. As described above, the Lagrange multiplier is empirically determined by a value selected by the user based on the clinical application at hand. The objective function of Eqn. (10) is subsequently produced in step 220 by adding the first, second, and third terms together. This objective function is then minimized, as indicated at step 222, using, for example, a nonlinear conjugate gradient method. The minimization process proceeds until a stopping criterion is satisfied, as described above.
With reference now to
As described above, σn2 is the noise variance, which is a parameter indicative of noise in the imaging system employed to acquire the relevant image data. For example, in an x-ray imaging system, the noise parameter, σn2, may be the noise variance associated with the nth x-ray detector. Alternatively, in an MRI system, the noise parameter, σn2, may be the estimated noise variance in the nth receiver coil. After the system noise matrix, D, has been produced, the following matrix multiplication is performed:
X
T
DX (19);
as indicated at step 312. The result of this operation is subsequently scaled by the Lagrange multiplier, as indicated at step 314. The objective function of Eqn. (13) is subsequently produced in step 320 by adding the first, second, and third terms together. This objective function is then minimized, as indicated at step 322, using, for example, a nonlinear conjugate gradient method. The minimization process proceeds until a stopping criterion is satisfied, as described above.
With the appropriate choice of a prior image, Ip, a target image, I, can be reconstructed with both high SNR and spatial resolution. In general, the SNR of the prior image, Ip, is imparted to the target image, I, by way of the aforementioned PICCS image reconstruction methods. Therefore, consideration of how the prior image, Ip, is produced can provide benefits to the overall image quality of the target image, I. For example, the quality of preexisting medical images can be improved by forming a low noise, low spatial resolution prior image from the existing images and appropriately resampling the existing images to produce image data from which higher quality target images may be reconstructed.
It is noted that the present invention provides an image reconstruction method applicable to a number of different imaging modalities including x-ray computed tomography (“CT”), x-ray C-arm imaging, x-ray tomosynthesis imaging, magnetic resonance imaging (“MRI”), positron emission tomography (“PET”), single photon emission computed tomography (“SPECT”), optical imaging, and ultrasound imaging. More specifically, the present invention provides an image reconstruction method that provides an increase in achievable signal-to-noise ratio (“SNR”) in reconstructed images, without a significant decrease in spatial resolution. This method is also applicable for post-processing applications, in which the SNR of an existing image is enhanced. Additionally, the method is applicable to other imaging modalities, including optical imaging using light in both the visible and other spectral ranges.
As described above, the presently described constrained image reconstruction method is applicable to many different medical imaging modalities and may be utilized in many different clinical applications. A number of such exemplary clinical applications are described below to illustrate the broad scope of the present invention. Such embodiments do not necessarily represent the full scope of the invention, however, and reference is made therefore to the claims and herein for interpreting the scope of the invention.
Referring now to
When the acquired image data from which the provided image was originally reconstructed is not available, the provided image is resampled, for example by forward projection, to produce synthesized image data, as indicated at step 406. The choice of resampling kernel may vary depending on the desired clinical application and on considerations of transform domain characteristics. Exemplary resampling kernels include those that transform the provided image into x-ray projection space, for example by Radon transform; k-space, for example by Fourier transform; and a wavelet transform domain, for example by wavelet transform. In some instances, the provided image may also be resampled in the image domain prior to performing a forward projection. For example, the provided image may be oversampled to increase its spatial resolution. This oversampling provides an additional guarantee that the spatial resolution is properly encoded in the synthesized image data. Because the provided image is resampled in this manner, it does not matter how the provided image was originally acquired. Therefore, the image reconstruction method is “universal,” in as much as it is applicable to any such previously obtained image.
Using the produced prior image and the synthesized image data, a high spatial resolution, high SNR target image is reconstructed using a prior image constrained compressed sensing (“PICCS”) reconstruction technique, as indicated at step 408. As noted, the reconstructed target image depicts the same subject as the provided image, but contains a higher SNR by virtue of the high SNR prior image, and similar or increased spatial resolution by virtue of the synthesized image data. Exemplary image reconstruction includes the minimization of an objective function that includes terms related to the prior image and synthesized image data. For example, the objective function in Eqn. (7) may be modified as follows:
α∥Ψ1(I−IP)∥pp+(1−α)∥Ψ2I∥pp+λ∥RI−RI0∥22 (20);
where R is a resampling operator, such as those discussed above, and may in some instances be selected as an identity matrix; I0 is the provided image; and RI0 is the synthesized data. When R is selected as the identity matrix, the implication is to allow the forward projection of the target image, I, and of the provided image, I0, to equal their respective selves. In such an instance, the objective function provided in Eqn. (20) may be rewritten as:
α∥Ψ1(I−IP)∥pp+(1−α)∥Ψ2I∥pp+λ∥I−I0∥22 (21).
Thus, the image reconstruction method can be used to the advantage of enhancing preexisting images as well as guiding future data acquisitions. For example, an x-ray computed tomography data acquisition can be altered such that a lower x-ray radiation dose is imparted to the subject. In this instance the SNR of an image reconstructed from such lower dose image data would
With initial reference to
The rotation of the gantry and the operation of the x-ray source 513 are governed by a control mechanism 520 of the CT system. The control mechanism 520 includes an x-ray controller 522 that provides power and timing signals to the x-ray source 513 and a gantry motor controller 523 that controls the rotational speed and position of the gantry 512. A data acquisition system (“DAS”) 524 in the control mechanism 520 samples analog data from detector elements 518 and converts the data to digital signals for subsequent processing. An image reconstructor 525, receives sampled and digitized x-ray data from the DAS 524 and performs high speed image reconstruction. The reconstructed image is applied as an input to a computer 526 which stores the image in a mass storage device 528.
The computer 526 also receives commands and scanning parameters from an operator via console 530 that has a keyboard. An associated display 532 allows the operator to observe the reconstructed image and other data from the computer 526. The operator supplied commands and parameters are used by the computer 526 to provide control signals and information to the DAS 524, the x-ray controller 522 and the gantry motor controller 523. In addition, computer 526 operates a table motor controller 534 which controls a motorized table 536 to position the patient 515 in the gantry 512.
Referring particularly now to
The pulse sequence server 610 functions in response to instructions downloaded from the workstation 602 to operate a gradient system 618 and a radiofrequency (“RF”) system 620. Gradient waveforms necessary to perform the prescribed scan are produced and applied to the gradient system 618, which excites gradient coils in an assembly 622 to produce the magnetic field gradients Gx, Gy, and Gz used for position encoding MR signals. The gradient coil assembly 622 forms part of a magnet assembly 624 that includes a polarizing magnet 626 and a whole-body RF coil 628.
RF excitation waveforms are applied to the RF coil 628, or a separate local coil (not shown in
The RF system 620 also includes one or more RF receiver channels. Each RF receiver channel includes an RF amplifier that amplifies the MR signal received by the coil 628 to which it is connected, and a detector that detects and digitizes the I and Q quadrature components of the received MR signal. The magnitude of the received MR signal may thus be determined at any sampled point by the square root of the sum of the squares of the I and Q components:
M=√{square root over (I2+Q2)} (22);
and the phase of the received MR signal may also be determined:
The pulse sequence server 610 also optionally receives patient data from a physiological acquisition controller 630. The controller 630 receives signals from a number of different sensors connected to the patient, such as electrocardiograph (“ECG”) signals from electrodes, or respiratory signals from a bellows or other respiratory monitoring device. Such signals are typically used by the pulse sequence server 610 to synchronize, or “gate,” the performance of the scan with the subject's heart beat or respiration.
The pulse sequence server 610 also connects to a scan room interface circuit 632 that receives signals from various sensors associated with the condition of the patient and the magnet system. It is also through the scan room interface circuit 632 that a patient positioning system 634 receives commands to move the patient to desired positions during the scan.
The digitized MR signal samples produced by the RF system 620 are received by the data acquisition server 612. The data acquisition server 612 operates in response to instructions downloaded from the workstation 602 to receive the real-time MR data and provide buffer storage, such that no data is lost by data overrun. In some scans, the data acquisition server 612 does little more than pass the acquired MR data to the data processor server 614. However, in scans that require information derived from acquired MR data to control the further performance of the scan, the data acquisition server 612 is programmed to produce such information and convey it to the pulse sequence server 610. For example, during prescans, MR data is acquired and used to calibrate the pulse sequence performed by the pulse sequence server 610. Also, navigator signals may be acquired during a scan and used to adjust the operating parameters of the RF system 620 or the gradient system 618, or to control the view order in which k-space is sampled. The data acquisition server 612 may also be employed to process MR signals used to detect the arrival of contrast agent in a magnetic resonance angiography (“MRA”) scan. In all these examples, the data acquisition server 612 acquires MR data and processes it in real-time to produce information that is used to control the scan.
The data processing server 614 receives MR data from the data acquisition server 612 and processes it in accordance with instructions downloaded from the workstation 602. Such processing may include, for example: Fourier transformation of raw k-space MR data to produce two or three-dimensional images; the application of filters to a reconstructed image; the performance of a backprojection image reconstruction of acquired MR data; the generation of functional MR images; and the calculation of motion or flow images.
Images reconstructed by the data processing server 614 are conveyed back to the workstation 602 where they are stored. Real-time images are stored in a data base memory cache (not shown in
The present invention has been described in terms of one or more preferred embodiments, and it should be appreciated that many equivalents, alternatives, variations, and modifications, aside from those expressly stated, are possible and within the scope of the invention.
This invention was made with United States government support under EB005712 awarded by the National Institutes of Health. The United States government has certain rights in this invention.