阅读说明：本技术 用于多材料分解的方法和系统 (Method and system for multi-material decomposition ) 是由 萨斯什·拉马尼 吴明烨 B·德梅恩 彼得·伊迪克 于 2021-05-13 设计创作，主要内容包括：本发明题为“用于多材料分解的方法和系统”。本发明提供了用于计算机断层摄影术的多材料分解的各种方法和系统。在一个实施方案中,一种方法包括：经由成像系统采集多个x射线光谱的投影数据,基于投影数据和成像系统的校准数据来估计多种材料的路径长度,基于从校准数据导出的线性化模型迭代地细化估计路径长度,以及由迭代细化的估计路径长度重建多种材料中的每种材料的材料密度图像。通过在不对成像系统的物理学建模的情况下以这种方式确定路径长度估计,可更快速地执行准确的材料分解并且对系统的物理学变化的敏感性较低,并且此外可扩展到两种以上的材料。(The invention provides a method and system for multi-material decomposition. Various methods and systems for multi-material decomposition for computed tomography are provided. In one embodiment, a method comprises: acquiring projection data for 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, iteratively refining the estimated path lengths based on a linearization model derived from the calibration data, and reconstructing a material density image for each material of the plurality of materials from the iteratively refined estimated path lengths. By determining 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 physical changes of the system, and furthermore can be extended to more than two materials.)

1. A method, comprising:

acquiring projection data for a plurality of x-ray spectra via an imaging system;

estimating path lengths of a plurality of materials based on the projection data and calibration data of the imaging system;

iteratively refining the estimated path length based on a linearized model derived from the calibration data; and

reconstructing a material density image for each material of the plurality of materials from the iteratively refined estimated path lengths.

2. The method of claim 1, wherein estimating path lengths of the plurality of materials based on the projection data and the calibration data of the imaging system comprises: estimating the path lengths of the plurality of materials based on the projection data and the calibration data without modeling physics of the imaging system, 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 path lengths of the plurality of materials based on the projection data and the calibration data of the imaging system comprises: performing an inverse function lookup of the calibration data to generate a first estimate of a path length of the plurality of materials corresponding to the projection data.

4. The method of claim 3, wherein estimating path lengths of the plurality of materials based on the projection data and the calibration data of the imaging system further comprises: 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 preliminary estimates of path lengths of the plurality of materials.

5. The method of claim 3, wherein performing the inverse function lookup of the calibration data to generate a first estimate of the path length of the plurality of materials corresponding to the projection data comprises: selecting, for each sinogram interval, a candidate vector of material path lengths in the calibration data that, when input to a forward model constructed from the calibration data, produces results within a threshold distance of the projection data, and calculating, for each sinogram interval, a first estimate of the path length from a weighted sum of the candidate vectors.

6. The method of claim 4, wherein generating the linear approximation of the forward model comprises: calculating a matrix of coefficients 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 comprises the known path lengths and the known projection measurements.

7. The method of claim 1, further comprising: at least two monochromatic sinograms 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 sinograms.

8. A method, comprising:

acquiring projection data for a plurality of x-ray spectra via an imaging system;

calculating preliminary estimates of path lengths of a plurality of materials based on the projection data and calibration data of the imaging system without modeling physics of the imaging system including the plurality of x-ray spectra and a response of a detector to the plurality of x-ray spectra;

iteratively updating the preliminary estimates of the path lengths of the plurality of materials to obtain final estimates of path lengths of the plurality of materials; and

reconstructing a material density image for each material of the plurality of materials from the final estimate of the path length 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: calculating a first estimate of path length based on an inverse function lookup of the calibration data, and solving a system of linear equations for the preliminary estimates of path length for the plurality of materials, the system of linear equations being constructed based on a linear approximation of a forward model of the calibration data.

10. The method of claim 8, wherein iteratively updating the 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: iteratively minimizing a statistical function, wherein the projection data and the calibration data are input, initialized with preliminary estimates of the path lengths of the plurality of materials to determine a final estimate of the path lengths of the plurality of materials.

Technical Field

Embodiments of the subject matter disclosed herein relate to Computed Tomography (CT) imaging, and more particularly to multi-material decomposition for CT imaging.

Background

A dual-or multi-energy spectrum Computed Tomography (CT) imaging system may display the density of different materials in an object and generate images corresponding to a plurality of monochromatic x-ray energy levels. CT imaging systems may derive behavior at different monochromatic energy levels based on signals from at least two regions of photon energy in the spectrum (e.g., a low energy portion and a 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 x-ray attenuation process is dominated by two physical processes: compton scattering and the photoelectric effect. The signals detected from the two energy regions provide sufficient information to resolve the energy dependence on the material being imaged. The signals detected from the two energy regions provide sufficient information to determine the relative composition of an object composed of two hypothetical materials.

Disclosure of Invention

In one embodiment, a method comprises: acquiring projection data for a plurality of x-ray spectra via an imaging system; estimating path lengths of a plurality of materials based on the acquired projection data and calibration data of the imaging system; and reconstructing a material density image for each of the plurality of materials from the estimated path lengths. In the following, the term "path length" is used to specify the line integral of the material density image along the line connecting the x-ray source to the individual detector elements. By determining 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 physical changes of the system, and furthermore can be extended to more than two materials.

It should be appreciated that the brief description above is provided to introduce in simplified form a selection of 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 defined uniquely 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.

Drawings

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

fig. 1 shows a pictorial view of an imaging system according to an embodiment;

FIG. 2 shows a block schematic diagram of an exemplary imaging system according to one embodiment;

FIG. 3 depicts a high level flow chart illustrating an exemplary method for multi-material decomposition in accordance with one embodiment;

FIG. 4 depicts a high level flow chart that illustrates an exemplary method for calculating preliminary estimates of path lengths for multiple materials based on calibration data, in accordance with one embodiment;

FIG. 5 shows a set of graphs depicting an exemplary multi-dimensional surface for initial estimation of path lengths of a plurality of materials, in accordance with one embodiment;

FIG. 6 shows a set of graphs depicting an exemplary linear model for refining an initial estimate of path length for a plurality of materials, in accordance with one embodiment;

FIG. 7 shows a graph illustrating an exemplary estimated sinogram and a true sinogram of a material according to one embodiment;

FIG. 8 shows a graph illustrating an enlarged view of the graph of FIG. 7;

figure 9 illustrates a graph showing the difference of the estimated sinogram of figure 7 from the true sinogram, in accordance with one embodiment; and is

FIG. 10 illustrates a high level flow chart showing an exemplary method for estimating a monochromatic sinogram, according to one embodiment.

Detailed Description

The following description relates to various embodiments of spectral Computed Tomography (CT) imaging. In particular, methods and systems for multi-material decomposition in spectral CT imaging are provided. An example of a CT imaging system that may be used to acquire images in accordance with the present techniques is shown in FIGS. 1 and 2. CT imaging systems may be configured with energy discriminating detectors, such as photon counting detectors, that provide fidelity in which materials can be distinguished via material decomposition. The material decomposition may be performed in the projection domain, the image domain, or in conjunction with the reconstruction. However, the joint material decomposition-reconstruction method is generally computationally expensive compared to both the projection domain method and the image domain method. Furthermore, projection domain methods for multi-material decomposition (i.e., material decomposition of two or more materials) typically rely on a priori 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 practice, such a priori knowledge is difficult to obtain with high accuracy, especially if the detector spectral response varies across the detector array and/or when the emitted x-ray spectrum varies over the x-ray tube lifetime. Methods for multi-material decomposition, such as the method depicted in fig. 3, overcome the challenges of these previous methods by estimating multi-material path lengths based solely on calibration data and without knowledge of the physics of the multi-energy CT imaging system. The method comprises a two-step approach, wherein a preliminary estimate is first obtained based on calibration data and acquired projection data, and then an iterative scheme is initialized using the preliminary estimate, which optimizes statistical criteria to further enhance the statistical accuracy of the multi-material path length. As shown in fig. 4, the method for obtaining a preliminary estimate of the multi-material path length includes a combination of a function back-finding and a local multi-linear fitting method. The process of calculating such preliminary estimates is illustrated in fig. 5 and 6 with simulated data. The results of further simulations are shown in fig. 7 to 9, indicating that the preliminary estimate is highly accurate, so that the iterative scheme can obtain a final estimate of the multi-material path length with even higher accuracy with only a few iterations. The method described herein for directly estimating multi-material path lengths may also be adapted to estimate monochromatic sinograms at multiple energies, as shown in FIG. 10.

Fig. 1 shows an exemplary CT system 100 configured for CT imaging. Specifically, the CT system 100 is configured to image a subject 112 (such as a patient, an inanimate object, one or more manufacturing components) and/or a foreign object (such as a dental implant, a stent, and/or a contrast agent present within 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 a beam of x-ray radiation 106 (see fig. 2) for imaging a subject 112 lying on a table 114. In particular, the x-ray source 104 is configured to project a beam of x-ray radiation 106 toward a detector array 108 positioned on an opposite side of the gantry 102. Although fig. 1 depicts only one x-ray source 104, 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 a patient. In some embodiments, the x-ray source 104 may enable dual-spectrum imaging by fast peak kilovoltage (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, where 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. Thus, it should be understood that the methods described herein may be implemented with a variety of multispectral acquisition techniques and are not limited to a particular described embodiment.

In certain embodiments, the CT system 100 further comprises 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), and so forth. As further described herein, in some examples, the image processor unit 110 may use an analytical image reconstruction method (such as FBP) in addition to the iterative image reconstruction method.

In some CT imaging system configurations, an X-ray source projects a cone-shaped beam of X-ray radiation that is collimated to lie within an X-Y-Z plane of a Cartesian coordinate system and is commonly referred to as an "imaging plane". The x-ray radiation beam passes through an object being imaged, such as a patient or subject. The x-ray radiation beam impinges upon an array of detector elements after being attenuated by the object. 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 measure of the x-ray beam attenuation at the detector location. Attenuation measurements from all detector elements are acquired separately to produce a transmission profile.

In some CT systems, a gantry is used to rotate an x-ray source and a detector array around an object to be imaged in an imaging plane such that the angle at which the radiation beam intersects the object constantly changes. 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 the object includes a set of views made at different gantry angles, or view angles, during one rotation of the x-ray source and detector. It is contemplated that the benefits of the methods described herein stem from medical imaging modalities other than CT, and thus, as used herein, the term "view" is not limited to the use described above with respect to projection data from one gantry angle. The term "view" is used to refer to one data acquisition whenever there are multiple data acquisitions from different angles, whether from a CT imaging system or another imaging system acquiring multispectral attenuation measurements.

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, an image corresponding to a three-dimensional representation of the object. One method for reconstructing an image from a set of projection data is known in the art as the filtered backprojection 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, as well as iterative reconstruction techniques. The method converts the attenuation measurements from the scan into "CT numbers" or "Hounsfield units" which are used to control the brightness of corresponding pixels on the display device.

To reduce the total scan time, a "helical" scan may be performed. To perform a "helical" scan, the patient is moved with the table 114 while data for a prescribed number of slices is acquired. Such systems produce a single helix from a cone beam helical scan. The helix mapped out by the cone beam yields projection data from which an image in each prescribed slice can be reconstructed.

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

Fig. 2 shows an exemplary imaging system 200 similar to CT system 100 of fig. 1. According to aspects of the present disclosure, the imaging system 200 is configured for imaging a subject 204 (e.g., the subject 112 of fig. 1). In one embodiment, the imaging system 200 includes a detector array 108 (see fig. 1). The detector array 108 also includes a plurality of detector elements 202 that together sense the x-ray radiation beam 106 (see fig. 2) passing through a subject 204, such as a patient, to acquire corresponding projection data. Thus, in one embodiment, detector array 108 is fabricated 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 acquiring projection data.

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

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 is then subject to pre-processing and calibration to adjust the data to represent the line integrals of the attenuation coefficients of the scanned subject 204. The processed data is commonly referred to as projections.

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

The acquired projection data set may be used for Basis 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-each corresponding map or image of the 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 underlying material (e.g., bone, soft tissue, and/or contrast agent) in the imaging volume.

Once reconstructed, the base material image produced by the imaging system 200 shows the internal features of the subject 204 in terms of the densities of the two base materials. A density image may be displayed to show these features. In conventional methods of diagnosing medical conditions (such as disease states), and more generally medical events, a radiologist or physician will consider a hard copy or display of density images to discern characteristic features of interest. Such features may include lesions, sizes and shapes of particular anatomical structures or organs, as well as other features that should be discernable in the image based on the skill and knowledge of the individual practitioner.

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 further includes an x-ray controller 210 configured to provide power and timing signals to the x-ray source 104. In addition, 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.

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

In addition, the computing device 216 provides commands and parameters to one or more of the DAS 214, x-ray controller 210, and gantry motor controller 212 to control system operations, such as data acquisition and/or processing. In certain embodiments, the computing device 216 controls system operation based on operator input. The computing device 216 receives operator input, including, for example, commands and/or scanning parameters, via an operator console 220 operatively 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 scanning parameters.

Although only one operator console 220 is shown in fig. 2, more than one operator console may be coupled to the imaging system 200, for example, for inputting or outputting system parameters, requesting examinations, drawing data, and/or viewing images. Further, in certain embodiments, the imaging system 200 may be coupled to a plurality of displays, printers, workstations, and/or the like, located locally or remotely, e.g., within an institution or hospital, or at disparate locations, via one or more configurable wired and/or wireless networks (such as the internet and/or virtual private networks, wireless telephone networks, wireless local area networks, wired local area networks, wireless wide area networks, wired wide area networks, and the like).

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.

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

As previously described, DAS 214 samples and digitizes projection data acquired by detector elements 202. Subsequently, the image reconstructor 230 performs a high speed reconstruction using the sampled and digitized x-ray data. Although fig. 2 illustrates the image reconstructor 230 as a separate entity, in certain embodiments, the image reconstructor 230 may form a portion of the computing device 216. Alternatively, the image reconstructor 230 may not be present in the imaging system 200, and the computing device 216 may instead perform one or more functions of the image reconstructor 230. Further, the image reconstructor 230 may be located locally or remotely and may be operatively 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.

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 available patient information for diagnosis and evaluation. In certain embodiments, the computing device 216 may transmit the reconstructed image and/or 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 images 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.

Various methods and processes described further herein, such as the methods described below with reference to fig. 3, 4, and 10, may be stored as executable instructions in 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 a non-transitory memory and may apply the methods described herein to reconstruct an image from scan data. In another embodiment, the computing device 216 may include instructions in 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 across the image reconstructor 230 and the computing device 216.

In one embodiment, display 232 allows the operator to assess the imaged anatomy. The display 232 may also allow an operator to select a volume of interest (VOI) and/or request patient information, e.g., via a Graphical User Interface (GUI), for subsequent scanning or processing.

Fig. 3 illustrates a high-level flow diagram showing an exemplary method 300 for multi-material decomposition, according to one embodiment. In particular, the method 300 involves performing multi-material decomposition in the projection domain based on calibration data rather than based on physics knowledge of the imaging system (such as models of incident x-ray spectra and/or detector responses). The method 300 is described with reference to the systems and components of fig. 1 and 2, but it should be understood that the method 300 may be implemented with other systems and components without departing from the scope of the present disclosure. The method 300 may be implemented, for example, as executable instructions in a non-transitory memory of the computing device 216 and/or the image reconstructor 230, and may be executed by a processor of the computing device 216 and/or the image reconstructor 230 to perform the acts described herein below.

The method 300 begins at 305. At 305, the method 300 performs a scan of the 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, such that the position of the x-ray source 104 relative to the scanned subject 204 is adjusted while generating x-rays. Method 300 also acquires projection data measured by energy discriminating detectors (including photon counting detectors such as detector elements 202 of detector array 108), e.g., via DAS 214. The projection data includes multi-energy photon count measurements, e.g., where photon measurements of different energies are classified into predefined energy intervals.

At 310, the method 300 calculates preliminary estimates of path lengths for the plurality of materials based on the calibration data and the projection data. The calibration data comprises calibration data for calibrating the imaging system, which is acquired with the imaging system during a calibration scan of a phantom comprising known materials and a known path length. The method 300 calculates a preliminary estimate of the path length of the plurality of materials based on the calibration data, which will result in observed multi-energy measurements of the projection data. To this end, the 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 observed multi-energy measurements, and refine the initial estimate by solving a system of linear equations constructed based on the initial estimate. An exemplary method for calculating preliminary estimates of path lengths for a plurality of materials based on calibration data is further described herein with reference to fig. 4.

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 as the initial estimate. For example, the method 300 may calculate the final estimate m by iteratively solving the following equation^final：

m∈Ω_mWherein Ω is_mComprising a set of (physically) feasible material path length vectors, F is a statistical criterion (e.g. log-likelihood) that, when optimized, yields a (physically) meaningful estimate of the multi-material path length, F^calIncluding a distinguishable forward model that maps known material path lengths to corresponding known projection values (p-values) or x-ray intensities (I-values), p-values, using interpolation or polynomial models^measIncluding the p value measured at 305, and I^measIncluding the x-ray intensity measured at 305. Although the above formula applies to the entire sinogram (i.e., all sinogram intervals i), it should be understood that the above is used to finally estimate m according to the component of F^finalThe equations of (a) can be separated in the sinogram interval i and therefore solved in parallel.

The above statistical problem may be understood as setting a maximum likelihood (or maximum a posteriori) estimator for a desired multi-material path length. In one example, the function F may comprise a poisson log-likelihood function such that:

where index i indexes sinogram intervals, index k indexes energy intervals, measured x-ray intensitiesCorresponding to the scanning of air in the absence of material,to map known material path lengths to known normalized x-ray intensitiesWherein the distinguishable forward model isAndcorresponding to the corrected intensity measurement and the air scan, respectively, in the absence of material. In another example, the function F may comprise a weighted least squares criterion function such that:

whereinAndin proportion, andis a distinguishable 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 the function measures the x-ray intensity I^measLower (or even zero in the case of photon counting detectors, or negative due to electronic noise in the case of energy integrating detectors) scenarios, which prevent the calculation of p^measApplication of (2) log (-), wherein p^measComponent (b) ofCorresponding to the sinogram interval index i and the energy interval index k, is computed as:

at 320, the method 300 reconstructs a material density image for each material based on the final estimate of path length. For example, the final estimate of the path length for each material includes a material-based projection, so the method 300 reconstructs a material density image for each material from the final path length estimate for the material. At 325, the method 300 outputs a material density image. For example, method 300 may output the material density image to a display device, such as 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.

At 330, method 300 generates a monochromatic image at multiple energies based on the material density image. For example, the method 300 may selectively combine material density images to generate a monochromatic image at a given energy that simulates acquisition using an x-ray source that emits photons only at the given energy. The method 300 may generate a plurality of monochromatic images from the material density image at different energies. At 335, method 300 outputs a monochrome image. For example, method 300 may output the monochrome image to a display device, such as 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. The method 300 then returns.

Fig. 4 illustrates a high-level flow diagram showing an exemplary method 400 for calculating preliminary estimates of path lengths for multiple materials based on calibration data, according to one embodiment. In particular, 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 the forward model at the first estimate. The method 400 is described with reference to the systems and components of fig. 1 and 2, but it should be understood that the method 400 may be implemented with other systems and components without departing from the scope of the present disclosure. The method 400 may be implemented, for example, as executable instructions in a non-transitory memory of the computing device 216 and/or the image reconstructor 230, and may be executed by a processor of the computing device 216 and/or the image reconstructor 230 to perform the acts described herein below. Method 400 may include a subroutine of method 300, and in particular may include an act 310 of calculating preliminary estimates of path lengths of the plurality of materials based on the calibration data and the projection data.

The method 400 begins at 405. At 405, the method 400 loads projection data. For example, the method 400 loads the projection data acquired at 305 as described above. The projection data includes multi-energy photon count measurements, e.g., where photon measurements of different energies are classified into predefined energy intervals.

Continuing at 410, the 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 to calibrate the imaging system using a phantom comprising known materials and known path lengths. For example, the calibration data may be provided as a set S^calStoring in a non-transitory memory:

wherein each vector isLx1 vectors of known path length of L material, each vectorThe Kx1 vector, i is the sinogram interval index, and j is the experiment index, for known p-values (log-normalized intensity values) corresponding to K-1, …, K-energy intervals, where different experiments correspond to different material combinations. Distinguishable forward model of the ith sinogram intervalIs based on a calibration set S^calConstructed by mapping the material path length m to a p-value, such that using interpolation values such as non-uniform rational B-splines or lagrange interpolation values,

or, alternatively, using a polynomial model,

then, the measurement vector is logarithmically normalized for a given K energy interval,

and a corresponding x-ray intensity measurement vector,

corresponding to a sinogram interval i (i.e., corresponding to one x-ray projection line), the method 400 derives calibration data S^calTo find possible candidate vectorsThe candidate vector generates a proximity measurement vectorP-value vector ofFor example, the method 400 may determine:

wherein the function 1_(.)Is an indicator function that is one (i.e., one), 3 ≦ K if the condition in brackets is satisfied or is otherwise zero_sub≤K,T_pIs thatToAnd 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 comprise an absolute difference or a squared difference function. The result of the resulting set of search results,

thus, including the set of material vectors, the method 400 obtains a rough estimate of the desired multi-material path length from the set of material vectors. For example, the method 400 calculates an initial estimate by performing a simple weighted sum of material vectors:

wherein w_nIncluding imparting a material path lengthHigher priority weights, their corresponding p-value vectorsIs closer in value toThe numerical proximity depends on the measured vectorP value of (1)Is then based on the intensity measurement vectorCorresponding x-ray intensity measurements ofThe higher the intensity, the higher the confidence. Thus, the weighted model provides a statistically relevant estimate by giving higher weights to energy intervals with higher detected x-ray intensities.

At 415, the method 400 generates a linear approximation of the forward model based on the first estimate of the path length. For example, the method 400 may be implemented by evaluating the first estimatePerforming a forward model on surrounding calibration dataFurther refining the first estimate by local multiple linear approximation modelingThe multiple linear approximation model produces a system of linear equations that can be expressed as:

wherein the coefficient matrixComprising coefficients of a linear model passing through a pair matrixIn a known group of layersSolving the linear equation set to obtain a matrix corresponding to the matrixWhere the known path lengths of the multi-linear basis vectors are stacked (i.e.,one path length for each material. Note that in the following paragraphs, the vector m_iandp_iNow the linearization around the operating point defined by the search step is shown. After determining the coefficient matrix, the method 400 then passes the measured p-value based on a first estimate of the path lengthSets a similar linear system to generate a linear approximation of the forward model such that:

continuing at 420, the method 400 solves a system of linear equations based on linear approximation to obtain a preliminary estimate of path length. For example, to solve the above as a preliminary estimate of path lengthThe method 400 calculates, for the established linear system:

wherein the weighting matrixAnd measured intensityProportionally to give relatively better weighting to energy intervals with higher measured x-ray intensities.

An exemplary example of a local linear model of the forward model constructed based on the first estimate of path length is further described herein with reference to fig. 6.

At 425, the method 400 outputs a preliminary estimate of the path length. For example, the method 400 may output a preliminary estimate of the path length to memory, such that the method 300 may iteratively calculate a final estimate of the path length for each material, with the preliminary estimate as the initial estimate, as described above. The method 400 then returns.

Fig. 5 shows a set of graphs 500 depicting an exemplary multi-dimensional surface for initial estimation of path lengths of multiple materials, in accordance with one embodiment. The set of graphs 500 includes a graph for each of a plurality of energy intervals, in particular, in the depicted example K-8 spectral energy intervals, 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 mapping of a material path length (e.g.,) Mapping to log-normalized calibration data (e.g.,) The forward direction model of (e.g.,) A graph of the surface of (a).

For example, the first graph 510 shows a first surface 511 of a forward model that maps material path lengths to log-normalized calibration data p for a first energy interval (e.g., k 1) of a photon-counting CT imaging system_ij1. Specifically, first surface 511 maps the material path lengths of the two materials, water (W) and iodine (I), to calibration data. Black dots on each surface including the first surface 511 correspond to measurements along the x-ray pathIncluding the measured values 513 (e.g.,)。[W，I]points 516 and 517 on the plane represent [ W, I ]]Possible estimation of path lengthThese estimates use, for example, S^lookupCriteria are collectively derived from graphs 510, 520, 530, 540, 550, 560, 570, and 580, and from which corresponding measurements are derivedOf material along the x-ray path [ W, I ]]A first estimate 515 of the path length (e.g.,) (i.e., an 8 x1 vector of p values in the depicted example), as described above. Point 519 indicates the actual or real path length of the material.

Preliminary estimate 515 at [ W, I]In a plane, and the representations correspond to measured valuesThe estimated path lengths of water and iodine. Thus, it is possible to provideFor illustrative purposes, the preliminary estimate 515 is replicated in all graphs 520, 530, 540, 550, 560, 570, and 580. It should be understood that two materials, water and iodine, are depicted for ease of illustration, as the dimensions of a forward model of three or more materials are difficult to show in a two-dimensional representation. However, the illustrated method is applicable to two or more materials.

To improve the first estimate obtained from the surface of the forward model shown in fig. 5, a local joint linear model of the surface may be constructed. As an illustrative and non-limiting example, fig. 6 shows a set of graphs 600 depicting an exemplary linear model for refining an initial estimate of path length for a plurality of materials depicted in the set of graphs 500 of fig. 5.

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

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

Alternatively, the first estimate (e.g.,) And to improve the estimation (e.g.,) Both can be modeled using the intensity domain forward directionObtained by appropriately modifying the method 400 described above to use the intensity domain calibration data accordinglyAnd measurement dataMapping known material path lengths to known x-ray intensity valuesIt may be advantageous to operate in the intensity domain, especially ifHas a poor signal-to-noise ratio, or ifIs non-positive.

The preliminary estimate thus obtained is used to initialize an iterative optimizer (e.g., as described above with reference to fig. 3) to further refine the path length estimate. In some examples, the preliminary estimate of path length thus obtained may be used directly to reconstruct a material-based image, since it is substantially close to the real case. That is, in some examples, iterative optimization of path length estimates initialized with preliminary estimates may be omitted from the multi-energy material decomposition. Furthermore, in examples where iterative optimization is performed, the iterative method converges on the final estimate of the path length with fewer iterations and greater accuracy. Since the step of performing an inverse function lookup to compute the first estimate and generating a linear approximation of the forward model based on the first estimate is not a computational burden, the overall computational complexity of the material decomposition for multiple energies is reduced. Further, according to the systems and methods provided herein, accurate material decomposition can be performed without the need to model or model the physics of the imaging system or the multiple energies and x-ray spectra of multiple materials.

To illustrate the efficacy of the techniques provided herein for multi-energy material decomposition, fig. 7 shows a graph 700 illustrating exemplary distributions 705 and 710 and a corresponding true distribution 715 in an estimated sinogram of a material according to one embodiment. In particular, the graph 700 depicts an exemplary distribution 705 in the estimated sinogram corresponding to a preliminary estimate (e.g., computed according to the method 400 as described above) and an exemplary distribution 710 in the estimated sinogram corresponding to an optimized estimate (e.g., iteratively computed at 315 as described above).

As shown, exemplary distributions 705, 710, and 715 in the estimated sinogram are similar enough to be indistinguishable in fig. 7. To illustrate the differences between exemplary distributions 705, 710, and 715 in the estimated sinogram, fig. 8 shows a graph 800 showing an enlarged view of region 725 of graph 700. In this view, the deviation of the preliminary estimate, depicted by 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 effectively indistinguishable from the corresponding distribution 715 in the true sinogram.

To quantify the differences between the exemplary distributions 705 and 710 in the estimated sinogram and the corresponding distribution 715 in the true sinogram, fig. 9 shows a graph 900 showing the differences of the exemplary distributions 705 and 710 in the estimated sinogram and the corresponding distribution 715 in the true sinogram. In particular, the graph 900 depicts a graph showing the difference 905 of an exemplary distribution 705 in the preliminary estimate sinogram from the corresponding distribution 715 in the true sinogram, and the difference 910 of an exemplary distribution 710 in the final estimate sinogram from 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.

The difference 910 is less than 0.01g/cm compared to the real case²This corresponds to the maximum integral with respect to the density of water (40 g/cm in this example)²) 0.025% error, and the difference 905 is lower than 0.16g/cm compared to the real case²This corresponds to an error of 0.4% with respect to the maximum integral of the density of water. Thus, while 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, and such image artifacts may be effectively eliminated using the optimized final estimate. Furthermore, since the preliminary estimate is substantially close to the real case, 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.

While the method described above directly estimates multi-material path length, the above method can be adjusted to estimate monochromatic sinograms at multiple monochromatic energies (keV) instead of multi-material path length. For example, the method described above may be adapted to estimate the monochromatic sinogram by simple variable substitution. As an illustrative example, fig. 10 shows a high-level flow diagram illustrating an illustrative method 1000 for estimating a monochromatic sinogram, according to one embodiment. In particular, the method 1000 involves estimating a monochromatic sinogram at multiple energies based on calibration data to adaptively reduce noise covariance. The method 1000 is described with reference to the systems and components of fig. 1 and 2, but it should be understood that the method 1000 may be implemented with other systems and components without departing from the scope of the present disclosure. The method 1000 may be implemented, for example, as executable instructions in a non-transitory memory of the computing device 216 and/or the image reconstructor 230, and may be executed by a processor of the computing device 216 and/or the image reconstructor 230 to perform the acts described herein below.

Method 1000 begins at 1005. At 1005, method 1000 performs a scan of the 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, such that the position of the x-ray source 104 relative to the scanned subject 204 is adjusted while generating x-rays. Method 1000 also acquires projection data measured by energy discriminating detectors (including photon counting detectors such as detector elements 202 of detector array 108) via, for example, DAS 214. The projection data includes multi-energy photon count measurements, e.g., where photon measurements of different energies are classified into predefined energy intervals.

At 1010, the method 1000 estimates monochromatic sinograms at a plurality of keV based on the projection data and the calibration data. For example, if the vector μ of the line integrals of the monochromatic attenuations is represented as:

μ＝[μ_keV1,…,μ_keVN]，

the variation in the variables in the calibration data can then be applied. For example, a forward model f of the calibration data^calCan be expressed as:

f^cal(μ)＝f^cal(Q×m)，

where μ ═ Q × m denotes the transformation from the material path length m to the monochromatic sinogram. Thus, the forward model f is used in either of the foregoing with respect to method 300 or 400^cal(m) the independent variables may all be substituted by equal amounts, so that f^cal(Q^-1X μ). Unknown in this case is the vector μ of the monochromatic sinogram. For example, it can be re-expressedIterative optimization to make the final estimation vector mu of the monochromatic sinogram^finalComprises the following steps:

μ∈Ω_μwhich may be solved as described above with respect to fig. 3. Once the final estimate vector μ of the monochromatic sinogram is obtained^finalThe final estimate m of the material path length can be obtained by calculating the following equation^final：

m^final＝Q^-1×μ^final。

Alternatively, the final estimate m of the material path length^finalCan be obtained by the following simpler optimization:

m∈Ω_mwherein W is_μIs a matrix dependent on μ that weights the different components of μ appropriately for statistical benefit, and the 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 calculated as μ as described above^finalCan thus be implemented in the image domain for m^finalThereby potentially removing the final estimate μ of the monochromatic sinogram^finalAny beam hardening errors in (2). For a system in which method 1000 solves for m in image space^finalThe prior function (m) may include a variety of image processing, machine learning, and penalty functions based on deep learning.

At 1015, method 1000 reconstructs a monochromatic image from the estimated monochromatic sinogram. For example, the method 1000 utilizes the final estimate vector μ of the monochromatic sinogram^finalImage reconstruction is performed to generate corresponding monochromatic images at a plurality of keV. Then, at 1020, method 1000 outputs a monochrome image. For example, the method 300 may output the monochrome image to a display device, such as a display device232. 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.

Further, in some examples, at 1020, method 1000 optionally converts the monochrome image to a material density image. For example, method 1000 may convert a monochrome image to a material density image via a direct linear transformation in the image domain. In such examples, at 1025, method 1000 may optionally output a material density image. For example, the method 1000 may output the material density image to a display device (such as display device 232), and/or mass storage 218 for storage, and/or PACS 224 for remote viewing. Thus, rather than obtaining and generating a monochromatic image from a material density image, the systems and methods provided herein can estimate a monochromatic image from projection data and calibration data, and can generate a material density image from the monochromatic image. Method 1000 then returns.

Technical effects of the present disclosure include reconstructing material-based images of a plurality of materials. Another technical effect of the present disclosure includes reducing 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 a material-based image.

In one embodiment, a method comprises: acquiring projection data for 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, iteratively refining the estimated path lengths based on a linearization model derived from the calibration data, and reconstructing a material density image for each material of the plurality of materials from the iteratively refined estimated path lengths.

In a first example of the method, estimating path lengths of the plurality of materials based on the projection data and 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, which includes the plurality of x-ray spectra and the spectral response of the detectors of the imaging system. In a second example of the method, which optionally includes the first example, estimating path lengths of the plurality of materials based on the projection data and calibration data of the imaging system includes: an inverse function lookup of the calibration data is performed to generate a first estimate of path lengths of the plurality of materials corresponding to the projection data. In a third example of the method, which optionally includes one or more of the first example and the second example, estimating path lengths of the plurality of materials based on the projection data and calibration data of the imaging system further comprises: 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 preliminary estimates of path lengths for the plurality of materials. In a fourth example of the method, which optionally includes one or more of the first through third examples, iteratively refining the estimated path length based on a linearized model derived from the calibration data comprises: 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 the 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, which optionally includes one or more of the first example through the fourth example, performing an inverse function lookup of the calibration data to generate a first estimate of path length of the plurality of materials corresponding to the projection data comprises: a candidate vector of material path lengths in the calibration data is selected for each sinogram interval that, when input to the forward model, produces results within a threshold distance of the projection data, and a first estimate of path length is calculated from a weighted sum of the candidate vectors for each sinogram interval. In a sixth example of the method, which optionally includes one or more of the first through fifth examples, the method further comprises selecting a weight for the weighted sum based on intensity measurements of the projection data, wherein a higher intensity corresponds to a higher weight. In a seventh example of the method, which optionally includes one or more of the first through sixth examples, generating the linear approximation of the forward model comprises: a matrix of coefficients is calculated which, when multiplied by a matrix of known path lengths for each material, yields 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, which optionally includes one or more of the first through seventh examples, solving a system of linear equations based on linear approximation to obtain preliminary estimates of path lengths of a plurality of materials comprises: a matrix of preliminary estimates of path lengths for the plurality of materials is calculated, which when multiplied by the coefficient matrix yields 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 through eighth examples, the method further comprises: at least two monochromatic sinograms 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 sinograms. In a tenth example of the method, which optionally includes one or more of the first through ninth examples, reconstructing the material density image for each of the plurality of materials comprises: a plurality of material density images are estimated from at least two monochromatic images reconstructed from at least two monochromatic sinograms using a pure image domain technique including at least one of a linear transformation of the at least two monochromatic sinograms and an iterative approach using image domain prior information.

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

In a first example of the method, calculating a preliminary estimate of the path lengths of the plurality of materials comprises: the method further 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 estimates of the path lengths for 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, which optionally includes the first example, iteratively updating the preliminary estimates of path lengths of the plurality of materials to obtain a final estimate of path lengths of the plurality of materials includes: the statistical function is iteratively minimized, with the projection data and calibration data as inputs, initialized with preliminary estimates of path lengths of the multiple materials to determine final estimates of path lengths of the multiple materials. In a third example of the method, which optionally includes one or more of the first example and the second example, the plurality of materials includes at least three materials.

In yet another embodiment, a system comprises: an x-ray source configured to generate an x-ray beam toward a subject; a detector array comprising a plurality of detector elements configured to detect an x-ray beam attenuated by a subject; and a computing device communicatively coupled to the x-ray source and the detector array, the computing device configured with instructions in a non-transitory memory that, when executed, cause the computing device to: controlling an x-ray source and a detector array to scan a subject with a plurality of x-ray beams at different energy levels and acquire projection data; estimating path lengths of the plurality of materials based on the projection data and calibration data of the x-ray source and the detector array; iteratively refining the estimated path length based on a linearized model derived from the calibration data; and reconstructing a material density image for each of the plurality of materials from the iteratively refined estimated path lengths.

In a first example of the system, the computing device is further configured with instructions in non-transitory memory that, when executed, cause the computing device to estimate path lengths of the plurality of materials based on the projection data and the calibration data without modeling physics of the x-ray source and the detector array. In a second example of the system, which optionally includes the first example, the computing device is further configured with instructions in the 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, which optionally includes one or more of the first example and the second example, the computing device is further configured with instructions in the non-transitory memory that, when executed, cause the computing device to iteratively calculate a final estimate of a path length for each of the plurality of materials, wherein the preliminary estimate of the path length is an initial estimate, wherein the estimated path length used to reconstruct the material density image comprises the final estimate of the path length. In a fourth example of the system, which optionally includes one or more of the first through third examples, the computing device is further configured with instructions in the non-transitory memory that, when executed, cause the computing device to estimate at least two monochromatic sinograms based on the projection data and the calibration data and reconstruct at least two monochromatic images from the at least two monochromatic sinograms. In a fifth example of the system, optionally including one or more of the first example through the fourth example, the system further comprises a display device communicatively coupled to the computing device, and the computing device is further configured with instructions in the non-transitory memory that, when executed, cause 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 through fifth examples, the computing device is further configured with instructions in the non-transitory memory that, when executed, cause the computing device to estimate a plurality of material density images from at least two monochromatic images, the at least two monochromatic images reconstructed from at least two monochromatic sinograms 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 sinograms and an iterative approach using image domain prior information.

As used herein, an element or step recited in the singular and proceeded with the word "a" or "an" should be understood as not excluding plural said elements or steps, unless such exclusion is explicitly recited. Furthermore, references to "one embodiment" of the present invention are not intended to be interpreted as excluding the existence of additional embodiments that also incorporate the recited features. Moreover, unless explicitly stated to the contrary, embodiments "comprising," "including," or "having" an element or a plurality of elements having a particular property may include additional such elements not having that property. The terms "including" and "in … are used as shorthand, language equivalents of the respective terms" comprising "and" wherein ". Furthermore, the terms "first," "second," and "third," etc. are used merely as labels, and are not intended to impose numerical requirements or a particular positional order on their objects.

This written description uses examples to disclose the invention, including the best mode, and also to enable any person skilled in the relevant art to practice the invention, including making and using any devices or systems and performing any incorporated methods. The patentable scope of the invention is defined by the claims, and may include other examples that occur to those of ordinary skill in the art. Such other examples are intended to be 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 insubstantial differences from the literal languages of the claims.