HK1179397B - Method for radiation dose reduction using prior image constrained image reconstruction - Google Patents

Description

Method of radiometric reduction for image reconstruction using prior image constraints

Cross Reference to Related Applications

This application claims the benefit of U.S. patent application serial No. 12/783,058, filed on 9/5/2010 and entitled "method for Radiation Dose Reduction Using the former Image structured Image recovery".

Statement regarding federally sponsored research in the United states

This application was made with U.S. government support by the following agencies: national institute of health, NIH EB 005721. The united states government has certain rights in the invention.

Background

The field of the invention is medical imaging systems and methods. More particularly, the present invention relates to a method of image reconstruction that allows maintaining a desired level of signal-to-noise ratio while balancing other imaging considerations, such as reduction in radiation dose and scan time.

In medical imaging, and other imaging techniques, the signal-to-noise ratio ("SNR") is used as a quantitative measure of image quality. In general, SNR is defined as the ratio of the average intensity value in an image to the root mean square ("RMS") noise σ. The term "net signal" refers to the difference of the average signal value over the image from the background value, while the term RMS noise refers to the standard deviation of the noise values in the image. As the SNR in medical images decreases, it becomes increasingly difficult to distinguish anatomical features from other clinical findings that are important to the clinician. Therefore, it is generally desirable to maintain a relatively high SNR in medical imaging applications.

In computed tomography systems, an x-ray source emits a fan-shaped beam that is collimated to lie within an x-y plane of a Cartesian coordinate system, termed an "image plane". The x-ray beam passes through an object being imaged, such as a medical patient, and impinges upon an array of radiation detectors. The intensity of the emitted radiation depends on the attenuation of the x-ray beam by the object, and each detector produces a separate electrical signal, which is a measure of the beam attenuation. Attenuation measurements from all detectors are obtained separately to produce a so-called "transmission curve" (profile), "attenuation curve", or "projection".

The source and detector arrays in conventional CT systems are rotated on a gantry in the imaging plane and around the object so that the angle at which the x-ray beam intersects the object continuously changes. The transmission curves from the detector array at a given angle are referred to as "views" and "scans" of the object, including a set of views made at different angular orientations during one selection of the x-ray source and detector. In a 2D scan, data is processed to construct a two-dimensional image corresponding to a two-dimensional slice acquired through the object. The primary method for reconstructing an image from 2D data is known in the art as the filtered back projection technique. This image reconstruction process converts the attenuation measurements acquired during the scan into integers called "CT numbers" or "Hounsfield units", which are used to control the brightness of the corresponding pixel on the display.

Magnetic resonance imaging ("MRI") utilizes the phenomenon of nuclear magnetic resonance ("NMR") to produce images. When a substance such as human tissue is subjected to a uniform magnetic field (polarizing field B)₀) The magnetic moments of the nuclei within the tissue attempt to align with the polarizing field, but precess around the polarizing field in random order at its characteristic larmor frequency. If the substance or tissue is subjected to a magnetic field in the x-y plane and close to the larmor frequency (excitation field B)₁) Then net aligning the magnetic moment M_zCan be rotated or "deflected" into the x-y plane to produce a net transverse magnetic moment M_xy. In the excitation signal B₁After being terminated, signals are transmitted through the excited nuclei or "spins" and may be received and processed to form images.

When using these "MR" signals to generate images, gradient magnetic fields (G) are used_x、G_yAnd G_z). Typically, the region to be imaged is scanned by a sequence of measurement cycles in which these gradients vary depending on the particular localization method used. The resulting set of received MR signals is digitized and processed to reconstruct an image using one of a number of well-known reconstruction techniques.

The measurement cycle for acquiring each MR signal is performed under the direction of the pulse sequence generated by the pulse sequencer. Clinically available MRI systems store a library of pulse sequences that can be specified to meet the needs of many different clinical applications. Research-type MRI systems include a library of clinically proven pulse sequences, and the system may also develop new pulse sequences.

The MR signals acquired by an MRI system are signal samples of an object under examination in fourier space (or often referred to in the art as "k-space"). Each MR measurement cycle, or pulse sequence, generally samples a portion of k-space along the sampling trajectory characteristic of this pulse sequence. Most pulse sequences are sampled in k-space in a raster scan-like pattern sometimes referred to as a "spin-twist," "fourier," "line," or "cartesian" scan. Spin-twist scanning techniques employ variable amplitude phase encoding magnetic field gradient pulses to phase encode spatial information in the direction of the gradient prior to acquisition of MR spin echo signals. For example, in a two-dimensional implementation ("2 DFT"), a phase encoding gradient (G) is applied in one direction by applying a phase encoding gradient along that direction_y) To encode spatial information and subsequently in the presence of a readout magnetic field gradient (G) in a direction orthogonal to the phase encoding direction_x) The spin echo signal is acquired. The readout gradient present during spin echo acquisition encodes spatial information in the orthogonal direction. In a typical 2DFT pulse sequence, a phase encoding gradient pulse G_yIs incremented in the measurement cycle and "view" sequence acquired during the scan to produce a set of k-space MR data from which the entire image can be reconstructed.

There are also many other k-space sampling modes used by MRI systems. These include "radial", or "projection reconstruction" scans, in which k-space is sampled as a set of radial sampling trajectories extending from the middle of k-space. The pulse sequences for radial scanning are characterized by the absence of phase encoding gradients and the presence of readout gradients that change the direction from one pulse sequence view to the next. There are also many k-space sampling methods that are closely related to radial scanning and sample along a curved k-space sampling trajectory rather than along a straight radial trajectory.

An image is reconstructed from the acquired k-space data by converting the k-space data set to an image space data set. There are many different methods for performing this task, and this is often determined by the technique used to acquire k-space data. Using cartesian grids of k-space data derived from 2D or 3D spin-twist acquisition, for example, the most common reconstruction method is the inverse fourier transform ("2 DFT" or "3 DFT") along each of the 2 or 3 axes of the data set. Using radial k-space datasets and their variables, the most common reconstruction methods include "re-meshing" the k-space samples to create a cartesian grid of k-space samples and then performing a 2DFT or 3DFT on the re-meshed k-space datasets. In an alternative approach, the radial k-space dataset may also be converted to Radon space by performing a 1DFT per radial projection view, and then the Radon space dataset is converted to image space by performing a filtered back-projection.

According to standard image reconstruction theory, in order to reconstruct an image without aliasing artifacts, the sampling rate for acquiring image data must satisfy the so-called Nyquist (Nyquist) criterion, which is set on Nyquist-Shannon sampling theory. In addition, in standard image reconstruction theory, no specific prior information about the image is required. On the other hand, when some prior information about the desired image is available and suitably incorporated into the image reconstruction step, the image can be accurately reconstructed even if the Nyquist criterion is violated. For example, if the desired image is known to be circularly symmetric and spatially uniform, only one view of the parallel beam projections (i.e., one projection view) is required to accurately reconstruct the linear attenuation coefficient of the object. As another example, if it is known that the desired image includes only a single point, two orthogonal projections that intersect at that point are required to accurately reconstruct that image point. Thus, if prior information about the desired image is known, such as if the desired image is a set of sparsely distributed points, the image may be reconstructed from a set of data acquired in a manner that does not satisfy the Nyquist criterion. More generally, knowledge about the sparsity of the desired image may be employed to relax the Nyquist criterion; however, generalizing these issues to build a strict image reconstruction theory is an unusual task.

The nyquist criterion is used as one of the most important bases in the field of information science. However, it may also play a paramount role in modern medical imaging modalities, such as magnetic resonance imaging ("MRI") and x-ray computed tomography ("CT"). When the number of data samples acquired by the imaging system is less than the requirements specified by the nyquist criterion, artifacts appear in the reconstructed image. In general, such image artifacts include aliasing and streaking artifacts. In practice, the nyquist criterion is often violated, either intentionally or by unavoidable circumstances. For example, in order to shorten the data acquisition time in time-resolved MR angiography studies, undersampled projection reconstructions, or radial, acquisition methods are often deliberately introduced.

The risks associated with exposure to ionizing radiation for medical imaging, including x-ray computed tomography ("CT") and core muscle perfusion imaging ("MPI"), have become increasingly important considerations in recent years as the number of CT and nuclear MPI studies has increased dramatically. The reported effective radiation dose from cardiac CT angiography procedures is about 5-20 millisievert ("mSv") for male patients and even higher for female patients. In addition to this dose, there is a smaller radiation dose for the calcium-integrated CT scan that is conventionally performed prior to intravenous contrast agent injection. To perform CT-MPI as part of a full cardiac CT study, it would be required to acquire images about 20-30 times on the same region of the heart, resulting in an increase in radiation dose of about twenty to thirty times, which is an unacceptable level of radiation exposure.

When parameters of an x-ray imaging study, such as tube current and the product of tube current and time, "mAs" are altered to increase the radiation dose applied to the subject, the signal-to-noise ratio ("SNR") of the resulting image suffers. For example, reducing the tube current produces a relative reduction in radiation dose; however, the signal present in the resulting image is increased, thereby affecting the SNR according to the following relationship:

where μ is the measured linear attenuation coefficient and σ is the RMS noise. Thus, if mAs is halved, the SNR will decreaseThis corresponds to a reduction in SNR of about 30%. Thus, while reducing mAs during an x-ray imaging study provides a beneficial reduction in the radiation dose applied to the subject being imaged, the resulting image suffers from increased noise, and therefore, reduced SNR. Such images have limited clinical value.

Depending on the technique used, many MR scans currently require many minutes to acquire the necessary data used to generate the medical images. Reducing this scan time is an important consideration because reduced scan time increases patient flow, improves patient comfort, and improves image quality by reducing motion artifacts. Many different strategies have been developed to shorten this scan time.

One such strategy is commonly referred to as "parallel MRI" ("pMRI"). Parallel MRI techniques use spatial information from a radio frequency ("RF") receiver coil array to replace spatial encoding that would otherwise have to be obtained in a sequential manner using RF pulses and magnetic field gradients, such as phase and frequency encoding gradients. Each of the spatially independent receiver coils of the array carries specific spatial information and has a different spatial sensitivity profile. This information is used to obtain a complete spatial encoding of the received MR signals, e.g. by combining the simultaneously acquired data received from each individual coil. Parallel MRI techniques allow undersampling of k-space by reducing the number of phase encoded k-space sampling lines acquired while maintaining the maximum range covered in fixed k-space. The combination of the independent MR signals generated by the independent receiver coils enables a reduction of the acquisition time required for the image, the reduced factor being related to the number of receiver coils, compared to conventional k-space data acquisition.

Although using parallel MRI reduces the amount of time required to image a subject without increasing gradient switching rates or RF power, parallel MRI methods suffer from signal-to-noise ratio ("SNR") loss. In general, the SNR of an image reconstructed using a parallel MRI method is reduced according to the following relationship:

where g is the so-called geometry factor, or "g-factor," and R is the acceleration factor, describing the degree of undersampling employed and associated with and limited by the number of receiver coils in the array. Thus, the parallel MRI approach is limited by the reduction in achievable SNR, offsetting the benefits provided by the reduced scan time requirements.

It is therefore desirable to provide a method of reconstructing an image of an object from medical image data such that a higher signal-to-noise ratio ("SNR") can be obtained than with currently available methods. It is further desirable to provide a method of reconstructing an image of an object in the above-described manner, whereby the trade-off between SNR and other considerations, such as radiation dose in x-ray imaging and scan time in magnetic resonance imaging, can be balanced without a significant loss in SNR.

Disclosure of Invention

The present invention overcomes the above-described disadvantages and provides a method of reconstructing an image of an object with a higher signal-to-noise ratio ("SNR") inherited from a previous image formed by averaging image data acquired in a slice direction, but without the lower spatial resolution often present in such a previous image. A plurality of images is reconstructed from the acquired image data, wherein the plurality of images correspond to a respective plurality of image slices. For each image slice position, a weighting value is calculated using the reconstructed images. For example, a comparison between adjacent image slices is used to calculate the desired weighting value for a given slice position. The weighted image data is formed by applying the calculated plurality of weights to the acquired image data, thereby filtering the acquired image data in a direction orthogonal to the orientation of the plurality of image slices (the so-called "slice direction"). A previous image is reconstructed from this weighted image data and then used in the PICCS reconstruction algorithm to reconstruct a target image of the subject that inherits the higher SNR of the previous image, but retains the higher spatial resolution of the acquired image data.

The present invention provides image reconstruction methods applicable to a plurality of different imaging modalities, including x-ray computed tomography ("CT"), x-ray c-arm imaging, magnetic resonance imaging ("MRI"), positron emission tomography ("PET"), and single photon emission computed tomography ("SPECT"). More particularly, the present invention provides an image reconstruction method that provides an increase in the achievable signal-to-noise ratio ("SNR") in the reconstructed image without a significant reduction in spatial resolution. In this way, tradeoffs can be made with respect to other imaging considerations, such as the radiation dose applied to the subject and the overall scan time, without a significant reduction in SNR compared to currently available image reconstruction methods.

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 preferred embodiments of the invention. These embodiments do not necessarily represent the full scope of the invention, however, and reference is therefore made to the claims and herein for interpreting the scope of the invention.

Drawings

FIG. 1 is a flow chart illustrating the steps of an exemplary image reconstruction method employed in practicing the present invention;

FIG. 2 is a flow chart illustrating the steps of another exemplary image reconstruction method employed in practicing the present invention;

FIG. 3 is a flow chart illustrating the steps of yet another exemplary image reconstruction method employed in practicing the present invention;

FIG. 4A is a pictorial diagram of an exemplary x-ray computed tomography ("CT") imaging system;

FIG. 4B is a block diagram of the CT imaging system of FIG. 4A;

FIG. 5A is a pictorial view of an exemplary C-arm x-ray imaging system;

FIG. 5B is a block diagram of the C-arm x-ray imaging system of FIG. 5A;

FIG. 6A is a pictorial view of an x-ray source and detector in the C-arm x-ray imaging system of FIG. 5A;

FIG. 6B is a pictorial view of a C-arm scan path employed by the C-arm x-ray imaging system of FIG. 5A;

FIG. 7 is a flow chart illustrating steps of an exemplary image reconstruction method according to the present invention employed when using the x-ray CT system of FIGS. 4A and 4B, or the x-ray C-arm imaging system of FIGS. 5A and 5B;

FIG. 8 is a graphical illustration of the weighted averaging and image reconstruction process of FIG. 7;

FIG. 9 is a block diagram of an exemplary magnetic resonance imaging ("MRI") system in which the present invention may be employed;

FIG. 10 is a diagrammatic view of an exemplary gradient echo pulse sequence for guiding the MRI system of FIG. 9 to acquire image data by sampling along a series of radial k-space trajectories in accordance with the present invention;

FIG. 11 is a graphical representation of a k-space sampling pattern generated using the pulse sequence of FIG. 10; and

FIG. 12 is a flowchart illustrating steps of an exemplary method for generating prior images and reconstructing an image in accordance with the present invention when an MRI system such as that shown in FIG. 9 is employed.

Detailed Description

Generally, a method of reconstructing an image from a set of data includes a series of steps to estimate a desired image I from measured data samples Y. More specifically, the image reconstruction should satisfy the following compatibility condition:

AI=Y （3）;

where A is the system matrix. In general, the system matrix, a, can be viewed as a forward projection operator (forward projection operator) that relates the desired image I to the acquired data samples Y. When dealing with computed tomography ("CT") imaging, the system matrix may include reprojection operations, whereas in magnetic resonance imaging ("MRI"), the system matrix may include fourier transform operations. In other words, the compatible condition of equation (3) indicates that the forward operation should substantially emulate the actual data acquisition steps when the image is faithfully reconstructed, thereby yielding a correct estimate of the measured projection data.

Turning now to the method of the present invention, a method for reconstructing an image of desired quality is provided. In general, the iterative image reconstruction method is constrained with "prior images", where the theory of compressed sensing ("CS") is used. For example, in addition to the sparse transform commonly used in CS, the image is further thinned out by subtracting the previous image from the desired image. In this way, the image reconstruction method is referred to as compressive sensing of prior image constraints, or "PICCS". Using PICCS, images can be accurately reconstructed using a much smaller number of samples than required by the CS method.

More particularly, given a preceding image I_PAnd a desired image I to be reconstructed, the method for image reconstruction of the invention is realized by minimizing the following objective function:

α||Ψ₁(I-I_P)||₁+(1-α)||Ψ₂I||₁（4）;

therein Ψ₁And Ψ₂Is sparse transformation, | |. | calculation₁Is L₁Modulo arithmetic, and α is a regularization parameter for the relative weight of two terms in the objective function of the control formula (4). Note that:

l representing an N-dimensional vector x₁-a mould. More generally, from true L₁Deviations of the modes are possible while still leaving sufficient image quality in the desired image. For example, the objective function of equation (4) may be generalized as:

α||Ψ₁(I-I_P)||_p+(1-α)||Ψ₂I||_p,（6）;

wherein | |_pIs a modulo operation, having the form:

as described above, preferably, p ═ 1.0; however, in alternative embodiments, different values of p are possible. It will be appreciated by those skilled in the art that the more the value of p deviates from p by 1.0, in general, more degradation will become apparent in the reconstructed desired image.

Sparse transform in equation (4), Ψ₁And Ψ₂Generally speaking, different; however, in alternative embodiments, Ψ₁And Ψ₂May be the same sparse transform. Exemplary sparse transforms include wavelet transforms, first order finite differences, second order finite differences, and discrete gradient transforms such as, for example, discrete gradient transforms having the form

Where the indices m and n indicate the location of the pixels in the image I. Is defined asIs often referred to as a "gradient image".

Both terms in the objective function of equation (4) are important. As a result of their importance, the choice of the regularization parameter α is used to control the overall image reconstruction process. The choice of the regularization parameter α will therefore depend on the preceding image I_PThere are also upcoming clinical applications of (a) to (b). For example, the second term (1- α) | | Ψ in the objective function of equation (4)₂I||₁Reduced streaking, or potentially from a preceding image I_POther artifacts inherited in (1). By way of further example, selecting a regularization parameter α ≈ 0.3-0.7 is generally sufficient for cardiac imaging applications.

To better incorporate the compatible condition of equation (3) into the overall image reconstruction, a method of lagrange multipliers is used. In this way, a compatibility condition is employed to add a further constraint on the minimization of the objective function set in equation (4). A new objective function is thus generated, having the form:

where λ is the Lagrangian multiplier, X is the difference matrix, andis L₂The square of the modulo operation, which for an N-dimensional vector, x, has the form:

the difference matrix in formula (9) takes into account the compatibility condition of formula (3) and has the following form:

X=AI-Y （11）。

the lagrange multiplier, λ, is empirically determined for the particular imaging system used in practicing the present invention. For example, by desired data compatibility requirements and with the preceding image I_pIs determined to determine the lagrange multiplier lambda. When a larger Lagrangian multiplier λ is selected, the reconstructed mapLike having low noise variation; however, this can be realized as a loss of the higher spatial resolution characteristic of the previous image. Similarly, when a smaller lagrange multiplier λ is used, the higher spatial resolution characteristics of the prior image are better maintained, although the noise variation may be higher in the desired image. Such a situation affects the contrast-to-noise ratio achievable by the imaging system used. However, as will be described below, a larger lagrangian multiplier λ is not required to provide an increase, or maintenance, of the existing desired level of signal-to-noise ratio ("SNR"). Instead, the preceding image I of higher SNR_pIs generated and the SNR is applied to the desired image I.

The objective function presented in equation (9) may be further varied to account for the noise of the imaging system. In this way, the following objective function is minimized:

α||Ψ₁(I-I_P)||₁+(1-α)||Ψ₂I||₁+λ(X^TDX) (12);

wherein X^TIs the transpose of the difference matrix X and D is the system noise matrix, which is a diagonal matrix with the following matrix elements:

whereinIs the noise variation and is a parametric indication of noise in the imaging system used in practicing the invention. For example, in an x-ray imaging system, noise parametersIs the noise variation associated with the nth x-ray detector. Optionally, in the MR imaging system, the noise parameterIs the estimated noise variation in the nth receiver coil.

In the method of the invention, a preceding image I_PPlaying several roles. First, it serves as a seed image in the iterative reconstruction, speeding up the overall image reconstruction method. Second, the previous image I_PIs used to further sparsify the desired image I and is thus used as yet another sparse transform. Further, as will be described in detail below, the previous image I_pSNR of (d) is applied to the desired image I, thereby allowing reconstruction of the desired image I with higher SNR and higher spatial resolution.

The following provides possible prior images I with reference to different imaging modalities_PA brief discussion of (1); however, it will be appreciated by those skilled in the art that the previous image I_POther than those explicitly described herein, may be used depending on the clinical application. As described herein, the previous image I_PIs an image of the object comprising a priori information representing a desired image to be reconstructed. Preceding image I_PMay be derived from previously performed imaging studies, or may be derived from the number of images acquired for a desired imageFrom which image data acquired in the same process is reconstructed. In general, the prior image I is acquired using the same imaging modality as the desired image_P(ii) a However, there are applications where the prior image I may be obtained from a different imaging modality than the desired image_PSuch as, for example, when a combined PET-CT system is employed.

Referring now to FIG. 1, one implementation of the method of the present invention employs the objective function of equation (4) and begins by initializing the regularization parameter α, as shown at step 100. The choice of the regularization parameter a determines a trade-off between the sparsity of the desired image and the effect of the previous image on the desired image. Thus, the value of the regularization parameter, α, will vary depending on the clinical setting to be performed. For example, values of α ≈ 0.3-0.7 are generally sufficient for cardiac imaging applications. Subsequently, the first and second terms in the objective function of equation (4) are initialized, as shown in steps 102 and 104, respectively. First term α | | Ψ₁(I-I_P)||₁Is initiated in step 106, wherein a previous picture I_PIs subtracted from the estimate of the desired image I to produce a "difference image". Preceding image I_PWill depend on the imaging modality and the particular clinical application. Accordingly, various alternative embodiments of these options will be discussed in detail below. Subsequently, as shown in step 108, by applying the sparse transform Ψ₁To sparsify the difference image. As described above, the sparse transform Ψ₁Any number of mathematical operations may be used including wavelet transforms, first order finite differences, second order finite differences, and discrete gradient transforms. Then, in step 110, L of this thinned difference image is calculated₁-a mould. The results of this process are then weighted by the regularization parameter α, as shown in step 112.

The second term (1- α) | | Ψ in the objective function of equation (4)₂I||₁Begins in step 114 by applying the sparse transformation Ψ₂To sparsify the estimate of the desired image I. Then, in step 116, this desired graph that is being thinned out is computedLike estimated L₁-a mould. When discrete gradient changesIs selected as the sparse transform Ψ₂Then, steps 114 and 116 may be considered as computing a total variation (totalvanion) TV of the desired image estimate, having the form:

l after calculating the sparse desired image estimate₁After modulo, this result is weighted by (1- α), as shown in step 118. The objective function of equation (4) is then generated in step 120 by adding the first and second terms together. However, this objective function is minimized using, for example, a non-linear conjugate gradient method, as shown in step 122. This minimization process continues until the stopping criteria are met. The stopping criteria include, for example, comparing a current estimate of the desired image with estimates of the desired image from previous iterations. Such a stopping criterion has the following form:

wherein the content of the first and second substances,is the value of the (k + 1) th estimate of the desired image at pixel location (i, j), andis the value of the k-th estimate of the desired image at pixel location (i, j).

Referring now to FIG. 2, another implementation of the method of the present invention employs the objective function of equation (9) and begins by initializing the regularization parameter α, as shown at step 200. Then, formula (II)(9) Is initialized as shown in steps 202 and 204, respectively. This process continues in the same manner as described above with reference to steps 102 and 104 in fig. 1. Now, however, the compatible conditions of formula (3) are combined to the third termThe third item is initialized in step 206. First, as shown in step 208, a difference matrix X is generated. As detailed above, the difference matrix X, corresponding to the compatible condition of equation (3), has the following form:

X=AI-Y (16)。

thus, the difference matrix is determined by applying the system matrix a to the estimate of the desired image I and then subtracting the acquired image data Y corresponding to the desired image. Next, L of the difference matrix X is calculated in step 210₂-the square of the mode. At L where a difference matrix X is generated₂After the square of the modulus, the lagrange multiplier λ is determined and the difference matrix X is weighted with this multiplier, as shown in step 212. As described above, the lagrangian multiplier is empirically determined and values selected by the user based on the clinical application to be performed. The objective function of equation (9) is then generated in step 220 by adding the first, second, and third terms together. However, using, for example, a non-linear conjugate gradient method, this objective function is minimized, as shown in step 222. As described above, this minimization process continues until the stopping criteria are met.

Referring now to FIG. 3, yet another implementation of the method of the present invention employs the objective function of equation (11) and begins by initializing the regularization parameter α, as shown at step 300. Subsequently, the first and second terms in the objective function of equation (11) are initialized, as shown in steps 302 and 304, respectively. This process continues in the same manner as described above with reference to steps 102 and 104 in fig. 1. Now, however, the compatible condition of equation (3) and the impression of noise in the imaging system are combined to the third term λ (X)^TDX), the third entry is initialized in step 306. HeadFirst, as shown in step 308, and as described with reference to step 208 in FIG. 2, a difference matrix X is generated. Next, as shown in step 310, a system noise matrix D is generated. The system noise matrix D is a diagonal matrix having matrix elements determined according to:

as has been described above, in the above-mentioned,is the noise variation and is a parametric indication of noise in the imaging system used in practicing the invention. For example, in an x-ray imaging system, noise parametersIs the noise variation associated with the nth x-ray detector. Optionally, in the MR imaging system, the noise parameterIs the estimated noise variation in the nth receiver coil. After generating the system noise matrix D, the following matrix multiplication is performed:

X^TDX (18);

as shown in step 312. The result of this operation is then scaled by the lagrangian multiplier, as shown in step 314. The objective function of equation (11) is then generated in step 320 by adding the first, second, and third terms together. However, this objective function is minimized, for example by a non-linear conjugate gradient method, as shown in step 322. As described above, this minimization process continues until the stopping criteria are met.

Following for the previous image I_pWith proper choice, a desired image I with higher SNR and spatial resolution can be reconstructed. In general, the preceding image I is reconstructed by the PICCS image reconstruction method described above_pIs applied to the desired image I. Thus, how to generate the previous image I_pMay provide benefits to how the image data is originally acquired. For example, in x-ray imaging, lower tube currents are used to reduce the radiation dose applied to the object being imaged. Such a reduction in tube current produces an increase in the amount of noise present in images reconstructed from data acquired with reduced tube current. However, this increase in image noise is mitigated by the method of the present invention, while maintaining the desired level of spatial resolution, thereby allowing for a reduction in the radiation dose applied to the subject with substantially minimal compromise to the quality of the resulting image.

It is noted that the present invention provides image reconstruction methods that are applicable to a plurality of different imaging modalities, including x-ray computed tomography ("CT"), x-ray c-arm imaging, magnetic resonance imaging ("MRI"), positron emission tomography ("PET"), and single photon emission computed tomography ("SPECT"). More particularly, the present invention provides an image reconstruction method that provides an increase in the achievable signal-to-noise ratio ("SNR") in the reconstructed image without a significant reduction in spatial resolution. In this way, tradeoffs can be made with respect to other imaging considerations, such as the radiation dose applied to the subject and the overall scan time, without a significant reduction in SNR compared to currently available image reconstruction methods.

As mentioned above, the present invention is applicable to many different medical imaging modalities and in many different clinical applications. Several such exemplary clinical applications are described below to illustrate the broad scope of the invention. Such embodiments do not necessarily represent the full scope of the invention, and reference is therefore made to the claims herein for interpreting the scope of the invention.

Referring first to fig. 4A and 4B, an x-ray computed tomography ("CT") imaging system 410 includes a gantry 412 representing a "third generation" CT scanner. Gantry 412 has an x-ray source 413 that projects a fan beam, or cone beam, of x-rays 414 toward a detector array 416 located on the opposite side of the gantry. The detector array 416 is formed by a plurality of detector elements 418 that together sense the projected x-rays that pass through the medical patient 415. Each detector element 418 produces an electrical signal that represents the intensity of an incident x-ray beam and the resulting attenuation of the beam as it passes through the patient. During a scan to acquire x-ray projection data, gantry 412 and the components mounted thereon rotate about a center of rotation 419 located within patient 415.

Rotation of the gantry and operation of x-ray source 413 are governed by a control mechanism 420 of the CT system. Control mechanism 420 includes an x-ray controller 422 that provides power and timing signals to x-ray source 413 and a gantry motor controller 423 that controls the rotational speed and position of gantry 412. A data acquisition system ("DAS") 424 in control mechanism 420 samples analog data from detector elements 418 and converts the data to digital signals for subsequent processing. An image reconstructor 425 receives sampled and digitized x-ray data from DAS 424 and performs high speed image reconstruction. The reconstructed image is applied as an input to a computer 426, which stores the image in a mass storage device 428.

The computer 426 receives commands and scanning parameters from an operator via a console 430 having a keyboard. An associated display 432 allows the operator to observe the reconstructed image and other data from computer 426. The operator supplied commands and parameters are used by the computer 426 to provide control signals and information to the DAS 424, the x-ray controller 422, and the gantry motor controller 423. In addition, computer 426 operates a table motor controller 434, which controls a motorized table 436 to place patient 415 in gantry 412.

Referring now specifically to fig. 5A and 5B, embodiments of the present invention employ an x-ray system designed for use with interventional surgical instruments. Characterized as a frame having a C-arm 510 carrying an x-ray source assembly 512 on one end thereof and an x-ray detector array assembly 514 on the other end thereof. The gantry enables the x-ray source 512 and the detector 514 to be oriented at different positions and angles around a patient positioned on a table 516, while allowing physical contact with the patient.

The gantry includes an L-shaped base 518 having a horizontal leg 520 extending below the table 516, and a vertical leg 522 extending upwardly from the end of the horizontal leg 520 spaced from the table 516. A support arm 524 is rotatably secured to the upper end of the vertical leg 522 for rotation about a horizontal pivot axis 526. The pivot shaft 526 is aligned with the centerline of the table 516, and the arms 524 extend radially outward from the pivot shaft 526 to support the C-arm drive assembly 527 at its outer extremities. The C-arm 510 is slidably secured to a drive assembly 527 and coupled to a drive motor (not shown) that causes the C-arm 510 to slide to orbit about a C-axis 528 as indicated by arrow 530. The pivot axis 526 and the C-axis 528 intersect each other at an isocenter (isocenter) 536 located on the table 516, and they are perpendicular to each other.

An x-ray source assembly 512 is mounted at one end of the C-arm 510 and a detector array assembly 514 is mounted to the other end thereof. As will be discussed in greater detail below, x-ray source 512 emits a cone-shaped x-ray beam that is directed toward detector array 514. Both assemblies 512 and 514 extend radially inward toward pivot axis 526 so that the central ray of this cone beam passes through the system isocenter 536. During acquisition of x-ray attenuation data from an object placed on the table 516, the central ray of the cone beam may thus be rotated about the system isocenter about the pivot axis 526 or the C-axis 528, or both.

As shown in FIG. 6A, x-ray source assembly 512 includes an x-ray source 532 that emits a cone-shaped x-ray beam 533 when energized. The central ray 534 passes through the system isocenter 536 and impinges on a two-dimensional flat-panel digital detector 538 housed in the detector assembly 514. The detector 538 is a two-dimensional array of 2048 by 2048 element detector elements, 41 cm by 41 cm in size. Each element produces an electrical signal that is representative of the intensity of incident x-rays and, therefore, represents the attenuation of the rays caused by the x-rays as they pass through the patient. During a scan, x-ray source 532 and detector array 538 acquire x-ray attenuation projection data from different angles around system isocenter 536. The detector array can acquire 30 projections or views per second, which is a limiting factor in determining how many views can be acquired for a predetermined scan path and speed.

With particular reference to FIG. 6B, the rotation of components 512 and 514 and the operation of x-ray source 532 are controlled by a control mechanism 540 of the CT system. Control mechanism 540 includes an x-ray controller 542 that provides power and timing signals to x-ray source 532. A data acquisition system ("DAS") in control mechanism 540 samples data from detector elements 538 and communicates the data to an image reconstructor 545. An image reconstructor 545 receives digitized x-ray data from DAS 544 and performs high speed image reconstruction in accordance with the methods of the present invention. The reconstructed image is applied as an input to a computer 546, which stores the image in a mass storage device 549 and further processes the image.

Control mechanism 540 also includes a pivot motor controller 547 and a C-axis motor controller 548. In response to motion commands from computer 546, motor controllers 547 and 548 provide power to motors in the x-ray system, resulting in rotation about pivot axis 526 and C-axis 528, respectively. A program executed by computer 546 generates motion commands to motor drivers 547 and 548 to move components 512 and 514 in predetermined scan paths.

Referring now to FIG. 7, a flowchart illustrating steps of an exemplary reconstruction method in accordance with the present invention is shown. The image reconstruction method is generally applicable to x-ray imaging systems and is described herein with reference to x-ray CT systems and x-ray C-arm systems, in which image data is acquired by simultaneously moving an x-ray source and a detector assembly around an object being imaged.

The method begins with the acquisition of image data, as shown in step 700. This image data is acquired, for example, as a plurality of projection space data sets, each acquired by operating the x-ray source and detector to rotate around the object being imaged. For example, each projection space data set corresponds to an image slice position, such that the plurality of projection space data sets together correspond to a plurality of images of the object associated with a respective plurality of image slice positions. Further, the acquired image data is undersampled, thereby significantly reducing the data acquisition time for each projection space data set, and indeed the image data as a whole. In addition, the number of measured x-ray projections is reduced due to the use of undersampling. This, in turn, achieves a significant reduction in the x-ray dose applied to the object. As will be described in detail below, the achievable reduction in x-ray dose without a significant reduction in signal-to-noise ratio ("SNR") in the desired image is proportional to the number of slice positions included when generating a prior image of the object.

From the acquired image data, a plurality of image slices are then reconstructed, as shown in step 702. This image reconstruction is a conventional image reconstruction such as a filtered back-projection type of image reconstruction. The image thus reconstructed will include a higher level of noise. However, although these images are rendered inefficiently for clinical use in view of the noise level, they are useful for generating prior images of the subject with high SNR, which in turn are applied to subsequently reconstructed images.

The generation of the previous image according to the invention starts with the calculation of a plurality of weights from the reconstructed image, as shown in step 704. By way of example, for N image slices available for averaging, the simplest weighting scheme that can be used to average the image slices is to apply the same 1/N weight to all N slices. Alternatively, however, the weights may be selected depending on the image pixel values at a given image slice. For example, for N image slices, the weight z for a given slice position is calculated from the information in the reconstructed image slice. The weights are calculated, for example, by comparing each pixel position in an image slice with the corresponding pixel position in the nearest neighboring image slice. For example, on a pixel-by-pixel basis, the following comparisons are made:

|I_z(m,n)-I_z±1(m,n)|≥T(19)

wherein I_z(m, n) is the pixel value of the z-th image slice at pixel position (m, n); i is_z±1(m, n) are pixel values in image slices adjacent to the z-th image slice at the same pixel position (m, n); and T is a threshold, such as 1000Hounsfield units ("HU"). When the comparison in equation (18) is satisfied, the weight w is then applied_z1/(N-1) to the z-th image slice; otherwise, weight w_z1/N to the image slice.

Using the above-described weighting scheme, an average image volume is produced with a lower noise level, but also with reduced spatial resolution in the z-direction. This image volume is used as a prior image in the PICCS reconstruction algorithm, along with the raw projection data, to reconstruct each individual image slice. As a result of the higher SNR of the previous image, the final image will have a lower noise level, but the z-resolution will be restored, since the original projection data and its early (prior) spatial resolution are used to reconstruct the image slices.

An alternative implementation of the above-described image averaging process is to implement the weighting scheme of the projection data directly. For example, using an equivalent weighting scheme, after selecting the number of image slices to be averaged, N, a weighting function of 1/N is applied directly to the projection data along the z-direction in the detector plane.

After the weights are calculated, they are applied to the respective projection space data sets, as shown in step 706. For example, the following calculation is performed: ...

WhereinIs obtained by mapping the image slice position z to_nN (r) projection space data set P (r)_n) Applying the nth weight w_nThe resulting nth weighted projection data set. As described above, the weight w for a given slice position_zAnd may be, for example, 1/(N-1) or 1/N. Further, it is apparent that other weight values may also be implemented. The processing of weighting the projection space data sets in this way achieves the effect of reducing the contribution of these projection space data sets which contain significantly different signal information than the other projection space data sets. Projection space data set that is then weightedAre combined to form an averaged image data set, as shown in step 708. For example, an averaged image data setIs generated according to:

the averaged image data set is referred to because the weighting and adding together of the projection space data sets as described above has the effect of being similar to averaging the projection space data sets together. From another perspective, this process is one process of filtering in a direction orthogonal to the image slice (such as along the z-axis when the image slice is a transverse, or axial, image).

From the averaged image data set, a previous image I of the object_pIs reconstructed as shown in step 710. This image reconstruction is a conventional image reconstruction such as a filtered back-projection type of reconstruction. Since this averaged image dataset comprises contributions of signal information from a plurality of different projection space datasets, the preceding image I_pIncluding a higher SNR. This previous image I is then used, as shown in step 712_pTo reconstruct the desired object of the object. For example, PICCS image reconstruction is generally performed according to any of the methods shown in fig. 1-3 by selecting an estimated image, generating a thinned-out image using the selected estimated image and an over-the-front image, and reconstructing a desired image of the object using the thinned-out image, the selected estimated image, and the originally acquired image data. Further details regarding PICCS image reconstruction are disclosed, for example, in co-pending U.S. patent application serial No.12/248,590, which is incorporated herein by reference in its entirety.

By way of example, and with particular reference now to FIG. 8, the image reconstruction process described above is illustrated in an exemplary data diagram. Acquired projection data 802 corresponding to, for example, five slice positions is reconstructed to produce an associated series of five image slices 804, which reconstruction is generally illustrated as arrow 806. From the reconstructed image slice 804, and in accordance with the method described above, a set of five weights 808 is calculated, the calculation of which is illustrated generally as arrow 810. As described above, for each image slice position z, a different weight value w is calculated_z. These weights are then applied to the acquired projection data 802 to produce an averaged projection data set from which the prior image 812 was reconstructed. The process of averaging the projection data and reconstructing the prior image is generally illustrated by arrow 814. Using the prior image, the PICCS algorithm is then implemented to reconstruct a series of target image slices 816, the target image slices 816 incorporating the higher SNR of the prior image 812 while exhibiting an appreciable spatial resolution not found in the prior image 812. The PICCS algorithm is generally illustrated as arrow 818.

By implementing a previous image, such as the one described above, in the PICCS image reconstruction framework, the higher SNR of the previous image is applied to the desired image reconstructed from the PICCS image reconstruction. In this way, images with higher SNR than can be achieved using PICCS alone can be achieved while still maintaining the spatial resolution benefits obtained by PICCS image reconstruction. That is, the noise variation σ of the desired image²With respect to preceding image I_pAnd the spatial resolution of the desired image is related to the selected estimated image, which is generally significantly higher than the previous image. Applying the SNR from the previous image while maintaining the spatial resolution of the estimate from the selected image breaks the traditional relationship between SNR and spatial resolution:

is for two-dimensional imaging, an

Is for three-dimensional imaging, where Δ x is the spatial resolution.

The method of the invention thus allows the reconstruction of images with a higher SNR than that obtainable by previous image reconstruction methods, such as conventional filtered backprojection, without suffering a reduction in spatial resolution. Due to this ability to effectively increase the SNR of the reconstructed image, the tube current and mAs of the x-ray CT system, the x-ray C-arm system, etc. can be reduced, so that an image of the object having a desired SNR and spatial resolution can be obtained while significantly reducing the x-ray dose applied to the object.

It is noted that in general the number N of projection space data sets used for generating the averaged image data set corresponds to the increase in SNR obtainable using the method of the invention. Furthermore, the number N of projection space data sets used to generate the averaged image data set is related to the extent to which the tube current can be reduced without a significant reduction in SNR compared to the SNR obtainable using conventional image reconstruction methods and without reducing the tube current. For example, it is conceivable that the relationship between SNR and tube current shown in equation (1) may be modified as follows:

with particular reference to FIG. 9, the preferred embodiment of the present invention is employed in a magnetic resonance imaging ("MRI") system. The MRI system includes a workstation 910 having a display 912 and a keyboard 914. Workstation 910 includes a processor 916, which is a commercially available programmable machine running a commercially available operating system. The workstation 910 provides an operator interface that allows scan specifications to be entered into the MRI system. The workstation 910 is coupled to four servers: a pulse sequence server 918; a data acquisition server 920; a data processing server 922, and a data storage server 923. Workstation 910 and each of servers 918, 920, 922, and 923 are connected to communicate with each other.

The pulse sequence server 918 is operative to operate the gradient system 924 and the radio frequency ("RF") system 926 in response to instructions downloaded from the workstation 910. Gradient waveforms required to perform a prescribed scan are generated and applied to a gradient system 924 which excites the gradient coils in assembly 928 to produce the magnetic field gradients Gx, Gy and Gz which are used to positionally encode the MR signals. The gradient coil assembly 928 forms part of a magnet assembly 930, the magnet assembly 930 including a polarizing magnet 932 and a whole-body RF coil 934.

RF excitation waveforms are applied to the RF coil 934 by the RF system 926 to perform a prescribed magnetic resonance pulse sequence. Responsive MR signals detected by the RF coil 934, or a separate local coil (not shown in fig. 9), are received, amplified, demodulated, filtered, and digitized by the RF system 926 under the direction of commands generated by the pulse sequence server 918. The RF system 926 includes an RF transmitter for generating the various RF pulses used in the MR pulse sequence. The RF transmitter responds to the scan prescription and instructions from the pulse sequence server 918 to generate RF pulses of the desired frequency, phase and pulse amplitude waveform. The generated RF pulses may be applied to a whole-body RF coil 934 or one or more local coils or coil arrays (not shown in fig. 9).

The RF system 926 also includes one or more RF receiver channels. Each RF receiver channel includes an RF amplifier that amplifies the MR signals received by the coil to which it is connected, and a detector that detects and digitizes the I and Q quadrature components of the received MR signals. The magnitude of the received MR signal can thus be determined at any sample point by the square root of the sum of the squares of the I and Q components:

and the phase of the received MR signals may also determine:

the pulse sequence server 918 can also optionally receive patient data from the physiological acquisition controller 936. The controller 936 receives signals from a number of different sensors connected to the patient, such as ECG signals from the electrodes or respiration signals from the bellows (bellows). These signals are typically used by the pulse sequence server 918 to synchronize or "gate" the execution of the scan with the subject's breathing or heartbeat.

The pulse sequence server 918 is also connected to a scan room interface circuit 938 that receives signals from various sensors associated with the condition of the patient and the magnet system. Additionally, through the scan room interface circuit 938, the patient positioning system 940 receives commands to move the patient to desired positions during the scan.

The digitized MR signal samples produced by the RF system 926 are received by the data acquisition server 920. The data acquisition server 920 operates in response to instructions downloaded from the workstation 910 to receive real-time MR data and provide buffer storage so that data is not lost due to data spills. In some scans, the data acquisition server 920 simply passes the acquired MR data to the data processing server 922. However, in scans where information derived from the acquired MR signals is needed to control further performance of the scan, the data acquisition server 920 is programmed to generate and transmit such information to the pulse sequence server 918. For example, during a pre-scan, MR data is acquired and used to calibrate the pulse sequence performed by the pulse sequence server 918. In addition, navigator signals may be acquired during the scan and used to adjust RF or gradient system operating parameters or to control the view sequence in which k-space is sampled. Also, the data acquisition server 920 may be used to process MR signals used to detect the arrival of contrast agent in a Magnetic Resonance Angiography (MRA) scan. In all of these examples, the data acquisition server 920 acquires MR data and processes the data in real time to generate information for controlling the scan.

The data processing server 922 receives MR data from the data acquisition server 920 and processes the data according to instructions downloaded by the workstation 910. Such processing may include, for example: fourier transforming the raw k-space MR data to produce a two or three dimensional image; applying a filter to the reconstructed image; performing a backprojection image reconstruction of the acquired MR data; calculating a functional MR image; and computing a motion or flow image.

The images reconstructed by the data processing server 922 are transmitted back to the workstation 910 where they are stored. The real-time image is stored in a database memory cache (not shown) from which it can be output to the operator display 912 or display 942 located near the magnet assembly 930 for use by the attending physician. Either the batch mode images or the selected real-time images are stored in a host database on disk storage 944. When these images have been reconstructed and transferred to storage, the data processing server 922 notifies the data storage server 923 on the workstation 910. The workstation 910 may be used by an operator to archive the images, produce films, or send the images to other facilities over a network.

An exemplary pulse sequence used to guide an MRI system to acquire image data in accordance with the present invention is illustrated in fig. 10. Such an exemplary pulse sequence is commonly referred to as gradient echo ("GRE") and samples k-space along a series of radial projections, such as extending outward from the center of k-space to the peripheral portion in k-space defined by a k-space radius R. The pulse sequence includes a radio frequency ("RF") excitation pulse 1000 that is generated when a slice-selected gradient 1002 occurs. By sampling in a single plane in k-space, a single two-dimensional slice can be acquired using this pulse sequence, and multiple planes in k-space can be employed using this pulse sequence.

For example, and referring now to FIG. 11, k-space may be sampled along a plurality of radial projections, or radial trajectories 1110, located in a first k-space plane 1102. As described above, the radial projection extends outward from the middle of k-space to the peripheral portion in k-space defined by the k-space radius R, as shown by dashed line 1104.

Referring again to fig. 10, when acquiring multiple two-dimensional slices, the gradient of slice selection 1002 is a lobe-selective gradient (slab-selective gradient), followed by a phase encoding gradient lobe 1010 and an inverted gradient lobe 1012, opposite in polarity to the phase encoding lobe 1010. A phase encoding gradient lobe 1010 and an inverted gradient lobe 1012 are made (play out) along the same gradient axis as the lobe selection gradient 1002. During subsequent iterations of the pulse sequence, the phase encoding gradient 1010 is incremented through (step through) a plurality of different values, thereby sampling k-space in each of a plurality of k-space planes.

In this way, different planes in k-space, such as planes 1106, 1108, 1110, and 1112 in FIG. 11, are sampled. Such planes are substantially parallel and are disposed along a common longitudinal axis 1114. The direction of the longitudinal axis 1114 is defined by the lobe selection gradient 1002. For example, if the lobe selection gradient 1002 is made along the Gz gradient axis, the longitudinal axis 1114 will lie along the z-axis and each sampled k-space plane will correspond to an image slice along the z-axis. That is, the image will be acquired as a transverse slice through the object. Those skilled in the artThe person will readily understand that the lobe-selection gradient 1002 may similarly be along G_xGradient axis, G_yThe gradient axis, or a combination of the three basic gradient axes, is made so that the image is acquired as a sagittal, coronal, or oblique slice, respectively.

Two in-plane read gradients 1014 and 1016 are made during acquisition of nuclear magnetic resonance ("NMR") echo signals 1018 to employ k-space in a given plane (1102, 1106, 1108, 1110, 1112) along the radial trajectory 1100. These in-plane gradients 1014 and 1016 are perpendicular to the lobe selection gradient 1002 and they are perpendicular to each other. During subsequent iterations of the pulse sequence, the in-plane readout gradients 1014 and 1016 are incremented through a series of values to rotate the view angle of the radial sampling trajectory 1100, as described in detail below. In this way, all of the series of values of the in-plane readout gradients 1014 and 1016 are employed before the phase encoding gradient lobe 1010 is incremented to subsequent values, such that a plane in k-space is sampled along multiple radial trajectories 1100 before moving to a different plane in k-space. Each of the in-plane readout gradients is preceded by a pre-phasing gradient lobe 1020 and 1022 and followed by an inverted gradient lobe 1024 and 1026.

It will be appreciated by those skilled in the art that several sampling trajectories other than the one described above may be employed. For example, a straight-line trajectory extending from one point on the k-space peripheral boundary, through the center of k-space, and to an opposite point on the k-space peripheral boundary may also be used. Another variant, similar to the straight-line projection reconstruction pulse sequence, is to sample along a curved path instead of a straight line. Such pulse sequences are described, for example, in F.E.Boada et al, "Fast Three Dimensional Sound Imaging," magnetic resonance in Medicine,1997; (37): 706-. It will be further appreciated by those skilled in the art that cartesian sampling patterns may also be employed to practice the present invention.

The above-described MRI systems may be used in a variety of clinical applications to acquire image data that is used to reconstruct one or more images. The image reconstruction method of the present invention is particularly useful when employing a so-called "parallel imaging" method, in which the signal-to-noise ratio ("SNR") of the resulting image is reduced according to:

where g is the so-called geometric form factor, or "g-factor", and R is the so-called "acceleration factor". In general, the g-factor is a coil (coil) dependent noise amplification factor. In general, when a parallel array of RF receiver coils is implemented to acquire image data, an acceleration R equal to the number of coils in the array can be achieved. This acceleration further corresponds to a reduction in the amount of k-space employed in the acquired image data. For example, an acceleration of R =4 may be achieved in a coil array having at least four receiver coils, and image data acquired with such an acceleration is acquired by undersampling 25% of the samples required by the k-space nyquist criterion. The increase in SNR of the SNR penalty using parallel imaging techniques can be compensated for with the method of the present invention, as described below.

Referring now to fig. 12, a flow chart illustrating steps of an exemplary magnetic resonance image reconstruction method in accordance with the present invention is shown. The method begins with the acquisition of image data, as shown in step 1200. This image data is acquired, for example, as a plurality of k-space datasets, each of which is acquired by directing the MRI system to perform a pulse sequence (such as the one illustrated in fig. 10) to sample k-space. Those skilled in the art will appreciate that k-space may be sampled in any of a variety of different ways, such as along a radial trajectory, a helical trajectory, or a cartesian trajectory. For example, each k-space data set corresponds to an image slice position, such that a plurality of k-space data sets are put together corresponding to a plurality of images of the object associated with a respective plurality of image slice positions. Further, the acquired image data is undersampled, thereby significantly reducing the data acquisition time for each k-space data set, and indeed the image data as a whole.

Note that image data is acquired using a radio frequency ("RF") receiver coil array. As a result, according to equation (27), the SNR obtainable by reconstructing an image using a conventional reconstruction method can be significantly reduced. However, by implementing the methods described herein, the SNR of the desired image is significantly increased without affecting the achievable spatial resolution.

From the acquired image data, a plurality of images are then reconstructed, as shown in step 1202. This image reconstruction is a conventional image reconstruction and depends on the way k-space is sampled. For example, when sampling k-space along a series of radial trajectories, image reconstruction may be performed using a filtered backprojection type method, or by re-meshing the radial samples into a Cartesian grid and performing a Fourier transform. Such reconstructed images include artifacts, such as aliasing, as a result of undersampling employed in acquiring the image data. Thus, although these images are rendered inefficiently for clinical use in view of artifacts, they are useful for generating prior images of subjects with high SNR, which in turn are applied to subsequently reconstructed images.

The generation of the previous image according to the invention starts with the calculation of a plurality of weights from the reconstructed image, as shown in step 1204. Determining weights for a given slice position means that weights are calculated from information in a given reconstructed image, but they are subsequently applied to the associated k-space data set.

As described above, the same 1/N weight may be applied to each of the N reconstructed image slices; however, the weights may also be calculated based on the comparison of pixels in the current image slice with their neighbors. For example, the weight z for a given slice position is calculated from the information in the N reconstructed image slices by comparing each pixel position in the image slice with the corresponding pixel position in the nearest neighboring image slice. For example, on a pixel-by-pixel basis, the following comparisons are made:

|I_z(m,n)-I_z±1(m,n)|≥T（28）

wherein I_z(m, n) is the pixel value of the z-th image segment at pixel position (m, n); i is_z±1(m, n) are pixel values in image segments adjacent to the z-th image segment at the same pixel position (m, n); and T is a threshold. When the comparison in equation (28) is satisfied, the weight w is then applied_z1/(N-1) to the z-th image slice; otherwise, weight w_z1/N to the image slice.

After the weights are calculated, they are applied to the respective k-space data sets, as shown in step 1206. For example, the following calculation is performed: ...

WhereinIs obtained by mapping the image slice position z to_nOf the nth k-space data set S (k)_n) Applying the nth weight w_nThe resulting nth weighted k-data set. The processing of weighting each k-space data set in this way achieves the effect of reducing the contribution of these k-space data sets containing significantly different signal information than the other k-space data sets. Projection space data set that is then weightedAre aggregated to form an averaged image data set, as shown in step 1208. For example, an averaged image dataset is generated according to

The averaged image data sets are so called because the effect of weighting the k-space data sets and adding them together as described above is similar to averaging the k-space data sets together. From another perspective, this process is one process of filtering in a direction orthogonal to the image segment (such as along the z-axis when the image segment is a landscape, or axial, image).

From the averaged image data set, a previous image I of the object_pIs reconstructed as shown in step 1210. This image reconstruction is a conventional image reconstruction and depends on the way k-space is sampled. For example, when sampling k-space along a series of radial trajectories, image reconstruction may be performed using a filtered backprojection type method, or by re-meshing the radial samples into a Cartesian grid and performing a Fourier transform. Since this averaged image dataset comprises contributions of signal information from a plurality of different k-space datasets, the preceding image I_pIncluding a higher SNR. This previous image I is then used, as shown in step 1212_pTo reconstruct the desired object of the object. For example, according to any of the methods shown in fig. 1-3, PICCS image reconstruction is generally performed by selecting an estimated image, generating a thinned-out image using the selected estimated image and a prior image, and reconstructing a desired image of the object using the thinned-out image, the selected estimated image, and the originally acquired image data. Further details regarding PICCS image reconstruction are disclosed, for example, in co-pending U.S. patent application serial No.12/248,590, which is incorporated herein by reference in its entirety.

By implementing a previous image, such as the one described above, in the PICCS image reconstruction framework, the higher SNR of the previous image is applied to the desired image reconstructed from the PICCS image reconstruction. In this way, a map with a higher SNR than can be achieved using PICCS alone can be achievedLike while still maintaining the spatial resolution benefits obtained through PICCS image reconstruction. That is, the noise variation σ of the desired image²With respect to preceding image I_pAnd the spatial resolution of the desired image is related to the selected estimated image, which is generally significantly higher than the previous image. The application of SNR from a previous image, while maintaining the spatial resolution of the estimate from the selected image, breaks the traditional relationship between SNR and spatial resolution:

is for two-dimensional imaging, an

Is for three-dimensional imaging, where Δ x is the spatial resolution.

Thus, the method of the present invention allows the reconstruction of images from undersampled image data having a higher SNR than obtainable by previous image reconstruction methods without sacrificing the reduction in spatial resolution. Due to the ability to effectively increase the SNR of the reconstructed images, a parallel MRI method can be implemented without a regular loss of SNR. In this way, the scan time can be significantly reduced while maintaining the desired image quality. It is noted that in general the number N of k-space data sets used for generating the averaged image data set corresponds to the increase in SNR obtainable using the method of the invention. For example, it is contemplated that the SNR relationship described in equation (27) above may be modified as follows:

the present invention has been described in terms of one or more preferred embodiments, and it is understood that many equivalents, alternatives, variations, modifications, and the like from among these direct representations are possible and are within the scope of the present invention.

Claims (13)

1. A method for reconstructing an image of an object with an x-ray Computed Tomography (CT) system, the steps of the method comprising:

a) acquiring image data from the subject with the CT system;

b) reconstructing a plurality of images from the acquired image data, the plurality of images corresponding to a respective plurality of image slices;

c) calculating a plurality of weights from the reconstructed plurality of images;

d) generating weighted image data by applying the calculated plurality of weights to the acquired image data, thereby filtering the acquired image data in a direction orthogonal to the orientation of the plurality of image slices;

e) reconstructing a prior image from the weighted image data;

f) selecting an estimated image of the object;

g) generating a thinned image of the object using the prior image and the estimated image; and

h) reconstructing a desired image of the object using the thinned image, the estimated image, and the acquired image data.

2. The method of claim 1, wherein the acquired image data includes a plurality of projection space data sets, each projection space data set corresponding to one of the respective plurality of image slices.

3. The method of claim 2, wherein each of the plurality of weights calculated in step c) corresponds to a respective one of the plurality of projection space data sets, and a selected number of weights are calculated.

4. The method of claim 3, wherein step d) includes applying each of the plurality of calculated weights to a corresponding one of the plurality of projection space data sets, thereby producing a respective plurality of weighted projection space data sets.

5. The method of claim 4, wherein step e) comprises generating an average image dataset by summing the plurality of weighted projection space datasets.

6. The method of claim 4, wherein a selected number of weights calculated corresponds to the number of projection space data sets filtered in step d).

7. The method of claim 3, wherein the selected number of weights calculated is proportional to a reduction in x-ray dose applied to the object during step a) obtained without significantly degrading signal-to-noise ratio in the desired image.

8. The method of claim 1, wherein step c) comprises comparing a difference between one of the reconstructed plurality of images and another of the reconstructed plurality of images adjacent to the one of the reconstructed plurality of images to a threshold.

9. The method of claim 8, wherein the reconstructed plurality of images includes a number of images, and when the difference is greater than or equal to the threshold, one of the plurality of weights corresponding to one of the reconstructed plurality of images is assigned to a value less than the number of images.

10. The method of claim 8, wherein the threshold is 1000Hounsfield units.

11. A method for reconstructing an image of a subject with a medical imaging system, the steps of the method comprising:

a) acquiring, with the medical imaging system, a plurality of image datasets from the subject, each of the plurality of image datasets corresponding to a respective one of a plurality of image slice locations;

b) generating an average image data set by averaging a plurality of image data sets in a direction orthogonal to the orientation of the plurality of image slice positions;

c) reconstructing a previous image of the object from the averaged image data set;

d) selecting an estimated image of the object;

e) generating a thinned image of the object using the prior image and the estimated image; and

f) reconstructing a desired image of the subject using the thinned image, the estimated image, and the acquired plurality of image data sets.

12. The method of claim 11, wherein the medical imaging system is at least one of the group of: x-ray computed tomography systems, x-ray C-arm imaging systems, positron emission tomography systems, single photon emission computed tomography systems, and magnetic resonance imaging systems.

13. The method of claim 11, wherein step b) comprises:

i) reconstructing a plurality of images from the acquired plurality of image data sets, each of the reconstructed plurality of images corresponding to a respective one of the plurality of image slice locations;

ii) calculating a plurality of weights from the reconstructed plurality of images;

iii) generating a plurality of weighted image datasets by applying the calculated plurality of weights to the acquired plurality of image datasets, thereby filtering the acquired image datasets in a direction orthogonal to the orientation of the plurality of image slice positions; and

iv) generating an average image data set by adding the generated plurality of weighted image data sets.