Methods and systems for estimation of basis and multi-basis image reconstruction of imaging data

Abstract

An image reconstruction system includes a memory configured to store image data obtained from an imaging system. The system also includes a processor operatively coupled to the memory and configured to formulate reconstruction of the image data as a non-convex optimization program. The processor also generates a solution to the non-convex optimization problem. The processor also reconstructs the image data to generate a multi-basis image based on the solution to the non-convex optimization problem. A pathlength estimation system includes a memory configured to store projection data obtained from a projection imaging system. The system also includes a processor operatively coupled to the memory and configured to formulate estimation of the pathlength data as a non-convex optimization program. The processor also generates a solution to the non-convex optimization problem. The processor also estimates the pathlength data to generate a multi-basis pathlength based on the solution to the non-convex optimization problem.

Description

Atty. Dkt. No. 05400-0082-PCTMETHODS AND SYSTEMS FOR ESTIMATION OF BASIS AND MULTI-BASIS IMAGE RECONSTRUCTION OF IMAGING DATA CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] The present application claims the priority benefit of U. S. Provisional Patent App. No. 63 / 749,090 filed on January724, 2025 and U. S. Provisional Patent App. No. 63 / 895,162 file on October 7, 2025, the entire disclosures of which are incorporated by reference herein.REFERENCE TO GOVERNMENT RIGHTS

[0002] This invention was made with government support under CA263660 and CA287302 awarded by the National Institutes of Health. The government has certain rights in the invention.BACKGROUND

[0003] Computed tomography (CT) refers to an imaging procedure that uses X-rays to generate cross-sectional images of a subject of interest, such as a human body. A CT scanner generally utilizes a motorized table to rotate an X-ray source around a patient, or other object being scanned. Detectors are used to record the X-rays that pass through the patient, and these recorded X-rays are used to generate cross-sectional images of the patient. These image slices can be stacked together to generate a 3-dimensional image of the patient. Other imaging systems that utilize x-rays includes standard x-ray systems, tomosynthesis imaging systems, laminographic imaging systems, etc.SUMMARY

[0004] An illustrative image reconstruction system includes a memory configured to store image data obtained from an imaging system. The system also includes a processor operatively coupled to the memory and configured to formulate reconstruction of the image data as a non-convex optimization program. The processor also generates a solution to the non-convex optimization problem. The processor also reconstructs the image data to generate a multi-basis image based on the solution to the non-convex optimization problem.

[0005] In one embodiment, the image data includes multiple polychromatic spectra. In another embodiment, the image data originates from an single energy x-ray imaging process or a dual energy x-ray imaging process. In one embodiment, the processor inverts a nonAtty. Dkt. No. 05400-0082-PCTlinear basis path length data model to estimate basis path lengths for the image data, and the solution to the non-convex optimization problem is based at least in part on the estimated basis path lengths. In another embodiment, the processor converts a volume conservation constraint into a constraint on the basis path lengths to estimate the basis path lengths.

[0006] In another embodiment, the processor converts a mass conservation constraint into a constraint on the basis path lengths to estimate the basis path lengths. In one embodiment, the non-convex optimization problem is based on the inverted non-linear basis path length data model. In another embodiment, the processor generates the solution by numerically converging the non-convex optimization problem. In another embodiment, the image data is collected from a single spectrum in an imaging process that utilizes x-rays. In another embodiment, the processor uses a basis-region technique to reduce a number of voxel values in the image data.

[0007] An illustrative method of performing image reconstruction includes storing, in a memory, image data obtained from an imaging system. The method also includes formulating, by a processor operatively coupled to the memory, reconstruction of the image data as a non-convex optimization program. The method also includes generating, by the processor, a solution to the non-convex optimization problem. The method also includes reconstructing, by the processor, the image data to generate a multi-basis image based on the solution to the non-convex optimization problem.

[0008] In one embodiment, the image data includes multiple polychromatic spectra. In another embodiment, the image data originates from a single energy x-ray imaging process or a dual energy x-ray imaging process. In another embodiment, the method includes inverting, by the processor, a non-linear basis path length data model and using the inverted non-linear basis path length data model to estimate basis path lengths for the image data, where the solution to the non-convex optimization problem is based at least in part on the estimated basis path lengths. In another embodiment, the method includes converting, by the processor, a volume conservation constraint into a constraint on the basis path lengths to estimate the basis path lengths. In another embodiment, the method includes converting, by the processor, a mass conservation constraint into a constraint on the basis path lengths to estimate the basis path lengths.

[0009] In one embodiment, the non-convex optimization problem is based on the inverted non-linear basis path length data model. In another embodiment, the method includesAtty. Dkt. No. 05400-0082-PCTgenerating, by the processor, the solution by numerically converging the non-convex optimization problem. In another embodiment, the method includes collecting, by the processor, the image data from a single spectrum in an imaging process that utilizes x-rays. In another embodiment, the method includes using, by the processor, a basis-region technique to reduce a number of voxel values in the image data.

[0010] Other principal features and advantages of the invention will become apparent to those skilled in the art upon review of the following drawings, the detailed description, and the appended claims.BRIEF DESCRIPTION OF THE DRAWINGS

[0011] Illustrative embodiments of the invention will hereafter be described with reference to the accompanying drawings, wherein like numerals denote like elements.

[0012] Fig. 1 A shows physical phantom images in accordance with an illustrative embodiment.

[0013] Fig. IB are graphs showing results at 80kV, 135kV, and De-80&135kV in accordance with an illustrative embodiment.

[0014] Fig. 1C is a chart showing results at 80kV, 135kV, and DE in accordance with an illustrative embodiment.

[0015] Fig. 2 depicts abdomen images in accordance with an illustrative embodiment.

[0016] Fig. 3 shows conventional images reconstructed, respectively, from 60-kVp and 120-kVp data in CBCT in accordance with an illustrative embodiment.

[0017] Fig. 4 depicts images obtained from standard DE data of two distinct spectra (60 and 120 kVp) in CBCT in accordance with an illustrative embodiment.

[0018] Fig. 5 shows physical quantities estimated from standard DE data of two distinct spectra (60 & 120 kVP) in CBCT in accordance with an illustrative embodiment.

[0019] Fig. 6 shows images obtained from data of a single spectral (60kVp) in CBCT in accordance with an illustrative embodiment.

[0020] Fig. 7 shows physical quantities obtained from data of a single spectral (60 kVp) in CBCT in accordance with an illustrative embodiment.Atty. Dkt. No. 05400-0082-PCT

[0021] Fig. 8 shows images obtained from data of a single spectral (120 kVp) in CBCT in accordance with an illustrative embodiment.

[0022] Fig. 9 shows physical quantities obtained from data of a single spectral (120 kVp) in CBCT in accordance with an illustrative embodiment.

[0023] Fig. 10 shows conventional images reconstructed, respectively, from 40-kVp and 100-kVp data in CBCT in accordance with an illustrative embodiment.

[0024] Fig. 11 shows images obtained from standard DE data of two distinct spectra (40 and 100 kVp) in CBCT in accordance with an illustrative embodiment.

[0025] Fig. 12 shows physical quantities estimated from standard DE data of 2 distinct spectra (40 and 100 kVP) in CBCT in accordance with an illustrative embodiment.

[0026] Fig. 13 shows images obtained from data of a single spectral (40kVp) in CBCT in accordance with an illustrative embodiment.

[0027] Fig. 14 shows physical quantities obtained from data of a single spectral (40 kVp) in CBCT in accordance with an illustrative embodiment.

[0028] Fig. 15 shows images obtained from data of a single spectral (100 kVp) in standard CBCT in accordance with an illustrative embodiment.

[0029] Fig. 16 shows physical quantities obtained from data of a single spectral (100 kVp) in CBCT in accordance with an illustrative embodiment.

[0030] Fig. 17 shows standard x-ray images from projection data of differing kVps in accordance with an illustrative embodiment.

[0031] Fig. 18 shows pathlength images from two sets of projection data (40 and 100 kVp) in accordance with an illustrative embodiment.

[0032] Fig. 19 shows pathlength images from two sets of projection data (40 and 100 kVp) in accordance with an illustrative embodiment.

[0033] Fig. 20 shows pathlength images from one set of projection data (40 kVp) in accordance with an illustrative embodiment.

[0034] Fig. 21 shows pathlength images from one set of projection data (40 kVp) in accordance with an illustrative embodiment.Atty. Dkt. No. 05400-0082-PCT

[0035] Fig. 22 shows pathlength images from one set of projection data (100 kVp) in accordance with an illustrative embodiment.

[0036] Fig. 23 shows pathlength images from one set of projection data (loO kVp) in accordance with an illustrative embodiment.

[0037] Fig. 24 depicts a set of 5 spatially complementary basis regions that partitions a full image array in accordance with an illustrative embodiment.

[0038] Fig. 25 depicts basis images of air. water, bone, 20-mg / ml iodine solution, titanium, and stainless steel reconstructed, respectively, from simulated 80- and 140-kV conventional data in accordance with an illustrative embodiment.

[0039] Fig. 26 displays the numerical convergence properties of the dNCPD algorithm in terms of convergence metrics as functions of iteration n from simulated conventional data generated with the 80-kV spectrum in accordance with an illustrative embodiment.

[0040] Fig. 27 depicts the basis regions used for a physical DE-phantom study in accordance with an illustrative embodiment.

[0041] Fig. 28 shows a set of selected basis regions that yields the basis images of visually minimal artifacts for use the final reconstruction of the basis images via the dNCPD algorithm in accordance with an illustrative embodiment.

[0042] Fig. 29 depicts pseudo-codes of the dNCPD algorithm for empirically solving Eq. (67) in accordance with an illustrative embodiment.

[0043] Fig. 30 depicts a computing system for performing image reconstruction in accordance with an illustrative embodiment.

[0044] Fig. 31 depicts a computing system for pathlength estimation in accordance with an illustrative embodiment.DETAILED DESCRIPTION

[0045] In current applications of standard X-ray imaging, data, i.e., X-ray images, are collected with a conventional energy -integrating detector (EID) and also often a single polychromatic spectrum. As used herein, standard X-ray imaging with conventional EID is simply referred to as X-ray imaging, and X-ray imaging with a single spectrum is referred to as single energy (SE) X-ray imaging. While SE X-ray imaging has found a wide array of medical and non-medical applications, it yields X-ray images that necessarily containAtty. Dkt. No. 05400-0082-PCToverlapping artifacts arising from the superimposition of different types of distinct basis materials, such as tissue and bone, and also the beam-hardening artifact resulting from the polychromatic X-ray spectrum. Both artifacts may dampen the utility of X-ray images.

[0046] Thus, there is a need for development of a technique that removes the overlapping artifacts from different types of basis materials and also the beam-hardening effect for improving the utility of SE X-ray imaging in various applications. For example, in chest SE X-ray imaging, effort exists in minimizing the overlapping artifacts of lung tissue and bony ribs for improved detection or visualization of possible lung nodules. Similarly, in SE mammography, there exists an increased interest in separating the contributions of fatty and fibrograndular / cancerous tissues to the mammogram, which can be used for estimating the breast density, an index that is used for assessment of breast cancer risk.

[0047] X-ray images are generally heavily processed, often with non-linear schemes for enhancing the level of visualization-task utility in practical products and applications of X-ray imaging. While methods have been developed for ad hoc estimation of physical quantities such as tissue density from a heavily processed X-ray image such as a clinical SE mammographic image, the estimated quantities can represent only qualitative surrogates of the corresponding physical quantities of interest, with a considerable degree of variability. Furthermore, because the processing schemes can be significantly different among various vendors and research laboratories, the quantities estimated from these images are often associated with a substantial degree of variability as well.

[0048] It is possible to minimize the overlapping artifacts caused by different types of basis materials, such as bone and tissue, by working on raw X-ray images prior to their being processed for visualization. For convenience without causing confusion, a raw X-ray image is referred to as data or simply as an X-ray image hereinafter. A physics-based approach to minimizing the overlapping artifacts from different types of basis materials and beamhardening effect in an X-ray image is to estimate the individual projection images of the basis materials such as tissue and bone within the subject. As used herein, the projection image of a basis material is referred to as the basis pathlength (BPL) image.

[0049] Dual-energy (DE) X-ray imaging has been developed for yielding the BPL images of the basis materials considered from DE data collected with two distinct spectra by use of either kV-switch schemes or dual / multi-layer (or photon-counting) detectors. While some of the DE techniques may yield the BPLs and virtual pathlengths (VPLs) of practical utility,Atty. Dkt. No. 05400-0082-PCTthey provide largely approximate solutions to their respective non-linear BPL-data models, and thus approximate estimations of the BPLs. Moreover, the workflows and / or hardware involved in the DE techniques can be considerably more complex and / or costly in terms of scanning workflow, time, dose, and / or detector, than those in SE X-ray imaging, which takes one scan of a single spectrum with a conventional EID.

[0050] In an effort to avoid the workflow and / or cost issues in DE X-ray imaging, algorithms have been investigated for estimating two BPL images of basis materials such as tissue and bone within the subject from collected SE data (i.e., SE X-ray image). As the algorithms generally approximately solve the relevant non-linear BPL-data models, they yield generally two approximate BPLs from SE data. There is no evidence that traditional algorithms can accurately and stably invert the corresponding non-linear BPL-data models.

[0051] The inventors have developed an algorithm that can accurately and stably invert the non-linear BPL-data models to estimate BPLs from data collected with multiple polychromatic spectra in X-ray imaging, which includes SE and DE X-ray imaging as special cases. Specifically, based upon the standard, non-linear data model of X-ray imaging, one can obtain the non-linear BPL-data models that relate BPLs to data in X-ray imaging. Accurate estimation of the BPLs from data is thus equivalent to accurately inverting the non-linear BPL-data models. The well-known volume-conservation (VC) and mass-conservation (MC) constraints in CT imaging can be converted into constraints on the BPLs, which can then be exploited for augmenting data in X-ray imaging to accurately and stably estimate the BPLs, as discussed below.

[0052] In an illustrative embodiment, the inventors formulated the inverse problem of a non-linear BPL-data model in X-ray imaging as a non-convex optimization program (or problem), which can include the VC, MC, and additional constraints on the BPLs for possibly improving the stability of the inverse problem. The inventors then developed an algorithm, referred to as the dynamic non-convex primal-dual (dNCPD) algorithm, to numerically converge the optimization program or equivalently, to accurately and stably invert the nonlinear BPL-data model for estimating the BPLs from data (i.e., multi-spectral X-ray images in X-ray imaging, which includes SE and DE X-ray imaging as special cases of practical significance).

[0053] Computer-simulation studies were designed and conducted first for verifying the correctness of the dNCPD algorithm and its computer implementation in terms of accuratelyAtty. Dkt. No. 05400-0082-PCTinverting the non-linear BPL-data model for accurate estimation of the BPLs. In the testing, ideal simulated data, i.e.. ideal X-ray images, are generated from digital phantoms of anatomies of clinical and application relevance. Once the dNCPD algorithm and its computer implementation were verified, the inventors conducted studies to evaluate the algorithm’s stability in BPL estimations from real data. The real data includes X-ray images, which are collected from physical phantoms and other subjects in real-world SE and DE X-ray imaging and thus are inconsistent with the non-linear BPL-data models.

[0054] As discussed above, the task of BPL estimations from X-ray images is equivalent to the inversion of the corresponding non-linear BPL-data model. As such, based upon a nonlinear BPL-data model, the inventors formulated the inverse problem (i.e., the BPL-estimation problem) as an optimization program and then developed an algorithm, referred to as the dNCPD algorithm, to accurately and stably invert the non-linear BPL-data model. Specifically, the algorithm estimates the BPLs through numerically converging the optimization program from data collected in X-ray imaging.

[0055] Data models in X-ray imaging are discussed below, including the standard, nonlinear data model. Data collected with ray j from a subject scanned in X-ray imaging can be modeled adequately by the standard, non-linear data model as:g^ = —In I dE exp dl g'(E, r) (1)where g'(E, r) denotes the linear attenuation coefficient (LAC) of energy E at location r, c / j (E) is the effective spectrum s over ray / , which is the product of the filtered X-ray source spectrum and energy response of the EID, and s = 1, 2,and Sj is the total number of distinct spectra (or, equivalently, of data sets acquired with distinct spectra) for ray j. The detector-energy response is proportional to energy E in a conventional EID, or can be of a specific shape for each detector layer of a multilayer detector or within each energy window in a photon-counting detector. The letter J is used to denote the total number of rays for which data are collected, i.e., J = 1, 2,..., J.

[0056] As LAC g'(E, r) is a function of two variables E and r, one can decompose / / (E, r) into a series of the products of two functions of single variables E and r, with a focus on basis material types or interaction types. In this work, a focus is on the LAC decomposition into basis material types.Atty. Dkt. No. 05400-0082-PCT

[0057] Basis images and LAC decomposition are discussed below. Two types of basis images are considered. These basis images bk(f) and bk(r) decompose LAC '(E, r) into the forms of:K(E,r) = ^ pk' (E)bk' (f), and (2) fc=l K= ^ pk(E)bk(r), (3)k=l

[0058] where / (E) and i-ik(E~) denote, respectively, the LAC and mass attenuation coefficient (MAC) of basis material k. and K is the total number of basis materials considered (e.g., tissue and bone). Additionally, / ( ) = Pk k(E and Pkis the mass density of basis material k. Equations (2) and (3) are referred to as the mass- and volume-based decompositions of LAC, respectively.

[0059] As discussed in more detail below, one can obtain the MC constraint on the basis images as:K K^b;)(f) ff(O = b^and^bk^ / Peff(.r) = bt(r); (4)k k and one can also obtain, under Assumption 1 (discussed in more detail below), the VC constraint on the basis images asK K^ bk(r) = bt(f) and ^ bk(r) / pk= bt(r), (5)k k where petf(r) denotes the effective mass density defined in Eq. (27); and b®(r) is known and satisfies b®(r) = 1, or b®(r) = 0, at voxels to which the constraints are, or are not, applied. Clearly, b(t)(r) = 0 at voxels containing no materials, including those outside the exterior boundary of the subject support at which neither VC nor MC constraint is applied.

[0060] BPLs and their constraints are described below. The precise definition of the BPLs depends upon the basis images chosen. For basis images bk(f) and bk(r) in Eqs. (2) and (3), one can define their corresponding BPLs as:dl hl',(r') and djkAtty. Dkt. No. 05400-0082-PCTwhere the line integration is carried out over ray j.

[0061] Again, one can, under Assumptions 1 and 2 (discussed below), obtain the MC constraints on the BPLs as:K KY 'Jdi'k~= dPand dJfc / Peff = 4^' k=l Peffk=1and also the VC constraints on the BPLs asK Kdjk= df* and djk / pk= dd\ (8) k= l k=l where dj^ denotes the total PL over ray j through b^(r), i.e.,dd)= dl b^(r). (9)

[0062] Because b^ r) is assumed to be known, dj^ is also known for j = 1, 2,...,]. Again, it can be observed in Eqs. (7) and (8) that if knowledge of pkand pe^ is available, there are two constraints on each of BPLs djkand dj'k. The constraints in Eqs. (7) and (8) can be used for data augmentation as discussed below. It is also noted that other types of BPLs can be introduced, as revealed by the examples discussed below;

[0063] The non-linear BPL-data model is described below One can first obtain nonlinear BPL-data models for an individual ray from which one can then obtain the non-linear BPL-data models for multiple rays. For simplicity, “non-linear BPL-data model” and “BPL-data model” are used interchangeably hereinafter. As discussed above, the explicit decomposition of LAC p'(E, r) depends upon the selection of basis images. The discussion below focuses on, under Assumptions 1 and 2, obtaining the BPL-data models for djk, with the understanding that the discussion can readily be extended to obtain the BPL-data models for djfc.

[0064] The standard BPL-data model for a single ray is now' described. Substituting Eqs. (3) and (6) into Eq. (1), one obtains—In J dEq^(E) exp[— (10)where d^(E} (i.e., the VPL) is given byAtty. Dkt. No. 05400-0082-PCTKd?\E} = ^ gk(E)djk. (11)fc=i

[0065] Introducing vector dj = (d1;, d2j djK)Tof size K, one can rewrite Eq. (10) as:M= -|nE (<* / )■ <12)m=l where m = 1, 2,..., M are the indices of energy bins. Additionally:( K \<13)fc=l / where gkm= gk(E) evaluated at energy E = Em.

[0066] Furthermore, introducing vector dj = (cZ1;, d2jdjKof size K. one can reexpress Eq. (12) as(14)fc=l where d;is chosen such that the Taylor series expansion of gj is performed with respect to dJk. Additionally:h(s)(dG)) =denotes the linear term of the Taylor expansion. Also, Ag^^dj, d;) includes the nonlinear component of the BPL-data model in Eq. (12), which is defined as:A^(d7, dy) = gf - <(dO))dyfc. (16)

[0067] The non-linear BPL-data models for a single ray, including the data model with no constraints and the data model with VC constraints, MC constraints, and both VC and MC constraints are considered below. The non-linear BPL-data model with no constraints is discussed first. Using vectors g}- = (g^,g^2.....g ^ )T,gj =(g^, gf*,..., g^,Atty. Dkt. No. 05400-0082-PCTand Agj = (Ag^, Ag^ Ag^ )Tof size Sj, one can write the BPL-data model in Eq. (14) in a vector-matrix form as follows:gj = gj + Agj(dp dj) (17)9 J = where matrix ^(d,) of size (S) + 1) X K is given by:h“(d,)... k»(d()\ / .'?(<!,) / .'?(<!;)... k®(d,) (18)

[0068] The non-linear BPL-data model with the VC constraint is also considered. It can be observed that the VC constraints in Eqs. (7) and (8) on the BPL of basis image bk(f) require knowledge of effective density peff of the subject and mass density pkof basis material k, respectively. Because pe^ is generally unknown and because pkis known often for chosen basis material k, included herein is a discussion of the application of the VC constraint in Eq. (8) to the BPL of basis image bk(r) for obtaining the non-linear BPL-data model below.

[0069] Enforcing the VC constraint on djk in Eq. (8), one can obtain for ray j an augmented non-linear BPL-data model in a matrix-vector form of Eq. (17), with vectors 9j = g^\gj2}. ■ ■ ■, gS}d^,gj= g^,gj(2),..., g(.S}d^)T. andAgj=( g^, Ag^ Ag^, 0)Tdenote of size Sj + 1 andM1>(4) '■“(4) - '>“(4)\ '■“(4) ^’(4) - C(4) ^(4) = (19) / >ty(d,.) / >ty(d;)... / .^>(d;) \ pp... pp J where = X / pk, and k — 1, 2,... K.Atty. Dkt. No. 05400-0082-PCT

[0070] The BPL-data model with the MC constraint is considered next. If the MC constraint on d]kin Eq. (7) is available (i.e.. if Assumption 2 is satisfied,) it can, similar to the VC constraint, be exploited also to augment data for accurate, stable BPL estimation. In this case, the BPL-data model with the MC constraint can readily be obtained from that with the VC constraint simply by replacing f3k^ = 1 / pkwith = l / peff in the last row of matrix ^(d,) in Eq. (19).

[0071] The BPL-data model with both VC and MC constraints is also considered. If both VC and MC constraints in Eqs. (7) and (8) on djkare available, they can then be applied simultaneously to augmenting data. In this case, the BPL-data model with both VC and MC constraints can be obtained in the form identical to that of Eq. (17), but with vectors gj = (.9^: 9^. 9? ■■ -.i = (J^L ’ ^9^. Adj, 0, 0)Tof size S, + 2 and matrix %,■(<!,■) of size (S,- + 2 j x K / .“(dj... / .“(dW '•li’A)... k®(d;) =, ■ ■, (20) ’(d,) '■fl id,.)... / £' (<• / )... / ?«\ pp’... / J';"" / where [3k^ = X / pkwith (3km^ = l / peff for the volume- and mass-based decompositions, and k = X,2,... K.

[0072] It is noted that the BPL-data model in Eq. (17) is for individual ray j, involving arbitrary K and Sj. Discussed below are the imaging cases with specific selections of K andSJ

[0073] Imaging condition with K < Sj. Under this imaging condition, a BPL-data model is given by Eqs. (17) and (18), and it can in general be inverted for accurate and stable estimation of K BPLs from data acquired with Sj distinct spectra, without invoking either of the VC and MC constraints. For example, it is w ell-established that in DE X-ray imaging, 2 BPLs can be accurately and stably estimated from DE data collected with 5 = 2, or Sj > 2, distinct spectra.Atty. Dkt. No. 05400-0082-PCT

[0074] Imaging condition with K > Sj. This imaging condition is of high practical interest because the number, i.e., K, of BPLs to be estimated is larger than the number of data sets, Sj. Three specific cases under this imaging condition are described below.

[0075] For the case with K = Sj + 1, one may accurately and stably invert the BPL-data model in Eq. (17) with a constraint, e.g.. the VC constraint in matrix Eq. (19) to estimate K BPLs only from data acquired with 5} distinct spectra. For example, it is possible to invert the BPL-data model with a VC constraint for accurate and stable estimation of K = 2 BPLs from SE data in SE X-ray imaging, as demonstrated in the numerical studies described below. However, it is unlikely to accurately estimate K = 2 BPLs only from SE data in SE X-ray imaging without invoking the VC or MC constraint.

[0076] For the case with K = Sj + 2, the BPL-data model in Eq. (17) can incorporate both VC and MC constraints, if available, for data augmentation so that the BPL-data model may accurately and stably be inverted for K BPLs from data collected only w ith Sj distinct spectra. In particular, if K = 3, it may be possible to accurately invert the BPL-data model in Eq. (17) incorporated with both VC and MC constraints, provided both constraints are available, to estimate K = 3 BPLs from only SE data collected in X-ray imaging. However, it is unlikely to accurately estimate K = 3 BPLs by invoking only either the VC or MC constraint.

[0077] For the case with K > Sj + 3, it remains unlikely that the BPL-data model of Eqs. (17) and (18) can always be inverted accurately and stably even if the VC and MC constraints can be applied. For example, it is generally unlikely to accurately estimate K = 4 BPLs from only SE data of SE X-ray imaging.

[0078] An alternative application of the constraints is described below. In the discussion above, the VC constraint is incorporated into the BPL-data model through the matrix in Eq. (19), while the number K of unknow n BPLs remains unchanged. On the other hand, one can exploit the constraints for, instead of augmenting data, reducing the number of unknow n BPLs. Included below is an alternative approach to incorporating the VC constraint into a BPL-data model. The approach can readily be extended to MC constraint or both VC and MC constraints.

[0079] Enforcing the VC constraint on cljkin Eq. (8), one can re-express Eq. (10) as:Atty. Dkt. No. 05400-0082-PCTK9jS)= ~lnJ dEq^s\E)exp -^' {ik(E')djk, (21) fc=i

[0080] where ^(E) = qf\E) exp[-pK(E)pKdf\pk(E') = pk(E) - pk(E^,and K = K — 1. It can be observed that Eq. (21) is in a form exactly identical to that of Eq. (10), thus reducing the number of unknown BPLs from K to K. Therefore, following the same strategy discussed above,, one can readily obtain a BPL-data model in the form of Eqs. (17) and (18) for K BPLs djk.

[0081] The non-linear BPL-data models for multiple rays are also considered. Based upon the BPL-data model in Eq. (17) for individual ray j, one can obtain the BPL-data model for multiple rays entailed in a single view or multiple views of X-ray imaging. Introducing vectors, g = (g[, gl.., gj )T, g = (gT, gj,.. gj )T, and Ag(d, d)T=(AgiCd d, )T, Ag2(d2, d2),..., Ag;(d;, d;)T)T, which are of size 2=15)-, one can combine the BPL-data models of Eq. (17) of all individual rays involved to form the BPL-data model of multiple rays as:g = g + Ag(d, d) (22) g = 7L(d)d,where vector d = (d^, d2,..., d )Tis of size K x J and are referred to as BPL images, and 7f(d) a block diagonal matrix of size (S;=i ■$) + / ) x (J x K) given by:(23)

[0082] It is noted that the BPL-data models for individual BPL djkof single ray j, and for multiple BPLs d of multiple rays, are based upon basis image bk(r). As mentioned, the same approach can straightforwardly be applied to obtaining the BPL-data models for BPLs dj'kof basis image bk(r).

[0083] Inversion of the non-linear BPL-data models is discussed below, along with optimization programs and development of the dNCPD algorithm. In an attempt to estimateAtty. Dkt. No. 05400-0082-PCTBPL images d from knowledge of through inverting the BPL-data model in Eq. (22), the inventors formulated an optimization program as:d* = arg min d>(g(M\g(d)) s.t. ‘P(d), (24)d where g(M)denotes the measured data for all of the rays, (g(M), g(d)) denotes an objective function, and (d) is a set of constraints.

[0084] One can choose, for example,1d* = arg min- ||g(M)- %(d)d - Ag(d, d)||^ s.t. TV(d) < y & d > 0, (25) d 2 where the objective function in Eq. (25) is non-convex due to the inclusion of Ag(d, d). The objective function is also non-convex due to a TV constraint on a BPL image in X-ray imaging or on multiple BPL images at multiple views in tomosythesis or CT imaging, as well as possibly on BPL images of the basis images. Also, the positivity constraint is on BPLs for each individual ray. Clearly, both constraints are convex, and additional, relevant constraints can be devised as needed.

[0085] The dNCPD algorithm was developed to numerically converge the non-convex optimization program in Eq. (25) to invert accurately and stably the BPL-data model in Eq. (22) for estimating BPLs in X-ray imaging. As Eq. (25) is non-convex because its objective is non-convex, this provides the basis to form a strategy to develop an algorithm to solve the non-convex optimization program in Eq. (25). The strategy was also devised to develop an algorithm for numerically converging a non-convex optimization program involved in reconstruction of basis images in CT imaging.

[0086] The PD algorithm Analysis of Eq. (25) reveals that its non-convexity results from the dependence of term Ag(d, d) on d, which is to be estimated. In the strategy mentioned, one can thus convexify Eq. (25) by replacing Ag(d, d) with constant Ag independent of d. Subsequently, one can readily obtain a primaldual (PD) algorithm to mathematically solve the convexified version of Eq. (25). While Eq. (25) and its convexified version mathematically differ from each other, the form of the latter is largely identical to that of the former. This observation provided motivation for the development of the dNCPD algorithm for solving Eq. (25) through tailoring the derived PD algorithm.

[0087] In an illustrative embodiment, the dNCPD algorithm can be obtained by making appropriate changes in the PD algorithm or in essence to its pseudo codes in which one canAtty. Dkt. No. 05400-0082-PCTuse d(n) to denote the estimation of d at iteration n. Specifically, the dNCPD algorithm and its pseudo-codes are obtained by making changes to the PD algorithm and its pseudo-codes, which include: (1) d, Jf (d) are replaced by d^ and (d^), and (ii) Ag is replaced byAg(d^n+1\ d(n)), which can be computed by using Eq. (22). It can be observed that the pseudo-codes of the PC and dNCPD algorithms entail virtually identical iterative procedures, suggesting that the numerical convergence characteristics of the dNCPD algorithm can be similar to the well-characterized characteristics of the PD algorithm. It is noted that the mathematical convergence of the dNCPD algorithm is not established in solving the non- convex optimization program in Eq. (25). However, as the numerical studies reveal, the dNCPD algorithm appears to numerically converge to Eq. (25). and thus invert the BPL-data model in Eq. (22), in terms of single- or double-precision floating-point format.

[0088] The inventors designed and conducted two types of studies to verify and evaluate the accuracy and stability of the NCPD algorithm, which was developed to invert the BPL- data model in Eq. (22), along with Eq. (23), for estimating the BPL images in X-ray imaging. Specifically, one can consider the inversion of the BPL-data model in Eq. (22), along with Eq. (23) that is specified by Eq. (18) with neither VC nor MC constraint or by Eq. (19) with the VC constraint. It is noted that the numerical studies can readily be performed if the MC constraint is available. In the verification studies, the inventors used the ideal simulated data that are generated from digital phantoms with anatomies of application relevance by using the data model in Eq. (10), without including any additional physical factors such as noise. Therefore, ideal simulated data are consistent completely with the data model and can be used to verify the accuracy of the dNCPD algorithm and its computer implementation in inverting the BPL-data model in Eq. (22), along with Eq. (23). Furthermore, in the real data studies, the inventors collected real data from physical phantoms and other subjects in X-ray imaging. As real data contain real-world physical factors such as noise, scatter, and other physical factors, which are not included in (i.e., are inconsistent with) the BPL-data model, they are used to evaluate the the stability of the dNCPD algorithm verified.

[0089] The study results in the simulated- and real-data studies reveal that quantitatively accurate BPLs and VPLs can be estimated through accurate and stable inversion of the nonlinear BPL-data model in in X-ray imaging.

[0090] Based upon the analysis of the VC and MC constraints and the derivation of the respective BPL-data models, the inventors formulated the inverse problem of the BPLdataAtty. Dkt. No. 05400-0082-PCTmodels as a non-convex optimization program. The inventors then developed the dNCPD algorithm to numerically converge the non-convex optimization program to accurately and stably invert the BPL-data models, yielding BPLs from data collected in X-ray imaging. In particular, exploiting the VC constraint, the dNCPD algorithm can accurately and stably invert the corresponding BPL-data models to yield two or three BPLs, and the VPLs, from SE or DE data in simulated X-ray imaging. While the dNCPD algorithm is developed for accurately inverting the non-linear BPL-data models, it is motivated by the dNCPD algorithm for reconstruction of multibasis images only from SE data in conventional CT.

[0091] In the verification study, SE and DE data of X-ray imaging were simulated using digital phantoms designed in the BPL-data models. The digital phantoms are therefore consistent with the BPL-data model, and they enable the verification of the correctness (i.e., accuracy) of the dNCPD algorithm and its implementation. Conversely, the inventors collected real SE and DE data from physical phantoms and other subjects in real-world X-ray imaging and used them to evaluate the stability of the dNCPD algorithm in estimation of the BPLs. Both verification and evaluation studies reveal that the dNCPD algorithm can accurately and stably invert the BPL-data models to estimate two, or three, BPLs, as well as VPLs, of interest, respectively, from SE or DE data in X-ray imaging.

[0092] While the evaluation studies demonstrate the stability of the dNCPD algorithm in real-data studies with data containing physical factors such as noise and scatter, the accuracy level of other physical quantities such as parameterand can also understandably impact the BPL estimations with the dNCPD algorithm. The inventors have conducted evaluation studies with different estimates of the physical quantities. These study results confirm the stability- of the dNCPD algorithm, relative to the variability of the physical quantities, in estimation of the BPLs and VPLs from SE or DE data in X-ray imaging.

[0093] As the dNCPD algorithm yields BPLs and VPLs without a beam-hardening effect, the BPLs estimated for all rays in tomosynthesis or CT imaging can subsequently be used for reconstructing of basis images and the virtual monochromatic images (VMIs) of the subject scanned. The dNCPD algorithm can also be applied to estimating BPLs from multi-spectral data collected with dual / multi-layer and photon-counting detectors.

[0094] As discussed above, the inventors have developed and evaluated the dNCPD algorithm to accurately and stably invert the BPL-data model for estimating BPLs in X-ray¬ imaging, without changing scan workflow, increasing imaging dose and time, or usingAtty. Dkt. No. 05400-0082-PCTadvanced detectors (dual / multi-layer and photon-counting detectors.) In particular, the dNCPD algorithm may be exploited for yielding BPL images of, e.g., tissue and bone, in chest X-ray or iodine contrast. The algorithm can also be used for normal breast tissue in contrast-enhanced mammography of fibrograndular / cancerous and adipose tissues in mammography, and for yielding BPL images of tissue, bone, iodine, etc. in DE X-ray imaging. As the BPL images and VPL images are indeed additions to, instead of the replacement of, the corresponding conventional X-ray images, a natural adoption of the approach in X-ray imaging may be plausible for yielding additional information of possible clinical or other application utility.

[0095] Using the BPLs, VPLs, and additional pathlengths of various parameter images such as, but not limited to, images of material densities, contrast concentrations, and physical quantities such as (effective) electron density and (effective) atomic numbers, one can readily obtain tomographic, tomosynthesis, and laminographic images. These images obtained from BPLs, VPLs, and other additional pathlengths are estimated by use of various existing and / or newly developed reconstruction algorithms from which one can readily yield various parameter images such as images of material densities, contrast concentrations, and physical quantities such as (effective) electron density and (effective) atomic numbers. The tomographic, tomosynthesis, laminographic images, and parameter images mentioned above can be 2D (with 2D spatial dimensions,) 3D (with 3D spatial dimensions,) 3D (with 2D spatial dimensions and a temporal dimension,) 4D (with 3D spatial dimensions and a temporal dimension,) and nD (with 2D or 3D spatial dimensions, a temporal dimension, and multiple (> I) parameter dimensions) dimensions. In addition to the dNCPD algorithm, it was also found that a non-convex ADMM solver or Al / deep learning approaches can be developed to to accurately and stably invert the BPL-data model for estimating BPLs in X-ray imaging, without changing scan workflow, increasing imaging dose and time, and using advanced detectors (dual / multi-layer and photoncounting detectors.)

[0096] Included below is discussion of figures that support the above-described embodiments. In the figures, various acronyms / abbreviations are used. The acronym CBCT refers to cone-beam computed tomography, CT refers to computed tomography, and DE refers to dual energy7. Dual energy7refers to the fact that in standard dual-energy7CT, data are collected with two distinct, effective X-ray spectra largely determined by the X-ray source voltage, often referred to as kVp or kV data, and images reconstructed from the data are referred to as kVp or kV images. The acronym LAC is linear attenuation coefficient, which isAtty. Dkt. No. 05400-0082-PCTa function of energy. The acronym VMI is a virtual monochromatic image, which represents the spatial distribution of LACs and is thus a function of energy as well, which can be written as VMI at some keV.

[0097] Figs. 1 and 2 show results of image reconstruction and quantitative estimation of LAC from data of diagnostic CT. Specifically, Fig. 1 A shows physical phantom images in accordance with an illustrative embodiment. Fig. IB includes graphs showing results at 80kV, 135kV, and De-80&135kV in accordance with an illustrative embodiment. Fig. 1C is a chart showing results at 80kV. 135kV. and DE in accordance with an illustrative embodiment. Fig. 2 depicts abdomen images in accordance with an illustrative embodiment.

[0098] Figs. 3-9 show results of image reconstruction and physical quantity estimation from phantom data of cone-beam CT. The physical quantities estimated include electron density and effective atomic number. Fig. 3 shows conventional images reconstructed, respectively, from 60-kVp and 120-kVp data in CBCT in accordance with an illustrative embodiment. Fig. 4 depicts images obtained from standard DE data of two distinct spectra (60 and 120 kVp) in CBCT in accordance with an illustrative embodiment. Fig. 5 shows physical quantities estimated from standard DE data of two distinct spectra (60 & 120 kVP) in CBCT in accordance with an illustrative embodiment. Fig. 6 shows images obtained from data of a single spectral (60kVp) in CBCT in accordance with an illustrative embodiment. Fig. 7 shows physical quantities obtained from data of a single spectral (60 kVp) in CBCT in accordance with an illustrative embodiment. Fig. 8 shows images obtained from data of a single spectral (120 kVp) in CBCT in accordance with an illustrative embodiment. Fig. 9 shows physical quantities obtained from data of a single spectral (120 kVp) in CBCT in accordance with an illustrative embodiment.

[0099] Figs. 10-16 show image reconstruction and physical quantity estimation from mouse data of cone-beam CT. Specifically, Fig. 10 shows conventional images reconstructed, respectively, from 40-kVp and 100-kVp data in CBCT in accordance with an illustrative embodiment. Fig. 11 shows images obtained from standard DE data of two distinct spectra (40 and 100 kVp) in CBCT in accordance with an illustrative embodiment. Fig. 12 shows physical quantities estimated from standard DE data of 2 distinct spectra (40 and 100 kVP) in CBCT in accordance with an illustrative embodiment. Fig. 13 shows images obtained from data of a single spectral (40kVp) in CBCT in accordance w ith an illustrative embodiment. Fig. 14 shows physical quantities obtained from data of a single spectral (40 kVp) in CBCTAtty. Dkt. No. 05400-0082-PCTin accordance with an illustrative embodiment. Fig. 15 shows images obtained from data of a single spectral (100 kVp) in standard CBCT in accordance with an illustrative embodiment. Fig. 16 shows physical quantities obtained from data of a single spectral (100 kVp) in CBCT in accordance with an illustrative embodiment.

[0100] Figs. 17-23 show pathlength images of basis images, VMIs, and physical quantities from mouse projection data. Specifically, Fig. 17 shows standard x-ray images from projection data of differing kVps in accordance with an illustrative embodiment. Fig. 18 shows pathlength images from two sets of projection data (40 and 100 kVp) in accordance with an illustrative embodiment. Fig. 19 shows pathlength images from two sets of projection data (40 and 100 kVp) in accordance with an illustrative embodiment. Fig. 20 shows pathlength images from one set of projection data (40 kVp) in accordance with an illustrative embodiment. Fig. 21 shows pathlength images from one set of projection data (40 kVp) in accordance with an illustrative embodiment. Fig. 22 shows pathlength images from one set of projection data (100 kVp) in accordance with an illustrative embodiment. Fig. 23 shows pathlength images from one set of projection data (loO kVp) in accordance with an illustrative embodiment.

[0101] Included below is additional discussion regarding the mathematical analysis of basis-image and BPL constraints. If one lets (E, r) denote the MAC of the compounds or intimate mixtures of elements, one can relate the MAC to the corresponding LAC through:'(E,r) = (E,r)pef{(r), (26) whereeff(r), the effective density of the compounds or element mixtures, is given by:Peff GO =anddM(r) = ^5mk(r). (27)

[0102] The term zlM(r) depicts the total mass within volume AV such as the volume of the image voxel considered, of the compounds or element mixtures, 5mk(r) is the mass of material k, at r, and K denotes the total number of basis materials involved.

[0103] Applying the Bragg rule to compounds or intimate mixtures of elements, one can write the MAC r(E, r) as:K(m) v(E,r) = ^pk(E)ak (f), (28)k where akmr)the mass fraction of basis material k given byAtty. Dkt. No. 05400-0082-PCT(m)z^ _5mk(f) _ 8mk(r) “k WW) S^mk(r)'( }

[0104] Mass-based LAC decompositions are also considered. Substituting. Eq. (28) into Eq. (26), one obtains a mass-based decomposition of LAC, given by:Kr) = pfe(E)a^m)(r)Peff (r). (30)fc=i

[0105] It is noted that once basis material k is selected, its density pkis a known constant independent of r, and one can rewrite Eq. (30) as:p'(E,r) = y ^(£)«tm)(^)Peff^, (31)pkwhich constitutes another mass-based decomposition of the LAC.

[0106] With respect to volume-based LAC decompositions, Assumption 1 is introduced as follows:K AV = ^ Sv^r), (32)k where 5vk(r) denotes the volume of material k within AV. e.g., of a voxel, at r. The mass densify of material k within 8vk(r) can be written as:8mk(r)Pk (33) 8vk^ ’which is intrinsic to basis material k and thus independent of r.

[0107] Using Eqs. (27) and (32) in Eq. (30). one obtains, under Assumption 1, a volumebased decomposition:K p'(E,r) = ^ pk(E)a^(f)pk, (34)k=l where akv)(r)pkdenotes the volume fraction of material k and can be expressed as:8vk(f) 3vk(r)(35) “k (r)~ AV ~ %8vk(rY

[0108] Furthermore, one can rewrite Eq. (34) as:Atty. Dkt. No. 05400-0082-PCTK H'(E,r) = ^ pk' (E}a^{f), (36) / c=i which also constitutes another volume-based decomposition of LAC under Assumption 1.

[0109] Using Eqs. (29) and (35), the MC and VC constraints can be represented, respectively, as:K K^ t4m)(r) = l and ^ ^v)(r) = l. (37) fc=i fc=i

[0110] Constraints on basis images are considered below. Comparing Eq. (31) with Eq. (2), and Eq. (30) with Eq. (3), one obtains:b' (O = aJm)(r)^^, (38)Pk andbk(r) = a^\r)peff(r\ (39) where b / ((r) is dimensionless and bk(r) is of dimension cm 'g.

[0111] Furthermore, under Assumption 1, comparing Eq. (36) with Eq. (2), and Eq. (34) with Eq. (3), provides:bk(r) = (40) andbk(f) = afcV)(r)pk, (41) where, again, bk' ( ) is dimensionless and bk(r) is of dimension cm 'g.

[0112] Applying the MC constraint to basis images in Eqs. (38) and (39), respectively, one can obtain the MC constraints on basis images bk(f) and bk(r) as:K Kybk(f)Pkf^ = bt(r) and V bk(r) / pe{{(r) = bt(r); (42) k Peff r) kwhereas the application of the VC constraint to basis images in Eqs. (40) and (41), respectively, yields the VC constraints on basis images bk'(r) and bk(f) as:Atty. Dkt. No. 05400-0082-PCTK K^ bk(r) = bt(r) and bk(r) / pk= Z?t(r), (43)k k where b^(r) = 1, or b^ r) = 0, at voxels to which one of the constraints is, or is not, applied. Clearly, b(t)(r) = 0 at voxels containing no materials, including those outside the exterior boundary of the subject support.

[0113] With respect to BPL constraints. Assumption 2 is that total density peff is independent of r. Under Assumptions 1 and 2, considering Eqs. (42) and (43), one can obtain the MC constraints on BPLs d!'kand dJkdefined in Eq. (6) as:K KXdi’kp^= and dJ'k / pe{{ = d?’ (44)fc=lPeffk=l and also the VC constraints on BPLs djkand djkasK Kdj'k= d^ and djk / pk= d?\ (45) fc=l k=l where d^ denotes the total PL over ray j through b^ (r), i.e.,dj(t)= [ dl b (r). (46) ^3

[0114] Because b^ (r) is assumed to be known, dj^ is thus known for j = 1, 2.,..., J. Again, it can be observed in Eqs. (44) and (45) that if knowledge of pkand peffis available, there are two constraints on each of BPLs and djkand dj'k. The constraints in Eqs. (44) and (45) can be used for data augmention as discussed below.

[0115] It is noted that other types of BPLs can also be defined, including(47)djk' = f dl bk(r) (48)lJd = [ dl bk(r) / peff(r) (49)Atty. Dkt. No. 05400-0082-PCTd^ = dl bk(r) / pk. (50)lJ

[0116] Comparison of Eqs. (6) and (48) indicates that dj'k= d^'. It is also possible to derive constraints on the BPLs and their data models, and the dNCPD algorithm may be tailored to invert these BPL-data models to estimate the BPLs.

[0117] Determination of constraint information is discussed below, including the approaches to determining constraint information, i.e., knowledge of b^\r) and d^.

[0118] Determination of Z ^(r) of a subject. (1) ^(r) can be determined from the subject’s 3D or 2D images, which are obtained from its CT, tomosynthesis, or laminogram scans with X-rays or gamma rays and / or from other scans such as MRI and optical scans. Terahertz imaging, mm-wave imaging, and X-ray backscatter imaging or a combination of the techniques / methods / systems mentioned above. (2) b(t)(f) can be measured directly with various techniques such as the mechanical profilometer scans, laser and optical scans, Terahertz scans, mm-wave scans, and X-ray backscatter scans or scans with combinations of the techniques / methods / systems mentioned of the subject surfaces; (3) b^(r) can be estimated by using simulation data of phantoms that are of similar anatomy of the subjects; and (4) b^(r) can be estimated from an existing database of 2D (with 2D spatial dimensions,) 3D (with 3D spatial dimensions,) 3D (with 2D spatial dimensions and a temporal dimension,) 4D (with 3D spatial dimensions and a temporal dimension,) and nD (with 2D or 3D spatial dimensions, a temporal dimension, and multiple (> 1) parameter dimensions) images of CT, MRI, tomosynthesis, laminogram, ultrasound, photoacoustic, Terahertz, any tomographic images, or the combined tomographic images of subjects by use of Al technologies.

[0119] Determination of d ^of a subject. (1) d ^can readily be determined from b^(r) of the subj ect that is determined / estimated already as described immediately above; (2) d^ can be measured directly with various techniques such as experimental measurements and / or laser and optical scans of the subject surfaces; (3) d^ can be estimated by using simulation data of phantoms that are of similar anatomy of the subjects; and (4) d ^can be estimated from an existing database of b^(r) and / or d^ by use of Al technologies.Atty. Dkt. No. 05400-0082-PCT

[0120] Obtaining tomographic, tomosynthesis, and laminographic images is discussed below. From the BPLs, VPLs, and additional pathlengths of various parameter images such as images of material densities, contrast concentrations, and physical quantities such as (effective) electron density and (effective) atomic numbers, one can readily obtain tomographic, tomosynthesis, and laminographic images. The tomographic, tomosynthesis, and laminographic images and also the parameter images mentioned above can be 2D (with 2D spatial dimensions,) 3D (with 3D spatial dimensions,) 3D (with 2D spatial dimensions and a temporal dimension,) 4D (with 3D spatial dimensions and a temporal dimension,) and nD (with 2D or 3D spatial dimensions, a temporal dimension, and multiple (> 1) parameter dimensions) dimensions.

[0121] Estimation of the effective concentrations and densities of basis materials is also considered. It is of interest to estimate effective concentration and effective mass density of basis material k at r. One can use:(r) = «fcV)(r)Q (r), (51) where Ck(r) and Ck(r) can denote the effective and actual concentration or mass density of basis material k.

[0122] For example, using Eq. (41), one can re-express Eq. (51) as:bk(r)Ck(r} = ^^Ck(fY (52)Pk

[0123] If Ck(r) is chosen as the concentration pk(r) of basis material k, one can then obtain the effective concentration of basis material k at r as:bk(r) T)k(f) = -?jk(r). (53)Pk

[0124] If one chooses Ck(r) as mass density pkof basis material k. one can then obtain the effective density of basis material k at r as:Pk(f) = bk(r). (54)

[0125] It is noted that basis image bk(r) defined in Eq. (41) can be interpreted simply as the effective mass density of basis material k at r (see Eq. (54).) If multiple basis materials are present at r, this means that 0 < a^(f) < 1 and 0 < pk(r) = bk(r) =Atty. Dkt. No. 05400-0082-PCTakr)pk< pk, whereas if only basis material k is at r, this means that ak(r) = 1 and p (r) — bk(r) — a r pk— pk. It can also be observed that:^k(r) = - i]k(r)- (55)Pk

[0126] One can also define a total effective quantity' of all the basis materials as:K Kbk(r) C(r)= =Ck(r). (56)k=l k=l Pk

[0127] If Ck(r) is chosen to be concentration rik(r) of basis material k. C(r) thus denotes the total effective concentration of all basis materials at r. However, if if Ck(r) is chosen mass density r / k(rj of basis material k. C(r) then denotes the total effective mass density of all basis materials at r.

[0128] Using K' to denote the number of basis materials in a subset of the entire set of the K basis materials, where K' < K. one can also define an effective quantity of the selected set of the basis materials as:C'(r) = Q(r) = ^ ^Ck(r), (57)k~i k~iPkwhere it is assumed that, without loss of generality, the indices of the subset of basis materials considered are 1, 2,, K’. If Ck(r) is chosen to be concentration pk(r) of basis material k, C’(r) thus denotes the total effective concentration of the basis materials in the subset at r. If Ck(r) is chosen as mass density r / k(r) of basis material k. C'(r) then denotes the total effective mass density’ of the basis materials in the subset at r.

[0129] Furthermore, it is straightforward to estimate, from knowledge of basis images reconstructed, the effective electron density pe(r) and the effective atomic number Zeff(r) within the subject imaged. One can also estimate maps (or images) of other physical quantities of interest from knowledge of the estimated basis images.

[0130] Estimation of spectra is also considered. Spectra can be obtained by using spectrometer-based direct measurement methods or indirect measurement methods using techniques such as the step wedge technique. For example, in a direct measurement method, an energy-resolving detector is positioned in front of the X-ray beam and connected to aAtty. Dkt. No. 05400-0082-PCTmulti-channel analyzer for recording the energies of individual photons. In this scenario, a generated histogram of photon counts is a function of energy. Also, for example, using the indirect step wedge technique, a calibrated aluminum step wedge can be placed in the X-ray beam path, and the transmitted intensity through each operation can be measured using a nonenergy-resolving detector. By comparing these transmitted signals to those of the unattenuated beam, the spectrum can then be estimated. The direct and indirect methods provide cross-validation of the measured spectra.

[0131] Spectra can also be estimated using software packages, including, but not limited to, TASMIP and SpekCalc. These packages generate spectra from target and filter materials, including, but not limited to, tungsten and aluminum along with appropriate input parameters including anode angle, tube voltage, filter material, and filter thickness.

[0132] Thus, described herein are systems and methods for estimating basis pathlengths (BPLs) in X-ray imaging. An algorithm has been developed to accurately and stably invert non-linear BPL-data models to estimate the BPLs in X-ray imaging, with an eye to the estimations of 2 BPLs or 3 BPLs from data collected with a single spectrum or two distinct spectra using a conventional digital X-ray detector. Virtual PLs (VPLs) can also be obtained from the BPLs at any X-ray energy.

[0133] The approach taken by the inventors was to formulate the inverse problem of a non-linear BPL-data model as a non-convex optimization problem and subsequently develop an algorithm, referred to as the dynamic non-convex primal-dual (dNCPD) algorithm, to numerically solve the non-convex optimization program and thus invert the non-linear BPL-data model for yielding BPLs from data collected in X-ray imaging. The inventors developed the dNCPD algorithm and then carried out numerical studies to verify and evaluate, respectively, the accuracy and stability of the dNCPD algorithm using simulated data and real data collected in real-world X-ray imaging. The results of the studies confirm that the dNCPD algorithm can, under conditions of practical relevance, accurately and stably invert the non-linear BPL-data models to yield BPLs and VPLs in X-ray imaging.

[0134] The dNCPD algorithm can accurately and stably invert the non-linear BPL-data models to estimate BPLs in X-ray imaging. It can be applied in particular to estimating 2 BPLs of e.g.. tissue and bone only from data collected with a single spectrum in X-ray imaging, without using data collected in multi-spectral scans or / and with dual / multi-layer or photon-counting detectors.Atty. Dkt. No. 05400-0082-PCT

[0135] Included below is additional background and description of the proposed dNCPD algorithm. Also included below are additional studies and considerations used to evaluate and confirm the efficacy of the algorithm.

[0136] As discussed above, in standard computed tomography (CT), the polychromatic nature of the X-rays involved leads to a non-linear data model that relates data collected to the image, i.e., the spatial distribution of LAC within the subject of interest. When the reconstruction of an image from CT data is based upon a linear model, it results in beamhardening (BH) artifacts that may plague the visualization, and bias the quantification, of the image reconstructed. Furthermore, as the subject is composed generally of multi-basis materials, it can thus be of practical interest to reconstruct the spatial distributions of the basis materials as basis images.

[0137] The current approach is to reconstruct multi (> 2)-basis images, which may also address the issue of BH correction, from two sets of data collected with two distinct spectra in dual-energy CT (DECT), which are referred to simply as dual energy (DE) data hereinafter. Image-domain-decomposition and data-domain-decomposition methods are used for obtaining basis images from DE data. Additionally, optimization-based algorithms have been investigated for reconstruction of two- or multi (> 2)-basis images directly from DE data through inverting the non-linear data model.

[0138] The proposed dNCPD algorithm is a one-step algorithm for directly inverting the non-linear data model to numerically accurately and stably reconstruct multi-basis images from conventional data collected in standard CT. The inverse problem, i.e., reconstruction problem of multi-basis images, from conventional data can be considerably ill-posed as the result of limited information contained in conventional data in standard CT. In an attempt to alleviate the degree of ill-posedness of the inverse problem, the proposed system exploits, on one hand, the basis-region technique for reducing the number of unknowns involved in basis images and, on the other hand, the volume-conservation (VC) constraint for augmenting conventional data. Using the technique and constraint, the inventors formulate the inverse problem, i.e., the reconstruction problem of multi-basis images, as a non-convex optimization program. The inventors subsequently developed the dNCPD algorithm to solve the optimization program empirically for achieving numerically accurate reconstruction of multibasis images from conventional data.Atty. Dkt. No. 05400-0082-PCT

[0139] The inventors performed numerical studies with simulation data of digital phantoms to verify that the dNCPD algorithm developed and its computer implementation can numerically accurately invert the non-linear data model and reconstruct multi-basis images from conventional data. Furthermore, the inventors conducted studies using real conventional data collected from physical and clinical phantoms of relevance to evaluate the stability’ of the dNCPD algorithm in reconstruction of multi-basis images and virtual monochromatic images (VMIs).

[0140] To facilitate the development of the dNCPD algorithm, included below is a summary of the well-established non-linear data model with the X-ray polychromaticity included in standard CT. One can use gj(b) to denote model data of ray j, where j = 1, 2,..., J and J is the total number of rays involved;mis the normalized effective spectrum for ray j at energy bin m, where m = 1, 2,..., M and M is the total number of energy bins; and pkm the mass attenuation coefficient at m of material k, where k = 1, 2,..., K, and K is the total number of basis-material types. Let vector bk of size I denote basis image k on the full image array of size I with entries bki, where i = 1, 2,..., I. For discussion convenience, one can use vector b of size KI to denote congregated basis image obtained by concatenating bk in the order of k.

[0141] In standard CT. one can express the well-established nonlinear data model in a discrete-to-discrete form as:= -In Xm=i Qjm exp(-A[ Xk=i ukmbk). (58)

[0142] where vector Aj of size I denotes row j of matrix A of size J x I, which depicts the discrete X-ray transform, with element aji representing the weight of the contribution of voxel i (i.e., column i) to model data of ray j (i.e., row j). The variable T indicates the transpose operation. It can be observed in data model Eq. (58) that basis images bk of interest are related non-linearly to model data gj(b) for polychromatic spectrum qjm. For a monochromatic spectrum, i.e., qjm has only one non-zero value among all energy bins, basis images bk are related then linearly to model data gj(b).

[0143] While it is possible that K > 3 basis images bk can be accurately and stably reconstructed directly from DE data collected with two distinct spectra, the developed algorithm can be used to invert the non-linear data model in Eq. (58) for numerically accurate and stable reconstruction of K > 3 basis images bk directly from conventional data collected in standard CT. As discussed above, the inverse problem from conventional data is generallyAtty. Dkt. No. 05400-0082-PCThighly ill-posed as conventional data contain an amount of information significantly less than that in DE data.

[0144] Improving the well-posedness of the data model is described below. In an attempt to alleviate the degree of ill-posedness of the inverse problem from conventional data, the inventors exploited, on one hand, the basis-region technique for reducing the number of unknowns involved in the reconstruction of multi-basis images and, on the other hand, the volume conservation (VC) constraint for augmenting conventional data.

[0145] Reduction of unknowns with basis-region images (basis regions and basis-region images) are discussed below. As the basis materials of the subject imaged are often confined within spatial regions smaller than the full image array, one can partition the full image array into L spatially complementary basis regions, and assume that basis region 1, where 1 = 1, 2,..., L, contains Ki known types of basis materials. It is specifically assumed herein that each of the basis regions contains KI < 2 known types of basis materials.

[0146] One can let Ri denote a diagonal matrix of size I x I, specifying basis region 1 inside or outside which voxel values are 1 or 0, respectively, and < Di denote a set of indices k’s of basis-material types contained in basis region 1. Basis region 1 can then be completely characterized by Ri and < Di. Using vector b'ik of size Ii < I, referred to as the basis-region image, to denote the spatial distribution of basis material k within basis region 1, it can be seen that Ii = tr(Ri) < I and Ki = card(Oi) < K, which are the total number of voxels of basisregion image bik and the total number of basis-material types contained in basis region 1.

[0147] Basis-region images, basis images, and VMIs are discussed below. One can write basis images bk on the full image array in terms of basis-region images b'ik as:bk=b'ik(59) where R'lis a matrix of size Ii x I obtained from diagonal matrix Ri by removing rows containing all zeros in Ri; and k is a set of indices 1’s of all the basis regions containing basis material k. The transpose of matrix R'i transforms basis-region images b'ik of size Ii within basis region 1 into an image on the full image array of size I.

[0148] One can use vector fm(b) of size I to denote the VMI at energy bin m, which can be written as fm(b) =∑Kk=1ukmbkin terms of basis images bk and is thus a function of basis image b. Using this relationship, along with Eq. (59), one can re-express the VMI at energy bin m as:Atty. Dkt. No. 05400-0082-PCTfm b ') ~ 2f=iSfce<r>iufcm^i7’ b[k, (60)

[0149] where vector fm(b') of size I denotes the VMI as a function of basis-region image b'ik. As mentioned, using Eq. (59), one can readily show fm(b') = fm(b).

[0150] Possible reduction of unknowns is discussed below. For discussion convenience, one can introduce congregated basis-region image b' by concatenating b'ik in the order of k and 1 and refer to b' simply as the basis-region image. It can be observed that the total number of voxel values, i.e., unknowns, in basis-region image b' is Nb = £f=1KtI which generally is smaller than KI, the total number of voxel values, i.e.. unknowns, in basis image b. As mentioned above, assuming in the w ork that each basis region contains Ki < 2 known types of basis materials, one has Nb < 21 that is significantly smaller than KI, the total number of voxel values, i.e., unknowns, in K> 3 basis images defined on the full image array.Therefore, the introduction of basis regions can effectively reduce the total number of unknowns involved in the multi (i.e., K > 3)-basis images.

[0151] Data augmentation with VC constraint is discussed below. It has been previously investigated and observed that, in the case of Ki < 2, the basis-region-based data model in Eq. (62) (below) can be inverted for accurate reconstruction of basis-region image b' directly from DE data. However, it is highly unlikely to directly invert Eq. (62) for numerically accurate reconstruction of b' only from conventional data without imposing additional, adequate constraints. Therefore, the inventors apply the VC constraint on the basis-region images within basis region L e 'Pvc containing Kic= 2 known types of basis materials as:VCie (b') =^lcbl'ck / pk- lIlc= o, (61 ) where 'Pvc denotes a set of indices L of basis regions to which the VC constraint is applied and Lc = card('Pvc) < L; pk depicts the known density of basis material k; and vector liic is of size lie with entries of value 1. It is noted that no VC constraint needs to be applied to a basis region containing only Ki = 1 known type of basis materials such as air or metal. The total Zic Itcexplicit constraints in Eq. (61) are exploited for augmenting conventional data for accurate and stable reconstruction of basis-region images.

[0152] A basis-region-based data model for image reconstruction is described below. In order to reconstruct from conventional data the basis region images w ith a reduced number of unknowns, a data model is needed that relates conventional data to basis-region images.Atty. Dkt. No. 05400-0082-PCTLetting g'j(b') denote model data in terms of basis region image b', and substituting Eq. (59) into Eq. (58). one obtains the basis-region-based data model with reduced unknowns as ^■(b') = -ZnZ"=iiV7.m(y), (62) whereNjm(.b ) X( = 1 ^lk)- (6->)

[0153] It is noted that summations Xt=i XiG’i'kand Xf=i Ske<t>(areequivalent, as two ways to traverse the basis-region images over index 1 or k first. Again, as mentioned, using Eqs. (58) and (59), one can readily show g'j(b') = gj(b).

[0154] For discussion convenience below, one can introduce vector g'(b') of size J with entries g'j(b') and then re-express the ray-based data model in Eq. (62) in a matrix-vector form as:g'^b') = + hg(b', b'), (64) where H^b'^b' is the linear term when g'(b') is expanded into the Taylor series at point b' selected; matrix H(b') of size J x Nb has elements, / " VA T,m=lukmNjm(Pr) vU / h^b)=(65)

[0155] The symbol r'lii' is the element at row i and column i' of matrix R'. It is noted that Eq. (64) can simply be re-ordered to yield the defining expression for g(b', b') as:hg(b', b'} = g'(b'^) — H(b'^b', (66) which includes the non-linear component of the data model in Eq. (62). Therefore, the task is to invert Eq. (64) (or, equivalently, Eq. (62)) for numerically accurate reconstruction of basisregion images b' from which basis images bk can readily be obtained by using Eq (59).

[0156] Optimization-based image reconstruction is discussed below. With the basisregion images of reduced number of unknowns and VC constraint for data augmentation discussed above, the algorithm development is in order for numerically inverting the basisregion-based data model in Eq. (64) (or, equivalently, Eq. (62)) to achieve accurate and stable reconstruction of basis-region images from conventional data collected with a single spectrum in standard CT.Atty. Dkt. No. 05400-0082-PCT

[0157] The inventors formulated the reconstruction problem of basis-region image b' (i.e., the inverse problem of Eq. (64)) from conventional data as a constrained optimization program:b’* = argb,min^ ||[^[M]- hg(b', b')] - H(b')b'\\22(Eq. 67) 5- 1- ll(|V / ^i(b')|)||1< || (| V / ^2(b') |) Hi < ym2VClc(b') = 0 for lc= b' > 0,where vector g| |of size J denotes conventional data measured with a single spectrum in a standard CT scan; ||(| Vfmi (b')|)| 11 the total-variation (TV) of VMI fmi at energy bin mi selected; Lcis the number of basis regions subject to the VC constraint; and b' is an input independent of b'.

[0158] The optimization program of interest in Eq. (67) is nonconvex as a result of using the non-linear data model (i.e., Eq. (64)) in the objective function (i.e., line 1 in Eq. (67)). Because it is unclear if an algorithm can be developed for mathematically exactly converging the non-convex optimization program in Eq. (67). the inventors took the approach, as described below, to developing an algorithm to empirically converge the non-convex optimization program in Eq. (67) for numerically accurate and stable reconstruction of basisregion images. Fig. 29 depicts pseudo-codes of the dNCPD algorithm for empirically solving Eq. (67) in accordance with an illustrative embodiment.

[0159] Results of extensive numerical studies, including those presented below, indicate that the dNCPD algorithm can numerically accurately converge to the solution of the non-convex optimization program in Eq. (67). In the numerical studies below, reconstruction of basis-region image b' is obtained when convergence conditions are satisfied numerically in terms of the computer single-precision floating-point error. Numerically converged solution b'* in the pseudo-codes of the dNCPD algorithm is used as basis-region image b' reconstructed from which the basis images and VMIs can then be obtained by use of Eqs. (59) and (60).

[0160] Numerical studies were carried out in this work that are composed of a verification study and an evaluation study on the numerical accuracy and stability of the dNCPD algorithm. In the verification study, simulated data are generated from digitalAtty. Dkt. No. 05400-0082-PCTphantoms by use of non-linear data model Eq. (62) and thus are consistent completely with the data model. They are thus used for verifying that the dNCPD algorithm can accurately invert the non-linear data model in Eq. (62) through empirically converging the non-linear optimization program in Eq. (67). Meanwhile, in the evaluation study, real data of physical phantoms, which are collected by use of clinical DECT scanners, are inconsistent with the non-linear data model in Eq. (62) as they contain noise, scatter, and other physical factors. Therefore, real data are thus used for evaluating the stability of dNCPD algorithm in inverting the non-linear data model in Eq. (62) through empirically converging the non-convex optimization program in Eq. (67).

[0161] The inventors carried out a verification study to demonstrate that the dNCPD algorithm can empirically converge the non-convex optimization program in Eq. (67) for inverting the non-linear data model in Eq. (62) to reach numerically accurate reconstruction of multi (K > 2)-basis images, only from conventional data of a digital chest phantom.

[0162] The digital chest phantom contains a total of K = 6 basis materials, i.e., air, water, bone, 20-mg / ml iodine solution, titanium, and stainless steel. A set of 5 spatially complementary basis regions, as depicted in Fig. 24, partitions the full image array of 200x256 0.14-cm square pixels. It is assumed that basis materials air, bone, 20-mg / ml iodine solution, titanium, and stainless steel are confined, respectively, within the 5 complementary basis regions shown in Fig. 24 from left to right, and basis material water distributes within the first three basis region, which is the full image array minus the regions containing titanium and air-steel. The titanium basis region (column 4) contains only the titanium basis material, while air is also allowed in the air-steel basis region (column 5), for the needle might contain air. Therefore, each basis region contains up to 2 basis materials. Furthermore, as the phantom includes K = 6 basis materials, it thus has K = 6 basis images on the full image array, referred to as the truth basis images, each of which contains only single basis material, as shown in column 1 of Fig. 25.

[0163] Generation of simulated conventional data is discussed below. Using the truth basis images Eq. (62), the inventors generated simulated, conventional data with each of two distinct spectra at 1440 views uniformly distributed over 2n with a circular fan-beam geometry’ that has a source-to-rotation (SOR) distance of 100 cm, a source-to-detector (SOD) distance of 150-cm, and a linear detector that includes 896 detector bins of 0.08-cm bin size. Specifically, the two spectra of 80 kV and 140 kV are produced by use of the TASMICAtty. Dkt. No. 05400-0082-PCTpackage (Hernandez and Boone, 2014), whereas the mass attenuation coefficients of the basis materials are either looked up from the NIST database or generated using the NIST XCOM tool.

[0164] Parameter selection. The dNCPD algorithm is used to reconstruct basis-region images from simulated conventional data of the digital chest phantom. In the reconstruction, geometrical and spectral information and truth basis regions used for the data generation are also used in the algorithm; whereas the TV constraint parameters. ymiandm2 are computed from the truth VMIs at 50 and 100 keV, which can readily be obtained from the truth basis images in column 1 of Fig. 25.

[0165] Quantitative convergence results. Fig. 26 displays the numerical convergence properties of the dNCPD algorithm in terms of convergence metrics as functions of iteration n from simulated conventional data generated with the 80-kV spectrum. It can be observed that the necessary convergence conditions are satisfied numerically in terms of the computer single-precision floating-point error, confirming that the dNCPD algorithm can empirically converge the non-convex optimization program in Eq. 67. Furthermore, as image-quality metric Dbin the last panel in Fig. 26 shows, the condition on the reconstruction of basis images is satisfied numerically, verifying that for the conditions in the study, the dNCPD algorithm can numerically accurately invert the non-linear data model in Eq. (62) (or, equivalently, Eq. (58)). Similar convergence results are obtained for simulated conventional data generated with the 140-kV spectrum.

[0166] Basis images reconstructed. Also shown in columns 2 and 3 of Fig. 25 are basis images of air, water, bone, 20-mg / ml iodine solution, titanium, and stainless steel reconstructed, respectively, from simulated 80- and 140-kV conventional data. It can be observed that the basis images are visually identical to the truth basis images (column 1), thus corroborating the convergence results in Fig. 26. The results of basis image reconstruction, along with the convergence evidence, verify that the dNCPD algorithm can numerically converge Eq. (67) for numerically accurate reconstruction of basis images, or, equivalently, numerically accurate inversion of the non-linear data model in Eq. (62) only from conventional data in standard CT.

[0167] Once the dNCPD algorithm is verified for numerically accurate reconstruction of multi-basis images (or, equivalently, numerically accurate inversion of the non-linear data model in Eq. (62)) from simulated conventional data consistent with the data model, one canAtty. Dkt. No. 05400-0082-PCTthen evaluate the stability of the dNCPD algorithm for reconstruction of multi-basis images from real conventional data, which contain necessarily physical factors, such as noise, scatter, and spectrum inaccuracy, that are not considered in, and thus inconsistent with, the non-linear data model in Eq. (62).

[0168] A physical DE phantom study was also conducted. In this real-data study, the inventors used the widely-used, physical DE phantom that includes a solid water disk with an approximate size of an average pelvis, 7 inserts of iodine contrast solutions at 2, 2.5, 5, 7.5, 10, 15, 20 mg / ml concentrations, and 7 inserts of calcium solutions at 50. 100, 200, 300, 400, 500, 600 mg / ml concentrations, as well as 8 air gap holes. The iodine and calcium inserts with various know n concentration levels can be used for quantitative analysis including estimation of LAC and iodine concentration. The basis images and VMIs of the physical DE phantom are reconstructed on a full image array of 512 x 512 0.08-cm square pixels.

[0169] Real conventional data of the physical DE phantom were collected with 80-kV and 135-kV spectra by use of a clinical CT scanner; and the scanner has multiple curved detector rows each of which includes 896 equal-angular detector bins. The central curved detector row has SOR and SOD distances of 60 cm and 107 cm, thus forming a fan angle of 49°. In the reconstruction studies below, the inventors used real conventional data collected on the central curved detector row in axial mode at 1200 views evenly distributed over 2rt. For the clinical scanner, the 80- and 135-kV spectra are estimated, respectively, for use in the image reconstruction.

[0170] Parameter selection. As discussed, the physical DE phantom includes four (K = 4) types of basis materials, including air, w aler. 600-mg / ml calcium solution, and 20-mg / ml iodine solution; and the inventors partitioned the full image array into three basis regions each of which contains two types of basis materials. It can be observed that the optimization program in Eq. (67) is specified by several parameters, including TV-constraint parameters ymi& ym2 and also the selection of basis regions. While knowledge of these parameters is known in the verification study because the truth digital chest phantom is known, in a real-data study, the parameters need to be estimated empirically because their respective ‘‘truths’' are generally unknown or even undefined precisely. Instead, given a set of real conventional data, one can perform multiple reconstructions with different sets of ymi& ym2 and basis regions selected and then choose the set that yields the basis images of visually minimalAtty. Dkt. No. 05400-0082-PCTartifacts for yielding the “optimal’' image reconstruction from real conventional data. In the physical DE-phantom study, the basis regions so selected are shown in Fig. 27.

[0171] Reconstruction of basis images and VMIs from conventional data. As described above, given each of conventional data sets collected with either 80-kV or 135-kV spectrum, the inventors selected the constraint parameters and the basis regions shown in Fig. 27, that are needed in the dNCPD algorithm for reconstruction of basis region images. Using the basis-region images reconstructed in Eqs. (59) and (60), one can readily obtain the basis images and then VMI. Fig. 1 A displays the basis images and VMI at 80 keV reconstructed, respectively, from 80-kV (column 1) and 135-kV (column 2) conventional data sets. It can be observed that the dNCPD algorithm can stably reconstruct multi-basis images and also VMIs directly from real conventional data collected with either 80-kV or 135-kV spectrum.

[0172] Conversely, the two sets of data collected with 80-kV and 135-kV spectra can also be used to form DE data, referred to as 80& 135-kV DE data, similar to that collected in DECT. Therefore, for comparison, the inventors also apply an existing algorithm to reconstructing basis images and VMIs from 80&135-kV DE data formed. In column 3 of Fig.1A, the basis images and VMI at 80 keV reconstructed from 80&135-kV DE data are displayed. It can be observed in Fig. 1 A that basis images and VMIs reconstructed from conventional data collected with either 80-kV or 135-kV spectrum appear to be visually comparable to the corresponding basis images and VMIs reconstructed from 80&135-kV DE data.

[0173] Quantitative analysis of the reconstructed images. In addition to visual inspection, the inventors also conducted quantitative analysis of the basis images and VMIs that were reconstructed. The inventors first quantitatively analyzed estimates of iodine concentrations (ICs) within 6 regions-of-interest (ROIs) selected in Fig. 1A, which center around solid water, iodine solutions of 2, 5, 10, and 20 mg / ml, and calcium solution of 50 mg / ml, respectively. Considering an affine relationship between the concentration level and the pixel values in the iodine basis image, one can obtain IC values at the pixels in the iodine-basis image within the ROIs indicated in Fig. 1A and then average the IC values over each of the respective ROIs. Using the average ICs obtained, one can then compute their biases and standard deviations (SDs) relative to the respective reference ICs provided by the manufacturer of the physical DE phantom, which are summarized in the table of Fig. IC.Atty. Dkt. No. 05400-0082-PCT

[0174] It can be observed that the biases and SDs of the IC estimates within the iodine basis images reconstructed from either 80-kV or 135-kV conventional data are comparable and that they are also comparable to those obtained from 80& 135-kV DE data. The results thus demonstrates quantitatively accuracy and stability of the dNCPD algorithm in reconstruction of the basis images from conventional data collected with a single spectrum in standard CT.

[0175] A quantitative analysis of the VMIs reconstructed over energy range 30 ~ 140 keV was also performed by computing average LACs, along with their respective SDs, over 5 ROIs, labeled 0, 1, 4, 5, and 6 in Fig. 1 A, containing solid water, iodine solutions of 2 and 20 mg / ml, and calcium solutions of 50 and 600 mg / ml. The LACs, along with their SDs, averaged over their respective ROIs, as displayed in Fig. IB, are estimated from VMIs reconstructed from 80-kV or 135-kV conventional data. The LACs estimated are observed to agree generally well with the respective reference values (solid curves in Fig. IB) computed from the manufacturer’s spec sheet. The results reveal quantitative accuracy and stability of the dNCPD algorithm in VMI reconstruction across the energy range from conventional data collected with a single spectrum in standard CT. It can also be observed that the results in the left and middle panels of Fig. IB obtained from the conventional data appear to be comparable to those in the right panel of Fig. IB obtained from the 80&135-kV DE data.

[0176] The clinical abdomen phantom. The inventors also carried out an additional real-data study using conventional data collected from a clinical abdomen phantom with 80-kV or 135-kV spectrum. The clinical abdomen phantom is chosen because it possesses realistic anatomic complexity similar to that of a human abdomen. Considering the anatomic composition of the clinical abdomen phantom, four (K = 4) basis materials, including air, adipose, bone, and 20-mg / ml iodine solution, are selected for reconstructing multi (K = 4)-basis images directly from 80-kV or 135-kV conventional data. This study is used for further investigating the stability of dNCPD algorithm with respect to different levels of anatomic complexity and clinical relevance. The basis images and VMIs of the phantom are reconstructed on a full image array of 432 × 656 0.08-cm square pixels.

[0177] Real conventional data collected. Real conventional data of the clinical abdomen phantom are collected with 80-kV and 135-kV spectra at 1200 views evenly distributed over 2TI by use of the same clinical CT scanner as that in the study with the physical DE phantom.Atty. Dkt. No. 05400-0082-PCTThe 80- and 135-kV spectra estimated above are used also in image reconstruction of the clinical abdomen phantom.

[0178] Parameter selection. Considering four (K = 4) types of basis materials, including air. adipose, bone, and 20-mg / ml iodine solution, one can select basis regions that partition the full image array and assume two types of the basis materials within each of basis regions. One can take the same approach as that described above to selecting TV-constraint parameters ymi& ym2 and basis images. Namely, giving each set of real conventional data collected with 80-kV and 135-kV spectra, one can perform multiple reconstructions with different sets of TV constraint parameters ymi& ym2 and basis regions, and then choose the set that yields the basis images of visually minimal artifacts for use the final reconstruction of the basis images. In the study, the basis regions so chosen are shown in Fig. 28.

[0179] Reconstruction of basis images and VMIs from conventional data. As described above, for each set of 80-kV and 135-kV conventional data, the inventors selected the constraint parameters, including the basis regions shown in Fig. 28. for reconstruction of basis-region images by use of the dNCPD algorithm. Using the basis-region images reconstructed in Eqs. (59) and (60), one can readily obtain the basis images and then VMI at 80 keV, which are displayed in Fig. 2. It can be observed that the dNCPD algorithm can stably reconstruct multi-basis images and also VMIs directly from real conventional data collected with either 80-kV or 135-kV spectrum. Because no knowledge of the ground truth of the clinical abdomen phantom is available, no quantitative analysis of the reconstruction in terms of IC and LAC estimations, similar to that discussed above, is performed.

[0180] Again, the two sets of data of the clinical abdomen phantom collected with 80-kV and 135-kV spectra can be used to form DE data, also referred to as 80&135-kV DE data, similar to that collected in DECT. Therefore, for comparison, one can also apply the same existing algorithm to reconstructing basis images and VMIs from 80&135-kV DE data formed. In column 3 of Fig. 2, the basis images and VMI at 80 keV reconstructed from 80&135-kV DE data are displayed. It can be observed in Fig. 2 that basis images and VMIs reconstructed from conventional data collected with either 80-kV or 135-kV spectrum appear to be visually comparable to the corresponding basis images and VMIs reconstructed from 80&135-kV DE data.

[0181] Thus, described herein is the development of the dNCPD algorithm to invert the non-linear data model in standard CT for numerically accurate and stable reconstruction ofAtty. Dkt. No. 05400-0082-PCTmulti (> 2)-basis images and VMIs directly from conventional data collected with a single spectrum in standard CT. The development of the dNCPD algorithm is enabled by the exploitation of the basis-region technique for reducing the number of voxels in basis images to be reconstructed and the VC constraint for effectively augmenting conventional data.

[0182] Numerical studies were conducted on the performance of the dNCPD algorithm by using a qualitative metric of visual inspection and then quantitative metrics of (spatial-averaging) bias and variance of estimation of LACs and IC concentrations. Specifically, the accuracy of the dNCPD algorithm and its computer implementation are verified numerically first in a quantitative study with simulated conventional data of digital phantoms. Following the verification study, the inventors evaluated and demonstrated the stability of the dNCPD algorithm for numerically accurate reconstruction of the basis images from real conventional data collected with a low- or a high-kVp spectrum in standard CT.

[0183] The problem of image reconstruction (or, equivalently, of the inversion of the non-linear data model) is formulated as an optimization program in Eq. (67), which is non-convex as a consequence of its inclusion of the non-linear data model in CT imaging. In an attempt to empirically converge the non-convex optimization program, the inventors first convexify it to obtain a convexified optimization program so that a PD algorithm can be derived for accurately solving it; and the inventors then tailored the PD algorithm derived to obtain the dNCPD algorithm for empirically solving the non-convex optimization program.

[0184] As such, the pseudo-code structures of the dNCPD and PD algorithms are identical, except for the additional non-linear correction step involving Ag(b', b') in the former. This indeed also provides a convenient empirical check of the convergence of the dNCPD algorithm when monochromatic X-rays are considered. In this case, the non-linear data model turns into the standard linear data model, the non-convex optimization program degenerates into a convex optimization program, and the dNCPD algorithm thus becomes the PD algorithm. Furthermore, the convergence properties of the dNCPD algorithm, in terms of the convergence metrics shown in Fig. 26, are also similar to those of the PD algorithm under the same data conditions.

[0185] The dNCPD algorithm developed can be applied to reconstructing numerically accurately and stably multi (> 2)-basis images only from conventional data assuming that each of the basis regions contains KI < 2 types of basis materials. This assumption is generally reasonable because photoelectric effect and Compton scattering are the twoAtty. Dkt. No. 05400-0082-PCTdominant interaction mechanisms in the X-ray-energy range of standard diagnostic CT, indicating two degrees of freedom in the decomposition that thus allow for up to two types of basis materials in the basis region. Meanwhile, in the presence of K-edge materials, a representative basis material is included, such as 20-mg / ml iodine solution in the work, and basis regions are selected such that the iodine solution basis material is contained with a low-attenuating basis material, such as water or adipose, to increase the expansion space afforded by the decomposition model in that basis region.

[0186] On the other hand, the KI < 2 assumption means that certain DECT imaging applications where an image voxel may need to be decomposed into 3 basis materials, such as fat, liver tissue, and blood for liver fat quantification, may not be replaced always appropriately with standard scans only using a single spectrum. Parameter selection is critical to the performance of any reconstruction algorithms. In the work, the basis regions and TV constraint parameters are selected as their respective truth in the simulated-data study and by surveying the parameter space based on the empirical metric of visualization in the real-data study. Some guidance in the parameter space search can be provided in order to reduce the search domain. For example, TV values can be computed from FBP-reconstructed images from either low-kV or high-kV conventional data and then used as reference values for searching the TV-constraint parameters. The FBP-reconstructed images can also be segmented using simple hard-thresholding to provide a reference for selection of the basisregion partition. However, unlike the case in DECT with the VC constraint, some prior knowledge of the distribution of K-edge materials needs to be incorporated, as it might be challenging to separate iodine and bone based only upon pixel values.

[0187] The dNCPD algorithm can readily be generalized to reconstruct basis images and VMIs from conventional data collected with standard cone-beam CT and be applied directly to conventional data collected at sparse views in fan- and cone-beam CT as the non-convex optimization program and the dNCPD algorithm are formulated in terms of the geometry of individual X-rays in fan- and conebeam CT. It would also be interesting to tailor the dNCPD algorithm by including constraints on, e.g., image directional TVs, to address the reconstruction problem of basis images and VMIs from conventional data collected only over a limited-angular range. Finally, the work can be of practical implication as it reveals the possibility of obtaining multi-basis images and VMIs only from conventional data in standard CT, instead of data collected in DE, multi-spectra, or photon counting CT, which generallyAtty. Dkt. No. 05400-0082-PCTinvolve either additional scanning effort and / or additional, unique hardware components / systems.

[0188] In an illustrative embodiment, any of the operations described herein can be performed by a computing system that includes a processor, a memory, a user interface, a transceiver (e.g., a receiver and transmitter), etc. The operations described herein can be implemented as computer-readable instructions that are stored in the memory’. Upon execution of the computer-readable instructions by the processor, the computing system performs the various operations. As an example, Fig. 30 depicts a computing system for performing image reconstruction in accordance with an illustrative embodiment.

[0189] The system of Fig. 30 includes a computing device 3000 that has a processor 3005, an operating system 3010, a memory 3015, an input / output (I / O) system 3020, a network interface 3025, and an image reconstruction application 3030. In alternative embodiments, the computing device 3000 may include fewer, additional, and / or different components. The components of the computing device 3000 communicate with one another via one or more buses or any other interconnect system. The computing device 3000 can be any type of computing device (e.g., smartphone, tablet, laptop, desktop, etc.), including a dedicated standalone computing system that is designed to perform the calculations and generate the images described herein.

[0190] The system of Fig. 30 also includes an imaging system 3040 that is in direct or indirect (i.e., via the network 3035) communication with the computing device 3000. The imaging system 3040 can be an x-ray imaging system (e.g., general radiographic imaging, fluoroscopic imaging, mammographic imaging, angiographic imaging, etc.), a CAT imaging system, a tomosynthesis imaging system, a laminographic imaging system, etc. In an alternative embodiment, the computing device 3000 may be incorporated into the imaging system 3040.

[0191] The processor 3005 can be in electrical communication with and used to control the imaging system 3040. Alternatively, the processor may receive information from the imaging system 3040 without controlling the imaging system 3040. The processor can also be used to execute the image reconstruction application 3030. The processor 3005 can be any type of computer processor known in the art, and can include a plurality’ of processors and / or a plurality of processing cores. The processor 3005 can include a controller, a microcontroller, an audio processor, a graphics processing unit, a hardware accelerator, aAtty. Dkt. No. 05400-0082-PCTdigital signal processor, etc. Additionally, the processor 3005 may be implemented as a complex instruction set computer processor, a reduced instruction set computer processor, an x86 instruction set computer processor, etc. The processor 3005 is used to run the operating system 3010, which can be any type of operating system.

[0192] The operating system 3010 is stored in the memory 3015, which is also used to store programs, received imaging data from the imaging system 3040, mathematical constants and algorithms, network and communications data, peripheral component data, the image reconstruction application 3030, and other operating instructions. The memory 3015 can be one or more memory systems that include various types of computer memory such as flash memory, random access memory (RAM), dynamic (RAM), static (RAM), a universal serial bus (USB) drive, an optical disk drive, a tape drive, an internal storage device, a nonvolatile storage device, a hard disk drive (HDD), a volatile storage device, etc. In some embodiments, at least a portion of the memory 3015 can be in the cloud to provide cloud storage for the system. Similarly, in one embodiment, any of the computing components described herein (e.g., the processor 3005, etc.) can be implemented in the cloud such that the system can be run and controlled through cloud computing.

[0193] The I / O system 3020 is the framework which enables users and peripheral devices to interact with the computing device 3000. The I / O system 3020 can also include one or more speakers, one or more displays, one or more microphones, a keyboard, a mouse, one or more buttons or other controls, etc. that allow the user to interact with and control the computing device 3000. The I / O system 3020 also includes circuitry and a bus structure to interface with peripheral computing devices such as power sources, universal service bus (USB) devices, data acquisition cards, peripheral component interconnect express (PCIe) devices, serial advanced technology attachment (SATA) devices, high definition multimedia interface (HDMI) devices, proprietary connection devices, etc.

[0194] The network interface 3025 includes transceiver circuitry (e.g., a transmitter and a receiver) that allows the computing device 3000 to transmit and receive datato / from other devices such as the imaging system 3040, remote computing systems, servers, websites, etc. The network interface 3025 enables communication through the network 3035. which can be one or more communication networks. The network 3035 can include a cable network, a fiber network, a cellular network, a wi-fi network, a landline telephone network, a microwaveAtty. Dkt. No. 05400-0082-PCTnetwork, a satellite network, etc. The network interface 3025 also includes circuitry' to allow device-to-device communication such as Bluetooth® communication.

[0195] The image reconstruction application 3030 can include software and algorithms in the form of computer-readable instructions which, upon execution by the processor 3005, performs any of the various operations described herein such as receiving image data from the imaging system 3040, controlling the imaging system 3040, solving equations, formulating reconstruction as a non-convex optimization problem, solving the non-convex optimization problem, determining constraints, determining basis path lengths, reducing voxels, generating images, etc. The image reconstruction application 3030 can utilize the processor 3005 and / or the memory 3015 as discussed above. In an alternative implementation, the image reconstruction application 3030 can be remote or independent from the computing device 3000, but in communication therewith.

[0196] Fig. 31 depicts a computing system for performing pathlength estimation in accordance with an illustrative embodiment. Similar to the system of Fig. 30. the system of Fig. 31 includes a computing device 3100 that has a processor 3105, an operating system 3110, a memory 3115, an input / output (I / O) system 3120, and a network interface 3125. The computing device 3100 also includes a pathlength estimation application 3130. In alternative embodiments, the computing device 3100 may include fewer, additional, and / or different components. The components of the computing device 3100 communicate with one another via one or more buses or any other interconnect system. The computing device 3100 can be any type of computing device (e.g., smartphone, tablet, laptop, desktop, etc.), including a dedicated standalone computing system that is designed to perform the calculations and pathlength estimations described herein.

[0197] The system of Fig. 31 also includes a pathlength (or proj ection) imaging system 3140 that is in direct or indirect (i.e., via the network 3135) communication with the computing device 3100. The pathlength imaging system 3140 can be any type of x-ray projection imaging system used to obtain chest x-rays, mammography scans, bone scans, head scans, scans of any other body parts, security scans, industrial scans, etc. In an alternative embodiment, the computing device 3100 may be incorporated into the pathlength imaging system 3140. In an illustrative embodiment, pathlength estimation can be performed from a single view of conventional projection data collected from the pathlength imaging system 3140.Atty. Dkt. No. 05400-0082-PCT

[0198] The processor 3105 can be in electrical communication with and used to control the pathlength imaging system 3140. Alternatively, the processor may receive information from the pathlength imaging system 3140 without controlling the pathlength imaging system 3140. The processor can also be used to execute the pathlength estimation application 3130. Similar to the processor 3005 of Fig. 30, the processor 3105 can be any type of computer processor known in the art, and can include a plurality of processors and / or a plurality of processing cores. The processor 3105 is used to run the operating system 3110. which can be any type of operating system.

[0199] The operating system 3110 is stored in the memory 3115, which is also used to store programs, received imaging data from the pathlength imaging system 3140, mathematical constants and algorithms, network and communications data, peripheral component data, the pathlength estimation application 3130, and other operating instructions. The memory 3115 can be one or more memory systems that include various types of computer memory as discussed above with respect to the memory 3015 of Fig. 30. In some embodiments, any of the computing components of Fig. 31 (e.g., the processor 3105, the memory 3115, etc.) can be implemented in the cloud such that the system can be run and controlled through cloud computing. The I / O system 3120 and the network interface 3125 can operate in the same fashion as the I / O system 3020 and the network interface 3025 described with reference to Fig. 30.

[0200] The pathlength estimation application 3130 can include software and algorithms in the form of computer-readable instructions which, upon execution by the processor 3105, performs any of the various operations described herein such as receiving image data from the pathlength imaging system 3140, controlling the pathlength imaging system 3140, solving equations, inverting a non-linear BPL data model, estimating basis path lengths, determining virtual pathlengths, etc. The pathlength estimation application 3130 can utilize the processor 3105 and / or the memory 3115 as discussed above. In an alternative implementation, the pathlength estimation application 3130 can be remote or independent from the computing device 3100, but in communication therewith.

[0201] The word "illustrative" is used herein to mean serving as an example, instance, or illustration. Any aspect or design described herein as "illustrative" is not necessarily to be construed as preferred or advantageous over other aspects or designs. Further, for the purposes of this disclosure and unless otherwise specified, "a" or "an" means "one or more.”Atty. Dkt. No. 05400-0082-PCT

[0202] The foregoing description of illustrative embodiments of the invention has been presented for purposes of illustration and of description. It is not intended to be exhaustive or to limit the invention to the precise form disclosed, and modifications and variations are possible in light of the above teachings or may be acquired from practice of the invention. The embodiments were chosen and described in order to explain the principles of the invention and as practical applications of the invention to enable one skilled in the art to utilize the invention in various embodiments and with various modifications as suited to the particular use contemplated. It is intended that the scope of the invention be defined by the claims appended hereto and their equivalents.

Claims

Atty. Dkt. No. 05400-0082-PCTWHAT IS CLAIMED IS:

1. An image reconstruction system comprising:a memory configured to store image data obtained from an imaging system; anda processor operatively coupled to the memory and configured to: formulate reconstruction of the image data as a non-convex optimization program;generate a solution to the non-convex optimization problem; and reconstruct the image data to generate a multi-basis image based on the solution to the non-convex optimization problem.

2. The system of claim 1, wherein the image data includes multiple polychromatic spectra.

3. The system of claim 2, wherein the image data originates from an single energy x-ray imaging process or a dual energy x-ray imaging process.

4. The system of claim 1, wherein the processor inverts a non-linear basis path length data model to estimate basis path lengths for the image data, and wherein the solution to the non-convex optimization problem is based at least in part on the estimated basis path lengths.

5. The system of claim 4, wherein the processor converts a volume conservation constraint into a constraint on the basis path lengths to estimate the basis path lengths.

6. The system of claim 4, wherein the processor converts a mass conservation constraint into a constraint on the basis path lengths to estimate the basis path lengths.

7. The system of claim 4, wherein the non-convex optimization problem is based on the inverted non-linear basis path length data model.

8. The system of claim 1, wherein the processor generates the solution by numerically converging the non-convex optimization problem.Atty. Dkt. No. 05400-0082-PCT9. The system of claim 1, wherein the image data is collected from a single spectrum in an imaging process that utilizes x-rays.

10. The system of claim 1, wherein the processor uses a basis-region technique to reduce a number of voxel values in the image data.

11. A method of performing image reconstruction, the method comprising:storing, in a memory, image data obtained from an imaging system; formulating, by a processor operatively coupled to the memory, reconstruction of the image data as a non-convex optimization program;generating, by the processor, a solution to the non-convex optimization problem; andreconstructing, by the processor, the image data to generate a multi-basis image based on the solution to the non-convex optimization problem.

12. The method of claim 1, wherein the image data includes multiple polychromatic spectra.

13. The method of claim 12, wherein the image data originates from a single energy x-ray imaging process or a dual energy x-ray imaging process.

14. The method of claim 1, further comprising inverting, by the processor, a non-linear basis path length data model and using the inverted non-linear basis path length data model to estimate basis path lengths for the image data, wherein the solution to the non-convex optimization problem is based at least in part on the estimated basis path lengths.

15. The method of claim 14, further comprising converting, by the processor, a volume conservation constraint into a constraint on the basis path lengths to estimate the basis path lengths.

16. The method of claim 14, further comprising converting, by the processor, a mass conservation constraint into a constraint on the basis path lengths to estimate the basis path lengths.Atty. Dkt. No. 05400-0082-PCT17. The method of claim 14, wherein the non-convex optimization problem is based on the inverted non-linear basis path length data model.

18. The method of claim 11, further comprising generating, by the processor, the solution by numerically converging the non-convex optimization problem.

19. The method of claim 11, further comprising collecting, by the processor, the image data from a single spectrum in an imaging process that utilizes x-rays.

20. The method of claim 11, further comprising using, by the processor, a basis-region technique to reduce a number of voxel values in the image data.

Images