Method and system for multi-material decomposition

Abstract

The present invention is entitled "Methods and Systems for Multi-Material Decomposition". The present invention provides various methods and systems for multi-material decomposition for computed tomography. In one embodiment, a method includes: acquiring projection data of multiple x-ray spectra via an imaging system, estimating path lengths for multiple materials based on the projection data and calibration data for the imaging system, iteratively refining the estimated path lengths based on a linearized model derived from the calibration data, and reconstructing a material density image for each of the multiple materials from the iteratively refined estimated path lengths. By determining the path length estimates in this manner without modeling the physics of the imaging system, accurate material decomposition can be performed more quickly and with less sensitivity to changes in the physics of the system, and further can be extended to more than two materials.

Description

Technical Field

[0001] Embodiments of the subject matter disclosed herein relate to computed tomography (CT) imaging, and more particularly to multi-material decomposition for CT imaging. Background Art

[0002] Dual- or multi-spectral computed tomography (CT) imaging systems can display the density of different materials in an object and generate images corresponding to multiple monochromatic x-ray energy levels. The CT imaging system can derive the behavior at different monochromatic energy levels based on signals from at least two regions of photon energy in the spectrum (e.g., the low-energy portion and the high-energy portion of the incident x-ray spectrum). In a given energy region of medical CT, where the scanned object is a patient, the following two physical processes dominate the x-ray attenuation process: Compton scattering and the photoelectric effect. The signals detected from the two energy regions provide sufficient information to resolve the energy dependence of the imaged material. The signals detected from the two energy regions provide sufficient information to determine the relative composition of an object composed of two hypothetical materials. Summary of the invention

[0003] In one embodiment, a method includes: acquiring projection data of a plurality of x-ray spectra via an imaging system; estimating path lengths for a plurality of materials based on the acquired projection data and calibration data for the imaging system; and reconstructing a material density image for each of the plurality of materials from the estimated path lengths. Hereinafter, the term "path length" is used to designate a line integral of the material density image along a line connecting an x-ray source to respective detector elements. By determining the path length estimate in this manner without modeling the physics of the imaging system, accurate material decomposition can be performed more quickly and with less sensitivity to changes in the physics of the system, and furthermore can be extended to more than two materials.

[0004] It should be understood that the above brief description is provided to introduce in a simplified form selected concepts that are further described in the detailed description. It is not meant to identify key or essential features of the claimed subject matter, the scope of which is solely defined by the claims that follow the detailed description. Furthermore, the claimed subject matter is not limited to implementations that solve any disadvantages noted above or in any part of this disclosure. BRIEF DESCRIPTION OF THE DRAWINGS

[0005] The present invention will be better understood by reading the following description of non-limiting embodiments with reference to the accompanying drawings, in which:

[0006] Figure 1 shows a pictorial view of an imaging system according to one embodiment;

[0007] Figure 2shows a block schematic diagram of an exemplary imaging system according to one embodiment;

[0008] Figure 3 A high-level flow chart illustrating an exemplary method for multi-material decomposition is shown according to one embodiment;

[0009] Figure 4 A high-level flowchart illustrating an exemplary method for calculating preliminary estimates of path lengths for a variety of materials based on calibration data is shown, according to one embodiment;

[0010] Figure 5 shows a set of graphs depicting an exemplary multi-dimensional surface for initially estimating path lengths for various materials according to one embodiment;

[0011] Figure 6 shows a set of graphs depicting an exemplary linear model for refining initial estimates of path lengths for various materials, according to one embodiment;

[0012] Figure 7 shows a graph showing an exemplary estimated sinusoid and a true sinusoid of a material according to one embodiment;

[0013] Figure 8 Shows the description Figure 7 a graph of an enlarged view of a graph of;

[0014] Fig. 9 An embodiment of the invention is shown Figure 7 a graph of the difference between the estimated sinusoid and the true sinusoid; and

[0015] Fig.10 A high-level flow chart illustrating an exemplary method for estimating a monochromatic sinusoidal pattern is shown, according to one embodiment. DETAILED DESCRIPTION

[0016] The following description relates to various embodiments of spectral computed tomography (CT) imaging. In particular, methods and systems for performing multi-material decomposition in spectral CT imaging are provided. Figure 1 and Figure 2An example of a CT imaging system that can be used to acquire images according to the present techniques is shown in . The CT imaging system may be configured with an energy discriminating detector, such as a photon counting detector, which provides a fidelity capable of distinguishing materials via material decomposition. Material decomposition may be performed in the projection domain, the image domain, or in conjunction with reconstruction. However, joint material decomposition-reconstruction methods are typically computationally expensive compared to both projection domain methods and image domain methods. In addition, projection domain methods for multi-material decomposition (i.e., material decomposition of two or more materials) typically depend on prior knowledge of the imaging system, such as knowledge of some model of the effective x-ray spectrum, which is typically a combination of the x-ray spectrum and the detector spectral response. In implementation, such prior knowledge is difficult to obtain with high accuracy, particularly if the detector spectral response varies across the detector array and / or when the emitted x-ray spectrum varies with the lifetime of the x-ray tube. Methods for multi-material decomposition, such as Figure 3 The method described in , overcomes the challenges of these previous methods by estimating the multi-material path length based only on calibration data and without knowledge of the physics of the multi-energy CT imaging system. The method comprises a two-step approach, in which a preliminary estimate is first obtained based on the calibration data and the acquired projection data, and then the preliminary estimate is used to initialize an iterative scheme that optimizes a statistical criterion to further enhance the statistical accuracy of the multi-material path length. Figure 4 As shown in FIG. 1 , a method for obtaining a preliminary estimate of the multi-material path length includes a combination of function back-lookup and local multilinear fitting methods. The process of calculating such a preliminary estimate is performed using simulated data in Figure 5 and Figure 6 The results of further simulations are shown in Figures 7 to 9 , indicating that the preliminary estimate is highly accurate, allowing the iterative scheme to obtain a final estimate of the multi-material path length with even higher accuracy in only a few iterations. The method described herein for directly estimating the multi-material path length can also be adapted to estimate monochromatic sinusoidal patterns at multiple energies, such as Fig.10 shown.

[0017] Figure 1 An exemplary CT system 100 configured for CT imaging is shown. Specifically, the CT system 100 is configured to image a subject 112 (such as a patient, an inanimate object, one or more manufactured parts) and / or a foreign object (such as a dental implant, a stent, and / or a contrast agent present in the body). In one embodiment, the CT system 100 includes a gantry 102, which in turn may also include at least one x-ray source 104 configured to project an x-ray radiation beam 106 (see Figure 2) for imaging a subject 112 lying on an examination table 114. Specifically, the x-ray source 104 is configured to project an x-ray radiation beam 106 toward a detector array 108 positioned on an opposite side of the gantry 102. Figure 1 Only one x-ray source 104 is depicted, but in certain embodiments, multiple x-ray sources and detectors may be employed to project multiple x-ray radiation beams 106 to acquire projection data at different energy levels corresponding to the patient. In some embodiments, the x-ray source 104 may achieve dual energy spectrum imaging by rapid peak kilovolt voltage (kVp) switching. In some embodiments, the x-ray detector employed is a photon counting detector capable of distinguishing x-ray photons of different energies. In other embodiments, two sets of x-ray sources and detectors are used to generate dual energy projections, wherein one set of x-ray sources and detectors operates at a low kVp and the other set of x-ray sources and detectors operates at a high kVp. Therefore, it should be understood that the methods described herein may be implemented with a variety of multi-spectral acquisition techniques and are not limited to the specific described embodiments.

[0018] In certain embodiments, the CT system 100 further includes an image processor unit 110 configured to reconstruct an image of a target volume of the subject 112 using an iterative or analytical image reconstruction method. For example, the image processor unit 110 may reconstruct an image of the target volume of the patient using an analytical image reconstruction method such as filtered back projection (FBP). As another example, the image processor unit 110 may reconstruct an image of the target volume of the subject 112 using an iterative image reconstruction method such as advanced statistical iterative reconstruction (ASIR), conjugate gradient (CG), maximum likelihood expectation maximization (MLEM), model-based iterative reconstruction (MBIR), etc. As further described herein, in some examples, in addition to iterative image reconstruction methods, the image processor unit 110 may also use analytical image reconstruction methods such as FBP.

[0019] In some CT imaging system configurations, an x-ray source projects a cone-shaped x-ray radiation beam that is collimated to lie in the XYZ plane of a Cartesian coordinate system and is often referred to as the "imaging plane". The x-ray radiation beam passes through an object being imaged, such as a patient or subject. The x-ray radiation beam, after being attenuated by the object, impinges upon an array of detector elements. The intensity of the attenuated x-ray radiation beam received at the detector array depends on the attenuation of the radiation beam by the object. Each detector element of the array produces a separate electrical signal that is a measurement of the x-ray beam attenuation at the detector location. The attenuation measurements from all detector elements are collected individually to produce a transmission profile.

[0020] In some CT systems, a gantry is used to rotate the x-ray source and detector array around the object to be imaged in the imaging plane so that the angle at which the radiation beam intersects the object is constantly changing. A set of x-ray radiation attenuation measurements (e.g., projection data) from the detector array at one gantry angle is referred to as a "view." A "scan" of an object includes a set of views made at different gantry angles or viewing angles during one rotation of the x-ray source and detector. It is contemplated that the benefits of the methods described herein derive from medical imaging modalities other than CT, and therefore, as used herein, the term "view" is not limited to the uses described above with respect to projection data from one gantry angle. The term "view" is used to represent a data acquisition whenever there are multiple data acquisitions from different angles, whether from a CT imaging system or from another imaging system that acquires multi-spectral attenuation measurements.

[0021] The projection data is processed to reconstruct an image corresponding to a two-dimensional slice acquired through the object, or in some examples where the projection data includes data from multiple detector rows or from multiple gantry rotations, to reconstruct an image corresponding to a three-dimensional volume representation of the object. A method for reconstructing an image from a set of projection data is known in the art as a filtered back projection technique. Transmission and emission tomography reconstruction techniques also include statistical iterative methods such as maximum likelihood expectation maximization (MLEM) and ordered subset expectation reconstruction techniques and iterative reconstruction techniques. The method converts attenuation measurements from the scan into "CT numbers" or "Hohn's units," which are used to control the brightness of corresponding pixels on a display device.

[0022] To reduce the total scan time, a "helical" scan may be performed. To perform a "helical" scan, the patient is moved using the examination table 114 while data for a specified number of slices is acquired. Such systems produce a single spiral from a cone beam helical scan. The spiral mapped out by the cone beam produces projection data from which an image in each specified slice can be reconstructed.

[0023] As used herein, the phrase "reconstructing an image" is not intended to exclude embodiments of the invention in which data representing an image is generated rather than a visible image. Thus, as used herein, the term "image" broadly refers to both a visible image and data representing a visible image. However, many embodiments generate (or are configured to generate) at least one visible image. As used herein, the term image may refer to a two-dimensional image as well as a three-dimensional or higher dimensional image volume.

[0024] Figure 2 Shows something like Figure 1 1. An exemplary imaging system 200 of the CT system 100 of FIG. According to aspects of the present disclosure, the imaging system 200 is configured to image a subject 204 (eg, Figure 1 In one embodiment, the imaging system 200 includes a detector array 108 (see Figure 1 ). The detector array 108 also includes a plurality of detector elements 202 that together sense the x-ray radiation beam 106 that passes through a subject 204 (such as a patient) (see Figure 2 ) to collect corresponding projection data. Therefore, in one embodiment, the detector array 108 is manufactured in a multi-slice configuration including multiple rows of cells or detector elements 202. In such a configuration, one or more additional rows of detector elements 202 are arranged in a parallel configuration for collecting projection data.

[0025] In certain embodiments, the imaging system 200 is configured to traverse different angular positions around the subject 204 to collect the desired projection data. Thus, the gantry 102 and components mounted thereon can be configured to rotate about the rotation center 206 to collect, for example, projection data at different energy levels. Alternatively, in embodiments where the projection angle relative to the subject 204 changes over time, the mounted components can be configured to move along a generally curved line rather than along a section of a circle.

[0026] Thus, as the x-ray source 104 and the detector array 108 rotate, the detector array 108 collects data of the attenuated x-ray beam. The data collected by the detector array 108 then undergoes pre-processing and calibration to condition the data to represent the line integral of the attenuation coefficient of the scanned subject 204. The processed data is often referred to as a projection.

[0027] In some examples, individual detectors or detector elements 202 in the detector array 108 may include energy discriminating photon counting detectors that register interactions of individual photons into one or more energy bins. It should be appreciated that the methods described herein may also be implemented using energy integrating detectors.

[0028] The acquired projection data sets may be used for base material decomposition (BMD). During BMD, the measured projections are converted into a set of material density projections. The material density projections may be reconstructed to form a pair or set of material density maps or images - a map or image of each corresponding base material (such as a bone, soft tissue, and / or contrast agent map). The density maps or images may then be correlated to form a volume rendering of the base material (e.g., bone, soft tissue, and / or contrast agent) in the imaging volume.

[0029] Once reconstructed, the base material image produced by the imaging system 200 shows the internal features of the subject 204 represented by the density of the two base materials. The density image can be displayed to show these features. In traditional methods of diagnosing medical conditions (such as disease states), and more generally diagnosing medical events, a radiologist or physician will consider a hard copy or display of a density image to discern characteristic features of interest. Such features can include lesions, size and shape of specific anatomical structures or organs, as well as other features that should be discernible in the image based on the skill and knowledge of the individual practitioner.

[0030] In one embodiment, the imaging system 200 includes a control mechanism 208 to control movement of components, such as rotation of the gantry 102 and operation of the x-ray source 104. In certain embodiments, the control mechanism 208 also includes an x-ray controller 210 configured to provide power and timing signals to the x-ray source 104. Additionally, the control mechanism 208 includes a gantry motor controller 212 configured to control the rotational speed and / or position of the gantry 102 based on imaging requirements.

[0031] In certain embodiments, the control mechanism 208 further includes a data acquisition system (DAS) 214 configured to sample analog data received from the detector elements 202 and convert the analog data into digital signals for subsequent processing. The DAS 214 may include a component of the control mechanism 208 or a separate component, such as Figure 2 2. The DAS 214 may also be configured to selectively aggregate analog data from a subset of the detector elements 202 into a so-called macro detector, as further described herein. The data sampled and digitized by the DAS 214 is transmitted to a computer or computing device 216. In one example, the computing device 216 stores the data in a storage device or mass storage device 218. For example, the storage device 218 may include a hard drive, a floppy disk drive, a compact disk-read / write (CD-R / W) drive, a digital versatile disc (DVD) drive, a flash drive, and / or a solid-state storage drive.

[0032] In addition, the computing device 216 provides commands and parameters to one or more of the DAS 214, the x-ray controller 210, and the gantry motor controller 212 to control system operations, such as data acquisition and / or processing. In certain embodiments, the computing device 216 controls system operations based on operator input. The computing device 216 receives operator input, which includes, for example, commands and / or scan parameters, via an operator console 220 operably coupled to the computing device 216. The operator console 220 may include a keyboard (not shown) or a touch screen to allow an operator to specify commands and / or scan parameters.

[0033] Although Figure 2 Only one operator console 220 is shown, but more than one operator console may be coupled to the imaging system 200, e.g., for inputting or outputting system parameters, requesting exams, plotting data, and / or viewing images. In addition, in certain embodiments, the imaging system 200 may be coupled to multiple displays, printers, workstations, and / or the like located locally or remotely, e.g., within an institution or hospital or at an entirely different location, via one or more configurable wired and / or wireless networks (such as the Internet and / or a virtual private network, a wireless telephone network, a wireless local area network, a wired local area network, a wireless wide area network, a wired wide area network, etc.).

[0034] In one embodiment, for example, the imaging system 200 includes or is coupled to a picture archiving and communication system (PACS) 224. In an exemplary embodiment, the PACS 224 is further coupled to a remote system (such as a radiology department information system, a hospital information system) and / or to an internal or external network (not shown) to allow operators at different locations to supply commands and parameters and / or gain access to image data.

[0035] The computing device 216 uses operator-provided and / or system-defined commands and parameters to operate a table motor controller 226, which in turn can control the table 114, which can be a motorized table. Specifically, the table motor controller 226 can move the table 114 to properly position the subject 204 in the gantry 102 to acquire projection data corresponding to a target volume of the subject 204.

[0036] As previously described, the DAS 214 samples and digitizes the projection data acquired by the detector elements 202. Subsequently, the image reconstructor 230 uses the sampled and digitized x-ray data to perform high-speed reconstruction. Figure 2The image reconstructor 230 is shown as a separate entity, but in certain embodiments, the image reconstructor 230 may form part of the computing device 216. Alternatively, the image reconstructor 230 may not be present in the imaging system 200, and instead the computing device 216 may perform one or more functions of the image reconstructor 230. Furthermore, the image reconstructor 230 may be located locally or remotely, and may be operably connected to the imaging system 200 using a wired or wireless network. In particular, one exemplary embodiment may use computing resources in a "cloud" network cluster for the image reconstructor 230.

[0037] In one embodiment, the image reconstructor 230 stores the reconstructed image in the storage device 218. Alternatively, the image reconstructor 230 may transmit the reconstructed image to the computing device 216 to generate usable patient information for diagnosis and evaluation. In certain embodiments, the computing device 216 may transmit the reconstructed image and / or the patient information to a display or display device 232 that is communicatively coupled to the computing device 216 and / or the image reconstructor 230. In some embodiments, the reconstructed image may be transmitted from the computing device 216 or the image reconstructor 230 to the storage device 218 for short-term or long-term storage.

[0038] Various methods and processes further described herein (such as those described below in reference Figure 3 , Figure 4 and Fig.10 The methods described herein may be stored as executable instructions in a non-transitory memory on a computing device (or controller) in the imaging system 200. In one embodiment, the image reconstructor 230 may include such executable instructions in the non-transitory memory and may apply the methods described herein to reconstruct an image from the scan data. In another embodiment, the computing device 216 may include instructions in the non-transitory memory and may apply the methods described herein, at least in part, to the reconstructed image after receiving the reconstructed image from the image reconstructor 230. In another embodiment, the methods and processes described herein may be distributed on the image reconstructor 230 and the computing device 216.

[0039] In one embodiment, the display 232 allows the operator to assess the imaged anatomy. The display 232 may also allow the operator to select a volume of interest (VOI) and / or request patient information, such as via a graphical user interface (GUI), for subsequent scanning or processing.

[0040] Figure 3A high-level flowchart illustrating an exemplary method 300 for multi-material decomposition is shown according to one embodiment. Specifically, the method 300 involves performing multi-material decomposition in the projection domain based on calibration data rather than based on knowledge of the physics of the imaging system (such as a model of the incident x-ray spectrum and / or detector response). Figure 1 and Figure 2 The method 300 is described with reference to the systems and components of the present disclosure, but it should be understood that the method 300 can be implemented with other systems and components without departing from the scope of the present disclosure. The method 300 can be implemented, for example, as executable instructions in a non-transitory memory of the computing device 216 and / or the image reconstructor 230, and can be executed by a processor of the computing device 216 and / or the image reconstructor 230 to perform the actions described below.

[0041] The method 300 starts at 305. At 305, the method 300 performs a scan of a subject to acquire intensity measurements and generate projection data. The scan includes a dual energy CT scan or a multi-energy CT scan of the subject. To this end, the method 300 controls the x-ray controller 210 to drive the x-ray source 104 to emit x-rays at two or more energy levels, for example, while also controlling the gantry motor controller 212 and the table motor controller 226 to adjust the positions of the gantry 102 and the table 114, respectively, so that the position of the x-ray source 104 relative to the scanned subject 204 is adjusted while the x-rays are generated. The method 300 also acquires projection data measured by an energy discriminating detector (including a photon counting detector, such as the detector element 202 of the detector array 108), for example, via the DAS 214. The projection data includes multi-energy photon counting measurements, for example, where the photon measurements of different energies are classified into predefined energy intervals.

[0042] At 310, method 300 calculates preliminary estimates of path lengths for a plurality of materials based on calibration data and projection data. The calibration data includes calibration data for calibrating the imaging system, which is acquired with the imaging system during a calibration scan of a phantom including known materials and known path lengths. Method 300 calculates preliminary estimates of path lengths for a plurality of materials based on the calibration data, which will result in observed multi-energy measurements of the projection data. To this end, method 300 may perform a two-step estimation that performs an inverse function lookup on the calibration data to obtain an initial estimate of the multi-material path length that will result in the observed multi-energy measurements, and refines the initial estimate by solving a system of linear equations constructed based on the initial estimate. References herein Figure 4 Further described are exemplary methods for calculating preliminary estimates of path lengths for a variety of materials based on calibration data.

[0043] At 315, the method 300 iteratively calculates a final estimate of the path length for each of the plurality of materials, with the preliminary estimate serving as an initial estimate. For example, the method 300 may calculate the final estimate m by iteratively solving the following equation: final :

[0044]

[0045] m∈Ω m , where Ω m includes a set of (physically) feasible material path length vectors, F is a statistical criterion (e.g., log-likelihood value) that produces (physically) meaningful estimates of multi-material path lengths when optimized, f cal A differentiable forward model is included that uses an interpolated value or a polynomial model to map a known material path length to a corresponding known projection value (p-value) or x-ray intensity (I-value), p meas includes the p-value measured at 305, and I meas includes the x-ray intensity measured at 305. Although the above formula applies to the entire sinogram (ie, all sinogram intervals i), it should be understood that the above formula for the final estimate of m is different depending on the components of F. final The equations for can be separated in the sinogram interval i and can therefore be solved in parallel.

[0046] The above statistical problem can be understood as setting a maximum likelihood (or maximizing a posteriori) estimator for the expected multi-material path length. In one example, the function F may include a Poisson log-likelihood function such that:

[0047]

[0048] The index i denotes the sinogram interval, the index k denotes the energy interval, and the measured x-ray intensity corresponds to an air scan when no material is present, To map a known material path length to a known normalized x-ray intensity A differentiable forward model of and Corresponding to the corrected intensity measurement and air scan in the absence of material, respectively. In another example, function F may include a weighted least squares criterion function such that:

[0049]

[0050] in and proportional, and is a differentiable forward model that maps known path lengths to known p values, and the index i indexes the sinogram interval. The advantage of using the Poisson log-likelihood function of F is that it is accurate over the measured x-ray intensity I meas is lower (or even zero in the case of photon counting detectors, or negative due to electronic noise in the case of energy integrating detectors), which prevents meas -log(·) is applied in meas The amount Corresponding to the sinogram bin index i and energy bin index k, is calculated as:

[0051]

[0052] At 320, the method 300 reconstructs a material density image for each material based on the final estimate of the path length. For example, the final estimate of the path length for each material includes a projection based on the material, so the method 300 reconstructs a material density image for each material from the final path length estimate of the material. At 325, the method 300 outputs the material density image. For example, the method 300 may output the material density image to a display device such as the display device 232. Additionally or alternatively, the method 300 may output the material density image to the mass storage device 218 for storage and / or the PACS 224 for remote viewing.

[0053] At 330, method 300 generates a monochromatic image at a plurality of energies based on the material density images. For example, method 300 may selectively combine the material density images to generate a monochromatic image at a given energy that simulates an acquisition using an x-ray source that emits photons only at a given energy. Method 300 may generate a plurality of monochromatic images at different energies from the material density images. At 335, method 300 outputs the monochromatic image. For example, method 300 may output the monochromatic image to a display device, such as display device 232. Additionally or alternatively, method 300 may output the material density image to mass storage device 218 for storage and / or PACS 224 for remote viewing. Method 300 then returns.

[0054] Figure 4 A high-level flowchart is shown that illustrates an exemplary method 400 for calculating a preliminary estimate of path length for a plurality of materials based on calibration data, according to one embodiment. Specifically, the method 400 involves determining a first estimate of the path length based on the calibration data, and refining the first estimate by solving a system of linear equations that corresponds to a local linear approximation of a forward model at the first estimate. Figure 1 and Figure 2Method 400 is described with reference to systems and components of, but it should be understood that method 400 may be implemented with other systems and components without departing from the scope of the present disclosure. Method 400 may be implemented, for example, as executable instructions in a non-transitory memory of computing device 216 and / or image reconstructor 230, and may be executed by a processor of computing device 216 and / or image reconstructor 230 to perform the actions described herein below. Method 400 may include a subroutine of method 300, and specifically may include an action 310 of calculating a preliminary estimate of the path length of a plurality of materials based on calibration data and projection data.

[0055] The method 400 starts at 405. At 405, the method 400 loads projection data. For example, the method 400 loads the projection data acquired as described above at 305. The projection data includes multi-energy photon counting measurements, for example, where photon measurements of different energies are categorized into predefined energy bins.

[0056] Continuing at 410, method 400 performs an inverse function lookup on the calibration data to obtain a first estimate of the path length corresponding to the projection data. The calibration data includes measurements acquired during a calibration scan using a phantom including known materials and known path lengths to calibrate the imaging system. For example, the calibration data may be provided as a set S cal Stored in non-transitory memory:

[0057]

[0058] Each vector is an Lx1 vector of known path lengths of L material, each vector is a Kx1 vector of known p-values ​​(log-normalized intensity values) corresponding to energy bins k=1,…,K, i is the sinogram bin index, and j is the experiment index, where different experiments correspond to different material combinations. Differentiable forward model for the ith sinogram bin is based on the calibration set S cal The calibration set maps material path lengths m to p values ​​such that using an interpolation such as non-uniform rational B-splines or Lagrangian interpolation,

[0059]

[0060] Or alternatively, using a polynomial model,

[0061]

[0062] Then, for a given K energy bin, log-normalize the measurement vector,

[0063]

[0064] and the corresponding x-ray intensity measurement vector,

[0065]

[0066] Corresponding to sinogram interval i (ie, corresponding to an x-ray projection line), method 400 calculates from calibration data S cal Find possible candidate vectors in This candidate vector produces a close measurement vector The p-value vector For example, method 400 may determine:

[0067]

[0068] Function 1 (.) is an indicator function that is unity (i.e., one) if the conditions in parentheses are met or it is otherwise zero, 3 ≤ K sub ≤K,T p yes to The function dist(a, b) is a function that measures the distance between its parameters a and b. In some examples, the function dist(a, b) may include an absolute difference or a squared difference function. The resulting set of search results,

[0069]

[0070] Thus, a set of material vectors is included from which method 400 obtains a rough estimate of the expected multi-material path length. For example, method 400 calculates an initial estimate by performing a simple weighted sum of the material vectors:

[0071]

[0072] where w n Including giving material path length Higher priority weights, corresponding to the p-value vector Numerically closer to The numerical proximity depends on the vector being measured. The p-value in The credibility of The corresponding x-ray intensity measurements in The weighted model thus provides a statistically relevant estimate by giving higher weights to energy bins with higher detected x-ray intensities.

[0073] At 415, method 400 generates a linear approximation of the forward model based on the first estimate of the path length. For example, method 400 may generate a linear approximation of the forward model based on the first estimate of the path length. Perform forward modeling around calibration data The local multilinear approximation model is used to further refine the first estimate The multilinear approximation model produces a system of linear equations that can be expressed as follows:

[0074]

[0075] The coefficient matrix Contains the coefficients of the linear model, which is obtained by adding the matrix A set of known Solving the above linear equations, we get the matrix corresponding to the matrix The known path lengths of the multilinear basis vectors (i.e., One path length for each material. Note that in the following paragraphs, the vector m i andp i Now is the linearization around the operating point defined by the search step. After determining the coefficient matrix, method 400 then calculates the measured p-value based on the first estimate of the path length. A similar linear system is set up to generate a linear approximation of the forward model, such that:

[0076]

[0077] Continuing at 420, method 400 solves the system of linear equations based on linear approximation to obtain a preliminary estimate of the path length. For example, to solve the above preliminary estimate of the path length The linear system established, method 400 calculates:

[0078]

[0079] The weight matrix With the measured intensity proportionally to give relatively better weighting to energy bins with higher measured x-ray intensities.

[0080] References Figure 6 An illustrative example of a local linear model of a forward model constructed based on a first estimate of the path length is further described.

[0081] At 425, method 400 outputs the preliminary estimate of the path length. For example, method 400 can output the preliminary estimate of the path length to a memory so that method 300 can iteratively calculate a final estimate of the path length for each material, with the preliminary estimate serving as an initial estimate, as described above. Method 400 then returns.

[0082] Figure 5 A set of graphs 500 depicting an exemplary multidimensional surface for initially estimating the path length of a plurality of materials is shown, according to one embodiment. The set of graphs 500 includes a graph for each of a plurality of energy intervals, specifically, K=8 spectral energy intervals in the depicted example, including a first graph 510 for a first energy interval, a second graph 520 for a second energy interval, a third graph 530 for a third energy interval, a fourth graph 540 for a fourth energy interval, a fifth graph 550 for a fifth energy interval, a sixth graph 560 for a sixth energy interval, a seventh graph 570 for a seventh energy interval, and an eighth graph 580 for an eighth energy interval. Each graph in the set of graphs 500 depicts a graph that converts a material path length (e.g., ) is mapped to the log-normalized calibration data (e.g., ) of the forward model (e.g., ) surface curve.

[0083] For example, a first graph 510 shows a first surface 511 of a forward model that maps material path lengths to logarithmically normalized calibration data p for a first energy bin (eg, k=1) of a photon counting CT imaging system. ij1 Specifically, the first surface 511 maps the material path lengths of two materials, water (W) and iodine (I), to calibration data. The black dots on each surface, including the first surface 511, correspond to the measured values ​​along the x-ray path. including measurements 513 (e.g., ). Points 516 and 517 on the [W, I] plane represent possible estimates of the [W, I] path length These estimates are for example using S lookup The criteria are derived from graphs 510, 520, 530, 540, 550, 560, 570 and 580 together, and the corresponding measured values ​​are derived therefrom. A first estimate 515 of the [W, I] path length of the material along the x-ray path (e.g., ) (ie, in the depicted example, an 8x1 vector of p-values), as described above. Point 519 indicates the true or actual path length of the material.

[0084] The preliminary estimate 515 is in the [W, I] plane and represents the value corresponding to the measurement 515 is a graph of the estimated path lengths of water and iodine. Therefore, for illustrative purposes, the preliminary estimate 515 is replicated in all of the graphs 520, 530, 540, 550, 560, 570, and 580. It should be understood that for ease of illustration, two materials, water and iodine, are depicted because the dimensions of the forward model for three or more materials are difficult to show in a two-dimensional representation. However, the method shown is applicable to two or more materials.

[0085] In order to improve Figure 5 The first estimate of the surface obtained by the forward model shown can be used to construct a local joint linear model of the surface. As an illustrative and non-limiting example, Figure 6 A set of graphs 600 are shown, which depict the Figure 5 An exemplary linear model for initial estimation of path length for various materials is depicted in group of graphs 500 .

[0086] The set of graphs 600 includes a first graph 610 for a first energy interval, a second graph 620 for a second energy interval, a third graph 630 for a third energy interval, a fourth graph 640 for a fourth energy interval, a fifth graph 650 for a fifth energy interval, a sixth graph 660 for a sixth energy interval, a seventh graph 670 for a seventh energy interval, and an eighth graph 680 for an eighth energy interval. Each graph in the set of graphs 600 depicts a forward model depicted in the set of graphs 500. surfaces, including a first surface 511.

[0087] Based on Figure 5 For each estimate of each energy interval obtained as shown, a local joint linear model of each surface is constructed, resulting in a system of equations that can be solved to obtain an improved estimate that is better than the first estimate. For example, based on the first estimate 515, a linear model 613 of the surface 511 is constructed to correspond to the region 611 of the [W, I] plane used to construct the linear model 613 as described above. The linear model 613, together with the linear models in 620, 630, 640, 650, 660, 670 and 680, produces a joint system of linear equations that can be solved to find a better preliminary estimate 615 of the path length relative to the first estimate 515 (e.g., ), as shown in graph 610. Note that preliminary estimate 615 is in the [W, I] plane and represents the estimated path lengths of water and iodine. Therefore, for illustrative purposes, preliminary estimate 615 is replicated in all graphs 620, 630, 640, 650, 660, 670, and 680. In particular, preliminary estimate 615 is closer to true condition 519 than first estimate 515.

[0088] Alternatively, a first estimate (e.g., ) and improve estimates (e.g., ) can be achieved by using the intensity domain forward model The model is obtained by appropriately modifying the method 400 described above to use the intensity domain calibration data accordingly. and measurement data Mapping known material path lengths to known x-ray intensity values Operating in the intensity domain may be advantageous, especially if has a poor signal-to-noise ratio, or if It is not right.

[0089] The initial estimate thus obtained is used to initialize the iterative optimizer (e.g., as described above with reference to Figure 3 The invention relates to a method for reconstructing a material-based image by using a method for reconstructing a material-based image. The method further improves the path length estimation by using a method for reconstructing a material-based image. In some examples, since the preliminary estimate of the path length thus obtained is substantially close to the actual situation, the preliminary estimate of the path length can be directly used to reconstruct a material-based image. That is, in some examples, the iterative optimization of the path length estimate initialized with the preliminary estimate can be omitted from the multi-energy material decomposition. In addition, in examples in which iterative optimization is performed, the iterative method converges to the final estimate of the path length with fewer iterations and higher accuracy. Since the steps of performing an inverse function lookup to calculate the first estimate and generating a linear approximation of the forward model based on the first estimate are not computationally burdensome, the overall computational complexity of the material decomposition for multiple energies is reduced. In addition, according to the systems and methods provided herein, accurate material decomposition can be performed without simulating or modeling the physics of the imaging system or the x-ray spectra of multiple energies and multiple materials.

[0090] To illustrate the efficacy of the techniques provided herein for multi-energetic material decomposition, Figure 7 A graph 700 is shown showing exemplary distributions 705 and 710 in an estimated sinogram of a material according to one embodiment and a corresponding true distribution 715. Specifically, graph 700 depicts an exemplary distribution 705 in an estimated sinogram corresponding to a preliminary estimate (e.g., calculated according to method 400 as described above) and an exemplary distribution 710 in an estimated sinogram corresponding to an optimized estimate (e.g., iteratively calculated at 315 as described above).

[0091] As shown, the exemplary distributions 705, 710, and 715 in the estimated sinograms are similar enough to be Figure 7 To illustrate the differences between the exemplary distributions 705, 710, and 715 in the estimated sinograms, Figure 8Graph 800 is shown, which shows a magnified view of region 725 of graph 700. In this view, the deviation of the preliminary estimate depicted by the exemplary distribution 705 in the estimated sinogram from the exemplary distribution 710 in the final optimized sinogram estimate is more visible. However, the final estimate depicted by the exemplary distribution 710 in the estimated sinogram is still virtually indistinguishable from the corresponding distribution 715 in the true sinogram.

[0092] To quantify the difference between the exemplary distributions 705 and 710 in the estimated sinograms and the corresponding distribution 715 in the true sinogram, Fig. 9 A graph 900 is shown showing the difference between the exemplary distributions 705 and 710 in the estimated sinogram and the corresponding distribution 715 in the true sinogram. Specifically, graph 900 depicts a graph showing the difference 905 between the exemplary distribution 705 in the preliminary estimated sinogram and the corresponding distribution 715 in the true sinogram, and the difference 910 between the exemplary distribution 710 in the final estimated sinogram and the corresponding distribution 715 in the true sinogram. As shown, the error or difference 905 of the preliminary estimate is greater than the error or difference 910 of the final estimate.

[0093] The difference 910 is less than 0.01g / cm compared with the actual situation 2 , which corresponds to the maximum integral with respect to the water density (40 g / cm in this example 2 )0.025% error, and the difference 905 is less than 0.16g / cm compared with the actual situation. 2 , which is equivalent to a maximum integrated error of 0.4% relative to the water density. Therefore, although the preliminary estimate is relatively effective, if the material-based image is reconstructed directly from the preliminary estimate, some image artifacts may be caused by the difference 905, while using the optimized final estimate can effectively eliminate such image artifacts. In addition, because the preliminary estimate is substantially close to the true situation, using the preliminary estimate as an initial estimate for iteratively calculating the path length estimate provides a substantial improvement over using a predetermined initial estimate or an initial estimate derived from a model-based approach.

[0094] Although the method described above estimates multi-material path length directly, the method described above can be adapted to estimate monochromatic sinograms at multiple monochromatic energies (keV) instead of multi-material path length. For example, the method described above can be adapted to estimate monochromatic sinograms by simple variable substitution. As an illustrative example, Fig.10 A high-level flowchart is shown that illustrates an exemplary method 1000 for estimating a monochromatic sinusoidal pattern according to one embodiment. Specifically, the method 1000 involves estimating a monochromatic sinusoidal pattern at multiple energies based on calibration data to adaptively reduce noise covariance. Figure 1 and Figure 2 The method 1000 is described with reference to the systems and components of the present disclosure, but it should be understood that the method 1000 can be implemented with other systems and components without departing from the scope of the present disclosure. The method 1000 can be implemented, for example, as executable instructions in a non-transitory memory of the computing device 216 and / or the image reconstructor 230, and can be executed by a processor of the computing device 216 and / or the image reconstructor 230 to perform the actions described below.

[0095] The method 1000 starts at 1005. At 1005, the method 1000 performs a scan of a subject to acquire projection data. The scan includes a dual energy CT scan or a multi-energy CT scan of the subject. To this end, the method 1000 controls the x-ray controller 210 to drive the x-ray source 104 to emit x-rays at two or more energy levels, for example, while also controlling the gantry motor controller 212 and the table motor controller 226 to adjust the positions of the gantry 102 and the table 114, respectively, so that the position of the x-ray source 104 relative to the scanned subject 204 is adjusted while generating x-rays. The method 1000 also acquires projection data measured by an energy discriminating detector (including a photon counting detector, such as the detector element 202 of the detector array 108) via, for example, the DAS 214. The projection data includes multi-energy photon counting measurements, for example, where the photon measurements of different energies are classified into predefined energy intervals.

[0096] At 1010, method 1000 estimates a monochromatic sinogram at multiple keV based on the projection data and calibration data. For example, if the vector μ of the line integral of the monochromatic attenuation is represented as:

[0097] μ=[μ keV1 ,…,μ keVN ],

[0098] Then the changes in the variables in the calibration data can be applied. For example, the forward model f of the calibration data cal It can be expressed as:

[0099] f cal (μ) = f cal (Q×m),

[0100] Wherein, μ = Q × m represents the transformation from the material path length m to the monochromatic sinusoid. Therefore, whether the forward model f is used in the above method 300 or 400, cal (m), the independent variable can be replaced by the same amount, so that f cal (Q -1 ×μ). In this case, the unknown is the vector μ of the monochromatic sinusoid. For example, the iterative optimization can be reformulated so that the final estimated vector μ of the monochromatic sinusoid is final for:

[0101]

[0102] μ∈Ω μ , which can be as described above about Figure 3 Once the final estimated vector μ of the monochromatic sinusoidal graph is obtained final , the final estimate of the material path length m can be obtained by calculating the following formula final :

[0103] m final =Q -1 ×μ final .

[0104] Alternatively, the final estimate of the material path length m final This can be obtained by the following simpler optimization:

[0105]

[0106] m∈Ω m , where W μ is a matrix dependent on μ that appropriately weights the different components of μ for statistical benefit, and prior(m) comprises a function that applies prior information about the material path length m. It should be understood that since the material decomposition process has been performed by calculating μ as described above, final can be implemented in the image domain for m final The above expression for , thereby potentially removing the final estimate μ of the monochromatic sinusoidal pattern final For which method 1000 solves m in image space final As an example, the prior function (m) can include various image processing, machine learning, and deep learning-based penalty functions.

[0107] At 1015, method 1000 reconstructs a monochrome image from the estimated monochrome sinusoidal graph. For example, method 1000 utilizes the final estimated vector μ of the monochrome sinusoidal graph final The method 1000 may then perform image reconstruction to generate a corresponding monochrome image at multiple keV. Then, at 1020, the method 1000 outputs the monochrome image. For example, the method 300 may output the monochrome image to a display device, such as the display device 232. Additionally or alternatively, the method 300 may output the material density image to the mass storage device 218 for storage and / or the PACS 224 for remote viewing.

[0108] In addition, in some examples, at 1020, method 1000 optionally converts the monochrome image into a material density image. For example, method 1000 may convert the monochrome image into a material density image via a direct linear transformation in the image domain. In such examples, at 1025, method 1000 may optionally output the material density image. For example, method 1000 may output the material density image to a display device (such as display device 232), and / or a large capacity storage device 218 for storage and / or PACS 224 for remote viewing. Therefore, the systems and methods provided herein are capable of estimating a monochrome image from projection data and calibration data, and are capable of generating a material density image from a monochrome image, rather than obtaining a material density image and generating a monochrome image from the material density image. Then, method 1000 returns.

[0109] A technical effect of the present disclosure includes reconstructing material-based images of multiple materials. Another technical effect of the present disclosure includes reducing the computational complexity for calculating material path lengths. Another technical effect of the present disclosure includes improving the accuracy of material decomposition of two or more materials. Yet another technical effect of the present disclosure includes displaying material-based images of three or more materials generated from projection data without using physics modeling. Another technical effect of the present disclosure includes reconstructing two or more monochromatic images directly from projection data rather than material-based images.

[0110] In one embodiment, a method includes acquiring projection data of multiple x-ray spectra via an imaging system, estimating path lengths of multiple materials based on the projection data and calibration data of the imaging system, iteratively refining the estimated path lengths based on a linearized model derived from the calibration data, and reconstructing a material density image for each of the multiple materials from the iteratively refined estimated path lengths.

[0111] In a first example of the method, estimating the path lengths of the plurality of materials based on the projection data and the calibration data of the imaging system includes: estimating the path lengths of the plurality of materials based on the projection data and the calibration data without modeling the physics of the imaging system, the physics of the imaging system including a plurality of x-ray spectra and a spectral response of a detector of the imaging system. In a second example of the method optionally including the first example, estimating the path lengths of the plurality of materials based on the projection data and the calibration data of the imaging system includes: performing an inverse function lookup on the calibration data to generate a first estimate of the path lengths of the plurality of materials corresponding to the projection data. In a third example of the method optionally including one or more of the first and second examples, estimating the path lengths of the plurality of materials based on the projection data and the calibration data of the imaging system also includes: generating a linear approximation of a forward model constructed from the calibration data, and solving a system of linear equations based on the linear approximation to obtain a preliminary estimate of the path lengths of the plurality of materials. In a fourth example of the method, optionally including one or more of the first to third examples, iteratively refining the estimated path length based on the linearized model derived from the calibration data includes: iteratively calculating a final estimate of the path length for each of the plurality of materials, wherein the preliminary estimate of the path length is used as an initial estimate, wherein the material density image is reconstructed from the final estimate of the path length for each of the plurality of materials. In a fifth example of the method, optionally including one or more of the first to fourth examples, performing an inverse function lookup on the calibration data to generate a first estimate of the path length for the plurality of materials corresponding to the projection data includes: selecting a candidate vector of material path lengths in the calibration data for each sinogram interval, the candidate vector producing a result within a threshold distance of the projection data when input to the forward model, and for each sinogram interval, calculating the first estimate of the path length from a weighted sum of the candidate vectors. In a sixth example of the method, optionally including one or more of the first to fifth examples, the method further includes selecting weights for the weighted sum based on intensity measurements of the projection data, wherein higher intensities correspond to higher weights. In a seventh example of the method, optionally including one or more of the first to sixth examples, generating a linear approximation of the forward model includes: calculating a coefficient matrix that, when multiplied by a matrix of known path lengths for each material, results in a corresponding matrix of known projection measurements, wherein the calibration data includes the known path lengths and the known projection measurements. In an eighth example of the method, optionally including one or more of the first to seventh examples, solving a system of linear equations based on the linear approximation to obtain preliminary estimates of path lengths for a plurality of materials includes: calculating a matrix of preliminary estimates of path lengths for a plurality of materials that, when multiplied by the coefficient matrix, results in a corresponding matrix of projection measurements for the projection data.In a ninth example of the method optionally including one or more of the first to eighth examples, the method further includes: estimating at least two monochromatic sinusoids based on the projection data and the calibration data, and reconstructing at least two monochromatic images from the at least two monochromatic sinusoids. In a tenth example of the method optionally including one or more of the first to ninth examples, reconstructing a material density image for each of a plurality of materials includes: estimating a plurality of material density images from the at least two monochromatic images, the at least two monochromatic images being reconstructed from the at least two monochromatic sinusoids by using a pure image domain technique, the pure image domain technique including at least one of a linear transformation of the at least two monochromatic sinusoids and an iterative scheme using image domain prior information.

[0112] In another embodiment, a method includes: acquiring projection data of multiple x-ray spectra via an imaging system; calculating preliminary estimates of path lengths for multiple materials based on the projection data and calibration data for the imaging system without modeling the physics of the imaging system, the physics of the imaging system including the multiple x-ray spectra and the response of the detector to the multiple x-ray spectra; iteratively updating the preliminary estimates of path lengths for the multiple materials to obtain final estimates of path lengths for the multiple materials; and reconstructing a material density image for each of the multiple materials from the estimated path lengths.

[0113] In a first example of the method, calculating a preliminary estimate of the path length of the plurality of materials includes: calculating a first estimate of the path length based on an inverse function lookup of the calibration data, and solving a system of linear equations for the preliminary estimate of the path length of the plurality of materials, the system of linear equations being constructed based on a linear approximation of a forward model of the calibration data. In a second example of the method optionally including the first example, iteratively updating the preliminary estimate of the path length of the plurality of materials to obtain a final estimate of the path length of the plurality of materials includes: iteratively minimizing a statistical function with projection data and calibration data as inputs, initialized with the preliminary estimate of the path length of the plurality of materials, to determine a final estimate of the path length of the plurality of materials. In a third example of the method optionally including one or more of the first and second examples, the plurality of materials includes at least three materials.

[0114] In yet another embodiment, a system includes: an x-ray source configured to generate an x-ray beam toward a subject; a detector array including a plurality of detector elements configured to detect the x-ray beam attenuated by the subject; and a computing device communicatively coupled to the x-ray source and the detector array, the computing device being configured with instructions in a non-volatile memory that, when executed, cause the computing device to: control the x-ray source and the detector array to scan the subject with a plurality of x-ray beams at different energy levels and to acquire projection data; estimate path lengths of a plurality of materials based on the projection data and calibration data of the x-ray source and the detector array; iteratively refine the estimated path length based on a linearized model derived from the calibration data; and reconstruct a material density image of each of the plurality of materials from the iteratively refined estimated path length.

[0115] In a first example of the system, the computing device is further configured with instructions in a non-transitory memory that, when executed, cause the computing device to estimate path lengths of a plurality of materials based on the projection data and the calibration data without modeling the physics of the x-ray source and the detector array. In a second example of the system that optionally includes the first example, the computing device is further configured with instructions in a non-transitory memory that, when executed, cause the computing device to perform an inverse function lookup on the calibration data to generate a first estimate of the path lengths of the plurality of materials corresponding to the projection data, generate a linear approximation of a forward model constructed from the calibration data, and solve a system of linear equations based on the linear approximation to obtain a preliminary estimate of the path lengths of the plurality of materials. In a third example of the system that optionally includes one or more of the first and second examples, the computing device is further configured with instructions in a non-transitory memory that, when executed, cause the computing device to iteratively calculate a final estimate of the path length of each of the plurality of materials, wherein the preliminary estimate of the path length is used as an initial estimate, wherein the estimated path length used to reconstruct the material density image includes the final estimate of the path length. In a fourth example of the system optionally including one or more of the first to third examples, the computing device is further configured with instructions in a non-transitory memory, which when executed causes the computing device to estimate at least two monochromatic sinusoids based on the projection data and the calibration data and reconstruct at least two monochromatic images from the at least two monochromatic sinusoids. In a fifth example of the system optionally including one or more of the first to fourth examples, the system further includes a display device communicatively coupled to the computing device, and the computing device is further configured with instructions in a non-transitory memory, which when executed causes the computing device to output the material density image to the display device for display. In a sixth example of the system optionally including one or more of the first to fifth examples, the computing device is further configured with instructions in a non-transitory memory, which when executed causes the computing device to estimate a plurality of material density images from at least two monochromatic images, the at least two monochromatic images being reconstructed from the at least two monochromatic sinusoids by using a pure image domain technique, the pure image domain technique comprising at least one of a linear transformation of the at least two monochromatic sinusoids and an iterative scheme using image domain prior information.

[0116] As used herein, the elements or steps listed in the singular and beginning with the word "one" or "a kind of" should be understood as not excluding a plurality of the elements or steps, unless such exclusion is explicitly stated. In addition, the reference to "an embodiment" of the present invention is not intended to be interpreted as excluding the existence of additional embodiments that also include the cited features. In addition, unless explicitly stated to the contrary, the embodiment of "including", "comprising" or "having" an element or multiple elements with a specific characteristic may include additional such elements without the characteristic. The terms "including" and "in..." are used as the concise language equivalents of the corresponding terms "including" and "wherein". In addition, the terms "first", "second" and "third" etc. are only used as marks, and are not intended to impose numerical requirements or specific positional order on their objects.

[0117] This written description uses examples to disclose the invention, including the best mode, and also to enable a person of ordinary skill in the relevant art to practice the invention, including making and using any devices or systems and performing any included methods. The patentable scope of the invention is defined by the claims, and may include other examples that occur to a person of ordinary skill in the art. Such other examples are intended to fall within the scope of the claims if they have structural elements that do not differ from the literal language of the claims, or if they include equivalent structural elements with insignificant differences from the literal language of the claims.

Claims

1. A method for multi-material decomposition, comprising: acquiring projection data of a plurality of x-ray spectra via an imaging system; estimating path lengths for a plurality of materials based on the projection data and calibration data for the imaging system acquired with the imaging system during a calibration scan of a phantom including known materials and known path lengths, wherein the path lengths are used to specify line integrals of a material density image along lines connecting an x-ray source to respective detector elements; iteratively refining the estimated path length based on a linearized model derived from the calibration data; as well as A material density image for each of the plurality of materials is reconstructed from the iteratively refined estimated path lengths.

2. The method of claim 1 , wherein estimating the path lengths of the plurality of materials based on the projection data and the calibration data of the imaging system comprises: The path lengths of the plurality of materials are estimated based on the projection data and the calibration data without modeling the physics of the imaging system, including the plurality of x-ray spectra and a spectral response of a detector of the imaging system.

3. The method of claim 1 , wherein estimating the path lengths of the plurality of materials based on the projection data and the calibration data of the imaging system comprises: An inverse function lookup is performed on the calibration data to generate first estimates of path lengths for the plurality of materials corresponding to the projection data.

4. The method of claim 3, wherein estimating the path lengths of the plurality of materials based on the projection data and the calibration data of the imaging system further comprises: A linear approximation of a forward model constructed from the calibration data is generated, and a system of linear equations is solved based on the linear approximation to obtain preliminary estimates of the path lengths of the plurality of materials.

5. The method of claim 3, wherein performing the inverse function lookup on the calibration data to generate a first estimate of the path lengths of the plurality of materials corresponding to the projection data comprises: A candidate vector for the material path length in the calibration data is selected for each sinogram interval, the candidate vector producing a result within a threshold distance of the projection data when input into a forward model constructed from the calibration data, and a first estimate of the path length is calculated for each sinogram interval from a weighted sum of the candidate vectors.

6. The method of claim 4, wherein generating the linear approximation of the forward model comprises: A coefficient matrix is ​​calculated that, when multiplied by a matrix of known path lengths for each material, results in a corresponding matrix of known projection measurements, wherein the calibration data includes the known path lengths and the known projection measurements.

7. The method according to claim 1, further comprising: At least two monochromatic sinusoids are estimated based on the projection data and the calibration data, and at least two monochromatic images are reconstructed from the at least two monochromatic sinusoids.

8. A method for multi-material decomposition, comprising: acquiring projection data of a plurality of x-ray spectra via an imaging system; calculating, without modeling the physics of the imaging system, preliminary estimates of path lengths for the plurality of materials based on the projection data and calibration data for the imaging system acquired with the imaging system during a calibration scan of a phantom including known materials and known path lengths, wherein the path lengths are used to specify line integrals of a material density image along lines connecting an x-ray source to respective detector elements; iteratively updating preliminary estimates of the path lengths for the plurality of materials to obtain final estimates of the path lengths for the plurality of materials; as well as A material density image for each of the plurality of materials is reconstructed from the final estimates of the path lengths for the plurality of materials.

9. The method of claim 8, wherein calculating a preliminary estimate of the path lengths of the plurality of materials comprises: A first estimate of path length is calculated based on an inverse function lookup of the calibration data, and a system of linear equations constructed based on a linear approximation of a forward model of the calibration data is solved for the preliminary estimate of path length for the plurality of materials.

10. The method of claim 8, wherein iteratively updating preliminary estimates of the path lengths of the plurality of materials to obtain final estimates of the path lengths of the plurality of materials comprises: A statistical function is iteratively minimized with the projection data and the calibration data as inputs, initialized with preliminary estimates of the path lengths for the plurality of materials, to determine final estimates of the path lengths for the plurality of materials.

Images